libMesh
Loading...
Searching...
No Matches
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 15 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.

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 an affine-skewed box. Unlike a unit cube, its
218 // non-unit quality values cannot be confused with the generic
219 // fallback value of 1.
220 _mesh->add_point(Point(0, 0, 0), 0);
221 _mesh->add_point(Point(2, 0, 0), 1);
222 _mesh->add_point(Point(2.25, 1, 0), 2);
223 _mesh->add_point(Point(0.25, 1, 0), 3);
224 _mesh->add_point(Point(0, 0, 1), 4);
225 _mesh->add_point(Point(2, 0, 1), 5);
226 _mesh->add_point(Point(2.25, 1, 1), 6);
227 _mesh->add_point(Point(0.25, 1, 1), 7);
228
229 // With some combinations of face triangulations, even a
230 // simple box has no tetrahedralization that doesn't have
231 // either interior discontinuities or a 0-volume pancake tet!
232 //
233 // The initial "natural" way to orient our sides is commented
234 // out here; permutations that fix diagonals for us are used
235 // instead.
236 // Listing both z sides first also tests nonconsecutive local
237 // side edges in the polyhedron edge lookup.
238 const std::vector<std::vector<unsigned int>> nodes_on_side =
239 { {0, 1, 2, 3}, // min z
240 {5, 6, 7, 4}, // max z
241 {0, 1, 5, 4}, // min y
242 // {1, 2, 6, 5}, // max x - bad
243 {2, 6, 5, 1}, // max x
244 {2, 3, 7, 6}, // max y
245 // {3, 0, 4, 7}, // min x - bad
246 {0, 4, 7, 3} }; // min x
247#endif
248
249 // Build all the sides.
250 std::vector<std::shared_ptr<Polygon>> sides(nodes_on_side.size());
251
252 for (auto s : index_range(nodes_on_side))
253 {
254 const auto & nodes_on_s = nodes_on_side[s];
255 sides[s] = std::make_shared<C0Polygon>(nodes_on_s.size());
256 for (auto i : index_range(nodes_on_s))
257 sides[s]->set_node(i, _mesh->node_ptr(nodes_on_s[i]));
258 }
259
260 std::unique_ptr<libMesh::Node> mid_elem_node;
261 std::unique_ptr<Elem> polyhedron = std::make_unique<C0Polyhedron>(sides, mid_elem_node);
262 _mesh->add_elem(std::move(polyhedron));
263 if (mid_elem_node)
264 _mesh->add_node(std::move(mid_elem_node));
265 _mesh->prepare_for_use();
266 }
267 else
268 {
269 const unsigned int dim = test_elem->dim();
270 const unsigned int use_x = dim > 0;
271 const unsigned int use_y = dim > 1;
272 const unsigned int use_z = dim > 2;
273
274 MeshTools::Generation::build_cube (*_mesh,
275 N*use_x, N*use_y, N*use_z,
276 minpos, maxpos,
277 minpos, use_y*maxpos,
278 minpos, use_z*maxpos,
279 elem_type);
280 }
281
282 // Use non-default subdomain ids so we can properly test their
283 // preservation.
284 //
285 // Use long subdomain names so we can test that we're not
286 // truncating as badly as we used to in ExodusII.
287 for (const auto & elem :
288 this->_mesh->element_ptr_range())
289 {
290 const subdomain_id_type sbdid = 10 + (elem->id() % 10);
291 elem->subdomain_id() = sbdid;
292 std::ostringstream sbdname;
293 sbdname <<
294 "a_very_long_subdomain_name_for_the_subdomain_with_number_" << sbdid;
295 this->_mesh->set_subdomain_name(sbdid, sbdname.str());
296 }
297
298 // Make sure our mesh's cache knows about them all for later
299 // test comparisons
300 this->_mesh->cache_elem_data();
301
302 // We may have updated only a portion of subdomain names if we're
303 // on a distributed mesh, but every processor should know about
304 // every name.
305 this->_mesh->sync_subdomain_name_map();
306 }
dim
unsigned int dim
Definition adaptivity_ex3.C:114
PerElemTest::_mesh
std::unique_ptr< Mesh > _mesh
Definition elem_test.h:26
libMesh::DofObject::id
dof_id_type id() const
Definition dof_object.h:819
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition elem.h:96
libMesh::Elem::set_node
virtual Node *& set_node(const unsigned int i)
Definition elem.h:2567
libMesh::Elem::build
static std::unique_ptr< Elem > build(const ElemType type, Elem *p=nullptr)
Definition elem.C:442
libMesh::Elem::subdomain_id
subdomain_id_type subdomain_id() const
Definition elem.h:2591
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition point.h:40
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:317
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::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition libmesh_common.h:144
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

References PerElemTest< elem_type >::_mesh, 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, libMesh::Elem::set_node(), and TestCommWorld.

◆ test_bounding_box()

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

Definition at line 17 of file elem_test.C.

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

References PerElemTest< elem_type >::_mesh, libMesh::BoundingBox::contains_point(), libMesh::BoundingBox::max(), libMesh::BoundingBox::min(), libMesh::BoundingBox::scale(), and libMesh::BoundingBox::union_with().

◆ test_center_node_on_side()

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

Definition at line 566 of file elem_test.C.

567 {
568 LOG_UNIT_TEST;
569
570 for (const auto & elem :
571 this->_mesh->active_local_element_ptr_range())
572 for (const auto s : elem->side_index_range())
573 {
574 if (elem->type() == EDGE2 || elem->type() == EDGE3 || elem->type() == EDGE4)
575 CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(s), elem->center_node_on_side(s));
576 else if (elem->type() == TRI6 || elem->type() == TRI7)
577 CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(s + 3), elem->center_node_on_side(s));
578 else if (elem->type() == QUAD8 || elem->type() == QUAD9 ||
579 elem->type() == QUADSHELL8 || elem->type() == QUADSHELL9)
580 CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(s + 4), elem->center_node_on_side(s));
581 else if (elem->type() == HEX27)
582 CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(s + 20), elem->center_node_on_side(s));
583 else if (elem->type() == PRISM18 && s >= 1 && s <= 3)
584 CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(s + 14), elem->center_node_on_side(s));
585 else if ((elem->type() == PRISM20 ||
586 elem->type() == PRISM21) && s >= 1 && s <= 3)
587 CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(s + 14), elem->center_node_on_side(s));
588 else if (elem->type() == PRISM20 ||
589 elem->type() == PRISM21)
590 CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(18 + (s == 4)), elem->center_node_on_side(s));
591 else if (elem->type() == PYRAMID14 && s == 4)
592 CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(13), elem->center_node_on_side(s));
593 else if (elem->type() == PYRAMID18)
594 {
595 if (s < 4)
596 CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(s + 14), elem->center_node_on_side(s));
597 else
598 CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(13), elem->center_node_on_side(s));
599 }
600 else
601 CPPUNIT_ASSERT_EQUAL(invalid_uint, elem->center_node_on_side(s));
602 }
603 }
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

References PerElemTest< elem_type >::_mesh, 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.

◆ test_contains_point_node()

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

Definition at line 291 of file elem_test.C.

292 {
293 LOG_UNIT_TEST;
294
295 for (const auto & elem :
296 this->_mesh->active_local_element_ptr_range())
297 {
298 if (elem->infinite())
299 continue;
300
301 for (const auto n : elem->node_index_range())
302#ifndef LIBMESH_ENABLE_EXCEPTIONS
303 // If this node has a singular Jacobian, we need exceptions in order
304 // to catch the failed inverse_map solve and return the singular
305 // master point. Therefore, if we don't have exceptions and we're
306 // at a singular node, we can't test this. As of the writing of
307 // this comment, this issue exists for only Pyramid elements at
308 // the apex.
309 if (elem->local_singular_node(elem->point(n), TOLERANCE*TOLERANCE) == invalid_uint)
310#endif
311 CPPUNIT_ASSERT(elem->contains_point(elem->point(n)));
312 }
313 }
libMesh::TOLERANCE
static constexpr Real TOLERANCE
Definition libmesh_common.h:151

References PerElemTest< elem_type >::_mesh, libMesh::invalid_uint, and libMesh::TOLERANCE.

◆ test_double_refinement()

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

Definition at line 874 of file elem_test.C.

875 {
876 LOG_UNIT_TEST;
877
878 test_n_refinements(2);
879 }
ElemTest::test_n_refinements
void test_n_refinements(unsigned int n)
Definition elem_test.C:660

References ElemTest< elem_type >::test_n_refinements().

◆ test_elem_side_builder()

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

Definition at line 638 of file elem_test.C.

639 {
640 LOG_UNIT_TEST;
641
642 ElemSideBuilder cache;
643 for (auto & elem : this->_mesh->active_local_element_ptr_range())
644 for (const auto s : elem->side_index_range())
645 {
646 const auto side = elem->build_side_ptr(s);
647
648 auto & cached_side = cache(*elem, s);
649 CPPUNIT_ASSERT_EQUAL(side->type(), cached_side.type());
650 for (const auto n : side->node_index_range())
651 CPPUNIT_ASSERT_EQUAL(side->node_ref(n), cached_side.node_ref(n));
652
653 const auto & const_cached_side = cache(const_cast<const Elem &>(*elem), s);
654 CPPUNIT_ASSERT_EQUAL(side->type(), const_cached_side.type());
655 for (const auto n : side->node_index_range())
656 CPPUNIT_ASSERT_EQUAL(side->node_ref(n), const_cached_side.node_ref(n));
657 }
658 }
libMesh::ElemSideBuilder
Helper for building element sides that minimizes the construction of new elements.
Definition elem_side_builder.h:51

References PerElemTest< elem_type >::_mesh.

◆ test_flip()

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

Definition at line 354 of file elem_test.C.

355 {
356 LOG_UNIT_TEST;
357
358 BoundaryInfo & boundary_info = this->_mesh->get_boundary_info();
359
360 for (const auto & elem :
361 this->_mesh->active_local_element_ptr_range())
362 {
363 if (elem->infinite())
364 continue;
365
366 if (elem->type() == C0POLYHEDRON)
367 continue;
368
369 const Point vertex_avg = elem->vertex_average();
370
371 const unsigned int n_sides = elem->n_sides();
372 std::vector<std::set<Point*>> side_nodes(n_sides);
373 std::vector<Elem*> neighbors(n_sides);
374 std::vector<std::vector<boundary_id_type>> bcids(n_sides);
375 for (auto s : make_range(n_sides))
376 {
377 for (auto n : elem->nodes_on_side(s))
378 side_nodes[s].insert(elem->node_ptr(n));
379 neighbors[s] = elem->neighbor_ptr(s);
380 boundary_info.boundary_ids(elem, s, bcids[s]);
381 }
382
383 elem->flip(&boundary_info);
384
385 // We should just be flipped, not twisted, so our map should
386 // still be affine.
387 // ... except for stupid singular pyramid maps
388 // ... or the polygons we're deliberately testing non-affine
389 if ((elem->dim() < 3 ||
390 elem->n_vertices() != 5) &&
391 elem_type != C0POLYGON)
392 CPPUNIT_ASSERT(elem->has_affine_map());
393 else if (elem_type == C0POLYGON)
394 CPPUNIT_ASSERT(!elem->has_affine_map());
395
396 // The neighbors and bcids should have flipped in a way
397 // consistently with the nodes (unless this is a 0D NodeElem)
398 bool something_changed = false;
399 for (auto s : make_range(n_sides))
400 {
401 std::set<Point*> new_side_nodes;
402 for (auto n : elem->nodes_on_side(s))
403 new_side_nodes.insert(elem->node_ptr(n));
404
405 std::vector<boundary_id_type> new_bcids;
406 boundary_info.boundary_ids(elem, s, new_bcids);
407
408 unsigned int old_side = libMesh::invalid_uint;
409 for (auto os : make_range(n_sides))
410 if (new_side_nodes == side_nodes[os])
411 old_side = os;
412
413 if (old_side != s)
414 something_changed = true;
415
416 CPPUNIT_ASSERT(old_side != libMesh::invalid_uint);
417
418 CPPUNIT_ASSERT(neighbors[old_side] ==
419 elem->neighbor_ptr(s));
420
421 CPPUNIT_ASSERT(bcids[old_side] == new_bcids);
422 }
423 CPPUNIT_ASSERT(!elem->dim() || something_changed);
424
425 const Point new_vertex_avg = elem->vertex_average();
426 for (const auto d : make_range(LIBMESH_DIM))
427 LIBMESH_ASSERT_FP_EQUAL(vertex_avg(d), new_vertex_avg(d),
428 TOLERANCE*TOLERANCE);
429 }
430 }
libMesh::BoundaryInfo
The BoundaryInfo class contains information relevant to boundary conditions including storing faces,...
Definition boundary_info.h:58
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:1389

References PerElemTest< elem_type >::_mesh, libMesh::BoundaryInfo::boundary_ids(), libMesh::C0POLYGON, libMesh::C0POLYHEDRON, libMesh::invalid_uint, libMesh::make_range(), and libMesh::TOLERANCE.

◆ test_is_internal()

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

Definition at line 881 of file elem_test.C.

882 {
883 LOG_UNIT_TEST;
884
885 for (const auto & elem :
886 this->_mesh->active_local_element_ptr_range())
887 for (const auto nd : elem->node_index_range())
888 {
889 if ((elem->type() == EDGE3 || elem->type() == EDGE4) && nd >= 2)
890 CPPUNIT_ASSERT(elem->is_internal(nd));
891 else if (elem->type() == HEX27 && nd == 26)
892 CPPUNIT_ASSERT(elem->is_internal(nd));
893 else if (elem->type() == PRISM21 && nd == 20)
894 CPPUNIT_ASSERT(elem->is_internal(nd));
895 else if ((elem->type() == QUAD9 || elem->type() == QUADSHELL9) && nd == 8)
896 CPPUNIT_ASSERT(elem->is_internal(nd));
897 else if (elem->type() == TRI7 && nd == 6)
898 CPPUNIT_ASSERT(elem->is_internal(nd));
899 else if (elem->type() == INFHEX18 && nd == 17)
900 CPPUNIT_ASSERT(elem->is_internal(nd));
901 else if (elem->type() == INFQUAD6 && nd == 5)
902 CPPUNIT_ASSERT(elem->is_internal(nd));
903 // Accomodate for mid-element node of C0 polyhedron
904 else if (elem->type() == C0POLYHEDRON && nd == elem->n_vertices())
905 CPPUNIT_ASSERT(elem->is_internal(nd));
906 else
907 CPPUNIT_ASSERT(!elem->is_internal(nd));
908 }
909 }

References PerElemTest< elem_type >::_mesh, libMesh::C0POLYHEDRON, libMesh::EDGE3, libMesh::EDGE4, libMesh::HEX27, libMesh::INFHEX18, libMesh::INFQUAD6, libMesh::PRISM21, libMesh::QUAD9, libMesh::QUADSHELL9, and libMesh::TRI7.

◆ test_maps()

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

Definition at line 231 of file elem_test.C.

232 {
233 LOG_UNIT_TEST;
234
235 for (const auto & elem :
236 this->_mesh->active_local_element_ptr_range())
237 {
238 for (const auto edge : elem->edge_index_range())
239 for (const auto side_on_edge : elem->sides_on_edge(edge))
240 for (const auto node_on_edge : elem->nodes_on_edge(edge))
241 CPPUNIT_ASSERT(elem->is_node_on_side(node_on_edge, side_on_edge));
242
243 for (const auto side : elem->side_index_range())
244 for (const auto node_on_side : elem->nodes_on_side(side))
245 CPPUNIT_ASSERT(elem->is_node_on_side(node_on_side, side));
246
247 for (const auto edge : elem->edge_index_range())
248 for (const auto node_on_edge : elem->nodes_on_edge(edge))
249 CPPUNIT_ASSERT(elem->is_node_on_edge(node_on_edge, edge));
250
251 for (const auto edge : elem->edge_index_range())
252 for (const auto side_on_edge : elem->sides_on_edge(edge))
253 CPPUNIT_ASSERT(elem->is_edge_on_side(edge, side_on_edge));
254 }
255 }

References PerElemTest< elem_type >::_mesh.

◆ test_n_refinements()

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

Definition at line 660 of file elem_test.C.

661 {
662#ifdef LIBMESH_ENABLE_AMR
663 // We don't support refinement of all element types
664 if (elem_type == EDGE4 ||
665 elem_type == PRISM20 ||
666 elem_type == PYRAMID5 ||
667 elem_type == PYRAMID13 ||
668 elem_type == PYRAMID14 ||
669 elem_type == PYRAMID18 ||
670 elem_type == C0POLYGON ||
671 elem_type == C0POLYHEDRON)
672 return;
673
674 auto refining_mesh = this->_mesh->clone();
675
676 MeshRefinement mr(*refining_mesh);
677 mr.uniformly_refine(n);
678
679 std::set<std::pair<dof_id_type, unsigned int>> parent_node_was_touched;
680 std::set<std::pair<dof_id_type, unsigned int>> parent_child_was_touched;
681
682 for (const Elem * elem : refining_mesh->active_element_ptr_range())
683 {
684 CPPUNIT_ASSERT_EQUAL(elem->level(), n);
685 CPPUNIT_ASSERT(!elem->ancestor());
686 CPPUNIT_ASSERT(elem->active());
687 CPPUNIT_ASSERT(!elem->subactive());
688 CPPUNIT_ASSERT(!elem->has_children());
689 CPPUNIT_ASSERT(!elem->has_ancestor_children());
690 CPPUNIT_ASSERT(!elem->interior_parent());
691
692 const Elem * parent = elem->parent();
693 CPPUNIT_ASSERT(parent);
694 CPPUNIT_ASSERT(parent->ancestor());
695 CPPUNIT_ASSERT(!parent->active());
696 CPPUNIT_ASSERT(!parent->subactive());
697 CPPUNIT_ASSERT(parent->has_children());
698 CPPUNIT_ASSERT(!parent->has_ancestor_children());
699 CPPUNIT_ASSERT(!parent->interior_parent());
700 if (n == 1)
701 {
702 CPPUNIT_ASSERT_EQUAL(parent, elem->top_parent());
703 CPPUNIT_ASSERT_EQUAL(parent, parent->top_parent());
704 }
705 else
706 {
707 CPPUNIT_ASSERT(parent != elem->top_parent());
708 CPPUNIT_ASSERT(parent != parent->top_parent());
709 CPPUNIT_ASSERT_EQUAL(elem->top_parent(), parent->top_parent());
710 }
711
712 CPPUNIT_ASSERT(parent->is_ancestor_of(elem));
713 const unsigned int c = parent->which_child_am_i(elem);
714 CPPUNIT_ASSERT(c < parent->n_children());
715 CPPUNIT_ASSERT_EQUAL(elem, parent->child_ptr(c));
716 parent_child_was_touched.emplace(parent->id(), c);
717
718 CPPUNIT_ASSERT_EQUAL(parent->n_nodes_in_child(c), elem->n_nodes());
719 for (auto n : make_range(elem->n_nodes()))
720 {
721 CPPUNIT_ASSERT_EQUAL(parent->is_vertex_on_child(c, n), elem->is_vertex(n));
722
723 auto pn = parent->as_parent_node(c, n);
724 CPPUNIT_ASSERT_EQUAL(pn, parent->get_node_index(elem->node_ptr(n)));
725 if (pn == libMesh::invalid_uint)
726 continue;
727 CPPUNIT_ASSERT_EQUAL(parent->is_vertex_on_parent(c, n), parent->is_vertex(pn));
728 parent_node_was_touched.emplace(parent->id(), pn);
729 }
730
731 for (auto s : make_range(parent->n_sides()))
732 {
733 if (parent->is_child_on_side(c,s))
734 {
735 auto parent_side = parent->build_side_ptr(s);
736
737 // Implicitly assuming here that s is the child side
738 // too - we support that now and hopefully won't have
739 // to change it later
740 auto child_side = elem->build_side_ptr(s);
741
742 // 2D Inf FE inverse_map not yet implemented?
743 if (!parent_side->infinite())
744 for (const Node & node : child_side->node_ref_range())
745 CPPUNIT_ASSERT(parent_side->contains_point(node));
746
747 // We can at least check for sharing some node with
748 // a side, on both finite and infinite elements, for
749 // those children that aren't just the middle elements
750 // of triangular sides.
751 bool shares_a_side_node = false;
752 bool shares_a_vertex = false;
753 for (const Node & node : child_side->node_ref_range())
754 if (parent_side->local_node(node.id()) != invalid_uint)
755 shares_a_side_node = true;
756 for (const Node & node : elem->node_ref_range())
757 if (parent->local_node(node.id()) < parent->n_vertices())
758 shares_a_vertex = true;
759
760 CPPUNIT_ASSERT(shares_a_side_node || !shares_a_vertex);
761
762 if (elem->neighbor_ptr(s) && !elem->neighbor_ptr(s)->is_remote())
763 CPPUNIT_ASSERT_EQUAL(parent->child_neighbor(elem->neighbor_ptr(s)), elem);
764 }
765 }
766
767 for (auto e : make_range(parent->n_edges()))
768 {
769 if (parent->is_child_on_edge(c,e))
770 {
771 auto parent_edge = parent->build_edge_ptr(e);
772
773 // Implicitly assuming here that e is the child edge
774 // too - we support that now and hopefully won't have
775 // to change it later
776 auto child_edge = elem->build_edge_ptr(e);
777
778 // 1D Inf FE inverse_map not yet implemented?
779 if (!parent_edge->infinite())
780 for (const Node & node : child_edge->node_ref_range())
781 CPPUNIT_ASSERT(parent_edge->contains_point(node));
782
783 // We can at least check for sharing some node with
784 // an edge, on both finite and infinite elements
785 bool shares_an_edge_node = false;
786 for (const Node & node : child_edge->node_ref_range())
787 if (parent_edge->local_node(node.id()) != invalid_uint)
788 shares_an_edge_node = true;
789
790 CPPUNIT_ASSERT(shares_an_edge_node);
791 }
792 }
793
794 if (parent->has_affine_map())
795 CPPUNIT_ASSERT(elem->has_affine_map());
796 }
797
798 // It's possible for a parent element on a distributed mesh to not
799 // have all its children available on any one processor
800 TestCommWorld->set_union(parent_child_was_touched);
801 TestCommWorld->set_union(parent_node_was_touched);
802
803 for (const Elem * elem : refining_mesh->local_element_ptr_range())
804 {
805 if (elem->active())
806 continue;
807
808 // With only one layer of refinement the family tree methods
809 // should have the full number of elements, even if some are
810 // remote.
811 if (n == 1)
812 {
813 std::vector<const Elem *> family;
814 elem->family_tree(family);
815 CPPUNIT_ASSERT_EQUAL(family.size(),
816 std::size_t(elem->n_children() + 1));
817
818 family.clear();
819 elem->total_family_tree(family);
820 CPPUNIT_ASSERT_EQUAL(family.size(),
821 std::size_t(elem->n_children() + 1));
822
823 family.clear();
824 elem->active_family_tree(family);
825 CPPUNIT_ASSERT_EQUAL(family.size(),
826 std::size_t(elem->n_children()));
827
828 for (auto s : make_range(elem->n_sides()))
829 {
830 family.clear();
831 elem->active_family_tree_by_side(family,s);
832 if (!elem->build_side_ptr(s)->infinite())
833 CPPUNIT_ASSERT_EQUAL(double(family.size()),
834 std::pow(2.0, int(elem->dim()-1)));
835 else
836 CPPUNIT_ASSERT_EQUAL(double(family.size()),
837 std::pow(2.0, int(elem->dim()-2)));
838 for (const Elem * child : family)
839 {
840 if (child->is_remote())
841 continue;
842
843 unsigned int c = elem->which_child_am_i(child);
844 CPPUNIT_ASSERT(elem->is_child_on_side(c, s));
845 }
846 }
847 }
848
849 if (elem->level() + 1 == n)
850 {
851 for (auto c : make_range(elem->n_children()))
852 {
853 auto it = parent_child_was_touched.find(std::make_pair(elem->id(), c));
854 CPPUNIT_ASSERT(it != parent_child_was_touched.end());
855 }
856
857 for (auto n : make_range(elem->n_nodes()))
858 {
859 auto it = parent_node_was_touched.find(std::make_pair(elem->id(), n));
860 CPPUNIT_ASSERT(it != parent_node_was_touched.end());
861 }
862 }
863 }
864#endif
865 }
libMesh::Elem::has_ancestor_children
bool has_ancestor_children() const
Definition elem.h:3010
libMesh::Elem::n_nodes_in_child
virtual unsigned int n_nodes_in_child(unsigned int) const
Definition elem.h:667
libMesh::Elem::is_child_on_side
virtual bool is_child_on_side(const unsigned int c, const unsigned int s) const =0
libMesh::Elem::active
bool active() const
Definition elem.h:2958
libMesh::Elem::is_ancestor_of
bool is_ancestor_of(const Elem *descendant) const
Definition elem.h:3026
libMesh::Elem::child_neighbor
Elem * child_neighbor(Elem *elem)
Definition elem.h:2657
libMesh::Elem::is_child_on_edge
virtual bool is_child_on_edge(const unsigned int c, const unsigned int e) const
Definition elem.C:2333
libMesh::Elem::has_children
bool has_children() const
Definition elem.h:2996
libMesh::Elem::parent
const Elem * parent() const
Definition elem.h:3047
libMesh::Elem::child_ptr
const Elem * child_ptr(unsigned int i) const
Definition elem.h:3180
libMesh::Elem::ancestor
bool ancestor() const
Definition elem.C:2019
libMesh::Elem::is_vertex_on_parent
virtual bool is_vertex_on_parent(unsigned int c, unsigned int n) const
Definition elem.C:3532
libMesh::Elem::which_child_am_i
unsigned int which_child_am_i(const Elem *e) const
Definition elem.h:3209
libMesh::Elem::is_vertex_on_child
virtual bool is_vertex_on_child(unsigned int, unsigned int n) const
Definition elem.h:773
libMesh::Elem::as_parent_node
virtual unsigned int as_parent_node(unsigned int c, unsigned int n) const
Definition elem.C:2413
libMesh::Elem::build_edge_ptr
virtual std::unique_ptr< Elem > build_edge_ptr(const unsigned int i)=0
libMesh::Elem::is_vertex
virtual bool is_vertex(const unsigned int i) const =0
libMesh::Elem::get_node_index
unsigned int get_node_index(const Node *node_ptr) const
Definition elem.h:2554
libMesh::Elem::interior_parent
const Elem * interior_parent() const
Definition elem.C:1160
libMesh::Elem::subactive
bool subactive() const
Definition elem.h:2976
libMesh::Elem::has_affine_map
virtual bool has_affine_map() const
Definition elem.h:1191
libMesh::Elem::top_parent
const Elem * top_parent() const
Definition elem.h:3073
libMesh::Elem::build_side_ptr
virtual std::unique_ptr< Elem > build_side_ptr(const unsigned int i)=0
libMesh::MeshRefinement
Implements (adaptive) mesh refinement algorithms for a MeshBase.
Definition mesh_refinement.h:62
libMesh::Node
A Node is like a Point, but with more information.
Definition node.h:55
TestCommWorld
Communicator * TestCommWorld
n_nodes
const dof_id_type n_nodes
Definition tecplot_io.C:67

References PerElemTest< elem_type >::_mesh, libMesh::Elem::active(), libMesh::Elem::ancestor(), libMesh::Elem::as_parent_node(), libMesh::Elem::build_edge_ptr(), libMesh::Elem::build_side_ptr(), libMesh::C0POLYGON, libMesh::C0POLYHEDRON, libMesh::Elem::child_neighbor(), libMesh::Elem::child_ptr(), libMesh::EDGE4, libMesh::Elem::get_node_index(), libMesh::Elem::has_affine_map(), libMesh::Elem::has_ancestor_children(), libMesh::Elem::has_children(), libMesh::DofObject::id(), libMesh::Elem::interior_parent(), libMesh::invalid_uint, libMesh::Elem::is_ancestor_of(), libMesh::Elem::is_child_on_edge(), libMesh::Elem::is_child_on_side(), libMesh::Elem::is_vertex(), libMesh::Elem::is_vertex_on_child(), libMesh::Elem::is_vertex_on_parent(), libMesh::Elem::local_node(), libMesh::make_range(), libMesh::Elem::n_edges(), libMesh::Elem::n_nodes_in_child(), libMesh::Elem::n_sides(), libMesh::Elem::n_vertices(), libMesh::Elem::parent(), libMesh::PRISM20, libMesh::PYRAMID13, libMesh::PYRAMID14, libMesh::PYRAMID18, libMesh::PYRAMID5, libMesh::Elem::subactive(), TestCommWorld, libMesh::Elem::top_parent(), libMesh::MeshRefinement::uniformly_refine(), and libMesh::Elem::which_child_am_i().

Referenced by ElemTest< elem_type >::test_double_refinement(), and ElemTest< elem_type >::test_refinement().

◆ test_node_edge_map_consistency()

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

Definition at line 911 of file elem_test.C.

912 {
913 LOG_UNIT_TEST;
914
915 for (const auto & elem : this->_mesh->active_local_element_ptr_range())
916 {
917 for (const auto nd : elem->node_index_range())
918 {
919 auto adjacent_edge_ids = elem->edges_adjacent_to_node(nd);
920
921 if (elem->dim() < 2)
922 {
923 // 0D elements don't have edges.
924 // 1D elements *are* "edges", but we don't consider
925 // them to *have* edges, so assert that here.
926 CPPUNIT_ASSERT(adjacent_edge_ids.empty());
927 }
928 else
929 {
930 // For 2D and 3D elements, on each edge which is
931 // claimed to be adjacent, check that "nd" is indeed
932 // on it
933 for (const auto & edge_id : adjacent_edge_ids)
934 {
935 auto node_ids_on_edge = elem->nodes_on_edge(edge_id);
936 CPPUNIT_ASSERT(std::find(node_ids_on_edge.begin(), node_ids_on_edge.end(), nd) != node_ids_on_edge.end());
937 }
938 }
939 }
940 }
941 }

References PerElemTest< elem_type >::_mesh.

◆ test_orient()

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

Definition at line 432 of file elem_test.C.

433 {
434 LOG_UNIT_TEST;
435
436 BoundaryInfo & boundary_info = this->_mesh->get_boundary_info();
437
438 for (const auto & elem :
439 this->_mesh->active_local_element_ptr_range())
440 {
441 if (elem->infinite())
442 continue;
443
444 const Point vertex_avg = elem->vertex_average();
445
446 const unsigned int n_sides = elem->n_sides();
447 std::vector<std::set<Point*>> side_nodes(n_sides);
448 std::vector<Elem*> neighbors(n_sides);
449 std::vector<std::vector<boundary_id_type>> bcids(n_sides);
450 for (auto s : make_range(n_sides))
451 {
452 for (auto n : elem->nodes_on_side(s))
453 side_nodes[s].insert(elem->node_ptr(n));
454 neighbors[s] = elem->neighbor_ptr(s);
455 boundary_info.boundary_ids(elem, s, bcids[s]);
456 }
457
458 CPPUNIT_ASSERT(!elem->is_flipped());
459
460 if (elem->id()%2)
461 {
462 elem->flip(&boundary_info);
463 CPPUNIT_ASSERT(elem->is_flipped());
464 }
465
466 elem->orient(&boundary_info);
467 CPPUNIT_ASSERT(!elem->is_flipped());
468
469 // Our map should still be affine.
470 // ... except for stupid singular pyramid maps
471 // ... or the polygons we're deliberately testing non-affine
472 // ... or the polyhedra we deliberately don't define "affine
473 // map" for.
474 if ((elem->dim() < 3 ||
475 elem->n_vertices() != 5) &&
476 elem_type != C0POLYGON &&
477 elem_type != C0POLYHEDRON)
478 CPPUNIT_ASSERT(elem->has_affine_map());
479 else if (elem_type == C0POLYGON &&
480 elem_type != C0POLYHEDRON)
481 CPPUNIT_ASSERT(!elem->has_affine_map());
482
483 // The neighbors and bcids should have flipped back to where
484 // they were.
485 for (auto s : make_range(n_sides))
486 {
487 std::set<Point*> new_side_nodes;
488 for (auto n : elem->nodes_on_side(s))
489 new_side_nodes.insert(elem->node_ptr(n));
490
491 std::vector<boundary_id_type> new_bcids;
492 boundary_info.boundary_ids(elem, s, new_bcids);
493
494 CPPUNIT_ASSERT(side_nodes[s] ==
495 new_side_nodes);
496
497 CPPUNIT_ASSERT(neighbors[s] ==
498 elem->neighbor_ptr(s));
499
500 CPPUNIT_ASSERT(bcids[s] == new_bcids);
501 }
502
503 const Point new_vertex_avg = elem->vertex_average();
504 for (const auto d : make_range(LIBMESH_DIM))
505 LIBMESH_ASSERT_FP_EQUAL(vertex_avg(d), new_vertex_avg(d),
506 TOLERANCE*TOLERANCE);
507 }
508 }

References PerElemTest< elem_type >::_mesh, libMesh::BoundaryInfo::boundary_ids(), libMesh::C0POLYGON, libMesh::C0POLYHEDRON, libMesh::make_range(), and libMesh::TOLERANCE.

◆ test_orient_elements()

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

Definition at line 510 of file elem_test.C.

511 {
512 LOG_UNIT_TEST;
513
514 const Mesh old_mesh {*this->_mesh};
515
516 BoundaryInfo & boundary_info = this->_mesh->get_boundary_info();
517 const BoundaryInfo & old_boundary_info = old_mesh.get_boundary_info();
518 CPPUNIT_ASSERT(&boundary_info != &old_boundary_info);
519
520 for (const auto & elem :
521 this->_mesh->active_local_element_ptr_range())
522 {
523 if (elem->infinite())
524 continue;
525
526 if (elem->id()%2)
527 {
528 elem->flip(&boundary_info);
529 CPPUNIT_ASSERT(elem->is_flipped());
530 }
531 }
532
533 MeshTools::Modification::orient_elements(*this->_mesh);
534
535 // I should really create a MeshBase::operator==()...
536 for (const auto & elem :
537 this->_mesh->active_local_element_ptr_range())
538 {
539 const Elem & old_elem = old_mesh.elem_ref(elem->id());
540
541 CPPUNIT_ASSERT(!elem->is_flipped());
542
543 // Elem::operator==() uses node ids to compare
544 CPPUNIT_ASSERT(*elem == old_elem);
545
546 const unsigned int n_sides = elem->n_sides();
547 for (auto s : make_range(n_sides))
548 {
549 std::vector<boundary_id_type> bcids, old_bcids;
550 boundary_info.boundary_ids(elem, s, bcids);
551 old_boundary_info.boundary_ids(&old_elem, s, old_bcids);
552 CPPUNIT_ASSERT(bcids == old_bcids);
553
554 if (elem->neighbor_ptr(s))
555 {
556 CPPUNIT_ASSERT(old_elem.neighbor_ptr(s));
557 CPPUNIT_ASSERT_EQUAL(elem->neighbor_ptr(s)->id(),
558 old_elem.neighbor_ptr(s)->id());
559 }
560 else
561 CPPUNIT_ASSERT(!old_elem.neighbor_ptr(s));
562 }
563 }
564 }
libMesh::Elem::neighbor_ptr
const Elem * neighbor_ptr(unsigned int i) const
Definition elem.h:2615
libMesh::Mesh
The Mesh class is a thin wrapper, around the ReplicatedMesh class by default.
Definition mesh.h:51
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:273

References PerElemTest< elem_type >::_mesh, libMesh::BoundaryInfo::boundary_ids(), libMesh::DofObject::id(), libMesh::make_range(), libMesh::Elem::neighbor_ptr(), and libMesh::MeshTools::Modification::orient_elements().

◆ test_permute()

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

Definition at line 315 of file elem_test.C.

316 {
317 LOG_UNIT_TEST;
318
319 for (const auto & elem :
320 this->_mesh->active_local_element_ptr_range())
321 {
322 if (elem->infinite())
323 continue;
324
325 const Point centroid = elem->true_centroid();
326 const Point vertex_avg = elem->vertex_average();
327 Point quasicc;
328 if (elem->dim() < 3)
329 quasicc = elem->quasicircumcenter();
330
331 for (const auto p : IntRange<unsigned int>(0, elem->n_permutations()))
332 {
333 elem->permute(p);
334 CPPUNIT_ASSERT(elem->has_invertible_map());
335 const Point new_centroid = elem->true_centroid();
336 const Point new_vertex_avg = elem->vertex_average();
337 Point new_quasicc;
338 if (elem->dim() < 3)
339 new_quasicc = elem->quasicircumcenter();
340 for (const auto d : make_range(LIBMESH_DIM))
341 {
342 // Getting a little FP error from Pyramid18
343 LIBMESH_ASSERT_FP_EQUAL(centroid(d), new_centroid(d),
344 TOLERANCE*std::sqrt(TOLERANCE));
345 LIBMESH_ASSERT_FP_EQUAL(vertex_avg(d), new_vertex_avg(d),
346 TOLERANCE*TOLERANCE);
347 LIBMESH_ASSERT_FP_EQUAL(quasicc(d), new_quasicc(d),
348 TOLERANCE*TOLERANCE);
349 }
350 }
351 }
352 }
int
void ErrorVector unsigned int
Definition adjoints_ex3.C:360
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:54

References PerElemTest< elem_type >::_mesh, libMesh::make_range(), and libMesh::TOLERANCE.

◆ test_quality()

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

Definition at line 53 of file elem_test.C.

54 {
55 LOG_UNIT_TEST;
56
57 for (const auto & elem : this->_mesh->active_local_element_ptr_range())
58 {
59 // EDGE_LENGTH_RATIO is one metric that is defined on all elements
60 const Real edge_length_ratio = elem->quality(EDGE_LENGTH_RATIO);
61
62 // We use "0" to mean infinity rather than inf or NaN, and
63 // every quality other than that should be 1 or larger (worse)
64 CPPUNIT_ASSERT_LESSEQUAL(edge_length_ratio, Real(1)); // 1 <= edge_length_ratio
65
66 // We're building isotropic meshes, where even elements
67 // dissected from cubes ought to have tolerable quality.
68 //
69 // Worst I see is 2 on tets and the skewed C0Polyhedron, but
70 // let's add a little tolerance in case we decide to play with
71 // rotated meshes here later.
72 CPPUNIT_ASSERT_LESSEQUAL(Real(2+TOLERANCE), edge_length_ratio); // edge_length_ratio <= 2
73
74 // The MIN_ANGLE and MAX_ANGLE quality metrics are also defined on all elements
75 const Real min_angle = elem->quality(MIN_ANGLE);
76 const Real max_angle = elem->quality(MAX_ANGLE);
77
78 // Reference Quads/Hexes have maximum internal angles of 90
79 // degrees, but pyramids actually have an obtuse interior
80 // angle of acos(-1/3) ~ 109.47 deg formed by edge pairs
81 // {(0, 4), (2, 4)} and {(1,4), (3,4)}.
82 if (elem->type() != C0POLYGON)
83 CPPUNIT_ASSERT_LESSEQUAL((std::acos(Real(-1)/3) * 180 / libMesh::pi) + TOLERANCE, max_angle);
84
85 // We're doing some polygon testing with a squashed pentagon
86 // that has a 135 degree angle
87 CPPUNIT_ASSERT_LESSEQUAL(135 + TOLERANCE, max_angle);
88
89 const auto assert_quality =
90 [elem](const ElemQuality quality,
91 const Real expected_value,
92 const Real expected_lower_bound,
93 const Real expected_upper_bound)
94 {
95 const auto bounds = elem->qual_bounds(quality);
96 LIBMESH_ASSERT_FP_EQUAL(expected_lower_bound, bounds.first, TOLERANCE);
97 LIBMESH_ASSERT_FP_EQUAL(expected_upper_bound, bounds.second, TOLERANCE);
98 LIBMESH_ASSERT_FP_EQUAL(expected_value, elem->quality(quality), TOLERANCE);
99 };
100
101 if (elem->type() == C0POLYGON)
102 {
103 const std::vector<ElemQuality> expected = {
104 EDGE_LENGTH_RATIO,
105 JACOBIAN,
106 SCALED_JACOBIAN,
107 MAX_ANGLE,
108 MIN_ANGLE
109 };
110
111 CPPUNIT_ASSERT(expected == Quality::valid(elem->type()));
112 CPPUNIT_ASSERT_EQUAL(
113 std::string("Edge Length Ratio"),
114 Quality::name(EDGE_LENGTH_RATIO));
115
116 assert_quality(EDGE_LENGTH_RATIO, std::sqrt(Real(2)), 1., 4.);
117 assert_quality(JACOBIAN, 0.5, 0.5, 1.);
118 assert_quality(
119 SCALED_JACOBIAN, Real(1) / std::sqrt(Real(2)), 0.5, 1.);
120 assert_quality(MAX_ANGLE, 135., 108., 144.);
121 assert_quality(MIN_ANGLE, 90., 54., 108.);
122 }
123 else if (elem->type() == C0POLYHEDRON)
124 {
125 const std::vector<ElemQuality> expected = {
126 EDGE_LENGTH_RATIO,
127 MAX_ANGLE,
128 MIN_ANGLE,
129 MAX_DIHEDRAL_ANGLE,
130 MIN_DIHEDRAL_ANGLE
131 };
132
133 CPPUNIT_ASSERT(expected == Quality::valid(elem->type()));
134
135 const Real min_polyhedron_angle =
136 std::acos(Real(1) / std::sqrt(Real(17))) *
137 Real(180) / libMesh::pi;
138 assert_quality(EDGE_LENGTH_RATIO, 2., 1., 4.);
139 assert_quality(
140 MAX_ANGLE, Real(180) - min_polyhedron_angle, 60., 180.);
141 assert_quality(MIN_ANGLE, min_polyhedron_angle, 30., 180.);
142 assert_quality(MAX_DIHEDRAL_ANGLE, 90., 60., 90.);
143 assert_quality(
144 MIN_DIHEDRAL_ANGLE, min_polyhedron_angle, 30., 90.);
145 }
146
147 // Notes on minimum angle we expect to see:
148 // 1.) 1D Elements don't have interior angles, so the base
149 // class implementation currently returns 0 for those
150 // elements.
151 // 2.) Reference triangles/tetrahedra have min interior angle
152 // of 45 degrees, however, here we are checking Tets in a
153 // build_cube() mesh which have angles as small as
154 // acos(2/sqrt(6)) so we use that as our lower bound here.
155 if (elem->dim() > 1)
156 CPPUNIT_ASSERT_GREATEREQUAL((std::acos(Real(2)/std::sqrt(Real(6))) * 180 / libMesh::pi) - TOLERANCE, min_angle);
157
158 // MIN,MAX_DIHEDRAL_ANGLE are implemented for all 3D elements
159 // tested here.
160 if (elem->dim() > 2)
161 {
162 const Real min_dihedral_angle = elem->quality(MIN_DIHEDRAL_ANGLE);
163 const Real max_dihedral_angle = elem->quality(MAX_DIHEDRAL_ANGLE);
164
165 // Debugging
166 // libMesh::out << "Elem type: " << Utility::enum_to_string(elem->type())
167 // << ", min_dihedral_angle = " << min_dihedral_angle
168 // << ", max_dihedral_angle = " << max_dihedral_angle
169 // << std::endl;
170
171 // Assert that we match expected values for build_cube() meshes.
172 // * build_cube() meshes of hexes, tetrahedra, and prisms have max_dihedral_angle == 90
173 // * build_cube() meshes of pyramids have max_dihedral_angle == 60 between adjacent triangular faces
174 CPPUNIT_ASSERT_LESSEQUAL (90 + TOLERANCE, max_dihedral_angle);
175
176 // * build_cube() meshes of hexes have min_dihedral_angle == 90
177 // * build_cube() meshes of prisms, tets, and pyramids have min_dihedral_angle == 45
178 // * For the InfPrism tests, we construct a single
179 // InfPrism by hand with interior angle ~53.13 deg. so that
180 // is the minimum dihedral angle we expect in that case.
181 CPPUNIT_ASSERT_GREATEREQUAL(45 - TOLERANCE, min_dihedral_angle);
182 }
183
184 // The generic JACOBIAN metrics require exactly dim() adjacent
185 // edges at a vertex, which is not guaranteed for arbitrary
186 // C0Polyhedron topology.
187 if (elem->type() != C0POLYHEDRON)
188 {
189 // The largest non-infinite elements in this test (Hexes and
190 // Prisms) have a nodal Jacobian (volume) of 8, e.g. the 2x2x2
191 // reference Hex. The infinite elements that we construct for
192 // this testing have a max nodal area of 64 since they are
193 // created (see setUp()) with max side lengths of 4 (4*4*4=64).
194 const Real jac = elem->quality(JACOBIAN);
195 if (elem->infinite())
196 CPPUNIT_ASSERT_LESSEQUAL (64 + TOLERANCE, jac);
197 else
198 CPPUNIT_ASSERT_LESSEQUAL (8 + TOLERANCE, jac);
199
200 // The smallest 2D/3D nodal areas for regular elements in this
201 // test are found in tetrahedra, which have a minimum value of
202 // 2.0. However, we return a default value of 1.0 for 0D and
203 // 1D elements here, and we have a custom distorted-pentagon
204 // C0Polygon with a 0.5 at 3 nodes, so we handle those cases
205 // too.
206 if (elem->dim() < 2)
207 CPPUNIT_ASSERT_GREATEREQUAL(1 - TOLERANCE, jac);
208 else if (elem->type() == C0POLYGON)
209 CPPUNIT_ASSERT_GREATEREQUAL(0.5 - TOLERANCE, jac);
210 else
211 CPPUNIT_ASSERT_GREATEREQUAL(2 - TOLERANCE, jac);
212
213 // The scaled Jacobian should always be <= 1. The minimum
214 // scaled Jacobian value I observed was ~0.408248 for the
215 // tetrahedral meshes. This is consistent with the way that
216 // we generate Tet meshes with build_cube(), as we don't
217 // simply refine the reference tetrahedron in that case.
218 const Real scaled_jac = elem->quality(SCALED_JACOBIAN);
219 CPPUNIT_ASSERT_LESSEQUAL (1 + TOLERANCE, scaled_jac);
220 CPPUNIT_ASSERT_GREATEREQUAL( Real(0.4), scaled_jac);
221
222 // Debugging
223 // libMesh::out << "elem->type() = " << Utility::enum_to_string(elem->type())
224 // << ", jac = " << jac
225 // << ", scaled_jac = " << scaled_jac
226 // << std::endl;
227 }
228 }
229 }
libMesh::Quality::name
std::string name(const ElemQuality q)
This function returns a string containing some name for q.
Definition elem_quality.C:42
libMesh::Quality::valid
std::vector< ElemQuality > valid(const ElemType t)
Definition elem_quality.C:355
libMesh::pi
const Real pi
.
Definition libmesh.h:292
libMesh::ElemQuality
ElemQuality
Defines an enum for element quality metrics.
Definition enum_elem_quality.h:34

References PerElemTest< elem_type >::_mesh, 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::Quality::name(), libMesh::pi, libMesh::Real, libMesh::SCALED_JACOBIAN, libMesh::TOLERANCE, and libMesh::Quality::valid().

◆ test_refinement()

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

Definition at line 867 of file elem_test.C.

868 {
869 LOG_UNIT_TEST;
870
871 test_n_refinements(1);
872 }

References ElemTest< elem_type >::test_n_refinements().

◆ test_side_subdomain()

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

Definition at line 615 of file elem_test.C.

616 {
617 LOG_UNIT_TEST;
618
619 for (const auto & elem :
620 this->_mesh->active_local_element_ptr_range())
621 {
622 std::unique_ptr<const Elem> side;
623
624 for (const auto s : elem->side_index_range())
625 {
626 CPPUNIT_ASSERT_EQUAL(elem->build_side_ptr(s)->subdomain_id(), elem->subdomain_id());
627
628 elem->build_side_ptr(side, s);
629 CPPUNIT_ASSERT_EQUAL(side->subdomain_id(), elem->subdomain_id());
630
631 // We don't require that the "lightweight" side_ptr work
632 // the same way
633 // CPPUNIT_ASSERT_EQUAL(elem->side_ptr(s)->subdomain_id(), elem->subdomain_id());
634 }
635 }
636 }

References PerElemTest< elem_type >::_mesh.

◆ test_side_type()

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

Definition at line 605 of file elem_test.C.

606 {
607 LOG_UNIT_TEST;
608
609 for (const auto & elem :
610 this->_mesh->active_local_element_ptr_range())
611 for (const auto s : elem->side_index_range())
612 CPPUNIT_ASSERT_EQUAL(elem->build_side_ptr(s)->type(), elem->side_type(s));
613 }

References PerElemTest< elem_type >::_mesh.

◆ test_static_data()

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

Definition at line 257 of file elem_test.C.

258 {
259 LOG_UNIT_TEST;
260
261 for (const auto & elem :
262 this->_mesh->active_local_element_ptr_range())
263 {
264 CPPUNIT_ASSERT(elem->n_nodes() <= Elem::max_n_nodes);
265
266 const ElemType etype = elem->type();
267
268 CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(elem->dim()),
269 Elem::type_to_dim_map[etype]);
270 CPPUNIT_ASSERT_EQUAL(elem->default_order(),
271 Elem::type_to_default_order_map[etype]);
272
273 // If we have an element type with topology defined solely by
274 // the type, then that should match the runtime topology. If
275 // not, then we should be aware of it.
276 if (!elem->runtime_topology())
277 {
278 CPPUNIT_ASSERT_EQUAL(elem->n_nodes(), Elem::type_to_n_nodes_map[etype]);
279 CPPUNIT_ASSERT_EQUAL(elem->n_sides(), Elem::type_to_n_sides_map[etype]);
280 CPPUNIT_ASSERT_EQUAL(elem->n_edges(), Elem::type_to_n_edges_map[etype]);
281 }
282 else
283 {
284 CPPUNIT_ASSERT_EQUAL(invalid_uint, Elem::type_to_n_nodes_map[etype]);
285 CPPUNIT_ASSERT_EQUAL(invalid_uint, Elem::type_to_n_sides_map[etype]);
286 CPPUNIT_ASSERT_EQUAL(invalid_uint, Elem::type_to_n_edges_map[etype]);
287 }
288 }
289 }
libMesh::Elem::type_to_n_nodes_map
static const unsigned int type_to_n_nodes_map[INVALID_ELEM]
This array maps the integer representation of the ElemType enum to the number of nodes in the element...
Definition elem.h:643
libMesh::Elem::type_to_n_edges_map
static const unsigned int type_to_n_edges_map[INVALID_ELEM]
This array maps the integer representation of the ElemType enum to the number of edges on the element...
Definition elem.h:742
libMesh::Elem::type_to_dim_map
static const unsigned int type_to_dim_map[INVALID_ELEM]
This array maps the integer representation of the ElemType enum to the geometric dimension of the ele...
Definition elem.h:628
libMesh::Elem::type_to_n_sides_map
static const unsigned int type_to_n_sides_map[INVALID_ELEM]
This array maps the integer representation of the ElemType enum to the number of sides on the element...
Definition elem.h:678
libMesh::Elem::max_n_nodes
static const unsigned int max_n_nodes
The maximum number of nodes any element can contain.
Definition elem.h:654
libMesh::Elem::type_to_default_order_map
static const Order type_to_default_order_map[INVALID_ELEM]
This array maps the integer representation of the ElemType enum to the default approximation order of...
Definition elem.h:990
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition enum_elem_type.h:33

References PerElemTest< elem_type >::_mesh, 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.

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.

Referenced by MeshPerElemTest< elem_type >::meshes_equal_enough(), PerElemTest< elem_type >::setUp(), ElemTest< elem_type >::test_bounding_box(), ElemTest< elem_type >::test_center_node_on_side(), ElemTest< elem_type >::test_contains_point_node(), ElemTest< elem_type >::test_elem_side_builder(), ElemTest< elem_type >::test_flip(), ElemTest< elem_type >::test_is_internal(), ElemTest< elem_type >::test_maps(), ElemTest< elem_type >::test_n_refinements(), ElemTest< elem_type >::test_node_edge_map_consistency(), ElemTest< elem_type >::test_orient(), ElemTest< elem_type >::test_orient_elements(), ElemTest< elem_type >::test_permute(), ElemTest< elem_type >::test_quality(), ElemTest< elem_type >::test_side_subdomain(), ElemTest< elem_type >::test_side_type(), ElemTest< elem_type >::test_static_data(), ExodusTest< elem_type >::test_write(), and XdrIOTest< elem_type >::test_write().

◆ 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 Sun Sep 6 2026 20:28:34 for libMesh by doxygen 1.9.8