357 if (multiple_subdomains && tangle_mesh)
358 libmesh_not_implemented_msg(
359 "Arbitrary mesh tangling with multiple subdomains is not supported.");
363 const auto dim = ref_elem->dim();
370 unsigned int n_elems_per_side = 4 / Elem::type_to_default_order_map[type];
379 mesh, n_elems_per_side, n_elems_per_side, 0., 1., 0., 1., type);
397 libmesh_error_msg(
"Unsupported dimension " <<
dim);
401 DistortHyperCube dh(
dim);
415 const Real dr = 1. / (n_elems_per_side * Elem::type_to_default_order_map[type]);
421 Real dist1_closest = std::numeric_limits<Real>::max();
423 Real dist2_closest = std::numeric_limits<Real>::max();
425 for (
const auto * node :
mesh.local_node_ptr_range())
427 const auto dist1 = (*node - p1).
norm();
428 if (dist1 < dist1_closest)
430 dist1_closest = dist1;
431 node1_id = node->id();
434 const auto dist2 = (*node - p2).
norm();
435 if (dist2 < dist2_closest)
437 dist2_closest = dist2;
438 node2_id = node->id();
443 unsigned int node1_rank;
447 unsigned int node2_rank;
452 if ((node1_id != DofObject::invalid_id) && (node2_id != DofObject::invalid_id))
459 const Point diff = node1 - node2;
463 node1 -= tangle_damping_factor * diff;
464 node2 += tangle_damping_factor * diff;
473 CPPUNIT_ASSERT(unsmoothed_info.mesh_is_tangled);
477 std::unordered_map<subdomain_id_type, Real> distorted_subdomain_volumes;
478 if (multiple_subdomains)
481 for (
auto * elem :
mesh.active_element_ptr_range())
483 unsigned int subdomain_id = 0;
485 if (elem->vertex_average()(d) > 0.5)
487 elem->subdomain_id() += subdomain_id;
493 for (
auto * elem :
mesh.active_element_ptr_range())
494 distorted_subdomain_volumes[elem->subdomain_id()] += elem->volume();
499 libmesh_error_msg_if(elem_orders.size() != 1,
500 "The variational smoother cannot be used for mixed-order meshes!");
504 const auto scale_factor = *elem_orders.begin() * (type_is_pyramid ? 2 * 4 : 1);
507 auto node_distortion_is = [&n_elems_per_side, &
dim, &boundary_info, &scale_factor, &type_is_prism](
510 std::vector<boundary_id_type> boundary_ids;
511 boundary_info.boundary_ids(&node, boundary_ids);
514 const auto num_dofs =
dim - boundary_ids.size();
534 std::size_t num_zero_or_one = 0;
536 bool distorted =
false;
539 const Real r = node(d);
540 const Real R = r * n_elems_per_side * scale_factor;
541 CPPUNIT_ASSERT_GREATER(-distortion_tol * distortion_tol, r);
542 CPPUNIT_ASSERT_GREATER(-distortion_tol * distortion_tol, 1 - r);
544 bool d_distorted = std::abs(R - std::round(R)) > distortion_tol;
545 if (type_is_prism && (scale_factor == 3))
550 const Real R_prism = R / scale_factor * 2;
551 const bool d_distorted_prism =
552 std::abs(R_prism - std::round(R_prism)) > distortion_tol;
553 d_distorted &= d_distorted_prism;
555 distorted |= d_distorted;
559 CPPUNIT_ASSERT_GREATEREQUAL(
dim - num_dofs, num_zero_or_one);
565 return distorted == distortion;
569 for (
auto node :
mesh.node_ptr_range())
570 CPPUNIT_ASSERT(node_distortion_is(*node,
true));
580 SquareToParallelogram stp;
583 else if (type_is_prism)
588 CubeToParallelepiped ctp;
599 ParallelogramToSquare pts;
602 else if (type_is_prism)
607 ParallelepipedToCube ptc;
611 if (multiple_subdomains)
616 std::unordered_map<subdomain_id_type, Real> smoothed_subdomain_volumes;
617 for (
auto * elem :
mesh.active_element_ptr_range())
618 smoothed_subdomain_volumes[elem->subdomain_id()] += elem->volume();
620 for (
const auto & pair : distorted_subdomain_volumes)
622 const auto & subdomain_id = pair.first;
625 libmesh_map_find(smoothed_subdomain_volumes, subdomain_id),
632 std::set<dof_id_type> nodes_checked;
633 for (
const auto * elem :
mesh.active_element_ptr_range())
635 for (
const auto local_node_id :
make_range(elem->n_nodes()))
637 const auto & node = elem->node_ref(local_node_id);
638 if (nodes_checked.find(node.id()) != nodes_checked.end())
641 nodes_checked.insert(node.id());
647 if (local_node_id > 8 && local_node_id < 13)
659 const auto & base = elem->node_ref(local_node_id - 9);
660 const auto & apex = elem->node_ref(4);
663 CPPUNIT_ASSERT(node.relative_fuzzy_equals(base + x * (apex - base), 1e-3));
666 else if (local_node_id > 13)
676 const auto & base1 = elem->node_ref(local_node_id - 14);
677 const auto & base2 = elem->node_ref((local_node_id - 13) % 4);
678 const auto & apex = elem->node_ref(4);
680 const auto node_approx = (0.3141064847 * base1 +
681 0.3141064847 * base2 +
682 0.3717870306 * apex);
683 CPPUNIT_ASSERT(node.relative_fuzzy_equals(node_approx, 1e-3));
688 CPPUNIT_ASSERT(node_distortion_is(node,
false, 1e-2));
const std::set< Order > & elem_default_orders() const
A Node is like a Point, but with more information.
void minloc(T &r, unsigned int &min_id) const
static constexpr Real TOLERANCE
const MeshQualityInfo & get_mesh_info() const
Getter for the _system's _mesh_info attribute.
const Parallel::Communicator & comm() const
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
void sync_dofobject_data_by_id(const Communicator &comm, const Iterator &range_begin, const Iterator &range_end, SyncFunctor &sync)
Request data about a range of ghost dofobjects uniquely identified by their id.
void broadcast(T &data, const unsigned int root_id=0, const bool identical_sizes=false) const
std::string enum_to_string(const T e)
bool absolute_fuzzy_equals(const T &var1, const T2 &var2, const Real tol=TOLERANCE *TOLERANCE)
Function to check whether two variables are equal within an absolute tolerance.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
virtual const Node & node_ref(const dof_id_type i) const
bool relative_fuzzy_equals(const T &var1, const T2 &var2, const Real tol=TOLERANCE *TOLERANCE)
Function to check whether two variables are equal within a relative tolerance.
virtual void setup()
Setup method that creates equation systems, system, and constraints, to be called just prior to smoot...
A Point defines a location in LIBMESH_DIM dimensional Real space.
const Elem & get(const ElemType type_in)
virtual void smooth() override
Redefinition of the smooth function from the base class.
