diff --git a/Makefile.in b/Makefile.in
index 043fb445fa8..615a533febc 100644
--- a/Makefile.in
+++ b/Makefile.in
@@ -341,12 +341,12 @@ am__libmesh_dbg_la_SOURCES_DIST = src/base/default_coupling.C \
 	src/geom/cell_prism20.C src/geom/cell_prism21.C \
 	src/geom/cell_prism6.C src/geom/cell_pyramid.C \
 	src/geom/cell_pyramid13.C src/geom/cell_pyramid14.C \
-	src/geom/cell_pyramid5.C src/geom/cell_tet.C \
-	src/geom/cell_tet10.C src/geom/cell_tet14.C \
-	src/geom/cell_tet4.C src/geom/edge.C src/geom/edge_edge2.C \
-	src/geom/edge_edge3.C src/geom/edge_edge4.C \
-	src/geom/edge_inf_edge2.C src/geom/elem.C \
-	src/geom/elem_cutter.C src/geom/elem_quality.C \
+	src/geom/cell_pyramid18.C src/geom/cell_pyramid5.C \
+	src/geom/cell_tet.C src/geom/cell_tet10.C \
+	src/geom/cell_tet14.C src/geom/cell_tet4.C src/geom/edge.C \
+	src/geom/edge_edge2.C src/geom/edge_edge3.C \
+	src/geom/edge_edge4.C src/geom/edge_inf_edge2.C \
+	src/geom/elem.C src/geom/elem_cutter.C src/geom/elem_quality.C \
 	src/geom/elem_refinement.C src/geom/elem_side_builder.C \
 	src/geom/face.C src/geom/face_inf_quad.C \
 	src/geom/face_inf_quad4.C src/geom/face_inf_quad6.C \
@@ -704,6 +704,7 @@ am__objects_1 = src/base/libmesh_dbg_la-default_coupling.lo \
 	src/geom/libmesh_dbg_la-cell_pyramid.lo \
 	src/geom/libmesh_dbg_la-cell_pyramid13.lo \
 	src/geom/libmesh_dbg_la-cell_pyramid14.lo \
+	src/geom/libmesh_dbg_la-cell_pyramid18.lo \
 	src/geom/libmesh_dbg_la-cell_pyramid5.lo \
 	src/geom/libmesh_dbg_la-cell_tet.lo \
 	src/geom/libmesh_dbg_la-cell_tet10.lo \
@@ -1102,12 +1103,12 @@ am__libmesh_devel_la_SOURCES_DIST = src/base/default_coupling.C \
 	src/geom/cell_prism20.C src/geom/cell_prism21.C \
 	src/geom/cell_prism6.C src/geom/cell_pyramid.C \
 	src/geom/cell_pyramid13.C src/geom/cell_pyramid14.C \
-	src/geom/cell_pyramid5.C src/geom/cell_tet.C \
-	src/geom/cell_tet10.C src/geom/cell_tet14.C \
-	src/geom/cell_tet4.C src/geom/edge.C src/geom/edge_edge2.C \
-	src/geom/edge_edge3.C src/geom/edge_edge4.C \
-	src/geom/edge_inf_edge2.C src/geom/elem.C \
-	src/geom/elem_cutter.C src/geom/elem_quality.C \
+	src/geom/cell_pyramid18.C src/geom/cell_pyramid5.C \
+	src/geom/cell_tet.C src/geom/cell_tet10.C \
+	src/geom/cell_tet14.C src/geom/cell_tet4.C src/geom/edge.C \
+	src/geom/edge_edge2.C src/geom/edge_edge3.C \
+	src/geom/edge_edge4.C src/geom/edge_inf_edge2.C \
+	src/geom/elem.C src/geom/elem_cutter.C src/geom/elem_quality.C \
 	src/geom/elem_refinement.C src/geom/elem_side_builder.C \
 	src/geom/face.C src/geom/face_inf_quad.C \
 	src/geom/face_inf_quad4.C src/geom/face_inf_quad6.C \
@@ -1464,6 +1465,7 @@ am__objects_2 = src/base/libmesh_devel_la-default_coupling.lo \
 	src/geom/libmesh_devel_la-cell_pyramid.lo \
 	src/geom/libmesh_devel_la-cell_pyramid13.lo \
 	src/geom/libmesh_devel_la-cell_pyramid14.lo \
+	src/geom/libmesh_devel_la-cell_pyramid18.lo \
 	src/geom/libmesh_devel_la-cell_pyramid5.lo \
 	src/geom/libmesh_devel_la-cell_tet.lo \
 	src/geom/libmesh_devel_la-cell_tet10.lo \
@@ -1859,12 +1861,12 @@ am__libmesh_oprof_la_SOURCES_DIST = src/base/default_coupling.C \
 	src/geom/cell_prism20.C src/geom/cell_prism21.C \
 	src/geom/cell_prism6.C src/geom/cell_pyramid.C \
 	src/geom/cell_pyramid13.C src/geom/cell_pyramid14.C \
-	src/geom/cell_pyramid5.C src/geom/cell_tet.C \
-	src/geom/cell_tet10.C src/geom/cell_tet14.C \
-	src/geom/cell_tet4.C src/geom/edge.C src/geom/edge_edge2.C \
-	src/geom/edge_edge3.C src/geom/edge_edge4.C \
-	src/geom/edge_inf_edge2.C src/geom/elem.C \
-	src/geom/elem_cutter.C src/geom/elem_quality.C \
+	src/geom/cell_pyramid18.C src/geom/cell_pyramid5.C \
+	src/geom/cell_tet.C src/geom/cell_tet10.C \
+	src/geom/cell_tet14.C src/geom/cell_tet4.C src/geom/edge.C \
+	src/geom/edge_edge2.C src/geom/edge_edge3.C \
+	src/geom/edge_edge4.C src/geom/edge_inf_edge2.C \
+	src/geom/elem.C src/geom/elem_cutter.C src/geom/elem_quality.C \
 	src/geom/elem_refinement.C src/geom/elem_side_builder.C \
 	src/geom/face.C src/geom/face_inf_quad.C \
 	src/geom/face_inf_quad4.C src/geom/face_inf_quad6.C \
@@ -2221,6 +2223,7 @@ am__objects_3 = src/base/libmesh_oprof_la-default_coupling.lo \
 	src/geom/libmesh_oprof_la-cell_pyramid.lo \
 	src/geom/libmesh_oprof_la-cell_pyramid13.lo \
 	src/geom/libmesh_oprof_la-cell_pyramid14.lo \
+	src/geom/libmesh_oprof_la-cell_pyramid18.lo \
 	src/geom/libmesh_oprof_la-cell_pyramid5.lo \
 	src/geom/libmesh_oprof_la-cell_tet.lo \
 	src/geom/libmesh_oprof_la-cell_tet10.lo \
@@ -2616,12 +2619,12 @@ am__libmesh_opt_la_SOURCES_DIST = src/base/default_coupling.C \
 	src/geom/cell_prism20.C src/geom/cell_prism21.C \
 	src/geom/cell_prism6.C src/geom/cell_pyramid.C \
 	src/geom/cell_pyramid13.C src/geom/cell_pyramid14.C \
-	src/geom/cell_pyramid5.C src/geom/cell_tet.C \
-	src/geom/cell_tet10.C src/geom/cell_tet14.C \
-	src/geom/cell_tet4.C src/geom/edge.C src/geom/edge_edge2.C \
-	src/geom/edge_edge3.C src/geom/edge_edge4.C \
-	src/geom/edge_inf_edge2.C src/geom/elem.C \
-	src/geom/elem_cutter.C src/geom/elem_quality.C \
+	src/geom/cell_pyramid18.C src/geom/cell_pyramid5.C \
+	src/geom/cell_tet.C src/geom/cell_tet10.C \
+	src/geom/cell_tet14.C src/geom/cell_tet4.C src/geom/edge.C \
+	src/geom/edge_edge2.C src/geom/edge_edge3.C \
+	src/geom/edge_edge4.C src/geom/edge_inf_edge2.C \
+	src/geom/elem.C src/geom/elem_cutter.C src/geom/elem_quality.C \
 	src/geom/elem_refinement.C src/geom/elem_side_builder.C \
 	src/geom/face.C src/geom/face_inf_quad.C \
 	src/geom/face_inf_quad4.C src/geom/face_inf_quad6.C \
@@ -2978,6 +2981,7 @@ am__objects_4 = src/base/libmesh_opt_la-default_coupling.lo \
 	src/geom/libmesh_opt_la-cell_pyramid.lo \
 	src/geom/libmesh_opt_la-cell_pyramid13.lo \
 	src/geom/libmesh_opt_la-cell_pyramid14.lo \
+	src/geom/libmesh_opt_la-cell_pyramid18.lo \
 	src/geom/libmesh_opt_la-cell_pyramid5.lo \
 	src/geom/libmesh_opt_la-cell_tet.lo \
 	src/geom/libmesh_opt_la-cell_tet10.lo \
@@ -3372,12 +3376,12 @@ am__libmesh_prof_la_SOURCES_DIST = src/base/default_coupling.C \
 	src/geom/cell_prism20.C src/geom/cell_prism21.C \
 	src/geom/cell_prism6.C src/geom/cell_pyramid.C \
 	src/geom/cell_pyramid13.C src/geom/cell_pyramid14.C \
-	src/geom/cell_pyramid5.C src/geom/cell_tet.C \
-	src/geom/cell_tet10.C src/geom/cell_tet14.C \
-	src/geom/cell_tet4.C src/geom/edge.C src/geom/edge_edge2.C \
-	src/geom/edge_edge3.C src/geom/edge_edge4.C \
-	src/geom/edge_inf_edge2.C src/geom/elem.C \
-	src/geom/elem_cutter.C src/geom/elem_quality.C \
+	src/geom/cell_pyramid18.C src/geom/cell_pyramid5.C \
+	src/geom/cell_tet.C src/geom/cell_tet10.C \
+	src/geom/cell_tet14.C src/geom/cell_tet4.C src/geom/edge.C \
+	src/geom/edge_edge2.C src/geom/edge_edge3.C \
+	src/geom/edge_edge4.C src/geom/edge_inf_edge2.C \
+	src/geom/elem.C src/geom/elem_cutter.C src/geom/elem_quality.C \
 	src/geom/elem_refinement.C src/geom/elem_side_builder.C \
 	src/geom/face.C src/geom/face_inf_quad.C \
 	src/geom/face_inf_quad4.C src/geom/face_inf_quad6.C \
@@ -3734,6 +3738,7 @@ am__objects_5 = src/base/libmesh_prof_la-default_coupling.lo \
 	src/geom/libmesh_prof_la-cell_pyramid.lo \
 	src/geom/libmesh_prof_la-cell_pyramid13.lo \
 	src/geom/libmesh_prof_la-cell_pyramid14.lo \
+	src/geom/libmesh_prof_la-cell_pyramid18.lo \
 	src/geom/libmesh_prof_la-cell_pyramid5.lo \
 	src/geom/libmesh_prof_la-cell_tet.lo \
 	src/geom/libmesh_prof_la-cell_tet10.lo \
@@ -5093,6 +5098,7 @@ am__depfiles_remade = src/apps/$(DEPDIR)/amr_dbg-amr.Po \
 	src/geom/$(DEPDIR)/libmesh_dbg_la-cell_pyramid.Plo \
 	src/geom/$(DEPDIR)/libmesh_dbg_la-cell_pyramid13.Plo \
 	src/geom/$(DEPDIR)/libmesh_dbg_la-cell_pyramid14.Plo \
+	src/geom/$(DEPDIR)/libmesh_dbg_la-cell_pyramid18.Plo \
 	src/geom/$(DEPDIR)/libmesh_dbg_la-cell_pyramid5.Plo \
 	src/geom/$(DEPDIR)/libmesh_dbg_la-cell_tet.Plo \
 	src/geom/$(DEPDIR)/libmesh_dbg_la-cell_tet10.Plo \
@@ -5155,6 +5161,7 @@ am__depfiles_remade = src/apps/$(DEPDIR)/amr_dbg-amr.Po \
 	src/geom/$(DEPDIR)/libmesh_devel_la-cell_pyramid.Plo \
 	src/geom/$(DEPDIR)/libmesh_devel_la-cell_pyramid13.Plo \
 	src/geom/$(DEPDIR)/libmesh_devel_la-cell_pyramid14.Plo \
+	src/geom/$(DEPDIR)/libmesh_devel_la-cell_pyramid18.Plo \
 	src/geom/$(DEPDIR)/libmesh_devel_la-cell_pyramid5.Plo \
 	src/geom/$(DEPDIR)/libmesh_devel_la-cell_tet.Plo \
 	src/geom/$(DEPDIR)/libmesh_devel_la-cell_tet10.Plo \
@@ -5217,6 +5224,7 @@ am__depfiles_remade = src/apps/$(DEPDIR)/amr_dbg-amr.Po \
 	src/geom/$(DEPDIR)/libmesh_oprof_la-cell_pyramid.Plo \
 	src/geom/$(DEPDIR)/libmesh_oprof_la-cell_pyramid13.Plo \
 	src/geom/$(DEPDIR)/libmesh_oprof_la-cell_pyramid14.Plo \
+	src/geom/$(DEPDIR)/libmesh_oprof_la-cell_pyramid18.Plo \
 	src/geom/$(DEPDIR)/libmesh_oprof_la-cell_pyramid5.Plo \
 	src/geom/$(DEPDIR)/libmesh_oprof_la-cell_tet.Plo \
 	src/geom/$(DEPDIR)/libmesh_oprof_la-cell_tet10.Plo \
@@ -5279,6 +5287,7 @@ am__depfiles_remade = src/apps/$(DEPDIR)/amr_dbg-amr.Po \
 	src/geom/$(DEPDIR)/libmesh_opt_la-cell_pyramid.Plo \
 	src/geom/$(DEPDIR)/libmesh_opt_la-cell_pyramid13.Plo \
 	src/geom/$(DEPDIR)/libmesh_opt_la-cell_pyramid14.Plo \
+	src/geom/$(DEPDIR)/libmesh_opt_la-cell_pyramid18.Plo \
 	src/geom/$(DEPDIR)/libmesh_opt_la-cell_pyramid5.Plo \
 	src/geom/$(DEPDIR)/libmesh_opt_la-cell_tet.Plo \
 	src/geom/$(DEPDIR)/libmesh_opt_la-cell_tet10.Plo \
@@ -5341,6 +5350,7 @@ am__depfiles_remade = src/apps/$(DEPDIR)/amr_dbg-amr.Po \
 	src/geom/$(DEPDIR)/libmesh_prof_la-cell_pyramid.Plo \
 	src/geom/$(DEPDIR)/libmesh_prof_la-cell_pyramid13.Plo \
 	src/geom/$(DEPDIR)/libmesh_prof_la-cell_pyramid14.Plo \
+	src/geom/$(DEPDIR)/libmesh_prof_la-cell_pyramid18.Plo \
 	src/geom/$(DEPDIR)/libmesh_prof_la-cell_pyramid5.Plo \
 	src/geom/$(DEPDIR)/libmesh_prof_la-cell_tet.Plo \
 	src/geom/$(DEPDIR)/libmesh_prof_la-cell_tet10.Plo \
@@ -7497,6 +7507,7 @@ libmesh_SOURCES = \
         src/geom/cell_pyramid.C \
         src/geom/cell_pyramid13.C \
         src/geom/cell_pyramid14.C \
+        src/geom/cell_pyramid18.C \
         src/geom/cell_pyramid5.C \
         src/geom/cell_tet.C \
         src/geom/cell_tet10.C \
@@ -8567,6 +8578,8 @@ src/geom/libmesh_dbg_la-cell_pyramid13.lo: src/geom/$(am__dirstamp) \
 	src/geom/$(DEPDIR)/$(am__dirstamp)
 src/geom/libmesh_dbg_la-cell_pyramid14.lo: src/geom/$(am__dirstamp) \
 	src/geom/$(DEPDIR)/$(am__dirstamp)
+src/geom/libmesh_dbg_la-cell_pyramid18.lo: src/geom/$(am__dirstamp) \
+	src/geom/$(DEPDIR)/$(am__dirstamp)
 src/geom/libmesh_dbg_la-cell_pyramid5.lo: src/geom/$(am__dirstamp) \
 	src/geom/$(DEPDIR)/$(am__dirstamp)
 src/geom/libmesh_dbg_la-cell_tet.lo: src/geom/$(am__dirstamp) \
@@ -9730,6 +9743,8 @@ src/geom/libmesh_devel_la-cell_pyramid13.lo: src/geom/$(am__dirstamp) \
 	src/geom/$(DEPDIR)/$(am__dirstamp)
 src/geom/libmesh_devel_la-cell_pyramid14.lo: src/geom/$(am__dirstamp) \
 	src/geom/$(DEPDIR)/$(am__dirstamp)
+src/geom/libmesh_devel_la-cell_pyramid18.lo: src/geom/$(am__dirstamp) \
+	src/geom/$(DEPDIR)/$(am__dirstamp)
 src/geom/libmesh_devel_la-cell_pyramid5.lo: src/geom/$(am__dirstamp) \
 	src/geom/$(DEPDIR)/$(am__dirstamp)
 src/geom/libmesh_devel_la-cell_tet.lo: src/geom/$(am__dirstamp) \
@@ -10830,6 +10845,8 @@ src/geom/libmesh_oprof_la-cell_pyramid13.lo: src/geom/$(am__dirstamp) \
 	src/geom/$(DEPDIR)/$(am__dirstamp)
 src/geom/libmesh_oprof_la-cell_pyramid14.lo: src/geom/$(am__dirstamp) \
 	src/geom/$(DEPDIR)/$(am__dirstamp)
+src/geom/libmesh_oprof_la-cell_pyramid18.lo: src/geom/$(am__dirstamp) \
+	src/geom/$(DEPDIR)/$(am__dirstamp)
 src/geom/libmesh_oprof_la-cell_pyramid5.lo: src/geom/$(am__dirstamp) \
 	src/geom/$(DEPDIR)/$(am__dirstamp)
 src/geom/libmesh_oprof_la-cell_tet.lo: src/geom/$(am__dirstamp) \
@@ -11930,6 +11947,8 @@ src/geom/libmesh_opt_la-cell_pyramid13.lo: src/geom/$(am__dirstamp) \
 	src/geom/$(DEPDIR)/$(am__dirstamp)
 src/geom/libmesh_opt_la-cell_pyramid14.lo: src/geom/$(am__dirstamp) \
 	src/geom/$(DEPDIR)/$(am__dirstamp)
+src/geom/libmesh_opt_la-cell_pyramid18.lo: src/geom/$(am__dirstamp) \
+	src/geom/$(DEPDIR)/$(am__dirstamp)
 src/geom/libmesh_opt_la-cell_pyramid5.lo: src/geom/$(am__dirstamp) \
 	src/geom/$(DEPDIR)/$(am__dirstamp)
 src/geom/libmesh_opt_la-cell_tet.lo: src/geom/$(am__dirstamp) \
@@ -13027,6 +13046,8 @@ src/geom/libmesh_prof_la-cell_pyramid13.lo: src/geom/$(am__dirstamp) \
 	src/geom/$(DEPDIR)/$(am__dirstamp)
 src/geom/libmesh_prof_la-cell_pyramid14.lo: src/geom/$(am__dirstamp) \
 	src/geom/$(DEPDIR)/$(am__dirstamp)
+src/geom/libmesh_prof_la-cell_pyramid18.lo: src/geom/$(am__dirstamp) \
+	src/geom/$(DEPDIR)/$(am__dirstamp)
 src/geom/libmesh_prof_la-cell_pyramid5.lo: src/geom/$(am__dirstamp) \
 	src/geom/$(DEPDIR)/$(am__dirstamp)
 src/geom/libmesh_prof_la-cell_tet.lo: src/geom/$(am__dirstamp) \
@@ -14956,6 +14977,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_dbg_la-cell_pyramid.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_dbg_la-cell_pyramid13.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_dbg_la-cell_pyramid14.Plo@am__quote@ # am--include-marker
+@AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_dbg_la-cell_pyramid18.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_dbg_la-cell_pyramid5.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_dbg_la-cell_tet.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_dbg_la-cell_tet10.Plo@am__quote@ # am--include-marker
@@ -15018,6 +15040,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_devel_la-cell_pyramid.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_devel_la-cell_pyramid13.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_devel_la-cell_pyramid14.Plo@am__quote@ # am--include-marker
+@AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_devel_la-cell_pyramid18.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_devel_la-cell_pyramid5.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_devel_la-cell_tet.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_devel_la-cell_tet10.Plo@am__quote@ # am--include-marker
@@ -15080,6 +15103,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_oprof_la-cell_pyramid.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_oprof_la-cell_pyramid13.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_oprof_la-cell_pyramid14.Plo@am__quote@ # am--include-marker
+@AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_oprof_la-cell_pyramid18.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_oprof_la-cell_pyramid5.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_oprof_la-cell_tet.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_oprof_la-cell_tet10.Plo@am__quote@ # am--include-marker
@@ -15142,6 +15166,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_opt_la-cell_pyramid.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_opt_la-cell_pyramid13.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_opt_la-cell_pyramid14.Plo@am__quote@ # am--include-marker
+@AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_opt_la-cell_pyramid18.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_opt_la-cell_pyramid5.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_opt_la-cell_tet.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_opt_la-cell_tet10.Plo@am__quote@ # am--include-marker
@@ -15204,6 +15229,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_prof_la-cell_pyramid.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_prof_la-cell_pyramid13.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_prof_la-cell_pyramid14.Plo@am__quote@ # am--include-marker
+@AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_prof_la-cell_pyramid18.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_prof_la-cell_pyramid5.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_prof_la-cell_tet.Plo@am__quote@ # am--include-marker
 @AMDEP_TRUE@@am__include@ @am__quote@src/geom/$(DEPDIR)/libmesh_prof_la-cell_tet10.Plo@am__quote@ # am--include-marker
@@ -17591,6 +17617,13 @@ src/geom/libmesh_dbg_la-cell_pyramid14.lo: src/geom/cell_pyramid14.C
 @AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
 @am__fastdepCXX_FALSE@	$(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_dbg_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_dbg_la_CXXFLAGS) $(CXXFLAGS) -c -o src/geom/libmesh_dbg_la-cell_pyramid14.lo `test -f 'src/geom/cell_pyramid14.C' || echo '$(srcdir)/'`src/geom/cell_pyramid14.C
 
+src/geom/libmesh_dbg_la-cell_pyramid18.lo: src/geom/cell_pyramid18.C
+@am__fastdepCXX_TRUE@	$(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_dbg_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_dbg_la_CXXFLAGS) $(CXXFLAGS) -MT src/geom/libmesh_dbg_la-cell_pyramid18.lo -MD -MP -MF src/geom/$(DEPDIR)/libmesh_dbg_la-cell_pyramid18.Tpo -c -o src/geom/libmesh_dbg_la-cell_pyramid18.lo `test -f 'src/geom/cell_pyramid18.C' || echo '$(srcdir)/'`src/geom/cell_pyramid18.C
+@am__fastdepCXX_TRUE@	$(AM_V_at)$(am__mv) src/geom/$(DEPDIR)/libmesh_dbg_la-cell_pyramid18.Tpo src/geom/$(DEPDIR)/libmesh_dbg_la-cell_pyramid18.Plo
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	$(AM_V_CXX)source='src/geom/cell_pyramid18.C' object='src/geom/libmesh_dbg_la-cell_pyramid18.lo' libtool=yes @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@	$(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_dbg_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_dbg_la_CXXFLAGS) $(CXXFLAGS) -c -o src/geom/libmesh_dbg_la-cell_pyramid18.lo `test -f 'src/geom/cell_pyramid18.C' || echo '$(srcdir)/'`src/geom/cell_pyramid18.C
+
 src/geom/libmesh_dbg_la-cell_pyramid5.lo: src/geom/cell_pyramid5.C
 @am__fastdepCXX_TRUE@	$(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_dbg_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_dbg_la_CXXFLAGS) $(CXXFLAGS) -MT src/geom/libmesh_dbg_la-cell_pyramid5.lo -MD -MP -MF src/geom/$(DEPDIR)/libmesh_dbg_la-cell_pyramid5.Tpo -c -o src/geom/libmesh_dbg_la-cell_pyramid5.lo `test -f 'src/geom/cell_pyramid5.C' || echo '$(srcdir)/'`src/geom/cell_pyramid5.C
 @am__fastdepCXX_TRUE@	$(AM_V_at)$(am__mv) src/geom/$(DEPDIR)/libmesh_dbg_la-cell_pyramid5.Tpo src/geom/$(DEPDIR)/libmesh_dbg_la-cell_pyramid5.Plo
@@ -20713,6 +20746,13 @@ src/geom/libmesh_devel_la-cell_pyramid14.lo: src/geom/cell_pyramid14.C
 @AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
 @am__fastdepCXX_FALSE@	$(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_devel_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_devel_la_CXXFLAGS) $(CXXFLAGS) -c -o src/geom/libmesh_devel_la-cell_pyramid14.lo `test -f 'src/geom/cell_pyramid14.C' || echo '$(srcdir)/'`src/geom/cell_pyramid14.C
 
+src/geom/libmesh_devel_la-cell_pyramid18.lo: src/geom/cell_pyramid18.C
+@am__fastdepCXX_TRUE@	$(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_devel_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_devel_la_CXXFLAGS) $(CXXFLAGS) -MT src/geom/libmesh_devel_la-cell_pyramid18.lo -MD -MP -MF src/geom/$(DEPDIR)/libmesh_devel_la-cell_pyramid18.Tpo -c -o src/geom/libmesh_devel_la-cell_pyramid18.lo `test -f 'src/geom/cell_pyramid18.C' || echo '$(srcdir)/'`src/geom/cell_pyramid18.C
+@am__fastdepCXX_TRUE@	$(AM_V_at)$(am__mv) src/geom/$(DEPDIR)/libmesh_devel_la-cell_pyramid18.Tpo src/geom/$(DEPDIR)/libmesh_devel_la-cell_pyramid18.Plo
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	$(AM_V_CXX)source='src/geom/cell_pyramid18.C' object='src/geom/libmesh_devel_la-cell_pyramid18.lo' libtool=yes @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@	$(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_devel_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_devel_la_CXXFLAGS) $(CXXFLAGS) -c -o src/geom/libmesh_devel_la-cell_pyramid18.lo `test -f 'src/geom/cell_pyramid18.C' || echo '$(srcdir)/'`src/geom/cell_pyramid18.C
+
 src/geom/libmesh_devel_la-cell_pyramid5.lo: src/geom/cell_pyramid5.C
 @am__fastdepCXX_TRUE@	$(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_devel_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_devel_la_CXXFLAGS) $(CXXFLAGS) -MT src/geom/libmesh_devel_la-cell_pyramid5.lo -MD -MP -MF src/geom/$(DEPDIR)/libmesh_devel_la-cell_pyramid5.Tpo -c -o src/geom/libmesh_devel_la-cell_pyramid5.lo `test -f 'src/geom/cell_pyramid5.C' || echo '$(srcdir)/'`src/geom/cell_pyramid5.C
 @am__fastdepCXX_TRUE@	$(AM_V_at)$(am__mv) src/geom/$(DEPDIR)/libmesh_devel_la-cell_pyramid5.Tpo src/geom/$(DEPDIR)/libmesh_devel_la-cell_pyramid5.Plo
@@ -23835,6 +23875,13 @@ src/geom/libmesh_oprof_la-cell_pyramid14.lo: src/geom/cell_pyramid14.C
 @AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
 @am__fastdepCXX_FALSE@	$(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_oprof_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_oprof_la_CXXFLAGS) $(CXXFLAGS) -c -o src/geom/libmesh_oprof_la-cell_pyramid14.lo `test -f 'src/geom/cell_pyramid14.C' || echo '$(srcdir)/'`src/geom/cell_pyramid14.C
 
+src/geom/libmesh_oprof_la-cell_pyramid18.lo: src/geom/cell_pyramid18.C
+@am__fastdepCXX_TRUE@	$(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_oprof_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_oprof_la_CXXFLAGS) $(CXXFLAGS) -MT src/geom/libmesh_oprof_la-cell_pyramid18.lo -MD -MP -MF src/geom/$(DEPDIR)/libmesh_oprof_la-cell_pyramid18.Tpo -c -o src/geom/libmesh_oprof_la-cell_pyramid18.lo `test -f 'src/geom/cell_pyramid18.C' || echo '$(srcdir)/'`src/geom/cell_pyramid18.C
+@am__fastdepCXX_TRUE@	$(AM_V_at)$(am__mv) src/geom/$(DEPDIR)/libmesh_oprof_la-cell_pyramid18.Tpo src/geom/$(DEPDIR)/libmesh_oprof_la-cell_pyramid18.Plo
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	$(AM_V_CXX)source='src/geom/cell_pyramid18.C' object='src/geom/libmesh_oprof_la-cell_pyramid18.lo' libtool=yes @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@	$(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_oprof_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_oprof_la_CXXFLAGS) $(CXXFLAGS) -c -o src/geom/libmesh_oprof_la-cell_pyramid18.lo `test -f 'src/geom/cell_pyramid18.C' || echo '$(srcdir)/'`src/geom/cell_pyramid18.C
+
 src/geom/libmesh_oprof_la-cell_pyramid5.lo: src/geom/cell_pyramid5.C
 @am__fastdepCXX_TRUE@	$(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_oprof_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_oprof_la_CXXFLAGS) $(CXXFLAGS) -MT src/geom/libmesh_oprof_la-cell_pyramid5.lo -MD -MP -MF src/geom/$(DEPDIR)/libmesh_oprof_la-cell_pyramid5.Tpo -c -o src/geom/libmesh_oprof_la-cell_pyramid5.lo `test -f 'src/geom/cell_pyramid5.C' || echo '$(srcdir)/'`src/geom/cell_pyramid5.C
 @am__fastdepCXX_TRUE@	$(AM_V_at)$(am__mv) src/geom/$(DEPDIR)/libmesh_oprof_la-cell_pyramid5.Tpo src/geom/$(DEPDIR)/libmesh_oprof_la-cell_pyramid5.Plo
@@ -26957,6 +27004,13 @@ src/geom/libmesh_opt_la-cell_pyramid14.lo: src/geom/cell_pyramid14.C
 @AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
 @am__fastdepCXX_FALSE@	$(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_opt_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_opt_la_CXXFLAGS) $(CXXFLAGS) -c -o src/geom/libmesh_opt_la-cell_pyramid14.lo `test -f 'src/geom/cell_pyramid14.C' || echo '$(srcdir)/'`src/geom/cell_pyramid14.C
 
+src/geom/libmesh_opt_la-cell_pyramid18.lo: src/geom/cell_pyramid18.C
+@am__fastdepCXX_TRUE@	$(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_opt_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_opt_la_CXXFLAGS) $(CXXFLAGS) -MT src/geom/libmesh_opt_la-cell_pyramid18.lo -MD -MP -MF src/geom/$(DEPDIR)/libmesh_opt_la-cell_pyramid18.Tpo -c -o src/geom/libmesh_opt_la-cell_pyramid18.lo `test -f 'src/geom/cell_pyramid18.C' || echo '$(srcdir)/'`src/geom/cell_pyramid18.C
+@am__fastdepCXX_TRUE@	$(AM_V_at)$(am__mv) src/geom/$(DEPDIR)/libmesh_opt_la-cell_pyramid18.Tpo src/geom/$(DEPDIR)/libmesh_opt_la-cell_pyramid18.Plo
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	$(AM_V_CXX)source='src/geom/cell_pyramid18.C' object='src/geom/libmesh_opt_la-cell_pyramid18.lo' libtool=yes @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@	$(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_opt_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_opt_la_CXXFLAGS) $(CXXFLAGS) -c -o src/geom/libmesh_opt_la-cell_pyramid18.lo `test -f 'src/geom/cell_pyramid18.C' || echo '$(srcdir)/'`src/geom/cell_pyramid18.C
+
 src/geom/libmesh_opt_la-cell_pyramid5.lo: src/geom/cell_pyramid5.C
 @am__fastdepCXX_TRUE@	$(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_opt_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_opt_la_CXXFLAGS) $(CXXFLAGS) -MT src/geom/libmesh_opt_la-cell_pyramid5.lo -MD -MP -MF src/geom/$(DEPDIR)/libmesh_opt_la-cell_pyramid5.Tpo -c -o src/geom/libmesh_opt_la-cell_pyramid5.lo `test -f 'src/geom/cell_pyramid5.C' || echo '$(srcdir)/'`src/geom/cell_pyramid5.C
 @am__fastdepCXX_TRUE@	$(AM_V_at)$(am__mv) src/geom/$(DEPDIR)/libmesh_opt_la-cell_pyramid5.Tpo src/geom/$(DEPDIR)/libmesh_opt_la-cell_pyramid5.Plo
@@ -30079,6 +30133,13 @@ src/geom/libmesh_prof_la-cell_pyramid14.lo: src/geom/cell_pyramid14.C
 @AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
 @am__fastdepCXX_FALSE@	$(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_prof_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_prof_la_CXXFLAGS) $(CXXFLAGS) -c -o src/geom/libmesh_prof_la-cell_pyramid14.lo `test -f 'src/geom/cell_pyramid14.C' || echo '$(srcdir)/'`src/geom/cell_pyramid14.C
 
+src/geom/libmesh_prof_la-cell_pyramid18.lo: src/geom/cell_pyramid18.C
+@am__fastdepCXX_TRUE@	$(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_prof_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_prof_la_CXXFLAGS) $(CXXFLAGS) -MT src/geom/libmesh_prof_la-cell_pyramid18.lo -MD -MP -MF src/geom/$(DEPDIR)/libmesh_prof_la-cell_pyramid18.Tpo -c -o src/geom/libmesh_prof_la-cell_pyramid18.lo `test -f 'src/geom/cell_pyramid18.C' || echo '$(srcdir)/'`src/geom/cell_pyramid18.C
+@am__fastdepCXX_TRUE@	$(AM_V_at)$(am__mv) src/geom/$(DEPDIR)/libmesh_prof_la-cell_pyramid18.Tpo src/geom/$(DEPDIR)/libmesh_prof_la-cell_pyramid18.Plo
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	$(AM_V_CXX)source='src/geom/cell_pyramid18.C' object='src/geom/libmesh_prof_la-cell_pyramid18.lo' libtool=yes @AMDEPBACKSLASH@
+@AMDEP_TRUE@@am__fastdepCXX_FALSE@	DEPDIR=$(DEPDIR) $(CXXDEPMODE) $(depcomp) @AMDEPBACKSLASH@
+@am__fastdepCXX_FALSE@	$(AM_V_CXX@am__nodep@)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_prof_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_prof_la_CXXFLAGS) $(CXXFLAGS) -c -o src/geom/libmesh_prof_la-cell_pyramid18.lo `test -f 'src/geom/cell_pyramid18.C' || echo '$(srcdir)/'`src/geom/cell_pyramid18.C
+
 src/geom/libmesh_prof_la-cell_pyramid5.lo: src/geom/cell_pyramid5.C
 @am__fastdepCXX_TRUE@	$(AM_V_CXX)$(LIBTOOL) $(AM_V_lt) --tag=CXX $(AM_LIBTOOLFLAGS) $(LIBTOOLFLAGS) --mode=compile $(CXX) $(DEFS) $(DEFAULT_INCLUDES) $(INCLUDES) $(libmesh_prof_la_CPPFLAGS) $(CPPFLAGS) $(libmesh_prof_la_CXXFLAGS) $(CXXFLAGS) -MT src/geom/libmesh_prof_la-cell_pyramid5.lo -MD -MP -MF src/geom/$(DEPDIR)/libmesh_prof_la-cell_pyramid5.Tpo -c -o src/geom/libmesh_prof_la-cell_pyramid5.lo `test -f 'src/geom/cell_pyramid5.C' || echo '$(srcdir)/'`src/geom/cell_pyramid5.C
 @am__fastdepCXX_TRUE@	$(AM_V_at)$(am__mv) src/geom/$(DEPDIR)/libmesh_prof_la-cell_pyramid5.Tpo src/geom/$(DEPDIR)/libmesh_prof_la-cell_pyramid5.Plo
@@ -34106,6 +34167,7 @@ distclean: distclean-recursive
 	-rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_pyramid.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_pyramid13.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_pyramid14.Plo
+	-rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_pyramid18.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_pyramid5.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_tet.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_tet10.Plo
@@ -34168,6 +34230,7 @@ distclean: distclean-recursive
 	-rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_pyramid.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_pyramid13.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_pyramid14.Plo
+	-rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_pyramid18.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_pyramid5.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_tet.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_tet10.Plo
@@ -34230,6 +34293,7 @@ distclean: distclean-recursive
 	-rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_pyramid.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_pyramid13.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_pyramid14.Plo
+	-rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_pyramid18.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_pyramid5.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_tet.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_tet10.Plo
@@ -34292,6 +34356,7 @@ distclean: distclean-recursive
 	-rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_pyramid.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_pyramid13.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_pyramid14.Plo
+	-rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_pyramid18.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_pyramid5.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_tet.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_tet10.Plo
@@ -34354,6 +34419,7 @@ distclean: distclean-recursive
 	-rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_pyramid.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_pyramid13.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_pyramid14.Plo
+	-rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_pyramid18.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_pyramid5.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_tet.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_tet10.Plo
@@ -36441,6 +36507,7 @@ maintainer-clean: maintainer-clean-recursive
 	-rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_pyramid.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_pyramid13.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_pyramid14.Plo
+	-rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_pyramid18.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_pyramid5.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_tet.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_dbg_la-cell_tet10.Plo
@@ -36503,6 +36570,7 @@ maintainer-clean: maintainer-clean-recursive
 	-rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_pyramid.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_pyramid13.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_pyramid14.Plo
+	-rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_pyramid18.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_pyramid5.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_tet.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_devel_la-cell_tet10.Plo
@@ -36565,6 +36633,7 @@ maintainer-clean: maintainer-clean-recursive
 	-rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_pyramid.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_pyramid13.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_pyramid14.Plo
+	-rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_pyramid18.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_pyramid5.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_tet.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_oprof_la-cell_tet10.Plo
@@ -36627,6 +36696,7 @@ maintainer-clean: maintainer-clean-recursive
 	-rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_pyramid.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_pyramid13.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_pyramid14.Plo
+	-rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_pyramid18.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_pyramid5.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_tet.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_opt_la-cell_tet10.Plo
@@ -36689,6 +36759,7 @@ maintainer-clean: maintainer-clean-recursive
 	-rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_pyramid.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_pyramid13.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_pyramid14.Plo
+	-rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_pyramid18.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_pyramid5.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_tet.Plo
 	-rm -f src/geom/$(DEPDIR)/libmesh_prof_la-cell_tet10.Plo
diff --git a/include/Makefile.in b/include/Makefile.in
index 328648c0577..3511610a9eb 100644
--- a/include/Makefile.in
+++ b/include/Makefile.in
@@ -717,6 +717,7 @@ include_HEADERS = \
         geom/cell_pyramid.h \
         geom/cell_pyramid13.h \
         geom/cell_pyramid14.h \
+        geom/cell_pyramid18.h \
         geom/cell_pyramid5.h \
         geom/cell_tet.h \
         geom/cell_tet10.h \
@@ -1074,6 +1075,7 @@ include_HEADERS = \
         utils/tree_node.h \
         utils/utility.h \
         utils/vectormap.h \
+        utils/win_gettimeofday.h \
         utils/xdr_cxx.h 
 
 
diff --git a/include/fe/fe.h b/include/fe/fe.h
index 34d64d7d7c3..8f34d12e566 100644
--- a/include/fe/fe.h
+++ b/include/fe/fe.h
@@ -1324,6 +1324,16 @@ OutputShape fe_fdm_deriv(const Elem * elem,
                             const unsigned int, const Point &,
                             const bool));
 
+template <typename OutputShape>
+OutputShape fe_fdm_deriv(const ElemType type,
+                         const Order order,
+                         const unsigned int i,
+                         const unsigned int j,
+                         const Point & p,
+                         OutputShape(*shape_func)
+                           (const ElemType, const Order,
+                            const unsigned int, const Point &));
+
 template <typename OutputShape>
 OutputShape
 fe_fdm_second_deriv(const Elem * elem,
@@ -1337,6 +1347,18 @@ fe_fdm_second_deriv(const Elem * elem,
                        const unsigned int, const unsigned int,
                        const Point &, const bool));
 
+template <typename OutputShape>
+OutputShape fe_fdm_second_deriv(const ElemType type,
+                                const Order order,
+                                const unsigned int i,
+                                const unsigned int j,
+                                const Point & p,
+                                OutputShape(*deriv_func)
+                                  (const ElemType, const Order,
+                                   const unsigned int,
+                                   const unsigned int,
+                                   const Point &));
+
 /**
  * Helper functions for Lagrange-based basis functions.
  */
diff --git a/include/geom/cell_prism20.h b/include/geom/cell_prism20.h
index ef4b8513800..b9d0d1ab50f 100644
--- a/include/geom/cell_prism20.h
+++ b/include/geom/cell_prism20.h
@@ -211,7 +211,8 @@ class Prism20 final : public Prism
                             std::vector<dof_id_type> & conn) const override;
 
   /**
-   * \returns 2 for all edge nodes and 4 for face nodes.
+   * \returns 2 for all edge nodes, 4 for quad face nodes, 6 for tri
+   * face nodes.
    */
   virtual unsigned int n_second_order_adjacent_vertices (const unsigned int) const override;
 
diff --git a/include/geom/cell_prism21.h b/include/geom/cell_prism21.h
index d31edd49825..55bd80fcaef 100644
--- a/include/geom/cell_prism21.h
+++ b/include/geom/cell_prism21.h
@@ -214,7 +214,8 @@ class Prism21 final : public Prism
                             std::vector<dof_id_type> & conn) const override;
 
   /**
-   * \returns 2 for all edge nodes and 4 for face nodes.
+   * \returns 2 for all edge nodes, 4 for quad face nodes, 3 for tri
+   * face nodes, 6 for the interior node.
    */
   virtual unsigned int n_second_order_adjacent_vertices (const unsigned int) const override;
 
diff --git a/include/geom/cell_pyramid18.h b/include/geom/cell_pyramid18.h
new file mode 100644
index 00000000000..55d293499a0
--- /dev/null
+++ b/include/geom/cell_pyramid18.h
@@ -0,0 +1,287 @@
+// The libMesh Finite Element Library.
+// Copyright (C) 2002-2022 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
+
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+
+
+
+#ifndef LIBMESH_CELL_PYRAMID18_H
+#define LIBMESH_CELL_PYRAMID18_H
+
+// Local includes
+#include "libmesh/cell_pyramid.h"
+
+namespace libMesh
+{
+
+/**
+ * The \p Pyramid18 is an element in 3D composed of 18 nodes, designed
+ * to interface with a QUAD9 element on the base and a TRI7 element on
+ * each of the triangular faces.  Cubit will generate hybrid meshes
+ * with linear pyramids, but as of version 18 will not export
+ * quadratic pyramids.  Paraview may support 13-node pyramids, but
+ * does not render 18-node pyramids correctly.  So even if this
+ * element works in libmesh, we are currently limited in what we can do
+ * with it outside the library...
+ *
+ * The node numbering for the pyramid18 is given below:
+ *
+ * \verbatim
+ *   PYRAMID18:
+ *                       o 4
+ *                     //|\
+ *                    // | \
+ *                   //  |  \
+ *                  //   |   \
+ *              12 o/    |    o 11
+ *                //     |     \
+ *               /o 9    o 10   \
+ *              //       |       \          zeta
+ *             //        |        \          ^   eta (into page)
+ *          3 o/.......o.|........o 2        | /
+ *           ./       7  |       /           |/
+ *          ./           |      /            o---> xi
+ *         ./            |     /
+ *        ./             |    /
+ *     8 o/       o      |   o 6
+ *      ./        13     |  /
+ *     ./                | /
+ *    ./                 |/
+ *    o--------o---------o
+ *    0        5         1
+ * \endverbatim
+ *
+ * And it also includes four triangle face nodes:
+ * Node 14, on side 0, equidistant from 0/1/4 or 5/9/10
+ * Node 15, on side 1, equidistant from 1/2/4 or 6/10/11
+ * Node 16, on side 2, equidistant from 2/3/4 or 7/11/12
+ * Node 17, on side 3, equidistant from 0/3/4 or 8/9/12
+ *
+ * (xi, eta, zeta): { zeta-1 <= xi   <= 1-zeta
+ *                  { zeta-1 <= eta  <= 1-zeta
+ *                  {      0 <= zeta <= 1
+ * are the reference element coordinates associated with the given
+ * numbering.
+ *
+ * \author Roy H. Stogner
+ * \date 2022
+ * \brief A 3D pyramid element with 18 nodes.
+ */
+class Pyramid18 final : public Pyramid
+{
+public:
+
+  /**
+   * Constructor.  By default this element has no parent.
+   */
+  explicit
+  Pyramid18 (Elem * p=nullptr) :
+    Pyramid(Pyramid18::n_nodes(), p, _nodelinks_data)
+  {}
+
+  Pyramid18 (Pyramid18 &&) = delete;
+  Pyramid18 (const Pyramid18 &) = delete;
+  Pyramid18 & operator= (const Pyramid18 &) = delete;
+  Pyramid18 & operator= (Pyramid18 &&) = delete;
+  virtual ~Pyramid18() = default;
+
+  /**
+   * \returns 18.
+   */
+  virtual unsigned int n_nodes() const override { return num_nodes; }
+
+  /**
+   * \returns \p PYRAMID18.
+   */
+  virtual ElemType type () const override { return PYRAMID18; }
+
+  /**
+   * FIXME: we don't yet have a refinement pattern for pyramids...
+   * \returns 1.
+   */
+  virtual unsigned int n_sub_elem() const override { return 1; }
+
+  /**
+   * \returns \p true if the specified (local) node number is a vertex.
+   */
+  virtual bool is_vertex(const unsigned int i) const override;
+
+  /**
+   * \returns \p true if the specified (local) node number is an edge.
+   */
+  virtual bool is_edge(const unsigned int i) const override;
+
+  /**
+   * \returns \p true if the specified (local) node number is a face.
+   */
+  virtual bool is_face(const unsigned int i) const override;
+
+  /**
+   * \returns \p true if the specified (local) node number is on the
+   * specified side.
+   */
+  virtual bool is_node_on_side(const unsigned int n,
+                               const unsigned int s) const override;
+
+  virtual std::vector<unsigned int> nodes_on_side(const unsigned int s) const override;
+
+  virtual std::vector<unsigned int> nodes_on_edge(const unsigned int e) const override;
+
+  /**
+   * \returns \p true if the specified (local) node number is on the
+   * specified edge.
+   */
+  virtual bool is_node_on_edge(const unsigned int n,
+                               const unsigned int e) const override;
+
+  /**
+   * \returns \p true if the element map is definitely affine within
+   * numerical tolerances.
+   */
+  virtual bool has_affine_map () const override;
+
+  /**
+   * \returns SECOND.
+   */
+  virtual Order default_order() const override;
+
+  /**
+   * Don't hide Pyramid::key() defined in the base class.
+   */
+  using Pyramid::key;
+
+  /**
+   * \returns An id associated with the \p s side of this element.
+   * The id is not necessarily unique, but should be close.  This is
+   * particularly useful in the \p MeshBase::find_neighbors() routine.
+   *
+   * We reimplement this method here for the \p Pyramid18 since we can
+   * use the center node of the base face to provide a perfect (unique)
+   * key.
+   */
+  virtual dof_id_type key (const unsigned int s) const override;
+
+  /**
+   * \returns \p Pyramid18::side_nodes_map[side][side_node] after doing some range checking.
+   */
+  virtual unsigned int local_side_node(unsigned int side,
+                                       unsigned int side_node) const override;
+
+  /**
+   * \returns \p Pyramid18::edge_nodes_map[edge][edge_node] after doing some range checking.
+   */
+  virtual unsigned int local_edge_node(unsigned int edge,
+                                       unsigned int edge_node) const override;
+
+  /**
+   * Builds a \p QUAD9 or \p TRI7 coincident with face i.
+   * The \p std::unique_ptr<Elem> handles the memory aspect.
+   */
+  virtual std::unique_ptr<Elem> build_side_ptr (const unsigned int i,
+                                                bool proxy=false) override;
+
+  /**
+   * Rebuilds a \p QUAD9 or \p TRI7 built coincident with face i.
+   */
+  virtual void build_side_ptr (std::unique_ptr<Elem> & elem,
+                               const unsigned int i) override;
+
+  /**
+   * Builds a \p EDGE3 coincident with edge i.
+   * The \p std::unique_ptr<Elem> handles the memory aspect.
+   */
+  virtual std::unique_ptr<Elem> build_edge_ptr (const unsigned int i) override;
+
+  /**
+   * Rebuilds a \p EDGE3 coincident with edge i.
+   */
+  virtual void build_edge_ptr (std::unique_ptr<Elem> & edge, const unsigned int i) override;
+
+  virtual void connectivity(const unsigned int sc,
+                            const IOPackage iop,
+                            std::vector<dof_id_type> & conn) const override;
+
+  /**
+   * \returns 2 for all edge nodes, 4 for quad face nodes, 3 for tri
+   * face nodes.
+   */
+  virtual unsigned int n_second_order_adjacent_vertices (const unsigned int n) const override;
+
+  /**
+   * \returns The element-local number of the \f$ v^{th} \f$ vertex
+   * that defines the \f$ n^{th} \f$ second-order node.
+   */
+  virtual unsigned short int second_order_adjacent_vertex (const unsigned int n,
+                                                           const unsigned int v) const override;
+
+  /**
+   * Geometric constants for Pyramid18.
+   */
+  static const int num_nodes = 18;
+  static const int num_sides = 5;
+  static const int num_edges = 8;
+  static const int num_children = 0; // not implemented
+  static const int nodes_per_side = 9;
+  static const int nodes_per_edge = 3;
+
+  /**
+   * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ side to
+   * element node numbers.
+   */
+  static const unsigned int side_nodes_map[num_sides][nodes_per_side];
+
+  /**
+   * This maps the \f$ j^{th} \f$ node of the \f$ i^{th} \f$ edge to
+   * element node numbers.
+   */
+  static const unsigned int edge_nodes_map[num_edges][nodes_per_edge];
+
+  virtual void permute(unsigned int perm_num) override final;
+
+  virtual void flip(BoundaryInfo *) override final;
+
+  unsigned int center_node_on_side(const unsigned short side) const override final;
+
+  ElemType side_type (const unsigned int s) const override final;
+
+protected:
+
+  /**
+   * Data for links to nodes.
+   */
+  Node * _nodelinks_data[num_nodes];
+
+
+
+#ifdef LIBMESH_ENABLE_AMR
+
+  /**
+   * Matrix used to create the elements children.
+   */
+  virtual Real embedding_matrix (const unsigned int,
+                                 const unsigned int,
+                                 const unsigned int) const override
+  { libmesh_not_implemented(); return 0.; }
+
+  LIBMESH_ENABLE_TOPOLOGY_CACHES;
+
+#endif // LIBMESH_ENABLE_AMR
+
+};
+
+} // namespace libMesh
+
+
+#endif // LIBMESH_CELL_PYRAMID18_H
diff --git a/include/include_HEADERS b/include/include_HEADERS
index 85b43d9d4f1..ac26f86e084 100644
--- a/include/include_HEADERS
+++ b/include/include_HEADERS
@@ -127,6 +127,7 @@ include_HEADERS =  \
         geom/cell_pyramid.h \
         geom/cell_pyramid13.h \
         geom/cell_pyramid14.h \
+        geom/cell_pyramid18.h \
         geom/cell_pyramid5.h \
         geom/cell_tet.h \
         geom/cell_tet10.h \
@@ -484,4 +485,5 @@ include_HEADERS =  \
         utils/tree_node.h \
         utils/utility.h \
         utils/vectormap.h \
+        utils/win_gettimeofday.h \
         utils/xdr_cxx.h 
diff --git a/include/libmesh/Makefile.am b/include/libmesh/Makefile.am
index 2d49e348066..66aa7993112 100644
--- a/include/libmesh/Makefile.am
+++ b/include/libmesh/Makefile.am
@@ -117,6 +117,7 @@ BUILT_SOURCES = \
         cell_pyramid.h \
         cell_pyramid13.h \
         cell_pyramid14.h \
+        cell_pyramid18.h \
         cell_pyramid5.h \
         cell_tet.h \
         cell_tet10.h \
@@ -911,6 +912,9 @@ cell_pyramid13.h: $(top_srcdir)/include/geom/cell_pyramid13.h
 cell_pyramid14.h: $(top_srcdir)/include/geom/cell_pyramid14.h
 	$(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@
 
+cell_pyramid18.h: $(top_srcdir)/include/geom/cell_pyramid18.h
+	$(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@
+
 cell_pyramid5.h: $(top_srcdir)/include/geom/cell_pyramid5.h
 	$(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@
 
diff --git a/include/libmesh/Makefile.in b/include/libmesh/Makefile.in
index b1b04d126ea..153bff28ac0 100644
--- a/include/libmesh/Makefile.in
+++ b/include/libmesh/Makefile.in
@@ -544,9 +544,9 @@ BUILT_SOURCES = auto_ptr.h default_coupling.h dirichlet_boundaries.h \
 	cell_inf_prism12.h cell_inf_prism6.h cell_prism.h \
 	cell_prism15.h cell_prism18.h cell_prism20.h cell_prism21.h \
 	cell_prism6.h cell_pyramid.h cell_pyramid13.h cell_pyramid14.h \
-	cell_pyramid5.h cell_tet.h cell_tet10.h cell_tet14.h \
-	cell_tet4.h compare_elems_by_level.h edge.h edge_edge2.h \
-	edge_edge3.h edge_edge4.h edge_inf_edge2.h elem.h \
+	cell_pyramid18.h cell_pyramid5.h cell_tet.h cell_tet10.h \
+	cell_tet14.h cell_tet4.h compare_elems_by_level.h edge.h \
+	edge_edge2.h edge_edge3.h edge_edge4.h edge_inf_edge2.h elem.h \
 	elem_cutter.h elem_hash.h elem_internal.h elem_quality.h \
 	elem_range.h elem_side_builder.h face.h face_inf_quad.h \
 	face_inf_quad4.h face_inf_quad6.h face_quad.h face_quad4.h \
@@ -1246,6 +1246,9 @@ cell_pyramid13.h: $(top_srcdir)/include/geom/cell_pyramid13.h
 cell_pyramid14.h: $(top_srcdir)/include/geom/cell_pyramid14.h
 	$(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@
 
+cell_pyramid18.h: $(top_srcdir)/include/geom/cell_pyramid18.h
+	$(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@
+
 cell_pyramid5.h: $(top_srcdir)/include/geom/cell_pyramid5.h
 	$(AM_V_GEN)rm -f $@ && $(LN_S) -f $< $@
 
diff --git a/include/utils/win_gettimeofday.h b/include/utils/win_gettimeofday.h
index d0a9de2cec9..450daa534c8 100644
--- a/include/utils/win_gettimeofday.h
+++ b/include/utils/win_gettimeofday.h
@@ -5,6 +5,9 @@
 #ifndef _WIN_GETTIMEOFDAY_H_
 #define _WIN_GETTIMEOFDAY_H_
 
+#include "libmesh/libmesh_config.h"
+
+#ifndef LIBMESH_HAVE_GETTIMEOFDAY
 #ifdef _MSC_VER
 #define NOMINMAX
 #define WIN32_LEAN_AND_MEAN
@@ -39,5 +42,6 @@ inline int gettimeofday(struct timeval * tp, struct timezone * tzp)
 #else // _MSC_VER
 #error "gettimeofday() is not implemented"
 #endif // _MSC_VER
+#endif // !LIBMESH_HAVE_GETTIMEOFDAY
 
 #endif // _WIN_GETTIMEOFDAY_H_
diff --git a/reference_elements/3D/one_pyramid18.xda b/reference_elements/3D/one_pyramid18.xda
new file mode 100644
index 00000000000..4e8b0c97381
--- /dev/null
+++ b/reference_elements/3D/one_pyramid18.xda
@@ -0,0 +1,33 @@
+libMesh-0.7.0+
+1	 # number of elements
+18	 # number of nodes
+n/a	 # boundary condition specification file
+.	 # subdomain id specification file
+n/a	 # processor id specification file
+n/a	 # p-level specification file
+1	 # n_elem at level 0, [ type sid (n0 ... nN-1) ]
+37 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
+-1.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00
+1.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00
+1.0000000000000000e+00 1.0000000000000000e+00 0.0000000000000000e+00
+-1.0000000000000000e+00 1.0000000000000000e+00 0.0000000000000000e+00
+0.0000000000000000e+00 0.0000000000000000e+00 1.0000000000000000e+00
+0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00
+1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+0.0000000000000000e+00 1.0000000000000000e+00 0.0000000000000000e+00
+-1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+-5.0000000000000000e-01 -5.0000000000000000e-01 5.0000000000000000e-01
+5.0000000000000000e-01 -5.0000000000000000e-01 5.0000000000000000e-01
+5.0000000000000000e-01 5.0000000000000000e-01 5.0000000000000000e-01
+-5.0000000000000000e-01 5.0000000000000000e-01 5.0000000000000000e-01
+0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00
+0.0000000000000000e+00 -6.6666666666666666e-01 3.3333333333333333e-01
+6.6666666666666666e-01 0.0000000000000000e+00 3.3333333333333333e-01
+0.0000000000000000e+00 6.6666666666666666e-01 3.3333333333333333e-01
+-6.6666666666666666e-01 0.0000000000000000e+00 3.3333333333333333e-01
+5	 # number of boundary conditions
+0 0 0 
+0 1 1 
+0 2 2 
+0 3 3 
+0 4 4 
diff --git a/src/fe/fe.C b/src/fe/fe.C
index 966eb6043cb..5234c371ac6 100644
--- a/src/fe/fe.C
+++ b/src/fe/fe.C
@@ -780,6 +780,57 @@ OutputShape fe_fdm_deriv(const Elem * elem,
 }
 
 
+template <typename OutputShape>
+OutputShape fe_fdm_deriv(const ElemType type,
+                         const Order order,
+                         const unsigned int i,
+                         const unsigned int j,
+                         const Point & p,
+                         OutputShape(*shape_func)
+                           (const ElemType, const Order,
+                            const unsigned int, const Point &))
+{
+  libmesh_assert_less (j, LIBMESH_DIM);
+
+  // cheat by using finite difference approximations:
+  const Real eps = 1.e-6;
+  Point pp = p, pm = p;
+
+  switch (j)
+    {
+      // d()/dxi
+    case 0:
+      {
+        pp(0) += eps;
+        pm(0) -= eps;
+        break;
+      }
+
+      // d()/deta
+    case 1:
+      {
+        pp(1) += eps;
+        pm(1) -= eps;
+        break;
+      }
+
+      // d()/dzeta
+    case 2:
+      {
+        pp(2) += eps;
+        pm(2) -= eps;
+        break;
+      }
+
+    default:
+      libmesh_error_msg("Invalid derivative index j = " << j);
+    }
+
+  return (shape_func(type, order, i, pp) -
+          shape_func(type, order, i, pm))/2./eps;
+}
+
+
 template <typename OutputShape>
 OutputShape
 fe_fdm_second_deriv(const Elem * elem,
@@ -863,6 +914,88 @@ fe_fdm_second_deriv(const Elem * elem,
 }
 
 
+template <typename OutputShape>
+OutputShape
+fe_fdm_second_deriv(const ElemType type,
+                    const Order order,
+                    const unsigned int i,
+                    const unsigned int j,
+                    const Point & p,
+                    OutputShape(*deriv_func)
+                      (const ElemType, const Order,
+                       const unsigned int, const unsigned int,
+                       const Point &))
+{
+  // cheat by using finite difference approximations:
+  const Real eps = 1.e-5;
+  Point pp = p, pm = p;
+  unsigned int deriv_j = 0;
+
+  switch (j)
+    {
+      //  d^2() / dxi^2
+    case 0:
+      {
+        pp(0) += eps;
+        pm(0) -= eps;
+        deriv_j = 0;
+        break;
+      }
+
+      // d^2() / dxi deta
+    case 1:
+      {
+        pp(1) += eps;
+        pm(1) -= eps;
+        deriv_j = 0;
+        break;
+      }
+
+      // d^2() / deta^2
+    case 2:
+      {
+        pp(1) += eps;
+        pm(1) -= eps;
+        deriv_j = 1;
+        break;
+      }
+
+      // d^2()/dxidzeta
+    case 3:
+      {
+        pp(2) += eps;
+        pm(2) -= eps;
+        deriv_j = 0;
+        break;
+      }                  // d^2()/deta^2
+
+      // d^2()/detadzeta
+    case 4:
+      {
+        pp(2) += eps;
+        pm(2) -= eps;
+        deriv_j = 1;
+        break;
+      }
+
+      // d^2()/dzeta^2
+    case 5:
+      {
+        pp(2) += eps;
+        pm(2) -= eps;
+        deriv_j = 2;
+        break;
+      }
+
+    default:
+      libmesh_error_msg("Invalid shape function derivative j = " << j);
+    }
+
+  return (deriv_func(type, order, i, deriv_j, pp) -
+          deriv_func(type, order, i, deriv_j, pm))/2./eps;
+}
+
+
 void rational_fe_weighted_shapes(const Elem * elem,
                                  const FEType underlying_fe_type,
                                  std::vector<std::vector<Real>> & shapes,
@@ -1228,6 +1361,13 @@ Real fe_fdm_deriv<Real>(const Elem *, const Order, const unsigned int,
                           (const Elem *, const Order, const unsigned int,
                            const Point &, const bool));
 
+template
+Real fe_fdm_deriv<Real>(const ElemType, const Order, const unsigned int,
+                        const unsigned int, const Point &,
+                        Real(*shape_func)
+                          (const ElemType, const Order, const unsigned int,
+                           const Point &));
+
 template
 RealGradient
 fe_fdm_deriv<RealGradient>(const Elem *, const Order, const unsigned int,
@@ -1236,6 +1376,14 @@ fe_fdm_deriv<RealGradient>(const Elem *, const Order, const unsigned int,
                              (const Elem *, const Order, const unsigned int,
                               const Point &, const bool));
 
+template
+Real
+fe_fdm_second_deriv<Real>(const ElemType, const Order, const unsigned int,
+                          const unsigned int, const Point &,
+                          Real(*shape_func)
+                            (const ElemType, const Order, const unsigned int,
+                             const unsigned int, const Point &));
+
 template
 Real
 fe_fdm_second_deriv<Real>(const Elem *, const Order, const unsigned int,
diff --git a/src/fe/fe_abstract.C b/src/fe/fe_abstract.C
index cf35923c34a..52388e60b43 100644
--- a/src/fe/fe_abstract.C
+++ b/src/fe/fe_abstract.C
@@ -498,6 +498,16 @@ void FEAbstract::get_refspace_nodes(const ElemType itemType, std::vector<Point>
         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
@@ -675,6 +685,7 @@ bool FEAbstract::on_reference_element(const Point & p, const ElemType t, const R
     case PYRAMID5:
     case PYRAMID13:
     case PYRAMID14:
+    case PYRAMID18:
       {
         // Check that the point is on the same side of all the faces
         // by testing whether:
diff --git a/src/fe/fe_interface.C b/src/fe/fe_interface.C
index 0c9f7d0fa7d..282fef14687 100644
--- a/src/fe/fe_interface.C
+++ b/src/fe/fe_interface.C
@@ -2363,6 +2363,8 @@ unsigned int FEInterface::max_order(const FEType & fe_t,
         case PYRAMID13:
         case PYRAMID14:
           return 2;
+        case PYRAMID18:
+          return 3;
         default:
           return unknown;
         }
@@ -2399,6 +2401,7 @@ unsigned int FEInterface::max_order(const FEType & fe_t,
         case PYRAMID5:
         case PYRAMID13:
         case PYRAMID14:
+        case PYRAMID18:
           return unlimited;
         default:
           return unknown;
@@ -2445,6 +2448,7 @@ unsigned int FEInterface::max_order(const FEType & fe_t,
         case PYRAMID5:
         case PYRAMID13:
         case PYRAMID14:
+        case PYRAMID18:
           return 0;
         default:
           return unknown;
@@ -2485,6 +2489,7 @@ unsigned int FEInterface::max_order(const FEType & fe_t,
         case PYRAMID5:
         case PYRAMID13:
         case PYRAMID14:
+        case PYRAMID18:
           return 0;
         default:
           return unknown;
@@ -2520,6 +2525,7 @@ unsigned int FEInterface::max_order(const FEType & fe_t,
         case PYRAMID5:
         case PYRAMID13:
         case PYRAMID14:
+        case PYRAMID18:
           return unlimited;
         default:
           return unknown;
@@ -2557,6 +2563,7 @@ unsigned int FEInterface::max_order(const FEType & fe_t,
         case PYRAMID5:
         case PYRAMID13:
         case PYRAMID14:
+        case PYRAMID18:
           return 0;
         default:
           return unknown;
@@ -2598,6 +2605,7 @@ unsigned int FEInterface::max_order(const FEType & fe_t,
         case PYRAMID5:
         case PYRAMID13:
         case PYRAMID14:
+        case PYRAMID18:
           return 0;
         default:
           return unknown;
@@ -2642,6 +2650,7 @@ unsigned int FEInterface::max_order(const FEType & fe_t,
         case PYRAMID5:
         case PYRAMID13:
         case PYRAMID14:
+        case PYRAMID18:
           return 0;
         default:
           return unknown;
@@ -2687,6 +2696,7 @@ unsigned int FEInterface::max_order(const FEType & fe_t,
         case PYRAMID5:
         case PYRAMID13:
         case PYRAMID14:
+        case PYRAMID18:
           return 0;
         default:
           return unknown;
diff --git a/src/fe/fe_l2_lagrange.C b/src/fe/fe_l2_lagrange.C
index aa6bf535c03..1a81df4d41c 100644
--- a/src/fe/fe_l2_lagrange.C
+++ b/src/fe/fe_l2_lagrange.C
@@ -85,6 +85,7 @@ unsigned int l2_lagrange_n_dofs(const ElemType t, const Order o)
           case PYRAMID5:
           case PYRAMID13:
           case PYRAMID14:
+          case PYRAMID18:
             return 5;
 
           case INVALID_ELEM:
@@ -147,6 +148,7 @@ unsigned int l2_lagrange_n_dofs(const ElemType t, const Order o)
 
           case PYRAMID5:
           case PYRAMID14:
+          case PYRAMID18:
             return 14;
 
           case INVALID_ELEM:
@@ -181,6 +183,9 @@ unsigned int l2_lagrange_n_dofs(const ElemType t, const Order o)
           case PRISM21:
             return 21;
 
+          case PYRAMID18:
+            return 18;
+
           case INVALID_ELEM:
             return 0;
 
diff --git a/src/fe/fe_lagrange.C b/src/fe/fe_lagrange.C
index daf41d37de1..49aa3176dee 100644
--- a/src/fe/fe_lagrange.C
+++ b/src/fe/fe_lagrange.C
@@ -251,32 +251,34 @@ void lagrange_nodal_soln(const Elem * elem,
               return;
             }
 
-          case PYRAMID13:
+          case PYRAMID18:
             {
-              libmesh_assert_equal_to (elem_soln.size(), 5);
-              libmesh_assert_equal_to (nodal_soln.size(), 13);
+              libmesh_assert_equal_to (nodal_soln.size(), 18);
 
-              nodal_soln[0]  = elem_soln[0];
-              nodal_soln[1]  = elem_soln[1];
-              nodal_soln[2]  = elem_soln[2];
-              nodal_soln[3]  = elem_soln[3];
-              nodal_soln[4]  = elem_soln[4];
-              nodal_soln[5]  = .5*(elem_soln[0] + elem_soln[1]);
-              nodal_soln[6]  = .5*(elem_soln[1] + elem_soln[2]);
-              nodal_soln[7]  = .5*(elem_soln[2] + elem_soln[3]);
-              nodal_soln[8]  = .5*(elem_soln[3] + elem_soln[0]);
-              nodal_soln[9]  = .5*(elem_soln[0] + elem_soln[4]);
-              nodal_soln[10] = .5*(elem_soln[1] + elem_soln[4]);
-              nodal_soln[11] = .5*(elem_soln[2] + elem_soln[4]);
-              nodal_soln[12] = .5*(elem_soln[3] + elem_soln[4]);
+              nodal_soln[14] = (elem_soln[0] + elem_soln[1] + elem_soln[4])/Real(3);
+              nodal_soln[15] = (elem_soln[1] + elem_soln[2] + elem_soln[4])/Real(3);
+              nodal_soln[16] = (elem_soln[2] + elem_soln[3] + elem_soln[4])/Real(3);
+              nodal_soln[17] = (elem_soln[0] + elem_soln[3] + elem_soln[4])/Real(3);
 
-              return;
+              libmesh_fallthrough();
             }
 
           case PYRAMID14:
+            {
+              if (type == PYRAMID14)
+                libmesh_assert_equal_to (nodal_soln.size(), 14);
+
+              nodal_soln[13] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[2] + elem_soln[3]);
+
+              libmesh_fallthrough();
+            }
+
+          case PYRAMID13:
             {
               libmesh_assert_equal_to (elem_soln.size(), 5);
-              libmesh_assert_equal_to (nodal_soln.size(), 14);
+
+              if (type == PYRAMID13)
+                libmesh_assert_equal_to (nodal_soln.size(), 13);
 
               nodal_soln[0]  = elem_soln[0];
               nodal_soln[1]  = elem_soln[1];
@@ -291,11 +293,9 @@ void lagrange_nodal_soln(const Elem * elem,
               nodal_soln[10] = .5*(elem_soln[1] + elem_soln[4]);
               nodal_soln[11] = .5*(elem_soln[2] + elem_soln[4]);
               nodal_soln[12] = .5*(elem_soln[3] + elem_soln[4]);
-              nodal_soln[13] = .25*(elem_soln[0] + elem_soln[1] + elem_soln[2] + elem_soln[3]);
 
               return;
             }
-
           default:
             {
               // By default the element solution _is_ nodal,
@@ -448,6 +448,7 @@ unsigned int lagrange_n_dofs(const ElemType t, const Order o)
           case PYRAMID5:
           case PYRAMID13:
           case PYRAMID14:
+          case PYRAMID18:
             return 5;
 
           case INVALID_ELEM:
@@ -503,6 +504,7 @@ unsigned int lagrange_n_dofs(const ElemType t, const Order o)
             return 13;
 
           case PYRAMID14:
+          case PYRAMID18:
             return 14;
 
           case INVALID_ELEM:
@@ -529,6 +531,9 @@ unsigned int lagrange_n_dofs(const ElemType t, const Order o)
           case PRISM21:
             return 21;
 
+          case PYRAMID18:
+            return 18;
+
           case TRI7:
             return 7;
 
@@ -681,6 +686,7 @@ unsigned int lagrange_n_dofs_at_node(const ElemType t,
           case PYRAMID5:
           case PYRAMID13:
           case PYRAMID14:
+          case PYRAMID18:
             {
               switch (n)
                 {
@@ -729,6 +735,9 @@ unsigned int lagrange_n_dofs_at_node(const ElemType t,
           case PRISM21:
             return (n < 18);
 
+          case PYRAMID18:
+            return (n < 14);
+
           case TRI7:
             return (n < 6);
 
@@ -751,6 +760,7 @@ unsigned int lagrange_n_dofs_at_node(const ElemType t,
           case EDGE4:
           case PRISM20:
           case PRISM21:
+          case PYRAMID18:
           case TRI7:
           case TET14:
             return 1;
diff --git a/src/fe/fe_lagrange_shape_3D.C b/src/fe/fe_lagrange_shape_3D.C
index c19d976efb6..deea82da3a0 100644
--- a/src/fe/fe_lagrange_shape_3D.C
+++ b/src/fe/fe_lagrange_shape_3D.C
@@ -412,6 +412,7 @@ Real fe_lagrange_3D_shape(const ElemType type,
           case PYRAMID5:
           case PYRAMID13:
           case PYRAMID14:
+          case PYRAMID18:
             {
               libmesh_assert_less (i, 5);
 
@@ -775,6 +776,7 @@ Real fe_lagrange_3D_shape(const ElemType type,
                                "High order on first order elements only supported for L2 families");
             libmesh_fallthrough();
           case PYRAMID14:
+          case PYRAMID18:
             {
               libmesh_assert_less (i, 14);
 
@@ -906,6 +908,108 @@ Real fe_lagrange_3D_shape(const ElemType type,
                       fe_lagrange_1D_quadratic_shape(i0[i], p1d));
             }
 
+            // Weird rational shape functions with weirder bubbles...
+          case PYRAMID18:
+            {
+              libmesh_assert_less (i, 18);
+
+              const Real xi   = p(0);
+              const Real eta  = p(1);
+              const Real zeta = p(2);
+              const Real eps  = 1.e-35;
+
+              // The "normalized coordinates" defined by Graglia.  These are
+              // the planes which define the faces of the pyramid.
+              const Real
+                p1 = 0.5*(1. - eta - zeta), // back
+                p2 = 0.5*(1. + xi  - zeta), // left
+                p3 = 0.5*(1. + eta - zeta), // front
+                p4 = 0.5*(1. - xi  - zeta); // right
+
+              // Denominators are perturbed by epsilon to avoid
+              // divide-by-zero issues.
+              const Real
+                den = (-1. + zeta + eps),
+                den2 = den*den;
+
+              // Bubble functions on triangular sides.  We actually
+              // have a degree of freedom to play with here, and I'm
+              // not certain how best to use it, so let's leave it
+              // as a variable in case we figure that out later.
+              constexpr Real alpha = 0.5;
+              const Real
+                bub_f1 = ((1-alpha)*(1-zeta) + alpha*(-eta)),
+                bub_f2 = ((1-alpha)*(1-zeta) + alpha*(xi)),
+                bub_f3 = ((1-alpha)*(1-zeta) + alpha*(eta)),
+                bub_f4 = ((1-alpha)*(1-zeta) + alpha*(-xi));
+
+              const Real
+                bub1 = bub_f1*p1*p2*p4*zeta/den2,
+                bub2 = bub_f2*p1*p2*p3*zeta/den2,
+                bub3 = bub_f3*p2*p3*p4*zeta/den2,
+                bub4 = bub_f4*p1*p3*p4*zeta/den2;
+
+              switch(i)
+                {
+                case 0:
+                  return p4*p1*(xi*eta - zeta + zeta*zeta)/den2 + 3*(bub1+bub4);
+
+                case 1:
+                  return -p1*p2*(xi*eta + zeta - zeta*zeta)/den2 + 3*(bub1+bub2);
+
+                case 2:
+                  return p2*p3*(xi*eta - zeta + zeta*zeta)/den2 + 3*(bub2+bub3);
+
+                case 3:
+                  return -p3*p4*(xi*eta + zeta - zeta*zeta)/den2 + 3*(bub3+bub4);
+
+                case 4:
+                  return zeta*(2.*zeta - 1.) + 3*(bub1+bub2+bub3+bub4);
+
+                case 5:
+                  return -4.*p2*p1*p4*eta/den2 - 12*bub1;
+
+                case 6:
+                  return 4.*p1*p2*p3*xi/den2 - 12*bub2;
+
+                case 7:
+                  return 4.*p2*p3*p4*eta/den2 - 12*bub3;
+
+                case 8:
+                  return -4.*p3*p4*p1*xi/den2 - 12*bub4;
+
+                case 9:
+                  return -4.*p1*p4*zeta/den - 12*(bub1+bub4);
+
+                case 10:
+                  return -4.*p2*p1*zeta/den - 12*(bub1+bub2);
+
+                case 11:
+                  return -4.*p3*p2*zeta/den - 12*(bub2+bub3);
+
+                case 12:
+                  return -4.*p4*p3*zeta/den - 12*(bub3+bub4);
+
+                case 13:
+                  return 16.*p1*p2*p3*p4/den2;
+
+                case 14:
+                  return 27*bub1;
+
+                case 15:
+                  return 27*bub2;
+
+                case 16:
+                  return 27*bub3;
+
+                case 17:
+                  return 27*bub4;
+
+                default:
+                  libmesh_error_msg("Invalid i = " << i);
+                }
+            }
+
             // quadratic Lagrange shape functions with cubic bubbles
           case TET14:
             {
@@ -1187,6 +1291,7 @@ Real fe_lagrange_3D_shape_deriv(const ElemType type,
           case PYRAMID5:
           case PYRAMID13:
           case PYRAMID14:
+          case PYRAMID18:
             {
               libmesh_assert_less (i, 5);
 
@@ -2111,6 +2216,7 @@ Real fe_lagrange_3D_shape_deriv(const ElemType type,
                                "High order on first order elements only supported for L2 families");
             libmesh_fallthrough();
           case PYRAMID14:
+          case PYRAMID18:
             {
               libmesh_assert_less (i, 14);
 
@@ -2630,6 +2736,12 @@ Real fe_lagrange_3D_shape_deriv(const ElemType type,
                   libmesh_error_msg("Invalid shape function derivative j = " << j);
                 }
             }
+
+          case PYRAMID18:
+            {
+              return fe_fdm_deriv(type, order, i, j, p, fe_lagrange_3D_shape<T>);
+            }
+
           default:
             libmesh_error_msg("ERROR: Unsupported 3D element type!: " << Utility::enum_to_string(type));
           }
@@ -2724,6 +2836,7 @@ Real fe_lagrange_3D_shape_second_deriv(const ElemType type,
           case PYRAMID5:
           case PYRAMID13:
           case PYRAMID14:
+          case PYRAMID18:
             {
               libmesh_assert_less (i, 5);
 
@@ -3568,6 +3681,7 @@ Real fe_lagrange_3D_shape_second_deriv(const ElemType type,
                                "High order on first order elements only supported for L2 families");
             libmesh_fallthrough();
           case PYRAMID14:
+          case PYRAMID18:
             {
               libmesh_assert_less (i, 14);
 
@@ -4692,6 +4806,13 @@ Real fe_lagrange_3D_shape_second_deriv(const ElemType type,
                   libmesh_error_msg("Invalid shape function derivative j = " << j);
                 }
             }
+
+          case PYRAMID18:
+            {
+              return fe_fdm_second_deriv(type, order, i, j, p,
+                                         fe_lagrange_3D_shape_deriv<T>);
+            }
+
           default:
             libmesh_error_msg("ERROR: Unsupported 3D element type!: " << Utility::enum_to_string(type));
           }
diff --git a/src/fe/fe_monomial.C b/src/fe/fe_monomial.C
index 17a92811a43..b018c262505 100644
--- a/src/fe/fe_monomial.C
+++ b/src/fe/fe_monomial.C
@@ -77,6 +77,7 @@ unsigned int monomial_n_dofs(const ElemType t, const Order o)
           case PYRAMID5:
           case PYRAMID13:
           case PYRAMID14:
+          case PYRAMID18:
             return 4;
 
           case INVALID_ELEM:
@@ -127,6 +128,7 @@ unsigned int monomial_n_dofs(const ElemType t, const Order o)
           case PYRAMID5:
           case PYRAMID13:
           case PYRAMID14:
+          case PYRAMID18:
             return 10;
 
           case INVALID_ELEM:
@@ -177,6 +179,7 @@ unsigned int monomial_n_dofs(const ElemType t, const Order o)
           case PYRAMID5:
           case PYRAMID13:
           case PYRAMID14:
+          case PYRAMID18:
             return 20;
 
           case INVALID_ELEM:
diff --git a/src/geom/cell_pyramid18.C b/src/geom/cell_pyramid18.C
new file mode 100644
index 00000000000..b224023da01
--- /dev/null
+++ b/src/geom/cell_pyramid18.C
@@ -0,0 +1,509 @@
+// The libMesh Finite Element Library.
+// Copyright (C) 2002-2022 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
+
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
+
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
+
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
+
+
+// Local includes
+#include "libmesh/side.h"
+#include "libmesh/cell_pyramid18.h"
+#include "libmesh/edge_edge3.h"
+#include "libmesh/face_tri7.h"
+#include "libmesh/face_quad9.h"
+#include "libmesh/enum_io_package.h"
+#include "libmesh/enum_order.h"
+
+namespace libMesh
+{
+
+
+
+
+// ------------------------------------------------------------
+// Pyramid18 class static member initializations
+const int Pyramid18::num_nodes;
+const int Pyramid18::num_sides;
+const int Pyramid18::num_edges;
+const int Pyramid18::num_children;
+const int Pyramid18::nodes_per_side;
+const int Pyramid18::nodes_per_edge;
+
+const unsigned int Pyramid18::side_nodes_map[Pyramid18::num_sides][Pyramid18::nodes_per_side] =
+  {
+    {0, 1, 4, 5, 10,  9, 14, 99, 99}, // Side 0 (front)
+    {1, 2, 4, 6, 11, 10, 15, 99, 99}, // Side 1 (right)
+    {2, 3, 4, 7, 12, 11, 16, 99, 99}, // Side 2 (back)
+    {3, 0, 4, 8,  9, 12, 17, 99, 99}, // Side 3 (left)
+    {0, 3, 2, 1,  8,  7,  6,  5, 13}  // Side 4 (base)
+  };
+
+const unsigned int Pyramid18::edge_nodes_map[Pyramid18::num_edges][Pyramid18::nodes_per_edge] =
+  {
+    {0, 1,  5}, // Edge 0
+    {1, 2,  6}, // Edge 1
+    {2, 3,  7}, // Edge 2
+    {0, 3,  8}, // Edge 3
+    {0, 4,  9}, // Edge 4
+    {1, 4, 10}, // Edge 5
+    {2, 4, 11}, // Edge 6
+    {3, 4, 12}  // Edge 7
+  };
+
+// ------------------------------------------------------------
+// Pyramid18 class member functions
+
+bool Pyramid18::is_vertex(const unsigned int i) const
+{
+  if (i < 5)
+    return true;
+  return false;
+}
+
+
+
+bool Pyramid18::is_edge(const unsigned int i) const
+{
+  if (i < 5)
+    return false;
+  if (i > 12)
+    return false;
+  return true;
+}
+
+
+
+bool Pyramid18::is_face(const unsigned int i) const
+{
+  if (i > 12)
+    return true;
+  return false;
+}
+
+
+
+bool Pyramid18::is_node_on_side(const unsigned int n,
+                                const unsigned int s) const
+{
+  libmesh_assert_less (s, n_sides());
+  return std::find(std::begin(side_nodes_map[s]),
+                   std::end(side_nodes_map[s]),
+                   n) != std::end(side_nodes_map[s]);
+}
+
+std::vector<unsigned>
+Pyramid18::nodes_on_side(const unsigned int s) const
+{
+  libmesh_assert_less(s, n_sides());
+  auto trim = (s == 4) ? 0 : 2;
+  return {std::begin(side_nodes_map[s]), std::end(side_nodes_map[s]) - trim};
+}
+
+std::vector<unsigned>
+Pyramid18::nodes_on_edge(const unsigned int e) const
+{
+  libmesh_assert_less(e, n_edges());
+  return {std::begin(edge_nodes_map[e]), std::end(edge_nodes_map[e])};
+}
+
+bool Pyramid18::is_node_on_edge(const unsigned int n,
+                                const unsigned int e) const
+{
+  libmesh_assert_less (e, n_edges());
+  return std::find(std::begin(edge_nodes_map[e]),
+                   std::end(edge_nodes_map[e]),
+                   n) != std::end(edge_nodes_map[e]);
+}
+
+
+
+bool Pyramid18::has_affine_map() const
+{
+  // TODO: If the base is a parallelogram and all the triangular faces are planar,
+  // the map should be linear, but I need to test this theory...
+  return false;
+}
+
+
+
+Order Pyramid18::default_order() const
+{
+  return THIRD;
+}
+
+
+
+dof_id_type Pyramid18::key (const unsigned int s) const
+{
+  libmesh_assert_less (s, this->n_sides());
+
+  switch (s)
+    {
+    case 0: // triangular face 1
+    case 1: // triangular face 2
+    case 2: // triangular face 3
+    case 3: // triangular face 4
+      return this->compute_key (this->node_id(s+14));
+
+    case 4:  // the quad face at z=0
+      return this->compute_key (this->node_id(13));
+
+    default:
+      libmesh_error_msg("Invalid side s = " << s);
+    }
+}
+
+
+
+unsigned int Pyramid18::local_side_node(unsigned int side,
+                                        unsigned int side_node) const
+{
+  libmesh_assert_less (side, this->n_sides());
+
+  // Never more than 9 nodes per side.
+  libmesh_assert_less(side_node, Pyramid18::nodes_per_side);
+
+  // Some sides have 7 nodes.
+  libmesh_assert(side == 4 || side_node < 7);
+
+  return Pyramid18::side_nodes_map[side][side_node];
+}
+
+
+
+unsigned int Pyramid18::local_edge_node(unsigned int edge,
+                                        unsigned int edge_node) const
+{
+  libmesh_assert_less(edge, this->n_edges());
+  libmesh_assert_less(edge_node, Pyramid18::nodes_per_edge);
+
+  return Pyramid18::edge_nodes_map[edge][edge_node];
+}
+
+
+
+std::unique_ptr<Elem> Pyramid18::build_side_ptr (const unsigned int i, bool proxy)
+{
+  libmesh_assert_less (i, this->n_sides());
+
+  std::unique_ptr<Elem> face;
+  if (proxy)
+    {
+#ifdef LIBMESH_ENABLE_DEPRECATED
+      libmesh_deprecated();
+      switch (i)
+        {
+        case 0:
+        case 1:
+        case 2:
+        case 3:
+          {
+            face = std::make_unique<Side<Tri7,Pyramid18>>(this,i);
+            break;
+          }
+
+        case 4:
+          {
+            face = std::make_unique<Side<Quad9,Pyramid18>>(this,i);
+            break;
+          }
+
+        default:
+          libmesh_error_msg("Invalid side i = " << i);
+        }
+#else
+      libmesh_error();
+#endif // LIBMESH_ENABLE_DEPRECATED
+    }
+
+  else
+    {
+      switch (i)
+        {
+        case 0: // triangular face 1
+        case 1: // triangular face 2
+        case 2: // triangular face 3
+        case 3: // triangular face 4
+          {
+            face = std::make_unique<Tri7>();
+            break;
+          }
+        case 4: // the quad face at z=0
+          {
+            face = std::make_unique<Quad9>();
+            break;
+          }
+        default:
+          libmesh_error_msg("Invalid side i = " << i);
+        }
+
+      // Set the nodes
+      for (auto n : face->node_index_range())
+        face->set_node(n) = this->node_ptr(Pyramid18::side_nodes_map[i][n]);
+    }
+
+#ifdef LIBMESH_ENABLE_DEPRECATED
+  if (!proxy) // proxy sides used to leave parent() set
+#endif
+    face->set_parent(nullptr);
+  face->set_interior_parent(this);
+
+  face->subdomain_id() = this->subdomain_id();
+  face->set_mapping_type(this->mapping_type());
+#ifdef LIBMESH_ENABLE_AMR
+  face->set_p_level(this->p_level());
+#endif
+
+  return face;
+}
+
+
+
+void Pyramid18::build_side_ptr (std::unique_ptr<Elem> & side,
+                                const unsigned int i)
+{
+  libmesh_assert_less (i, this->n_sides());
+
+  switch (i)
+    {
+    case 0: // triangular face 1
+    case 1: // triangular face 2
+    case 2: // triangular face 3
+    case 3: // triangular face 4
+      {
+        if (!side.get() || side->type() != TRI7)
+          {
+            side = this->build_side_ptr(i, false);
+            return;
+          }
+        break;
+      }
+    case 4:  // the quad face at z=0
+      {
+        if (!side.get() || side->type() != QUAD9)
+          {
+            side = this->build_side_ptr(i, false);
+            return;
+          }
+        break;
+      }
+    default:
+      libmesh_error_msg("Invalid side i = " << i);
+    }
+
+  side->subdomain_id() = this->subdomain_id();
+  side->set_mapping_type(this->mapping_type());
+
+  // Set the nodes
+  for (auto n : side->node_index_range())
+    side->set_node(n) = this->node_ptr(Pyramid18::side_nodes_map[i][n]);
+}
+
+
+
+std::unique_ptr<Elem> Pyramid18::build_edge_ptr (const unsigned int i)
+{
+  return this->simple_build_edge_ptr<Edge3,Pyramid18>(i);
+}
+
+
+
+void Pyramid18::build_edge_ptr (std::unique_ptr<Elem> & edge, const unsigned int i)
+{
+  this->simple_build_edge_ptr<Pyramid18>(edge, i, EDGE3);
+}
+
+
+
+void Pyramid18::connectivity(const unsigned int libmesh_dbg_var(sc),
+                             const IOPackage iop,
+                             std::vector<dof_id_type> & /*conn*/) const
+{
+  libmesh_assert(_nodes);
+  libmesh_assert_less (sc, this->n_sub_elem());
+  libmesh_assert_not_equal_to (iop, INVALID_IO_PACKAGE);
+
+  switch (iop)
+    {
+    case TECPLOT:
+      {
+        // TODO
+        libmesh_not_implemented();
+      }
+
+    case VTK:
+      {
+        // TODO
+        libmesh_not_implemented();
+      }
+
+    default:
+      libmesh_error_msg("Unsupported IO package " << iop);
+    }
+}
+
+
+
+unsigned int Pyramid18::n_second_order_adjacent_vertices (const unsigned int n) const
+{
+  switch (n)
+    {
+    case 5:
+    case 6:
+    case 7:
+    case 8:
+    case 9:
+    case 10:
+    case 11:
+    case 12:
+      return 2;
+
+    case 13:
+      return 4;
+
+    case 14:
+    case 15:
+    case 16:
+    case 17:
+      return 3;
+
+    default:
+      libmesh_error_msg("Invalid node n = " << n);
+    }
+}
+
+
+unsigned short int Pyramid18::second_order_adjacent_vertex (const unsigned int n,
+                                                            const unsigned int v) const
+{
+  libmesh_assert_greater_equal (n, this->n_vertices());
+  libmesh_assert_less (n, this->n_nodes());
+
+  switch (n)
+    {
+    case 5:
+    case 6:
+    case 7:
+    case 8:
+    case 9:
+    case 10:
+    case 11:
+    case 12:
+      {
+        libmesh_assert_less (v, 2);
+
+        // This is the analog of the static, const arrays
+        // {Hex,Prism,Tet10}::_second_order_adjacent_vertices
+        // defined in the respective source files... possibly treat
+        // this similarly once the Pyramid13 has been added?
+        constexpr unsigned short node_list[8][2] =
+          {
+            {0,1},
+            {1,2},
+            {2,3},
+            {0,3},
+            {0,4},
+            {1,4},
+            {2,4},
+            {3,4}
+          };
+
+        return node_list[n-5][v];
+      }
+
+      // mid-face node on bottom
+    case 13:
+      {
+        libmesh_assert_less (v, 4);
+
+        // The vertex nodes surrounding node 13 are 0, 1, 2, and 3.
+        // Thus, the v'th node is simply = v.
+        return cast_int<unsigned short>(v);
+      }
+
+      // mid-face nodes on triangles
+    case 14:
+    case 15:
+    case 16:
+    case 17:
+      {
+        libmesh_assert_less (v, 3);
+
+        constexpr unsigned short node_list[4][3] =
+          {
+            {0,1,4},
+            {1,2,4},
+            {2,3,4},
+            {0,3,4}
+          };
+
+        return node_list[n-14][v];
+      }
+
+    default:
+      libmesh_error_msg("Invalid n = " << n);
+
+    }
+}
+
+
+
+void Pyramid18::permute(unsigned int perm_num)
+{
+  libmesh_assert_less (perm_num, 4);
+
+  for (unsigned int i = 0; i != perm_num; ++i)
+    {
+      swap4nodes(0,1,2,3);
+      swap4nodes(5,6,7,8);
+      swap4nodes(9,10,11,12);
+      swap4nodes(14,15,16,17);
+      swap4neighbors(0,1,2,3);
+    }
+}
+
+
+void Pyramid18::flip(BoundaryInfo * boundary_info)
+{
+  swap2nodes(0,1);
+  swap2nodes(2,3);
+  swap2nodes(6,8);
+  swap2nodes(9,10);
+  swap2nodes(11,12);
+  swap2nodes(15,17);
+  swap2neighbors(1,3);
+  swap2boundarysides(1,3,boundary_info);
+  swap2boundaryedges(1,3,boundary_info);
+  swap2boundaryedges(4,5,boundary_info);
+  swap2boundaryedges(6,7,boundary_info);
+}
+
+
+unsigned int Pyramid18::center_node_on_side(const unsigned short side) const
+{
+  libmesh_assert_less (side, Pyramid18::num_sides);
+  return side == 4 ? 13 : side+14;
+}
+
+
+ElemType Pyramid18::side_type (const unsigned int s) const
+{
+  libmesh_assert_less (s, 5);
+  if (s < 4)
+    return TRI7;
+  return QUAD9;
+}
+
+
+} // namespace libMesh
diff --git a/src/geom/elem.C b/src/geom/elem.C
index bec33a4a2ee..2dec756daee 100644
--- a/src/geom/elem.C
+++ b/src/geom/elem.C
@@ -58,6 +58,7 @@
 #include "libmesh/cell_pyramid5.h"
 #include "libmesh/cell_pyramid13.h"
 #include "libmesh/cell_pyramid14.h"
+#include "libmesh/cell_pyramid18.h"
 #include "libmesh/fe_base.h"
 #include "libmesh/mesh_base.h"
 #include "libmesh/quadrature_nodal.h"
@@ -336,6 +337,8 @@ std::unique_ptr<Elem> Elem::build(const ElemType type,
       return std::make_unique<Pyramid13>(p);
     case PYRAMID14:
       return std::make_unique<Pyramid14>(p);
+    case PYRAMID18:
+      return std::make_unique<Pyramid18>(p);
 
 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
       // 1D infinite elements
@@ -2647,6 +2650,7 @@ ElemType Elem::first_order_equivalent_type (const ElemType et)
     case PYRAMID5:
     case PYRAMID13:
     case PYRAMID14:
+    case PYRAMID18:
       return PYRAMID5;
 
 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
@@ -2774,6 +2778,8 @@ ElemType Elem::second_order_equivalent_type (const ElemType et,
       return PRISM20;
     case PRISM21:
       return PRISM21;
+    case PYRAMID18:
+      return PYRAMID18;
 
     case PYRAMID5:
     case PYRAMID13:
diff --git a/src/geom/elem_quality.C b/src/geom/elem_quality.C
index 63aaa62c4df..8effff28968 100644
--- a/src/geom/elem_quality.C
+++ b/src/geom/elem_quality.C
@@ -420,6 +420,7 @@ std::vector<ElemQuality> Quality::valid(const ElemType t)
     case PYRAMID5:
     case PYRAMID13:
     case PYRAMID14:
+    case PYRAMID18:
       {
         // None yet
         break;
diff --git a/src/geom/reference_elem.C b/src/geom/reference_elem.C
index 02852ce4af2..7c50d81a881 100644
--- a/src/geom/reference_elem.C
+++ b/src/geom/reference_elem.C
@@ -203,6 +203,7 @@ void init_ref_elem_table()
     ref_elem_file[PYRAMID5] = ElemDataStrings::one_pyramid;
     ref_elem_file[PYRAMID13] = ElemDataStrings::one_pyramid13;
     ref_elem_file[PYRAMID14] = ElemDataStrings::one_pyramid14;
+    ref_elem_file[PYRAMID18] = ElemDataStrings::one_pyramid18;
   }
 
   // Read'em
diff --git a/src/geom/reference_elem.data b/src/geom/reference_elem.data
index 118e91736fd..bccc91946da 100644
--- a/src/geom/reference_elem.data
+++ b/src/geom/reference_elem.data
@@ -1,7 +1,7 @@
 // THIS FILE IS AUTOMATICALLY GENERATED -- DO NOT EDIT
 // src/geom/reference_elem.data automatically created by
 // contrib/bin/generate_reference_elem.sh reference_elements/ src/geom/reference_elem.data
-// in directory /home/stogrh/libmesh/fresh
+// in directory /home/stogrh/libmesh/git
 
   typedef std::map<ElemType, std::string> FileMapType;
   FileMapType ref_elem_file;
@@ -585,6 +585,43 @@ static const char* one_pyramid14 =
 ;
 
 
+static const char* one_pyramid18 =
+"libMesh-0.7.0+\n"
+"1	 # number of elements\n"
+"18	 # number of nodes\n"
+"n/a	 # boundary condition specification file\n"
+".	 # 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 sid (n0 ... nN-1) ]\n"
+"37 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17\n"
+"-1.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00\n"
+"1.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00\n"
+"1.0000000000000000e+00 1.0000000000000000e+00 0.0000000000000000e+00\n"
+"-1.0000000000000000e+00 1.0000000000000000e+00 0.0000000000000000e+00\n"
+"0.0000000000000000e+00 0.0000000000000000e+00 1.0000000000000000e+00\n"
+"0.0000000000000000e+00 -1.0000000000000000e+00 0.0000000000000000e+00\n"
+"1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00\n"
+"0.0000000000000000e+00 1.0000000000000000e+00 0.0000000000000000e+00\n"
+"-1.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00\n"
+"-5.0000000000000000e-01 -5.0000000000000000e-01 5.0000000000000000e-01\n"
+"5.0000000000000000e-01 -5.0000000000000000e-01 5.0000000000000000e-01\n"
+"5.0000000000000000e-01 5.0000000000000000e-01 5.0000000000000000e-01\n"
+"-5.0000000000000000e-01 5.0000000000000000e-01 5.0000000000000000e-01\n"
+"0.0000000000000000e+00 0.0000000000000000e+00 0.0000000000000000e+00\n"
+"0.0000000000000000e+00 -6.6666666666666666e-01 3.3333333333333333e-01\n"
+"6.6666666666666666e-01 0.0000000000000000e+00 3.3333333333333333e-01\n"
+"0.0000000000000000e+00 6.6666666666666666e-01 3.3333333333333333e-01\n"
+"-6.6666666666666666e-01 0.0000000000000000e+00 3.3333333333333333e-01\n"
+"5	 # number of boundary conditions\n"
+"0 0 0\n"
+"0 1 1\n"
+"0 2 2\n"
+"0 3 3\n"
+"0 4 4\n"
+;
+
+
 static const char* one_tet =
 "libMesh-0.7.0+\n"
 "1	 # number of elements\n"
diff --git a/src/libmesh_SOURCES b/src/libmesh_SOURCES
index 2b2b05a2b7c..f1d2d5217a7 100644
--- a/src/libmesh_SOURCES
+++ b/src/libmesh_SOURCES
@@ -141,6 +141,7 @@ libmesh_SOURCES =  \
         src/geom/cell_pyramid.C \
         src/geom/cell_pyramid13.C \
         src/geom/cell_pyramid14.C \
+        src/geom/cell_pyramid18.C \
         src/geom/cell_pyramid5.C \
         src/geom/cell_tet.C \
         src/geom/cell_tet10.C \
diff --git a/src/mesh/exodusII_io_helper.C b/src/mesh/exodusII_io_helper.C
index c51c47bb79a..aea3ccd3a17 100644
--- a/src/mesh/exodusII_io_helper.C
+++ b/src/mesh/exodusII_io_helper.C
@@ -542,6 +542,13 @@ void ExodusII_IO_Helper::init_conversion_map()
     conv.n_nodes = 14;
     conv.exodus_type = "PYRAMID14";
   }
+  {
+    auto & conv = conversion_map[PYRAMID18];
+    conv.libmesh_type = PYRAMID18;
+    conv.dim = 3;
+    conv.n_nodes = 18;
+    conv.exodus_type = "PYRAMID18";
+  }
 }
 
 
@@ -652,15 +659,12 @@ void ExodusII_IO_Helper::init_element_equivalence_map()
   // PRISM21 equivalences
   element_equivalence_map["WEDGE21"] = PRISM21;
 
-  // PYRAMID5 equivalences
+  // PYRAMID equivalences
   element_equivalence_map["PYRAMID"]  = PYRAMID5;
   element_equivalence_map["PYRAMID5"] = PYRAMID5;
-
-  // PYRAMID13 equivalences
   element_equivalence_map["PYRAMID13"] = PYRAMID13;
-
-  // PYRAMID14 equivalences
   element_equivalence_map["PYRAMID14"] = PYRAMID14;
+  element_equivalence_map["PYRAMID18"] = PYRAMID18;
 }
 
 const ExodusII_IO_Helper::Conversion &
diff --git a/src/mesh/mesh_generation.C b/src/mesh/mesh_generation.C
index 3360b6276aa..f1cc6235004 100644
--- a/src/mesh/mesh_generation.C
+++ b/src/mesh/mesh_generation.C
@@ -134,6 +134,7 @@ unsigned int idx(const ElemType type,
     case PYRAMID5: // PYRAMID5's are created from an initial HEX27 discretization
     case PYRAMID13:
     case PYRAMID14:
+    case PYRAMID18:
     case PRISM15:
     case PRISM18:
     case PRISM20:
@@ -960,6 +961,7 @@ void MeshTools::Generation::build_cube(UnstructuredMesh & mesh,
           case PYRAMID5: // PYRAMIDs are created from an initial HEX27 discretization
           case PYRAMID13:
           case PYRAMID14:
+          case PYRAMID18:
             {
               mesh.reserve_elem(nx*ny*nz);
               break;
@@ -1001,6 +1003,7 @@ void MeshTools::Generation::build_cube(UnstructuredMesh & mesh,
           case PYRAMID5: // PYRAMIDs are created from an initial HEX27 discretization
           case PYRAMID13:
           case PYRAMID14:
+          case PYRAMID18:
           case PRISM15:
           case PRISM18:
             {
@@ -1083,6 +1086,7 @@ void MeshTools::Generation::build_cube(UnstructuredMesh & mesh,
           case PYRAMID5: // PYRAMIDs are created from an initial HEX27 discretization
           case PYRAMID13:
           case PYRAMID14:
+          case PYRAMID18:
           case PRISM15:
           case PRISM18:
           case PRISM20:
@@ -1263,6 +1267,7 @@ void MeshTools::Generation::build_cube(UnstructuredMesh & mesh,
           case PYRAMID5: // PYRAMIDs are created from an initial HEX27 discretization
           case PYRAMID13:
           case PYRAMID14:
+          case PYRAMID18:
             {
               for (unsigned int k=0; k<(2*nz); k += 2)
                 for (unsigned int j=0; j<(2*ny); j += 2)
@@ -1293,7 +1298,8 @@ void MeshTools::Generation::build_cube(UnstructuredMesh & mesh,
                       elem->set_node(19) = mesh.node_ptr(idx(type,nx,ny,i,  j+1,k+2));
 
                       if ((type == HEX27) || (type == TET4) || (type == TET10) || (type == TET14) ||
-                          (type == PYRAMID5) || (type == PYRAMID13) || (type == PYRAMID14))
+                          (type == PYRAMID5) || (type == PYRAMID13) || (type == PYRAMID14) ||
+                          (type == PYRAMID18))
                         {
                           elem->set_node(20) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k)  );
                           elem->set_node(21) = mesh.node_ptr(idx(type,nx,ny,i+1,j,  k+1));
@@ -1496,7 +1502,8 @@ void MeshTools::Generation::build_cube(UnstructuredMesh & mesh,
             (type == TET14) ||
             (type == PYRAMID5) ||
             (type == PYRAMID13) ||
-            (type == PYRAMID14))
+            (type == PYRAMID14) ||
+            (type == PYRAMID18))
           {
             // Temporary storage for new elements. (24 tets per hex, 6 pyramids)
             std::vector<std::unique_ptr<Elem>> new_elements;
@@ -1552,7 +1559,7 @@ void MeshTools::Generation::build_cube(UnstructuredMesh & mesh,
                           }
                       } // end if ((type == TET4) || (type == TET10) || (type == TET14))
 
-                    else // type==PYRAMID5 || type==PYRAMID13 || type==PYRAMID14
+                    else // type==PYRAMID*
                       {
                         // Build 1 sub-pyramid per side.
                         new_elements.push_back( Elem::build(PYRAMID5) );
@@ -1573,7 +1580,7 @@ void MeshTools::Generation::build_cube(UnstructuredMesh & mesh,
                         // 4 (the square base) with the same b_id.
                         if (b_id != BoundaryInfo::invalid_id)
                           boundary_info.add_side(sub_elem.get(), 4, b_id);
-                      } // end else type==PYRAMID5 || type==PYRAMID13 || type==PYRAMID14
+                      } // end else type==PYRAMID*
                   }
               }
 
@@ -1592,7 +1599,7 @@ void MeshTools::Generation::build_cube(UnstructuredMesh & mesh,
                 mesh.add_elem( std::move(new_elements[i]) );
               }
 
-          } // end if (type == TET4,TET10,TET14,PYRAMID5,PYRAMID13,PYRAMID14
+          } // end if (type == TET*,PYRAMID*)
 
 
         // Use all_second_order to convert the TET4's to TET10's or PYRAMID5's to PYRAMID14's
@@ -1602,7 +1609,7 @@ void MeshTools::Generation::build_cube(UnstructuredMesh & mesh,
         else if (type == PYRAMID13)
           mesh.all_second_order(/*full_ordered=*/false);
 
-        else if (type == TET14)
+        else if ((type == TET14) || (type == PYRAMID18))
           mesh.all_complete_order();
 
 
diff --git a/src/mesh/unstructured_mesh.C b/src/mesh/unstructured_mesh.C
index e2f8851c5d7..6dbec0afd59 100644
--- a/src/mesh/unstructured_mesh.C
+++ b/src/mesh/unstructured_mesh.C
@@ -1631,6 +1631,7 @@ void UnstructuredMesh::all_complete_order ()
           lo_elem->type() != TET14 &&
           lo_elem->type() != PRISM20 &&
           lo_elem->type() != PRISM21 &&
+          lo_elem->type() != PYRAMID18 &&
           lo_elem->default_order() != FIRST)
         continue;
 
@@ -1650,6 +1651,10 @@ void UnstructuredMesh::all_complete_order ()
                lo_elem->type() == PRISM15 ||
                lo_elem->type() == PRISM6)
         hi_type = PRISM21;
+      else if (lo_elem->type() == PYRAMID14 ||
+               lo_elem->type() == PYRAMID13 ||
+               lo_elem->type() == PYRAMID5)
+        hi_type = PYRAMID18;
       else
         libmesh_assert(lo_elem->type() == TRI6 || lo_elem->type() == TRI3);
 
diff --git a/src/quadrature/quadrature_conical_3D.C b/src/quadrature/quadrature_conical_3D.C
index b01ed469c15..44f9d89e683 100644
--- a/src/quadrature/quadrature_conical_3D.C
+++ b/src/quadrature/quadrature_conical_3D.C
@@ -38,6 +38,7 @@ void QConical::init_3D(const ElemType, unsigned int)
     case PYRAMID5:
     case PYRAMID13:
     case PYRAMID14:
+    case PYRAMID18:
       {
         this->conical_product_pyramid();
         return;
diff --git a/src/quadrature/quadrature_gauss_3D.C b/src/quadrature/quadrature_gauss_3D.C
index d7409d9e7ce..448f2bad411 100644
--- a/src/quadrature/quadrature_gauss_3D.C
+++ b/src/quadrature/quadrature_gauss_3D.C
@@ -531,6 +531,7 @@ void QGauss::init_3D(const ElemType, unsigned int)
     case PYRAMID5:
     case PYRAMID13:
     case PYRAMID14:
+    case PYRAMID18:
       {
         // We compute the Pyramid rule as a conical product of a
         // Jacobi rule with alpha==2 on the interval [0,1] two 1D
diff --git a/src/quadrature/quadrature_grid_3D.C b/src/quadrature/quadrature_grid_3D.C
index b0f1147dd1a..c5a8e6fa2a4 100644
--- a/src/quadrature/quadrature_grid_3D.C
+++ b/src/quadrature/quadrature_grid_3D.C
@@ -111,6 +111,7 @@ void QGrid::init_3D(const ElemType, unsigned int)
     case PYRAMID5:
     case PYRAMID13:
     case PYRAMID14:
+    case PYRAMID18:
       {
         const unsigned int np = (_order+1)*(_order+2)*(_order+3)/6;
         _points.resize(np);
diff --git a/src/quadrature/quadrature_nodal_3D.C b/src/quadrature/quadrature_nodal_3D.C
index 09f89f85494..3b1d1d23aa8 100644
--- a/src/quadrature/quadrature_nodal_3D.C
+++ b/src/quadrature/quadrature_nodal_3D.C
@@ -202,6 +202,33 @@ void QNodal::init_3D(const ElemType, unsigned int)
         return;
       }
 
+    case PYRAMID18:
+      {
+        libmesh_error_msg_if(!allow_nodal_pyramid_quadrature,
+                             "Nodal quadrature on Pyramid elements is not allowed by default since\n"
+                             "the Jacobian of the inverse element map is not well-defined at the Pyramid apex.\n"
+                             "Set the QBase::allow_nodal_pyramid_quadrature flag to true to ignore skip this check.");
+
+        _points =
+          {
+            Point(-1,-1,0.), Point(+1,-1,0.), Point(+1,+1,0.),
+            Point(-1,+1,0.), Point(0.,0.,+1), Point(0.,-1,0.),
+            Point(+1,0.,0.), Point(0.,+1,0.), Point(-1,0.,0.),
+            Point(-.5,-.5,.5), Point(-.5,.5,.5), Point(.5,.5,.5),
+            Point(-.5,.5,.5), Point(0.,0.,0.), Point(0.,-2/Real(3),1/Real(3)),
+            Point(2/Real(3),0.,1/Real(3)), Point(0.,2/Real(3),1/Real(3)),
+            Point(-2/Real(3),0.,1/Real(3))
+          };
+
+        // Even with triangle faces to play with, I can't seem to get
+        // exact integrals of any higher order functions without
+        // negative weights.  So I punt and just use QTrap weights.
+        _weights = {1/Real(4), 1/Real(4), 1/Real(4), 1/Real(4), 1/Real(3),
+                    0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
+
+        return;
+      }
+
     case TET14:
       {
         _points.resize(14);
diff --git a/src/quadrature/quadrature_simpson_3D.C b/src/quadrature/quadrature_simpson_3D.C
index 9e3e14bc2b1..bd77b1e7d4a 100644
--- a/src/quadrature/quadrature_simpson_3D.C
+++ b/src/quadrature/quadrature_simpson_3D.C
@@ -139,6 +139,7 @@ void QSimpson::init_3D(const ElemType, unsigned int)
     case PYRAMID5:
     case PYRAMID13:
     case PYRAMID14:
+    case PYRAMID18:
       {
         libmesh_error_msg_if(!allow_nodal_pyramid_quadrature,
                              "Nodal quadrature on Pyramid elements is not allowed by default since\n"
diff --git a/src/quadrature/quadrature_trap_3D.C b/src/quadrature/quadrature_trap_3D.C
index 32e7fd4210e..4126491a67e 100644
--- a/src/quadrature/quadrature_trap_3D.C
+++ b/src/quadrature/quadrature_trap_3D.C
@@ -91,6 +91,7 @@ void QTrap::init_3D(const ElemType, unsigned int)
     case PYRAMID5:
     case PYRAMID13:
     case PYRAMID14:
+    case PYRAMID18:
       {
         libmesh_error_msg_if(!allow_nodal_pyramid_quadrature,
                              "Nodal quadrature on Pyramid elements is not allowed by default since\n"
diff --git a/tests/fe/fe_lagrange_test.C b/tests/fe/fe_lagrange_test.C
index 4863bf14511..acc3fa6550d 100644
--- a/tests/fe/fe_lagrange_test.C
+++ b/tests/fe/fe_lagrange_test.C
@@ -38,4 +38,5 @@ INSTANTIATE_FETEST(THIRD, LAGRANGE, PRISM21);
 INSTANTIATE_FETEST(FIRST, LAGRANGE, PYRAMID5);
 INSTANTIATE_FETEST(SECOND, LAGRANGE, PYRAMID13);
 INSTANTIATE_FETEST(SECOND, LAGRANGE, PYRAMID14);
+INSTANTIATE_FETEST(THIRD, LAGRANGE, PYRAMID18);
 #endif
diff --git a/tests/fe/fe_test.h b/tests/fe/fe_test.h
index e705c2ac37a..2adb76413de 100644
--- a/tests/fe/fe_test.h
+++ b/tests/fe/fe_test.h
@@ -254,7 +254,8 @@ Gradient rational_test_grad (const Point& p,
 #define FE_CAN_TEST_CUBIC \
   (((family != LAGRANGE && family != L2_LAGRANGE) || \
     (elem_type != TRI7 && elem_type != TET14 && \
-     elem_type != PRISM20 && elem_type != PRISM21)) && order > 2)
+     elem_type != PRISM20 && elem_type != PRISM21 && \
+     elem_type != PYRAMID18)) && order > 2)
 
 
 template <Order order, FEFamily family, ElemType elem_type, unsigned int build_nx>
diff --git a/tests/geom/elem_test.C b/tests/geom/elem_test.C
index 2795448558c..c4a2d6b8275 100644
--- a/tests/geom/elem_test.C
+++ b/tests/geom/elem_test.C
@@ -253,8 +253,9 @@ public:
               new_quasicc = elem->quasicircumcenter();
             for (const auto d : make_range(LIBMESH_DIM))
               {
+                // Getting a little FP error from Pyramid18
                 LIBMESH_ASSERT_FP_EQUAL(centroid(d), new_centroid(d),
-                                        TOLERANCE*TOLERANCE);
+                                        TOLERANCE*std::sqrt(TOLERANCE));
                 LIBMESH_ASSERT_FP_EQUAL(vertex_avg(d), new_vertex_avg(d),
                                         TOLERANCE*TOLERANCE);
                 LIBMESH_ASSERT_FP_EQUAL(quasicc(d), new_quasicc(d),
@@ -474,6 +475,13 @@ public:
             CPPUNIT_ASSERT_EQUAL((unsigned int)(18 + (s == 4)), elem->center_node_on_side(s));
           else if (elem->type() == PYRAMID14 && s == 4)
             CPPUNIT_ASSERT_EQUAL((unsigned int)(13), elem->center_node_on_side(s));
+          else if (elem->type() == PYRAMID18)
+            {
+              if (s < 4)
+                CPPUNIT_ASSERT_EQUAL((unsigned int)(s + 14), elem->center_node_on_side(s));
+              else
+                CPPUNIT_ASSERT_EQUAL((unsigned int)(13), elem->center_node_on_side(s));
+            }
           else
             CPPUNIT_ASSERT_EQUAL(invalid_uint, elem->center_node_on_side(s));
         }
@@ -577,6 +585,7 @@ INSTANTIATE_ELEMTEST(PRISM21);
 INSTANTIATE_ELEMTEST(PYRAMID5);
 INSTANTIATE_ELEMTEST(PYRAMID13);
 INSTANTIATE_ELEMTEST(PYRAMID14);
+INSTANTIATE_ELEMTEST(PYRAMID18);
 
 #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
 INSTANTIATE_ELEMTEST(INFHEX8);
diff --git a/tests/geom/which_node_am_i_test.C b/tests/geom/which_node_am_i_test.C
index 0521fa0e23f..6bd48e96b52 100644
--- a/tests/geom/which_node_am_i_test.C
+++ b/tests/geom/which_node_am_i_test.C
@@ -76,6 +76,17 @@ public:
     for (unsigned int edge=4; edge<8; ++edge)
       CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(4),
                            pyr14.local_edge_node(edge, /*node=*/1));
+
+    const Elem & pyr18 = ReferenceElem::get(PYRAMID18);
+    // The central node of the base should be node 13
+    CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(13), pyr18.local_side_node(/*side=*/4, /*node=*/8));
+    // The central node of the sides should be nodes 14-17
+    for (unsigned int side=0; side<4; ++side)
+      CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(14+side), pyr18.local_side_node(side, /*node=*/6));
+    // Second node of edges 4-7 is the apex node (4).
+    for (unsigned int edge=4; edge<8; ++edge)
+      CPPUNIT_ASSERT_EQUAL(static_cast<unsigned int>(4),
+                           pyr18.local_edge_node(edge, /*node=*/1));
   }
 
 
diff --git a/tests/mesh/mesh_generation_test.C b/tests/mesh/mesh_generation_test.C
index e0c1092e4fc..987fedd71d4 100644
--- a/tests/mesh/mesh_generation_test.C
+++ b/tests/mesh/mesh_generation_test.C
@@ -188,8 +188,8 @@ public:
         CPPUNIT_ASSERT_EQUAL(mesh.n_nodes(),
                              cast_int<dof_id_type>((2*n+1)*(2*n+1)*(2*n+1) + 14*n*n*n + 4*3*(n+1)*n*n));
         break;
-      case 18: // Second-order prisms and hexes use the same nodes
-      case 27:
+      case 18:
+      case 27: // Second-order prisms and hexes use the same nodes
         CPPUNIT_ASSERT_EQUAL(mesh.n_nodes(),
                              cast_int<dof_id_type>((2*n+1)*(2*n+1)*(2*n+1)));
         break;
diff --git a/tests/quadrature/quadrature_test.C b/tests/quadrature/quadrature_test.C
index 13632badfec..18e185abaec 100644
--- a/tests/quadrature/quadrature_test.C
+++ b/tests/quadrature/quadrature_test.C
@@ -180,7 +180,8 @@ private:
     // We are testing integration on the reference elements here, so
     // the element map is not relevant, and we can safely use nodal
     // Pyramid quadrature.
-    if (elem_type == PYRAMID5 || elem_type == PYRAMID13 || elem_type == PYRAMID14 || elem_type == PYRAMID18)
+    if (elem_type == PYRAMID5 || elem_type == PYRAMID13 ||
+        elem_type == PYRAMID14 || elem_type == PYRAMID18)
       qrule->allow_nodal_pyramid_quadrature = true;
 
     qrule->init (elem_type);
@@ -333,7 +334,7 @@ public:
        {TET4, TET10, TET14},
        {HEX8, HEX20, HEX27},
        {PRISM6, PRISM15, PRISM18, PRISM20, PRISM21},
-       {PYRAMID5, PYRAMID13, PYRAMID14}};
+       {PYRAMID5, PYRAMID13, PYRAMID14, PYRAMID18}};
 
     const std::function<Real(int,int,int)> true_values[] =
       {edge_integrals,
@@ -364,6 +365,11 @@ public:
               elem_type == TET14 || elem_type == PRISM20)
             exactorder--;
 
+          // And one partially-cubic element can only do linears.
+          // Seriously.
+          if (elem_type == PYRAMID18)
+            exactorder=1;
+
           testPolynomials(QNODAL, order, elem_type, true_values[i], exactorder);
         }
   }