libMesh
mesh_smoother_test.C
Go to the documentation of this file.
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> // LIBMESH_HAVE_SOLVER define
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"
20 
21 #include "test_comm.h"
22 #include "libmesh_cppunit.h"
23 
24 #include <random>
25 
26 namespace
27 {
28 using namespace libMesh;
29 using Utility::pow;
30 
31 // Distortion function that doesn't distort boundary nodes
32 // 2D only, use for LaplaceMeshSmoother
33 class DistortSquare : public FunctionBase<Real>
34 {
35  std::unique_ptr<FunctionBase<Real>> clone() const override
36  {
37  return std::make_unique<DistortSquare>();
38  }
39 
40  Real operator()(const Point &, const Real = 0.) override
41  {
42  libmesh_not_implemented();
43  } // scalar-only API
44 
45  // Skew inward based on a cubic function
46  void operator()(const Point & p, const Real, DenseVector<Real> & output)
47  {
48  output.resize(3);
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));
53  output(2) = 0;
54  }
55 };
56 
57 // Distortion function supporting 1D, 2D, and 3D
58 // Use for VariationalMeshSmoother
59 class DistortHyperCube : public FunctionBase<Real>
60 {
61 public:
62  DistortHyperCube(const unsigned int dim) : _dim(dim) {}
63 
64 private:
65  std::unique_ptr<FunctionBase<Real>> clone() const override
66  {
67  return std::make_unique<DistortHyperCube>(_dim);
68  }
69 
70  Real operator()(const Point &, const Real = 0.) override { libmesh_not_implemented(); }
71 
72  void operator()(const Point & p, const Real, DenseVector<Real> & output) override
73  {
74  output.resize(3);
75  output.zero();
76 
77  // Count how many coordinates are exactly on the boundary (0 or 1)
78  std::array<bool, 3> is_on_boundary = {false, false, false};
79  for (unsigned int i = 0; i < _dim; ++i)
80  if (std::abs(p(i)) < TOLERANCE || std::abs(p(i) - 1.) < TOLERANCE)
81  is_on_boundary[i] = true;
82 
83  // Distort only those directions not fixed on the boundary
84  for (unsigned int i = 0; i < _dim; ++i)
85  {
86  if (!is_on_boundary[i]) // only distort free dimensions
87  {
88  // This value controls the strength of the distortion
89  Real modulation = 0.3;
90  Real xi = 2. * p(i) - 1.;
91  for (unsigned int j = 0; j < _dim; ++j)
92  {
93  if (j != i)
94  {
95  Real pj = std::clamp(p(j), 0., 1.); // ensure numeric safety
96  modulation *= (pj - 0.5) * (pj - 0.5) * 4.; // quadratic bump centered at 0.5
97  }
98  }
99  const auto delta = (pow<3>(xi) - xi) * modulation;
100  // Check for delta = 0 to make sure we perturb every point
101  output(i) = (delta > TOLERANCE) ? p(i) + delta : 1.03 * p(i);
102  }
103  else
104  output(i) = p(i); // dimension on boundary remains unchanged
105  }
106  }
107 
108  const unsigned int _dim;
109 };
110 
111 class SquareToParallelogram : public FunctionBase<Real>
112 {
113  std::unique_ptr<FunctionBase<Real>> clone() const override
114  {
115  return std::make_unique<SquareToParallelogram>();
116  }
117 
118  Real operator()(const Point &, const Real = 0.) override
119  {
120  libmesh_not_implemented();
121  } // scalar-only API
122 
123  // Has the effect that a square, meshed into right triangles with diagonals
124  // rising from lower-left to upper-right, is transformed into a left-leaning
125  // parallelogram of eqilateral triangles.
126  void operator()(const Point & p, const Real, DenseVector<Real> & output)
127  {
128  output.resize(3);
129  output(0) = p(0) - 0.5 * p(1);
130  output(1) = 0.5 * std::sqrt(Real(3)) * p(1);
131  output(2) = p(2);
132  }
133 };
134 
135 class ParallelogramToSquare : public FunctionBase<Real>
136 {
137  std::unique_ptr<FunctionBase<Real>> clone() const override
138  {
139  return std::make_unique<ParallelogramToSquare>();
140  }
141 
142  Real operator()(const Point &, const Real = 0.) override
143  {
144  libmesh_not_implemented();
145  } // scalar-only API
146 
147  // Has the effect that a left-leaning parallelogram of equilateral triangles is transformed
148  // into a square of right triangle with diagonals rising from lower-left to upper-right.
149  // This is the inversion of the SquareToParallelogram mapping.
150  void operator()(const Point & p, const Real, DenseVector<Real> & output)
151  {
152  output.resize(3);
153  output(0) = p(0) + p(1) / std::sqrt(Real(3));
154  output(1) = (2. / std::sqrt(Real(3))) * p(1);
155  output(2) = p(2);
156  }
157 };
158 
159 class CubeToParallelepiped : public FunctionBase<Real>
160 {
161  std::unique_ptr<FunctionBase<Real>> clone() const override
162  {
163  return std::make_unique<CubeToParallelepiped>();
164  }
165 
166  Real operator()(const Point &, const Real = 0.) override
167  {
168  libmesh_not_implemented();
169  } // scalar-only API
170 
171  // Has the effect that a cube, meshed into right triangular prisms with diagonals
172  // rising from lower-right to upper-left, is transformed into a right-leaning
173  // parallelepiped of equilateral triangular prisms. Additionally, the z direction
174  // is scaled to ensure element height to base aspect ratios match the target element.
175  // Without the correct aspect ratio, the smoothed mesh nodes are off and the
176  // distortion_is assertions fail.
177  void operator()(const Point & p, const Real, DenseVector<Real> & output)
178  {
179  output.resize(3);
180  output(0) = p(0) + 0.5 * p(1);
181  output(1) = 0.5 * std::sqrt(3) * p(1);
182  // Adjusting z to get the triangular area to element height aspect ratio To
183  // match the target element
184  output(2) = p(2) * 0.25 * std::sqrt(3.);
185  }
186 };
187 
188 class ParallelepipedToCube : public FunctionBase<Real>
189 {
190  std::unique_ptr<FunctionBase<Real>> clone() const override
191  {
192  return std::make_unique<ParallelepipedToCube>();
193  }
194 
195  Real operator()(const Point &, const Real = 0.) override
196  {
197  libmesh_not_implemented();
198  } // scalar-only API
199 
200  // Has the effect that a right-leaning parallelepiped of equilaterial triangular
201  // prisms with is transformed into a cube of right triangular prisms with
202  // diagonals rising from lower-right to upper-left. Additionally, the z direction
203  // is scaled to ensure element heights align with the values expected by the
204  // distortion_is function.
205  // This is the inversion of the CubeToParallelepiped mapping.
206  void operator()(const Point & p, const Real, DenseVector<Real> & output)
207  {
208  output.resize(3);
209  output(0) = p(0) - p(1) / std::sqrt(3.);
210  output(1) = (2. / std::sqrt(3)) * p(1);
211  // Undoing the aspect ratio adjustment to get the z divisions to work with
212  // distortion_is
213  output(2) = p(2) * 4. / std::sqrt(3.);
214  }
215 };
216 }
217 
218 using namespace libMesh;
219 
220 class MeshSmootherTest : public CppUnit::TestCase
221 {
226 public:
227  LIBMESH_CPPUNIT_TEST_SUITE(MeshSmootherTest);
228 
229 #if LIBMESH_DIM > 2
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);
236 
237  CPPUNIT_TEST(testVariationalQuad4);
238  CPPUNIT_TEST(testVariationalQuad4MultipleSubdomains);
239  CPPUNIT_TEST(testVariationalQuad8);
240  CPPUNIT_TEST(testVariationalQuad9);
241  CPPUNIT_TEST(testVariationalQuad9MultipleSubdomains);
242  CPPUNIT_TEST(testVariationalQuad4Tangled);
243 
244  CPPUNIT_TEST(testVariationalTri3);
245  CPPUNIT_TEST(testVariationalTri6);
246  CPPUNIT_TEST(testVariationalTri6MultipleSubdomains);
247 
248  CPPUNIT_TEST(testVariationalHex8);
249  CPPUNIT_TEST(testVariationalHex20);
250  CPPUNIT_TEST(testVariationalHex27);
251  CPPUNIT_TEST(testVariationalHex27Tangled);
252  CPPUNIT_TEST(testVariationalHex27MultipleSubdomains);
253 
254  CPPUNIT_TEST(testVariationalPyramid5);
255  CPPUNIT_TEST(testVariationalPyramid13);
256  CPPUNIT_TEST(testVariationalPyramid14);
257  CPPUNIT_TEST(testVariationalPyramid18);
258  CPPUNIT_TEST(testVariationalPyramid18MultipleSubdomains);
259 
260  CPPUNIT_TEST(testVariationalPrism6);
261  CPPUNIT_TEST(testVariationalPrism15);
262  CPPUNIT_TEST(testVariationalPrism18);
263  CPPUNIT_TEST(testVariationalPrism20);
264  CPPUNIT_TEST(testVariationalPrism21);
265  CPPUNIT_TEST(testVariationalPrism21MultipleSubdomains);
266 
267 #if defined(LIBMESH_HAVE_GZSTREAM)
268  CPPUNIT_TEST(testVariationalMixed2D);
269  CPPUNIT_TEST(testVariationalMixed3D);
270 #endif
271 
272 #endif // LIBMESH_ENABLE_VSMOOTHER
273 #endif
274 
275  CPPUNIT_TEST_SUITE_END();
276 
277 public:
278  void setUp() {}
279 
280  void tearDown() {}
281 
283  {
284  LOG_UNIT_TEST;
285 
286  constexpr unsigned int n_elems_per_side = 4;
287 
289  mesh, n_elems_per_side, n_elems_per_side, 0., 1., 0., 1., type);
290 
291  // Move it around so we have something that needs smoothing
292  DistortSquare ds;
294 
295  // Assert the distortion is as expected
296  auto center_distortion_is =
297  [](const Node & node, int d, bool distortion, Real distortion_tol = TOLERANCE) {
298  const Real r = node(d);
299  const Real R = r * n_elems_per_side;
300  CPPUNIT_ASSERT_GREATER(-TOLERANCE * TOLERANCE, r);
301  CPPUNIT_ASSERT_GREATER(-TOLERANCE * TOLERANCE, 1 - r);
302 
303  // If we're at the boundaries we should *never* be distorted
304  if (std::abs(node(0)) < TOLERANCE * TOLERANCE ||
305  std::abs(node(0) - 1) < TOLERANCE * TOLERANCE)
306  {
307  const Real R1 = node(1) * n_elems_per_side;
308  CPPUNIT_ASSERT_LESS(TOLERANCE * TOLERANCE, std::abs(R1 - std::round(R1)));
309  return true;
310  }
311 
312  if (std::abs(node(1)) < TOLERANCE * TOLERANCE ||
313  std::abs(node(1) - 1) < TOLERANCE * TOLERANCE)
314  {
315  const Real R0 = node(0) * n_elems_per_side;
316  CPPUNIT_ASSERT_LESS(TOLERANCE * TOLERANCE, std::abs(R0 - std::round(R0)));
317 
318  return true;
319  }
320 
321  // If we're at the center we're fine
322  if (std::abs(r - 0.5) < TOLERANCE * TOLERANCE)
323  return true;
324 
325  return ((std::abs(R - std::round(R)) > distortion_tol) == distortion);
326  };
327 
328  for (auto node : mesh.node_ptr_range())
329  {
330  CPPUNIT_ASSERT(center_distortion_is(*node, 0, true));
331  CPPUNIT_ASSERT(center_distortion_is(*node, 1, true));
332  }
333 
334  // Enough iterations to mostly fix us up. Laplace seems to be at 1e-3
335  // tolerance by iteration 6, so hopefully everything is there on any
336  // system by 8.
337  for (unsigned int i = 0; i != 8; ++i)
338  smoother.smooth();
339 
340  // Make sure we're not too distorted anymore.
341  for (auto node : mesh.node_ptr_range())
342  {
343  CPPUNIT_ASSERT(center_distortion_is(*node, 0, false, 1e-3));
344  CPPUNIT_ASSERT(center_distortion_is(*node, 1, false, 1e-3));
345  }
346  }
347 
349  VariationalMeshSmoother & smoother,
350  const ElemType type,
351  const bool multiple_subdomains = false,
352  const bool tangle_mesh=false,
353  const Real tangle_damping_factor=1.0)
354  {
355  LOG_UNIT_TEST;
356 
357  if (multiple_subdomains && tangle_mesh)
358  libmesh_not_implemented_msg(
359  "Arbitrary mesh tangling with multiple subdomains is not supported.");
360 
361  // Get mesh dimension, determine whether type is triangular
362  const auto * ref_elem = &(ReferenceElem::get(type));
363  const auto dim = ref_elem->dim();
364  const bool type_is_tri = Utility::enum_to_string(type).compare(0, 3, "TRI") == 0;
365  const bool type_is_pyramid = Utility::enum_to_string(type).compare(0, 7, "PYRAMID") == 0;
366  const bool type_is_prism = Utility::enum_to_string(type).compare(0, 5, "PRISM") == 0;
367 
368  // Used fewer elems for higher order types, as extra midpoint nodes will add
369  // enough complexity
370  unsigned int n_elems_per_side = 4 / Elem::type_to_default_order_map[type];
371 
372  switch (dim)
373  {
374  case 1:
375  MeshTools::Generation::build_line(mesh, n_elems_per_side, 0., 1., type);
376  break;
377  case 2:
379  mesh, n_elems_per_side, n_elems_per_side, 0., 1., 0., 1., type);
380  break;
381 
382  case 3:
384  n_elems_per_side,
385  n_elems_per_side,
386  n_elems_per_side,
387  0.,
388  1.,
389  0.,
390  1.,
391  0.,
392  1.,
393  type);
394  break;
395 
396  default:
397  libmesh_error_msg("Unsupported dimension " << dim);
398  }
399 
400  // Move it around so we have something that needs smoothing
401  DistortHyperCube dh(dim);
403 
404  const auto & boundary_info = mesh.get_boundary_info();
405 
406  if (tangle_mesh)
407  {
408  // We tangle the mesh by (partially) swapping the locations of 2 nodes.
409  // If the n-dimentional hypercube is represented as a grid with
410  // (undistorted) divisions occuring at integer numbers, we will
411  // (partially) swap the nodes nearest p1 = (1,1,1) and p2 = (2,1,2).
412 
413  // Define p1 and p2 given the integer indices
414  // Undistorted elem side length
415  const Real dr = 1. / (n_elems_per_side * Elem::type_to_default_order_map[type]);
416  const Point p1 = Point(dr, dim > 1 ? dr : 0, dim > 2 ? dr : 0);
417  const Point p2 = Point(2. * dr, dim > 1 ? dr : 0, dim > 2 ? 2. * dr : 0);
418 
419  // ids and distances of mesh nodes nearest p1 and p2
421  Real dist1_closest = std::numeric_limits<Real>::max();
423  Real dist2_closest = std::numeric_limits<Real>::max();
424 
425  for (const auto * node : mesh.local_node_ptr_range())
426  {
427  const auto dist1 = (*node - p1).norm();
428  if (dist1 < dist1_closest)
429  {
430  dist1_closest = dist1;
431  node1_id = node->id();
432  }
433 
434  const auto dist2 = (*node - p2).norm();
435  if (dist2 < dist2_closest)
436  {
437  dist2_closest = dist2;
438  node2_id = node->id();
439  }
440  }
441 
442  // parallel communication
443  unsigned int node1_rank;
444  mesh.comm().minloc(dist1_closest, node1_rank);
445  mesh.comm().broadcast(node1_id, node1_rank);
446 
447  unsigned int node2_rank;
448  mesh.comm().minloc(dist2_closest, node2_rank);
449  mesh.comm().broadcast(node2_id, node2_rank);
450 
451  // modify the nodes
452  if ((node1_id != DofObject::invalid_id) && (node2_id != DofObject::invalid_id))
453  {
454  libmesh_assert(node1_id != node2_id);
455 
456  auto & node1 = mesh.node_ref(node1_id);
457  auto & node2 = mesh.node_ref(node2_id);
458 
459  const Point diff = node1 - node2;
460  // Note that a damping factor of 1 will swap the nodes and one
461  // of 0.5 will move the nodes to the mean of their original
462  // locations.
463  node1 -= tangle_damping_factor * diff;
464  node2 += tangle_damping_factor * diff;
465  }
466 
467  SyncNodalPositions sync_object(mesh);
468  Parallel::sync_dofobject_data_by_id (mesh.comm(), mesh.nodes_begin(), mesh.nodes_end(), sync_object);
469 
470  // Make sure the mesh is tangled
471  smoother.setup(); // Need this to create the system we are about to access
472  const auto & unsmoothed_info = smoother.get_mesh_info();
473  CPPUNIT_ASSERT(unsmoothed_info.mesh_is_tangled);
474  }
475 
476  // Add multiple subdomains if requested
477  std::unordered_map<subdomain_id_type, Real> distorted_subdomain_volumes;
478  if (multiple_subdomains)
479  {
480  // Modify the subdomain ids in an interesting way
481  for (auto * elem : mesh.active_element_ptr_range())
482  {
483  unsigned int subdomain_id = 0;
484  for (const auto d : make_range(dim))
485  if (elem->vertex_average()(d) > 0.5)
486  ++subdomain_id;
487  elem->subdomain_id() += subdomain_id;
488  }
489 
490  // This loop should NOT be combined with the one above because we need to
491  // finish checking and updating subdomain ids for all elements before
492  // recording the final subdomain volumes.
493  for (auto * elem : mesh.active_element_ptr_range())
494  distorted_subdomain_volumes[elem->subdomain_id()] += elem->volume();
495  }
496 
497  // Get the mesh order
498  const auto & elem_orders = mesh.elem_default_orders();
499  libmesh_error_msg_if(elem_orders.size() != 1,
500  "The variational smoother cannot be used for mixed-order meshes!");
501  // For pyramids, the factor 2 accounts for the account that cubes of side
502  // length s are divided into pyramids of base side length s and height s/2.
503  // The factor 4 lets us catch triangular face midpoint nodes for PYRAMID18 elements.
504  const auto scale_factor = *elem_orders.begin() * (type_is_pyramid ? 2 * 4 : 1);
505 
506  // Function to assert the node distortion is as expected
507  auto node_distortion_is = [&n_elems_per_side, &dim, &boundary_info, &scale_factor, &type_is_prism](
508  const Node & node, bool distortion, Real distortion_tol = TOLERANCE) {
509  // Get boundary ids associated with the node
510  std::vector<boundary_id_type> boundary_ids;
511  boundary_info.boundary_ids(&node, boundary_ids);
512 
513  // This tells us what type of node we are: internal, sliding, or fixed
514  const auto num_dofs = dim - boundary_ids.size();
515  /*
516  * The following cases of num_dofs are possible, ASSUMING all boundaries
517  * are non-overlapping
518  * 3D: 3-0, 3-1, 3-2, 3-3
519  * = 3 2 1 0
520  * internal sliding, sliding, fixed
521  * 2D: 2-0, 2-1, 2-2
522  * = 2 1 0
523  * internal sliding, fixed
524  * 1D: 1-0, 1-1
525  * = 1 0
526  * internal fixed
527  *
528  * We expect that R is an integer in [0, n_elems_per_side] for
529  * num_dofs of the node's cooridinantes, while the remaining coordinates
530  * are fixed to the boundary with value 0 or 1. In other words, at LEAST
531  * dim - num_dofs coordinantes should be 0 or 1.
532  */
533 
534  std::size_t num_zero_or_one = 0;
535 
536  bool distorted = false;
537  for (const auto d : make_range(dim))
538  {
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);
543 
544  bool d_distorted = std::abs(R - std::round(R)) > distortion_tol;
545  if (type_is_prism && (scale_factor == 3))
546  {
547  // Adjustment required for triangular face nodes of PRISM20/21 elements.
548  // These elements have fe_order 3. This takes care of nodes occuring at
549  // thirds, but not halves.
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;
554  }
555  distorted |= d_distorted;
556  num_zero_or_one += (absolute_fuzzy_equals(r, 0.) || absolute_fuzzy_equals(r, 1.));
557  }
558 
559  CPPUNIT_ASSERT_GREATEREQUAL(dim - num_dofs, num_zero_or_one);
560 
561  // We can never expect a fixed node to be distorted
562  if (num_dofs == 0)
563  // if (num_dofs < dim)
564  return true;
565  return distorted == distortion;
566  };
567 
568  // Make sure our DistortHyperCube transformation has distorted the mesh
569  for (auto node : mesh.node_ptr_range())
570  CPPUNIT_ASSERT(node_distortion_is(*node, true));
571 
572  // Transform the square mesh of triangles to a parallelogram mesh of
573  // triangles. This will allow the Variational Smoother to smooth the mesh
574  // to the optimal case of equilateral triangles
575  if (type_is_tri)
576  {
577  // Transform the square mesh of triangles to a parallelogram mesh of
578  // triangles. This will allow the Variational Smoother to smooth the mesh
579  // to the optimal case of equilateral triangles
580  SquareToParallelogram stp;
582  }
583  else if (type_is_prism)
584  {
585  // Transform the cube mesh of prisms to a parallelepiped mesh of
586  // prisms. This will allow the Variational Smoother to smooth the mesh
587  // to the optimal case of equilateral triangular prisms
588  CubeToParallelepiped ctp;
590  }
591 
592  smoother.smooth();
593 
594  if (type_is_tri)
595  {
596  // Transform the parallelogram mesh back to a square mesh. In the case of
597  // the Variational Smoother, equilateral triangules will be transformed
598  // into right triangules that align with the original undistorted mesh.
599  ParallelogramToSquare pts;
601  }
602  else if (type_is_prism)
603  {
604  // Transform the parallelepiped mesh back to a cube mesh. In the case of
605  // the Variational Smoother, equilateral triangular prisms will be transformed
606  // into right triangular prisms that align with the original undistorted mesh.
607  ParallelepipedToCube ptc;
609  }
610 
611  if (multiple_subdomains)
612  {
613  // Make sure the subdomain volumes didn't change. Although this doesn't
614  // guarantee that the subdomain didn't change, it is a good indicator
615  // and is friedly to sliding subdomain boundary nodes.
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();
619 
620  for (const auto & pair : distorted_subdomain_volumes)
621  {
622  const auto & subdomain_id = pair.first;
623  CPPUNIT_ASSERT(
624  relative_fuzzy_equals(libmesh_map_find(distorted_subdomain_volumes, subdomain_id),
625  libmesh_map_find(smoothed_subdomain_volumes, subdomain_id),
626  TOLERANCE));
627  }
628  }
629  else
630  {
631  // Make sure we're not too distorted anymore
632  std::set<dof_id_type> nodes_checked;
633  for (const auto * elem : mesh.active_element_ptr_range())
634  {
635  for (const auto local_node_id : make_range(elem->n_nodes()))
636  {
637  const auto & node = elem->node_ref(local_node_id);
638  if (nodes_checked.find(node.id()) != nodes_checked.end())
639  continue;
640 
641  nodes_checked.insert(node.id());
642 
643  // Check for special case where pyramidal base-to-apex midpoint
644  // nodes are not smoothed to the actual midpoints.
645  if (type_is_pyramid)
646  {
647  if (local_node_id > 8 && local_node_id < 13)
648  {
649  // Base-Apex midpoint nodes
650  //
651  // Due to the nature of the pyramidal target element and the
652  // cubic nature of the test mesh, for a dilation weight o
653  // 0.5, smoothed base-apex midpoint nodes are smoothed to
654  // the point base + x * (apex - base) instead of
655  // base + 0.5 (apex - base), where x is in (0,1) and
656  // depends on the number of nodes in the element.
657  // We hard-code this check below.
658 
659  const auto & base = elem->node_ref(local_node_id - 9);
660  const auto & apex = elem->node_ref(4);
661  const Real x = (type == PYRAMID18) ? 0.569332 : 0.549876;
662 
663  CPPUNIT_ASSERT(node.relative_fuzzy_equals(base + x * (apex - base), 1e-3));
664  continue;
665  }
666  else if (local_node_id > 13)
667  {
668  // Triangular face midpoint nodes
669  //
670  // Due to the nature of the pyramidal target element and the
671  // cubic nature of the test mesh, for a dilation weight o
672  // 0.5, smoothed triangular face midpoint nodes are
673  // smoothed a weighted average of the constituent
674  // vertices instead of an unweighted average.
675  // We hard-code this check below.
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);
679 
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));
684  continue;
685  }
686  }
687 
688  CPPUNIT_ASSERT(node_distortion_is(node, false, 1e-2));
689 
690  }
691  }
692  }
693  }
694 
695  // Function that distorts and smooths a mesh and checks that the result is
696  // equal to the original "gold" mesh
698  {
699  // Make a copy of the gold mesh to distort and then smooth
700  ReplicatedMesh mesh(gold_mesh);
701 
702  // Move it around so we have something that needs smoothing
703  DistortHyperCube dh(mesh.mesh_dimension());
705 
706  // Function to assert the distortion is as expected
707  auto mesh_distortion_is = [](const ReplicatedMesh & gold_mesh,
708  const ReplicatedMesh & mesh,
709  const bool distortion,
710  Real distortion_tol = TOLERANCE) {
711 
712  CPPUNIT_ASSERT(gold_mesh.n_nodes() == mesh.n_nodes());
713  for (const auto node_id : make_range(gold_mesh.n_nodes()))
714  {
715  const auto & gold_node = gold_mesh.node_ref(node_id);
716  const auto & node = mesh.node_ref(node_id);
717  for (const auto d : make_range(gold_mesh.mesh_dimension()))
718  {
719  const bool d_distorted = std::abs(gold_node(d) - node(d)) > distortion_tol;
720  if (d_distorted)
721  // mesh is distorted from gold_mesh
722  return distortion;
723  }
724  }
725 
726  // mesh is not distorted from gold_mesh
727  return !distortion;
728  };
729 
730  // Make sure the mesh has been distorted
731  CPPUNIT_ASSERT(mesh_distortion_is(gold_mesh, mesh, true));
732 
733  // Turn off subdomain boundary preservation because DistortHyperCube
734  // does not preserve subdomain boundaries
735  // Also set dilation coefficient to 0 because the reference volume of the
736  // distorted mesh won't be equal to the reference volume of the smoothed
737  // mesh and will smooth to a slightly different solution.
738  VariationalMeshSmoother smoother(mesh, 0.0, false);
739  smoother.smooth();
740 
741  // Make sure the mesh has been smoothed to the gold mesh
742  CPPUNIT_ASSERT(mesh_distortion_is(gold_mesh, mesh, false));
743  }
744 
746  {
748  LaplaceMeshSmoother laplace(mesh);
749 
750  testLaplaceSmoother(mesh, laplace, QUAD4);
751  }
752 
754  {
756  LaplaceMeshSmoother laplace(mesh);
757 
758  testLaplaceSmoother(mesh, laplace, TRI3);
759  }
760 
761 #ifdef LIBMESH_ENABLE_VSMOOTHER
763  {
765  VariationalMeshSmoother variational(mesh);
766 
767  testVariationalSmoother(mesh, variational, EDGE2);
768  }
769 
771  {
773  VariationalMeshSmoother variational(mesh);
774 
775  testVariationalSmoother(mesh, variational, EDGE3);
776  }
777 
779  {
781  VariationalMeshSmoother variational(mesh);
782 
783  testVariationalSmoother(mesh, variational, EDGE3, true);
784  }
785 
787  {
789  VariationalMeshSmoother variational(mesh);
790 
791  testVariationalSmoother(mesh, variational, QUAD4);
792  }
793 
795  {
797  VariationalMeshSmoother variational(mesh);
798 
799  testVariationalSmoother(mesh, variational, QUAD4, true);
800  }
801 
803  {
805  VariationalMeshSmoother variational(mesh);
806 
807  testVariationalSmoother(mesh, variational, QUAD8);
808  }
809 
811  {
813  VariationalMeshSmoother variational(mesh);
814 
815  testVariationalSmoother(mesh, variational, QUAD9);
816  }
817 
819  {
821  VariationalMeshSmoother variational(mesh);
822 
823  testVariationalSmoother(mesh, variational, QUAD9, true);
824  }
825 
827  {
829  VariationalMeshSmoother variational(mesh);
830 
831  testVariationalSmoother(mesh, variational, QUAD4, false, true, 0.65);
832  }
833 
835  {
837  VariationalMeshSmoother variational(mesh);
838 
839  testVariationalSmoother(mesh, variational, TRI3);
840  }
841 
843  {
845  VariationalMeshSmoother variational(mesh);
846 
847  testVariationalSmoother(mesh, variational, TRI6);
848  }
849 
851  {
853  VariationalMeshSmoother variational(mesh);
854 
855  testVariationalSmoother(mesh, variational, TRI3, true);
856  }
857 
859  {
861  VariationalMeshSmoother variational(mesh);
862 
863  testVariationalSmoother(mesh, variational, HEX8);
864  }
865 
867  {
869  VariationalMeshSmoother variational(mesh);
870 
871  testVariationalSmoother(mesh, variational, HEX20);
872  }
873 
875  {
877  VariationalMeshSmoother variational(mesh);
878 
879  testVariationalSmoother(mesh, variational, HEX27);
880  }
881 
883  {
885  VariationalMeshSmoother variational(mesh);
886 
887  testVariationalSmoother(mesh, variational, HEX27, true);
888  }
889 
891  {
893  VariationalMeshSmoother variational(mesh);
894 
895  testVariationalSmoother(mesh, variational, HEX27, false, true, 0.55);
896  }
897 
899  {
901  VariationalMeshSmoother variational(mesh);
902 
903  testVariationalSmoother(mesh, variational, PYRAMID5);
904  }
905 
907  {
909  VariationalMeshSmoother variational(mesh);
910 
911  testVariationalSmoother(mesh, variational, PYRAMID13);
912  }
913 
915  {
917  VariationalMeshSmoother variational(mesh);
918 
919  testVariationalSmoother(mesh, variational, PYRAMID14);
920  }
921 
923  {
925  VariationalMeshSmoother variational(mesh);
926 
927  testVariationalSmoother(mesh, variational, PYRAMID18);
928  }
929 
931  {
933  VariationalMeshSmoother variational(mesh);
934 
935  testVariationalSmoother(mesh, variational, PYRAMID18, true);
936  }
937 
939  {
941  VariationalMeshSmoother variational(mesh);
942 
943  testVariationalSmoother(mesh, variational, PRISM6);
944  }
945 
947  {
949  VariationalMeshSmoother variational(mesh);
950 
951  testVariationalSmoother(mesh, variational, PRISM15);
952  }
953 
955  {
957  VariationalMeshSmoother variational(
958  mesh, 0.5, true, TOLERANCE * TOLERANCE, 10 * TOLERANCE * TOLERANCE);
959 
960  testVariationalSmoother(mesh, variational, PRISM18);
961  }
962 
964  {
966  VariationalMeshSmoother variational(
967  mesh, 0.5, true, TOLERANCE * TOLERANCE, 10 * TOLERANCE * TOLERANCE);
968 
969  testVariationalSmoother(mesh, variational, PRISM20);
970  }
971 
973  {
975  VariationalMeshSmoother variational(
976  mesh, 0.5, true, TOLERANCE * TOLERANCE, 10 * TOLERANCE * TOLERANCE);
977 
978  testVariationalSmoother(mesh, variational, PRISM21);
979  }
980 
982  {
984  mesh.allow_renumbering(false);
985  VariationalMeshSmoother variational(mesh);
986 
987  testVariationalSmoother(mesh, variational, PRISM21, true);
988  }
989 
991  {
993  mesh.read("meshes/quad4_tri3_smoothed.xda.gz");
994 
995  testVariationalSmootherRegression(mesh);
996  }
997 
999  {
1001  mesh.read("meshes/hex8_prism6_smoothed.xda.gz");
1002 
1003  testVariationalSmootherRegression(mesh);
1004  }
1005 #endif // LIBMESH_ENABLE_VSMOOTHER
1006 };
1007 
virtual dof_id_type n_nodes() const override final
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...
Definition: elem.h:992
ElemType
Defines an enum for geometric element types.
const std::set< Order > & elem_default_orders() const
Definition: mesh_base.h:289
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.
Definition: node.h:52
virtual void zero() override final
Set every element in the vector to 0.
Definition: dense_vector.h:420
void testLaplaceSmoother(ReplicatedMesh &mesh, MeshSmoother &smoother, ElemType type)
void allow_renumbering(bool allow)
If false is passed in then this mesh will no longer be renumbered when being prepared for use...
Definition: mesh_base.h:1196
void minloc(T &r, unsigned int &min_id) const
libMesh::Parallel::Communicator * TestCommWorld
Definition: driver.C:171
static constexpr Real TOLERANCE
unsigned int dim
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:396
MeshBase & mesh
const MeshQualityInfo & get_mesh_info() const
Getter for the _system&#39;s _mesh_info attribute.
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 build_square(UnstructuredMesh &mesh, const unsigned int nx, const unsigned int ny, const Real xmin=0., const Real xmax=1., const Real ymin=0., const Real ymax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
A specialized build_cube() for 2D meshes.
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.
Definition: mesh_base.h:165
virtual void smooth()=0
Function which actually performs the smoothing operations.
This class defines the data structures necessary for Laplace smoothing.
T pow(const T &x)
Definition: utility.h:328
static constexpr dof_id_type invalid_id
An invalid id to distinguish an uninitialized DofObject.
Definition: dof_object.h:482
void testVariationalQuad4Tangled()
libmesh_assert(ctx)
void testVariationalQuad9MultipleSubdomains()
CPPUNIT_TEST_SUITE_REGISTRATION(MeshSmootherTest)
auto norm(const T &a) -> decltype(std::abs(a))
Definition: tensor_tools.h:74
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.
Definition: fuzzy_equals.h:64
void testVariationalHex27MultipleSubdomains()
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void testVariationalSmootherRegression(const ReplicatedMesh &gold_mesh)
This is an implementation of Larisa Branets&#39; smoothing algorithms.
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...
Definition: int_range.h:140
unsigned int mesh_dimension() const
Definition: mesh_base.C:372
void build_line(UnstructuredMesh &mesh, const unsigned int nx, const Real xmin=0., const Real xmax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
A specialized build_cube() for 1D meshes.
void redistribute(MeshBase &mesh, const FunctionBase< Real > &mapfunc)
Deterministically perturb the nodal locations.
virtual const Node & node_ref(const dof_id_type i) const
Definition: mesh_base.h:596
This class provides the necessary interface for mesh smoothing.
Definition: mesh_smoother.h:38
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.
Definition: fuzzy_equals.h:78
Defines a dense vector for use in Finite Element-type computations.
Definition: dof_map.h:74
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.
Definition: point.h:39
const Elem & get(const ElemType type_in)
void build_cube(UnstructuredMesh &mesh, const unsigned int nx=0, const unsigned int ny=0, const unsigned int nz=0, const Real xmin=0., const Real xmax=1., const Real ymin=0., const Real ymax=1., const Real zmin=0., const Real zmax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
Builds a (elements) cube.
virtual dof_id_type n_nodes() const =0
virtual void smooth() override
Redefinition of the smooth function from the base class.
uint8_t dof_id_type
Definition: id_types.h:67
void testVariationalEdge3MultipleSubdomains()