Extending a unit test to include pyramids

fem/fe/fe_base.cpp

@@ -1993,6 +1993,57 @@ void Poly_1D::CalcDBinomTerms(const int p, const double x, const double y,
    }
 }
 
+void Poly_1D::CalcJacobi(const int p, const double alpha, const double beta,
+                         const double x, double *u)
+{
+   // use the recursive definition for [-1,1]:
+   // 2(n+1)(n+a+b+1)(2n+a+b)P^{(a,b)}_{n+1}(z) =
+   //    (2n+a+b+1)((2n+a+b+2)(2n+a+b)z+a^2-b^2)P^{(a,b)}_n(z)
+   //    -2(n+a)(n+b)(2n+a+b+2)P^{(a,b)}_{n-1}(z)
+   u[0] = 1.;
+   if (p == 0) { return; }
+   double z = 2.*x - 1.;
+   u[1] = 0.5 * ((alpha + beta + 2.) * z + alpha - beta);
+   for (int n = 1; n < p; n++)
+   {
+      double c2 = 2.*(1. + n)*(1. + n + alpha + beta)*(2.*n + alpha + beta);
+      double c1 = (1. + 2.*n + alpha + beta) *
+                  ((2. + 2.*n + alpha + beta)*
+                   (2.*n + alpha + beta)*z + pow(alpha,2) - pow(beta,2));
+      double c0 = 2.*(alpha + n)*(beta + n)*(2. + 2.*n + alpha + beta);
+
+      u[n+1] = (c1 * u[n] - c0 * u[n-1]) / c2;
+   }
+}
+
+void Poly_1D::CalcJacobi(const int p, const double alpha, const double beta,
+                         const double x, double *u, double *d)
+{
+   // use the recursive definition for [-1,1]:
+   // 2(n+1)(n+a+b+1)(2n+a+b)P^{(a,b)}_{n+1}(z) =
+   //    (2n+a+b+1)((2n+a+b+2)(2n+a+b)z+a^2-b^2)P^{(a,b)}_n(z)
+   //    -2(n+a)(n+b)(2n+a+b+2)P^{(a,b)}_{n-1}(z)
+   u[0] = 1.;
+   d[0] = 0.;
+   if (p == 0) { return; }
+   double z = 2.*x - 1.;
+   u[1] = 0.5 * ((alpha + beta + 2.) * z + alpha - beta);
+   d[1] = alpha + beta + 2.;
+   for (int n = 1; n < p; n++)
+   {
+      double c2 = 2.*(1. + n)*(1. + n + alpha + beta)*(2.*n + alpha + beta);
+      double c1 = (1. + 2.*n + alpha + beta) *
+                  ((2. + 2.*n + alpha + beta)*
+                   (2.*n + alpha + beta)*z + pow(alpha,2) - pow(beta,2));
+      double c0 = 2.*(alpha + n)*(beta + n)*(2. + 2.*n + alpha + beta);
+      double dc1 = 2.*(1. + 2.*n + alpha + beta) *
+                   (2. + 2.*n + alpha + beta)*(2.*n + alpha + beta);
+
+      u[n+1] = (c1 * u[n] - c0 * u[n-1]) / c2;
+      d[n+1] = (c1 * d[n] + dc1 * u[n] - c0 * d[n-1]) / c2;
+   }
+}
+
 void Poly_1D::CalcLegendre(const int p, const double x, double *u)
 {
    // use the recursive definition for [-1,1]:

fem/fe/fe_base.hpp

@@ -1112,7 +1112,28 @@ public:
    static void CalcBernstein(const int p, const double x, double *u, double *d)
    { CalcBinomTerms(p, x, 1. - x, u, d); }
 
+   // Evaluate the values of the Jacobi polynomials, P^{(\alpha,\beta)}_n,
+   // on [0,1] with degrees ranging from 0 to p and parameters alpha and beta
+   // at a point x.  Internally the point is transformed to the traditional
+   // domain of [-1,1].
+   static void CalcJacobi(const int p, const double alpha, const double beta,
+                          const double x, double *u);
+
+   // Evaluate the values and derivatives of the Jacobi polynomials,
+   // P^{(\alpha,\beta)}_n, on [0,1] with degrees ranging from 0 to p and
+   // parameters alpha and beta at a point x.  Internally the point is
+   // transformed to the traditional domain of [-1,1].
+   static void CalcJacobi(const int p, const double alpha, const double beta,
+                          const double x, double *u, double *d);
+
+   // Evaluate the values of the Legendre polynomials, P_n, on [0,1] with
+   // degrees ranging from 0 to p at a point x.  Internally the point is
+   // transformed to the traditional domain of [-1,1].
    static void CalcLegendre(const int p, const double x, double *u);
+
+   // Evaluate the values and derivatives of the Legendre polynomials, P_n,
+   // on [0,1] with degrees ranging from 0 to p at a point x.  Internally the
+   // point is transformed to the traditional domain of [-1,1].
    static void CalcLegendre(const int p, const double x, double *u, double *d);
 
    ~Poly_1D();

fem/fe/fe_h1.cpp

@@ -1040,4 +1040,315 @@ void H1_WedgeElement::CalcDShape(const IntegrationPoint &ip,
    }
 }
 
+H1_PyramidElement::H1_PyramidElement(const int p, const int btype)
+   : NodalFiniteElement(3, Geometry::PYRAMID,
+                        (p + 1) * (p + 2) * (2 * p + 3) / 6,
+                        p, FunctionSpace::Qk)
+{
+   const double *cp = poly1d.ClosedPoints(p, VerifyNodal(VerifyClosed(btype)));
+
+   const double **mcp = new const double*[p];
+   mcp[0] = NULL;
+   for (int k=1; k<p; k++)
+   {
+      mcp[k] = poly1d.ClosedPoints(k, VerifyNodal(VerifyClosed(btype)));
+   }
+
+#ifndef MFEM_THREAD_SAFE
+   shape_x.SetSize(p + 1);
+   shape_y.SetSize(p + 1);
+   shape_z.SetSize(p + 1);
+   dshape_x.SetSize(p + 1);
+   dshape_y.SetSize(p + 1);
+   dshape_z.SetSize(p + 1);
+   ddshape_x.SetSize(p + 1);
+   ddshape_y.SetSize(p + 1);
+   ddshape_z.SetSize(p + 1);
+   u.SetSize(dof);
+   du.SetSize(dof, dim);
+   ddu.SetSize(dof, (dim * (dim + 1)) / 2);
+#else
+   Vector shape_x(p + 1), shape_y(p + 1), shape_z(p + 1);
+#endif
+
+   // vertices
+   Nodes.IntPoint(0).Set3(cp[0], cp[0], cp[0]);
+   Nodes.IntPoint(1).Set3(cp[p], cp[0], cp[0]);
+   Nodes.IntPoint(2).Set3(cp[p], cp[p], cp[0]);
+   Nodes.IntPoint(3).Set3(cp[0], cp[p], cp[0]);
+   Nodes.IntPoint(4).Set3(cp[0], cp[0], cp[p]);
+
+   // edges
+   int o = 5;
+   for (int i = 1; i < p; i++)  // (0,1)
+   {
+      Nodes.IntPoint(o++).Set3(cp[i], cp[0], cp[0]);
+   }
+   for (int i = 1; i < p; i++)  // (1,2)
+   {
+      Nodes.IntPoint(o++).Set3(cp[p], cp[i], cp[0]);
+   }
+   for (int i = 1; i < p; i++)  // (3,2)
+   {
+      Nodes.IntPoint(o++).Set3(cp[i], cp[p], cp[0]);
+   }
+   for (int i = 1; i < p; i++)  // (0,3)
+   {
+      Nodes.IntPoint(o++).Set3(cp[0], cp[i], cp[0]);
+   }
+   for (int i = 1; i < p; i++)  // (0,4)
+   {
+      Nodes.IntPoint(o++).Set3(cp[0], cp[0], cp[i]);
+   }
+   for (int i = 1; i < p; i++)  // (1,4)
+   {
+      Nodes.IntPoint(o++).Set3(cp[p-i], cp[0], cp[i]);
+   }
+   for (int i = 1; i < p; i++)  // (2,4)
+   {
+      Nodes.IntPoint(o++).Set3(cp[p-i], cp[p-i], cp[i]);
+   }
+   for (int i = 1; i < p; i++)  // (3,4)
+   {
+      Nodes.IntPoint(o++).Set3(cp[0], cp[p-i], cp[i]);
+   }
+
+   // quadrilateral face
+   for (int j = 1; j < p; j++)
+   {
+      for (int i = 1; i < p; i++)
+      {
+         Nodes.IntPoint(o++).Set3(cp[i], cp[j], cp[0]);
+      }
+   }
+
+   // triangular faces
+   for (int j = 1; j < p; j++)
+      for (int i = 1; i + j < p; i++)  // (0,1,4)
+      {
+         double w = cp[i] + cp[j] + cp[p-i-j];
+         Nodes.IntPoint(o++).Set3(cp[i]/w, cp[0], cp[j]/w);
+      }
+   for (int j = 1; j < p; j++)
+      for (int i = 1; i + j < p; i++)  // (1,2,4)
+      {
+         double w = cp[i] + cp[j] + cp[p-i-j];
+         Nodes.IntPoint(o++).Set3(1.0 - cp[j]/w, cp[i]/w, cp[j]/w);
+      }
+   for (int j = 1; j < p; j++)
+      for (int i = 1; i + j < p; i++)  // (3,4,2)
+      {
+         double w = cp[i] + cp[j] + cp[p-i-j];
+         Nodes.IntPoint(o++).Set3(cp[j]/w, 1.0 - cp[i]/w, cp[i]/w);
+      }
+   for (int j = 1; j < p; j++)
+      for (int i = 1; i + j < p; i++)  // (0,4,3)
+      {
+         double w = cp[i] + cp[j] + cp[p-i-j];
+         Nodes.IntPoint(o++).Set3(cp[0], cp[j]/w, cp[i]/w);
+      }
+   // interior
+   for (int k = 1; k < p - 1; k++)
+   {
+      double wk = 1.0 - cp[k];
+      for (int j = 1; j < p - k; j++)
+      {
+         for (int i = 1; i < p - k; i++)
+         {
+            Nodes.IntPoint(o++).Set3(mcp[p-k][i] * wk,
+                                     mcp[p-k][j] * wk,
+                                     cp[k]);
+         }
+      }
+   }
+
+   MFEM_ASSERT(o == dof,
+               "Number of nodes does not match the "
+               "number of degrees of freedom");
+   DenseMatrix T(dof);
+
+   for (int m = 0; m < dof; m++)
+   {
+      const IntegrationPoint &ip = Nodes.IntPoint(m);
+
+      double oz = 1.0 - ip.z;
+
+      double tol = 1e-6;
+
+      double x = (ip.z < 1.0) ? (ip.x / (1.0 - ip.z)) : 0.0;
+      double y = (ip.z < 1.0) ? (ip.y / (1.0 - ip.z)) : 0.0;
+      double z = ip.z;
+
+      poly1d.CalcLegendre(p, x, shape_x);
+      poly1d.CalcLegendre(p, y, shape_y);
+
+      o = 0;
+      for (int i = 0; i <= p; i++)
+      {
+         for (int j = 0; j <= p; j++)
+         {
+            int maxij = std::max(i, j);
+            poly1d.CalcJacobi(p - maxij, 2.0 * (maxij + 1.0), 0.0, z, shape_z);
+
+            for (int k = 0; k <= p - maxij; k++)
+            {
+               if (oz <= tol)
+               {
+                  if (maxij == 0)
+                  {
+                     T(o++, m) = shape_z(k);
+                  }
+                  else
+                  {
+                     T(o++, m) = 0.0;
+                  }
+               }
+               else
+               {
+                  T(o++, m) = shape_x(i) * shape_y(j) * shape_z(k) *
+                              pow(1.0 - ip.z, maxij);
+               }
+            }
+         }
+      }
+   }
+   Ti.Factor(T);
+
+   delete [] mcp;
+}
+
+void H1_PyramidElement::CalcShape(const IntegrationPoint &ip,
+                                  Vector &shape) const
+{
+   const int p = order;
+
+#ifdef MFEM_THREAD_SAFE
+   Vector shape_x(order+1);
+   Vector shape_y(order+1);
+   Vector shape_z(order+1);
+   Vector u(dof);
+#endif
+
+   double oz = 1.0 - ip.z;
+
+   double tol = 1e-6;
+
+   double x = (ip.z < 1.0) ? (ip.x / (1.0 - ip.z)) : 0.0;
+   double y = (ip.z < 1.0) ? (ip.y / (1.0 - ip.z)) : 0.0;
+   double z = ip.z;
+
+   poly1d.CalcLegendre(p, x, shape_x);
+   poly1d.CalcLegendre(p, y, shape_y);
+
+   for (int o = 0, i = 0; i <= p; i++)
+   {
+      for (int j = 0; j <= p; j++)
+      {
+         int maxij = std::max(i, j);
+         poly1d.CalcJacobi(p - maxij, 2.0 * (maxij + 1.0), 0.0, z, shape_z);
+         for (int k = 0; k <= p - maxij; k++)
+         {
+            if (oz <= tol)
+            {
+               if (maxij == 0)
+               {
+                  u(o++) = shape_z(k);
+               }
+               else
+               {
+                  u(o++) = 0.0;
+               }
+            }
+            else
+            {
+               u(o++) =
+                  shape_x(i) * shape_y(j) * shape_z(k) * pow(1.0 - ip.z, maxij);
+            }
+         }
+      }
+   }
+
+   Ti.Mult(u, shape);
+}
+
+void H1_PyramidElement::CalcDShape(const IntegrationPoint &ip,
+                                   DenseMatrix &dshape) const
+{
+   const int p = order;
+
+#ifdef MFEM_THREAD_SAFE
+   Vector  shape_x(p + 1),  shape_y(p + 1),  shape_z(p + 1);
+   Vector dshape_x(p + 1), dshape_y(p + 1), dshape_z(p + 1);
+   DenseMatrix du(dof, dim);
+#endif
+
+   double oz = 1.0 - ip.z;
+
+   double tol = 1e-6;
+
+   if (oz <= tol)
+   {
+      mfem::out << "H1_PyramidElement::CalcDShape near apex" << std::endl;
+   }
+
+   double x = (ip.z < 1.0) ? (ip.x / (1.0 - ip.z)) : 0.0;
+   double y = (ip.z < 1.0) ? (ip.y / (1.0 - ip.z)) : 0.0;
+   double z = ip.z;
+
+   poly1d.CalcLegendre(p, x, shape_x, dshape_x);
+   poly1d.CalcLegendre(p, y, shape_y, dshape_y);
+
+   for (int o = 0, i = 0; i <= p; i++)
+   {
+      for (int j = 0; j <= p; j++)
+      {
+         int maxij = std::max(i, j);
+         poly1d.CalcJacobi(p - maxij, 2.0 * (maxij + 1.0), 0.0, z,
+                           shape_z, dshape_z);
+         for (int k = 0; k <= p - maxij; k++)
+         {
+            if (oz > tol)
+            {
+               du(o, 0) = dshape_x(i) * shape_y(j) * shape_z(k) *
+                          pow(1.0 - ip.z, maxij - 1);
+               du(o, 1) = shape_x(i) * dshape_y(j) * shape_z(k) *
+                          pow(1.0 - ip.z, maxij - 1);
+               du(o, 2) = (ip.x * dshape_x(i) * shape_y(j) * shape_z(k) +
+                           ip.y * shape_x(i) * dshape_y(j) * shape_z(k) -
+                           shape_x(i) * shape_y(j) * (1.0 - ip.z) *
+                           (shape_z(k) * maxij -
+                            (1.0 - ip.z) * dshape_z(k))
+                          ) * pow(1.0 - ip.z, maxij - 2);
+            }
+            else
+            {
+               // The gradient is multivalued at the apex. The following
+               // values were produced by setting x and y equal to zero
+               // and computing the limit as z->1.
+
+               du(o, 0) = 0.0;
+               du(o, 1) = 0.0;
+               du(o, 2) = 0.0;
+
+               if (i == 0 && j == 0)
+               {
+                  du(o, 2) = k * (k + 1) * (k + 2) * (k + 3) / 6;
+               }
+               else if (i == 0 && j == 1)
+               {
+                  du(o, 1) = -(k + 1) * (k + 2) * (k + 3) * (k + 4) / 24;
+               }
+               else if (i == 1 && j == 0)
+               {
+                  du(o, 0) = -(k + 1) * (k + 2) * (k + 3) * (k + 4) / 24;
+               }
+            }
+            o++;
+         }
+      }
+   }
+
+   Ti.Mult(du, dshape);
+}
+
 }

fem/fe/fe_h1.hpp

@@ -148,6 +148,31 @@ public:
                            DenseMatrix &dshape) const;
 };
 
+class H1_PyramidElement : public NodalFiniteElement
+{
+private:
+#ifndef MFEM_THREAD_SAFE
+   // mutable Vector shape_0, shape_1, shape_2;
+   // mutable Vector dshape_0_0, dshape_1_0, dshape_2_0;
+   // mutable Vector dshape_0_1, dshape_1_1, dshape_2_1;
+   // mutable Vector u;
+   // mutable DenseMatrix du;
+
+   mutable Vector shape_x, shape_y, shape_z;
+   mutable Vector dshape_x, dshape_y, dshape_z, u;
+   mutable Vector ddshape_x, ddshape_y, ddshape_z;
+   mutable DenseMatrix du, ddu;
+#endif
+   DenseMatrixInverse Ti;
+
+public:
+   H1_PyramidElement(const int p,
+                     const int btype = BasisType::GaussLobatto);
+   virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
+   virtual void CalcDShape(const IntegrationPoint &ip,
+                           DenseMatrix &dshape) const;
+};
+
 } // namespace mfem
 
 #endif

fem/fe_coll.cpp

@@ -1749,7 +1749,7 @@ H1_FECollection::H1_FECollection(const int p, const int dim, const int btype)
          H1_dof[Geometry::TETRAHEDRON] = (TriDof*pm3)/3;
          H1_dof[Geometry::CUBE] = QuadDof*pm1;
          H1_dof[Geometry::PRISM] = TriDof*pm1;
-         H1_dof[Geometry::PYRAMID] = 0;
+         H1_dof[Geometry::PYRAMID] = pm1*pm2*(2*p-3)/6;
          if (b_type == BasisType::Positive)
          {
             H1_Elements[Geometry::TETRAHEDRON] = new H1Pos_TetrahedronElement(p);
@@ -1763,7 +1763,7 @@ H1_FECollection::H1_FECollection(const int p, const int dim, const int btype)
             H1_Elements[Geometry::CUBE] = new H1_HexahedronElement(p, btype);
             H1_Elements[Geometry::PRISM] = new H1_WedgeElement(p, btype);
          }
-         H1_Elements[Geometry::PYRAMID] = new LinearPyramidFiniteElement;
+         H1_Elements[Geometry::PYRAMID] = new H1_PyramidElement(p, btype);
 
          const int &TetDof = H1_dof[Geometry::TETRAHEDRON];
          TetDofOrd[0] = new int[24*TetDof];
@@ -1861,16 +1861,7 @@ H1_FECollection::H1_FECollection(const int p, const int dim, const int btype)
 const FiniteElement *
 H1_FECollection::FiniteElementForGeometry(Geometry::Type GeomType) const
 {
-   if (GeomType != Geometry::PYRAMID || this->GetOrder() == 1)
-   {
-      return H1_Elements[GeomType];
-   }
-   else
-   {
-      MFEM_ABORT("H1 Pyramid basis functions are not yet supported "
-                 "for order > 1.");
-      return NULL;
-   }
+   return H1_Elements[GeomType];
 }
 
 const int *H1_FECollection::DofOrderForOrientation(Geometry::Type GeomType,

fem/intrules.cpp

@@ -1675,7 +1675,7 @@ IntegrationRule *IntegrationRules::PyramidIntegrationRule(int Order)
       ipp.x = ipc.x * (1.0 - ipc.z);
       ipp.y = ipc.y * (1.0 - ipc.z);
       ipp.z = ipc.z;
-      ipp.weight = ipc.weight / 3.0;
+      ipp.weight = 4.0 * ipc.weight * pow(1.0 - ipc.z, 2);
    }
    return PyramidIntRules[Order];
 }

tests/unit/fem/test_eigs.cpp

@@ -38,8 +38,9 @@ enum MeshType
    WEDGE2 = 12,
    TETRAHEDRA = 13,
    WEDGE4 = 14,
-   MIXED3D6 = 15,
-   MIXED3D8 = 16
+   PYRAMID6 = 15,
+   MIXED3D6 = 16,
+   MIXED3D8 = 17
 };
 
 Mesh * GetMesh(MeshType type);
@@ -70,6 +71,7 @@ TEST_CASE("Laplacian Eigenvalues",
           mt == MeshType::WEDGE2     ||
           mt == MeshType::TETRAHEDRA ||
           mt == MeshType::WEDGE4     ||
+          mt == MeshType::PYRAMID6   ||
           mt == MeshType::MIXED3D8 )
       {
          mesh->UniformRefinement();
@@ -580,6 +582,40 @@ Mesh * GetMesh(MeshType type)
          v[0] = 3; v[1] = 0; v[2] = 4; v[3] = 8; v[4] = 5; v[5] = 9;
          mesh->AddWedge(v);
          break;
+      case PYRAMID6:
+         mesh = new Mesh(3, 9, 6);
+         c[0] = 0.0; c[1] = 0.0; c[2] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = a_; c[1] = 0.0; c[2] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = a_; c[1] = b_; c[2] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = 0.0; c[1] = b_; c[2] = 0.0;
+         mesh->AddVertex(c);
+         c[0] = 0.0; c[1] = 0.0; c[2] = c_;
+         mesh->AddVertex(c);
+         c[0] = a_; c[1] = 0.0; c[2] = c_;
+         mesh->AddVertex(c);
+         c[0] = a_; c[1] = b_; c[2] = c_;
+         mesh->AddVertex(c);
+         c[0] = 0.0; c[1] = b_; c[2] = c_;
+         mesh->AddVertex(c);
+         c[0] = 0.5 * a_; c[1] = 0.5 * b_; c[2] = 0.5 * c_;
+         mesh->AddVertex(c);
+
+         v[0] = 0; v[1] = 1; v[2] = 2; v[3] = 3; v[4] = 8;
+         mesh->AddPyramid(v);
+         v[0] = 0; v[1] = 4; v[2] = 5; v[3] = 1; v[4] = 8;
+         mesh->AddPyramid(v);
+         v[0] = 1; v[1] = 5; v[2] = 6; v[3] = 2; v[4] = 8;
+	 mesh->AddPyramid(v);
+         v[0] = 2; v[1] = 6; v[2] = 7; v[3] = 3; v[4] = 8;
+	 mesh->AddPyramid(v);
+         v[0] = 3; v[1] = 7; v[2] = 4; v[3] = 0; v[4] = 8;
+	 mesh->AddPyramid(v);
+         v[0] = 7; v[1] = 6; v[2] = 5; v[3] = 4; v[4] = 8;
+	 mesh->AddPyramid(v);
+         break;
       case MIXED3D6:
          mesh = new Mesh(3, 12, 6);
          c[0] = 0.0; c[1] = 0.0; c[2] = 0.0;