diff --git a/include/geom/cell_prism.h b/include/geom/cell_prism.h
index fbe7688280..20bfa80560 100644
--- a/include/geom/cell_prism.h
+++ b/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
diff --git a/include/geom/cell_pyramid.h b/include/geom/cell_pyramid.h
index fe347d1953..7f8a9c4f39 100644
--- a/include/geom/cell_pyramid.h
+++ b/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
 
diff --git a/include/geom/face_tri.h b/include/geom/face_tri.h
index bfc249a0cb..595ee04014 100644
--- a/include/geom/face_tri.h
+++ b/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
diff --git a/src/fe/fe_abstract.C b/src/fe/fe_abstract.C
index 20d435ae44..741e0cc089 100644
--- a/src/fe/fe_abstract.C
+++ b/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);
 }
 
 
diff --git a/src/geom/cell_prism.C b/src/geom/cell_prism.C
index 618d681051..180720f34c 100644
--- a/src/geom/cell_prism.C
+++ b/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] =
diff --git a/src/geom/cell_pyramid.C b/src/geom/cell_pyramid.C
index 04bdf5298a..1169693a8d 100644
--- a/src/geom/cell_pyramid.C
+++ b/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] =
diff --git a/src/geom/face_tri.C b/src/geom/face_tri.C
index f63a640954..b6ebd6a879 100644
--- a/src/geom/face_tri.C
+++ b/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] =
diff --git a/tests/fe/fe_test.h b/tests/fe/fe_test.h
index 468d982d5f..d2a6e49836 100644
--- a/tests/fe/fe_test.h
+++ b/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;
 
diff --git a/tests/geom/elem_test.C b/tests/geom/elem_test.C
index 41ced6f6b7..943e4f8554 100644
--- a/tests/geom/elem_test.C
+++ b/tests/geom/elem_test.C
@@ -50,6 +50,31 @@ public:
       }
   }
 
+  void test_ref_elem()
+  {
+    LOG_UNIT_TEST;
+
+    for (const auto & elem : this->_mesh->active_local_element_ptr_range())
+      {
+        // We don't currently support reference infinite elements
+        if (elem->infinite())
+          continue;
+
+        // We'll probably never support reference polytopes
+        if (elem->runtime_topology())
+          continue;
+
+        const Elem & ref_elem = *elem->reference_elem();
+
+        for (auto i : elem->node_index_range())
+          {
+            LIBMESH_ASSERT_REALVEC_EQUAL(elem->master_point(i),
+                                         ref_elem.point(i),
+                                         TOLERANCE*TOLERANCE);
+          }
+      }
+  }
+
   void test_quality()
   {
     LOG_UNIT_TEST;
@@ -944,6 +969,7 @@ public:
 
 #define ELEMTEST                                \
   CPPUNIT_TEST( test_bounding_box );            \
+  CPPUNIT_TEST( test_ref_elem );                \
   CPPUNIT_TEST( test_quality );                 \
   CPPUNIT_TEST( test_node_edge_map_consistency ); \
   CPPUNIT_TEST( test_maps );                    \