diff --git a/fem/fe/fe_h1.cpp b/fem/fe/fe_h1.cpp
index 2a0474fbdb..3fc55ea2bc 100644
--- a/fem/fe/fe_h1.cpp
+++ b/fem/fe/fe_h1.cpp
@@ -1107,7 +1107,7 @@ H1_FuentesPyramidElement::H1_FuentesPyramidElement(const int p, const int btype)
    {
       for (int i = 1; i < p; i++)
       {
-         Nodes.IntPoint(o++).Set3(cp[i], cp[j], cp[0]);
+         Nodes.IntPoint(o++).Set3(cp[i], cp[p-j], cp[0]);
       }
    }
 
diff --git a/fem/fe/fe_nd.cpp b/fem/fe/fe_nd.cpp
index 348b194c1f..fd1fe15e8d 100644
--- a/fem/fe/fe_nd.cpp
+++ b/fem/fe/fe_nd.cpp
@@ -1684,7 +1684,7 @@ ND_FuentesPyramidElement::ND_FuentesPyramidElement(const int p,
    for (int j = 1; j < p; j++)
       for (int i = 0; i < p; i++)
       {
-         Nodes.IntPoint(o).Set3(qop[i], qcp[j], 0.);
+         Nodes.IntPoint(o).Set3(qop[i], qcp[p-j], 0.);
          dof2tk[o++] = 0; // (1 0 0)
       }
 
@@ -1692,7 +1692,7 @@ ND_FuentesPyramidElement::ND_FuentesPyramidElement(const int p,
    for (int j = 0; j < p; j++)
       for (int i = 1; i < p; i++)
       {
-         Nodes.IntPoint(o).Set3(qcp[i], qop[j], 0.);
+         Nodes.IntPoint(o).Set3(qcp[i], qop[p-1-j], 0.);
          dof2tk[o++] = 7; // (0 -1 0)
       }
 
diff --git a/fem/fe/fe_rt.cpp b/fem/fe/fe_rt.cpp
index 613f6d0792..7f69eb1c9d 100644
--- a/fem/fe/fe_rt.cpp
+++ b/fem/fe/fe_rt.cpp
@@ -1317,7 +1317,7 @@ RT_FuentesPyramidElement::RT_FuentesPyramidElement(const int p)
    for (int j = 0; j <= p; j++)
       for (int i = 0; i <= p; i++)  // (3,2,1,0)
       {
-         Nodes.IntPoint(o).Set3(bop[i], bop[j], 0.);
+         Nodes.IntPoint(o).Set3(bop[i], bop[p-j], 0.);
          dof2nk[o++] = 0;
       }
    // triangular faces
diff --git a/mesh/mesh.cpp b/mesh/mesh.cpp
index 66df625fce..f9ac5faa09 100644
--- a/mesh/mesh.cpp
+++ b/mesh/mesh.cpp
@@ -3583,6 +3583,7 @@ void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
    }
 
 #undef VTX
+#undef VTXP
 
 #if 0
    ofstream test_stream("debug.mesh");
diff --git a/tests/unit/fem/test_3d_bilininteg.cpp b/tests/unit/fem/test_3d_bilininteg.cpp
index 948bd5c1bb..ab62f922ea 100644
--- a/tests/unit/fem/test_3d_bilininteg.cpp
+++ b/tests/unit/fem/test_3d_bilininteg.cpp
@@ -3348,7 +3348,7 @@ TEST_CASE("3D Bilinear Weak Gradient Integrators",
    VectorFunctionCoefficient  df3_coef(dim, Grad_f3);
    VectorFunctionCoefficient dqf3_coef(dim, Grad_qf3);
 
-   for (int type = first_3D_et; type <= last_3D_et-1; type++) // MLS: Fix this
+   for (int type = first_3D_et; type <= last_3D_et; type++)
    {
       Mesh mesh =
          Mesh::MakeCartesian3D(n, n, n, (Element::Type)type, 2.0, 3.0, 5.0);
@@ -3555,7 +3555,7 @@ TEST_CASE("3D Bilinear Scalar Weak Divergence Integrators",
    VectorFunctionCoefficient  Vf3_coef(dim, Vf3);
    FunctionCoefficient       dVf3_coef(Div_Vf3);
 
-   for (int type = first_3D_et; type <= last_3D_et-1; type++) // MLS: Fix this
+   for (int type = first_3D_et; type <= last_3D_et; type++)
    {
       Mesh mesh =
          Mesh::MakeCartesian3D(n, n, n, (Element::Type)type, 2.0, 3.0, 5.0);
@@ -3708,7 +3708,7 @@ TEST_CASE("3D Bilinear Weak Divergence Integrators",
    FunctionCoefficient      dDF3_coef(Div_DF3);
    FunctionCoefficient      dMF3_coef(Div_MF3);
 
-   for (int type = first_3D_et; type <= last_3D_et-1; type++) // MLS: Fix this
+   for (int type = first_3D_et; type <= last_3D_et; type++)
    {
       Mesh mesh =
          Mesh::MakeCartesian3D(n, n, n, (Element::Type)type, 2.0, 3.0, 5.0);
@@ -4043,7 +4043,7 @@ TEST_CASE("3D Bilinear Weak Curl Integrators",
    VectorFunctionCoefficient dDF3_coef(dim, Curl_DF3);
    VectorFunctionCoefficient dMF3_coef(dim, Curl_MF3);
 
-   for (int type = first_3D_et; type <= last_3D_et-1; type++) // MLS: Fix this
+   for (int type = first_3D_et; type <= last_3D_et; type++)
    {
       Mesh mesh =
          Mesh::MakeCartesian3D(n, n, n, (Element::Type)type, 2.0, 3.0, 5.0);
@@ -4365,7 +4365,7 @@ TEST_CASE("3D Bilinear Weak Div Cross Integrators",
    VectorFunctionCoefficient  VF3_coef(dim, VcrossF3);
    FunctionCoefficient       dVF3_coef(Div_VcrossF3);
 
-   for (int type = first_3D_et; type <= last_3D_et-1; type++) // MLS: Fix this
+   for (int type = first_3D_et; type <= last_3D_et; type++)
    {
       Mesh mesh =
          Mesh::MakeCartesian3D(n, n, n, (Element::Type)type, 2.0, 3.0, 5.0);
@@ -4500,7 +4500,7 @@ TEST_CASE("3D Bilinear Weak Curl Cross Integrators",
    VectorFunctionCoefficient  VxF3_coef(dim, VcrossF3);
    VectorFunctionCoefficient dVxF3_coef(dim, Curl_VcrossF3);
 
-   for (int type = first_3D_et; type <= last_3D_et-1; type++) // MLS: Fix this
+   for (int type = first_3D_et; type <= last_3D_et; type++)
    {
       Mesh mesh =
          Mesh::MakeCartesian3D(n, n, n, (Element::Type)type, 2.0, 3.0, 5.0);
@@ -4632,7 +4632,7 @@ TEST_CASE("3D Bilinear Weak Grad Dot Product Integrators",
    FunctionCoefficient       VdotF3_coef(VdotF3);
    VectorFunctionCoefficient   dVF3_coef(dim, GradVdotF3);
 
-   for (int type = first_3D_et; type <= last_3D_et-1; type++) // MLS: Fix this
+   for (int type = first_3D_et; type <= last_3D_et; type++)
    {
       Mesh mesh =
          Mesh::MakeCartesian3D(n, n, n, (Element::Type)type, 2.0, 3.0, 5.0);
@@ -4769,7 +4769,7 @@ TEST_CASE("3D Bilinear Grad Div Integrators",
    VectorFunctionCoefficient dVdf3_coef(dim, GradVdotGrad_f3);
    FunctionCoefficient       dVdF3_coef(DivVDivF3);
 
-   for (int type = first_3D_et; type <= last_3D_et-1; type++) // MLS: Fix this
+   for (int type = first_3D_et; type <= last_3D_et; type++)
    {
       Mesh mesh =
          Mesh::MakeCartesian3D(n, n, n, (Element::Type)type, 2.0, 3.0, 5.0);
@@ -5344,7 +5344,7 @@ TEST_CASE("3D Bilinear Curl Curl Integrators",
    VectorFunctionCoefficient  qdF3_coef(dim, qCurlF3);
    VectorFunctionCoefficient dqdF3_coef(dim, Curl_qCurlF3);
 
-   for (int type = first_3D_et; type <= last_3D_et-1; type++) // MLS: Fix this
+   for (int type = first_3D_et; type <= last_3D_et; type++)
    {
       Mesh mesh =
          Mesh::MakeCartesian3D(n, n, n, (Element::Type)type, 2.0, 3.0, 5.0);
@@ -5489,7 +5489,7 @@ TEST_CASE("3D Bilinear Mixed Curl Curl Integrators",
    VectorFunctionCoefficient dDdF3_coef(dim, Curl_DCurlF3);
    VectorFunctionCoefficient dMdF3_coef(dim, Curl_MCurlF3);
 
-   for (int type = first_3D_et; type <= last_3D_et-1; type++) // MLS: Fix this
+   for (int type = first_3D_et; type <= last_3D_et; type++)
    {
       Mesh mesh =
          Mesh::MakeCartesian3D(n, n, n, (Element::Type)type, 2.0, 3.0, 5.0);
@@ -5645,7 +5645,7 @@ TEST_CASE("3D Bilinear Mixed Cross Curl Curl Integrators",
    VectorFunctionCoefficient  VdF3_coef(dim, VcrossCurlF3);
    VectorFunctionCoefficient dVdF3_coef(dim, Curl_VcrossCurlF3);
 
-   for (int type = first_3D_et; type <= last_3D_et-1; type++) // MLS: Fix this
+   for (int type = first_3D_et; type <= last_3D_et; type++)
    {
       Mesh mesh =
          Mesh::MakeCartesian3D(n, n, n, (Element::Type)type, 2.0, 3.0, 5.0);
@@ -5716,7 +5716,7 @@ TEST_CASE("3D Bilinear Mixed Cross Grad Curl Integrators",
    VectorFunctionCoefficient  Vdf3_coef(dim, VcrossGrad_f3);
    VectorFunctionCoefficient dVdf3_coef(dim, Curl_VcrossGrad_f3);
 
-   for (int type = first_3D_et; type <= last_3D_et-1; type++) // MLS: Fix this
+   for (int type = first_3D_et; type <= last_3D_et; type++)
    {
       Mesh mesh =
          Mesh::MakeCartesian3D(n, n, n, (Element::Type)type, 2.0, 3.0, 5.0);
@@ -5801,7 +5801,7 @@ TEST_CASE("3D Bilinear Div Div Integrators",
    FunctionCoefficient         qdF3_coef(qDivF3);
    VectorFunctionCoefficient  dqdF3_coef(dim, Grad_qDivF3);
 
-   for (int type = first_3D_et; type <= last_3D_et-1; type++) // MLS: Fix this
+   for (int type = first_3D_et; type <= last_3D_et; type++)
    {
       Mesh mesh =
          Mesh::MakeCartesian3D(n, n, n, (Element::Type)type, 2.0, 3.0, 5.0);