1 #include <libmesh/boundary_info.h> 2 #include <libmesh/elem.h> 3 #include <libmesh/enum_elem_type.h> 4 #include <libmesh/function_base.h> 5 #include <libmesh/libmesh.h> 6 #include <libmesh/mesh_generation.h> 7 #include <libmesh/mesh_modification.h> 8 #include <libmesh/mesh_smoother_laplace.h> 9 #include <libmesh/mesh_smoother_vsmoother.h> 10 #include <libmesh/mesh_tools.h> 11 #include <libmesh/node.h> 12 #include <libmesh/reference_elem.h> 13 #include <libmesh/replicated_mesh.h> 14 #include <libmesh/system.h> 15 #include "libmesh/face_tri.h" 16 #include "libmesh/utility.h" 17 #include "libmesh/enum_to_string.h" 18 #include "libmesh/parallel_ghost_sync.h" 34 std::unique_ptr<FunctionBase<Real>> clone()
const override 36 return std::make_unique<DistortSquare>();
41 libmesh_not_implemented();
48 const Real eta = 2 * p(0) - 1;
49 const Real zeta = 2 * p(1) - 1;
50 output(0) = p(0) + (pow<3>(eta) - eta) * p(1) * (1 - p(1));
51 output(1) = p(1) + (pow<3>(zeta) - zeta) * p(0) * (1 - p(0));
61 DistortHyperCube(
const unsigned int dim) : _dim(
dim) {}
64 std::unique_ptr<FunctionBase<Real>> clone()
const override 66 return std::make_unique<DistortHyperCube>(_dim);
69 Real operator()(
const Point &,
const Real = 0.)
override { libmesh_not_implemented(); }
77 std::array<bool, 3> is_on_boundary = {
false,
false,
false};
78 for (
unsigned int i = 0; i < _dim; ++i)
80 is_on_boundary[i] =
true;
83 for (
unsigned int i = 0; i < _dim; ++i)
85 if (!is_on_boundary[i])
88 Real modulation = 0.3;
89 Real xi = 2. * p(i) - 1.;
90 for (
unsigned int j = 0; j < _dim; ++j)
94 Real pj = std::clamp(p(j), 0., 1.);
95 modulation *= (pj - 0.5) * (pj - 0.5) * 4.;
98 const auto delta = (pow<3>(xi) - xi) * modulation;
100 output(i) = (delta >
TOLERANCE) ? p(i) + delta : 1.03 * p(i);
107 const unsigned int _dim;
112 std::unique_ptr<FunctionBase<Real>> clone()
const override 114 return std::make_unique<SquareToParallelogram>();
119 libmesh_not_implemented();
128 output(0) = p(0) - 0.5 * p(1);
129 output(1) = 0.5 * std::sqrt(
Real(3)) * p(1);
136 std::unique_ptr<FunctionBase<Real>> clone()
const override 138 return std::make_unique<ParallelogramToSquare>();
143 libmesh_not_implemented();
152 output(0) = p(0) + p(1) / std::sqrt(
Real(3));
153 output(1) = (2. / std::sqrt(
Real(3))) * p(1);
171 CPPUNIT_TEST(testLaplaceQuad);
172 CPPUNIT_TEST(testLaplaceTri);
173 #if defined(LIBMESH_ENABLE_VSMOOTHER) && defined(LIBMESH_HAVE_SOLVER) 174 CPPUNIT_TEST(testVariationalEdge2);
175 CPPUNIT_TEST(testVariationalEdge3);
176 CPPUNIT_TEST(testVariationalEdge3MultipleSubdomains);
178 CPPUNIT_TEST(testVariationalQuad);
179 CPPUNIT_TEST(testVariationalQuadMultipleSubdomains);
180 CPPUNIT_TEST(testVariationalQuadTangled);
182 CPPUNIT_TEST(testVariationalTri3);
183 CPPUNIT_TEST(testVariationalTri6);
184 CPPUNIT_TEST(testVariationalTri6MultipleSubdomains);
186 CPPUNIT_TEST(testVariationalHex8);
187 CPPUNIT_TEST(testVariationalHex20);
188 CPPUNIT_TEST(testVariationalHex27);
189 CPPUNIT_TEST(testVariationalHex27Tangled);
190 CPPUNIT_TEST(testVariationalHex27MultipleSubdomains);
192 CPPUNIT_TEST(testVariationalPyramid5);
193 CPPUNIT_TEST(testVariationalPyramid13);
194 CPPUNIT_TEST(testVariationalPyramid14);
195 CPPUNIT_TEST(testVariationalPyramid18);
196 CPPUNIT_TEST(testVariationalPyramid18MultipleSubdomains);
197 #endif // LIBMESH_ENABLE_VSMOOTHER 200 CPPUNIT_TEST_SUITE_END();
211 constexpr
unsigned int n_elems_per_side = 4;
214 mesh, n_elems_per_side, n_elems_per_side, 0., 1., 0., 1., type);
221 auto center_distortion_is =
223 const Real r = node(d);
224 const Real R = r * n_elems_per_side;
232 const Real R1 = node(1) * n_elems_per_side;
240 const Real R0 = node(0) * n_elems_per_side;
250 return ((std::abs(R - std::round(R)) > distortion_tol) == distortion);
253 for (
auto node :
mesh.node_ptr_range())
255 CPPUNIT_ASSERT(center_distortion_is(*node, 0,
true));
256 CPPUNIT_ASSERT(center_distortion_is(*node, 1,
true));
262 for (
unsigned int i = 0; i != 8; ++i)
266 for (
auto node :
mesh.node_ptr_range())
268 CPPUNIT_ASSERT(center_distortion_is(*node, 0,
false, 1e-3));
269 CPPUNIT_ASSERT(center_distortion_is(*node, 1,
false, 1e-3));
276 const bool multiple_subdomains =
false,
277 const bool tangle_mesh=
false,
278 const Real tangle_damping_factor=1.0)
282 if (multiple_subdomains && tangle_mesh)
283 libmesh_not_implemented_msg(
284 "Arbitrary mesh tangling with multiple subdomains is not supported.");
288 const auto dim = ref_elem->dim();
303 mesh, n_elems_per_side, n_elems_per_side, 0., 1., 0., 1., type);
321 libmesh_error_msg(
"Unsupported dimension " <<
dim);
325 DistortHyperCube dh(
dim);
345 Real dist1_closest = std::numeric_limits<Real>::max();
347 Real dist2_closest = std::numeric_limits<Real>::max();
349 for (
const auto * node :
mesh.local_node_ptr_range())
351 const auto dist1 = (*node - p1).
norm();
352 if (dist1 < dist1_closest)
354 dist1_closest = dist1;
355 node1_id = node->id();
358 const auto dist2 = (*node - p2).
norm();
359 if (dist2 < dist2_closest)
361 dist2_closest = dist2;
362 node2_id = node->id();
367 unsigned int node1_rank;
371 unsigned int node2_rank;
383 const Point diff = node1 - node2;
387 node1 -= tangle_damping_factor * diff;
388 node2 += tangle_damping_factor * diff;
397 CPPUNIT_ASSERT(unsmoothed_info.mesh_is_tangled);
401 std::unordered_map<subdomain_id_type, Real> distorted_subdomain_volumes;
402 if (multiple_subdomains)
405 for (
auto * elem :
mesh.active_element_ptr_range())
407 unsigned int subdomain_id = 0;
409 if (elem->vertex_average()(d) > 0.5)
411 elem->subdomain_id() += subdomain_id;
417 for (
auto * elem :
mesh.active_element_ptr_range())
418 distorted_subdomain_volumes[elem->subdomain_id()] += elem->volume();
423 libmesh_error_msg_if(elem_orders.size() != 1,
424 "The variational smoother cannot be used for mixed-order meshes!");
428 const auto scale_factor = *elem_orders.begin() * (type_is_pyramid ? 2 * 4 : 1);
431 auto node_distortion_is = [&n_elems_per_side, &
dim, &boundary_info, &scale_factor](
434 std::vector<boundary_id_type> boundary_ids;
435 boundary_info.boundary_ids(&node, boundary_ids);
438 const auto num_dofs =
dim - boundary_ids.size();
458 std::size_t num_zero_or_one = 0;
460 bool distorted =
false;
463 const Real r = node(d);
464 const Real R = r * n_elems_per_side * scale_factor;
465 CPPUNIT_ASSERT_GREATER(-distortion_tol * distortion_tol, r);
466 CPPUNIT_ASSERT_GREATER(-distortion_tol * distortion_tol, 1 - r);
468 const bool d_distorted = std::abs(R - std::round(R)) > distortion_tol;
469 distorted |= d_distorted;
473 CPPUNIT_ASSERT_GREATEREQUAL(
dim - num_dofs, num_zero_or_one);
479 return distorted == distortion;
483 for (
auto node :
mesh.node_ptr_range())
484 CPPUNIT_ASSERT(node_distortion_is(*node,
true));
491 SquareToParallelogram stp;
503 ParallelogramToSquare pts;
507 if (multiple_subdomains)
512 std::unordered_map<subdomain_id_type, Real> smoothed_subdomain_volumes;
513 for (
auto * elem :
mesh.active_element_ptr_range())
514 smoothed_subdomain_volumes[elem->subdomain_id()] += elem->volume();
516 for (
const auto & pair : distorted_subdomain_volumes)
518 const auto & subdomain_id = pair.first;
521 libmesh_map_find(smoothed_subdomain_volumes, subdomain_id),
528 std::set<dof_id_type> nodes_checked;
529 for (
const auto * elem :
mesh.active_element_ptr_range())
531 for (
const auto local_node_id :
make_range(elem->n_nodes()))
533 const auto & node = elem->node_ref(local_node_id);
534 if (nodes_checked.find(node.id()) != nodes_checked.end())
537 nodes_checked.insert(node.id());
543 if (local_node_id > 8 && local_node_id < 13)
555 const auto & base = elem->node_ref(local_node_id - 9);
556 const auto & apex = elem->node_ref(4);
559 CPPUNIT_ASSERT(node.relative_fuzzy_equals(base + x * (apex - base), 1e-3));
562 else if (local_node_id > 13)
572 const auto & base1 = elem->node_ref(local_node_id - 14);
573 const auto & base2 = elem->node_ref((local_node_id - 13) % 4);
574 const auto & apex = elem->node_ref(4);
576 const auto node_approx = (0.3141064847 * base1 +
577 0.3141064847 * base2 +
578 0.3717870306 * apex);
579 CPPUNIT_ASSERT(node.relative_fuzzy_equals(node_approx, 1e-3));
584 CPPUNIT_ASSERT(node_distortion_is(node,
false, 1e-2));
596 testLaplaceSmoother(
mesh, laplace,
QUAD4);
604 testLaplaceSmoother(
mesh, laplace,
TRI3);
607 #ifdef LIBMESH_ENABLE_VSMOOTHER 613 testVariationalSmoother(
mesh, variational,
EDGE2);
621 testVariationalSmoother(
mesh, variational,
EDGE3);
629 testVariationalSmoother(
mesh, variational,
EDGE3,
true);
637 testVariationalSmoother(
mesh, variational,
QUAD4);
645 testVariationalSmoother(
mesh, variational,
QUAD4,
true);
653 testVariationalSmoother(
mesh, variational,
QUAD4,
false,
true, 0.65);
661 testVariationalSmoother(
mesh, variational,
TRI3);
669 testVariationalSmoother(
mesh, variational,
TRI6);
677 testVariationalSmoother(
mesh, variational,
TRI3,
true);
685 testVariationalSmoother(
mesh, variational,
HEX8);
693 testVariationalSmoother(
mesh, variational,
HEX20);
701 testVariationalSmoother(
mesh, variational,
HEX27);
709 testVariationalSmoother(
mesh, variational,
HEX27,
true);
717 testVariationalSmoother(
mesh, variational,
HEX27,
false,
true, 0.55);
759 #endif // LIBMESH_ENABLE_VSMOOTHER void testVariationalHex27Tangled()
static const Order type_to_default_order_map[INVALID_ELEM]
This array maps the integer representation of the ElemType enum to the default approximation order of...
ElemType
Defines an enum for geometric element types.
void testVariationalQuadTangled()
const std::set< Order > & elem_default_orders() const
The ReplicatedMesh class is derived from the MeshBase class, and is used to store identical copies of...
A Node is like a Point, but with more information.
virtual void zero() override final
Set every element in the vector to 0.
void testVariationalTri3()
void testLaplaceSmoother(ReplicatedMesh &mesh, MeshSmoother &smoother, ElemType type)
void minloc(T &r, unsigned int &min_id) const
libMesh::Parallel::Communicator * TestCommWorld
static constexpr Real TOLERANCE
void resize(const unsigned int n)
Resize the vector.
void testVariationalQuad()
const MeshQualityInfo & get_mesh_info() const
Getter for the _system's _mesh_info attribute.
void testVariationalEdge2()
void testVariationalPyramid18MultipleSubdomains()
const Parallel::Communicator & comm() const
void testVariationalSmoother(ReplicatedMesh &mesh, VariationalMeshSmoother &smoother, const ElemType type, const bool multiple_subdomains=false, const bool tangle_mesh=false, const Real tangle_damping_factor=1.0)
The libMesh namespace provides an interface to certain functionality in the library.
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
virtual void smooth()=0
Function which actually performs the smoothing operations.
void testVariationalHex27()
This class defines the data structures necessary for Laplace smoothing.
void testVariationalEdge3()
void testVariationalPyramid5()
void testVariationalHex8()
static const dof_id_type invalid_id
An invalid id to distinguish an uninitialized DofObject.
CPPUNIT_TEST_SUITE_REGISTRATION(MeshSmootherTest)
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)
void testVariationalHex20()
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.
void testVariationalHex27MultipleSubdomains()
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void testVariationalPyramid14()
This is an implementation of Larisa Branets' smoothing algorithms.
void testVariationalQuadMultipleSubdomains()
void testVariationalTri6MultipleSubdomains()
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...
void testVariationalPyramid18()
void testVariationalTri6()
void testVariationalPyramid13()
virtual const Node & node_ref(const dof_id_type i) const
This class provides the necessary interface for mesh smoothing.
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.
Defines a dense vector for use in Finite Element-type computations.
virtual void setup()
Setup method that creates equation systems, system, and constraints, to be called just prior to smoot...
Base class for functors that can be evaluated at a point and (optionally) time.
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.
void testVariationalEdge3MultipleSubdomains()