• 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 pyramidal finite elements currently implemented in Montjoie. The reference element is the symmetric pyramid, look at the figure below for the conventions used

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

// for a new pyramidal finite element (1 : scalar)
class MyFiniteElement : public PyramidReference<1>
{
   public :
   
   // number of dofs on each edge, quadrangle, triangle, pyramid, 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;

};

PyramidGeomReference

The class PyramidGeomReference implements shape functions methods for pyramidal elements. These shape functions are computed as a linear combination of orthogonal polynomials :

where ψ_j form an orthogonal basis of the pyramidal finite element space. We choose the following basis:

where P_i^{α, β} are Jacobi polynomials orthogonal with respect to weight (1-x)^α(1+x)^β. The coefficients c_{i, j, k} are chosen such that the rational functions ψ_{i, j, k} are actually orthonormal :

The coefficients a_{i, j} are computed by inverting a Vandermonde's matrix:

Most of the methods of this class come from the parent class ElementGeomReference. Below, we detail the methods and attributes specific to the class PyramidGeomReference.

GetNbPointsNodalInside	returns the number of nodal points inside the hexahedron.
type_interpolation_nodal	choice of nodal points
GetLegendrePolynomial	returns Legendre polynomials stored in the class
GetCoefLegendre	returns the coefficients used to normalize Legendre polynomialsj
GetEvenJacobiPolynomial	returns the even Jacobi polynomials
GetCoefEvenJacobi	returns the coefficients used to normalize even Jacobi 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

PyramidReference

The class PyramidReference is the base class for all pyramidal finite elements. The methods come from the parent class ElementReference.

PyramidClassical

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

The following orthogonal basis is used :

where are Jacobi polynomials on [-1,1] orthogonal with respect to the weight The basis functions are based on interpolation points x_j and computed by inverting the Vandermonde matrix :

The basis functions are then equal to

The interpolation points x_j are chosen to coincide with Gauss-Lobatto points on the quadrilateral face, and with nodal points of TriangleClassical on triangular faces. The interior points are placed on parallel planes z = cte with scaled Gauss-Lobatto points. The used quadrature formulas are exact for (1-z)² Q_2r+1 over the cube (QUADRATURE_JACOBI2) (the integrals can be transformed into integrals over the unit cube with Duffy transformation). The class PyramidClassical allows numerical computations over hybrid meshes, the corresponding classes for tetrahedra, prisms and hexahedra are TetrahedronClassical, WedgeClassical, HexahedronGauss. The class PyramidClassical is the class used for pyramids if you have TypeElement = TETRAHEDRON_CLASSICAL in the datafile.

PyramidLobatto

This class is the same as PyramidClassical, except that boundary integrals over the quadrilateral base are performed by using Gauss-Lobatto points instead of Gauss points. The class PyramidLobatto allows numerical computations over hybrid meshes, the corresponding classes for tetrahedra, prisms and hexahedra are TetrahedronClassical, WedgeLobatto, HexahedronLobatto. The class PyramidLobatto is the class used for pyramids if you have TypeElement = TETRAHEDRON_LOBATTO or TypeElement = HEXAHEDRON_LOBATTO in the datafile.

PyramidDgOrtho

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 PyramidClassical. The basis functions used for this class are the orthogonal functions

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

PyramidDgLegendre

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 pyramids). As a result, it can be used for a mesh with affine pyramids (the base of each pyramid must be a parallelogram). The orthogonal functions are equal to

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

PyramidHierarchic

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 PyramidClassical. 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_JACOBI2) for volume integrals. The class PyramidHierarchic allows numerical computations over hybrid meshes, the corresponding classes for tetrahedra, prisms and hexahedra are TetrahedronHierarchic, WedgeHierarchic, HexahedronHierarchic. The class PyramidHierarchic is the class used for pyramids if you have TypeElement = TETRAHEDRON_HIERARCHIC in the datafile.

PyramidHcurlFirstFamily

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

If we denote T the transformation from the symmetric cube to the symmetric pyramid

The finite element space on the cube is equal to

This space is generated with nearly orthogonal basis functions ψ_i :

Each degree of freedom is associated with a position x and a tangent t. The Vandermonde matrix is then equal to

The interpolatory basis functions are then equal to

For each triangular face the nodal points and tangents are chosen to coincide with those of the class TriangleHcurlFirstFamily (e.g. Gauss points along edges). Degrees of freedom for the quadrilateral base coincide with those of the class QuadrangleHcurlFirstFamily. For the degrees of freedom inside the element, nodal points are chosen as interior points of the nodal pyramid (PyramidClassical) of order r+1. These interpolatory basis functions have currently the major drawback of a large condition number (due to the loss of orthogonality of functions ψ_i ?) when r is high. Their use for r greater or equal to 6 should be avoided. The class PyramidHcurlFirstFamily allows numerical computations over hybrid meshes, the corresponding classes for tetrahedra, prisms and hexahedra are TetrahedronHcurlFirstFamily, WedgeHcurlFirstFamily, HexahedronHcurlFirstFamily. The class PyramidHcurlFirstFamily is the class used for pyramids if you have TypeElement = TETRAHEDRON_FIRST_FAMILY in the datafile.

PyramidHcurlOptimalFirstFamily

This finite element class is similar to the class PyramidHcurlFirstFamily. The finite element space has the same expression with q going until r+1 instead of r. Interior points along edges are chosen as interior Gauss-Lobatto points. It allows numerical computations over hybrid meshes, the corresponding classes for tetrahedra, prisms and hexahedra are TetrahedronHcurlFirstFamily, WedgeHcurlOptimalFirstFamily, HexahedronHcurlOptimalFirstFamily. The class PyramidHcurlOptimalFirstFamily is the class used for pyramids if you have TypeElement = TETRAHEDRON_OPTIMAL_FIRST_FAMILY in the datafile.

PyramidHcurlHpFirstFamily

This class implements hierarchical basis functions for edge elements of Nedelec's first family. It can be used to discretize H(curl). The finite element space is the same as for PyramidHcurlFirstFamily. 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, prisms and hexahedra are TetrahedronHcurlOptimalHpFirstFamily, WedgeHcurlHpFirstFamily, HexahedronHcurlHpFirstFamily. The class PyramidHcurlHpFirstFamily is the class used for pyramids if you have TypeElement = TETRAHEDRON_HP_FIRST_FAMILY in the datafile.

PyramidHcurlOptimalHpFirstFamily

This class implements hierarchical basis functions for edge elements of Nedelec's first family. It can be used to discretize H(curl). The finite element space is the same as for PyramidHcurlOptimalFirstFamily. This class implements hierarchical basis functions for optimal edge pyramidal edge elements of Nedelec's first family. By using such 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, prisms and hexahedra are TetrahedronHcurlOptimalHpFirstFamily, WedgeHcurlOptimalHpFirstFamily, HexahedronHcurlOptimalHpFirstFamily. The class PyramidHcurlOptimalHpFirstFamily is the class used for pyramids if you have TypeElement = TETRAHEDRON_OPTIMAL_HP_FIRST_FAMILY in the datafile.

PyramidHdivFirstFamily

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

This space is generated with nearly orthogonal basis functions ψ_i :

with the relations

Each degree of freedom is associated with a position x and a tangent t. The Vandermonde matrix is then equal to

The interpolatory basis functions are then equal to

For each triangular face the interpolation points and tangents are chosen to coincide with interior points of TriangleClassical of order r+2. Degrees of freedom for the quadrilateral base are Gauss points. For the degrees of freedom inside the element, interpolation points are chosen as interior points of the nodal pyramid (PyramidClassical) of order r+2. These interpolatory basis functions have currently the major drawback of a large condition number (due to the loss of orthogonality of functions ψ_i ?) when r is high. Their use for r greater or equal to 6 should be avoided. The class PyramidHdivFirstFamily allows numerical computations over hybrid meshes, the corresponding classes for tetrahedra, prisms and hexahedra are TetrahedronHdivFirstFamily, WedgeHdivFirstFamily, HexahedronHdivFirstFamily. The class PyramidHdivFirstFamily is the class used for pyramids if you have TypeElement = TETRAHEDRON_FIRST_FAMILY in the datafile.

PyramidHdivOptimalFirstFamily

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

This space is generated with nearly orthogonal basis functions ψ_i :

with the relations

Each degree of freedom is associated with a position x and a tangent t. The Vandermonde matrix is then equal to

The interpolatory basis functions are then equal to

For each triangular face the interpolation points and tangents are chosen to coincide with interior points of TriangleClassical of order r+2. Degrees of freedom for the quadrilateral base are Gauss points. For the degrees of freedom inside the element, interpolation points are chosen as interior points of the nodal pyramid (PyramidClassical) of order r+2. These interpolatory basis functions have currently the major drawback of a large condition number (due to the loss of orthogonality of functions ψ_i ?) when r is high. Their use for r greater or equal to 6 should be avoided. The class PyramidHdivOptimalFirstFamily allows numerical computations over hybrid meshes, the corresponding classes for tetrahedra, prisms and hexahedra are TetrahedronHdivFirstFamily, WedgeHdivOptimalFirstFamily, HexahedronHdivOptimalFirstFamily. The class PyramidHdivOptimalFirstFamily is the class used for pyramids if you have TypeElement = TETRAHEDRON_OPTIMAL_FIRST_FAMILY in the datafile.

PyramidHdivHpFirstFamily

This class implements hierarchical basis functions for facet elements of Nedelec's first family. It can be used to discretize H(div). The finite element space is the same as for PyramidHdivFirstFamily. 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, prisms and hexahedra are TetrahedronHdivOptimalHpFirstFamily, WedgeHdivHpFirstFamily, HexahedronHdivHpFirstFamily. The class PyramidHdivHpFirstFamily is the class used for pyramids if you have TypeElement = TETRAHEDRON_HP_FIRST_FAMILY in the datafile.

PyramidHdivOptimalHpFirstFamily

This class implements hierarchical basis functions for 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 PyramidHdivOptimalFirstFamily. 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, prisms and hexahedra are TetrahedronHdivOptimalHpFirstFamily, WedgeHdivOptimalHpFirstFamily, HexahedronHdivOptimalHpFirstFamily. The class PyramidHdivOptimalHpFirstFamily is the class used for pyramids if you have TypeElement = TETRAHEDRON_OPTIMAL_HP_FIRST_FAMILY in the datafile.

GetNbPointsNodalInside

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

Example :

int r = 4; // polynomial degree
PyramidGeomReference pyr;
pyr.ConstructFiniteElement(r);

// number of nodal points inside the pyramid ?
int nb_nodes_inside = pyr.GetNbPointsNodalInside(); // should give 5 here

Location :

FiniteElement/Pyramid/PyramidGeomReference.hxx   FiniteElement/Pyramid/PyramidGeomReferenceInline.cxx

type_interpolation_nodal

This public attribute stores the family of nodal points to use when ConstructFiniteElement is called. You can choose between the following values:

• LOBATTO_BASIS : the nodal points will be constructed to coincide with Gauss-Lobatto points on edges and quadrilateral points. The interior points are located on parallel planes.
• REGULAR_BASIS: the nodal are regularly spaced. This choice leads to a large Lebesgue constant, it implies a Runge phenomenon if nodal points are used for interpolation

The default value is LOBATTO_BASIS. The last value REGULAR_BASIS should be used cautiously for moderate polynomial degree.

Example :

int r = 4; // polynomial degree
PyramidGeomReference pyr;
pyr.type_interpolation_nodal = pyr.REGULAR_BASIS;
pyr.ConstructFiniteElement(r);

Location :

FiniteElement/Pyramid/PyramidGeomReference.hxx

GetLegendrePolynomial

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

Example :

int r = 4; // polynomial degree
PyramidGeomReference pyr;
pyr.ConstructFiniteElement(r);

// Legendre polynomials
const Matrix<Real_wp>& pol_leg = pyr.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/Pyramid/PyramidGeomReference.hxx   FiniteElement/Pyramid/PyramidGeomReferenceInline.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
PyramidGeomReference pyr;
pyr.ConstructFiniteElement(r);

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

// and coefficients of normalization
const VectReal_wp& coef_leg = pyr.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/Pyramid/PyramidGeomReference.hxx   FiniteElement/Pyramid/PyramidGeomReferenceInline.cxx

GetEvenJacobiPolynomial

This method returns the jacobi polynomials P_j^{2i+2, 0} stored in the class.

Example :

int r = 4; // polynomial degree
PyramidGeomReference pyr;
pyr.ConstructFiniteElement(r);

// jacobi polynomials for even numbers
const Vector<Matrix<Real_wp> >& jac = pyr.GetEvenJacobiPolynomial();

// coefficients of normalization
const Matrix<Real_wp&>& coef = pyr.GetCoefEvenJacobi();

// evaluating polynomials P_j^{2i+2, 0} for a given i
int i = 2; Real_wp x = 0.4; VectReal_wp Pj;
EvaluateJacobiPolynomial(jac(i), r+1-i, 2.0*x - 1, Pj);

// multiplication by weights
for (int j = 0; j < Pj.GetM(); j++)
  Pj(j) *= coef(i, j);

Location :

FiniteElement/Pyramid/PyramidGeomReference.hxx   FiniteElement/Pyramid/PyramidGeomReferenceInline.cxx

GetCoefEvenJacobi

This method returns the coefficient that normalizes jacobi polynomials P_j^{2i+2, 0} stored in the class :

Example :

int r = 4; // polynomial degree
PyramidGeomReference pyr;
pyr.ConstructFiniteElement(r);

// jacobi polynomials for even numbers
const Vector<Matrix<Real_wp> >& jac = pyr.GetEvenJacobiPolynomial();

// coefficients of normalization
const Matrix<Real_wp&>& coef = pyr.GetCoefEvenJacobi();

// evaluating polynomials P_j^{2i+2, 0} for a given i
int i = 2; Real_wp x = 0.4; VectReal_wp Pj;
EvaluateJacobiPolynomial(jac(i), r+1-i, 2.0*x - 1, Pj);

// multiplication by weights
for (int j = 0; j < Pj.GetM(); j++)
  Pj(j) *= coef(i, j);

Location :

FiniteElement/Pyramid/PyramidGeomReference.hxx   FiniteElement/Pyramid/PyramidGeomReferenceInline.cxx

ConstructRegularPoints

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

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; j++)
    for (int k = 0; k ≤ r-max(i,j); 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;
PyramidGeomReference::ConstructRegularPoints(r, points1d, points2d_tri, NumNodes2D_tri,
                                             points2d_quad, NumNodes2D_quad, points3d, NumNodes3D);

Location :

FiniteElement/Pyramid/PyramidGeomReference.hxx   FiniteElement/Pyramid/PyramidGeomReference.cxx

ConstructLobattoPoints

This method computes the nodal points of the pyramid 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 base.

Example :

int r = 3; // polynomial degree
Matrix<int> NumNodes2D_quad, coor;  
MeshNumbering<Dimension2>::ConstructQuadrilateralNumbering(r, NumNodes2D_quad, coor);
  
// then we compute nodal points for the required order
VectReal_wp points1d; VectR2 points2d_tri, points2d_quad; VectR3 points3d;
PyramidGeomReference::ConstructLobattoPoints(r, points1d, points2d_tri, points2d_quad, NumNodes2D_quad, points3d);

Location :

FiniteElement/Pyramid/PyramidGeomReference.hxx   FiniteElement/Pyramid/PyramidGeomReference.cxx

ComputeValuesPhiOrthoRef

This method evaluates the orthogonal functions ψ_j (as defined in the section devoted to the class PyramidGeomReference).

Example :

int r = 3; // polynomial degree
PyramidGeomReference pyr;
pyr.ConstructFiniteElement(r);

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

Location :

FiniteElement/Pyramid/PyramidGeomReference.hxx   FiniteElement/Pyramid/PyramidGeomReference.cxx

ComputeGradientPhiOrthoRef

This method evaluates the gradient of orthogonal polynomials ψ_j (as defined in the section devoted to the class PyramidGeomReference).

Example :

int r = 3; // polynomial degree
PyramidGeomReference pyr;
pyr.ConstructFiniteElement(r);

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

Location :

FiniteElement/Pyramid/PyramidGeomReference.hxx   FiniteElement/Pyramid/PyramidGeomReference.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;

PyramidGeomReference pyr;
// constructing orthogonal functions
pyr.ComputeOrthogonalFunctions(r);

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

Location :

FiniteElement/Pyramid/PyramidGeomReference.hxx   FiniteElement/Pyramid/PyramidGeomReference.cxx