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

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

// for a new triangular finite element (1 : scalar)
class MyFiniteElement : public TriangleReference<1>
{
   public :
   
   // number of dofs on each edge, triangle, tetrahedron, 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 (Dunavant's formulas for triangles)
   // redge, type_edge : quadrature formulas for edges
   void ConstructFiniteElement(int r, int rgeom = 0, int rquad = 0, int type_quad = -1,
                               int redge = 0, int type_edge = -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 R2& x, VectReal_wp& phi) const;

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

};

TriangleGeomReference

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

where ψ_j form an orthogonal basis of polynomials of degree r. We choose the Dubiner's basis:

where P_i^{α, β} are Jacobi polynomials orthogonal with respect to weight (1-x)^α(1+x)^β. The coefficients c_{i, j} are chosen such that the polynomials 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 TriangleGeomReference.

GetNbPointsNodalInside	returns the number of nodal points inside the triangle.
type_interpolation_nodal	choice of nodal points
GetLegendrePolynomial	returns Legendre polynomials stored in the class
GetInverseWeightPolynomial	returns the coefficients ci, j applied to obtain orthonormal polynomials ψj
GetCoefficientLegendre	returns the coefficients used to normalize Legendre polynomialsj
GetNumOrtho2D	returns the orthogonal polynomial number k from the tuple (i, j)
GetOddJacobiPolynomial	returns the odd Jacobi polynomials
GetCoefficientOddJacobi	returns the coefficients used to normalize odd Jacobi polynomialsj
GetInverseVandermonde	returns the inverse of the Vandermonde matrix
SetNodalPoints	constructs shapes functions from a set of nodal points
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

Methods and attributes specific to TriangleReference (inherited from ElementReference)

The class TriangleReference is the base class for all triangular finite elements. Below we list methods and attributes specific to this class (other methods come from the parent class ElementReference).

type_interpolation	sets the type of interpolation points
ConstructQuadrature	constructs the quadrature rules for 2-D and boundary integrals.
ConstructFiniteElement1D	constructs the boundary finite element (trace on an edge)

TriangleClassical

This class implements a scalar finite element with interpolatory basis functions. It can be used to discretize H¹. The finite element space considered consists of polynomials of total degree less or equal to r :

Similarly to TriangleGeomReference, the following orthogonal basis is used (Dubiner's basis) :

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 default interpolation points x_j are chosen as points coming from an electrostatic problem as described by Hestaven (1998). These points include Gauss-Lobatto points along each edge, we have displayed them for r=3, 6, and 9 in the picture below :

You can select another set of points by changing the attribute type_interpolation. You can also provide directly interpolation points as illustrated in the example below:

  TriangleClassical tri;
  // you construct the finite element with default interpolation points
  tri.ConstructFiniteElement(4);

  // then you can provide your own set of points
  VectR2 my_points(tri.GetNbDof());
  // first : three vertices of the unit triangle
  my_points(0).Init(0.0, 0.0); my_points(1).Init(0.0, 0.0); my_points(2).Init(0.0, 0.0);
  // then : points on edges of the unit triangle
  my_points(3).Init(0.2, 0.0); my_points(4).Init(0.5, 0.0); my_points(5).Init(0.8, 0.0);
  my_points(6).Init(0.8, 0.2); my_points(7).Init(0.5, 0.5); my_points(8).Init(0.2, 0.8);
  my_points(9).Init(0.0, 0.8); my_points(10).Init(0.0, 0.5); my_points(11).Init(0.0, 0.2);
  // finally : interior points
  my_points(12).Init(0.2, 0.2); my_points(13).Init(0.6, 0.2); my_points(14).Init(0.2, 0.6);

  tri.SetInterpolationPoints(my_points); // basis functions will be recomputed with the points provided

If you call SetInterpolationPoints, the provided points must be ordered as given in the example above. In the class TriangleClassical, the used quadrature formulas are Dunavant's formulas (QUADRATURE_GAUSS of class TriangleQuadrature) for volume integrals and Gauss-Legendre formulas for boundary integrals. The class TriangleClassical allows numerical computations over hybrid meshes, the corresponding class for quadrilaterals is QuadrangleGauss. The class TriangleClassical is the class used for triangles if you have TypeElement = TRIANGLE_CLASSICAL or TypeElement = QUADRANGLE_DG_GAUSS in the datafile.

TriangleLobatto

This class implements a scalar finite element with interpolatory basis functions. It can be used to discretize H¹. The basis functions are the same for TriangleLobatto and TriangleClassical. The only difference between the two classes is that boundary integrals are evaluated with Gauss-Lobatto points with TriangleLobatto. Besides, the class TriangleLobatto allows numerical computations over hybrid meshes, the corresponding class for quadrilaterals is QuadrangleLobatto. The class TriangleLobatto is the class used for triangles if you have TypeElement = TRIANGLE_LOBATTO or TypeElement = QUADRANGLE_LOBATTO in the datafile.

TriangleDgOrtho

This class implements a scalar finite element with orthogonal basis functions. It can be used to discretize L² (with a discontinuous Galerkin formulation). 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 therefore the quadrature rules obtained with Duffy transformation (QUADRATURE_TENSOR). Boundary integrals are evaluated with Gauss-Legendre rules. This class allows numerical computations over hybrid meshes, the corresponding class for quadrilaterals is QuadrangleDgGauss. The class TriangleDgOrtho is the class used for triangles if you have TypeElement = TRIANGLE_DG_ORTHO in the datafile.

TriangleHierarchic

This class implements a scalar finite element with hierarchical basis functions. It can be used to discretize H¹. The basis functions for this class are hierarchical, allowing the use of a variable order in a mesh. Two sets of basis functions are proposed : a set of basis functions with a tensor structure when expressed in the square (TENSOR_BASIS), a set of basis functions such that degrees of freedom associated with edges are "invariant" by rotation (INVARIANT_BASIS). This last choice is used by default because it is the restriction of 3-D hierarchical finite elements on triangular faces. To select the tensor basis, you can proceed as follows:

  TriangleHierarchic tri;

  tri.SetBasisType(tri.TENSOR_BASIS); // the choice is set here

  // then you construct the finite element
  tri.ConstructFiniteElement(4);

The basis functions (default choice : INVARIANT_BASIS) are the following ones :

If you select TENSOR_BASIS (by setting type_basis before calling ConstructFiniteElement), the following basis functions are implemented

The tensor structure is currently not exploited to obtain faster computations. As a result, symmetric quadrature rules (QUADRATURE_GAUSS) are used for evaluating volume integrals. This class allows numerical computations over hybrid meshes, the corresponding class for quadrilaterals is QuadrangleHierarchic. The class TriangleHierarchic is the class used for triangles if you have TypeElement = TRIANGLE_HIERARCHIC in the datafile.

TriangleHcurlFirstFamily

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

where

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

where

Each degree of freedom is associated with a position x and a tangente t, in the figures below, we display the positions and orientations of the degrees of freedom for r=2,4 and 6

The Vandermonde matrix is then equal to

The interpolatory basis functions are then equal to

Along the edges, the nodal points are Gauss points, whereas interior points are the interior Hesthaven's points (of order r+1). You can select regular interpolation points (instead of Hesthaven's points), in that case the basis functions are not computed through a Vandermonde, but directly as a product of monomials (as described in a paper of Graglia). You can also choose interior Gauss-Lobatto points, the choice is made as follows:

  TriangeHcurlFirstFamily tri;

  // you can choose between NODAL_REGULAR, NODAL_GAUSS or NODAL_LOBATTO
  tri.type_nodal_basis = NODAL_REGULAR;

  // then construct the finite element
  tri.ConstructFiniteElement(4);

The quadrature formulas used are symmetric (QUADRATURE_GAUSS). The class TriangleHcurlFirstFamily allows numerical computations over hybrid meshes, the corresponding class for quadrilaterals is QuadrangleHcurlFirstFamily. The class TriangleHcurlFirstFamily is the class used for triangles if you have TypeElement = TRIANGLE_FIRST_FAMILY or TypeElement = QUADRANGLE_FIRST_FAMILY in the datafile.

TriangleHcurlOptimalFirstFamily

This finite element class implements interpolatory basis functions for edge elements of Nedelec's first family. It can be used to discretize H(curl) space. The interpolatory basis functions of this class are constructed with the same technique as for TriangleHcurlFirstFamily, the only difference is that interior Gauss-Lobatto points are used by default for interpolation points along edges instead of Gauss points. The class TriangleHcurlOptimalFirstFamily allows numerical computations over hybrid meshes, the corresponding class for quadrilaterals is QuadrangleHcurlOptimalFirstFamily. The class TriangleHcurlFirstFamily is the class used for triangles if you have TypeElement = TRIANGLE_OPTIMAL_FIRST_FAMILY in the datafile.

TriangleHcurlSecondFamily

This finite element class implements interpolatory basis functions for edge elements of Nedelec's second family. It can be used to discretize H(curl) space. The associated finite element space is . It uses the same interpolatory basis functions as for TriangleClassical, with an orientation for each degree of freedom. This class allows numerical computations over hybrid meshes, the corresponding class for quadrilaterals is QuadrangleHcurlLobatto. It should be noticed that Nedelec's second family for quadrilaterals produces spurious modes, that's why the user is advised to prefer class TriangleHcurlFirstFamily for computation on hybrid meshes. The class TriangleHcurlSecondFamily is the class used for triangles if you have TypeElement = TRIANGLE_SECOND_FAMILY or TypeElement = QUADRANGLE_HCURL_LOBATTO in the data file.

TriangleHcurlOptimalHpFirstFamily

This finite element class implements hierarchical basis functions for edge elements of Nedelec's first family. It can be used to discretize H(curl) space. The finite element space is the same as for TriangleHcurlFirstFamily. The chosen basis functions are "invariant" by rotation and are equal to the restriction of 3-D basis functions for the tetrahedra over triangular faces.

The class TriangleHcurlOptimalHpFirstFamily allows numerical computations over hybrid meshes, the corresponding class for quadrilaterals is QuadrangleHcurlHpFirstFamily or QuadrangleHcurlOptimalHpFirstFamily(for optimal rate of convergence) . The class TriangleHcurlOptimalHpFirstFamily is the class used for triangles if you have TypeElement = TRIANGLE_HP_FIRST_FAMILY or TypeElement = TRIANGLE_OPTIMAL_HP_FIRST_FAMILY in the data file.

TriangleHdivFirstFamily

This finite element class implements interpolatory basis functions for facet elements (also known as Raviart Thomas elements). It can be used to discretize H(div) space. The finite element space considered here is equal to

where

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

where

Each degree of freedom is associated with a position x and a tangente t, in the figures below, we display the positions and orientations of the degrees of freedom for r=2,4 and 6

The Vandermonde matrix is then equal to

The interpolatory basis functions are then equal to

Along the edges, the nodal points are Gauss points, whereas interior points are the interior Hesthaven's points (of order r+1). You can select regular interpolation points (instead of Hesthaven's points) or interior Gauss-Lobatto points, the choice is made as follows:

  TriangeHdivFirstFamily tri;

  // you can choose between NODAL_REGULAR, NODAL_GAUSS or NODAL_LOBATTO
  tri.type_nodal_basis = NODAL_REGULAR;

  // then construct the finite element
  tri.ConstructFiniteElement(4);

The quadrature formulas used are symmetric (QUADRATURE_GAUSS). The class TriangleHdivFirstFamily allows numerical computations over hybrid meshes, the corresponding class for quadrilaterals is QuadrangleHdivFirstFamily. The class TriangleHdivFirstFamily is the class used for triangles if you have TypeElement = TRIANGLE_FIRST_FAMILY or TypeElement = QUADRANGLE_FIRST_FAMILY in the datafile.

TriangleHdivOptimalFirstFamily

This finite element class implements interpolatory basis functions for facet elements. It can be used to discretize H(div) space. The interpolatory basis functions of this class are constructed with the same technique as for TriangleHdivFirstFamily, the only difference is that interior Gauss-Lobatto points are used by default for interpolation points along edges instead of Gauss points. The class TriangleHdivOptimalFirstFamily allows numerical computations over hybrid meshes, the corresponding class for quadrilaterals is QuadrangleHdivOptimalFirstFamily. The class TriangleHdivFirstFamily is the class used for triangles if you have TypeElement = TRIANGLE_OPTIMAL_FIRST_FAMILY in the datafile.

TriangleHdivOptimalHpFirstFamily

This finite element class implements hierarchical basis functions for facet elements of Nedelec's first family. It can be used to discretize H(div) space. The finite element space is the same as for TriangleHdivFirstFamily. The chosen basis functions are given as

The class TriangleHdivOptimalHpFirstFamily allows numerical computations over hybrid meshes, the corresponding class for quadrilaterals is QuadrangleHdivHpFirstFamily or QuadrangleHdivOptimalHpFirstFamily(for optimal rate of convergence) . The class TriangleHdivOptimalHpFirstFamily is the class used for triangles if you have TypeElement = TRIANGLE_HP_FIRST_FAMILY or TypeElement = TRIANGLE_OPTIMAL_HP_FIRST_FAMILY in the data file.

ConstructRegularPoints

This method computes the nodal points of the triangle with a regular spacing. The nodes numbering is given in the last argument.

Example :

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

// then we compute nodal points with the provided numbering
VectReal_wp points1d; VectR2 points2d;
TriangleGeomReference::ConstructRegularPoints(r, points1d, points2d, NumNodes2D);

Location :

FiniteElement/Triangle/TriangleGeomReference.hxx   FiniteElement/Triangle/TriangleGeomReference.cxx

ConstructLobattoPoints

This method computes the nodal points of the triangle with a small Lebesgue constant. Depending on type_interpolation_nodal, the constructed points are Hesthaven's points (based on an anology with electrostatic's equilibrium) or Fekete points. The 1-D nodal points will be equal to Gauss-Lobatto points since the two families of points coincide with Gauss-Lobatto points on the edges of the triangle.

Example :

int r = 3; // polynomial degree

// then we compute nodal points for the required order
VectReal_wp points1d; VectR2 points2d;
TriangleGeomReference::ConstructLobattoPoints(r, points1d, points2d);

Location :

FiniteElement/Triangle/TriangleGeomReference.hxx   FiniteElement/Triangle/TriangleGeomReference.cxx

ComputeValuesPhiOrthoRef

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

Example :

int r = 3; // polynomial degree
TriangleGeomReference tri;
tri.ConstructFiniteElement(r);

// then we can evaluate orthogonal polynomials on a point of the unit triangle
R2 pt(0.2, 0.4); VectReal_wp psi;
tri.ComputeValuesPhiOrthoRef(pt, psi);

Location :

FiniteElement/Triangle/TriangleGeomReference.hxx   FiniteElement/Triangle/TriangleGeomReference.cxx

ComputeGradientPhiOrthoRef

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

Example :

int r = 3; // polynomial degree
TriangleGeomReference tri;
tri.ConstructFiniteElement(r);

// then we can evaluate the gradient orthogonal polynomials on a point of the unit triangle
R2 pt(0.2, 0.4); VectR2 grad_psi;
tri.ComputeGradientPhiOrthoRef(pt, grad_psi);

Location :

FiniteElement/Triangle/TriangleGeomReference.hxx   FiniteElement/Triangle/TriangleGeomReference.cxx

GetNbPointsNodalInside

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

Example :

int r = 4; // polynomial degree
TriangleGeomReference tri;
tri.ConstructFiniteElement(r);

// number of nodal points inside the triangle ?
int nb_nodes_inside = tri.GetNbPointsNodalInside(); // should give 3 here

Location :

FiniteElement/Triangle/TriangleGeomReference.hxx   FiniteElement/Triangle/TriangleGeomReferenceInline.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 correspond to Hesthaven's points based on an analogy with electrostatic's equilibrium
• FEKETE_BASIS : the nodal points will correspond to Fekete points that maximize the determinant of the Vandermonde matrix.
• 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
TriangleGeomReference tri;
tri.type_interpolation_nodal = tri.FEKETE_BASIS; // Fekete points will be considered
tri.ConstructFiniteElement(r);

Location :

FiniteElement/Triangle/TriangleGeomReference.hxx

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
TriangleGeomReference tri;
tri.ConstructFiniteElement(r);

// Legendre polynomials
const Matrix<Real_wp>& pol_leg = tri.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/Triangle/TriangleGeomReference.hxx   FiniteElement/Triangle/TriangleGeomReferenceInline.cxx

GetCoefficientLegendre

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

Example :

int r = 4; // polynomial degree
TriangleGeomReference tri;
tri.ConstructFiniteElement(r);

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

// and coefficients of normalization
const VectReal_wp& coef_leg = tri.GetCoefficientLegendre();

// 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/Triangle/TriangleGeomReference.hxx   FiniteElement/Triangle/TriangleGeomReferenceInline.cxx

GetInverseWeightPolynomial

This method returns the coefficient used to normalize orthogonal polynomials ψ_j stored in the class. With these coefficients, you have the orthonormality of orthogonal polynomials in the unit triangle :

The coefficients c_n are returned by GetInverseWeightPolynomial.

Example :

int r = 4; // polynomial degree
TriangleGeomReference tri;
tri.ConstructFiniteElement(r);

// coefficients of normalization
const VectReal_wp& coef_psi = tri.GetInverseWeightPolynomial();

// ComputeValuesPhiOrthoRef fills c_m psi_m
R2 pt(0.2, 0.1); VectReal_wp psi;
tri.ComputeValuesPhiOrthoRef(pt, psi);

Location :

FiniteElement/Triangle/TriangleGeomReference.hxx   FiniteElement/Triangle/TriangleGeomReferenceInline.cxx

GetNumOrtho2D

This method returns the numbering for orthogonal polynomials.

Example :

int r = 4; // polynomial degree
TriangleGeomReference tri;
tri.ConstructFiniteElement(r);

// coefficients of normalization
const VectReal_wp& coef_psi = tri.GetInverseWeightPolynomial();

// numbers
const Matrix<int>& num = tri.GetNumOrtho2D();

// double loop
for (int i = 0; i ≤ r; i++)
  for (int j = 0; j < r-i; j++)
    {
      // number for orthogonal polynomial
      int k = num(i, j);
    }

Location :

FiniteElement/Triangle/TriangleGeomReference.hxx   FiniteElement/Triangle/TriangleGeomReferenceInline.cxx

GetOddJacobiPolynomial

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

Example :

int r = 4; // polynomial degree
TriangleGeomReference tri;
tri.ConstructFiniteElement(r);

// jacobi polynomials for odd numbers
const Vector<Matrix<Real_wp> >& jac = tri.GetOddJacobiPolynomial();

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

// evaluating polynomials P_j^{2i+1, 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/Triangle/TriangleGeomReference.hxx   FiniteElement/Triangle/TriangleGeomReferenceInline.cxx

GetCoefficientOddJacobi

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

Example :

int r = 4; // polynomial degree
TriangleGeomReference tri;
tri.ConstructFiniteElement(r);

// jacobi polynomials for odd numbers
const Vector<Matrix<Real_wp> >& jac = tri.GetOddJacobiPolynomial();

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

// evaluating polynomials P_j^{2i+1, 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/Triangle/TriangleGeomReference.hxx   FiniteElement/Triangle/TriangleGeomReferenceInline.cxx

GetInverseVandermonde

This method returns the inverse of Vandermonde matrix stored in the class. This matrix is used to compute shape functions from orthogonal polynomials.

Example :

int r = 4; // polynomial degree
TriangleGeomReference tri;
tri.ConstructFiniteElement(r);

// retrieving VDM^{-1}
const Matrix<Real_wp>& invVDM = tri.GetInverseVandermonde();

Location :

FiniteElement/Triangle/TriangleGeomReference.hxx   FiniteElement/Triangle/TriangleGeomReferenceInline.cxx

SetNodalPoints

This method constructs the shape functions from nodal points provided by the user. It can be used in replacement of ConstructFiniteElement if you want to use your own set of points. The number of points should be equal to (r+1)(r+2)/2 with r+1 nodal points on each edge.

Example :

int r = 3;
VectR2 pts((r+1)*(r+2)/2);
pts(0).Init(0.0, 0.0);  pts(1).Init(1.0, 0.0);  pts(2).Init(0.0, 1.0);
// etc

TriangleGeomReference tri;
// constructing orthogonal polynomials
tri.ComputeOrthogonalFunctions(r);
// then you can provide your own set of nodal points
tri.SetNodalPoints(r, pts);

// evaluating shape functions
R2 pt_loc(0.2, 0.4); VectReal_wp phi;
tri.ComputeValuesPhiNodalRef(pt_loc, phi);

Location :

FiniteElement/Triangle/TriangleGeomReference.hxx   FiniteElement/Triangle/TriangleGeomReference.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;

TriangleGeomReference tri;
// constructing orthogonal polynomials
tri.ComputeOrthogonalFunctions(r);

// evaluating orthogonal functions
R2 pt_loc(0.2, 0.4); VectReal_wp psi;
tri.ComputeValuesPhiOrthoRef(pt_loc, psi);

Location :

FiniteElement/Triangle/TriangleGeomReference.hxx   FiniteElement/Triangle/TriangleGeomReference.cxx

type_interpolation

This public attribute stores the family of interpolation points to use when ConstructFiniteElement is called. This attribute is relevant only for interpolatory elements such as TriangleClassical. The list of available points is detailed in the description of type_interpolation_nodal.The default value is LOBATTO_BASIS.

Example :

int r = 4; // polynomial degree
TriangleClassical tri;
tri.type_interpolation = TriangleGeomReference::FEKETE_BASIS; // Fekete points will be considered
tri.ConstructFiniteElement(r);

Location :

FiniteElement/Triangle/TriangleReference.hxx

ConstructQuadrature

This method computes the quadrature rules to use for the interior and for the boundary. This method should not be called in a regular use, since it is already called by ConstructFiniteElement. If you define a new finite element that derives from TriangleReference, you can call this method to construct quadrature rules for the triangle. If r is the order for quadrature rules, the rules should integrate exactly polynomials of degree 2r or 2r-1 such that the elementary matrices are integrated accurately.

Example :

class MyTriangle : public TriangleReference<1>
{
public:
  void ConstructFiniteElement(int r, int rgeom = 0, int rquad = 0, int type_quad = -1,
                                int rsurf = 0, int type_surf = -1)
  {
    // constructing quadrature rules :
    this->ConstructQuadrature(r, type_quad, rsurf, type_surf);
  }

};

Location :

FiniteElement/Triangle/TriangleReference.hxx   FiniteElement/Triangle/TriangleReference.cxx

ConstructFiniteElement1D

This virtual method computes the 1-D finite element which is the trace of the triangular finite element. The default constructed 1-D element is EdgeGaussReference, which corresponds to the trace of TriangleClassical. Since this method is virtual, you can overload it without overloading ConstructFiniteElement.

Example :

class MyTriangle : public TriangleReference<1>
{
public:
  void ConstructFiniteElement1D(int r, int rgeom, int rquad, int type_quad)
  {
    // constructing 1-D finite element, e.g. EdgeLobattoReference
    if (this->element_surface != NULL)
      delete this->element_surface;

    EdgeLobattoReference* edge = new EdgeLobattoReference();
    edge->ConstructFiniteElement(r, rgeom, rquad, type_quad);    
    this->element_surface = edge;    
  }

};

Location :

FiniteElement/Triangle/TriangleReference.hxx   FiniteElement/Triangle/TriangleReference.cxx