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 30 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().

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

◆ 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:26
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:303
PerElemTest::_mesh
std::unique_ptr< Mesh > _mesh
Definition: elem_test.h:26

◆ 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:303
libMesh::TOLERANCE
static constexpr Real TOLERANCE
Definition: libmesh_common.h:151
PerElemTest::_mesh
std::unique_ptr< Mesh > _mesh
Definition: elem_test.h:26

◆ test_double_refinement()

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

Definition at line 809 of file elem_test.C.

810  {
811  LOG_UNIT_TEST;
812 
813  test_n_refinements(2);
814  }
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:26
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:303
libMesh::TOLERANCE
static constexpr Real TOLERANCE
Definition: libmesh_common.h:151
PerElemTest::_mesh
std::unique_ptr< Mesh > _mesh
Definition: elem_test.h:26
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:1348
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:176
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 816 of file elem_test.C.

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

817  {
818  LOG_UNIT_TEST;
819 
820  for (const auto & elem :
821  this->_mesh->active_local_element_ptr_range())
822  for (const auto nd : elem->node_index_range())
823  {
824  if ((elem->type() == EDGE3 || elem->type() == EDGE4) && nd >= 2)
825  CPPUNIT_ASSERT(elem->is_internal(nd));
826  else if (elem->type() == HEX27 && nd == 26)
827  CPPUNIT_ASSERT(elem->is_internal(nd));
828  else if (elem->type() == PRISM21 && nd == 20)
829  CPPUNIT_ASSERT(elem->is_internal(nd));
830  else if ((elem->type() == QUAD9 || elem->type() == QUADSHELL9) && nd == 8)
831  CPPUNIT_ASSERT(elem->is_internal(nd));
832  else if (elem->type() == TRI7 && nd == 6)
833  CPPUNIT_ASSERT(elem->is_internal(nd));
834  else if (elem->type() == INFHEX18 && nd == 17)
835  CPPUNIT_ASSERT(elem->is_internal(nd));
836  else if (elem->type() == INFQUAD6 && nd == 5)
837  CPPUNIT_ASSERT(elem->is_internal(nd));
838  // Accomodate for mid-element node of C0 polyhedron
839  else if (elem->type() == C0POLYHEDRON && nd == elem->n_vertices())
840  CPPUNIT_ASSERT(elem->is_internal(nd));
841  else
842  CPPUNIT_ASSERT(!elem->is_internal(nd));
843  }
844  }
PerElemTest::_mesh
std::unique_ptr< Mesh > _mesh
Definition: elem_test.h:26

◆ 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:26

◆ 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  // We can at least check for sharing some node with
683  // a side, on both finite and infinite elements, for
684  // those children that aren't just the middle elements
685  // of triangular sides.
686  bool shares_a_side_node = false;
687  bool shares_a_vertex = false;
688  for (const Node & node : child_side->node_ref_range())
689  if (parent_side->local_node(node.id()) != invalid_uint)
690  shares_a_side_node = true;
691  for (const Node & node : elem->node_ref_range())
692  if (parent->local_node(node.id()) < parent->n_vertices())
693  shares_a_vertex = true;
694 
695  CPPUNIT_ASSERT(shares_a_side_node || !shares_a_vertex);
696 
697  if (elem->neighbor_ptr(s) && !elem->neighbor_ptr(s)->is_remote())
698  CPPUNIT_ASSERT_EQUAL(parent->child_neighbor(elem->neighbor_ptr(s)), elem);
699  }
700  }
701 
702  for (auto e : make_range(parent->n_edges()))
703  {
704  if (parent->is_child_on_edge(c,e))
705  {
706  auto parent_edge = parent->build_edge_ptr(e);
707 
708  // Implicitly assuming here that e is the child edge
709  // too - we support that now and hopefully won't have
710  // to change it later
711  auto child_edge = elem->build_edge_ptr(e);
712 
713  // 1D Inf FE inverse_map not yet implemented?
714  if (!parent_edge->infinite())
715  for (const Node & node : child_edge->node_ref_range())
716  CPPUNIT_ASSERT(parent_edge->contains_point(node));
717 
718  // We can at least check for sharing some node with
719  // an edge, on both finite and infinite elements
720  bool shares_an_edge_node = false;
721  for (const Node & node : child_edge->node_ref_range())
722  if (parent_edge->local_node(node.id()) != invalid_uint)
723  shares_an_edge_node = true;
724 
725  CPPUNIT_ASSERT(shares_an_edge_node);
726  }
727  }
728 
729  if (parent->has_affine_map())
730  CPPUNIT_ASSERT(elem->has_affine_map());
731  }
732 
733  // It's possible for a parent element on a distributed mesh to not
734  // have all its children available on any one processor
735  TestCommWorld->set_union(parent_child_was_touched);
736  TestCommWorld->set_union(parent_node_was_touched);
737 
738  for (const Elem * elem : refining_mesh->local_element_ptr_range())
739  {
740  if (elem->active())
741  continue;
742 
743  // With only one layer of refinement the family tree methods
744  // should have the full number of elements, even if some are
745  // remote.
746  if (n == 1)
747  {
748  std::vector<const Elem *> family;
749  elem->family_tree(family);
750  CPPUNIT_ASSERT_EQUAL(family.size(),
751  std::size_t(elem->n_children() + 1));
752 
753  family.clear();
754  elem->total_family_tree(family);
755  CPPUNIT_ASSERT_EQUAL(family.size(),
756  std::size_t(elem->n_children() + 1));
757 
758  family.clear();
759  elem->active_family_tree(family);
760  CPPUNIT_ASSERT_EQUAL(family.size(),
761  std::size_t(elem->n_children()));
762 
763  for (auto s : make_range(elem->n_sides()))
764  {
765  family.clear();
766  elem->active_family_tree_by_side(family,s);
767  if (!elem->build_side_ptr(s)->infinite())
768  CPPUNIT_ASSERT_EQUAL(double(family.size()),
769  std::pow(2.0, int(elem->dim()-1)));
770  else
771  CPPUNIT_ASSERT_EQUAL(double(family.size()),
772  std::pow(2.0, int(elem->dim()-2)));
773  for (const Elem * child : family)
774  {
775  if (child->is_remote())
776  continue;
777 
778  unsigned int c = elem->which_child_am_i(child);
779  CPPUNIT_ASSERT(elem->is_child_on_side(c, s));
780  }
781  }
782  }
783 
784  if (elem->level() + 1 == n)
785  {
786  for (auto c : make_range(elem->n_children()))
787  {
788  auto it = parent_child_was_touched.find(std::make_pair(elem->id(), c));
789  CPPUNIT_ASSERT(it != parent_child_was_touched.end());
790  }
791 
792  for (auto n : make_range(elem->n_nodes()))
793  {
794  auto it = parent_node_was_touched.find(std::make_pair(elem->id(), n));
795  CPPUNIT_ASSERT(it != parent_node_was_touched.end());
796  }
797  }
798  }
799 #endif
800  }
libMesh::Elem::parent
const Elem * parent() const
Definition: elem.h:3044
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:303
TestCommWorld
libMesh::Parallel::Communicator * TestCommWorld
Definition: driver.C:218
PerElemTest::_mesh
std::unique_ptr< Mesh > _mesh
Definition: elem_test.h:26
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:296
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:176
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 846 of file elem_test.C.

847  {
848  LOG_UNIT_TEST;
849 
850  for (const auto & elem : this->_mesh->active_local_element_ptr_range())
851  {
852  for (const auto nd : elem->node_index_range())
853  {
854  auto adjacent_edge_ids = elem->edges_adjacent_to_node(nd);
855 
856  if (elem->dim() < 2)
857  {
858  // 0D elements don't have edges.
859  // 1D elements *are* "edges", but we don't consider
860  // them to *have* edges, so assert that here.
861  CPPUNIT_ASSERT(adjacent_edge_ids.empty());
862  }
863  else
864  {
865  // For 2D and 3D elements, on each edge which is
866  // claimed to be adjacent, check that "nd" is indeed
867  // on it
868  for (const auto & edge_id : adjacent_edge_ids)
869  {
870  auto node_ids_on_edge = elem->nodes_on_edge(edge_id);
871  CPPUNIT_ASSERT(std::find(node_ids_on_edge.begin(), node_ids_on_edge.end(), nd) != node_ids_on_edge.end());
872  }
873  }
874  }
875  }
876  }
PerElemTest::_mesh
std::unique_ptr< Mesh > _mesh
Definition: elem_test.h:26

◆ 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:26
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:1348
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:176
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:26
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:1348
libMesh::DofObject::id
dof_id_type id() const
Definition: dof_object.h:819
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:271
libMesh::Elem::neighbor_ptr
const Elem * neighbor_ptr(unsigned int i) const
Definition: elem.h:2612
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:176
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:26
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:176
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:26
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:292

◆ test_refinement()

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

Definition at line 802 of file elem_test.C.

803  {
804  LOG_UNIT_TEST;
805 
806  test_n_refinements(1);
807  }
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:26

◆ 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:26

◆ 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:303
PerElemTest::_mesh
std::unique_ptr< Mesh > _mesh
Definition: elem_test.h:26

Member Data Documentation

◆ _mesh

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

Definition at line 26 of file elem_test.h.

◆ libmesh_suite_name

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

Definition at line 27 of file elem_test.h.

The documentation for this class was generated from the following file:
Generated on Wed Jun 24 2026 08:04:21 for libMesh by   doxygen 1.8.14