diff --git a/fem/fe/fe_nd.cpp b/fem/fe/fe_nd.cpp index 2f8eaa0b32..1527c91672 100644 --- a/fem/fe/fe_nd.cpp +++ b/fem/fe/fe_nd.cpp @@ -2130,7 +2130,7 @@ void ND_FuentesPyramidElement::calcBasis(const int p, } // Family IV - // Re-using mu from Family I + // Re-using mu from Family II dmu = grad_mu0(z); phi_Q(p, mu01(z, xy, 2), mu01(z, xy, 1), phi_Q2_ij); for (int j=2; j<=p; j++) diff --git a/fem/fe/fe_pos.hpp b/fem/fe/fe_pos.hpp index 29e4773359..a06ede84e0 100644 --- a/fem/fe/fe_pos.hpp +++ b/fem/fe/fe_pos.hpp @@ -258,6 +258,13 @@ public: /// Arbitrary order H1 elements in 3D utilizing the Bernstein basis on a pyramid +/// +/// The pyramid affine-related coordinates $\lambda_i$ for $i=1,\ldots,5$ can +/// be used to define a positive H1 basis by noting that $\lambda_i \ge 0$ +/// inside the pyramid for all $i$ and that $\sum_{i=1}^5\lambda_i=1$. This +/// leads to $1 = (\sum_{i=1}^5\lambda_i)^p$. The terms of this product, +/// expanded as a polynomial in the $\lambda_i$, can be used as a Bernstein +/// basis of order $p$ on a pyramid. class H1Pos_PyramidElement : public PositiveFiniteElement, FuentesPyramid { protected: diff --git a/fem/fe/fe_pyramid.hpp b/fem/fe/fe_pyramid.hpp index 1f04e233bd..7ed48898bc 100644 --- a/fem/fe/fe_pyramid.hpp +++ b/fem/fe/fe_pyramid.hpp @@ -28,6 +28,12 @@ namespace mfem shape functions for the exact sequence elements of all shapes" by Federico Fuentes, Brendan Keith, Leszek Demkowicz, and Sriram Nagaraj, see https://doi.org/10.1016/j.camwa.2015.04.027. + + @note Many of the functions below, e.g. lam1, lam2, etc. and related + functions, are singular or multi-valued at the apex of the pyramid. The + values returned near the apex are computed in the limit z->1 using + (x, y, z) = ((1-z)/2, (1-z)/2, z) i.e. along the line from the center + of the base to the apex. */ class FuentesPyramid { @@ -115,6 +121,7 @@ public: static DenseMatrix grad_lam35(real_t x, real_t y, real_t z); static DenseMatrix grad_lam45(real_t x, real_t y, real_t z); + /// Computes $\lambda_i \nabla \lambda_5 - \lambda_5 \nabla \lambda_i$ static Vector lam15_grad_lam15(real_t x, real_t y, real_t z); static Vector lam25_grad_lam25(real_t x, real_t y, real_t z); static Vector lam35_grad_lam35(real_t x, real_t y, real_t z); @@ -135,6 +142,13 @@ public: { return Vector({lam1(x, y, z), lam4(x, y, z), lam5(x, y, z)}); } /// Vector functions related to the normals to the triangular faces + /// + /// Computes + /// $ + /// \lambda_i \nabla\lambda_j \times \nabla \lambda_5 + /// + \lambda_j \nabla\lambda_5 \times \nabla \lambda_i + /// + \lambda_5 \nabla\lambda_i \times \nabla \lambda_j + /// $ static Vector lam125_grad_lam125(real_t x, real_t y, real_t z); static Vector lam235_grad_lam235(real_t x, real_t y, real_t z); static Vector lam345_grad_lam345(real_t x, real_t y, real_t z); @@ -408,11 +422,6 @@ public: void phi_T(int p, Vector nu, DenseMatrix &u) const; void phi_T(int p, Vector nu, const DenseMatrix &grad_nu, DenseMatrix &u, DenseTensor &grad_u) const; - /* - static void E_E(int p, real_t s0, real_t s1, real_t *u, - real_t *duds0, real_t *duds1); - */ - // static void E_E(int p, real_t s0, real_t s1, real_t *u); /** This is a vector-valued function associated with an edge of a pyramid diff --git a/fem/fe/fe_rt.cpp b/fem/fe/fe_rt.cpp index 225deb3e62..2fc7332f23 100644 --- a/fem/fe/fe_rt.cpp +++ b/fem/fe/fe_rt.cpp @@ -1643,7 +1643,6 @@ void RT_FuentesPyramidElement::calcBasis(const int p, E_Q(p, mu01(z, xy, 2), grad_mu01(z, xy, 2), mu01(z, xy, 1), grad_mu01(z, xy, 1), E_ijk, dE_ijk); // Re-using phi_E from Family I - // phi_E(p, mu01(z), grad_mu01(z), phi_k, dphi_k); const real_t muz = mu0(z); const Vector dmuz(grad_mu0(z)); diff --git a/fem/fe_coll.cpp b/fem/fe_coll.cpp index e24b1d1135..5c81fa273d 100644 --- a/fem/fe_coll.cpp +++ b/fem/fe_coll.cpp @@ -1884,7 +1884,6 @@ H1_FECollection::H1_FECollection(const int p, const int dim, const int btype, else { H1_dof[Geometry::PYRAMID] = (p-1)*(p-2)/2; - // H1_dof[Geometry::PYRAMID] = 0; } if (b_type == BasisType::Positive) { @@ -1908,7 +1907,6 @@ H1_FECollection::H1_FECollection(const int p, const int dim, const int btype, H1_Elements[Geometry::PYRAMID] = new H1_FuentesPyramidElement(p, btype); } } - // H1_Elements[Geometry::PYRAMID] = new LinearPyramidFiniteElement; const int &TetDof = H1_dof[Geometry::TETRAHEDRON]; TetDofOrd[0] = (TetDof > 0) ? new int[24*TetDof] : nullptr; @@ -2257,7 +2255,6 @@ L2_FECollection::L2_FECollection(const int p, const int dim, const int btype, L2_Elements[Geometry::PYRAMID] = new L2_FuentesPyramidElement(p, btype); } } - // L2_Elements[Geometry::PYRAMID] = new P0PyrFiniteElement; L2_Elements[Geometry::TETRAHEDRON]->SetMapType(map_type); L2_Elements[Geometry::CUBE]->SetMapType(map_type); @@ -2490,7 +2487,6 @@ RT_FECollection::RT_FECollection(const int order, const int dim, RT_Elements[Geometry::PRISM] = new RT_WedgeElement(p); RT_dof[Geometry::PRISM] = p*pp1*(3*p + 4)/2; - // RT_Elements[Geometry::PYRAMID] = new RT0PyrFiniteElement(false); RT_Elements[Geometry::PYRAMID] = new RT_FuentesPyramidElement(p); RT_dof[Geometry::PYRAMID] = 3*p*pp1*pp1; } @@ -2907,8 +2903,6 @@ ND_FECollection::ND_FECollection(const int p, const int dim, ND_Elements[Geometry::PRISM] = new ND_WedgeElement(p); ND_dof[Geometry::PRISM] = p*pm1*(3*p-4)/2; - // ND_Elements[Geometry::PYRAMID] = new Nedelec1PyrFiniteElement; - // ND_dof[Geometry::PYRAMID] = 0; ND_Elements[Geometry::PYRAMID] = new ND_FuentesPyramidElement(p); ND_dof[Geometry::PYRAMID] = 3*p*pm1*pm1; } diff --git a/tests/unit/fem/test_poly1d.cpp b/tests/unit/fem/test_poly1d.cpp index 932c791645..703afef7e9 100644 --- a/tests/unit/fem/test_poly1d.cpp +++ b/tests/unit/fem/test_poly1d.cpp @@ -13,107 +13,7 @@ #include "catch.hpp" using namespace mfem; -/* -TEST_CASE("1D Jacobi Polynomials (alpha=0, beta=0)","[Poly_1D]") -{ - // For alpha = beta = 0 the Jacobi polynomials should equal the Legendre - // polynomials - int p = 6; - int order = 2 * p + 1; - const IntegrationRule &ir = IntRules.Get(Geometry::SEGMENT, order); - - DenseMatrix J(p+1, ir.GetNPoints()); - DenseMatrix L(p+1, ir.GetNPoints()); - - for (int i=0; i