Merge pull request #4292 from GiudGiud/PR_vertex_ave_2

Move side vertex average normal implementation from edge2 to edge

Author: roystgnr

src/geom/edge_edge2.C

@@ -151,6 +151,7 @@ Edge2::side_vertex_average_normal(const unsigned int s) const
   return (s == 0) ? v : -v;
 }
 
+
 dof_id_type Edge2::key () const
 {
   return this->compute_key(this->node_id(0),

src/geom/reference_elem.C

@@ -178,6 +178,9 @@ void init_ref_elem_table()
   {
     ref_elem_file.clear();
 
+    // 0D elements
+    ref_elem_file[NODEELEM] = ElemDataStrings::one_nodeelem;
+
     // 1D elements
     ref_elem_file[EDGE2]    = ElemDataStrings::one_edge;
     ref_elem_file[EDGE3]    = ElemDataStrings::one_edge3;

src/geom/reference_elem.data

@@ -6,6 +6,20 @@
   typedef std::map<ElemType, std::string> FileMapType;
   FileMapType ref_elem_file;
 
+static const char* one_nodeelem =
+"libMesh-0.7.0+\n"
+"1	 # number of elements\n"
+"1	 # number of nodes\n"
+".	 # boundary condition specification file\n"
+"n/a	 # subdomain id specification file\n"
+"n/a	 # processor id specification file\n"
+"n/a	 # p-level specification file\n"
+"1	 # n_elem at level 0, [ type (n0 ... nN-1) ]\n"
+"27 0\n"
+"0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00\n"
+"1	 # number of boundary conditions\n"
+"0 0 0\n"
+;
 
 static const char* one_edge =
 "libMesh-0.7.0+\n"

tests/geom/side_vertex_average_normal_test.C

@@ -30,6 +30,7 @@ class SideVertexAverageNormalTest : public CppUnit::TestCase
 public:
   LIBMESH_CPPUNIT_TEST_SUITE( SideVertexAverageNormalTest );
   CPPUNIT_TEST( testEdge2 );
+  CPPUNIT_TEST( testEdge3 );
   CPPUNIT_TEST( testTri3 );
   CPPUNIT_TEST( testQuad4 );
   CPPUNIT_TEST( testPyramid5 );
@@ -79,6 +80,43 @@ public:
     }
   }
 
+  void testEdge3()
+  {
+    LOG_UNIT_TEST;
+
+    {
+      // Reference
+      const Elem & edge3 = ReferenceElem::get(EDGE3);
+      const Point n1 = edge3.side_vertex_average_normal(0);
+      LIBMESH_ASSERT_FP_EQUAL(-1, n1(0), TOLERANCE*TOLERANCE);
+      LIBMESH_ASSERT_FP_EQUAL(0, n1(1), TOLERANCE*TOLERANCE);
+      LIBMESH_ASSERT_FP_EQUAL(0, n1(2), TOLERANCE*TOLERANCE);
+      const Point n2 = edge3.side_vertex_average_normal(1);
+      LIBMESH_ASSERT_FP_EQUAL(1, n2(0), TOLERANCE*TOLERANCE);
+      LIBMESH_ASSERT_FP_EQUAL(0, n2(1), TOLERANCE*TOLERANCE);
+      LIBMESH_ASSERT_FP_EQUAL(0, n2(2), TOLERANCE*TOLERANCE);
+    }
+
+    {
+      // Oriented, checked with the FE construction
+      std::vector<Point> pts = {Point(1, 0, 0), Point(1, 3, 0), Point(2.2344, 1.210293, 0)};
+      auto [edge3, nodes] = this->construct_elem(pts, EDGE3);
+      std::unique_ptr<libMesh::FEBase> fe(libMesh::FEBase::build(1, libMesh::FEType(1)));
+      libMesh::QGauss qface(0, libMesh::CONSTANT);
+      const std::vector<Point> & normals = fe->get_normals();
+      for (const auto s : make_range(edge3->n_sides()))
+      {
+        const std::unique_ptr<const Elem> face = edge3->build_side_ptr(s);
+        fe->attach_quadrature_rule(&qface);
+        fe->reinit(edge3.get(), s, TOLERANCE);
+        const Point n1 = edge3->side_vertex_average_normal(s);
+        LIBMESH_ASSERT_FP_EQUAL(normals[0](0), n1(0), TOLERANCE*TOLERANCE);
+        LIBMESH_ASSERT_FP_EQUAL(normals[0](1), n1(1), TOLERANCE*TOLERANCE);
+        LIBMESH_ASSERT_FP_EQUAL(normals[0](2), n1(2), TOLERANCE*TOLERANCE);
+      }
+    }
+  }
+
   void testTri3()
   {
     LOG_UNIT_TEST;