• Introduction
• Algebra
• Quadrature
• Common functions
• Mesh
• Finite Elements
  • Fi transform
  • Edges
  • Triangles
  • Quadrilaterals
  • Tetrahedra
  • Pyramids
  • Wedges
  • Hexahedra
• Computations
• Boundary conditions
• Post-processing
• Solved equations
• Time schemes

In this page, we detail the different prismatic finite elements currently implemented in Montjoie. The reference element is the unit prism, look at the figure below for the conventions used

If you want to implement a new prismatic element, you have to derive the class WedgeReference. This class is the base class for all prismatic elements implemented in Montjoie. It depends on the type of element (1: scalar finite element, 2: edge elements, 3: facet elements).

// for a new prismatic finite element (1 : scalar)
class MyFiniteElement : public WedgeReference<1>
{
   public :
   
   // number of dofs on each edge, quadrangle, triangle, wedge, etc ?
   void ConstructNumberMap(NumberMap& map, int dg) const;

   // main method used to construct the finite element object
   // r : order of approximation for basis functions
   // rgeom : order of approximation for geometry (shape functions)
   // rquad : order of quadrature
   // if rgeom is equal to 0 (and similarly rquad equal to 0)
   // we are considering that rgeom = r, and rquad = r
   // if type_quad is equal to -1, we use the default quadrature
   // formulas (Gauss-Legendre formulas)
   // rsurf_tri, rsurf_quad : orders for boundary integrals
   // type_surf_tri, type_surf_quad : type of integration for triangles and quadrangles
   void ConstructFiniteElement(int r, int rgeom = 0, int rquad = 0, int type_quad = -1,
                               int rsurf_tri = 0, int rsurf_quad = 0, int type_surf_tri = -1,
                               int type_surf_quad = -1, int gauss_z = -1);

   // computation of phi_i(x) for a given point x
   // the values phi_i(x) have to be placed in the array phi
   void ComputeValuesPhiRef(const R3& x, VectReal_wp& phi) const;

   // computation of  \nabla phi_i(x) for a given point x
   void ComputeGradientPhiRef(const R3& x, VectR3& grad_phi) const;

};

WedgeGeomReference

The class WedgeGeomReference implements shape functions methods for prismatic elements. These shape functions are computed by tensorization of unit triangle with interval [0, 1]:

where are shape functions of the triangle (as defined in TriangleGeomReference) and are 1-D Lagrange interpolation polynomials (as defined in Globatto). Most of the methods of the class WedgeGeomReference come from the parent class ElementGeomReference. Below, we detail the methods and attributes specific to the class WedgeGeomReference.

GetNbPointsNodalInside	returns the number of nodal points inside the hexahedron.
GetNumNodesTri	returns the nodes numbering (index k from coordinates (i, j)).
GetLegendrePolynomial	returns Legendre polynomials stored in the class
GetCoefLegendre	returns the coefficients used to normalize Legendre polynomialsj
ConstructRegularPoints	computes nodal points regularly spaced
ConstructLobattoPoints	computes nodal points with a small Lebesgue constant
ComputeOrthogonalFunctions	constructs orthogonal functions only
ComputeValuesPhiOrthoRef	Evaluates orthonormal polynomials at a given point
ComputeGradientPhiOrthoRef	Evaluates the gradient of orthonormal polynomials at a given point

WedgeReference

The class WedgeReference is the base class for all prismatic finite elements. The methods come from the parent class ElementReference.

WedgeClassical

This class implements a scalar finite element with interpolatory basis functions. It can be used to discretize H¹. The considered finite element space is equal to

The nodal basis functions are obtained by a tensor product between nodal functions of TriangleClassical and Gauss-Lobatto interpolation functions

The used quadrature formulas are also obtained as a tensor product between rules for the unit triangle and for the unit interval (QUADRATURE_GAUSS) for volume integrals. The class WedgeClassical allows numerical computations over hybrid meshes, the corresponding classes for tetrahedra, pyramids and hexahedra are TetrahedronClassical, PyramidClassical, HexahedronGauss. The class WedgeClassical is the class used for wedges if you have TypeElement = TETRAHEDRON_CLASSICAL in the datafile.

WedgeLobatto

This class is the same as WedgeClassical, except that boundary integrals over quadrilateral faces are evaluated with Gauss-Lobatto points. The volume integrals are computed with Gauss-Lobatto points (instead of Gauss-Legendre for WedgeClassical) in z. As a result, there is a partial mass lumping in z, the elementary mass matrix is block diagonal. It allows numerical computations over hybrid meshes, the corresponding classes for tetrahedra, pyramids and hexahedra are TetrahedronClassical, PyramidLobatto, HexahedronLobatto. The class WedgeLobatto is the class used for wedges if you have TypeElement = TETRAHEDRON_LOBATTO or TypeElement = HEXAHEDRON_LOBATTO in the datafile.

WedgeDgClassical

This class implements a scalar finite element with interpolatory basis functions. It can be used to discretize L² (with Discontinuous Galerkin formulation). The finite element space is the same as for WedgeClassical. The basis functions are given as

where ψ_j are Lagrange interpolation polynomials based on Gauss-Legendre points (instead of Gauss-Lobatto points for WedgeClassical). Similarly to WedgeLobatto, t the elementary mass matrix is block diagonal (partial mass lumping). It allows numerical computations over hybrid meshes, the corresponding classes for tetrahedra, pyramids and hexahedra are TetrahedronClassical, PyramidClassical, HexahedronDgGauss. The class WedgeDgClassical is the class used for wedges if you have TypeElement = TETRAHEDRON_DG_CLASSICAL or TypeElement = HEXAHEDRON_DG_GAUSS in the datafile.

WedgeDgOrtho

This class implements a scalar finite element with orthogonal basis functions. It can be used to discretize L² (with discontinuous Galerkin formulation). The finite element space is the same as for WedgeClassical. The basis functions used for this class are the orthogonal polynomials

where P_i^α,β are Jacobi polynomials. Because of tensorized structure, fast computations are obtained for a large order of approximation. The quadrature rules are obtained with Duffy transformation (QUADRATURE_TENSOR). Boundary integrals are evaluated with classical rules (not tensorized) for triangular faces and with Gauss points for quadrilateral faces. This class allows numerical computations over hybrid meshes, the corresponding classes for tetrahedra, pyramids and hexahedra are TetrahedronDgOrtho, PyramidDgOrtho, HexahedronDgGauss. The class WedgeDgOrtho is the class used for wedges if you have TypeElement = TETRAHEDRON_DG_ORTHO in the datafile.

WedgeDgLegendre

This class implements a scalar finite element with orthogonal basis functions. It can be used to discretize L² (with discontinuous Galerkin formulation). The finite element space is here P_r instead of the optimal space (leading to non-consistent approximation for non-affine wedges). As a result, it can be used for a mesh with affine wedges (the quadrilateral faces of wedges must be parallelograms). The basis functions used for this class are the orthogonal functions

It allows numerical computations over hybrid meshes, the corresponding classes for tetrahedra, pyramids and hexahedra are TetrahedronDgOrtho, PyramidDgLegendre, HexahedronDgLegendre. The class WedgeDgLegendre is the class used for wedges if you have TypeElement = TETRAHEDRON_DG_LEGENDRE in the datafile.

WedgeHierarchic

This class implements a scalar finite element with hierarchical basis functions. It can be used to discretize H¹. The finite element space is the same as for WedgeClassical. The basis functions for this class are hierarchical, allowing the use of a variable order in a mesh. The basis functions are equal to :

The tensor structure of this basis is exploited to obtain faster computations for large values of r. Tensorized quadrature rules are therefore used (QUADRATURE_TENSOR) for volume integrals whereas symmetric rules are used for boundary integrals. The class WedgeHierarchic allows numerical computations over hybrid meshes, the corresponding classes for tetrahedra, pyramids and hexahedra are TetrahedronHierarchic, PyramidHierarchic, HexahedronHierarchic. The class WedgeHierarchic is the class used for wedges if you have TypeElement = TETRAHEDRON_HIERARCHIC in the datafile.

WedgeHcurlFirstFamily

This class implements interpolatory basis functions for edge elements of Nedelec's first family. It can be used to discretize H(curl). The finite element space considered here is equal to

where is the finite element space for triangular edge elements (as defined in TriangleHcurlFirstFamily). Interpolatory basis functions are constructed as a tensor product between basis functions of TriangleHcurlFirstFamily and Gauss-Lobatto points, and a tensor product between interpolatory functions of TriangleClassical and Gauss points.

The class WedgeHcurlFirstFamily allows numerical computations over hybrid meshes, the corresponding classes for tetrahedra, pyramids and hexahedra are TetrahedronHcurlFirstFamily, PyramidHcurlFirstFamily, HexahedronHcurlFirstFamily. The class WedgeHcurlFirstFamily is the class used for wedges if you have TypeElement = TETRAHEDRON_FIRST_FAMILY in the datafile.

WedgeHcurlOptimalFirstFamily

This finite element class is similar to the class WedgeHcurlFirstFamily, but the considered finite element space is equal to :

Interpolation points along edges are chosen as interior Gauss-Lobatto points (instead of Gauss points), so that this element is compatible with pyramids, prisms and hexahedra. It allows numerical computations over hybrid meshes, the corresponding classes for tetrahedra, pyramids and hexahedra are TetrahedronHcurlFirstFamily, PyramidHcurlOptimalFirstFamily, HexahedronHcurlOptimalFirstFamily. The class WedgeHcurlOptimalFirstFamily is the class used for wedges if you have TypeElement = TETRAHEDRON_OPTIMAL_FIRST_FAMILY in the datafile.

WedgeHcurlHpFirstFamily

This class implements hierarchical basis functions for prismatic edge elements of Nedelec's first family. It can be used to discretize H(curl). The finite element space is the same as for WedgeHcurlFirstFamily. By using hierarchical functions, a variable order can be set for each element of the mesh. The basis functions are equal to :

This class allows numerical computations over hybrid meshes, the corresponding classes for tetrahedra, pyramids and hexahedra are TetrahedronHcurlOptimalHpFirstFamily, PyramidHcurlHpFirstFamily, HexahedronHcurlHpFirstFamily. The class WedgeHcurlHpFirstFamily is the class used for wedges if you have TypeElement = TETRAHEDRON_HP_FIRST_FAMILY in the datafile.

WedgeHcurlOptimalHpFirstFamily

This class implements hierarchical basis functions for prismatic optimal edge elements of Nedelec's first family. It can be used to discretize H(curl). The finite element space is the same as for WedgeHcurlOptimalFirstFamily. By using hierarchical functions, a variable order can be set for each element of the mesh. The basis functions are equal to :

This class allows numerical computations over hybrid meshes, the corresponding classes for tetrahedra, pyramids and hexahedra are TetrahedronHcurlOptimalHpFirstFamily, PyramidHcurlOptimalHpFirstFamily, HexahedronHcurlOptimalHpFirstFamily. The class PyramidHcurlOptimalHpFirstFamily is the class used for pyramids if you have TypeElement = TETRAHEDRON_OPTIMAL_HP_FIRST_FAMILY in the datafile.

WedgeHdivFirstFamily

This class implements interpolatory basis functions for facet elements of Nedelec's first family. It can be used to discretize H(div). The finite element space considered here is equal to

where is the finite element space for triangular facet elements (as defined in TriangleHdivFirstFamily). Interpolatory basis functions are constructed as a tensor product between basis functions of TriangleHdivFirstFamily and Gauss points, and a tensor product between interpolatory functions of the triangle and Gauss-Lobatto points.

The scalar basis functions of the triangle are based on interior points of TriangleClassical of order r+2. The class WedgeHdivFirstFamily allows numerical computations over hybrid meshes, the corresponding classes for tetrahedra, pyramids and hexahedra are TetrahedronHdivFirstFamily, PyramidHdivFirstFamily, HexahedronHdivFirstFamily. The class WedgeHdivFirstFamily is the class used for wedges if you have TypeElement = TETRAHEDRON_FIRST_FAMILY in the datafile.

WedgeHdivOptimalFirstFamily

This class implements interpolatory basis functions for facet elements of optimal Nedelec's first family. It can be used to discretize H(div). The finite element space considered here is equal to

where is the finite element space for triangular facet elements (as defined in TriangleHdivFirstFamily). Interpolatory basis functions are constructed as a tensor product between basis functions of TriangleHdivFirstFamily and Gauss points, and a tensor product between interpolatory functions of the triangle and Gauss-Lobatto points.

The scalar basis functions of the triangle are based on interior points of TriangleClassical of order r+2. The class WedgeHdivFirstFamily allows numerical computations over hybrid meshes, the corresponding classes for tetrahedra, pyramids and hexahedra are TetrahedronHdivFirstFamily, PyramidHdivOptimalFirstFamily, HexahedronHdivOptimalFirstFamily. The class WedgeHdivOptimalFirstFamily is the class used for wedges if you have TypeElement = TETRAHEDRON_OPTIMAL_FIRST_FAMILY in the datafile.

WedgeHdivHpFirstFamily

This class implements hierarchical basis functions for prismatic facet elements of Nedelec's first family. It can be used to discretize H(div). The finite element space is the same as for WedgeHdivFirstFamily. By using hierarchical functions, a variable order can be set for each element of the mesh. The basis functions are equal to :

This class allows numerical computations over hybrid meshes, the corresponding classes for tetrahedra, pyramids and hexahedra are TetrahedronHdivOptimalHpFirstFamily, PyramidHdivHpFirstFamily, HexahedronHdivHpFirstFamily. The class WedgeHdivHpFirstFamily is the class used for wedges if you have TypeElement = TETRAHEDRON_HP_FIRST_FAMILY in the datafile.

WedgeHdivOptimalHpFirstFamily

This class implements hierarchical basis functions for prismatic facet elements of optimal Nedelec's first family. It can be used to discretize H(div). The finite element space is the same as for WedgeHdivOptimalFirstFamily. By using hierarchical functions, a variable order can be set for each element of the mesh. The basis functions are equal to :

This class allows numerical computations over hybrid meshes, the corresponding classes for tetrahedra, pyramids and hexahedra are TetrahedronHdivOptimalHpFirstFamily, PyramidHdivHpFirstFamily, HexahedronHdivHpFirstFamily. The class WedgeHdivHpFirstFamily is the class used for wedges if you have TypeElement = TETRAHEDRON_HP_FIRST_FAMILY in the datafile.

GetNbPointsNodalInside

This method returns the number of nodal points strictly inside the wedge. It should be equal to (r-1)(r-2)(r-1)/2 where r is the polynomial degree.

Example :

int r = 4; // polynomial degree
WedgeGeomReference wed;
wed.ConstructFiniteElement(r);

// number of nodal points inside the wedge ?
int nb_nodes_inside = wed.GetNbPointsNodalInside(); // should give 9 here

Location :

FiniteElement/Wedge/WedgeGeomReference.hxx   FiniteElement/Wedge/WedgeGeomReferenceInline.cxx

GetNumNodesTri

This method returns the matrix that links the 2-D numbers (i, j) of a shape function to the 1-D numbers k. The number i is the number of the triangular shape function, and j is the number of the 1-D shape function (in z).

Example :

int r = 4; // polynomial degree
WedgeGeomReference wed;
wed.ConstructFiniteElement(r);

// shape functions satisfy the relation phi_k(x, y, z) = phi_i^tri(x, y) phi_j^1D(z)
// phi_i^tri are shape functions on the triangle
// phi_j^1D are shape functions on the unit interval
// if you want to retrieve the number k from numbers (i, j), use GetNumNodesTri 
const Matrix<int>& num = wed.GetNumNodesTri();

int i = 6, j = 2;
int k = num(i, j);

Location :

FiniteElement/Wedge/WedgeGeomReference.hxx   FiniteElement/Wedge/WedgeGeomReferenceInline.cxx

GetLegendrePolynomial

This method returns the Legendre polynomials stored in the class (used for orthogonal polynomials ψ_j). These Legendre polynomials have been constructed with GetJacobiPolynomial. You can call EvaluateJacobiPolynomial to evaluate these polynomials.

Example :

int r = 4; // polynomial degree
WedgeGeomReference wed;
wed.ConstructFiniteElement(r);

// Legendre polynomials
const Matrix<Real_wp>& pol_leg = wed.GetLegendrePolynomial();

// use EvaluateJacobiPolynomial to evaluate Pn(2*x-1)
VectReal_wp Pn; Real_wp x = 0.2; // x in [0, 1]
EvaluateJacobiPolynomial(pol_leg, r+1, 2.0*x - 1.0, Pn);

Location :

FiniteElement/Wedge/WedgeGeomReference.hxx   FiniteElement/Wedge/WedgeGeomReferenceInline.cxx

GetCoefLegendre

This method returns the coefficient used to normalize Legendre polynomials stored in the class (used for orthogonal functions ψ_j). With these coefficients, you have the orthonormality of Legendre polynomials in [0, 1] :

Example :

int r = 4; // polynomial degree
WedgeGeomReference wed;
wed.ConstructFiniteElement(r);

// Legendre polynomials
const Matrix<Real_wp>& pol_leg = wed.GetLegendrePolynomial();

// and coefficients of normalization
const VectReal_wp& coef_leg = wed.GetCoefLegendre();

// use EvaluateJacobiPolynomial to evaluate Pn(2*x-1)
VectReal_wp Pn; Real_wp x = 0.2; // x in [0, 1]
EvaluateJacobiPolynomial(pol_leg, r+1, 2.0*x - 1.0, Pn);

// multiplying by coefficient
for (int i = 0; i ≤ r; i++)
  Pn(i) *= coef_leg(i);

Location :

FiniteElement/Wedge/WedgeGeomReference.hxx   FiniteElement/Wedge/WedgeGeomReferenceInline.cxx

ConstructRegularPoints

This method computes the nodal points of the wedge with a regular spacing. The nodes numbering is given for the triangle, quadrilateral and wedge.

Example :

int r = 3; // polynomial degree
Matrix<int> NumNodes2D_tri(r+1, r+1);
NumNodes2D_tri.Fill(-1);
// filling the node numbers : NumNodes2D_tri(i, j) is the number of the node that will located at x_i, y_j on the triangle
int num = 0;
for (int i = 0; i ≤ r; i++)
  for (int j = 0; j ≤ r-i; j++)
    NumNodes2D_tri(i, j) = num++;

Matrix<int> NumNodes2D_quad(r+1, r+1);
NumNodes2D_quad.Fill(-1);
// filling the node numbers : NumNodes2D_quad(i, j) is the number of the node that will located at x_i, y_j on the square
num = 0;
for (int i = 0; i ≤ r; i++)
  for (int j = 0; j ≤ r; j++)
    NumNodes2D_quad(i, j) = num++;

Array3D<int> NumNodes3D(r+1, r+1, r+1);
NumNodes3D.Fill(-1);
// filling the node numbers : NumNodes3D(i, j, k) is the number of the node that will located at x_i, y_j, z_k
num = 0;
for (int i = 0; i ≤ r; i++)
  for (int j = 0; j ≤ r-i; j++)
    for (int k = 0; k ≤ r; k++)
      NumNodes3D(i, j, k) = num++;

// then we compute nodal points with the provided numbering
VectReal_wp points1d; VectR2 points2d_tri, points2d_quad; VectR3 points3d;
WedgeGeomReference::ConstructRegularPoints(r, points1d, points2d_tri, NumNodes2D_tri,
                                             points2d_quad, NumNodes2D_quad, points3d, NumNodes3D);

Location :

FiniteElement/Wedge/WedgeGeomReference.hxx   FiniteElement/Wedge/WedgeGeomReference.cxx

ConstructLobattoPoints

This method computes the nodal points of the wedge with a small Lebesgue constant. The constructed points are based on Gauss-Lobatto points. The 1-D nodal points will be equal to Gauss-Lobatto points and 2-D nodal points will coincide with Hesthaven's points of the unit triangle and Gauss-Lobatto points for the quadrilateral faces.

Example :

int r = 3; // polynomial degree
Matrix<int> NumNodes2D_quad, NumNodes2D_tri, NumNodes3D, coor;  
MeshNumbering<Dimension2>::ConstructTriangularNumbering(r, NumNodes2D_tri, coor);
MeshNumbering<Dimension2>::ConstructQuadrilateralNumbering(r, NumNodes2D_quad, coor);
MeshNumbering<Dimension3>::ConstructPrismaticNumbering(r, NumNodes3D, coor);

// then we compute nodal points for the required order
VectReal_wp points1d; VectR2 points2d_tri, points2d_quad; VectR3 points3d;
WedgeGeomReference::ConstructLobattoPoints(r, points1d, points2d_tri, NumNodes2D_tri,
                                           points2d_quad, NumNodes2D_quad, points3d, NumNodes3D);

Location :

FiniteElement/Wedge/WedgeGeomReference.hxx   FiniteElement/Wedge/WedgeGeomReference.cxx

ComputeValuesPhiOrthoRef

This method evaluates the orthogonal functions ψ_j. They are obtained as a product of orthogonal functions of the triangle with Legendre polynomials :

where are orthogonal polynomials of the triangle (as defined in TriangleGeomReference), and P_j^{0, 0} are Legendre polynomials.

Example :

int r = 3; // polynomial degree
WedgeGeomReference wed;
wed.ConstructFiniteElement(r);

// then we can evaluate orthogonal polynomials on a point of the unit prism
R3 pt(0.2, 0.4, 0.1); VectReal_wp psi;
wed.ComputeValuesPhiOrthoRef(pt, psi);

Location :

FiniteElement/Wedge/WedgeGeomReference.hxx   FiniteElement/Wedge/WedgeGeomReference.cxx

ComputeGradientPhiOrthoRef

This method evaluates the gradient of orthogonal polynomials ψ_j.

Example :

int r = 3; // polynomial degree
WedgeGeomReference wed;
wed.ConstructFiniteElement(r);

// then we can evaluate the gradient orthogonal polynomials on a point of the unit prism
R3 pt(0.2, 0.4, 0.1); VectR3 grad_psi;
wed.ComputeGradientPhiOrthoRef(pt, grad_psi);

Location :

FiniteElement/Wedge/WedgeGeomReference.hxx   FiniteElement/Wedge/WedgeGeomReference.cxx

ComputeOrthogonalFunctions

This method constructs the orthogonal functions ψ_j. If you call ConstructFiniteElement, this method is alreay called. If you call only ComputeOrthogonalFunctions, the shape functions will not be available.

Example :

int r = 3;

WedgeGeomReference wed;
// constructing orthogonal functions
wed.ComputeOrthogonalFunctions(r);

// evaluating orthogonal functions
R3 pt_loc(0.2, 0.4, 0.1); VectReal_wp psi;
wed.ComputeValuesPhiOrthoRef(pt_loc, psi);

Location :

FiniteElement/Wedge/WedgeGeomReference.hxx   FiniteElement/Wedge/WedgeGeomReference.cxx