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/utility.h"
17 #include "libmesh/enum_to_string.h"
18 #include "libmesh/parallel_ghost_sync.h"
19 
20 #include "test_comm.h"
21 #include "libmesh_cppunit.h"
22 
23 #include <random>
24 
25 namespace
26 {
27 using namespace libMesh;
28 using Utility::pow;
29 
30 // Distortion function that doesn't distort boundary nodes
31 // 2D only, use for LaplaceMeshSmoother
32 class DistortSquare : public FunctionBase<Real>
33 {
34  std::unique_ptr<FunctionBase<Real>> clone() const override
35  {
36  return std::make_unique<DistortSquare>();
37  }
38 
39  Real operator()(const Point &, const Real = 0.) override
40  {
41  libmesh_not_implemented();
42  } // scalar-only API
43 
44  // Skew inward based on a cubic function
45  void operator()(const Point & p, const Real, DenseVector<Real> & output)
46  {
47  output.resize(3);
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));
52  output(2) = 0;
53  }
54 };
55 
56 // Distortion function supporting 1D, 2D, and 3D
57 // Use for VariationalMeshSmoother
58 class DistortHyperCube : public FunctionBase<Real>
59 {
60 public:
61  DistortHyperCube(const unsigned int dim) : _dim(dim) {}
62 
63 private:
64  std::unique_ptr<FunctionBase<Real>> clone() const override
65  {
66  return std::make_unique<DistortHyperCube>(_dim);
67  }
68 
69  Real operator()(const Point &, const Real = 0.) override { libmesh_not_implemented(); }
70 
71  void operator()(const Point & p, const Real, DenseVector<Real> & output) override
72  {
73  output.resize(3);
74  output.zero();
75 
76  // Count how many coordinates are exactly on the boundary (0 or 1)
77  std::array<bool, 3> is_on_boundary = {false, false, false};
78  for (unsigned int i = 0; i < _dim; ++i)
79  if (std::abs(p(i)) < TOLERANCE || std::abs(p(i) - 1.) < TOLERANCE)
80  is_on_boundary[i] = true;
81 
82  // Distort only those directions not fixed on the boundary
83  for (unsigned int i = 0; i < _dim; ++i)
84  {
85  if (!is_on_boundary[i]) // only distort free dimensions
86  {
87  // This value controls the strength of the distortion
88  Real modulation = 0.3;
89  Real xi = 2. * p(i) - 1.;
90  for (unsigned int j = 0; j < _dim; ++j)
91  {
92  if (j != i)
93  {
94  Real pj = std::clamp(p(j), 0., 1.); // ensure numeric safety
95  modulation *= (pj - 0.5) * (pj - 0.5) * 4.; // quadratic bump centered at 0.5
96  }
97  }
98  const auto delta = (pow<3>(xi) - xi) * modulation;
99  // Check for delta = 0 to make sure we perturb every point
100  output(i) = (delta > TOLERANCE) ? p(i) + delta : 1.03 * p(i);
101  }
102  else
103  output(i) = p(i); // dimension on boundary remains unchanged
104  }
105  }
106 
107  const unsigned int _dim;
108 };
109 
110 class SquareToParallelogram : public FunctionBase<Real>
111 {
112  std::unique_ptr<FunctionBase<Real>> clone() const override
113  {
114  return std::make_unique<SquareToParallelogram>();
115  }
116 
117  Real operator()(const Point &, const Real = 0.) override
118  {
119  libmesh_not_implemented();
120  } // scalar-only API
121 
122  // Has the effect that a square, meshed into right triangles with diagonals
123  // rising from lower-left to upper-right, is transformed into a left-leaning
124  // parallelogram of eqilateral triangles.
125  void operator()(const Point & p, const Real, DenseVector<Real> & output)
126  {
127  output.resize(3);
128  output(0) = p(0) - 0.5 * p(1);
129  output(1) = 0.5 * std::sqrt(Real(3)) * p(1);
130  output(2) = 0;
131  }
132 };
133 
134 class ParallelogramToSquare : public FunctionBase<Real>
135 {
136  std::unique_ptr<FunctionBase<Real>> clone() const override
137  {
138  return std::make_unique<ParallelogramToSquare>();
139  }
140 
141  Real operator()(const Point &, const Real = 0.) override
142  {
143  libmesh_not_implemented();
144  } // scalar-only API
145 
146  // Has the effect that a left-leaning parallelogram of equilateral triangles is transformed
147  // into a square of right triangle with diagonals rising from lower-left to upper-right.
148  // This is the inversion of the SquareToParallelogram mapping.
149  void operator()(const Point & p, const Real, DenseVector<Real> & output)
150  {
151  output.resize(3);
152  output(0) = p(0) + p(1) / std::sqrt(Real(3));
153  output(1) = (2. / std::sqrt(Real(3))) * p(1);
154  output(2) = 0;
155  }
156 };
157 }
158 
159 using namespace libMesh;
160 
161 class MeshSmootherTest : public CppUnit::TestCase
162 {
167 public:
168  LIBMESH_CPPUNIT_TEST_SUITE(MeshSmootherTest);
169 
170 #if LIBMESH_DIM > 2
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);
177 
178  CPPUNIT_TEST(testVariationalQuad);
179  CPPUNIT_TEST(testVariationalQuadMultipleSubdomains);
180  CPPUNIT_TEST(testVariationalQuadTangled);
181 
182  CPPUNIT_TEST(testVariationalTri3);
183  CPPUNIT_TEST(testVariationalTri6);
184  CPPUNIT_TEST(testVariationalTri6MultipleSubdomains);
185 
186  CPPUNIT_TEST(testVariationalHex8);
187  CPPUNIT_TEST(testVariationalHex20);
188  CPPUNIT_TEST(testVariationalHex27);
189  CPPUNIT_TEST(testVariationalHex27Tangled);
190  CPPUNIT_TEST(testVariationalHex27MultipleSubdomains);
191 
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
198 #endif
199 
200  CPPUNIT_TEST_SUITE_END();
201 
202 public:
203  void setUp() {}
204 
205  void tearDown() {}
206 
208  {
209  LOG_UNIT_TEST;
210 
211  constexpr unsigned int n_elems_per_side = 4;
212 
214  mesh, n_elems_per_side, n_elems_per_side, 0., 1., 0., 1., type);
215 
216  // Move it around so we have something that needs smoothing
217  DistortSquare ds;
219 
220  // Assert the distortion is as expected
221  auto center_distortion_is =
222  [](const Node & node, int d, bool distortion, Real distortion_tol = TOLERANCE) {
223  const Real r = node(d);
224  const Real R = r * n_elems_per_side;
225  CPPUNIT_ASSERT_GREATER(-TOLERANCE * TOLERANCE, r);
226  CPPUNIT_ASSERT_GREATER(-TOLERANCE * TOLERANCE, 1 - r);
227 
228  // If we're at the boundaries we should *never* be distorted
229  if (std::abs(node(0)) < TOLERANCE * TOLERANCE ||
230  std::abs(node(0) - 1) < TOLERANCE * TOLERANCE)
231  {
232  const Real R1 = node(1) * n_elems_per_side;
233  CPPUNIT_ASSERT_LESS(TOLERANCE * TOLERANCE, std::abs(R1 - std::round(R1)));
234  return true;
235  }
236 
237  if (std::abs(node(1)) < TOLERANCE * TOLERANCE ||
238  std::abs(node(1) - 1) < TOLERANCE * TOLERANCE)
239  {
240  const Real R0 = node(0) * n_elems_per_side;
241  CPPUNIT_ASSERT_LESS(TOLERANCE * TOLERANCE, std::abs(R0 - std::round(R0)));
242 
243  return true;
244  }
245 
246  // If we're at the center we're fine
247  if (std::abs(r - 0.5) < TOLERANCE * TOLERANCE)
248  return true;
249 
250  return ((std::abs(R - std::round(R)) > distortion_tol) == distortion);
251  };
252 
253  for (auto node : mesh.node_ptr_range())
254  {
255  CPPUNIT_ASSERT(center_distortion_is(*node, 0, true));
256  CPPUNIT_ASSERT(center_distortion_is(*node, 1, true));
257  }
258 
259  // Enough iterations to mostly fix us up. Laplace seems to be at 1e-3
260  // tolerance by iteration 6, so hopefully everything is there on any
261  // system by 8.
262  for (unsigned int i = 0; i != 8; ++i)
263  smoother.smooth();
264 
265  // Make sure we're not too distorted anymore.
266  for (auto node : mesh.node_ptr_range())
267  {
268  CPPUNIT_ASSERT(center_distortion_is(*node, 0, false, 1e-3));
269  CPPUNIT_ASSERT(center_distortion_is(*node, 1, false, 1e-3));
270  }
271  }
272 
274  VariationalMeshSmoother & smoother,
275  const ElemType type,
276  const bool multiple_subdomains = false,
277  const bool tangle_mesh=false,
278  const Real tangle_damping_factor=1.0)
279  {
280  LOG_UNIT_TEST;
281 
282  if (multiple_subdomains && tangle_mesh)
283  libmesh_not_implemented_msg(
284  "Arbitrary mesh tangling with multiple subdomains is not supported.");
285 
286  // Get mesh dimension, determine whether type is triangular
287  const auto * ref_elem = &(ReferenceElem::get(type));
288  const auto dim = ref_elem->dim();
289  const bool type_is_tri = Utility::enum_to_string(type).compare(0, 3, "TRI") == 0;
290  const bool type_is_pyramid = Utility::enum_to_string(type).compare(0, 7, "PYRAMID") == 0;
291 
292  // Used fewer elems for higher order types, as extra midpoint nodes will add
293  // enough complexity
294  unsigned int n_elems_per_side = 4 / Elem::type_to_default_order_map[type];
295 
296  switch (dim)
297  {
298  case 1:
299  MeshTools::Generation::build_line(mesh, n_elems_per_side, 0., 1., type);
300  break;
301  case 2:
303  mesh, n_elems_per_side, n_elems_per_side, 0., 1., 0., 1., type);
304  break;
305 
306  case 3:
308  n_elems_per_side,
309  n_elems_per_side,
310  n_elems_per_side,
311  0.,
312  1.,
313  0.,
314  1.,
315  0.,
316  1.,
317  type);
318  break;
319 
320  default:
321  libmesh_error_msg("Unsupported dimension " << dim);
322  }
323 
324  // Move it around so we have something that needs smoothing
325  DistortHyperCube dh(dim);
327 
328  const auto & boundary_info = mesh.get_boundary_info();
329 
330  if (tangle_mesh)
331  {
332  // We tangle the mesh by (partially) swapping the locations of 2 nodes.
333  // If the n-dimentional hypercube is represented as a grid with
334  // (undistorted) divisions occuring at integer numbers, we will
335  // (partially) swap the nodes nearest p1 = (1,1,1) and p2 = (2,1,2).
336 
337  // Define p1 and p2 given the integer indices
338  // Undistorted elem side length
339  const Real dr = 1. / (n_elems_per_side * Elem::type_to_default_order_map[type]);
340  const Point p1 = Point(dr, dim > 1 ? dr : 0, dim > 2 ? dr : 0);
341  const Point p2 = Point(2. * dr, dim > 1 ? dr : 0, dim > 2 ? 2. * dr : 0);
342 
343  // ids and distances of mesh nodes nearest p1 and p2
345  Real dist1_closest = std::numeric_limits<Real>::max();
347  Real dist2_closest = std::numeric_limits<Real>::max();
348 
349  for (const auto * node : mesh.local_node_ptr_range())
350  {
351  const auto dist1 = (*node - p1).norm();
352  if (dist1 < dist1_closest)
353  {
354  dist1_closest = dist1;
355  node1_id = node->id();
356  }
357 
358  const auto dist2 = (*node - p2).norm();
359  if (dist2 < dist2_closest)
360  {
361  dist2_closest = dist2;
362  node2_id = node->id();
363  }
364  }
365 
366  // parallel communication
367  unsigned int node1_rank;
368  mesh.comm().minloc(dist1_closest, node1_rank);
369  mesh.comm().broadcast(node1_id, node1_rank);
370 
371  unsigned int node2_rank;
372  mesh.comm().minloc(dist2_closest, node2_rank);
373  mesh.comm().broadcast(node2_id, node2_rank);
374 
375  // modify the nodes
376  if ((node1_id != DofObject::invalid_id) && (node2_id != DofObject::invalid_id))
377  {
378  libmesh_assert(node1_id != node2_id);
379 
380  auto & node1 = mesh.node_ref(node1_id);
381  auto & node2 = mesh.node_ref(node2_id);
382 
383  const Point diff = node1 - node2;
384  // Note that a damping factor of 1 will swap the nodes and one
385  // of 0.5 will move the nodes to the mean of their original
386  // locations.
387  node1 -= tangle_damping_factor * diff;
388  node2 += tangle_damping_factor * diff;
389  }
390 
391  SyncNodalPositions sync_object(mesh);
392  Parallel::sync_dofobject_data_by_id (mesh.comm(), mesh.nodes_begin(), mesh.nodes_end(), sync_object);
393 
394  // Make sure the mesh is tangled
395  smoother.setup(); // Need this to create the system we are about to access
396  const auto & unsmoothed_info = smoother.get_mesh_info();
397  CPPUNIT_ASSERT(unsmoothed_info.mesh_is_tangled);
398  }
399 
400  // Add multiple subdomains if requested
401  std::unordered_map<subdomain_id_type, Real> distorted_subdomain_volumes;
402  if (multiple_subdomains)
403  {
404  // Modify the subdomain ids in an interesting way
405  for (auto * elem : mesh.active_element_ptr_range())
406  {
407  unsigned int subdomain_id = 0;
408  for (const auto d : make_range(dim))
409  if (elem->vertex_average()(d) > 0.5)
410  ++subdomain_id;
411  elem->subdomain_id() += subdomain_id;
412  }
413 
414  // This loop should NOT be combined with the one above because we need to
415  // finish checking and updating subdomain ids for all elements before
416  // recording the final subdomain volumes.
417  for (auto * elem : mesh.active_element_ptr_range())
418  distorted_subdomain_volumes[elem->subdomain_id()] += elem->volume();
419  }
420 
421  // Get the mesh order
422  const auto & elem_orders = mesh.elem_default_orders();
423  libmesh_error_msg_if(elem_orders.size() != 1,
424  "The variational smoother cannot be used for mixed-order meshes!");
425  // For pyramids, the factor 2 accounts for the account that cubes of side
426  // length s are divided into pyramids of base side length s and height s/2.
427  // The factor 4 lets us catch triangular face midpoint nodes for PYRAMID18 elements.
428  const auto scale_factor = *elem_orders.begin() * (type_is_pyramid ? 2 * 4 : 1);
429 
430  // Function to assert the node distortion is as expected
431  auto node_distortion_is = [&n_elems_per_side, &dim, &boundary_info, &scale_factor](
432  const Node & node, bool distortion, Real distortion_tol = TOLERANCE) {
433  // Get boundary ids associated with the node
434  std::vector<boundary_id_type> boundary_ids;
435  boundary_info.boundary_ids(&node, boundary_ids);
436 
437  // This tells us what type of node we are: internal, sliding, or fixed
438  const auto num_dofs = dim - boundary_ids.size();
439  /*
440  * The following cases of num_dofs are possible, ASSUMING all boundaries
441  * are non-overlapping
442  * 3D: 3-0, 3-1, 3-2, 3-3
443  * = 3 2 1 0
444  * internal sliding, sliding, fixed
445  * 2D: 2-0, 2-1, 2-2
446  * = 2 1 0
447  * internal sliding, fixed
448  * 1D: 1-0, 1-1
449  * = 1 0
450  * internal fixed
451  *
452  * We expect that R is an integer in [0, n_elems_per_side] for
453  * num_dofs of the node's cooridinantes, while the remaining coordinates
454  * are fixed to the boundary with value 0 or 1. In other words, at LEAST
455  * dim - num_dofs coordinantes should be 0 or 1.
456  */
457 
458  std::size_t num_zero_or_one = 0;
459 
460  bool distorted = false;
461  for (const auto d : make_range(dim))
462  {
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);
467 
468  const bool d_distorted = std::abs(R - std::round(R)) > distortion_tol;
469  distorted |= d_distorted;
470  num_zero_or_one += (absolute_fuzzy_equals(r, 0.) || absolute_fuzzy_equals(r, 1.));
471  }
472 
473  CPPUNIT_ASSERT_GREATEREQUAL(dim - num_dofs, num_zero_or_one);
474 
475  // We can never expect a fixed node to be distorted
476  if (num_dofs == 0)
477  // if (num_dofs < dim)
478  return true;
479  return distorted == distortion;
480  };
481 
482  // Make sure our DistortHyperCube transformation has distorted the mesh
483  for (auto node : mesh.node_ptr_range())
484  CPPUNIT_ASSERT(node_distortion_is(*node, true));
485 
486  // Transform the square mesh of triangles to a parallelogram mesh of
487  // triangles. This will allow the Variational Smoother to smooth the mesh
488  // to the optimal case of equilateral triangles
489  if (type_is_tri)
490  {
491  SquareToParallelogram stp;
493  }
494 
495  smoother.smooth();
496 
497  // Transform the parallelogram mesh back to a square mesh. In the case of
498  // the Variational Smoother, equilateral triangular elements will be
499  // transformed into right triangular elements that align with the original
500  // undistorted mesh.
501  if (type_is_tri)
502  {
503  ParallelogramToSquare pts;
505  }
506 
507  if (multiple_subdomains)
508  {
509  // Make sure the subdomain volumes didn't change. Although this doesn't
510  // guarantee that the subdomain didn't change, it is a good indicator
511  // and is friedly to sliding subdomain boundary nodes.
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();
515 
516  for (const auto & pair : distorted_subdomain_volumes)
517  {
518  const auto & subdomain_id = pair.first;
519  CPPUNIT_ASSERT(
520  relative_fuzzy_equals(libmesh_map_find(distorted_subdomain_volumes, subdomain_id),
521  libmesh_map_find(smoothed_subdomain_volumes, subdomain_id),
522  TOLERANCE));
523  }
524  }
525  else
526  {
527  // Make sure we're not too distorted anymore
528  std::set<dof_id_type> nodes_checked;
529  for (const auto * elem : mesh.active_element_ptr_range())
530  {
531  for (const auto local_node_id : make_range(elem->n_nodes()))
532  {
533  const auto & node = elem->node_ref(local_node_id);
534  if (nodes_checked.find(node.id()) != nodes_checked.end())
535  continue;
536 
537  nodes_checked.insert(node.id());
538 
539  // Check for special case where pyramidal base-to-apex midpoint
540  // nodes are not smoothed to the actual midpoints.
541  if (type_is_pyramid)
542  {
543  if (local_node_id > 8 && local_node_id < 13)
544  {
545  // Base-Apex midpoint nodes
546  //
547  // Due to the nature of the pyramidal target element and the
548  // cubic nature of the test mesh, for a dilation weight o
549  // 0.5, smoothed base-apex midpoint nodes are smoothed to
550  // the point base + x * (apex - base) instead of
551  // base + 0.5 (apex - base), where x is in (0,1) and
552  // depends on the number of nodes in the element.
553  // We hard-code this check below.
554 
555  const auto & base = elem->node_ref(local_node_id - 9);
556  const auto & apex = elem->node_ref(4);
557  const Real x = (type == PYRAMID18) ? 0.569332 : 0.549876;
558 
559  CPPUNIT_ASSERT(node.relative_fuzzy_equals(base + x * (apex - base), 1e-3));
560  continue;
561  }
562  else if (local_node_id > 13)
563  {
564  // Triangular face midpoint nodes
565  //
566  // Due to the nature of the pyramidal target element and the
567  // cubic nature of the test mesh, for a dilation weight o
568  // 0.5, smoothed triangular face midpoint nodes are
569  // smoothed a weighted average of the constituent
570  // vertices instead of an unweighted average.
571  // We hard-code this check below.
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);
575 
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));
580  continue;
581  }
582  }
583 
584  CPPUNIT_ASSERT(node_distortion_is(node, false, 1e-2));
585 
586  }
587  }
588  }
589  }
590 
592  {
594  LaplaceMeshSmoother laplace(mesh);
595 
596  testLaplaceSmoother(mesh, laplace, QUAD4);
597  }
598 
600  {
602  LaplaceMeshSmoother laplace(mesh);
603 
604  testLaplaceSmoother(mesh, laplace, TRI3);
605  }
606 
607 #ifdef LIBMESH_ENABLE_VSMOOTHER
609  {
611  VariationalMeshSmoother variational(mesh);
612 
613  testVariationalSmoother(mesh, variational, EDGE2);
614  }
615 
617  {
619  VariationalMeshSmoother variational(mesh);
620 
621  testVariationalSmoother(mesh, variational, EDGE3);
622  }
623 
625  {
627  VariationalMeshSmoother variational(mesh);
628 
629  testVariationalSmoother(mesh, variational, EDGE3, true);
630  }
631 
633  {
635  VariationalMeshSmoother variational(mesh);
636 
637  testVariationalSmoother(mesh, variational, QUAD4);
638  }
639 
641  {
643  VariationalMeshSmoother variational(mesh);
644 
645  testVariationalSmoother(mesh, variational, QUAD4, true);
646  }
647 
649  {
651  VariationalMeshSmoother variational(mesh);
652 
653  testVariationalSmoother(mesh, variational, QUAD4, false, true, 0.65);
654  }
655 
657  {
659  VariationalMeshSmoother variational(mesh);
660 
661  testVariationalSmoother(mesh, variational, TRI3);
662  }
663 
665  {
667  VariationalMeshSmoother variational(mesh);
668 
669  testVariationalSmoother(mesh, variational, TRI6);
670  }
671 
673  {
675  VariationalMeshSmoother variational(mesh);
676 
677  testVariationalSmoother(mesh, variational, TRI3, true);
678  }
679 
681  {
683  VariationalMeshSmoother variational(mesh);
684 
685  testVariationalSmoother(mesh, variational, HEX8);
686  }
687 
689  {
691  VariationalMeshSmoother variational(mesh);
692 
693  testVariationalSmoother(mesh, variational, HEX20);
694  }
695 
697  {
699  VariationalMeshSmoother variational(mesh);
700 
701  testVariationalSmoother(mesh, variational, HEX27);
702  }
703 
705  {
707  VariationalMeshSmoother variational(mesh);
708 
709  testVariationalSmoother(mesh, variational, HEX27, true);
710  }
711 
713  {
715  VariationalMeshSmoother variational(mesh);
716 
717  testVariationalSmoother(mesh, variational, HEX27, false, true, 0.55);
718  }
719 
721  {
723  VariationalMeshSmoother variational(mesh);
724 
725  testVariationalSmoother(mesh, variational, PYRAMID5);
726  }
727 
729  {
731  VariationalMeshSmoother variational(mesh);
732 
733  testVariationalSmoother(mesh, variational, PYRAMID13);
734  }
735 
737  {
739  VariationalMeshSmoother variational(mesh);
740 
741  testVariationalSmoother(mesh, variational, PYRAMID14);
742  }
743 
745  {
747  VariationalMeshSmoother variational(mesh);
748 
749  testVariationalSmoother(mesh, variational, PYRAMID18);
750  }
751 
753  {
755  VariationalMeshSmoother variational(mesh);
756 
757  testVariationalSmoother(mesh, variational, PYRAMID18, true);
758  }
759 #endif // LIBMESH_ENABLE_VSMOOTHER
760 };
761 
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:991
ElemType
Defines an enum for geometric element types.
const std::set< Order > & elem_default_orders() const
Definition: mesh_base.h:289
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 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 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
libmesh_assert(ctx)
static const dof_id_type invalid_id
An invalid id to distinguish an uninitialized DofObject.
Definition: dof_object.h:482
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
This is an implementation of Larisa Branets&#39; 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...
Definition: int_range.h:140
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...
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 void smooth() override
Redefinition of the smooth function from the base class.
uint8_t dof_id_type
Definition: id_types.h:67
void testVariationalEdge3MultipleSubdomains()