libMesh
Loading...
Searching...
No Matches
side_vertex_average_normal_test.C
Go to the documentation of this file.
1// libmesh includes
2#include <libmesh/cell_c0polyhedron.h>
3#include <libmesh/elem.h>
4#include <libmesh/enum_elem_type.h>
5#include <libmesh/face_c0polygon.h>
6#include <libmesh/mesh_generation.h>
7#include <libmesh/mesh_modification.h>
8#include <libmesh/mesh.h>
9#include <libmesh/reference_elem.h>
10#include <libmesh/replicated_mesh.h>
11#include <libmesh/node.h>
12#include <libmesh/enum_to_string.h>
13#include <libmesh/tensor_value.h>
14#include <libmesh/enum_elem_quality.h>
15#include <libmesh/fe_base.h>
16#include <libmesh/quadrature_gauss.h>
17
18// unit test includes
19#include "test_comm.h"
20#include "libmesh_cppunit.h"
21
22// C++ includes
23#include <iomanip>
24
25using namespace libMesh;
26
27class SideVertexAverageNormalTest : public CppUnit::TestCase
28{
29
30public:
31 LIBMESH_CPPUNIT_TEST_SUITE( SideVertexAverageNormalTest );
32 CPPUNIT_TEST( testEdge2 );
33 CPPUNIT_TEST( testEdge3 );
34 CPPUNIT_TEST( testTris );
35 CPPUNIT_TEST( testQuads );
36 CPPUNIT_TEST( testTets );
37 CPPUNIT_TEST( testPyramids );
38 CPPUNIT_TEST( testPrisms );
39 CPPUNIT_TEST( testHexes );
40 CPPUNIT_TEST( testC0Polygon );
41 CPPUNIT_TEST( testC0Polyhedron );
42#ifdef LIBMESH_ENABLE_AMR
43 CPPUNIT_TEST( testEdge3PRefine );
44#endif // LIBMESH_ENABLE_AMR
45 CPPUNIT_TEST_SUITE_END();
46
47public:
48 const RealVectorValue unit_x{1};
49 const RealVectorValue unit_y{0,1};
50 const RealVectorValue unit_z{0,0,1};
51
52 void setUp()
53 {
54 }
55
56 void tearDown()
57 {
58 }
59
60 void testEdge2()
61 {
62 LOG_UNIT_TEST;
63
64 {
65 // Reference
66 const Elem & edge2 = ReferenceElem::get(EDGE2);
67 const Point n1 = edge2.side_vertex_average_normal(0);
68 LIBMESH_ASSERT_REALVEC_EQUAL(-unit_x, n1, TOLERANCE*TOLERANCE);
69 const Point n2 = edge2.side_vertex_average_normal(1);
70 LIBMESH_ASSERT_REALVEC_EQUAL(unit_x, n2, TOLERANCE*TOLERANCE);
71 }
72 {
73 // Oriented
74 std::vector<Point> pts = {Point(1, 0, 0), Point(1, 3, 0)};
75 auto [edge2, nodes] = this->construct_elem(pts, EDGE2);
76 const Point n1 = edge2->side_vertex_average_normal(0);
77 LIBMESH_ASSERT_REALVEC_EQUAL(-unit_y, n1, TOLERANCE*TOLERANCE);
78 const Point n2 = edge2->side_vertex_average_normal(1);
79 LIBMESH_ASSERT_REALVEC_EQUAL(unit_y, n2, TOLERANCE*TOLERANCE);
80 }
81 }
82
83 void testEdge3()
84 {
85 LOG_UNIT_TEST;
86
87 {
88 // Reference
89 const Elem & edge3 = ReferenceElem::get(EDGE3);
90 const Point n1 = edge3.side_vertex_average_normal(0);
91 LIBMESH_ASSERT_REALVEC_EQUAL(-unit_x, n1, TOLERANCE*TOLERANCE);
92 const Point n2 = edge3.side_vertex_average_normal(1);
93 LIBMESH_ASSERT_REALVEC_EQUAL(unit_x, n2, TOLERANCE*TOLERANCE);
94 }
95
96 {
97 // Oriented, checked with the FE construction
98 std::vector<Point> pts = {Point(1, 0, 0), Point(1, 3, 0), Point(2.2344, 1.210293, 0)};
99 auto [edge3, nodes] = this->construct_elem(pts, EDGE3);
100 std::unique_ptr<libMesh::FEBase> fe(libMesh::FEBase::build(1, libMesh::FEType(1)));
101 libMesh::QGauss qface(0, libMesh::CONSTANT);
102 const std::vector<Point> & normals = fe->get_normals();
103 for (const auto s : make_range(edge3->n_sides()))
104 {
105 const std::unique_ptr<const Elem> face = edge3->build_side_ptr(s);
106 fe->attach_quadrature_rule(&qface);
107 fe->reinit(edge3.get(), s, TOLERANCE);
108 const Point n1 = edge3->side_vertex_average_normal(s);
109 LIBMESH_ASSERT_REALVEC_EQUAL(normals[0], n1, TOLERANCE*TOLERANCE);
110 }
111 }
112 }
113
114
115#ifdef LIBMESH_ENABLE_AMR
116 void testEdge3PRefine()
117 {
118 LOG_UNIT_TEST;
119
120 {
121 // Constructed so we can tweak the p level
122 std::vector<Point> pts = {Point(0, 1, 0), Point(0, 2, 0), Point(0, 1.5, 0)};
123 auto [edge3, nodes] = this->construct_elem(pts, EDGE3);
124 edge3->set_p_level(1);
125 const Point n1 = edge3->side_vertex_average_normal(0);
126 LIBMESH_ASSERT_REALVEC_EQUAL(-unit_y, n1, TOLERANCE*TOLERANCE);
127 const Point n2 = edge3->side_vertex_average_normal(1);
128 LIBMESH_ASSERT_REALVEC_EQUAL(unit_y, n2, TOLERANCE*TOLERANCE);
129 }
130 }
131#endif // LIBMESH_ENABLE_AMR
132
133
134 void testTris()
135 {
136 LOG_UNIT_TEST;
137
138 // Reference elements
139 const Elem * tris[] = {&ReferenceElem::get(TRI3),
140 &ReferenceElem::get(TRI6),
141 &ReferenceElem::get(TRI7)};
142
143 for (const Elem * tri : tris)
144 {
145 const Point n1 = tri->side_vertex_average_normal(0);
146 LIBMESH_ASSERT_REALVEC_EQUAL(-unit_y, n1, TOLERANCE*TOLERANCE);
147 const Point n2 = tri->side_vertex_average_normal(1);
148 LIBMESH_ASSERT_REALVEC_EQUAL((unit_x+unit_y).unit(), n2, TOLERANCE*TOLERANCE);
149 const Point n3 = tri->side_vertex_average_normal(2);
150 LIBMESH_ASSERT_REALVEC_EQUAL(-unit_x, n3, TOLERANCE*TOLERANCE);
151 }
152
153 {
154 // General shape
155 std::vector<Point> pts = {Point(1, 0, 0), Point(1, 1, 0), Point(0, 3, 1)};
156 auto [tri3, nodes] = this->construct_elem(pts, TRI3);
157 const Point n1 = tri3->side_vertex_average_normal(0);
158 LIBMESH_ASSERT_REALVEC_EQUAL((unit_x-unit_z).unit(), n1, TOLERANCE*TOLERANCE);
159 const Point n2 = tri3->side_vertex_average_normal(1);
160 LIBMESH_ASSERT_REALVEC_EQUAL((unit_x+unit_y-unit_z).unit(), n2, TOLERANCE*TOLERANCE);
161 const Point n3 = tri3->side_vertex_average_normal(2);
162 LIBMESH_ASSERT_REALVEC_EQUAL((-3*unit_x-2*unit_y+3*unit_z).unit(), n3, TOLERANCE*TOLERANCE);
163 }
164 }
165
166 void testQuads()
167 {
168 LOG_UNIT_TEST;
169
170 // Reference
171 const Elem * quads[] = {&ReferenceElem::get(QUAD4),
172 &ReferenceElem::get(QUAD8),
173 &ReferenceElem::get(QUAD9)};
174
175 for (const Elem * quad : quads)
176 {
177 const Point n1 = quad->side_vertex_average_normal(0);
178 LIBMESH_ASSERT_REALVEC_EQUAL(-unit_y, n1, TOLERANCE*TOLERANCE);
179 const Point n2 = quad->side_vertex_average_normal(1);
180 LIBMESH_ASSERT_REALVEC_EQUAL(unit_x, n2, TOLERANCE*TOLERANCE);
181 const Point n3 = quad->side_vertex_average_normal(2);
182 LIBMESH_ASSERT_REALVEC_EQUAL(unit_y, n3, TOLERANCE*TOLERANCE);
183 const Point n4 = quad->side_vertex_average_normal(3);
184 LIBMESH_ASSERT_REALVEC_EQUAL(-unit_x, n4, TOLERANCE*TOLERANCE);
185 }
186
187 {
188 // Planar, general shape
189 std::vector<Point> pts = {Point(1, 0, 0), Point(1, 3, 0), Point(-1, -1, 0), Point(0, -1, 0)};
190 auto [quad4, nodes] = this->construct_elem(pts, QUAD4);
191 const Point n1 = quad4->side_vertex_average_normal(0);
192 LIBMESH_ASSERT_REALVEC_EQUAL(unit_x, n1, TOLERANCE*TOLERANCE);
193 const Point n2 = quad4->side_vertex_average_normal(1);
194 LIBMESH_ASSERT_REALVEC_EQUAL((-2*unit_x+unit_y).unit(), n2, TOLERANCE*TOLERANCE);
195 const Point n3 = quad4->side_vertex_average_normal(2);
196 LIBMESH_ASSERT_REALVEC_EQUAL(-unit_y, n3, TOLERANCE*TOLERANCE);
197 const Point n4 = quad4->side_vertex_average_normal(3);
198 LIBMESH_ASSERT_REALVEC_EQUAL((unit_x-unit_y).unit(), n4, TOLERANCE*TOLERANCE);
199 }
200
201 {
202 // Non-planar
203 std::vector<Point> pts = {Point(0, 0, 0), Point(1, 0, 1), Point(1, 1, 0), Point(0, 1, 1)};
204 auto [quad4, nodes] = this->construct_elem(pts, QUAD4);
205 const Point n1 = quad4->side_vertex_average_normal(0);
206 LIBMESH_ASSERT_REALVEC_EQUAL(-unit_y, n1, TOLERANCE*TOLERANCE);
207 const Point n2 = quad4->side_vertex_average_normal(1);
208 LIBMESH_ASSERT_REALVEC_EQUAL(unit_x, n2, TOLERANCE*TOLERANCE);
209 const Point n3 = quad4->side_vertex_average_normal(2);
210 LIBMESH_ASSERT_REALVEC_EQUAL(unit_y, n3, TOLERANCE*TOLERANCE);
211 const Point n4 = quad4->side_vertex_average_normal(3);
212 LIBMESH_ASSERT_REALVEC_EQUAL(-unit_x, n4, TOLERANCE*TOLERANCE);
213 }
214
215 {
216 // Non-planar, general
217 std::vector<Point> pts = {Point(0, 0, 0), Point(1, -2, 3), Point(1, 1, 0), Point(0, 2, -10.5)};
218 auto [quad4, nodes] = this->construct_elem(pts, QUAD4);
219 for (const auto s : make_range(quad4->n_sides()))
220 {
221 const Point n1 = quad4->side_vertex_average_normal(s);
222 const Point normal = quad4->Elem::side_vertex_average_normal(s);
223 LIBMESH_ASSERT_REALVEC_EQUAL(normal, n1, TOLERANCE*TOLERANCE);
224 }
225 }
226 }
227
228
229 void testTets()
230 {
231 LOG_UNIT_TEST;
232
233 // Reference
234 const Elem * tets[] = {&ReferenceElem::get(TET4),
235 &ReferenceElem::get(TET10),
236 &ReferenceElem::get(TET14)};
237
238 for (const Elem * tet : tets)
239 {
240 const Point n1 = tet->side_vertex_average_normal(0);
241 LIBMESH_ASSERT_REALVEC_EQUAL(-unit_z, n1, TOLERANCE*TOLERANCE);
242 const Point n2 = tet->side_vertex_average_normal(1);
243 LIBMESH_ASSERT_REALVEC_EQUAL(-unit_y, n2, TOLERANCE*TOLERANCE);
244 const Point n3 = tet->side_vertex_average_normal(2);
245 LIBMESH_ASSERT_REALVEC_EQUAL((unit_x+unit_y+unit_z).unit(), n3, TOLERANCE*TOLERANCE);
246 const Point n4 = tet->side_vertex_average_normal(3);
247 LIBMESH_ASSERT_REALVEC_EQUAL(-unit_x, n4, TOLERANCE*TOLERANCE);
248 }
249 }
250
251 void testPyramids()
252 {
253 LOG_UNIT_TEST;
254
255 // Reference
256 const Elem * pyrs[] = {&ReferenceElem::get(PYRAMID5),
257 &ReferenceElem::get(PYRAMID13),
258 &ReferenceElem::get(PYRAMID14),
259 &ReferenceElem::get(PYRAMID18)};
260
261 for (const Elem * pyr : pyrs)
262 {
263 const Point n1 = pyr->side_vertex_average_normal(0);
264 LIBMESH_ASSERT_REALVEC_EQUAL((-unit_y+unit_z).unit(), n1, TOLERANCE*TOLERANCE);
265 const Point n2 = pyr->side_vertex_average_normal(1);
266 LIBMESH_ASSERT_REALVEC_EQUAL((unit_x+unit_z).unit(), n2, TOLERANCE*TOLERANCE);
267 const Point n3 = pyr->side_vertex_average_normal(2);
268 LIBMESH_ASSERT_REALVEC_EQUAL((unit_y+unit_z).unit(), n3, TOLERANCE*TOLERANCE);
269 const Point n4 = pyr->side_vertex_average_normal(3);
270 LIBMESH_ASSERT_REALVEC_EQUAL((-unit_x+unit_z).unit(), n4, TOLERANCE*TOLERANCE);
271 const Point n5 = pyr->side_vertex_average_normal(4);
272 LIBMESH_ASSERT_REALVEC_EQUAL(-unit_z, n5, TOLERANCE*TOLERANCE);
273 }
274 }
275
276 void testPrisms()
277 {
278 LOG_UNIT_TEST;
279
280 // Reference
281 const Elem * pris[] = {&ReferenceElem::get(PRISM6),
282 &ReferenceElem::get(PRISM15),
283 &ReferenceElem::get(PRISM18),
284 &ReferenceElem::get(PRISM20),
285 &ReferenceElem::get(PRISM21)};
286
287 for (const Elem * pri : pris)
288 {
289 const Point n1 = pri->side_vertex_average_normal(0);
290 LIBMESH_ASSERT_REALVEC_EQUAL(-unit_z, n1, TOLERANCE*TOLERANCE);
291 const Point n2 = pri->side_vertex_average_normal(1);
292 LIBMESH_ASSERT_REALVEC_EQUAL(-unit_y, n2, TOLERANCE*TOLERANCE);
293 const Point n3 = pri->side_vertex_average_normal(2);
294 LIBMESH_ASSERT_REALVEC_EQUAL((unit_x+unit_y).unit(), n3, TOLERANCE*TOLERANCE);
295 const Point n4 = pri->side_vertex_average_normal(3);
296 LIBMESH_ASSERT_REALVEC_EQUAL(-unit_x, n4, TOLERANCE*TOLERANCE);
297 const Point n5 = pri->side_vertex_average_normal(4);
298 LIBMESH_ASSERT_REALVEC_EQUAL(unit_z, n5, TOLERANCE*TOLERANCE);
299 }
300 }
301
302 void testHexes()
303 {
304 LOG_UNIT_TEST;
305
306 const Elem * hexes[] = {&ReferenceElem::get(HEX8),
307 &ReferenceElem::get(HEX20),
308 &ReferenceElem::get(HEX27)};
309
310 for (const Elem * hex : hexes)
311 {
312 const Point n1 = hex->side_vertex_average_normal(0);
313 LIBMESH_ASSERT_REALVEC_EQUAL(-unit_z, n1, TOLERANCE*TOLERANCE);
314 const Point n2 = hex->side_vertex_average_normal(1);
315 LIBMESH_ASSERT_REALVEC_EQUAL(-unit_y, n2, TOLERANCE*TOLERANCE);
316 const Point n3 = hex->side_vertex_average_normal(2);
317 LIBMESH_ASSERT_REALVEC_EQUAL(unit_x, n3, TOLERANCE*TOLERANCE);
318 const Point n4 = hex->side_vertex_average_normal(3);
319 LIBMESH_ASSERT_REALVEC_EQUAL(unit_y, n4, TOLERANCE*TOLERANCE);
320 const Point n5 = hex->side_vertex_average_normal(4);
321 LIBMESH_ASSERT_REALVEC_EQUAL(-unit_x, n5, TOLERANCE*TOLERANCE);
322 const Point n6 = hex->side_vertex_average_normal(5);
323 LIBMESH_ASSERT_REALVEC_EQUAL(unit_z, n6, TOLERANCE*TOLERANCE);
324 }
325
326 {
327 // Non-planar quad faces
328 std::vector<Point> pts = {Point(0, 0, 0), Point(1, 0, 3), Point(1, 1, 0), Point(0, 2, 1),
329 Point(-2.2, 0, 4), Point(-1, 0, 5), Point(-0.5, 1, 4), Point(-2.2, 2, 6)};
330 auto [hex8, nodes] = this->construct_elem(pts, HEX8);
331 std::unique_ptr<libMesh::FEBase> fe(libMesh::FEBase::build(3, libMesh::FEType(1)));
332 libMesh::QGauss qface(2, libMesh::CONSTANT);
333 const std::vector<Point> & normals = fe->get_normals();
334 for (const auto s : make_range(hex8->n_sides()))
335 {
336 const std::unique_ptr<const Elem> face = hex8->build_side_ptr(s);
337 fe->attach_quadrature_rule(&qface);
338 fe->reinit(hex8.get(), s, TOLERANCE);
339 const Point n1 = hex8->side_vertex_average_normal(s);
340 LIBMESH_ASSERT_REALVEC_EQUAL(normals[0], n1, TOLERANCE*TOLERANCE);
341 }
342 }
343 }
344
345 void testC0Polygon()
346 {
347 LOG_UNIT_TEST;
348 {
349 // Square
350 std::vector<Point> pts = {Point(0, 0, 0), Point(1, 0, 0), Point(1, 1, 0), Point(0, 1, 0)};
351 auto [square, nodes] = this->construct_elem(pts, C0POLYGON);
352 const Point n1 = square->side_vertex_average_normal(0);
353 LIBMESH_ASSERT_REALVEC_EQUAL(-unit_y, n1, TOLERANCE*TOLERANCE);
354 const Point n2 = square->side_vertex_average_normal(1);
355 LIBMESH_ASSERT_REALVEC_EQUAL(unit_x, n2, TOLERANCE*TOLERANCE);
356 const Point n3 = square->side_vertex_average_normal(2);
357 LIBMESH_ASSERT_REALVEC_EQUAL(unit_y, n3, TOLERANCE*TOLERANCE);
358 const Point n4 = square->side_vertex_average_normal(3);
359 LIBMESH_ASSERT_REALVEC_EQUAL(-unit_x, n4, TOLERANCE*TOLERANCE);
360 }
361 {
362 // Hexagon (not the ref one but an easy one to draw)
363 std::vector<Point> pts = {Point(0, 0, 0), Point(1, 0, 0), Point(1.5, 1, 0), Point(1, 2, 0), Point(0, 2, 0) , Point(-0.5, 1, 0)};
364 auto [hexagon, nodes] = this->construct_elem(pts, C0POLYGON);
365 const Point n1 = hexagon->side_vertex_average_normal(0);
366 LIBMESH_ASSERT_REALVEC_EQUAL(-unit_y, n1, TOLERANCE*TOLERANCE);
367 const Point n2 = hexagon->side_vertex_average_normal(1);
368 LIBMESH_ASSERT_REALVEC_EQUAL((2*unit_x-unit_y).unit(), n2, TOLERANCE*TOLERANCE);
369 const Point n3 = hexagon->side_vertex_average_normal(2);
370 LIBMESH_ASSERT_REALVEC_EQUAL((2*unit_x+unit_y).unit(), n3, TOLERANCE*TOLERANCE);
371 const Point n4 = hexagon->side_vertex_average_normal(3);
372 LIBMESH_ASSERT_REALVEC_EQUAL(unit_y, n4, TOLERANCE*TOLERANCE);
373 const Point n5 = hexagon->side_vertex_average_normal(4);
374 LIBMESH_ASSERT_REALVEC_EQUAL((-2*unit_x+unit_y).unit(), n5, TOLERANCE*TOLERANCE);
375 const Point n6 = hexagon->side_vertex_average_normal(5);
376 LIBMESH_ASSERT_REALVEC_EQUAL((-2*unit_x-unit_y).unit(), n6, TOLERANCE*TOLERANCE);
377 }
378 {
379 // Non-planar quad (see quad unit test)
380 // Note that we don't necessarily allow non-planar polys at this time
381 std::vector<Point> pts = {Point(0, 0, 0), Point(1, 0, 3), Point(1, 1, 0), Point(0, 2, 1)};
382 auto [poly4, nodes] = this->construct_elem(pts, C0POLYGON);
383 // Use a quad4 to get the normals
384 auto [quad4, nodes2] = this->construct_elem(pts, QUAD4);
385 const std::unique_ptr<const Elem> face = quad4->build_side_ptr(0);
386 std::unique_ptr<libMesh::FEBase> fe(libMesh::FEBase::build(2, libMesh::FEType(1)));
387 libMesh::QGauss qface(1, libMesh::CONSTANT);
388 fe->attach_quadrature_rule(&qface);
389 const std::vector<Point> & normals = fe->get_normals();
390 for (const auto s : make_range(quad4->n_sides()))
391 {
392 const std::unique_ptr<const Elem> face = quad4->build_side_ptr(s);
393 fe->attach_quadrature_rule(&qface);
394 fe->reinit(quad4.get(), s, TOLERANCE);
395 const Point n1 = quad4->side_vertex_average_normal(s);
396 LIBMESH_ASSERT_REALVEC_EQUAL(normals[0], n1, TOLERANCE*TOLERANCE);
397 }
398 }
399 }
400
401 void testC0Polyhedron()
402 {
403 LOG_UNIT_TEST;
404 {
405 // Cube
406 std::vector<Point> points = {Point(0, 0, 0), Point(1, 0, 0), Point(1, 1, 0), Point(0, 1, 0),
407 Point(0, 0, 1), Point(1, 0, 1), Point(1, 1, 1), Point(0, 1, 1)};
408
409 // See notes in elem_test.h
410 const std::vector<std::vector<unsigned int>> nodes_on_side =
411 { {0, 1, 2, 3}, // min z
412 {0, 1, 5, 4}, // min y
413 {2, 6, 5, 1}, // max x
414 {2, 3, 7, 6}, // max y
415 {0, 4, 7, 3}, // min x
416 {5, 6, 7, 4} }; // max z
417 // Note: we don't test the alternative triangulation for side vertex average normal
418 // as we don't expect any difference, the sides are triangulated polygons in both cases
419
420 // Create Nodes
421 std::vector<std::unique_ptr<Node>> nodes(points.size());
422 for (const auto i : index_range(points))
423 nodes[i] = Node::build(points[i], /*id*/ i);
424
425 // Build all the sides of the cube
426 std::vector<std::shared_ptr<Polygon>> sides(nodes_on_side.size());
427
428 for (auto s : index_range(nodes_on_side))
429 {
430 const auto & nodes_on_s = nodes_on_side[s];
431 sides[s] = std::make_shared<C0Polygon>(nodes_on_s.size());
432 for (auto i : index_range(nodes_on_s))
433 sides[s]->set_node(i, nodes[nodes_on_s[i]].get());
434 }
435
436 std::unique_ptr<libMesh::Node> mid_elem_node;
437 std::unique_ptr<Elem> polyhedron = std::make_unique<C0Polyhedron>(sides, mid_elem_node);
438 const Point n1 = polyhedron->side_vertex_average_normal(0);
439 LIBMESH_ASSERT_REALVEC_EQUAL(-unit_z, n1, TOLERANCE*TOLERANCE);
440 const Point n2 = polyhedron->side_vertex_average_normal(1);
441 LIBMESH_ASSERT_REALVEC_EQUAL(-unit_y, n2, TOLERANCE*TOLERANCE);
442 const Point n3 = polyhedron->side_vertex_average_normal(2);
443 LIBMESH_ASSERT_REALVEC_EQUAL(unit_x, n3, TOLERANCE*TOLERANCE);
444 const Point n4 = polyhedron->side_vertex_average_normal(3);
445 LIBMESH_ASSERT_REALVEC_EQUAL(unit_y, n4, TOLERANCE*TOLERANCE);
446 const Point n5 = polyhedron->side_vertex_average_normal(4);
447 LIBMESH_ASSERT_REALVEC_EQUAL(-unit_x, n5, TOLERANCE*TOLERANCE);
448 const Point n6 = polyhedron->side_vertex_average_normal(5);
449 LIBMESH_ASSERT_REALVEC_EQUAL(unit_z, n6, TOLERANCE*TOLERANCE);
450 }
451
452 {
453 // Parallepiped with planar faces
454 std::vector<Point> pts = {Point(0, 0, 0), Point(1, 0, 0), Point(1, 1, 0), Point(0, 1, 0),
455 Point(0, 0, 4), Point(1, 0, 4), Point(1, 1, 5), Point(0, 1, 5)};
456
457 // See notes in elem_test.h
458 const std::vector<std::vector<unsigned int>> nodes_on_side =
459 { {0, 1, 2, 3}, // min z
460 {0, 1, 5, 4}, // min y
461 {2, 6, 5, 1}, // max x
462 {2, 3, 7, 6}, // max y
463 {0, 4, 7, 3}, // min x
464 {5, 6, 7, 4} }; // max z
465
466 // Create Nodes
467 std::vector<std::unique_ptr<Node>> nodes(pts.size());
468 for (const auto i : index_range(pts))
469 nodes[i] = Node::build(pts[i], /*id*/ i);
470
471 // Build all the sides of the cube
472 std::vector<std::shared_ptr<Polygon>> sides(nodes_on_side.size());
473
474 for (auto s : index_range(nodes_on_side))
475 {
476 const auto & nodes_on_s = nodes_on_side[s];
477 sides[s] = std::make_shared<C0Polygon>(nodes_on_s.size());
478 for (auto i : index_range(nodes_on_s))
479 sides[s]->set_node(i, nodes[nodes_on_s[i]].get());
480 }
481
482 std::unique_ptr<libMesh::Node> mid_elem_node;
483 std::unique_ptr<Elem> polyhedron = std::make_unique<C0Polyhedron>(sides, mid_elem_node);
484
485 // Get a hex8 for normal comparisons
486 auto [hex8, nodes2] = this->construct_elem(pts, HEX8);
487 for (const auto s : make_range(hex8->n_sides()))
488 {
489 const Point n1 = polyhedron->side_vertex_average_normal(s);
490 const Point normal = hex8->Elem::side_vertex_average_normal(s);
491 LIBMESH_ASSERT_REALVEC_EQUAL(normal, n1, TOLERANCE*TOLERANCE);
492 }
493 }
494 }
495
496protected:
497
498 // Helper function that is called by test_elem_invertible() to build an Elem
499 // of the requested elem_type from the provided Points. Note: the
500 // Nodes which are constructed in order to construct the Elem are
501 // also returned since
502 std::pair<std::unique_ptr<Elem>, std::vector<std::unique_ptr<Node>>>
503 construct_elem(const std::vector<Point> & pts,
504 ElemType elem_type)
505 {
506 const unsigned int n_points = pts.size();
507
508 // Create Nodes
509 std::vector<std::unique_ptr<Node>> nodes(n_points);
510 for (unsigned int i=0; i<n_points; i++)
511 nodes[i] = Node::build(pts[i], /*id*/ i);
512
513 // Create Elem, assign nodes
514 std::unique_ptr<Elem> elem;
515 if (elem_type != C0POLYGON)
516 elem = Elem::build(elem_type, /*parent*/ nullptr);
517 else
518 elem = std::make_unique<C0Polygon>(n_points);
519
520 // Make sure we were passed consistent input to build this type of Elem
521 libmesh_error_msg_if(elem->n_nodes() != n_points,
522 "Wrong number of points "
523 << n_points
524 << " provided to build a "
525 << Utility::enum_to_string(elem_type));
526
527 for (unsigned int i=0; i<n_points; i++)
528 elem->set_node(i, nodes[i].get());
529
530 // Return Elem and Nodes we created
531 return std::make_pair(std::move(elem), std::move(nodes));
532 }
533};
534
535CPPUNIT_TEST_SUITE_REGISTRATION( SideVertexAverageNormalTest );
SideVertexAverageNormalTest::LIBMESH_CPPUNIT_TEST_SUITE
LIBMESH_CPPUNIT_TEST_SUITE(SideVertexAverageNormalTest)
SideVertexAverageNormalTest::construct_elem
std::pair< std::unique_ptr< Elem >, std::vector< std::unique_ptr< Node > > > construct_elem(const std::vector< Point > &pts, ElemType elem_type)
Definition side_vertex_average_normal_test.C:503
SideVertexAverageNormalTest::CPPUNIT_TEST
CPPUNIT_TEST(testEdge3PRefine)
SideVertexAverageNormalTest::CPPUNIT_TEST
CPPUNIT_TEST(testC0Polyhedron)
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition elem.h:96
libMesh::Elem::build
static std::unique_ptr< Elem > build(const ElemType type, Elem *p=nullptr)
Definition elem.C:442
libMesh::Elem::side_vertex_average_normal
virtual Point side_vertex_average_normal(const unsigned int s) const
Definition elem.C:3519
libMesh::FEGenericBase::build
static std::unique_ptr< FEGenericBase > build(const unsigned int dim, const FEType &type)
Builds a specific finite element type.
libMesh::FEType
class FEType hides (possibly multiple) FEFamily and approximation orders, thereby enabling specialize...
Definition fe_type.h:197
libMesh::Node::build
static std::unique_ptr< Node > build(const Node &n)
Definition node.h:315
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition point.h:40
libMesh::QGauss
This class implements specific orders of Gauss quadrature.
Definition quadrature_gauss.h:40
libmesh_cppunit.h
libMesh::ReferenceElem::get
const Elem & get(const ElemType type_in)
Definition reference_elem.C:254
libMesh::Utility::enum_to_string
std::string enum_to_string(const T e)
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
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::ElemType
ElemType
Defines an enum for geometric element types.
Definition enum_elem_type.h:33
libMesh::TOLERANCE
static constexpr Real TOLERANCE
Definition libmesh_common.h:151
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
CPPUNIT_TEST_SUITE_REGISTRATION
CPPUNIT_TEST_SUITE_REGISTRATION(SideVertexAverageNormalTest)
test_comm.h
Generated on Fri Aug 21 2026 13:32:31 for libMesh by doxygen 1.9.8