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_ref_elem ()
 
void test_static_topology ()
 
void test_higher_order_node_placement ()
 
void test_quality ()
 
void test_quality_bounds ()
 
void test_quality_permutation_invariance ()
 
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 ()
 
std::vector< const Node * > oriented_face_nodes (const Elem &elem, const unsigned int face)
 The nodes at the vertices of face, in the order Elem::oriented_face_vertices() gives.
 
void test_orientation ()
 
void test_orientation_conformity ()
 
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 17 of file elem_test.C.

Member Function Documentation

◆ oriented_face_nodes()

template<ElemType elem_type>
std::vector< const Node * > ElemTest< elem_type >::oriented_face_nodes ( const Elem &  elem,
const unsigned int  face 
)
inline

The nodes at the vertices of face, in the order Elem::oriented_face_vertices() gives.

Definition at line 828 of file elem_test.C.

829 {
830 std::vector<const Node *> sequence;
831 for (const auto v : elem.oriented_face_vertices(face))
832 sequence.push_back(elem.node_ptr(v));
833
834 return sequence;
835 }

References libMesh::Elem::node_ptr(), and libMesh::Elem::oriented_face_vertices().

Referenced by ElemTest< elem_type >::test_orientation_conformity().

◆ 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:833
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:2688
libMesh::Elem::build
static std::unique_ptr< Elem > build(const ElemType type, Elem *p=nullptr)
Definition elem.C:382
libMesh::Elem::subdomain_id
subdomain_id_type subdomain_id() const
Definition elem.h:2712
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition point.h:41
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:148
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 19 of file elem_test.C.

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

912 {
913 LOG_UNIT_TEST;
914
915 for (const auto & elem :
916 this->_mesh->active_local_element_ptr_range())
917 for (const auto s : elem->side_index_range())
918 {
919 if (elem->type() == EDGE2 || elem->type() == EDGE3 || elem->type() == EDGE4)
920 CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(s), elem->center_node_on_side(s));
921 else if (elem->type() == TRI6 || elem->type() == TRI7)
922 CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(s + 3), elem->center_node_on_side(s));
923 else if (elem->type() == QUAD8 || elem->type() == QUAD9 ||
924 elem->type() == QUADSHELL8 || elem->type() == QUADSHELL9)
925 CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(s + 4), elem->center_node_on_side(s));
926 else if (elem->type() == HEX27)
927 CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(s + 20), elem->center_node_on_side(s));
928 else if (elem->type() == PRISM18 && s >= 1 && s <= 3)
929 CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(s + 14), elem->center_node_on_side(s));
930 else if ((elem->type() == PRISM20 ||
931 elem->type() == PRISM21) && s >= 1 && s <= 3)
932 CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(s + 14), elem->center_node_on_side(s));
933 else if (elem->type() == PRISM20 ||
934 elem->type() == PRISM21)
935 CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(18 + (s == 4)), elem->center_node_on_side(s));
936 else if (elem->type() == PYRAMID14 && s == 4)
937 CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(13), elem->center_node_on_side(s));
938 else if (elem->type() == PYRAMID18)
939 {
940 if (s < 4)
941 CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(s + 14), elem->center_node_on_side(s));
942 else
943 CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(13), elem->center_node_on_side(s));
944 }
945 else
946 CPPUNIT_ASSERT_EQUAL(invalid_uint, elem->center_node_on_side(s));
947 }
948 }
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 552 of file elem_test.C.

553 {
554 LOG_UNIT_TEST;
555
556 for (const auto & elem :
557 this->_mesh->active_local_element_ptr_range())
558 {
559 if (elem->infinite())
560 continue;
561
562 for (const auto n : elem->node_index_range())
563#ifndef LIBMESH_ENABLE_EXCEPTIONS
564 // If this node has a singular Jacobian, we need exceptions in order
565 // to catch the failed inverse_map solve and return the singular
566 // master point. Therefore, if we don't have exceptions and we're
567 // at a singular node, we can't test this. As of the writing of
568 // this comment, this issue exists for only Pyramid elements at
569 // the apex.
570 if (elem->local_singular_node(elem->point(n), TOLERANCE*TOLERANCE) == invalid_uint)
571#endif
572 CPPUNIT_ASSERT(elem->contains_point(elem->point(n)));
573 }
574 }
libMesh::TOLERANCE
static constexpr Real TOLERANCE
Definition libmesh_common.h:155

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 1219 of file elem_test.C.

1220 {
1221 LOG_UNIT_TEST;
1222
1223 test_n_refinements(2);
1224 }
ElemTest::test_n_refinements
void test_n_refinements(unsigned int n)
Definition elem_test.C:1005

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 983 of file elem_test.C.

984 {
985 LOG_UNIT_TEST;
986
987 ElemSideBuilder cache;
988 for (auto & elem : this->_mesh->active_local_element_ptr_range())
989 for (const auto s : elem->side_index_range())
990 {
991 const auto side = elem->build_side_ptr(s);
992
993 auto & cached_side = cache(*elem, s);
994 CPPUNIT_ASSERT_EQUAL(side->type(), cached_side.type());
995 for (const auto n : side->node_index_range())
996 CPPUNIT_ASSERT_EQUAL(side->node_ref(n), cached_side.node_ref(n));
997
998 const auto & const_cached_side = cache(const_cast<const Elem &>(*elem), s);
999 CPPUNIT_ASSERT_EQUAL(side->type(), const_cached_side.type());
1000 for (const auto n : side->node_index_range())
1001 CPPUNIT_ASSERT_EQUAL(side->node_ref(n), const_cached_side.node_ref(n));
1002 }
1003 }
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 615 of file elem_test.C.

616 {
617 LOG_UNIT_TEST;
618
619 BoundaryInfo & boundary_info = this->_mesh->get_boundary_info();
620
621 for (const auto & elem :
622 this->_mesh->active_local_element_ptr_range())
623 {
624 if (elem->infinite())
625 continue;
626
627 if (elem->type() == C0POLYHEDRON)
628 continue;
629
630 const Point vertex_avg = elem->vertex_average();
631
632 const unsigned int n_sides = elem->n_sides();
633 std::vector<std::set<Point*>> side_nodes(n_sides);
634 std::vector<Elem*> neighbors(n_sides);
635 std::vector<std::vector<boundary_id_type>> bcids(n_sides);
636 for (auto s : make_range(n_sides))
637 {
638 for (auto n : elem->nodes_on_side(s))
639 side_nodes[s].insert(elem->node_ptr(n));
640 neighbors[s] = elem->neighbor_ptr(s);
641 boundary_info.boundary_ids(elem, s, bcids[s]);
642 }
643
644 elem->flip(&boundary_info);
645
646 // We should just be flipped, not twisted, so our map should
647 // still be affine.
648 // ... except for stupid singular pyramid maps
649 // ... or the polygons we're deliberately testing non-affine
650 if ((elem->dim() < 3 ||
651 elem->n_vertices() != 5) &&
652 elem_type != C0POLYGON)
653 CPPUNIT_ASSERT(elem->has_affine_map());
654 else if (elem_type == C0POLYGON)
655 CPPUNIT_ASSERT(!elem->has_affine_map());
656
657 // The neighbors and bcids should have flipped in a way
658 // consistently with the nodes (unless this is a 0D NodeElem)
659 bool something_changed = false;
660 for (auto s : make_range(n_sides))
661 {
662 std::set<Point*> new_side_nodes;
663 for (auto n : elem->nodes_on_side(s))
664 new_side_nodes.insert(elem->node_ptr(n));
665
666 std::vector<boundary_id_type> new_bcids;
667 boundary_info.boundary_ids(elem, s, new_bcids);
668
669 unsigned int old_side = libMesh::invalid_uint;
670 for (auto os : make_range(n_sides))
671 if (new_side_nodes == side_nodes[os])
672 old_side = os;
673
674 if (old_side != s)
675 something_changed = true;
676
677 CPPUNIT_ASSERT(old_side != libMesh::invalid_uint);
678
679 CPPUNIT_ASSERT(neighbors[old_side] ==
680 elem->neighbor_ptr(s));
681
682 CPPUNIT_ASSERT(bcids[old_side] == new_bcids);
683 }
684 CPPUNIT_ASSERT(!elem->dim() || something_changed);
685
686 const Point new_vertex_avg = elem->vertex_average();
687 for (const auto d : make_range(LIBMESH_DIM))
688 LIBMESH_ASSERT_FP_EQUAL(vertex_avg(d), new_vertex_avg(d),
689 TOLERANCE*TOLERANCE);
690 }
691 }
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_higher_order_node_placement()

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

Definition at line 146 of file elem_test.C.

147 {
148 LOG_UNIT_TEST;
149
150 // The side/edge map consistency is a static_assert in each element
151 // header now; what's left to check at runtime is where the
152 // higher-order nodes sit. Every non-vertex reference node is at
153 // the centroid of the vertices of its edge, its face, or the whole
154 // element, except for EDGE4's nodes, which trisect it.
155 for (const auto & elem : this->_mesh->active_local_element_ptr_range())
156 {
157 if (elem->infinite() || elem->runtime_topology())
158 continue;
159
160 const ElemType type = elem->type();
161
162 for (const auto i : elem->node_index_range())
163 {
164 if (elem->is_vertex(i))
165 continue;
166
167 // EDGE4's interior nodes trisect it instead
168 if (type == EDGE4)
169 {
170 LIBMESH_ASSERT_REALVEC_EQUAL(Point(i == 2 ? Real(-1)/3 : Real(1)/3),
171 elem->master_point(i),
172 TOLERANCE*TOLERANCE);
173 continue;
174 }
175
176 // Find the smallest subentity the node belongs to: an edge
177 // if it sits on one, else a face, else the element itself.
178 // Its vertices are what the node should be the centroid of.
179 std::vector<unsigned int> subentity_vertices;
180 if (elem->dim() > 1)
181 for (const auto e : elem->edge_index_range())
182 {
183 const auto nodes = elem->nodes_on_edge(e);
184 if (std::find(nodes.begin(), nodes.end(), i) != nodes.end())
185 subentity_vertices = {nodes[0], nodes[1]};
186 }
187 if (subentity_vertices.empty())
188 for (const auto s : elem->side_index_range())
189 {
190 const auto nodes = elem->nodes_on_side(s);
191 if (std::find(nodes.begin(), nodes.end(), i) != nodes.end())
192 for (auto n : nodes)
193 if (elem->is_vertex(n))
194 subentity_vertices.push_back(n);
195 }
196 if (subentity_vertices.empty())
197 for (const auto n : elem->node_index_range())
198 if (elem->is_vertex(n))
199 subentity_vertices.push_back(n);
200
201 Point centroid;
202 for (auto v : subentity_vertices)
203 centroid += elem->master_point(v);
204 centroid /= Real(subentity_vertices.size());
205
206 LIBMESH_ASSERT_REALVEC_EQUAL(centroid, elem->master_point(i),
207 TOLERANCE*TOLERANCE);
208 }
209 }
210 }
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition enum_elem_type.h:33

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

◆ test_is_internal()

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

Definition at line 1226 of file elem_test.C.

1227 {
1228 LOG_UNIT_TEST;
1229
1230 for (const auto & elem :
1231 this->_mesh->active_local_element_ptr_range())
1232 for (const auto nd : elem->node_index_range())
1233 {
1234 if ((elem->type() == EDGE3 || elem->type() == EDGE4) && nd >= 2)
1235 CPPUNIT_ASSERT(elem->is_internal(nd));
1236 else if (elem->type() == HEX27 && nd == 26)
1237 CPPUNIT_ASSERT(elem->is_internal(nd));
1238 else if (elem->type() == PRISM21 && nd == 20)
1239 CPPUNIT_ASSERT(elem->is_internal(nd));
1240 else if ((elem->type() == QUAD9 || elem->type() == QUADSHELL9) && nd == 8)
1241 CPPUNIT_ASSERT(elem->is_internal(nd));
1242 else if (elem->type() == TRI7 && nd == 6)
1243 CPPUNIT_ASSERT(elem->is_internal(nd));
1244 else if (elem->type() == INFHEX18 && nd == 17)
1245 CPPUNIT_ASSERT(elem->is_internal(nd));
1246 else if (elem->type() == INFQUAD6 && nd == 5)
1247 CPPUNIT_ASSERT(elem->is_internal(nd));
1248 // Accomodate for mid-element node of C0 polyhedron
1249 else if (elem->type() == C0POLYHEDRON && nd == elem->n_vertices())
1250 CPPUNIT_ASSERT(elem->is_internal(nd));
1251 else
1252 CPPUNIT_ASSERT(!elem->is_internal(nd));
1253 }
1254 }

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 492 of file elem_test.C.

493 {
494 LOG_UNIT_TEST;
495
496 for (const auto & elem :
497 this->_mesh->active_local_element_ptr_range())
498 {
499 for (const auto edge : elem->edge_index_range())
500 for (const auto side_on_edge : elem->sides_on_edge(edge))
501 for (const auto node_on_edge : elem->nodes_on_edge(edge))
502 CPPUNIT_ASSERT(elem->is_node_on_side(node_on_edge, side_on_edge));
503
504 for (const auto side : elem->side_index_range())
505 for (const auto node_on_side : elem->nodes_on_side(side))
506 CPPUNIT_ASSERT(elem->is_node_on_side(node_on_side, side));
507
508 for (const auto edge : elem->edge_index_range())
509 for (const auto node_on_edge : elem->nodes_on_edge(edge))
510 CPPUNIT_ASSERT(elem->is_node_on_edge(node_on_edge, edge));
511
512 for (const auto edge : elem->edge_index_range())
513 for (const auto side_on_edge : elem->sides_on_edge(edge))
514 CPPUNIT_ASSERT(elem->is_edge_on_side(edge, side_on_edge));
515 }
516 }

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 1005 of file elem_test.C.

1006 {
1007#ifdef LIBMESH_ENABLE_AMR
1008 // We don't support refinement of all element types
1009 if (elem_type == EDGE4 ||
1010 elem_type == PRISM20 ||
1011 elem_type == PYRAMID5 ||
1012 elem_type == PYRAMID13 ||
1013 elem_type == PYRAMID14 ||
1014 elem_type == PYRAMID18 ||
1015 elem_type == C0POLYGON ||
1016 elem_type == C0POLYHEDRON)
1017 return;
1018
1019 auto refining_mesh = this->_mesh->clone();
1020
1021 MeshRefinement mr(*refining_mesh);
1022 mr.uniformly_refine(n);
1023
1024 std::set<std::pair<dof_id_type, unsigned int>> parent_node_was_touched;
1025 std::set<std::pair<dof_id_type, unsigned int>> parent_child_was_touched;
1026
1027 for (const Elem * elem : refining_mesh->active_element_ptr_range())
1028 {
1029 CPPUNIT_ASSERT_EQUAL(elem->level(), n);
1030 CPPUNIT_ASSERT(!elem->ancestor());
1031 CPPUNIT_ASSERT(elem->active());
1032 CPPUNIT_ASSERT(!elem->subactive());
1033 CPPUNIT_ASSERT(!elem->has_children());
1034 CPPUNIT_ASSERT(!elem->has_ancestor_children());
1035 CPPUNIT_ASSERT(!elem->interior_parent());
1036
1037 const Elem * parent = elem->parent();
1038 CPPUNIT_ASSERT(parent);
1039 CPPUNIT_ASSERT(parent->ancestor());
1040 CPPUNIT_ASSERT(!parent->active());
1041 CPPUNIT_ASSERT(!parent->subactive());
1042 CPPUNIT_ASSERT(parent->has_children());
1043 CPPUNIT_ASSERT(!parent->has_ancestor_children());
1044 CPPUNIT_ASSERT(!parent->interior_parent());
1045 if (n == 1)
1046 {
1047 CPPUNIT_ASSERT_EQUAL(parent, elem->top_parent());
1048 CPPUNIT_ASSERT_EQUAL(parent, parent->top_parent());
1049 }
1050 else
1051 {
1052 CPPUNIT_ASSERT(parent != elem->top_parent());
1053 CPPUNIT_ASSERT(parent != parent->top_parent());
1054 CPPUNIT_ASSERT_EQUAL(elem->top_parent(), parent->top_parent());
1055 }
1056
1057 CPPUNIT_ASSERT(parent->is_ancestor_of(elem));
1058 const unsigned int c = parent->which_child_am_i(elem);
1059 CPPUNIT_ASSERT(c < parent->n_children());
1060 CPPUNIT_ASSERT_EQUAL(elem, parent->child_ptr(c));
1061 parent_child_was_touched.emplace(parent->id(), c);
1062
1063 CPPUNIT_ASSERT_EQUAL(parent->n_nodes_in_child(c), elem->n_nodes());
1064 for (auto n : make_range(elem->n_nodes()))
1065 {
1066 CPPUNIT_ASSERT_EQUAL(parent->is_vertex_on_child(c, n), elem->is_vertex(n));
1067
1068 auto pn = parent->as_parent_node(c, n);
1069 CPPUNIT_ASSERT_EQUAL(pn, parent->get_node_index(elem->node_ptr(n)));
1070 if (pn == libMesh::invalid_uint)
1071 continue;
1072 CPPUNIT_ASSERT_EQUAL(parent->is_vertex_on_parent(c, n), parent->is_vertex(pn));
1073 parent_node_was_touched.emplace(parent->id(), pn);
1074 }
1075
1076 for (auto s : make_range(parent->n_sides()))
1077 {
1078 if (parent->is_child_on_side(c,s))
1079 {
1080 auto parent_side = parent->build_side_ptr(s);
1081
1082 // Implicitly assuming here that s is the child side
1083 // too - we support that now and hopefully won't have
1084 // to change it later
1085 auto child_side = elem->build_side_ptr(s);
1086
1087 // 2D Inf FE inverse_map not yet implemented?
1088 if (!parent_side->infinite())
1089 for (const Node & node : child_side->node_ref_range())
1090 CPPUNIT_ASSERT(parent_side->contains_point(node));
1091
1092 // We can at least check for sharing some node with
1093 // a side, on both finite and infinite elements, for
1094 // those children that aren't just the middle elements
1095 // of triangular sides.
1096 bool shares_a_side_node = false;
1097 bool shares_a_vertex = false;
1098 for (const Node & node : child_side->node_ref_range())
1099 if (parent_side->local_node(node.id()) != invalid_uint)
1100 shares_a_side_node = true;
1101 for (const Node & node : elem->node_ref_range())
1102 if (parent->local_node(node.id()) < parent->n_vertices())
1103 shares_a_vertex = true;
1104
1105 CPPUNIT_ASSERT(shares_a_side_node || !shares_a_vertex);
1106
1107 if (elem->neighbor_ptr(s) && !elem->neighbor_ptr(s)->is_remote())
1108 CPPUNIT_ASSERT_EQUAL(parent->child_neighbor(elem->neighbor_ptr(s)), elem);
1109 }
1110 }
1111
1112 for (auto e : make_range(parent->n_edges()))
1113 {
1114 if (parent->is_child_on_edge(c,e))
1115 {
1116 auto parent_edge = parent->build_edge_ptr(e);
1117
1118 // Implicitly assuming here that e is the child edge
1119 // too - we support that now and hopefully won't have
1120 // to change it later
1121 auto child_edge = elem->build_edge_ptr(e);
1122
1123 // 1D Inf FE inverse_map not yet implemented?
1124 if (!parent_edge->infinite())
1125 for (const Node & node : child_edge->node_ref_range())
1126 CPPUNIT_ASSERT(parent_edge->contains_point(node));
1127
1128 // We can at least check for sharing some node with
1129 // an edge, on both finite and infinite elements
1130 bool shares_an_edge_node = false;
1131 for (const Node & node : child_edge->node_ref_range())
1132 if (parent_edge->local_node(node.id()) != invalid_uint)
1133 shares_an_edge_node = true;
1134
1135 CPPUNIT_ASSERT(shares_an_edge_node);
1136 }
1137 }
1138
1139 if (parent->has_affine_map())
1140 CPPUNIT_ASSERT(elem->has_affine_map());
1141 }
1142
1143 // It's possible for a parent element on a distributed mesh to not
1144 // have all its children available on any one processor
1145 TestCommWorld->set_union(parent_child_was_touched);
1146 TestCommWorld->set_union(parent_node_was_touched);
1147
1148 for (const Elem * elem : refining_mesh->local_element_ptr_range())
1149 {
1150 if (elem->active())
1151 continue;
1152
1153 // With only one layer of refinement the family tree methods
1154 // should have the full number of elements, even if some are
1155 // remote.
1156 if (n == 1)
1157 {
1158 std::vector<const Elem *> family;
1159 elem->family_tree(family);
1160 CPPUNIT_ASSERT_EQUAL(family.size(),
1161 std::size_t(elem->n_children() + 1));
1162
1163 family.clear();
1164 elem->total_family_tree(family);
1165 CPPUNIT_ASSERT_EQUAL(family.size(),
1166 std::size_t(elem->n_children() + 1));
1167
1168 family.clear();
1169 elem->active_family_tree(family);
1170 CPPUNIT_ASSERT_EQUAL(family.size(),
1171 std::size_t(elem->n_children()));
1172
1173 for (auto s : make_range(elem->n_sides()))
1174 {
1175 family.clear();
1176 elem->active_family_tree_by_side(family,s);
1177 if (!elem->build_side_ptr(s)->infinite())
1178 CPPUNIT_ASSERT_EQUAL(double(family.size()),
1179 std::pow(2.0, int(elem->dim()-1)));
1180 else
1181 CPPUNIT_ASSERT_EQUAL(double(family.size()),
1182 std::pow(2.0, int(elem->dim()-2)));
1183 for (const Elem * child : family)
1184 {
1185 if (child->is_remote())
1186 continue;
1187
1188 unsigned int c = elem->which_child_am_i(child);
1189 CPPUNIT_ASSERT(elem->is_child_on_side(c, s));
1190 }
1191 }
1192 }
1193
1194 if (elem->level() + 1 == n)
1195 {
1196 for (auto c : make_range(elem->n_children()))
1197 {
1198 auto it = parent_child_was_touched.find(std::make_pair(elem->id(), c));
1199 CPPUNIT_ASSERT(it != parent_child_was_touched.end());
1200 }
1201
1202 for (auto n : make_range(elem->n_nodes()))
1203 {
1204 auto it = parent_node_was_touched.find(std::make_pair(elem->id(), n));
1205 CPPUNIT_ASSERT(it != parent_node_was_touched.end());
1206 }
1207 }
1208 }
1209#endif
1210 }
libMesh::Elem::has_ancestor_children
bool has_ancestor_children() const
Definition elem.h:3131
libMesh::Elem::n_nodes_in_child
virtual unsigned int n_nodes_in_child(unsigned int) const
Definition elem.h:748
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:3079
libMesh::Elem::is_ancestor_of
bool is_ancestor_of(const Elem *descendant) const
Definition elem.h:3147
libMesh::Elem::child_neighbor
Elem * child_neighbor(Elem *elem)
Definition elem.h:2778
libMesh::Elem::is_child_on_edge
virtual bool is_child_on_edge(const unsigned int c, const unsigned int e) const
Definition elem.C:2435
libMesh::Elem::has_children
bool has_children() const
Definition elem.h:3117
libMesh::Elem::parent
const Elem * parent() const
Definition elem.h:3168
libMesh::Elem::child_ptr
const Elem * child_ptr(unsigned int i) const
Definition elem.h:3301
libMesh::Elem::ancestor
bool ancestor() const
Definition elem.C:2121
libMesh::Elem::is_vertex_on_parent
virtual bool is_vertex_on_parent(unsigned int c, unsigned int n) const
Definition elem.C:3935
libMesh::Elem::which_child_am_i
unsigned int which_child_am_i(const Elem *e) const
Definition elem.h:3330
libMesh::Elem::is_vertex_on_child
virtual bool is_vertex_on_child(unsigned int, unsigned int n) const
Definition elem.h:872
libMesh::Elem::as_parent_node
virtual unsigned int as_parent_node(unsigned int c, unsigned int n) const
Definition elem.C:2515
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:2675
libMesh::Elem::interior_parent
const Elem * interior_parent() const
Definition elem.C:1100
libMesh::Elem::subactive
bool subactive() const
Definition elem.h:3097
libMesh::Elem::has_affine_map
virtual bool has_affine_map() const
Definition elem.h:1312
libMesh::Elem::top_parent
const Elem * top_parent() const
Definition elem.h:3194
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 1256 of file elem_test.C.

1257 {
1258 LOG_UNIT_TEST;
1259
1260 for (const auto & elem : this->_mesh->active_local_element_ptr_range())
1261 {
1262 for (const auto nd : elem->node_index_range())
1263 {
1264 auto adjacent_edge_ids = elem->edges_adjacent_to_node(nd);
1265
1266 if (elem->dim() < 2)
1267 {
1268 // 0D elements don't have edges.
1269 // 1D elements *are* "edges", but we don't consider
1270 // them to *have* edges, so assert that here.
1271 CPPUNIT_ASSERT(adjacent_edge_ids.empty());
1272 }
1273 else
1274 {
1275 // For 2D and 3D elements, on each edge which is
1276 // claimed to be adjacent, check that "nd" is indeed
1277 // on it
1278 for (const auto & edge_id : adjacent_edge_ids)
1279 {
1280 auto node_ids_on_edge = elem->nodes_on_edge(edge_id);
1281 CPPUNIT_ASSERT(std::find(node_ids_on_edge.begin(), node_ids_on_edge.end(), nd) != node_ids_on_edge.end());
1282 }
1283 }
1284 }
1285 }
1286 }

References PerElemTest< elem_type >::_mesh.

◆ test_orient()

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

Definition at line 693 of file elem_test.C.

694 {
695 LOG_UNIT_TEST;
696
697 BoundaryInfo & boundary_info = this->_mesh->get_boundary_info();
698
699 for (const auto & elem :
700 this->_mesh->active_local_element_ptr_range())
701 {
702 if (elem->infinite())
703 continue;
704
705 const Point vertex_avg = elem->vertex_average();
706
707 const unsigned int n_sides = elem->n_sides();
708 std::vector<std::set<Point*>> side_nodes(n_sides);
709 std::vector<Elem*> neighbors(n_sides);
710 std::vector<std::vector<boundary_id_type>> bcids(n_sides);
711 for (auto s : make_range(n_sides))
712 {
713 for (auto n : elem->nodes_on_side(s))
714 side_nodes[s].insert(elem->node_ptr(n));
715 neighbors[s] = elem->neighbor_ptr(s);
716 boundary_info.boundary_ids(elem, s, bcids[s]);
717 }
718
719 CPPUNIT_ASSERT(!elem->is_flipped());
720
721 if (elem->id()%2)
722 {
723 elem->flip(&boundary_info);
724 CPPUNIT_ASSERT(elem->is_flipped());
725 }
726
727 elem->orient(&boundary_info);
728 CPPUNIT_ASSERT(!elem->is_flipped());
729
730 // Our map should still be affine.
731 // ... except for stupid singular pyramid maps
732 // ... or the polygons we're deliberately testing non-affine
733 // ... or the polyhedra we deliberately don't define "affine
734 // map" for.
735 if ((elem->dim() < 3 ||
736 elem->n_vertices() != 5) &&
737 elem_type != C0POLYGON &&
738 elem_type != C0POLYHEDRON)
739 CPPUNIT_ASSERT(elem->has_affine_map());
740 else if (elem_type == C0POLYGON &&
741 elem_type != C0POLYHEDRON)
742 CPPUNIT_ASSERT(!elem->has_affine_map());
743
744 // The neighbors and bcids should have flipped back to where
745 // they were.
746 for (auto s : make_range(n_sides))
747 {
748 std::set<Point*> new_side_nodes;
749 for (auto n : elem->nodes_on_side(s))
750 new_side_nodes.insert(elem->node_ptr(n));
751
752 std::vector<boundary_id_type> new_bcids;
753 boundary_info.boundary_ids(elem, s, new_bcids);
754
755 CPPUNIT_ASSERT(side_nodes[s] ==
756 new_side_nodes);
757
758 CPPUNIT_ASSERT(neighbors[s] ==
759 elem->neighbor_ptr(s));
760
761 CPPUNIT_ASSERT(bcids[s] == new_bcids);
762 }
763
764 const Point new_vertex_avg = elem->vertex_average();
765 for (const auto d : make_range(LIBMESH_DIM))
766 LIBMESH_ASSERT_FP_EQUAL(vertex_avg(d), new_vertex_avg(d),
767 TOLERANCE*TOLERANCE);
768 }
769 }

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 771 of file elem_test.C.

772 {
773 LOG_UNIT_TEST;
774
775 const Mesh old_mesh {*this->_mesh};
776
777 BoundaryInfo & boundary_info = this->_mesh->get_boundary_info();
778 const BoundaryInfo & old_boundary_info = old_mesh.get_boundary_info();
779 CPPUNIT_ASSERT(&boundary_info != &old_boundary_info);
780
781 for (const auto & elem :
782 this->_mesh->active_local_element_ptr_range())
783 {
784 if (elem->infinite())
785 continue;
786
787 if (elem->id()%2)
788 {
789 elem->flip(&boundary_info);
790 CPPUNIT_ASSERT(elem->is_flipped());
791 }
792 }
793
794 MeshTools::Modification::orient_elements(*this->_mesh);
795
796 // I should really create a MeshBase::operator==()...
797 for (const auto & elem :
798 this->_mesh->active_local_element_ptr_range())
799 {
800 const Elem & old_elem = old_mesh.elem_ref(elem->id());
801
802 CPPUNIT_ASSERT(!elem->is_flipped());
803
804 // Elem::operator==() uses node ids to compare
805 CPPUNIT_ASSERT(*elem == old_elem);
806
807 const unsigned int n_sides = elem->n_sides();
808 for (auto s : make_range(n_sides))
809 {
810 std::vector<boundary_id_type> bcids, old_bcids;
811 boundary_info.boundary_ids(elem, s, bcids);
812 old_boundary_info.boundary_ids(&old_elem, s, old_bcids);
813 CPPUNIT_ASSERT(bcids == old_bcids);
814
815 if (elem->neighbor_ptr(s))
816 {
817 CPPUNIT_ASSERT(old_elem.neighbor_ptr(s));
818 CPPUNIT_ASSERT_EQUAL(elem->neighbor_ptr(s)->id(),
819 old_elem.neighbor_ptr(s)->id());
820 }
821 else
822 CPPUNIT_ASSERT(!old_elem.neighbor_ptr(s));
823 }
824 }
825 }
libMesh::Elem::neighbor_ptr
const Elem * neighbor_ptr(unsigned int i) const
Definition elem.h:2736
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_orientation()

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

Definition at line 837 of file elem_test.C.

838 {
839 LOG_UNIT_TEST;
840
841 for (const auto & elem : this->_mesh->active_local_element_ptr_range())
842 {
843 if (elem->infinite() || elem->type() == C0POLYHEDRON)
844 continue;
845
846 // Permuting renumbers the vertices without moving them
847 for (const auto p : make_range(elem->n_permutations()))
848 {
849 elem->permute(p);
850
851 for (const auto e : make_range(elem->n_edges()))
852 {
853 const bool positive = elem->point(elem->local_edge_node(e, 0)) >
854 elem->point(elem->local_edge_node(e, 1));
855
856 CPPUNIT_ASSERT_EQUAL(positive, elem->positive_edge_orientation(e));
857 }
858
859 for (const auto f : make_range(elem->n_faces()))
860 {
861 const unsigned int n_vertices =
862 Elem::type_to_n_sides_map[elem->side_type(f)];
863 const std::vector<unsigned int> nodes = elem->nodes_on_side(f);
864
865 // The least vertex of the face and the direction from it
866 unsigned int least = 0;
867 for (const auto v : make_range(1u, n_vertices))
868 if (elem->point(nodes[v]) < elem->point(nodes[least]))
869 least = v;
870
871 const bool positive =
872 elem->point(nodes[(least + n_vertices - 1) % n_vertices]) <
873 elem->point(nodes[(least + 1) % n_vertices]);
874
875 const unsigned int orientation = elem->face_orientation(f);
876
877 CPPUNIT_ASSERT_LESS(2 * n_vertices, orientation);
878 CPPUNIT_ASSERT_EQUAL(2 * least + static_cast<unsigned int>(positive),
879 orientation);
880 CPPUNIT_ASSERT_EQUAL(positive, elem->positive_face_orientation(f));
881 }
882 }
883 }
884 }
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:759

References PerElemTest< elem_type >::_mesh, libMesh::C0POLYHEDRON, libMesh::make_range(), and libMesh::Elem::type_to_n_sides_map.

◆ test_orientation_conformity()

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

Definition at line 886 of file elem_test.C.

887 {
888 LOG_UNIT_TEST;
889
890 // Elements sharing a face may have different orientation indices for it, but the indices
891 // must encode the same vertex sequence.
892 for (const auto & elem : this->_mesh->active_local_element_ptr_range())
893 {
894 if (elem->infinite() || elem->type() == C0POLYHEDRON)
895 continue;
896
897 for (const auto f : make_range(elem->n_faces()))
898 {
899 const Elem * neighbor = elem->neighbor_ptr(f);
900 if (!neighbor || neighbor == remote_elem || neighbor->infinite())
901 continue;
902
903 const unsigned int f_neighbor = neighbor->which_neighbor_am_i(elem);
904
905 CPPUNIT_ASSERT(oriented_face_nodes(*elem, f) ==
906 oriented_face_nodes(*neighbor, f_neighbor));
907 }
908 }
909 }
ElemTest::oriented_face_nodes
std::vector< const Node * > oriented_face_nodes(const Elem &elem, const unsigned int face)
The nodes at the vertices of face, in the order Elem::oriented_face_vertices() gives.
Definition elem_test.C:828
libMesh::Elem::which_neighbor_am_i
unsigned int which_neighbor_am_i(const Elem *e) const
This function tells you which neighbor e is.
Definition elem.h:3057
libMesh::Elem::infinite
virtual bool infinite() const =0
libMesh::remote_elem
const RemoteElem * remote_elem
Definition remote_elem.C:57

References PerElemTest< elem_type >::_mesh, libMesh::C0POLYHEDRON, libMesh::Elem::infinite(), libMesh::make_range(), libMesh::Elem::neighbor_ptr(), ElemTest< elem_type >::oriented_face_nodes(), libMesh::remote_elem, and libMesh::Elem::which_neighbor_am_i().

◆ test_permute()

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

Definition at line 576 of file elem_test.C.

577 {
578 LOG_UNIT_TEST;
579
580 for (const auto & elem :
581 this->_mesh->active_local_element_ptr_range())
582 {
583 if (elem->infinite())
584 continue;
585
586 const Point centroid = elem->true_centroid();
587 const Point vertex_avg = elem->vertex_average();
588 Point quasicc;
589 if (elem->dim() < 3)
590 quasicc = elem->quasicircumcenter();
591
592 for (const auto p : IntRange<unsigned int>(0, elem->n_permutations()))
593 {
594 elem->permute(p);
595 CPPUNIT_ASSERT(elem->has_invertible_map());
596 const Point new_centroid = elem->true_centroid();
597 const Point new_vertex_avg = elem->vertex_average();
598 Point new_quasicc;
599 if (elem->dim() < 3)
600 new_quasicc = elem->quasicircumcenter();
601 for (const auto d : make_range(LIBMESH_DIM))
602 {
603 // Getting a little FP error from Pyramid18
604 LIBMESH_ASSERT_FP_EQUAL(centroid(d), new_centroid(d),
605 TOLERANCE*std::sqrt(TOLERANCE));
606 LIBMESH_ASSERT_FP_EQUAL(vertex_avg(d), new_vertex_avg(d),
607 TOLERANCE*TOLERANCE);
608 LIBMESH_ASSERT_FP_EQUAL(quasicc(d), new_quasicc(d),
609 TOLERANCE*TOLERANCE);
610 }
611 }
612 }
613 }
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 212 of file elem_test.C.

213 {
214 LOG_UNIT_TEST;
215
216 for (const auto & elem : this->_mesh->active_local_element_ptr_range())
217 {
218 // EDGE_LENGTH_RATIO is one metric that is defined on all elements
219 const Real edge_length_ratio = elem->quality(EDGE_LENGTH_RATIO);
220
221 // We use "0" to mean infinity rather than inf or NaN, and
222 // every quality other than that should be 1 or larger (worse)
223 CPPUNIT_ASSERT_LESSEQUAL(edge_length_ratio, Real(1)); // 1 <= edge_length_ratio
224
225 // We're building isotropic meshes, where even elements
226 // dissected from cubes ought to have tolerable quality.
227 //
228 // Worst I see is 2 on tets and the skewed C0Polyhedron, but
229 // let's add a little tolerance in case we decide to play with
230 // rotated meshes here later.
231 CPPUNIT_ASSERT_LESSEQUAL(Real(2+TOLERANCE), edge_length_ratio); // edge_length_ratio <= 2
232
233 // The MIN_ANGLE and MAX_ANGLE quality metrics are also defined on all elements
234 const Real min_angle = elem->quality(MIN_ANGLE);
235 const Real max_angle = elem->quality(MAX_ANGLE);
236
237 // Reference Quads/Hexes have maximum internal angles of 90
238 // degrees, but pyramids actually have an obtuse interior
239 // angle of acos(-1/3) ~ 109.47 deg formed by edge pairs
240 // {(0, 4), (2, 4)} and {(1,4), (3,4)}.
241 if (elem->type() != C0POLYGON)
242 CPPUNIT_ASSERT_LESSEQUAL((std::acos(Real(-1)/3) * 180 / libMesh::pi) + TOLERANCE, max_angle);
243
244 // We're doing some polygon testing with a squashed pentagon
245 // that has a 135 degree angle
246 CPPUNIT_ASSERT_LESSEQUAL(135 + TOLERANCE, max_angle);
247
248 const auto assert_quality =
249 [elem](const ElemQuality quality,
250 const Real expected_value,
251 const Real expected_lower_bound,
252 const Real expected_upper_bound)
253 {
254 const auto bounds = elem->qual_bounds(quality);
255 LIBMESH_ASSERT_FP_EQUAL(expected_lower_bound, bounds.first, TOLERANCE);
256 LIBMESH_ASSERT_FP_EQUAL(expected_upper_bound, bounds.second, TOLERANCE);
257 LIBMESH_ASSERT_FP_EQUAL(expected_value, elem->quality(quality), TOLERANCE);
258 };
259
260 if (elem->type() == C0POLYGON)
261 {
262 const std::vector<ElemQuality> expected = {
263 EDGE_LENGTH_RATIO,
264 JACOBIAN,
265 SCALED_JACOBIAN,
266 MAX_ANGLE,
267 MIN_ANGLE
268 };
269
270 CPPUNIT_ASSERT(expected == Quality::valid(elem->type()));
271 CPPUNIT_ASSERT_EQUAL(
272 std::string("Edge Length Ratio"),
273 Quality::name(EDGE_LENGTH_RATIO));
274
275 assert_quality(EDGE_LENGTH_RATIO, std::sqrt(Real(2)), 1., 4.);
276 assert_quality(JACOBIAN, 0.5, 0.5, 1.);
277 assert_quality(
278 SCALED_JACOBIAN, Real(1) / std::sqrt(Real(2)), 0.5, 1.);
279 assert_quality(MAX_ANGLE, 135., 108., 144.);
280 assert_quality(MIN_ANGLE, 90., 54., 108.);
281 }
282 else if (elem->type() == C0POLYHEDRON)
283 {
284 const std::vector<ElemQuality> expected = {
285 EDGE_LENGTH_RATIO,
286 MAX_ANGLE,
287 MIN_ANGLE,
288 MAX_DIHEDRAL_ANGLE,
289 MIN_DIHEDRAL_ANGLE
290 };
291
292 CPPUNIT_ASSERT(expected == Quality::valid(elem->type()));
293
294 const Real min_polyhedron_angle =
295 std::acos(Real(1) / std::sqrt(Real(17))) *
296 Real(180) / libMesh::pi;
297 assert_quality(EDGE_LENGTH_RATIO, 2., 1., 4.);
298 assert_quality(
299 MAX_ANGLE, Real(180) - min_polyhedron_angle, 60., 180.);
300 assert_quality(MIN_ANGLE, min_polyhedron_angle, 30., 180.);
301 assert_quality(MAX_DIHEDRAL_ANGLE, 90., 60., 90.);
302 assert_quality(
303 MIN_DIHEDRAL_ANGLE, min_polyhedron_angle, 30., 90.);
304 }
305
306 // Notes on minimum angle we expect to see:
307 // 1.) 1D Elements don't have interior angles, so the base
308 // class implementation currently returns 0 for those
309 // elements.
310 // 2.) Reference triangles/tetrahedra have min interior angle
311 // of 45 degrees, however, here we are checking Tets in a
312 // build_cube() mesh which have angles as small as
313 // acos(2/sqrt(6)) so we use that as our lower bound here.
314 if (elem->dim() > 1)
315 CPPUNIT_ASSERT_GREATEREQUAL((std::acos(Real(2)/std::sqrt(Real(6))) * 180 / libMesh::pi) - TOLERANCE, min_angle);
316
317 // MIN,MAX_DIHEDRAL_ANGLE are implemented for all 3D elements
318 // tested here.
319 if (elem->dim() > 2)
320 {
321 const Real min_dihedral_angle = elem->quality(MIN_DIHEDRAL_ANGLE);
322 const Real max_dihedral_angle = elem->quality(MAX_DIHEDRAL_ANGLE);
323
324 // Debugging
325 // libMesh::out << "Elem type: " << Utility::enum_to_string(elem->type())
326 // << ", min_dihedral_angle = " << min_dihedral_angle
327 // << ", max_dihedral_angle = " << max_dihedral_angle
328 // << std::endl;
329
330 // Assert that we match expected values for build_cube() meshes.
331 // * build_cube() meshes of hexes, tetrahedra, and prisms have max_dihedral_angle == 90
332 // * build_cube() meshes of pyramids have max_dihedral_angle == 60 between adjacent triangular faces
333 CPPUNIT_ASSERT_LESSEQUAL (90 + TOLERANCE, max_dihedral_angle);
334
335 // * build_cube() meshes of hexes have min_dihedral_angle == 90
336 // * build_cube() meshes of prisms, tets, and pyramids have min_dihedral_angle == 45
337 // * For the InfPrism tests, we construct a single
338 // InfPrism by hand with interior angle ~53.13 deg. so that
339 // is the minimum dihedral angle we expect in that case.
340 CPPUNIT_ASSERT_GREATEREQUAL(45 - TOLERANCE, min_dihedral_angle);
341 }
342
343 // The generic JACOBIAN metrics require exactly dim() adjacent
344 // edges at a vertex, which is not guaranteed for arbitrary
345 // C0Polyhedron topology.
346 if (elem->type() != C0POLYHEDRON)
347 {
348 // The largest non-infinite elements in this test (Hexes and
349 // Prisms) have a nodal Jacobian (volume) of 8, e.g. the 2x2x2
350 // reference Hex. The infinite elements that we construct for
351 // this testing have a max nodal area of 64 since they are
352 // created (see setUp()) with max side lengths of 4 (4*4*4=64).
353 const Real jac = elem->quality(JACOBIAN);
354 if (elem->infinite())
355 CPPUNIT_ASSERT_LESSEQUAL (64 + TOLERANCE, jac);
356 else
357 CPPUNIT_ASSERT_LESSEQUAL (8 + TOLERANCE, jac);
358
359 // The smallest 2D/3D nodal areas for regular elements in this
360 // test are found in tetrahedra, which have a minimum value of
361 // 2.0. However, we return a default value of 1.0 for 0D and
362 // 1D elements here, and we have a custom distorted-pentagon
363 // C0Polygon with a 0.5 at 3 nodes, so we handle those cases
364 // too.
365 if (elem->dim() < 2)
366 CPPUNIT_ASSERT_GREATEREQUAL(1 - TOLERANCE, jac);
367 else if (elem->type() == C0POLYGON)
368 CPPUNIT_ASSERT_GREATEREQUAL(0.5 - TOLERANCE, jac);
369 else
370 CPPUNIT_ASSERT_GREATEREQUAL(2 - TOLERANCE, jac);
371
372 // The scaled Jacobian should always be <= 1. The minimum
373 // scaled Jacobian value I observed was ~0.408248 for the
374 // tetrahedral meshes. This is consistent with the way that
375 // we generate Tet meshes with build_cube(), as we don't
376 // simply refine the reference tetrahedron in that case.
377 const Real scaled_jac = elem->quality(SCALED_JACOBIAN);
378 CPPUNIT_ASSERT_LESSEQUAL (1 + TOLERANCE, scaled_jac);
379 CPPUNIT_ASSERT_GREATEREQUAL( Real(0.4), scaled_jac);
380
381 // Debugging
382 // libMesh::out << "elem->type() = " << Utility::enum_to_string(elem->type())
383 // << ", jac = " << jac
384 // << ", scaled_jac = " << scaled_jac
385 // << std::endl;
386 }
387 }
388 }
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:377
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_quality_bounds()

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

Definition at line 390 of file elem_test.C.

391 {
392 LOG_UNIT_TEST;
393
394 // Every quality metric that libMesh considers valid for a type
395 // should evaluate to a value inside its suggested qual_bounds() on
396 // a reasonably-shaped element of that type.
397 for (const auto & elem :
398 this->_mesh->active_local_element_ptr_range())
399 {
400 // qual_bounds() isn't well defined for infinite elements.
401 if (elem->infinite())
402 continue;
403
404 for (const ElemQuality q : Quality::valid(elem->type()))
405 {
406 const std::pair<Real, Real> bounds = elem->qual_bounds(q);
407
408 // A metric may be listed as valid for a type without having
409 // suggested bounds defined; qual_bounds() returns (-1, -1)
410 // as a sentinel in that case, which we skip.
411 if (bounds.first == -1. && bounds.second == -1.)
412 continue;
413
414 const Real value = elem->quality(q);
415
416 // build_cube()'s particular tet dissection (splitting each
417 // cube into 6 tets) is a bit more elongated than the
418 // "suggested" SCALED_JACOBIAN range assumes: its measured
419 // floor is ~0.408 (see the SCALED_JACOBIAN check in
420 // test_quality() above), just under the suggested 0.5.
421 // That's a real, known property of this specific mesh, not
422 // a quality regression, so widen the tolerance for exactly
423 // this combination.
424 const bool is_tet = (elem->type() == TET4 ||
425 elem->type() == TET10 ||
426 elem->type() == TET14);
427 const Real tol = (q == SCALED_JACOBIAN && is_tet) ?
428 Real(0.1) : TOLERANCE;
429
430 std::ostringstream msg;
431 msg << "Quality metric " << Utility::enum_to_string(q)
432 << " on a " << Utility::enum_to_string(elem->type())
433 << " evaluated to " << value
434 << ", outside its suggested bounds ["
435 << bounds.first << ", " << bounds.second << "]";
436
437 CPPUNIT_ASSERT_MESSAGE(msg.str(),
438 value >= bounds.first - tol &&
439 value <= bounds.second + tol);
440 }
441 }
442 }
libMesh::Utility::enum_to_string
std::string enum_to_string(const T e)
value
static const bool value
Definition xdr_io.C:55

References PerElemTest< elem_type >::_mesh, libMesh::Utility::enum_to_string(), libMesh::Real, libMesh::SCALED_JACOBIAN, libMesh::TET10, libMesh::TET14, libMesh::TET4, libMesh::TOLERANCE, libMesh::Quality::valid(), and value.

◆ test_quality_permutation_invariance()

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

Definition at line 444 of file elem_test.C.

445 {
446 LOG_UNIT_TEST;
447
448 // A genuine element-shape quality metric cannot depend on how the
449 // element's nodes happen to be numbered: permuting the nodes (i.e.
450 // relabeling which physical corner sits at which local index, via
451 // one of the reference element's symmetries) leaves the physical
452 // shape -- and therefore every quality metric -- unchanged.
453 for (const auto & elem :
454 this->_mesh->active_local_element_ptr_range())
455 {
456 // Permutations of infinite elements aren't well defined.
457 if (elem->infinite())
458 continue;
459
460 const std::vector<ElemQuality> valid_metrics =
461 Quality::valid(elem->type());
462
463 if (valid_metrics.empty())
464 continue;
465
466 std::vector<Real> baseline(valid_metrics.size());
467 for (std::size_t i = 0; i != valid_metrics.size(); ++i)
468 baseline[i] = elem->quality(valid_metrics[i]);
469
470 for (const auto p : IntRange<unsigned int>(0, elem->n_permutations()))
471 {
472 elem->permute(p);
473
474 for (std::size_t i = 0; i != valid_metrics.size(); ++i)
475 {
476 const Real permuted_value = elem->quality(valid_metrics[i]);
477
478 std::ostringstream msg;
479 msg << "Quality metric " << Utility::enum_to_string(valid_metrics[i])
480 << " on a " << Utility::enum_to_string(elem->type())
481 << " changed from " << baseline[i]
482 << " to " << permuted_value
483 << " after applying permutation " << p;
484
485 CPPUNIT_ASSERT_DOUBLES_EQUAL_MESSAGE
486 (msg.str(), baseline[i], permuted_value, TOLERANCE);
487 }
488 }
489 }
490 }

References PerElemTest< elem_type >::_mesh, libMesh::Utility::enum_to_string(), libMesh::Real, libMesh::TOLERANCE, and libMesh::Quality::valid().

◆ test_ref_elem()

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

Definition at line 55 of file elem_test.C.

56 {
57 LOG_UNIT_TEST;
58
59 for (const auto & elem : this->_mesh->active_local_element_ptr_range())
60 {
61 // We don't currently support reference infinite elements
62 if (elem->infinite())
63 continue;
64
65 // We'll probably never support reference polytopes
66 if (elem->runtime_topology())
67 continue;
68
69 const Elem & ref_elem = *elem->reference_elem();
70
71 for (auto i : elem->node_index_range())
72 {
73 LIBMESH_ASSERT_REALVEC_EQUAL(elem->master_point(i),
74 ref_elem.point(i),
75 TOLERANCE*TOLERANCE);
76 }
77 }
78 }
libMesh::Elem::point
const Point & point(const unsigned int i) const
Definition elem.h:2583
libMesh::Elem::reference_elem
const Elem * reference_elem() const
Definition elem.C:508

References PerElemTest< elem_type >::_mesh, libMesh::Elem::point(), libMesh::Elem::reference_elem(), and libMesh::TOLERANCE.

◆ test_refinement()

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

Definition at line 1212 of file elem_test.C.

1213 {
1214 LOG_UNIT_TEST;
1215
1216 test_n_refinements(1);
1217 }

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 960 of file elem_test.C.

961 {
962 LOG_UNIT_TEST;
963
964 for (const auto & elem :
965 this->_mesh->active_local_element_ptr_range())
966 {
967 std::unique_ptr<const Elem> side;
968
969 for (const auto s : elem->side_index_range())
970 {
971 CPPUNIT_ASSERT_EQUAL(elem->build_side_ptr(s)->subdomain_id(), elem->subdomain_id());
972
973 elem->build_side_ptr(side, s);
974 CPPUNIT_ASSERT_EQUAL(side->subdomain_id(), elem->subdomain_id());
975
976 // We don't require that the "lightweight" side_ptr work
977 // the same way
978 // CPPUNIT_ASSERT_EQUAL(elem->side_ptr(s)->subdomain_id(), elem->subdomain_id());
979 }
980 }
981 }

References PerElemTest< elem_type >::_mesh.

◆ test_side_type()

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

Definition at line 950 of file elem_test.C.

951 {
952 LOG_UNIT_TEST;
953
954 for (const auto & elem :
955 this->_mesh->active_local_element_ptr_range())
956 for (const auto s : elem->side_index_range())
957 CPPUNIT_ASSERT_EQUAL(elem->build_side_ptr(s)->type(), elem->side_type(s));
958 }

References PerElemTest< elem_type >::_mesh.

◆ test_static_data()

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

Definition at line 518 of file elem_test.C.

519 {
520 LOG_UNIT_TEST;
521
522 for (const auto & elem :
523 this->_mesh->active_local_element_ptr_range())
524 {
525 CPPUNIT_ASSERT(elem->n_nodes() <= Elem::max_n_nodes);
526
527 const ElemType etype = elem->type();
528
529 CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(elem->dim()),
530 Elem::type_to_dim_map[etype]);
531 CPPUNIT_ASSERT_EQUAL(elem->default_order(),
532 Elem::type_to_default_order_map[etype]);
533
534 // If we have an element type with topology defined solely by
535 // the type, then that should match the runtime topology. If
536 // not, then we should be aware of it.
537 if (!elem->runtime_topology())
538 {
539 CPPUNIT_ASSERT_EQUAL(elem->n_nodes(), Elem::type_to_n_nodes_map[etype]);
540 CPPUNIT_ASSERT_EQUAL(elem->n_sides(), Elem::type_to_n_sides_map[etype]);
541 CPPUNIT_ASSERT_EQUAL(elem->n_edges(), Elem::type_to_n_edges_map[etype]);
542 }
543 else
544 {
545 CPPUNIT_ASSERT_EQUAL(invalid_uint, Elem::type_to_n_nodes_map[etype]);
546 CPPUNIT_ASSERT_EQUAL(invalid_uint, Elem::type_to_n_sides_map[etype]);
547 CPPUNIT_ASSERT_EQUAL(invalid_uint, Elem::type_to_n_edges_map[etype]);
548 }
549 }
550 }
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:841
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:649
libMesh::Elem::max_n_nodes
static const unsigned int max_n_nodes
The maximum number of nodes any element can contain.
Definition elem.h:735
libMesh::Elem::type_to_n_nodes_map
static constexpr 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:664
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:1089

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.

◆ test_static_topology()

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

Definition at line 83 of file elem_test.C.

84 {
85 LOG_UNIT_TEST;
86
87 for (const auto & elem : this->_mesh->active_local_element_ptr_range())
88 {
89 const ElemType type = elem->type();
90
91 for (const auto s : elem->side_index_range())
92 {
93 const ElemType side_type = Elem::side_type(type, s);
94 CPPUNIT_ASSERT_EQUAL(elem->side_type(s), side_type);
95
96 // A polytope's side type follows from its element type even
97 // though its side count does not, but it has no static node
98 // map to check; neither are the infinite elements' maps read
99 // by these lookups, so both stop at the side type.
100 if (elem->runtime_topology() || elem->infinite())
101 continue;
102
103 const auto nodes = elem->nodes_on_side(s);
104 CPPUNIT_ASSERT_EQUAL(std::size_t(Elem::type_to_n_nodes_map[side_type]),
105 nodes.size());
106 for (auto n : index_range(nodes))
107 {
108 CPPUNIT_ASSERT_EQUAL(elem->local_side_node(s, n),
109 Elem::local_side_node(type, s, n));
110 CPPUNIT_ASSERT_EQUAL(nodes[n],
111 Elem::local_side_node(type, s, n));
112 }
113 }
114
115 // 1D elements have no edges. A 2D element's edges are its sides,
116 // but nodes_on_edge() and local_edge_node() are not the APIs the
117 // loop above exercised, so those are checked here too.
118 if (elem->infinite() || elem->dim() < 2)
119 continue;
120
121 for (const auto e : elem->edge_index_range())
122 {
123 // There is no virtual edge_type() to check against, but the
124 // element classes state the same fact when they build an
125 // edge, so compare with that instead
126 CPPUNIT_ASSERT_EQUAL(elem->build_edge_ptr(e)->type(),
127 Elem::edge_type(type));
128
129 if (elem->runtime_topology())
130 continue;
131
132 const auto nodes = elem->nodes_on_edge(e);
133 CPPUNIT_ASSERT_EQUAL(std::size_t(Elem::type_to_n_nodes_map[Elem::edge_type(type)]),
134 nodes.size());
135 for (auto n : index_range(nodes))
136 {
137 CPPUNIT_ASSERT_EQUAL(elem->local_edge_node(e, n),
138 Elem::local_edge_node(type, e, n));
139 CPPUNIT_ASSERT_EQUAL(nodes[n],
140 Elem::local_edge_node(type, e, n));
141 }
142 }
143 }
144 }
libMesh::Elem::local_side_node
virtual unsigned int local_side_node(unsigned int side, unsigned int side_node) const =0
libMesh::Elem::edge_type
static ElemType edge_type(const ElemType t)
Definition elem.C:3380
libMesh::Elem::local_edge_node
virtual unsigned int local_edge_node(unsigned int edge, unsigned int edge_node) const =0
Similar to Elem::local_side_node(), but instead of a side id, takes an edge id and a node id on that ...
libMesh::Elem::side_type
virtual ElemType side_type(const unsigned int s) const =0

References PerElemTest< elem_type >::_mesh, libMesh::Elem::edge_type(), libMesh::index_range(), libMesh::Elem::local_edge_node(), libMesh::Elem::local_side_node(), libMesh::Elem::side_type(), and libMesh::Elem::type_to_n_nodes_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_higher_order_node_placement(), 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_orientation(), ElemTest< elem_type >::test_orientation_conformity(), ElemTest< elem_type >::test_permute(), ElemTest< elem_type >::test_quality(), ElemTest< elem_type >::test_quality_bounds(), ElemTest< elem_type >::test_quality_permutation_invariance(), ElemTest< elem_type >::test_ref_elem(), ElemTest< elem_type >::test_side_subdomain(), ElemTest< elem_type >::test_side_type(), ElemTest< elem_type >::test_static_data(), ElemTest< elem_type >::test_static_topology(), 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 Tue Nov 10 2026 04:57:41 for libMesh by doxygen 1.9.8