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/cell_prism.h" 17 #include "libmesh/utility.h" 18 #include "libmesh/enum_to_string.h" 19 #include "libmesh/parallel_ghost_sync.h" 35 std::unique_ptr<FunctionBase<Real>> clone()
const override 37 return std::make_unique<DistortSquare>();
42 libmesh_not_implemented();
49 const Real eta = 2 * p(0) - 1;
50 const Real zeta = 2 * p(1) - 1;
51 output(0) = p(0) + (pow<3>(eta) - eta) * p(1) * (1 - p(1));
52 output(1) = p(1) + (pow<3>(zeta) - zeta) * p(0) * (1 - p(0));
62 DistortHyperCube(
const unsigned int dim) : _dim(
dim) {}
65 std::unique_ptr<FunctionBase<Real>> clone()
const override 67 return std::make_unique<DistortHyperCube>(_dim);
70 Real operator()(
const Point &,
const Real = 0.)
override { libmesh_not_implemented(); }
78 std::array<bool, 3> is_on_boundary = {
false,
false,
false};
79 for (
unsigned int i = 0; i < _dim; ++i)
81 is_on_boundary[i] =
true;
84 for (
unsigned int i = 0; i < _dim; ++i)
86 if (!is_on_boundary[i])
89 Real modulation = 0.3;
90 Real xi = 2. * p(i) - 1.;
91 for (
unsigned int j = 0; j < _dim; ++j)
95 Real pj = std::clamp(p(j), 0., 1.);
96 modulation *= (pj - 0.5) * (pj - 0.5) * 4.;
99 const auto delta = (pow<3>(xi) - xi) * modulation;
101 output(i) = (delta >
TOLERANCE) ? p(i) + delta : 1.03 * p(i);
108 const unsigned int _dim;
113 std::unique_ptr<FunctionBase<Real>> clone()
const override 115 return std::make_unique<SquareToParallelogram>();
120 libmesh_not_implemented();
129 output(0) = p(0) - 0.5 * p(1);
130 output(1) = 0.5 * std::sqrt(
Real(3)) * p(1);
137 std::unique_ptr<FunctionBase<Real>> clone()
const override 139 return std::make_unique<ParallelogramToSquare>();
144 libmesh_not_implemented();
153 output(0) = p(0) + p(1) / std::sqrt(
Real(3));
154 output(1) = (2. / std::sqrt(
Real(3))) * p(1);
161 std::unique_ptr<FunctionBase<Real>> clone()
const override 163 return std::make_unique<CubeToParallelepiped>();
168 libmesh_not_implemented();
180 output(0) = p(0) + 0.5 * p(1);
181 output(1) = 0.5 * std::sqrt(3) * p(1);
184 output(2) = p(2) * 0.25 * std::sqrt(3.);
190 std::unique_ptr<FunctionBase<Real>> clone()
const override 192 return std::make_unique<ParallelepipedToCube>();
197 libmesh_not_implemented();
209 output(0) = p(0) - p(1) / std::sqrt(3.);
210 output(1) = (2. / std::sqrt(3)) * p(1);
213 output(2) = p(2) * 4. / std::sqrt(3.);
230 CPPUNIT_TEST(testLaplaceQuad);
231 CPPUNIT_TEST(testLaplaceTri);
232 #if defined(LIBMESH_ENABLE_VSMOOTHER) && defined(LIBMESH_HAVE_SOLVER) 233 CPPUNIT_TEST(testVariationalEdge2);
234 CPPUNIT_TEST(testVariationalEdge3);
235 CPPUNIT_TEST(testVariationalEdge3MultipleSubdomains);
237 CPPUNIT_TEST(testVariationalQuad4);
238 CPPUNIT_TEST(testVariationalQuad4MultipleSubdomains);
239 CPPUNIT_TEST(testVariationalQuad8);
240 CPPUNIT_TEST(testVariationalQuad9);
241 CPPUNIT_TEST(testVariationalQuad9MultipleSubdomains);
242 CPPUNIT_TEST(testVariationalQuad4Tangled);
244 CPPUNIT_TEST(testVariationalTri3);
245 CPPUNIT_TEST(testVariationalTri6);
246 CPPUNIT_TEST(testVariationalTri6MultipleSubdomains);
248 CPPUNIT_TEST(testVariationalHex8);
249 CPPUNIT_TEST(testVariationalHex20);
250 CPPUNIT_TEST(testVariationalHex27);
251 CPPUNIT_TEST(testVariationalHex27Tangled);
252 CPPUNIT_TEST(testVariationalHex27MultipleSubdomains);
254 CPPUNIT_TEST(testVariationalPyramid5);
255 CPPUNIT_TEST(testVariationalPyramid13);
256 CPPUNIT_TEST(testVariationalPyramid14);
257 CPPUNIT_TEST(testVariationalPyramid18);
258 CPPUNIT_TEST(testVariationalPyramid18MultipleSubdomains);
260 CPPUNIT_TEST(testVariationalPrism6);
261 CPPUNIT_TEST(testVariationalPrism15);
262 CPPUNIT_TEST(testVariationalPrism18);
263 CPPUNIT_TEST(testVariationalPrism20);
264 CPPUNIT_TEST(testVariationalPrism21);
265 CPPUNIT_TEST(testVariationalPrism21MultipleSubdomains);
267 #if defined(LIBMESH_HAVE_GZSTREAM) 268 CPPUNIT_TEST(testVariationalMixed2D);
269 CPPUNIT_TEST(testVariationalMixed3D);
272 #endif // LIBMESH_ENABLE_VSMOOTHER 275 CPPUNIT_TEST_SUITE_END();
286 constexpr
unsigned int n_elems_per_side = 4;
289 mesh, n_elems_per_side, n_elems_per_side, 0., 1., 0., 1., type);
296 auto center_distortion_is =
298 const Real r = node(d);
299 const Real R = r * n_elems_per_side;
307 const Real R1 = node(1) * n_elems_per_side;
315 const Real R0 = node(0) * n_elems_per_side;
325 return ((std::abs(R - std::round(R)) > distortion_tol) == distortion);
328 for (
auto node :
mesh.node_ptr_range())
330 CPPUNIT_ASSERT(center_distortion_is(*node, 0,
true));
331 CPPUNIT_ASSERT(center_distortion_is(*node, 1,
true));
337 for (
unsigned int i = 0; i != 8; ++i)
341 for (
auto node :
mesh.node_ptr_range())
343 CPPUNIT_ASSERT(center_distortion_is(*node, 0,
false, 1e-3));
344 CPPUNIT_ASSERT(center_distortion_is(*node, 1,
false, 1e-3));
351 const bool multiple_subdomains =
false,
352 const bool tangle_mesh=
false,
353 const Real tangle_damping_factor=1.0)
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();
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);
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;
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));
709 const bool distortion,
715 const auto & gold_node = gold_mesh.
node_ref(node_id);
719 const bool d_distorted = std::abs(gold_node(d) - node(d)) > distortion_tol;
731 CPPUNIT_ASSERT(mesh_distortion_is(gold_mesh,
mesh,
true));
742 CPPUNIT_ASSERT(mesh_distortion_is(gold_mesh,
mesh,
false));
750 testLaplaceSmoother(
mesh, laplace,
QUAD4);
758 testLaplaceSmoother(
mesh, laplace,
TRI3);
761 #ifdef LIBMESH_ENABLE_VSMOOTHER 767 testVariationalSmoother(
mesh, variational,
EDGE2);
775 testVariationalSmoother(
mesh, variational,
EDGE3);
783 testVariationalSmoother(
mesh, variational,
EDGE3,
true);
791 testVariationalSmoother(
mesh, variational,
QUAD4);
799 testVariationalSmoother(
mesh, variational,
QUAD4,
true);
807 testVariationalSmoother(
mesh, variational,
QUAD8);
815 testVariationalSmoother(
mesh, variational,
QUAD9);
823 testVariationalSmoother(
mesh, variational,
QUAD9,
true);
831 testVariationalSmoother(
mesh, variational,
QUAD4,
false,
true, 0.65);
839 testVariationalSmoother(
mesh, variational,
TRI3);
847 testVariationalSmoother(
mesh, variational,
TRI6);
855 testVariationalSmoother(
mesh, variational,
TRI3,
true);
863 testVariationalSmoother(
mesh, variational,
HEX8);
871 testVariationalSmoother(
mesh, variational,
HEX20);
879 testVariationalSmoother(
mesh, variational,
HEX27);
887 testVariationalSmoother(
mesh, variational,
HEX27,
true);
895 testVariationalSmoother(
mesh, variational,
HEX27,
false,
true, 0.55);
943 testVariationalSmoother(
mesh, variational,
PRISM6);
951 testVariationalSmoother(
mesh, variational,
PRISM15);
960 testVariationalSmoother(
mesh, variational,
PRISM18);
969 testVariationalSmoother(
mesh, variational,
PRISM20);
978 testVariationalSmoother(
mesh, variational,
PRISM21);
987 testVariationalSmoother(
mesh, variational,
PRISM21,
true);
993 mesh.
read(
"meshes/quad4_tri3_smoothed.xda.gz");
995 testVariationalSmootherRegression(
mesh);
1001 mesh.
read(
"meshes/hex8_prism6_smoothed.xda.gz");
1003 testVariationalSmootherRegression(
mesh);
1005 #endif // LIBMESH_ENABLE_VSMOOTHER virtual dof_id_type n_nodes() const override final
void testVariationalHex27Tangled()
void testVariationalQuad8()
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.
const std::set< Order > & elem_default_orders() const
virtual void read(const std::string &name, void *mesh_data=nullptr, bool skip_renumber_nodes_and_elements=false, bool skip_find_neighbors=false)=0
Interfaces for reading/writing a mesh to/from a file.
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 testVariationalPrism18()
void allow_renumbering(bool allow)
If false is passed in then this mesh will no longer be renumbered when being prepared for use...
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 testVariationalQuad4()
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)
void testVariationalPrism21MultipleSubdomains()
void testVariationalMixed2D()
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()
static constexpr dof_id_type invalid_id
An invalid id to distinguish an uninitialized DofObject.
void testVariationalPrism21()
void testVariationalQuad4Tangled()
void testVariationalPyramid5()
void testVariationalHex8()
void testVariationalQuad9()
void testVariationalQuad9MultipleSubdomains()
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
void testVariationalPrism6()
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 testVariationalPrism20()
void testVariationalSmootherRegression(const ReplicatedMesh &gold_mesh)
void testVariationalPyramid14()
This is an implementation of Larisa Branets' smoothing algorithms.
void testVariationalPrism15()
void testVariationalTri6MultipleSubdomains()
void testVariationalMixed3D()
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...
unsigned int mesh_dimension() const
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...
void testVariationalQuad4MultipleSubdomains()
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 dof_id_type n_nodes() const =0
virtual void smooth() override
Redefinition of the smooth function from the base class.
void testVariationalEdge3MultipleSubdomains()