libMesh
Public Member Functions | Protected Attributes | List of all members
ElemTest< elem_type > Class Template Reference
Inheritance diagram for ElemTest< elem_type >:
[legend]

Public Member Functions

void test_bounding_box ()
 
void test_quality ()
 
void test_maps ()
 
void test_static_data ()
 
void test_contains_point_node ()
 
void test_permute ()
 
void test_flip ()
 
void test_orient ()
 
void test_orient_elements ()
 
void test_center_node_on_side ()
 
void test_side_type ()
 
void test_side_subdomain ()
 
void test_elem_side_builder ()
 
void test_n_refinements (unsigned int n)
 
void test_refinement ()
 
void test_double_refinement ()
 
void test_is_internal ()
 
void test_node_edge_map_consistency ()
 
void setUp ()
 

Protected Attributes

std::unique_ptr< Mesh > _mesh
 
std::string libmesh_suite_name
 

Detailed Description

template<ElemType elem_type>
class ElemTest< elem_type >

Definition at line 14 of file elem_test.C.

Member Function Documentation

◆ setUp()

template<ElemType elem_type>
void PerElemTest< elem_type >::setUp ( )
inlineinherited

Definition at line 26 of file elem_test.h.

References libMesh::Elem::build(), libMesh::MeshTools::Generation::build_cube(), libMesh::C0POLYGON, libMesh::C0POLYHEDRON, dim, libMesh::index_range(), libMesh::INFHEX16, libMesh::INFHEX18, libMesh::INFHEX8, libMesh::INFPRISM12, libMesh::INFPRISM6, libMesh::INFQUAD4, libMesh::INFQUAD6, libMesh::make_range(), libMesh::Real, and libMesh::Elem::set_node().

27  {
28  const Real minpos = 1.5, maxpos = 5.5;
29  const unsigned int N = 2;
30 
31  _mesh = std::make_unique<Mesh>(*TestCommWorld);
32 
33  std::unique_ptr<Elem> test_elem;
34 
35  if (elem_type != C0POLYHEDRON)
36  test_elem = Elem::build(elem_type);
37 
38 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
39 #if LIBMESH_DIM > 1
40  if (test_elem.get() && test_elem->infinite())
41  {
42  Elem * elem = _mesh->add_elem(std::move(test_elem));
43 
44  const auto add_point =
45  [this, elem](const unsigned int i,
46  const Real x,
47  const Real y,
48  const Real
49 #if LIBMESH_DIM == 3
50  z
51 #endif
52  )
53  {
54 #if LIBMESH_DIM == 2
55  auto node = _mesh->add_point(Point(x, y), i);
56 #else
57  auto node = _mesh->add_point(Point(x, y, z), i);
58 #endif
59  elem->set_node(i, node);
60  };
61 
62  const Real halfpos = (minpos + maxpos) / 2.;
63 
64  if (elem_type == INFQUAD4 || elem_type == INFQUAD6 ||
65  elem_type == INFHEX8 || elem_type == INFHEX16 || elem_type == INFHEX18)
66  {
67  const bool is_quad = (elem_type == INFQUAD4 || elem_type == INFQUAD6);
68 
69  add_point(0, minpos, minpos, minpos);
70  add_point(1, maxpos, minpos, minpos);
71  add_point(2+is_quad, maxpos, maxpos, minpos);
72  add_point(3-is_quad, minpos, maxpos, minpos);
73 
74  if (elem_type == INFQUAD6)
75  {
76  add_point(4, halfpos, minpos, minpos);
77  add_point(5, halfpos, maxpos, minpos);
78  }
79  }
80  if (elem_type == INFHEX8 || elem_type == INFHEX16 || elem_type == INFHEX18)
81  {
82  add_point(4, minpos, minpos, maxpos);
83  add_point(5, maxpos, minpos, maxpos);
84  add_point(6, maxpos, maxpos, maxpos);
85  add_point(7, minpos, maxpos, maxpos);
86 
87  if (elem_type == INFHEX16 || elem_type == INFHEX18)
88  {
89  add_point(8, halfpos, minpos, minpos);
90  add_point(9, maxpos, halfpos, minpos);
91  add_point(10, halfpos, maxpos, minpos);
92  add_point(11, minpos, halfpos, minpos);
93  add_point(12, halfpos, minpos, maxpos);
94  add_point(13, maxpos, halfpos, maxpos);
95  add_point(14, halfpos, maxpos, maxpos);
96  add_point(15, minpos, halfpos, maxpos);
97  }
98  if (elem_type == INFHEX18)
99  {
100  add_point(16, halfpos, halfpos, minpos);
101  add_point(17, halfpos, halfpos, maxpos);
102  }
103  }
104  if (elem_type == INFPRISM6 || elem_type == INFPRISM12)
105  {
106  add_point(0, minpos, minpos, minpos);
107  add_point(1, maxpos, minpos, minpos);
108  add_point(2, halfpos, maxpos, minpos);
109  add_point(3, minpos, minpos, maxpos);
110  add_point(4, maxpos, minpos, maxpos);
111  add_point(5, halfpos, maxpos, maxpos);
112 
113  if (elem_type == INFPRISM12)
114  {
115  add_point(6, halfpos, minpos, minpos);
116  add_point(7, (halfpos + maxpos) / 2., halfpos, minpos);
117  add_point(8, (halfpos + minpos) / 2., halfpos, minpos);
118  add_point(9, halfpos, minpos, maxpos);
119  add_point(10, (halfpos + maxpos) / 2., halfpos, maxpos);
120  add_point(11, (halfpos + minpos) / 2., halfpos, maxpos);
121  }
122  }
123 
124  _mesh->prepare_for_use();
125  }
126  else
127 #endif // LIBMESH_DIM > 1
128 #endif // LIBMESH_ENABLE_INFINITE_ELEMENTS
129  if (elem_type == C0POLYGON)
130  {
131  // We're not going to implement build_square for e.g.
132  // pentagons any time soon.
133  //
134  // We should probably implement build_dual_mesh, though, and
135  // run it on a perturbed or triangle input so we don't just
136  // get quads in the output...
137 
138  _mesh->add_point(Point(0, 0), 0);
139  _mesh->add_point(Point(1, 0), 1);
140  _mesh->add_point(Point(1.5, 0.5), 2);
141  _mesh->add_point(Point(1, 1), 3);
142  _mesh->add_point(Point(0, 1), 4);
143 
144  std::unique_ptr<Elem> polygon = std::make_unique<C0Polygon>(5);
145  for (auto i : make_range(5))
146  polygon->set_node(i, _mesh->node_ptr(i));
147  polygon->set_id() = 0;
148 
149  _mesh->add_elem(std::move(polygon));
150  _mesh->prepare_for_use();
151  }
152  else if (elem_type == C0POLYHEDRON)
153  {
154  // There's even less point in having a build_cube for general
155  // polyhedra, so again we'll hand-make one for testing and
156  // we'll plan on making a build_dual_mesh for the future.
157 
158 #if 0
159  // Try a truncated cube for relatively simple verification.
160 
161  // We have some differing orientations on the top sides, to
162  // test handling of that.
163  // Lower octagon points, counterclockwise
164  _mesh->add_point(Point(1/Real(3), 0, 0), 0);
165  _mesh->add_point(Point(2/Real(3), 0, 0), 1);
166  _mesh->add_point(Point(1, 1/Real(3), 0), 2);
167  _mesh->add_point(Point(1, 2/Real(3), 0), 3);
168  _mesh->add_point(Point(2/Real(3), 1, 0), 4);
169  _mesh->add_point(Point(1/Real(3), 1, 0), 5);
170  _mesh->add_point(Point(0, 2/Real(3), 0), 6);
171  _mesh->add_point(Point(0, 1/Real(3), 0), 7);
172 
173  // Lower-middle points, counterclockwise
174  _mesh->add_point(Point(0, 0, 1/Real(3)), 8);
175  _mesh->add_point(Point(1, 0, 1/Real(3)), 9);
176  _mesh->add_point(Point(1, 1, 1/Real(3)), 10);
177  _mesh->add_point(Point(0, 1, 1/Real(3)), 11);
178 
179  // Upper-middle points, counterclockwise
180  _mesh->add_point(Point(0, 0, 2/Real(3)), 12);
181  _mesh->add_point(Point(1, 0, 2/Real(3)), 13);
182  _mesh->add_point(Point(1, 1, 2/Real(3)), 14);
183  _mesh->add_point(Point(0, 1, 2/Real(3)), 15);
184 
185  // Upper octagon points, counterclockwise
186  _mesh->add_point(Point(1/Real(3), 0, 1), 16);
187  _mesh->add_point(Point(2/Real(3), 0, 1), 17);
188  _mesh->add_point(Point(1, 1/Real(3), 1), 18);
189  _mesh->add_point(Point(1, 2/Real(3), 1), 19);
190  _mesh->add_point(Point(2/Real(3), 1, 1), 20);
191  _mesh->add_point(Point(1/Real(3), 1, 1), 21);
192  _mesh->add_point(Point(0, 2/Real(3), 1), 22);
193  _mesh->add_point(Point(0, 1/Real(3), 1), 23);
194 
195  const std::vector<std::vector<unsigned int>> nodes_on_side =
196  { {0, 1, 2, 3, 4, 5, 6, 7}, // min z
197  {0, 1, 9, 13, 17, 16, 12, 8}, // min y
198  {2, 3, 10, 14, 19, 18, 13, 9}, // max x
199  {4, 5, 11, 15, 21, 20, 14, 10}, // max y
200  {6, 7, 8, 12, 23, 22, 15, 11}, // min x
201  {16, 17, 18, 19, 20, 21, 22, 23}, // max z
202  {7, 0, 8}, // max nothing
203  {1, 2, 9}, // max x
204  {3, 4, 10}, // max xy
205  {5, 6, 11}, // max y
206  {23, 16, 12}, // max z
207  {17, 18, 13}, // max xz
208  {19, 20, 14}, // max xyz
209  {21, 22, 15} }; // max yz
210 #endif
211 
212 #if 1
213  // Or try a plain cube, for simpler debugging
214  _mesh->add_point(Point(0, 0, 0), 0);
215  _mesh->add_point(Point(1, 0, 0), 1);
216  _mesh->add_point(Point(1, 1, 0), 2);
217  _mesh->add_point(Point(0, 1, 0), 3);
218  _mesh->add_point(Point(0, 0, 1), 4);
219  _mesh->add_point(Point(1, 0, 1), 5);
220  _mesh->add_point(Point(1, 1, 1), 6);
221  _mesh->add_point(Point(0, 1, 1), 7);
222 
223  // With some combinations of face triangulations, even a
224  // simple cube has no tetrahedralization that doesn't have
225  // either interior discontinuities or a 0-volume pancake tet!
226  //
227  // The initial "natural" way to orient our sides is commented
228  // out here; permutations that fix diagonals for us are used
229  // instead.
230  const std::vector<std::vector<unsigned int>> nodes_on_side =
231  { {0, 1, 2, 3}, // min z
232  {0, 1, 5, 4}, // min y
233  // {1, 2, 6, 5}, // max x - bad
234  {2, 6, 5, 1}, // max x
235  {2, 3, 7, 6}, // max y
236  // {3, 0, 4, 7}, // min x - bad
237  {0, 4, 7, 3}, // min x
238  // {4, 5, 6, 7} }; // max z - bad
239  {5, 6, 7, 4} }; // max z
240 #endif
241 
242  // Build all the sides.
243  std::vector<std::shared_ptr<Polygon>> sides(nodes_on_side.size());
244 
245  for (auto s : index_range(nodes_on_side))
246  {
247  const auto & nodes_on_s = nodes_on_side[s];
248  sides[s] = std::make_shared<C0Polygon>(nodes_on_s.size());
249  for (auto i : index_range(nodes_on_s))
250  sides[s]->set_node(i, _mesh->node_ptr(nodes_on_s[i]));
251  }
252 
253  std::unique_ptr<Elem> polyhedron = std::make_unique<C0Polyhedron>(sides);
254  _mesh->add_elem(std::move(polyhedron));
255  _mesh->prepare_for_use();
256  }
257  else
258  {
259  const unsigned int dim = test_elem->dim();
260  const unsigned int use_x = dim > 0;
261  const unsigned int use_y = dim > 1;
262  const unsigned int use_z = dim > 2;
263 
264  MeshTools::Generation::build_cube (*_mesh,
265  N*use_x, N*use_y, N*use_z,
266  minpos, maxpos,
267  minpos, use_y*maxpos,
268  minpos, use_z*maxpos,
269  elem_type);
270  }
271 
272  // Use non-default subdomain ids so we can properly test their
273  // preservation
274  for (const auto & elem :
275  this->_mesh->element_ptr_range())
276  elem->subdomain_id() = 10 + (elem->id() % 10);
277 
278  // And make sure our mesh's cache knows about them all for later
279  // tests comparisons
280  this->_mesh->cache_elem_data();
281  }
libMesh::Elem::set_node
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:2558
PerElemTest::_mesh
std::unique_ptr< Mesh > _mesh
Definition: elem_test.h:22
dim
unsigned int dim
Definition: adaptivity_ex3.C:113
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
libMesh::DofObject::id
dof_id_type id() const
Definition: dof_object.h:828
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:144
libMesh::Elem::subdomain_id
subdomain_id_type subdomain_id() const
Definition: elem.h:2582
libMesh::make_range
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
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
libMesh::index_range
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:117
libMesh::MeshTools::Generation::build_cube
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.
Definition: mesh_generation.C:313

◆ test_bounding_box()

template<ElemType elem_type>
void ElemTest< elem_type >::test_bounding_box ( )
inline

Definition at line 16 of file elem_test.C.

References libMesh::BoundingBox::contains_point().

17  {
18  LOG_UNIT_TEST;
19 
20  for (const auto & elem :
21  this->_mesh->active_local_element_ptr_range())
22  {
23  const BoundingBox bbox = elem->loose_bounding_box();
24 
25  // The "loose" bounding box should actually be pretty tight
26  // in most of these cases, but for weirdly aligned triangles
27  // (such as occur in pyramid elements) it won't be, so we'll
28  // just test against a widened bounding box.
29  BoundingBox wide_bbox(elem->point(0), elem->point(0));
30 
31  for (unsigned int n = 0; n != elem->n_nodes(); ++n)
32  {
33  const Point & p = elem->point(n);
34 
35  if (!elem->infinite())
36  CPPUNIT_ASSERT(bbox.contains_point(p));
37 
38  wide_bbox.union_with
39  (BoundingBox(elem->point(n), elem->point(n)));
40  }
41 
42  wide_bbox.scale(1. / 3.);
43 
44  if (!elem->infinite() && elem->dim())
45  {
46  CPPUNIT_ASSERT(!bbox.contains_point(wide_bbox.min()));
47  CPPUNIT_ASSERT(!bbox.contains_point(wide_bbox.max()));
48  }
49  }
50  }
libMesh::BoundingBox::contains_point
bool contains_point(const Point &) const
Definition: bounding_box.C:35
PerElemTest::_mesh
std::unique_ptr< Mesh > _mesh
Definition: elem_test.h:22
libMesh::BoundingBox
Defines a Cartesian bounding box by the two corner extremum.
Definition: bounding_box.h:40
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39

◆ test_center_node_on_side()

template<ElemType elem_type>
void ElemTest< elem_type >::test_center_node_on_side ( )
inline

Definition at line 501 of file elem_test.C.

References libMesh::EDGE2, libMesh::EDGE3, libMesh::EDGE4, libMesh::HEX27, libMesh::invalid_uint, libMesh::PRISM18, libMesh::PRISM20, libMesh::PRISM21, libMesh::PYRAMID14, libMesh::PYRAMID18, libMesh::QUAD8, libMesh::QUAD9, libMesh::QUADSHELL8, libMesh::QUADSHELL9, libMesh::TRI6, and libMesh::TRI7.

502  {
503  LOG_UNIT_TEST;
504 
505  for (const auto & elem :
506  this->_mesh->active_local_element_ptr_range())
507  for (const auto s : elem->side_index_range())
508  {
509  if (elem->type() == EDGE2 || elem->type() == EDGE3 || elem->type() == EDGE4)
510  CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(s), elem->center_node_on_side(s));
511  else if (elem->type() == TRI6 || elem->type() == TRI7)
512  CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(s + 3), elem->center_node_on_side(s));
513  else if (elem->type() == QUAD8 || elem->type() == QUAD9 ||
514  elem->type() == QUADSHELL8 || elem->type() == QUADSHELL9)
515  CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(s + 4), elem->center_node_on_side(s));
516  else if (elem->type() == HEX27)
517  CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(s + 20), elem->center_node_on_side(s));
518  else if (elem->type() == PRISM18 && s >= 1 && s <= 3)
519  CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(s + 14), elem->center_node_on_side(s));
520  else if ((elem->type() == PRISM20 ||
521  elem->type() == PRISM21) && s >= 1 && s <= 3)
522  CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(s + 14), elem->center_node_on_side(s));
523  else if (elem->type() == PRISM20 ||
524  elem->type() == PRISM21)
525  CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(18 + (s == 4)), elem->center_node_on_side(s));
526  else if (elem->type() == PYRAMID14 && s == 4)
527  CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(13), elem->center_node_on_side(s));
528  else if (elem->type() == PYRAMID18)
529  {
530  if (s < 4)
531  CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(s + 14), elem->center_node_on_side(s));
532  else
533  CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(13), elem->center_node_on_side(s));
534  }
535  else
536  CPPUNIT_ASSERT_EQUAL(invalid_uint, elem->center_node_on_side(s));
537  }
538  }
libMesh::invalid_uint
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value for an unsigned int...
Definition: libmesh.h:310
PerElemTest::_mesh
std::unique_ptr< Mesh > _mesh
Definition: elem_test.h:22

◆ test_contains_point_node()

template<ElemType elem_type>
void ElemTest< elem_type >::test_contains_point_node ( )
inline

Definition at line 226 of file elem_test.C.

References libMesh::invalid_uint, and libMesh::TOLERANCE.

227  {
228  LOG_UNIT_TEST;
229 
230  for (const auto & elem :
231  this->_mesh->active_local_element_ptr_range())
232  {
233  if (elem->infinite())
234  continue;
235 
236  for (const auto n : elem->node_index_range())
237 #ifndef LIBMESH_ENABLE_EXCEPTIONS
238  // If this node has a singular Jacobian, we need exceptions in order
239  // to catch the failed inverse_map solve and return the singular
240  // master point. Therefore, if we don't have exceptions and we're
241  // at a singular node, we can't test this. As of the writing of
242  // this comment, this issue exists for only Pyramid elements at
243  // the apex.
244  if (elem->local_singular_node(elem->point(n), TOLERANCE*TOLERANCE) == invalid_uint)
245 #endif
246  CPPUNIT_ASSERT(elem->contains_point(elem->point(n)));
247  }
248  }
libMesh::invalid_uint
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value for an unsigned int...
Definition: libmesh.h:310
libMesh::TOLERANCE
static constexpr Real TOLERANCE
Definition: libmesh_common.h:151
PerElemTest::_mesh
std::unique_ptr< Mesh > _mesh
Definition: elem_test.h:22

◆ test_double_refinement()

template<ElemType elem_type>
void ElemTest< elem_type >::test_double_refinement ( )
inline

Definition at line 785 of file elem_test.C.

786  {
787  LOG_UNIT_TEST;
788 
789  test_n_refinements(2);
790  }
ElemTest::test_n_refinements
void test_n_refinements(unsigned int n)
Definition: elem_test.C:595

◆ test_elem_side_builder()

template<ElemType elem_type>
void ElemTest< elem_type >::test_elem_side_builder ( )
inline

Definition at line 573 of file elem_test.C.

574  {
575  LOG_UNIT_TEST;
576 
577  ElemSideBuilder cache;
578  for (auto & elem : this->_mesh->active_local_element_ptr_range())
579  for (const auto s : elem->side_index_range())
580  {
581  const auto side = elem->build_side_ptr(s);
582 
583  auto & cached_side = cache(*elem, s);
584  CPPUNIT_ASSERT_EQUAL(side->type(), cached_side.type());
585  for (const auto n : side->node_index_range())
586  CPPUNIT_ASSERT_EQUAL(side->node_ref(n), cached_side.node_ref(n));
587 
588  const auto & const_cached_side = cache(const_cast<const Elem &>(*elem), s);
589  CPPUNIT_ASSERT_EQUAL(side->type(), const_cached_side.type());
590  for (const auto n : side->node_index_range())
591  CPPUNIT_ASSERT_EQUAL(side->node_ref(n), const_cached_side.node_ref(n));
592  }
593  }
PerElemTest::_mesh
std::unique_ptr< Mesh > _mesh
Definition: elem_test.h:22
libMesh::ElemSideBuilder
Helper for building element sides that minimizes the construction of new elements.
Definition: elem_side_builder.h:50

◆ test_flip()

template<ElemType elem_type>
void ElemTest< elem_type >::test_flip ( )
inline

Definition at line 289 of file elem_test.C.

References libMesh::BoundaryInfo::boundary_ids(), libMesh::C0POLYGON, libMesh::C0POLYHEDRON, libMesh::invalid_uint, libMesh::make_range(), and libMesh::TOLERANCE.

290  {
291  LOG_UNIT_TEST;
292 
293  BoundaryInfo & boundary_info = this->_mesh->get_boundary_info();
294 
295  for (const auto & elem :
296  this->_mesh->active_local_element_ptr_range())
297  {
298  if (elem->infinite())
299  continue;
300 
301  if (elem->type() == C0POLYHEDRON)
302  continue;
303 
304  const Point vertex_avg = elem->vertex_average();
305 
306  const unsigned int n_sides = elem->n_sides();
307  std::vector<std::set<Point*>> side_nodes(n_sides);
308  std::vector<Elem*> neighbors(n_sides);
309  std::vector<std::vector<boundary_id_type>> bcids(n_sides);
310  for (auto s : make_range(n_sides))
311  {
312  for (auto n : elem->nodes_on_side(s))
313  side_nodes[s].insert(elem->node_ptr(n));
314  neighbors[s] = elem->neighbor_ptr(s);
315  boundary_info.boundary_ids(elem, s, bcids[s]);
316  }
317 
318  elem->flip(&boundary_info);
319 
320  // We should just be flipped, not twisted, so our map should
321  // still be affine.
322  // ... except for stupid singular pyramid maps
323  // ... or the polygons we're deliberately testing non-affine
324  if ((elem->dim() < 3 ||
325  elem->n_vertices() != 5) &&
326  elem_type != C0POLYGON)
327  CPPUNIT_ASSERT(elem->has_affine_map());
328  else if (elem_type == C0POLYGON)
329  CPPUNIT_ASSERT(!elem->has_affine_map());
330 
331  // The neighbors and bcids should have flipped in a way
332  // consistently with the nodes (unless this is a 0D NodeElem)
333  bool something_changed = false;
334  for (auto s : make_range(n_sides))
335  {
336  std::set<Point*> new_side_nodes;
337  for (auto n : elem->nodes_on_side(s))
338  new_side_nodes.insert(elem->node_ptr(n));
339 
340  std::vector<boundary_id_type> new_bcids;
341  boundary_info.boundary_ids(elem, s, new_bcids);
342 
343  unsigned int old_side = libMesh::invalid_uint;
344  for (auto os : make_range(n_sides))
345  if (new_side_nodes == side_nodes[os])
346  old_side = os;
347 
348  if (old_side != s)
349  something_changed = true;
350 
351  CPPUNIT_ASSERT(old_side != libMesh::invalid_uint);
352 
353  CPPUNIT_ASSERT(neighbors[old_side] ==
354  elem->neighbor_ptr(s));
355 
356  CPPUNIT_ASSERT(bcids[old_side] == new_bcids);
357  }
358  CPPUNIT_ASSERT(!elem->dim() || something_changed);
359 
360  const Point new_vertex_avg = elem->vertex_average();
361  for (const auto d : make_range(LIBMESH_DIM))
362  LIBMESH_ASSERT_FP_EQUAL(vertex_avg(d), new_vertex_avg(d),
363  TOLERANCE*TOLERANCE);
364  }
365  }
libMesh::invalid_uint
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value for an unsigned int...
Definition: libmesh.h:310
libMesh::TOLERANCE
static constexpr Real TOLERANCE
Definition: libmesh_common.h:151
PerElemTest::_mesh
std::unique_ptr< Mesh > _mesh
Definition: elem_test.h:22
libMesh::BoundaryInfo::boundary_ids
void boundary_ids(const Node *node, std::vector< boundary_id_type > &vec_to_fill) const
Fills a user-provided std::vector with the boundary ids associated with Node node.
Definition: boundary_info.C:1333
libMesh::BoundaryInfo
The BoundaryInfo class contains information relevant to boundary conditions including storing faces...
Definition: boundary_info.h:57
libMesh::make_range
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
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39

◆ test_is_internal()

template<ElemType elem_type>
void ElemTest< elem_type >::test_is_internal ( )
inline

Definition at line 792 of file elem_test.C.

References libMesh::EDGE3, libMesh::EDGE4, libMesh::HEX27, libMesh::INFHEX18, libMesh::INFQUAD6, libMesh::PRISM21, libMesh::QUAD9, libMesh::QUADSHELL9, and libMesh::TRI7.

793  {
794  LOG_UNIT_TEST;
795 
796  for (const auto & elem :
797  this->_mesh->active_local_element_ptr_range())
798  for (const auto nd : elem->node_index_range())
799  {
800  if ((elem->type() == EDGE3 || elem->type() == EDGE4) && nd >= 2)
801  CPPUNIT_ASSERT(elem->is_internal(nd));
802  else if (elem->type() == HEX27 && nd == 26)
803  CPPUNIT_ASSERT(elem->is_internal(nd));
804  else if (elem->type() == PRISM21 && nd == 20)
805  CPPUNIT_ASSERT(elem->is_internal(nd));
806  else if ((elem->type() == QUAD9 || elem->type() == QUADSHELL9) && nd == 8)
807  CPPUNIT_ASSERT(elem->is_internal(nd));
808  else if (elem->type() == TRI7 && nd == 6)
809  CPPUNIT_ASSERT(elem->is_internal(nd));
810  else if (elem->type() == INFHEX18 && nd == 17)
811  CPPUNIT_ASSERT(elem->is_internal(nd));
812  else if (elem->type() == INFQUAD6 && nd == 5)
813  CPPUNIT_ASSERT(elem->is_internal(nd));
814  else
815  CPPUNIT_ASSERT(!elem->is_internal(nd));
816  }
817  }
PerElemTest::_mesh
std::unique_ptr< Mesh > _mesh
Definition: elem_test.h:22

◆ test_maps()

template<ElemType elem_type>
void ElemTest< elem_type >::test_maps ( )
inline

Definition at line 166 of file elem_test.C.

167  {
168  LOG_UNIT_TEST;
169 
170  for (const auto & elem :
171  this->_mesh->active_local_element_ptr_range())
172  {
173  for (const auto edge : elem->edge_index_range())
174  for (const auto side_on_edge : elem->sides_on_edge(edge))
175  for (const auto node_on_edge : elem->nodes_on_edge(edge))
176  CPPUNIT_ASSERT(elem->is_node_on_side(node_on_edge, side_on_edge));
177 
178  for (const auto side : elem->side_index_range())
179  for (const auto node_on_side : elem->nodes_on_side(side))
180  CPPUNIT_ASSERT(elem->is_node_on_side(node_on_side, side));
181 
182  for (const auto edge : elem->edge_index_range())
183  for (const auto node_on_edge : elem->nodes_on_edge(edge))
184  CPPUNIT_ASSERT(elem->is_node_on_edge(node_on_edge, edge));
185 
186  for (const auto edge : elem->edge_index_range())
187  for (const auto side_on_edge : elem->sides_on_edge(edge))
188  CPPUNIT_ASSERT(elem->is_edge_on_side(edge, side_on_edge));
189  }
190  }
PerElemTest::_mesh
std::unique_ptr< Mesh > _mesh
Definition: elem_test.h:22

◆ test_n_refinements()

template<ElemType elem_type>
void ElemTest< elem_type >::test_n_refinements ( unsigned int  n)
inline

Definition at line 595 of file elem_test.C.

References libMesh::C0POLYGON, libMesh::C0POLYHEDRON, libMesh::EDGE4, libMesh::invalid_uint, libMesh::make_range(), libMesh::Elem::parent(), libMesh::Utility::pow(), libMesh::PRISM20, libMesh::PYRAMID13, libMesh::PYRAMID14, libMesh::PYRAMID18, libMesh::PYRAMID5, TIMPI::Communicator::set_union(), TestCommWorld, and libMesh::MeshRefinement::uniformly_refine().

596  {
597 #ifdef LIBMESH_ENABLE_AMR
598  // We don't support refinement of all element types
599  if (elem_type == EDGE4 ||
600  elem_type == PRISM20 ||
601  elem_type == PYRAMID5 ||
602  elem_type == PYRAMID13 ||
603  elem_type == PYRAMID14 ||
604  elem_type == PYRAMID18 ||
605  elem_type == C0POLYGON ||
606  elem_type == C0POLYHEDRON)
607  return;
608 
609  auto refining_mesh = this->_mesh->clone();
610 
611  MeshRefinement mr(*refining_mesh);
612  mr.uniformly_refine(n);
613 
614  std::set<std::pair<dof_id_type, unsigned int>> parent_node_was_touched;
615  std::set<std::pair<dof_id_type, unsigned int>> parent_child_was_touched;
616 
617  for (const Elem * elem : refining_mesh->active_element_ptr_range())
618  {
619  CPPUNIT_ASSERT_EQUAL(elem->level(), n);
620  CPPUNIT_ASSERT(!elem->ancestor());
621  CPPUNIT_ASSERT(elem->active());
622  CPPUNIT_ASSERT(!elem->subactive());
623  CPPUNIT_ASSERT(!elem->has_children());
624  CPPUNIT_ASSERT(!elem->has_ancestor_children());
625  CPPUNIT_ASSERT(!elem->interior_parent());
626 
627  const Elem * parent = elem->parent();
628  CPPUNIT_ASSERT(parent);
629  CPPUNIT_ASSERT(parent->ancestor());
630  CPPUNIT_ASSERT(!parent->active());
631  CPPUNIT_ASSERT(!parent->subactive());
632  CPPUNIT_ASSERT(parent->has_children());
633  CPPUNIT_ASSERT(!parent->has_ancestor_children());
634  CPPUNIT_ASSERT(!parent->interior_parent());
635  if (n == 1)
636  {
637  CPPUNIT_ASSERT_EQUAL(parent, elem->top_parent());
638  CPPUNIT_ASSERT_EQUAL(parent, parent->top_parent());
639  }
640  else
641  {
642  CPPUNIT_ASSERT(parent != elem->top_parent());
643  CPPUNIT_ASSERT(parent != parent->top_parent());
644  CPPUNIT_ASSERT_EQUAL(elem->top_parent(), parent->top_parent());
645  }
646 
647  CPPUNIT_ASSERT(parent->is_ancestor_of(elem));
648  const unsigned int c = parent->which_child_am_i(elem);
649  CPPUNIT_ASSERT(c < parent->n_children());
650  CPPUNIT_ASSERT_EQUAL(elem, parent->child_ptr(c));
651  parent_child_was_touched.emplace(parent->id(), c);
652 
653  CPPUNIT_ASSERT_EQUAL(parent->n_nodes_in_child(c), elem->n_nodes());
654  for (auto n : make_range(elem->n_nodes()))
655  {
656  CPPUNIT_ASSERT_EQUAL(parent->is_vertex_on_child(c, n), elem->is_vertex(n));
657 
658  auto pn = parent->as_parent_node(c, n);
659  CPPUNIT_ASSERT_EQUAL(pn, parent->get_node_index(elem->node_ptr(n)));
660  if (pn == libMesh::invalid_uint)
661  continue;
662  CPPUNIT_ASSERT_EQUAL(parent->is_vertex_on_parent(c, n), parent->is_vertex(pn));
663  parent_node_was_touched.emplace(parent->id(), pn);
664  }
665 
666  for (auto s : make_range(parent->n_sides()))
667  {
668  if (parent->is_child_on_side(c,s))
669  {
670  auto parent_side = parent->build_side_ptr(s);
671 
672  // Implicitly assuming here that s is the child side
673  // too - we support that now and hopefully won't have
674  // to change it later
675  auto child_side = elem->build_side_ptr(s);
676 
677  // 2D Inf FE inverse_map not yet implemented?
678  if (!parent_side->infinite())
679  for (const Node & node : child_side->node_ref_range())
680  CPPUNIT_ASSERT(parent_side->contains_point(node));
681 
682  if (elem->neighbor_ptr(s) && !elem->neighbor_ptr(s)->is_remote())
683  CPPUNIT_ASSERT_EQUAL(parent->child_neighbor(elem->neighbor_ptr(s)), elem);
684  }
685  }
686 
687  for (auto e : make_range(parent->n_edges()))
688  {
689  if (parent->is_child_on_edge(c,e))
690  {
691  auto parent_edge = parent->build_edge_ptr(e);
692 
693  // Implicitly assuming here that e is the child edge
694  // too - we support that now and hopefully won't have
695  // to change it later
696  auto child_edge = elem->build_edge_ptr(e);
697 
698  // 1D Inf FE inverse_map not yet implemented?
699  if (!parent_edge->infinite())
700  for (const Node & node : child_edge->node_ref_range())
701  CPPUNIT_ASSERT(parent_edge->contains_point(node));
702  }
703  }
704 
705  if (parent->has_affine_map())
706  CPPUNIT_ASSERT(elem->has_affine_map());
707  }
708 
709  // It's possible for a parent element on a distributed mesh to not
710  // have all its children available on any one processor
711  TestCommWorld->set_union(parent_child_was_touched);
712  TestCommWorld->set_union(parent_node_was_touched);
713 
714  for (const Elem * elem : refining_mesh->local_element_ptr_range())
715  {
716  if (elem->active())
717  continue;
718 
719  // With only one layer of refinement the family tree methods
720  // should have the full number of elements, even if some are
721  // remote.
722  if (n == 1)
723  {
724  std::vector<const Elem *> family;
725  elem->family_tree(family);
726  CPPUNIT_ASSERT_EQUAL(family.size(),
727  std::size_t(elem->n_children() + 1));
728 
729  family.clear();
730  elem->total_family_tree(family);
731  CPPUNIT_ASSERT_EQUAL(family.size(),
732  std::size_t(elem->n_children() + 1));
733 
734  family.clear();
735  elem->active_family_tree(family);
736  CPPUNIT_ASSERT_EQUAL(family.size(),
737  std::size_t(elem->n_children()));
738 
739  for (auto s : make_range(elem->n_sides()))
740  {
741  family.clear();
742  elem->active_family_tree_by_side(family,s);
743  if (!elem->build_side_ptr(s)->infinite())
744  CPPUNIT_ASSERT_EQUAL(double(family.size()),
745  std::pow(2.0, int(elem->dim()-1)));
746  else
747  CPPUNIT_ASSERT_EQUAL(double(family.size()),
748  std::pow(2.0, int(elem->dim()-2)));
749  for (const Elem * child : family)
750  {
751  if (child->is_remote())
752  continue;
753 
754  unsigned int c = elem->which_child_am_i(child);
755  CPPUNIT_ASSERT(elem->is_child_on_side(c, s));
756  }
757  }
758  }
759 
760  if (elem->level() + 1 == n)
761  {
762  for (auto c : make_range(elem->n_children()))
763  {
764  auto it = parent_child_was_touched.find(std::make_pair(elem->id(), c));
765  CPPUNIT_ASSERT(it != parent_child_was_touched.end());
766  }
767 
768  for (auto n : make_range(elem->n_nodes()))
769  {
770  auto it = parent_node_was_touched.find(std::make_pair(elem->id(), n));
771  CPPUNIT_ASSERT(it != parent_node_was_touched.end());
772  }
773  }
774  }
775 #endif
776  }
libMesh::Elem::parent
const Elem * parent() const
Definition: elem.h:3030
libMesh::Node
A Node is like a Point, but with more information.
Definition: node.h:52
libMesh::invalid_uint
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value for an unsigned int...
Definition: libmesh.h:310
TestCommWorld
libMesh::Parallel::Communicator * TestCommWorld
Definition: driver.C:171
PerElemTest::_mesh
std::unique_ptr< Mesh > _mesh
Definition: elem_test.h:22
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
libMesh::Utility::pow
T pow(const T &x)
Definition: utility.h:328
libMesh::MeshRefinement
Implements (adaptive) mesh refinement algorithms for a MeshBase.
Definition: mesh_refinement.h:61
libMesh::make_range
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
TIMPI::Communicator::set_union
void set_union(T &data, const unsigned int root_id) const

◆ test_node_edge_map_consistency()

template<ElemType elem_type>
void ElemTest< elem_type >::test_node_edge_map_consistency ( )
inline

Definition at line 819 of file elem_test.C.

820  {
821  LOG_UNIT_TEST;
822 
823  for (const auto & elem : this->_mesh->active_local_element_ptr_range())
824  {
825  for (const auto nd : elem->node_index_range())
826  {
827  auto adjacent_edge_ids = elem->edges_adjacent_to_node(nd);
828 
829  if (elem->dim() < 2)
830  {
831  // 0D elements don't have edges.
832  // 1D elements *are* "edges", but we don't consider
833  // them to *have* edges, so assert that here.
834  CPPUNIT_ASSERT(adjacent_edge_ids.empty());
835  }
836  else
837  {
838  // For 2D and 3D elements, on each edge which is
839  // claimed to be adjacent, check that "nd" is indeed
840  // on it
841  for (const auto & edge_id : adjacent_edge_ids)
842  {
843  auto node_ids_on_edge = elem->nodes_on_edge(edge_id);
844  CPPUNIT_ASSERT(std::find(node_ids_on_edge.begin(), node_ids_on_edge.end(), nd) != node_ids_on_edge.end());
845  }
846  }
847  }
848  }
849  }
PerElemTest::_mesh
std::unique_ptr< Mesh > _mesh
Definition: elem_test.h:22

◆ test_orient()

template<ElemType elem_type>
void ElemTest< elem_type >::test_orient ( )
inline

Definition at line 367 of file elem_test.C.

References libMesh::BoundaryInfo::boundary_ids(), libMesh::C0POLYGON, libMesh::C0POLYHEDRON, libMesh::make_range(), and libMesh::TOLERANCE.

368  {
369  LOG_UNIT_TEST;
370 
371  BoundaryInfo & boundary_info = this->_mesh->get_boundary_info();
372 
373  for (const auto & elem :
374  this->_mesh->active_local_element_ptr_range())
375  {
376  if (elem->infinite())
377  continue;
378 
379  const Point vertex_avg = elem->vertex_average();
380 
381  const unsigned int n_sides = elem->n_sides();
382  std::vector<std::set<Point*>> side_nodes(n_sides);
383  std::vector<Elem*> neighbors(n_sides);
384  std::vector<std::vector<boundary_id_type>> bcids(n_sides);
385  for (auto s : make_range(n_sides))
386  {
387  for (auto n : elem->nodes_on_side(s))
388  side_nodes[s].insert(elem->node_ptr(n));
389  neighbors[s] = elem->neighbor_ptr(s);
390  boundary_info.boundary_ids(elem, s, bcids[s]);
391  }
392 
393  CPPUNIT_ASSERT(!elem->is_flipped());
394 
395  if (elem->id()%2)
396  {
397  elem->flip(&boundary_info);
398  CPPUNIT_ASSERT(elem->is_flipped());
399  }
400 
401  elem->orient(&boundary_info);
402  CPPUNIT_ASSERT(!elem->is_flipped());
403 
404  // Our map should still be affine.
405  // ... except for stupid singular pyramid maps
406  // ... or the polygons we're deliberately testing non-affine
407  // ... or the polyhedra we deliberately don't define "affine
408  // map" for.
409  if ((elem->dim() < 3 ||
410  elem->n_vertices() != 5) &&
411  elem_type != C0POLYGON &&
412  elem_type != C0POLYHEDRON)
413  CPPUNIT_ASSERT(elem->has_affine_map());
414  else if (elem_type == C0POLYGON &&
415  elem_type != C0POLYHEDRON)
416  CPPUNIT_ASSERT(!elem->has_affine_map());
417 
418  // The neighbors and bcids should have flipped back to where
419  // they were.
420  for (auto s : make_range(n_sides))
421  {
422  std::set<Point*> new_side_nodes;
423  for (auto n : elem->nodes_on_side(s))
424  new_side_nodes.insert(elem->node_ptr(n));
425 
426  std::vector<boundary_id_type> new_bcids;
427  boundary_info.boundary_ids(elem, s, new_bcids);
428 
429  CPPUNIT_ASSERT(side_nodes[s] ==
430  new_side_nodes);
431 
432  CPPUNIT_ASSERT(neighbors[s] ==
433  elem->neighbor_ptr(s));
434 
435  CPPUNIT_ASSERT(bcids[s] == new_bcids);
436  }
437 
438  const Point new_vertex_avg = elem->vertex_average();
439  for (const auto d : make_range(LIBMESH_DIM))
440  LIBMESH_ASSERT_FP_EQUAL(vertex_avg(d), new_vertex_avg(d),
441  TOLERANCE*TOLERANCE);
442  }
443  }
libMesh::TOLERANCE
static constexpr Real TOLERANCE
Definition: libmesh_common.h:151
PerElemTest::_mesh
std::unique_ptr< Mesh > _mesh
Definition: elem_test.h:22
libMesh::BoundaryInfo::boundary_ids
void boundary_ids(const Node *node, std::vector< boundary_id_type > &vec_to_fill) const
Fills a user-provided std::vector with the boundary ids associated with Node node.
Definition: boundary_info.C:1333
libMesh::BoundaryInfo
The BoundaryInfo class contains information relevant to boundary conditions including storing faces...
Definition: boundary_info.h:57
libMesh::make_range
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
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39

◆ test_orient_elements()

template<ElemType elem_type>
void ElemTest< elem_type >::test_orient_elements ( )
inline

Definition at line 445 of file elem_test.C.

References libMesh::BoundaryInfo::boundary_ids(), libMesh::DofObject::id(), libMesh::make_range(), libMesh::Elem::neighbor_ptr(), and libMesh::MeshTools::Modification::orient_elements().

446  {
447  LOG_UNIT_TEST;
448 
449  const Mesh old_mesh {*this->_mesh};
450 
451  BoundaryInfo & boundary_info = this->_mesh->get_boundary_info();
452  const BoundaryInfo & old_boundary_info = old_mesh.get_boundary_info();
453  CPPUNIT_ASSERT(&boundary_info != &old_boundary_info);
454 
455  for (const auto & elem :
456  this->_mesh->active_local_element_ptr_range())
457  {
458  if (elem->infinite())
459  continue;
460 
461  if (elem->id()%2)
462  {
463  elem->flip(&boundary_info);
464  CPPUNIT_ASSERT(elem->is_flipped());
465  }
466  }
467 
468  MeshTools::Modification::orient_elements(*this->_mesh);
469 
470  // I should really create a MeshBase::operator==()...
471  for (const auto & elem :
472  this->_mesh->active_local_element_ptr_range())
473  {
474  const Elem & old_elem = old_mesh.elem_ref(elem->id());
475 
476  CPPUNIT_ASSERT(!elem->is_flipped());
477 
478  // Elem::operator==() uses node ids to compare
479  CPPUNIT_ASSERT(*elem == old_elem);
480 
481  const unsigned int n_sides = elem->n_sides();
482  for (auto s : make_range(n_sides))
483  {
484  std::vector<boundary_id_type> bcids, old_bcids;
485  boundary_info.boundary_ids(elem, s, bcids);
486  old_boundary_info.boundary_ids(&old_elem, s, old_bcids);
487  CPPUNIT_ASSERT(bcids == old_bcids);
488 
489  if (elem->neighbor_ptr(s))
490  {
491  CPPUNIT_ASSERT(old_elem.neighbor_ptr(s));
492  CPPUNIT_ASSERT_EQUAL(elem->neighbor_ptr(s)->id(),
493  old_elem.neighbor_ptr(s)->id());
494  }
495  else
496  CPPUNIT_ASSERT(!old_elem.neighbor_ptr(s));
497  }
498  }
499  }
PerElemTest::_mesh
std::unique_ptr< Mesh > _mesh
Definition: elem_test.h:22
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
libMesh::BoundaryInfo::boundary_ids
void boundary_ids(const Node *node, std::vector< boundary_id_type > &vec_to_fill) const
Fills a user-provided std::vector with the boundary ids associated with Node node.
Definition: boundary_info.C:1333
libMesh::DofObject::id
dof_id_type id() const
Definition: dof_object.h:828
libMesh::BoundaryInfo
The BoundaryInfo class contains information relevant to boundary conditions including storing faces...
Definition: boundary_info.h:57
libMesh::MeshTools::Modification::orient_elements
void orient_elements(MeshBase &mesh)
Redo the nodal ordering of each element as necessary to give the element Jacobian a positive orientat...
Definition: mesh_modification.C:262
libMesh::Elem::neighbor_ptr
const Elem * neighbor_ptr(unsigned int i) const
Definition: elem.h:2598
libMesh::make_range
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
libMesh::Mesh
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition: mesh.h:50

◆ test_permute()

template<ElemType elem_type>
void ElemTest< elem_type >::test_permute ( )
inline

Definition at line 250 of file elem_test.C.

References libMesh::make_range(), and libMesh::TOLERANCE.

251  {
252  LOG_UNIT_TEST;
253 
254  for (const auto & elem :
255  this->_mesh->active_local_element_ptr_range())
256  {
257  if (elem->infinite())
258  continue;
259 
260  const Point centroid = elem->true_centroid();
261  const Point vertex_avg = elem->vertex_average();
262  Point quasicc;
263  if (elem->dim() < 3)
264  quasicc = elem->quasicircumcenter();
265 
266  for (const auto p : IntRange<unsigned int>(0, elem->n_permutations()))
267  {
268  elem->permute(p);
269  CPPUNIT_ASSERT(elem->has_invertible_map());
270  const Point new_centroid = elem->true_centroid();
271  const Point new_vertex_avg = elem->vertex_average();
272  Point new_quasicc;
273  if (elem->dim() < 3)
274  new_quasicc = elem->quasicircumcenter();
275  for (const auto d : make_range(LIBMESH_DIM))
276  {
277  // Getting a little FP error from Pyramid18
278  LIBMESH_ASSERT_FP_EQUAL(centroid(d), new_centroid(d),
279  TOLERANCE*std::sqrt(TOLERANCE));
280  LIBMESH_ASSERT_FP_EQUAL(vertex_avg(d), new_vertex_avg(d),
281  TOLERANCE*TOLERANCE);
282  LIBMESH_ASSERT_FP_EQUAL(quasicc(d), new_quasicc(d),
283  TOLERANCE*TOLERANCE);
284  }
285  }
286  }
287  }
libMesh::IntRange
The IntRange templated class is intended to make it easy to loop over integers which are indices of a...
Definition: int_range.h:53
libMesh::TOLERANCE
static constexpr Real TOLERANCE
Definition: libmesh_common.h:151
PerElemTest::_mesh
std::unique_ptr< Mesh > _mesh
Definition: elem_test.h:22
libMesh::make_range
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
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39

◆ test_quality()

template<ElemType elem_type>
void ElemTest< elem_type >::test_quality ( )
inline

Definition at line 52 of file elem_test.C.

References libMesh::C0POLYGON, libMesh::C0POLYHEDRON, libMesh::EDGE_LENGTH_RATIO, libMesh::JACOBIAN, libMesh::MAX_ANGLE, libMesh::MAX_DIHEDRAL_ANGLE, libMesh::MIN_ANGLE, libMesh::MIN_DIHEDRAL_ANGLE, libMesh::pi, libMesh::Real, libMesh::SCALED_JACOBIAN, and libMesh::TOLERANCE.

53  {
54  LOG_UNIT_TEST;
55 
56  for (const auto & elem : this->_mesh->active_local_element_ptr_range())
57  {
58  // EDGE_LENGTH_RATIO is one metric that is defined on all elements
59  const Real edge_length_ratio = elem->quality(EDGE_LENGTH_RATIO);
60 
61  // We use "0" to mean infinity rather than inf or NaN, and
62  // every quality other than that should be 1 or larger (worse)
63  CPPUNIT_ASSERT_LESSEQUAL(edge_length_ratio, Real(1)); // 1 <= edge_length_ratio
64 
65  // We're building isotropic meshes, where even elements
66  // dissected from cubes ought to have tolerable quality.
67  //
68  // Worst I see is 2 on tets, but let's add a little tolerance
69  // in case we decide to play with rotated meshes here later
70  CPPUNIT_ASSERT_LESSEQUAL(Real(2+TOLERANCE), edge_length_ratio); // edge_length_ratio <= 2
71 
72  // The MIN_ANGLE and MAX_ANGLE quality metrics are also defined on all elements
73  const Real min_angle = elem->quality(MIN_ANGLE);
74  const Real max_angle = elem->quality(MAX_ANGLE);
75 
76  // Reference Quads/Hexes have maximum internal angles of 90
77  // degrees, but pyramids actually have an obtuse interior
78  // angle of acos(-1/3) ~ 109.47 deg formed by edge pairs
79  // {(0, 4), (2, 4)} and {(1,4), (3,4)}.
80  if (elem->type() != C0POLYGON)
81  CPPUNIT_ASSERT_LESSEQUAL((std::acos(Real(-1)/3) * 180 / libMesh::pi) + TOLERANCE, max_angle);
82 
83  // We're doing some polygon testing with a squashed pentagon
84  // that has a 135 degree angle
85  CPPUNIT_ASSERT_LESSEQUAL(135 + TOLERANCE, max_angle);
86 
87  // Notes on minimum angle we expect to see:
88  // 1.) 1D Elements don't have interior angles, so the base
89  // class implementation currently returns 0 for those
90  // elements.
91  // 2.) Reference triangles/tetrahedra have min interior angle
92  // of 45 degrees, however, here we are checking Tets in a
93  // build_cube() mesh which have angles as small as
94  // acos(2/sqrt(6)) so we use that as our lower bound here.
95  if (elem->dim() > 1)
96  CPPUNIT_ASSERT_GREATEREQUAL((std::acos(Real(2)/std::sqrt(Real(6))) * 180 / libMesh::pi) - TOLERANCE, min_angle);
97 
98  // MIN,MAX_DIHEDRAL_ANGLE are implemented for all 3D elements
99  if (elem->dim() > 2)
100  {
101  const Real min_dihedral_angle = elem->quality(MIN_DIHEDRAL_ANGLE);
102  const Real max_dihedral_angle = elem->quality(MAX_DIHEDRAL_ANGLE);
103 
104  // Debugging
105  // libMesh::out << "Elem type: " << Utility::enum_to_string(elem->type())
106  // << ", min_dihedral_angle = " << min_dihedral_angle
107  // << ", max_dihedral_angle = " << max_dihedral_angle
108  // << std::endl;
109 
110  // Assert that we match expected values for build_cube() meshes.
111  // * build_cube() meshes of hexes, tetrahedra, and prisms have max_dihedral_angle == 90
112  // * build_cube() meshes of pyramids have max_dihedral_angle == 60 between adjacent triangular faces
113  CPPUNIT_ASSERT_LESSEQUAL (90 + TOLERANCE, max_dihedral_angle);
114 
115  // * build_cube() meshes of hexes have min_dihedral_angle == 90
116  // * build_cube() meshes of prisms, tets, and pyramids have min_dihedral_angle == 45
117  // * For the InfPrism tests, we construct a single
118  // InfPrism by hand with interior angle ~53.13 deg. so that
119  // is the minimum dihedral angle we expect in that case.
120  CPPUNIT_ASSERT_GREATEREQUAL(45 - TOLERANCE, min_dihedral_angle);
121  }
122 
123  // The JACOBIAN and SCALED_JACOBIAN quality metrics are also defined on all elements
124 
125  // The largest non-infinite elements in this test (Hexes and
126  // Prisms) have a nodal Jacobian (volume) of 8, e.g. the 2x2x2
127  // reference Hex. The infinite elements that we construct for
128  // this testing have a max nodal area of 64 since they are
129  // created (see setUp()) with max side lengths of 4 (4*4*4=64).
130  const Real jac = elem->quality(JACOBIAN);
131  if (elem->infinite())
132  CPPUNIT_ASSERT_LESSEQUAL (64 + TOLERANCE, jac);
133  else
134  CPPUNIT_ASSERT_LESSEQUAL (8 + TOLERANCE, jac);
135 
136  // The smallest 2D/3D nodal areas for regular elements in this
137  // test are found in tetrahedra, which have a minimum value of
138  // 2.0. However, we return a default value of 1.0 for 0D and
139  // 1D elements here, and we have a custom distorted-pentagon
140  // C0Polygon with a 0.5 at 2 nodes and a custom C0Polyhedron
141  // with a 1, so we handle those cases too.
142  if (elem->dim() < 2 || elem->type() == C0POLYHEDRON)
143  CPPUNIT_ASSERT_GREATEREQUAL(1 - TOLERANCE, jac);
144  else if (elem->type() == C0POLYGON)
145  CPPUNIT_ASSERT_GREATEREQUAL(0.5 - TOLERANCE, jac);
146  else
147  CPPUNIT_ASSERT_GREATEREQUAL(2 - TOLERANCE, jac);
148 
149  // The scaled Jacobian should always be <= 1. The minimum
150  // scaled Jacobian value I observed was ~0.408248 for the
151  // tetrahedral meshes. This is consistent with the way that
152  // we generate Tet meshes with build_cube(), as we don't
153  // simply refine the reference tetrahedron in that case.
154  const Real scaled_jac = elem->quality(SCALED_JACOBIAN);
155  CPPUNIT_ASSERT_LESSEQUAL (1 + TOLERANCE, scaled_jac);
156  CPPUNIT_ASSERT_GREATEREQUAL( Real(0.4), scaled_jac);
157 
158  // Debugging
159  // libMesh::out << "elem->type() = " << Utility::enum_to_string(elem->type())
160  // << ", jac = " << jac
161  // << ", scaled_jac = " << scaled_jac
162  // << std::endl;
163  }
164  }
libMesh::TOLERANCE
static constexpr Real TOLERANCE
Definition: libmesh_common.h:151
PerElemTest::_mesh
std::unique_ptr< Mesh > _mesh
Definition: elem_test.h:22
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:144
libMesh::pi
const Real pi
.
Definition: libmesh.h:299

◆ test_refinement()

template<ElemType elem_type>
void ElemTest< elem_type >::test_refinement ( )
inline

Definition at line 778 of file elem_test.C.

779  {
780  LOG_UNIT_TEST;
781 
782  test_n_refinements(1);
783  }
ElemTest::test_n_refinements
void test_n_refinements(unsigned int n)
Definition: elem_test.C:595

◆ test_side_subdomain()

template<ElemType elem_type>
void ElemTest< elem_type >::test_side_subdomain ( )
inline

Definition at line 550 of file elem_test.C.

551  {
552  LOG_UNIT_TEST;
553 
554  for (const auto & elem :
555  this->_mesh->active_local_element_ptr_range())
556  {
557  std::unique_ptr<const Elem> side;
558 
559  for (const auto s : elem->side_index_range())
560  {
561  CPPUNIT_ASSERT_EQUAL(elem->build_side_ptr(s)->subdomain_id(), elem->subdomain_id());
562 
563  elem->build_side_ptr(side, s);
564  CPPUNIT_ASSERT_EQUAL(side->subdomain_id(), elem->subdomain_id());
565 
566  // We don't require that the "lightweight" side_ptr work
567  // the same way
568  // CPPUNIT_ASSERT_EQUAL(elem->side_ptr(s)->subdomain_id(), elem->subdomain_id());
569  }
570  }
571  }
PerElemTest::_mesh
std::unique_ptr< Mesh > _mesh
Definition: elem_test.h:22

◆ test_side_type()

template<ElemType elem_type>
void ElemTest< elem_type >::test_side_type ( )
inline

Definition at line 540 of file elem_test.C.

541  {
542  LOG_UNIT_TEST;
543 
544  for (const auto & elem :
545  this->_mesh->active_local_element_ptr_range())
546  for (const auto s : elem->side_index_range())
547  CPPUNIT_ASSERT_EQUAL(elem->build_side_ptr(s)->type(), elem->side_type(s));
548  }
PerElemTest::_mesh
std::unique_ptr< Mesh > _mesh
Definition: elem_test.h:22

◆ test_static_data()

template<ElemType elem_type>
void ElemTest< elem_type >::test_static_data ( )
inline

Definition at line 192 of file elem_test.C.

References libMesh::invalid_uint, libMesh::Elem::max_n_nodes, libMesh::Elem::type_to_default_order_map, libMesh::Elem::type_to_dim_map, libMesh::Elem::type_to_n_edges_map, libMesh::Elem::type_to_n_nodes_map, and libMesh::Elem::type_to_n_sides_map.

193  {
194  LOG_UNIT_TEST;
195 
196  for (const auto & elem :
197  this->_mesh->active_local_element_ptr_range())
198  {
199  CPPUNIT_ASSERT(elem->n_nodes() <= Elem::max_n_nodes);
200 
201  const ElemType etype = elem->type();
202 
203  CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(elem->dim()),
204  Elem::type_to_dim_map[etype]);
205  CPPUNIT_ASSERT_EQUAL(elem->default_order(),
206  Elem::type_to_default_order_map[etype]);
207 
208  // If we have an element type with topology defined solely by
209  // the type, then that should match the runtime topology. If
210  // not, then we should be aware of it.
211  if (!elem->runtime_topology())
212  {
213  CPPUNIT_ASSERT_EQUAL(elem->n_nodes(), Elem::type_to_n_nodes_map[etype]);
214  CPPUNIT_ASSERT_EQUAL(elem->n_sides(), Elem::type_to_n_sides_map[etype]);
215  CPPUNIT_ASSERT_EQUAL(elem->n_edges(), Elem::type_to_n_edges_map[etype]);
216  }
217  else
218  {
219  CPPUNIT_ASSERT_EQUAL(invalid_uint, Elem::type_to_n_nodes_map[etype]);
220  CPPUNIT_ASSERT_EQUAL(invalid_uint, Elem::type_to_n_sides_map[etype]);
221  CPPUNIT_ASSERT_EQUAL(invalid_uint, Elem::type_to_n_edges_map[etype]);
222  }
223  }
224  }
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition: enum_elem_type.h:33
libMesh::invalid_uint
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value for an unsigned int...
Definition: libmesh.h:310
PerElemTest::_mesh
std::unique_ptr< Mesh > _mesh
Definition: elem_test.h:22

Member Data Documentation

◆ _mesh

template<ElemType elem_type>
std::unique_ptr<Mesh> PerElemTest< elem_type >::_mesh
protectedinherited

Definition at line 22 of file elem_test.h.

◆ libmesh_suite_name

template<ElemType elem_type>
std::string PerElemTest< elem_type >::libmesh_suite_name
protectedinherited

Definition at line 23 of file elem_test.h.

The documentation for this class was generated from the following file:
Generated on Fri Jul 18 2025 13:27:39 for libMesh by   doxygen 1.8.14