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