Refspace fixes

include/geom/cell_prism.h

@@ -211,7 +211,7 @@ class Prism : public Cell
   /**
    * Master element node locations
    */
-  static const Real _master_points[18][3];
+  static const Real _master_points[21][3];
 
   /**
    * This maps the \f$ j^{th} \f$ Prism vertex to the (in this case) 3

include/geom/cell_pyramid.h

@@ -211,7 +211,7 @@ class Pyramid : public Cell
   /**
    * Master element node locations
    */
-  static const Real _master_points[14][3];
+  static const Real _master_points[18][3];
 
 #ifdef LIBMESH_ENABLE_AMR
 

include/geom/face_tri.h

@@ -205,7 +205,7 @@ class Tri : public Face
   /**
    * Master element node locations
    */
-  static const Real _master_points[6][3];
+  static const Real _master_points[7][3];
 
   /**
    * This maps the \f$ j^{th} \f$ node to the (in this case) 2 side

src/fe/fe_abstract.C

@@ -31,6 +31,7 @@
 #include "libmesh/periodic_boundary.h"
 #include "libmesh/quadrature.h"
 #include "libmesh/quadrature_gauss.h"
+#include "libmesh/reference_elem.h"
 #include "libmesh/remote_elem.h"
 #include "libmesh/tensor_value.h"
 #include "libmesh/threads.h"
@@ -405,230 +406,13 @@ void FEAbstract::get_refspace_nodes(const ElemType itemType, std::vector<Point>
                       Utility::enum_to_string(itemType));
 
   nodes.resize(n_nodes);
-  switch(itemType)
-    {
-    case NODEELEM:
-      {
-        nodes[0] = Point (0.,0.,0.);
-        return;
-      }
-    case EDGE3:
-      {
-        nodes[2] = Point (0.,0.,0.);
-        libmesh_fallthrough();
-      }
-    case EDGE2:
-      {
-        nodes[0] = Point (-1.,0.,0.);
-        nodes[1] = Point (1.,0.,0.);
-        return;
-      }
-    case EDGE4: // not nested with EDGE3
-      {
-        nodes[0] = Point (-1.,0.,0.);
-        nodes[1] = Point (1.,0.,0.);
-        nodes[2] = Point (-1./3.,0.,0.);
-        nodes[3] - Point (1./3.,0.,0.);
-        return;
-      }
-    case TRI7:
-      {
-        nodes[6] = Point (1./3.,1./3.,0.);
-        libmesh_fallthrough();
-      }
-    case TRI6:
-      {
-        nodes[3] = Point (.5,0.,0.);
-        nodes[4] = Point (.5,.5,0.);
-        nodes[5] = Point (0.,.5,0.);
-        libmesh_fallthrough();
-      }
-    case TRI3:
-    case TRISHELL3:
-      {
-        nodes[0] = Point (0.,0.,0.);
-        nodes[1] = Point (1.,0.,0.);
-        nodes[2] = Point (0.,1.,0.);
-        return;
-      }
-    case QUAD9:
-    case QUADSHELL9:
-      {
-        nodes[8] = Point (0.,0.,0.);
-        libmesh_fallthrough();
-      }
-    case QUAD8:
-    case QUADSHELL8:
-      {
-        nodes[4] = Point (0.,-1.,0.);
-        nodes[5] = Point (1.,0.,0.);
-        nodes[6] = Point (0.,1.,0.);
-        nodes[7] = Point (-1.,0.,0.);
-        libmesh_fallthrough();
-      }
-    case QUAD4:
-    case QUADSHELL4:
-      {
-        nodes[0] = Point (-1.,-1.,0.);
-        nodes[1] = Point (1.,-1.,0.);
-        nodes[2] = Point (1.,1.,0.);
-        nodes[3] = Point (-1.,1.,0.);
-        return;
-      }
-    case TET14:
-      {
-        nodes[10] = Point (1/Real(3),1/Real(3),0.);
-        nodes[11] = Point (1/Real(3),0.,1/Real(3));
-        nodes[12] = Point (1/Real(3),1/Real(3),1/Real(3));
-        nodes[13] = Point (0.,1/Real(3),1/Real(3));
-        libmesh_fallthrough();
-      }
-    case TET10:
-      {
-        nodes[4] = Point (.5,0.,0.);
-        nodes[5] = Point (.5,.5,0.);
-        nodes[6] = Point (0.,.5,0.);
-        nodes[7] = Point (0.,0.,.5);
-        nodes[8] = Point (.5,0.,.5);
-        nodes[9] = Point (0.,.5,.5);
-        libmesh_fallthrough();
-      }
-    case TET4:
-      {
-        nodes[0] = Point (0.,0.,0.);
-        nodes[1] = Point (1.,0.,0.);
-        nodes[2] = Point (0.,1.,0.);
-        nodes[3] = Point (0.,0.,1.);
-        return;
-      }
-    case HEX27:
-      {
-        nodes[20] = Point (0.,0.,-1.);
-        nodes[21] = Point (0.,-1.,0.);
-        nodes[22] = Point (1.,0.,0.);
-        nodes[23] = Point (0.,1.,0.);
-        nodes[24] = Point (-1.,0.,0.);
-        nodes[25] = Point (0.,0.,1.);
-        nodes[26] = Point (0.,0.,0.);
-        libmesh_fallthrough();
-      }
-    case HEX20:
-      {
-        nodes[8] = Point (0.,-1.,-1.);
-        nodes[9] = Point (1.,0.,-1.);
-        nodes[10] = Point (0.,1.,-1.);
-        nodes[11] = Point (-1.,0.,-1.);
-        nodes[12] = Point (-1.,-1.,0.);
-        nodes[13] = Point (1.,-1.,0.);
-        nodes[14] = Point (1.,1.,0.);
-        nodes[15] = Point (-1.,1.,0.);
-        nodes[16] = Point (0.,-1.,1.);
-        nodes[17] = Point (1.,0.,1.);
-        nodes[18] = Point (0.,1.,1.);
-        nodes[19] = Point (-1.,0.,1.);
-        libmesh_fallthrough();
-      }
-    case HEX8:
-      {
-        nodes[0] = Point (-1.,-1.,-1.);
-        nodes[1] = Point (1.,-1.,-1.);
-        nodes[2] = Point (1.,1.,-1.);
-        nodes[3] = Point (-1.,1.,-1.);
-        nodes[4] = Point (-1.,-1.,1.);
-        nodes[5] = Point (1.,-1.,1.);
-        nodes[6] = Point (1.,1.,1.);
-        nodes[7] = Point (-1.,1.,1.);
-        return;
-      }
-    case PRISM21:
-      {
-        nodes[20] = Point (1/Real(3),1/Real(3),0);
-        libmesh_fallthrough();
-      }
-    case PRISM20:
-      {
-        nodes[18] = Point (1/Real(3),1/Real(3),-1);
-        nodes[19] = Point (1/Real(3),1/Real(3),1);
-        libmesh_fallthrough();
-      }
-    case PRISM18:
-      {
-        nodes[15] = Point (.5,0.,0.);
-        nodes[16] = Point (.5,.5,0.);
-        nodes[17] = Point (0.,.5,0.);
-        libmesh_fallthrough();
-      }
-    case PRISM15:
-      {
-        nodes[6] = Point (.5,0.,-1.);
-        nodes[7] = Point (.5,.5,-1.);
-        nodes[8] = Point (0.,.5,-1.);
-        nodes[9] = Point (0.,0.,0.);
-        nodes[10] = Point (1.,0.,0.);
-        nodes[11] = Point (0.,1.,0.);
-        nodes[12] = Point (.5,0.,1.);
-        nodes[13] = Point (.5,.5,1.);
-        nodes[14] = Point (0.,.5,1.);
-        libmesh_fallthrough();
-      }
-    case PRISM6:
-      {
-        nodes[0] = Point (0.,0.,-1.);
-        nodes[1] = Point (1.,0.,-1.);
-        nodes[2] = Point (0.,1.,-1.);
-        nodes[3] = Point (0.,0.,1.);
-        nodes[4] = Point (1.,0.,1.);
-        nodes[5] = Point (0.,1.,1.);
-        return;
-      }
-    case PYRAMID18:
-      {
-        // triangle centers
-        nodes[14] = Point (-2/Real(3),0.,1/Real(3));
-        nodes[15] = Point (0.,2/Real(3),1/Real(3));
-        nodes[16] = Point (2/Real(3),0.,1/Real(3));
-        nodes[17] = Point (0.,-2/Real(3),1/Real(3));
-
-        libmesh_fallthrough();
-      }
-    case PYRAMID14:
-      {
-        // base center
-        nodes[13] = Point (0.,0.,0.);
 
-        libmesh_fallthrough();
-      }
-    case PYRAMID13:
-      {
-        // base midedge
-        nodes[5] = Point (0.,-1.,0.);
-        nodes[6] = Point (1.,0.,0.);
-        nodes[7] = Point (0.,1.,0.);
-        nodes[8] = Point (-1,0.,0.);
-
-        // lateral midedge
-        nodes[9] = Point (-.5,-.5,.5);
-        nodes[10] = Point (.5,-.5,.5);
-        nodes[11] = Point (.5,.5,.5);
-        nodes[12] = Point (-.5,.5,.5);
-
-        libmesh_fallthrough();
-      }
-    case PYRAMID5:
-      {
-        // base corners
-        nodes[0] = Point (-1.,-1.,0.);
-        nodes[1] = Point (1.,-1.,0.);
-        nodes[2] = Point (1.,1.,0.);
-        nodes[3] = Point (-1.,1.,0.);
-        // apex
-        nodes[4] = Point (0.,0.,1.);
-        return;
-      }
+  // For any element with a well-defined n_nodes we should also have a
+  // well-defined reference element.
+  const Elem & ref_elem = ReferenceElem::get(itemType);
 
-    default:
-      libmesh_error_msg("ERROR: Unknown element type " << Utility::enum_to_string(itemType));
-    }
+  for (auto i : make_range(n_nodes))
+    nodes[i] = ref_elem.point(i);
 }
 
 

src/geom/cell_prism.C

@@ -34,7 +34,7 @@ const int Prism::num_sides;
 const int Prism::num_edges;
 const int Prism::num_children;
 
-const Real Prism::_master_points[18][3] =
+const Real Prism::_master_points[21][3] =
   {
     {0, 0, -1},
     {1, 0, -1},
@@ -53,7 +53,10 @@ const Real Prism::_master_points[18][3] =
     {0, 0.5, 1},
     {0.5, 0, 0},
     {0.5, 0.5, 0},
-    {0, 0.5, 0}
+    {0, 0.5, 0},
+    {Real(1)/3, Real(1)/3, -1},
+    {Real(1)/3, Real(1)/3, 1},
+    {Real(1)/3, Real(1)/3, 0}
   };
 
 const unsigned int Prism::edge_sides_map[9][2] =

src/geom/cell_pyramid.C

@@ -33,7 +33,7 @@ const int Pyramid::num_sides;
 const int Pyramid::num_edges;
 const int Pyramid::num_children;
 
-const Real Pyramid::_master_points[14][3] =
+const Real Pyramid::_master_points[18][3] =
   {
     {-1, -1, 0},
     {1, -1, 0},
@@ -44,11 +44,15 @@ const Real Pyramid::_master_points[14][3] =
     {1, 0, 0},
     {0, 1, 0},
     {-1, 0, 0},
-    {0, -0.5, 0.5},
-    {0.5, 0, 0.5},
-    {0, 0.5, 0.5},
-    {-0.5, 0, 0.5},
-    {0, 0, 0}
+    {-0.5, -0.5, 0.5},
+    {0.5, -0.5, 0.5},
+    {0.5, 0.5, 0.5},
+    {-0.5, 0.5, 0.5},
+    {0, 0, 0},
+    {0, -Real(2)/3, Real(1)/3},
+    { Real(2)/3, 0, Real(1)/3},
+    {0,  Real(2)/3, Real(1)/3},
+    {-Real(2)/3, 0, Real(1)/3}
   };
 
 const unsigned int Pyramid::edge_sides_map[8][2] =

src/geom/face_tri.C

@@ -35,14 +35,15 @@ const int Tri::num_children;
 
 // Note: we can omit initialization of the third entry of each row because
 // static variables are automatically zero-initialized.
-const Real Tri::_master_points[6][3] =
+const Real Tri::_master_points[7][3] =
   {
     {0, 0},
     {1, 0},
     {0, 1},
     {0.5, 0},
     {0.5, 0.5},
-    {0, 0.5}
+    {0, 0.5},
+    {Real(1)/3, Real(1)/3}
   };
 
 const unsigned int Tri::adjacent_sides_map[/*num_vertices*/3][/*n_adjacent_sides*/2] =

tests/fe/fe_test.h

@@ -25,6 +25,7 @@
 
 #define FETEST                                  \
   CPPUNIT_TEST( testFEInterface );              \
+  CPPUNIT_TEST( testRefspaceNodes );            \
   CPPUNIT_TEST( testU );                        \
   CPPUNIT_TEST( testPartitionOfUnity );         \
   CPPUNIT_TEST( testGradU );                    \
@@ -668,8 +669,33 @@ class FETest : public FETestBase<order, family, elem_type, 1, CaseName> {
 #endif // LIBMESH_ENABLE_EXCEPTIONS
   }
 
+  void testRefspaceNodes()
+    {
+      LOG_UNIT_TEST;
+
+      if (!this->_elem)
+        return;
+
+      // We don't currently support reference infinite elements
+      if (this->_elem->infinite())
+        return;
+
+      // We'll probably never support reference polytopes
+      if (this->_elem->runtime_topology())
+        return;
+
+      std::vector<Point> nodes;
+      this->_fe->get_refspace_nodes(this->_elem->type(), nodes);
+      for (auto n : index_range(nodes))
+        LIBMESH_ASSERT_REALVEC_EQUAL(nodes[n],
+                                     this->_elem->master_point(n),
+                                     TOLERANCE*TOLERANCE);
+    }
+
   void testPartitionOfUnity()
     {
+      LOG_UNIT_TEST;
+
       if (!this->_elem)
         return;