Finite elements

The finiteElements library concerns all the stuff related to the elementary objects of a finite elements discretization. When these objects are provided by a mesh tool, they are GeomElement (for geometric element) and when they also carry interpolation information, they are Element. In few words, a Mesh is a collection of GeomElement objects and a finite element Space is a collection of Element objects. In the framework of the finite elements, it is common (and useful) to link any element to a reference element, named in this library: GeomRefElement and RefElement. GeomRefElement is a class whom inherit fundamental geometric elements GeomSegment, GeomTriangle, … and RefElement is a class whom inherit all the finite element types LagrangeSegment, LagrangeTriangle, LagrangeTetrahedron, …

The Interpolation class carries some general information on the finite element such that the finite element type and its conforming space. Besides, this library provides integration method stuff (IntegrationMethod, Quadrature, QuadratureRule).

The Interpolation class

The purpose of the Interpolation class is to store data related to finite element interpolation. It concerns only general description of an interpolation: type of the finite element interpolation (Lagrange, Hermite, …), a subtype among standard, Gauss-Lobatto (only for Lagrange) and first or second family (only for Nedelec), the “order” of the interpolation for Lagrange family and the space conformity (L2, H1, Hdiv, …). For instance, a Lagrange standard P1 with H1 conforming is the classic P1 Lagrange interpolation but a Lagrange standard P1 with L2 conforming means a discontinuous interpolation (like Galerkin discontinuous approximation) where the vertices of an element are not shared. For the moment, this class is only intended to support the following choices:

type	subtype	order
Lagrange	standard, Gauss-Lobatto	any order
Hermite	standard
Morley	standard	order 2
Argyris	standard	order 5
Crouzet-Raviart	standard	order 1
Raviart-Thomas	standard	any order
Nedelec	first or second family	any order
Nedelec face	first or second family	any order
Nedelec edge	first or second family	any order

Hint

Nedelec face and Nedelec edge are specific to 3D. Nedelec face corresponds to the 2D Raviart-Thomas family and Nedelec edge corresponds to the 2D Nedelec family.

Hint

The Interpolation class is intended to support a lot of finite element family, but all are not available. See child classes of RefElement class to know which families are actually available.

This class is mainly used with the GeomElement class (definition of the geometric transformation) and the Space class (definition of the finite element approximation related to the discretized space).

The class Interpolation has the following attributes:

class Interpolation
  {
   private:
    const FEType type_;       //interpolation type
    const FESubType subtype_; //interpolation subtype
    const number_t numtype_;  //additional number type (degree for Lagrange)
    SobolevType conformSpace_;//conforming space
    bool isoparametric_;      //isoparametric flag
    String name_;             //name of the interpolation type
    String subname_;          //interpolation sub_name
    ...
  };

FEType, FESubType and SpaceType are enumerations defined in the config.hpp file as follows:

enum SobolevType{L2=0,H1,Hdiv,Hcurl,Hrot=Hcurl,H2,Hinf};
enum FEType {Lagrange,Hermite,CrouzeixRaviart,Nedelec,RaviartThomas,NedelecFace, NedelecEdge,BuffaChristiansen, BuffaChristiansen,_Morley,_Argyris};
enum FESubType  {standard=0,gaussLobatto,firstFamily,secondFamily};

In addition, there is a static attribute to manage a unique list of instantiated Interpolation objects.

static std::vector<Interpolation*> theInterpolations;

This class proposes some access member functions and functions returning some interpolation properties:

FEType type() const;
FESubType subtype() const;
number_t numtype() const;
SobolevType conformSpace() const;
String name() const;
String subName() const;
String conformSpaceName();
bool isIsoparametric() const;
void isIsoparametric(bool is) ;
bool isContinuous();
bool areDofsAllScalar();
number_t maximumDegree();

some error member functions:

void badType() const;
void badSubType() const;
void badSpace() const;
void badType(const number_t) const;
void badSubType(const number_t) const;
void badDegree(const number_t) const;

There is only one explicit constructor (no copy constructor):

Interpolation(FEType, FESubType, number_t, SobolevType);

that initializes the attributes and add the created object in the static list theInterpolations.

The process creating an Interpolation object is governed by other classes (in particular the Space class) using the function:

Interpolation* findInterpolation(FEType, FESubType, number_t, SpaceType);
Interpolation& interpolation(FEType, FESubType, number_t, SobolevType);

These functions try to find an existing Interpolation object corresponding to the given parameters, if it does not exist one, a new object is created by calling the constructor.

It is possible to work with a collection of Interpolation objects that are handles as a vector by the Interpolations class (see the definition of PCollection). It may be used as following:

Interpolations ints(3);
for(Number i=1;i<=3;i++) ints(i)=interpolation(Lagrange,_standard,i,H1);

See also

• library= finiteElements,

• header= Interpolation.hpp,

• implementation=Interpolation.cpp,

• test=test_Interpolation.cpp,

• header dep= utils.h, config.h

The RefDof class

The RefDof class is a class representing a generalized degree of freedom in a reference element.

It is defined by:

• A support that can be a point, a side, a side of side or the whole element. When it is a point, we also have its coordinates.

• A dimension, that means the number of components of shape functions of the DoF.

• Projections data (type and direction) when the DoF is a projection DoF. The projection type is 0 for a general projection, 1 for a \(u\dot n\) DoF, 2 for a \(u\times n\) DoF.

• Derivative data (order and direction) when the DoF is a derivative DoF. The order is 0 when the DoF is defined by an ‘integral’ over its support, or \(n\) for a Hermite DoF derivative of order \(n>0\) or a moment DoF of order \(n>0\).

• A shareable property: A DoF is shareable when it is shared between adjacent elements according to space conformity. Such a DoF can be shared or not according to interpolation. For example, every DoFs of a Lagrange H1-confirming element on a side is shared, while DoFs in discontinuous interpolation are not.

class RefDof
{
 private:
  bool sharable_;              //D.o.F is shared according to space conformity
  DofLocalization where_;      //hierarchic localization of DoF
  number_t supportNum_;        //support localization
  number_t index_;             //rank of the DoF in its localization
  dimen_t supportDim_;         //dimension of the D.o.F geometric support
  number_t nodeNum_;           //node number when a point DoF
  dimen_t dim_;                //number of components of shape functions of D.o.F
  std::vector<real_t> coords_; //coordinates of D.o.F support if supportDim_=0
  number_t order_;             //order of a derivative D.o.F or a moment D.o.F
  std::vector<real_t> derivativeVector_; //direction vector(s) of a derivative
  ProjectionType projectionType_;        //type of projection
  std::vector<real_t> projectionVector_; //direction vector(s) of a projection
  DiffOpType diffop_;                    //DoF differential operator (_id,_dx,_dy, ...)
  string_t name_;                        //D.o.F-type name for documentation
}

This class proposes mainly a full constructor and accessors to member data.

Hint

For non-nodal DoFs (supportDim_=0), virtual coordinates are also defined. They are useful when dealing with essential conditions.

Hint

Up to now, derivative vectors are not used and the normal vector involved in \(n.\nabla\) DoFs is stored in the projection vector.

In order to sort DoFs efficiently, the following class is offered:

class DofKey
{
 public:
  DofLocalization where;
  number_t v1, v2;           //used for vertex, edge, element DoFs
  std::vector<number_t> vs;  //used for face DoFs
  number_t locIndex;         //local DoF index
}

associated to the comparison operator:

bool operator<(const DofKey& k1, const DofKey& k2);

See also

• library= finiteElements,

• header= RefDof.hpp,

• implementation=RefDof.cpp,

• test=test_segment.cpp,test_triangle.cpp, test_quadrangle.cpp, test_tetrahedron.cpp, test_hexahedron.cpp, test_prism.cpp.

• header dep= utils.h, config.h

Geometric data of a reference element

The GeomRefElement class

The GeomRefElement object carries geometric data of a reference element: type, dimension, number of vertices, sides and side of sides, measure, and coordinates of vertices. Furthermore, to help to find geometric data on sides, it also carries type, vertex numbers of each side on the first hand, and vertex numbers and relative number to parent side of each side of side on the second hand.

class GeomRefElement
{
 protected:
  ShapeType shapeType_;     //element shape
  const Dimen dim_;         //element dimension
  const Number nbVertices_, nbSides_, nbSideOfSides_; //number of vertices, ...
  const Real measure_;                       //length or area or volume
  vector<Real> centroid_;                    //coordinates of centroid
  vector<Real> vertices_;                    //coordinates of vertices
  vector<ShapeType> sideShapeTypes_;         //shape type of each side
  vector<vector<Number> > sideVertexNumbers_;
  vector<vector<Number> > sideOfSideVertexNumbers_;
  vector<vector<Number> > sideOfSideNumbers_;

 public:
  static vector<GeomRefElement*> theGeomRefElements; //vector carrying all GeomReferenceElements
};

It offers:

• few constructors (the default one, a general one with dimension, measure, centroid, and number of vertices, edges and faces and one more specific for each dimension):

  GeomRefElement(); //default constructor
  GeomRefElement(ShapeType, const Dimen d, const Real m, const Real c, const Number v, const Number e, const Number f);        //general
  GeomRefElement(ShapeType, const Real m = 1.0, const Real c = 0.5);     //1D
  GeomRefElement(ShapeType, const Real m, const Real c, const Number v); //2D
  GeomRefElement(ShapeType, const Real m, const Real c, const Number v, const Number e);                                        //3D
  

• a private copy constructor and a private assign operator,

• some public access functions to data:

  Dimen dim() const;
  Number nbSides() const;
  Number nbSideOfSides() const;
  Real measure() const;
  vector<Real>::const_iterator centroid() const;
  vector<Real>::const_iterator vertices() const;
  vector<ShapeType> sideShapeTypes() const;
  const std::vector<std::vector<Number> >& sideVertexNumbers() const;
  const std::vector<std::vector<Number> >& sideOfSideVertexNumbers() const;
  

• some public methods to get data on sides and side of sides:

  String shape(const Number sideNum = 0) const; //element (side) shape as a string
  ShapeType shapeType(const Number sideNum = 0) const; /element (side) shape as a ShapeType
  bool isSimplex() const;                              //true if a simplex
  Number nbVertices(const Number sideNum = 0) const;   //number of element (side) vertices
  Number sideVertex(const Number id, const Number sideNum = 0) const;//number of element (side)
  std::vector<Real>::const_iterator vertex(const Number vNum) const; //coordinates of n-th vertex
  Number sideVertexNumber(const Number vNum, const Number sideNum) const;
  Number sideOfSideVertexNumber(const Number vNum, const Number sideOfSideNum) const;
  Number sideOfSideNumber(const Number i, const Number sideNum) const; //local number of i-th edge
  void rotateVertices(const number_t newFirst,vector<number_t>&) const;// circular permutation
  vector<real_t> sideToElt(number_t,vector<real_t>::const_iterator) const;
  virtual Real measure(const Dimen dim, const Number sideNum = 0) const = 0;
  // projection of a point onto ref element
  virtual std::vector<real_t> projection(const vector<real_t>&, real_t&) const;
  // test if a point belongs to current element
  virtual bool contains(std::vector<real_t>& p, real_t tol= theTolerance) const;
  

• some protected append methods for vertices:

  //build vertex 1d/2d/3d coordinates
  void vertex(vector<Real>::iterator& it, const Real x1);
  void vertex(vector<Real>::iterator& it, const Real x1, const Real x2);
  void vertex(vector<Real>::iterator& it, const Real x1, const Real x2, const Real x3);
  

• some public error handlers:

  void noSuchSideOfSide(const Number) const;
  void noSuchSide(const Number) const;
  void noSuchSide(const Number, const Number, const Number = 0, const Number = 0) const;
  void noSuchFunction(const String& s) const;
  

• an external search function in the static run-time list of GeomRefElement:

  GeomRefElement* findGeomRefElement(ShapeType);
  

See also

• library= finiteElements,

• header= GeomElement.hpp,

• implementation=GeomElement.cpp,

• test=test_segment.cpp, test_triangle.cpp, test_quadrangle.cpp, test_tetrahedron.cpp, test_hexahedron.cpp, test_prism.cpp.

• header dep= utils.h, GeomElement.hpp

Child classes of GeomRefElement

The child classes of GeomRefElement are GeomRefSegment in 1D, GeomRefTriangle and GeomRefQuadrangle in 2D and GeomRefHexahedron, GeomRefTetrahedron, GeomRefPrism and GeomRefPyramid in 3D.

Except implementation of virtual functions, these child classes provide 2 new member functions used in class constructors:

void sideNumbering();       //local numbers of vertices on sides
void sideOfSideNumbering(); //local numbers of vertices on side of sides

The last one is not defined for GeomRefSegment.

See also

• library= finiteElements,

• header= GeomRefSegment.hpp,

• implementation=GeomRefSegment.cpp,

• test=test_segment.cpp, .

• header dep= config.h, GeomElement.hpp

See also

• library= finiteElements,

• header= GeomRefTriangle.hpp,

• implementation=GeomRefTriangle.cpp,

• test=test_triangle.cpp, .

• header dep= config.h, GeomElement.hpp

See also

• library= finiteElements,

• header= GeomRefQuadrangle.hpp,

• implementation=GeomRefQuadrangle.cpp,

• test=test_quadrangle.cpp, .

• header dep= config.h, GeomElement.hpp

See also

• library= finiteElements,

• header= GeomRefTetrahedron.hpp,

• implementation=GeomRefTetrahedron.cpp,

• test=test_tetrahedron.cpp, .

• header dep= config.h, GeomElement.hpp

See also

• library= finiteElements,

• header= GeomRefHexahedron.hpp,

• implementation=GeomRefHexahedron.cpp,

• test=test_hexahedron.cpp, .

• header dep= config.h, GeomElement.hpp

See also

• library= finiteElements,

• header= GeomRefPrism.hpp,

• implementation=GeomRefPrism.cpp,

• test=test_prism.cpp, .

• header dep= config.h, GeomElement.hpp

See also

• library= finiteElements,

• header= GeomRefPyramid.hpp,

• implementation=GeomRefPyramid.cpp,

• test=test_pyramid.cpp, .

• header dep= config.h, GeomElement.hpp

Reference elements

The ShapeValues class

A ShapeValues object stores the evaluation of a shape function at a given point. Thus, it carries 2 real vectors of:

class ShapeValues
{
 public:
  vector<Real> w;          //shape functions at a  point
  vector<vector<Real>> dw; //first derivatives (dx,dy,[dz])
  vector<vector<Real>> d2w;//2nd derivatives   (dxx,dyy,dxy,[dzz,dxz, dyz])
 ...
};

It offers:

• basic constructors:

  ShapeValues();      //!< default constructor
  ShapeValues(const ShapeValues&);                     //copy constructor
  ShapeValues& operator=(const ShapeValues&);  // assignment operator=
  ShapeValues(const RefElement&);              //constructor with associated RefElement
  ShapeValues(const RefElement&, const Dimen); //constructor with associated RefElement
  

• two public size accessors and empty query:

  Dimen dim() const;
  Number nbDofs() const;
  bool isEmpty() const;
  

• some functions addressing the data:

  void resize(const RefElement&, const Dimen); // resize
  void set(const real_t);            // set shape functions to const value
  void assign(const ShapeValues&);   // fast assign of shapevalues into current
  

• some mapping functions

  void extendToVector(dimen_t d);  //extend scalar shape function to vector
  void map(const ShapeValues&, const GeomMapData&, bool der1, bool der2);  //standard mapping
  void contravariantPiolaMap(const ShapeValues&, const GeomMapData&, bool der1, bool der2);
  void covariantPiolaMap(const ShapeValues&, const GeomMapData&, bool der1, bool der2);
  void Morley2dMap(const ShapeValues&, const GeomMapData&, bool der1, bool der2);
  void Argyris2dMap(const ShapeValues& rsv, const GeomMapData& gd, bool der1, bool der2);
  void changeSign(const std::vector<real_t>&, dimen_t);  //change sign of shape functions according to a sign vector
  

The contravariant and covariant Piola maps are respectively used for Hdiv and Hcurl finite element families. They preserve respectively the normal and tangential traces. If \(J\) denotes the jacobian of the mapping from the reference element to any element, the covariant map is \(J^{-t}\mathbf{\hat{u}}\) and the contravariant map \(J/|J|\mathbf{\hat{u}}\) where \(\mathbf{\hat{u}}\) is a vector field in the reference space. The Morley and Argyris map are specific to these elements, they preserve first and second derivatives involved in such elements.

See also

• library= finiteElements,

• header= ShapeValues.hpp,

• implementation=ShapeValues.cpp,

• test=test_pyramid.cpp, .

• header dep= config.h, RefElement.hpp, utils.hpp

The RefElement class

The RefElement is the main class of the finiteElement library. First, it has an object of each type previously seen in this chapter, namely Interpolation, GeomRefElement, ShapeValues and RefDof. Next, It is the mother class of the wide range of true reference elements, identified by shape and interpolation type. That is why it is a complex class collecting a lot of information and providing many numbering functions such as the number of DoFs, the same number on vertices, on sides, on side of sides, …

To define completely a reference element, we also need data on sides and on side of sides, such as DoF numbers and reference elements. Side and side of side reference elements will be built only when they are needed, equivalent to say only they have meaning. It is the case for H1 finite elements, but not for Hcurl and Hdiv finite elements.

So, the RefElement attributes are:

class RefElement
{
 public:
  GeomRefElement* geomRefElem_p;        //pointer to geometric reference element
  const Interpolation* interpolation_p; //interpolation parameters
  vector<RefDof*> refDofs;          //local reference Degrees of Freedom
  FEMapType mapType;                //type of map applied to shape functions of reference element
  DofCompatibility DoFCompatibility;//compatibility rule to applied to side DoFs
  dimen_t dimShapeFunction;         //dimension of shape function
  bool rotateDof;                   //if true, apply rotation (given by children) to shape values
  number_t maxDegree;               //maximum degree of shape functions

 protected:
  String name_;         //reference element name (for documentation)
  Number nbDofs_; //number of Degrees Of Freedom
  Number nbPts_;  //number of Points (for local coordinates of points)
  Number nbDofsOnVertices_;   //nb of sharable D.o.F on vertices (first D.o.F)
  Number nbDofsInSideOfSides_;//nb of sharable D.o.F on side of sides (edges)
  Number nbDofsInSides_;      //nb of sharable D.o.F on sides (faces or edges)
  Number nbInternalDofs_;     //nb of non-sharable D.o.F's
  vector<RefElement*> sideRefElems_;       //pointers to side reference elements
  vector<RefElement*> sideOfSideRefElems_; //pointers to side of side reference elements

 public:
  vector<vector<Number>> sideDofNumbers_;       //DoF numbers for each side
  vector<vector<Number>> sideOfSideDofNumbers_; //DoF numbers for each side of side
  PolynomialsBasis Wk;                     //shape functions basis as polynomials
  vector<PolynomialsBasis> dWk;            //derivatives of shape functions basis as polynomials
  bool hasShapeValues;                     //ref element has shapevalues, generally true
  map<Quadrature*,vector<ShapeValues>> qshvs_;   //temporary structure to store shape values
  map<Quadrature*,vector<ShapeValues>> qshvs_aux;//other temporary structure to store shape values

  static vector<RefElement*> theRefElements; //to store run-time RefElement pointers
  static const bool isSharable_ = true;
};

The following enumeration collects the types of affine mapping from reference element to any element. It depends on finite element type:

enum FEMapType {_standardMap,_contravariantPiolaMap,_covariantPiolaMap,_MorleyMap,_ArgyrisMap};

It offers:

• Some constructors, the default one and a constructor by shape type and interpolation:

  RefElement(); //default constructor
  RefElement(ShapeType, const Interpolation* ); //constructor by shape & interpolation
  

  Note that the copy constructor and the assignment operator are private.

• Some public accessors:

  const Interpolation& interpolation() const;
  String name() const;
  Number nbPts() const;
  Number nbDofsOnVertices() const;    //nb of D.o.F supported by vertices
  Number nbDofsInSides() const;       //nb of D.o.F strictly supported by side
  Number nbDofsInSideOfSides() const; //nb of D.o.F strictly supported by side of side
  bool isSimplex() const ;            //true if shape is a simplex
  dimen_t dim() const;                //element dimension
  
  const RefElement* refElement(Number sideNum = 0) const; //reference element of side
  const GeomRefElement* geomRefElement(Number sideNum = 0) const; //GeomRefElement of side
  Number nbDofs(const Number sideNum = 0, const Dimen sideDim = 0) const;//number of  D.o.F by side
  Number nbInternalDofs(const Number sideNum = 0, const Dimen sideDim = 0) const;
  ShapeType shapeType() const;   //shape of element
  

• Some build functions of side and side of side reference elements:

  void sideOfSideRefElement();      //reference element constructors on element edges
  void sideRefElement();            //reference element constructors on element faces
  

• Some functions related to shape function computation

  void buildPolynomialTree();       //build tree representation of polynomial shape functions
  number_t shapeValueSize() const;
  virtual void computeShapeFunctions()  //compute shape functions
  virtual void computeShapeValues(vector<Real>::const_iterator it_pt, bool der1 = true, bool der2=false)  const = 0; //compute shape functions at point
  

  Shape function of high order element are boring to find. So, XLiFE++ provides tools to get them in a formal way, using Polynomials class. The virtual function computeShapeFunctions() computes them and the buildPolynomialTree() function builds a tree representation of the formal representation of shape functions, in order to make faster their computations.

• Some virtual functions to guarantee correct matching of DoFs:

  virtual vector<number_t> DoFsMap(const number_t& i, const number_t& j, const number_t& k=0) const;
  virtual number_t sideDofsMap(const number_t& n, const number_t& i, const number_t& j, const number_ & k=0) const;
  virtual number_t sideofsideDofsMap(const number_t& n, const number_t& i, const number_t& j=0)const;
  number_t sideOf(number_t) const;       //side number of a DoF
  number_t sideOfSideOf(number_t) const; //return side of side number of a DoF
  

• Some general build functions for DoFs:

  void LagrangeRefDofs(const int, const int, const int, const Dimen);
  

  As the inheritance diagram is based on shape, these functions specify the interpolation type

• For output purpose, a RefElement may be split in first order elements, either simplices or of same shape, at every order:

  virtual vector<vector<number_t>> splitP1() const;
  virtual vector<pair<ShapeType,vector<number_t>>> splitO1() const;
  

• Some external functions to find a reference element in the static list of RefElement and create it if not found:

  RefElement* findRefElement(ShapeType, const Interpolation*);
  RefElement* selectRefSegment(const Interpolation*);
  RefElement* selectRefTriangle(const Interpolation*);
  RefElement* selectRefQuadrangle(const Interpolation*);
  RefElement* selectRefTetrahedron(const Interpolation*);
  RefElement* selectRefPrism(const Interpolation*);
  RefElement* selectRefHexahedron(const Interpolation*);
  RefElement* selectRefPyramid(const Interpolation*);
  

See also

• library= finiteElements,

• header= RefElement.hpp,

• implementation=RefElement.cpp,

• test=test_segment.cpp, test_triangle.cpp, test_quadrangle.cpp, test_tetrahedron.cpp, test_hexahedron.cpp, test_prism.cpp, test_pyramid.cpp,

• header dep= config.h, utils.h, Interpolation.hpp, ShapeValues.hpp, GeomRefElement.hpp, RefDof.hpp.

Child classes for segment

The RefSegment class has three main child classes :

• LagrangeStdSegment which implements standard Lagrange element at any order

• LagrangeGLSegment which implements Lagrange element based on Gauss-Lobatto points at any order

• HermiteStdSegment which implements standard Hermite element, currently only P3 Hermite

Child classes for triangle

The RefTriangle class provides the following types of element:

• standard Lagrange element at any order LagrangeStdTriangle and LagrangeStdTrianglePk; the template LagrangeStdTriangle class is defined only for low order (up to 6) with better performance and LagrangeStdTrianglePk class is designed for any order. It uses a general representation of shape functions on xy-monoms basis get from the solving of a linear system. This method is not stable for very high order.

• HermiteStdTriangle which implements standard Hermite element, currently only P3 Hermite (not H2 conforming)

• CrouzeixRaviartStdTriangle which implements the Crouzeix-Raviart element, currently only P1 (not H1 conforming)

• RaviartThomasTriangle which implements the Raviart-Thomas elements (Hdiv conforming): :cpp:class:` RaviartThomasStdTriangleP1` for Raviart-Thomas standard P1 element and RaviartThomasStdTrianglePk for Raviart-Thomas at any order.

• NedelecTriangle which implements the Nedelec elements (Hcurl conforming): NedelecFirstTriangleP1 for Nedelec first family of order 1 element, NedelecFirstTrianglePk for Nedelec first family of any order and NedelecSecondTrianglePk for Nedelec second family of any order.

• MorleyTriangle and ArgyrisTriangle elements designed for H2-approximation (non H2 conforming and H2 conforming).

Child classes for quadrangle

Child classes for tetrahedron

The RefTetrahedron class provides the following types of element:

• LagrangeStdTetrahedron class is defined only for low order (up to 6) with better performance and LagrangeStdTetrahedronPk class is designed for any order. It uses a general representation of shape functions on xyz-monoms basis get from the solving of a linear system. This method is not stable for very high order.

• CrouzeixRaviartTetrahedron which implements the Crouzeix-Raviart element, currently only P1 (not H1 conforming)

• NedelecFaceTetrahedron which implements the Nedelec elements Hdiv conforming: NedelecFaceFirstTetrahedronPk for Nedelec first family of any order and NedelecFaceSecondTetrahedronPk for Nedelec second family of any order. The second family is not yet available.

• NedelecEdgeTetrahedron which implements the Nedelec elements Hrot conforming: NedelecEdgeFirstTetrahedronPk for Nedelec first family of any order and NedelecEdgeSecondTetrahedronPk for Nedelec second family of any order. The second family is not yet available.

Short identifiers of face/edge elements are defined in the following enumerations:

enum FeFaceType
{
 _RT_1 = 1, RT_1 = _RT_1, _NF1_1 = _RT_1, NF1_1 = _RT_1,
 _RT_2 ,    RT_2 = _RT_2, _NF1_2 = _RT_2, NF1_2 = _RT_2,
 _RT_3 ,    RT_3 = _RT_3, _NF1_3 = _RT_3, NF1_3 = _RT_3,
 _RT_4 ,    RT_4 = _RT_4, _NF1_4 = _RT_4, NF1_4 = _RT_4,
 _RT_5 ,    RT_5 = _RT_5, _NF1_5 = _RT_5, NF1_5 = _RT_5,
 _BDM_1 ,    BDM_1 = _BDM_1, _NF2_1 = _BDM_1, NF2_1 = _BDM_1,
 _BDM_2 ,    BDM_2 = _BDM_2, _NF2_2 = _BDM_2, NF2_2 = _BDM_2,
 _BDM_3 ,    BDM_3 = _BDM_3, _NF2_3 = _BDM_3, NF2_3 = _BDM_3,
 _BDM_4 ,    BDM_4 = _BDM_4, _NF2_4 = _BDM_4, NF2_4 = _BDM_4,
 _BDM_5 ,    BDM_5 = _BDM_5, _NF2_5 = _BDM_5, NF2_5 = _BDM_5
};

enum FeEdgeType
{
 _N1_1 = 1, N1_1 = _N1_1, _NE1_1 = _N1_1, NE1_1 = _N1_1,
 _N1_2 ,    N1_2 = _N1_2, _NE1_2 = _N1_2, NE1_2 = _N1_2,
 _N1_3 ,    N1_3 = _N1_3, _NE1_3 = _N1_3, NE1_3 = _N1_3,
 _N1_4 ,    N1_4 = _N1_4, _NE1_4 = _N1_4, NE1_4 = _N1_4,
 _N1_5 ,    N1_5 = _N1_5, _NE1_5 = _N1_5, NE1_5 = _N1_5,
 _N2_1 ,    N2_1 = _N2_1, _NE2_1 = _N2_1, NE2_1 = _N2_1,
 _N2_2 ,    N2_2 = _N2_2, _NE2_2 = _N2_2, NE2_2 = _N2_2,
 _N2_3 ,    N2_3 = _N2_3, _NE2_3 = _N2_3, NE2_3 = _N2_3,
 _N2_4 ,    N2_4 = _N2_4, _NE2_4 = _N2_4, NE2_4 = _N2_4,
 _N2_5 ,    N2_5 = _N2_5, _NE2_5 = _N2_5, NE2_5 = _N2_5
};

The naming convention is based on the FE periodic table. Note that there is an equivalence between NF1_k and RT_k (Raviart Thomas), NF2_k BDM_k (Brezzi Douglas Marini), N1_k and NE1_k and, N2_k and NE2_k.

Hint

NedelecEdgeTetrahedron and NedelecFaceTetrahedron classes use a general process to build shape functions as polynomials related to moment DoFs. To match DoFs on shared edge or shared face, the ascending numbering of vertices is implicitly used even if it is not in fact. This method is particular touchy in case of face DoFs of Nedelec Edge element. Indeed, such DoFs depend on the choice of two tangent vectors, generally two edges of faces of the reference tetrahedron that are mapped to some edges of physical element in a non-trivial way. To ensure a correct matching of such DoFs on a shared face, a rotation has to be applied to shape values to guarantee that they correspond to the same tangent vectors. This is the role of the member function NedelecEdgeFirstTetrahedronPk::rotateDofs.

Note that if the element vertex numbering and face vertex numbering are ascending, rotateDofs does nothing. This trick is commonly used to overcome this difficulty but as XLiFE++ makes no assumption on geometry, this DoFs rotation is mandatory.

More details on edge/face elements are provided in a specific paper.

Child classes for hexahedron

Child classes for prism

Child classes for pyramid

Examples of elements

Examples of triangles

Fig. 100 Triangle P1 Fig. 101 Triangle P1	Fig. 102 Triangle P2 Fig. 103 Triangle P2
Fig. 104 Triangle P3 Fig. 105 Triangle P3	Fig. 106 Triangle P4 Fig. 107 Triangle P4
Fig. 108 Triangle P5 Fig. 109 Triangle P5	Fig. 110 Triangle P6 Fig. 111 Triangle P6

Examples of quadrangles

Fig. 112 Quadrangle Q1 Fig. 113 Quadrangle Q1	Fig. 114 Quadrangle Q2 Fig. 115 Quadrangle Q2
Fig. 116 Quadrangle Q3 Fig. 117 Quadrangle Q3	Fig. 118 Quadrangle Q4 Fig. 119 Quadrangle Q4
Fig. 120 Quadrangle Q5 Fig. 121 Quadrangle Q5	Fig. 122 Quadrangle Q6 Fig. 123 Quadrangle Q6

Examples of tetrahedra

Fig. 124 Tetrahedron P1 Fig. 125 Tetrahedron P1	Fig. 126 Tetrahedron P2 Fig. 127 Tetrahedron P2
Fig. 128 Tetrahedron P3 Fig. 129 Tetrahedron P3	Fig. 130 Tetrahedron P4 Fig. 131 Tetrahedron P4

Examples of hexahedra

Fig. 132 Hexahedron Q1: Fig. 133 Hexahedron Q1:	Fig. 134 Hexahedron Q2: Fig. 135 Hexahedron Q2:
Fig. 136 Hexahedron Q3: Fig. 137 Hexahedron Q3:

Examples of prisms

Fig. 138 Prism P1Q1: Fig. 139 Prism P1Q1:

Fig. 140 Prism P2Q2: Fig. 141 Prism P2Q2:

Examples of pyramids

Fig. 142 Pyramid P1Q1: Fig. 143 Pyramid P1Q1:

Fig. 144 Pyramid P2Q2: Fig. 145 Pyramid P2Q2:

Integration methods

To perform computation of integrals over elements, the library provides IntegrationMethod class which is an abstract class. Two abstract classes inherit from it : the SingleIM class to deal with single integral and the DoubleIM class to deal with double integral. All the classes have to inherit from these two classes.

Among integration methods, quadrature methods are general and well adapted to regular functions to be integrated. XLiFE++ provides the QuadratureIM class, inherited from SingleIM class, which encapsulates the usual quadrature formulae defined in QuadratureIM class. Besides, the ProductIM class, inherited from doubleIM class, manages a pair of SingleIM objects and may be used to take into double quadrature method (adapted to compute double integral with regular kernel).

The inheritance diagram looks like

Table 3 :widths: 20 5 20 5 40 :header-rows: 0

IntegrationMethod (abstact)	\(\longrightarrow\)	SingleIM (abstract)	\(\longrightarrow\)	QuadratureIM (based on Quadrature class) PolynomialIM (exact intg. of polynoms) LenoirSalles2dIR (analytic method for IR, Laplace P0, P1) LenoirSalles3dIR (analytic method for IR, Laplace triangle P0, P1) FilonIM (1D oscillatory integrals intg_0_T f(t)exp(-iat))
	\(\longrightarrow\)	DoubleIM (abstract)	\(\longrightarrow\)	ProductIM (product of single intg methods) LenoirSalles2DIM (analytic method for BEM Laplace P0 and P1) LenoirSalles3DIM (analytic method for BEM Laplace triangle P0, P1) SauterSchwabIM (quadrature for BEM Laplace, Helmholtz triangle) SauterSchwabSymIM (Sauter-Schwab method for symmetrical kernels) DuffyIM (quadrature for BEM Laplace, Helmholtz segment) CollinoIM (semi-analytic method for BEM Maxwell RT0)

As some computations, in particular BEM computations, require more than one integration method, XLiFE++ provides the IntegrationMethods class that collects IntegrationMethod objects with additional parameters (IntgMeth class).

The IntegrationMethod, SingleIM, DoubleIM classes

The IntegrationMethod class

The IntegrationMethod class, basis class of all integration methods, manages a name, a type of integral and two booleans:

class IntegrationMethod
{
  public :
   String name;
   IntegrationMethodType imType;
   bool requireRefElement;
   bool requirePhyElement;
};

Up to now, the types of integral defined in enumeration are:

enum IntegrationMethodType
  {
   _undefIM, _quadratureIM, _polynomialIM, _productIM, _LenoirSalles2dIM, Lenoir_Salles_2d =_LenoirSalles2dIM, _LenoirSalles3dIM, Lenoir_Salles_3d=_LenoirSalles3dIM,  _LenoirSalles2dIR, _LenoirSalles3dIR, _SauterSchwabIM, Sauter_Schwab=_SauterSchwabIM, _SauterSchwabSymIM, Sauter_Schwab_sym=_SauterSchwabSymIM, _DuffyIM, Duffy = _DuffyIM, _DuffySymIM, Duffy_sym = _DuffySymIM, _HMatrixIM, H_Matrix= _HMatrixIM, _CollinoIM, _FilonIM
  };

The _HMatrixIM type is a very particular case because it is used to choose the HMatrix representation and compression methods which can be seen in some way as an integration method!

This class provides few general functions, most of them being virtual ones

IntegrationMethod(IntegrationMethodType imt=_undefIM, bool ref=false, bool phy=false, const String& na="")
virtual ~IntegrationMethod(){};
virtual bool isSingleIM() const =0;
virtual bool isDoubleIM() const =0;
virtual void print(std::ostream& os) const;
IntegrationMethodType type() const;
virtual bool useQuadraturePoints() const;

std::ostream& operator<<(std::ostream&, const IntegrationMethod&);

The SingleIM, DoubleIM classes

The SingleIM, DoubleIM derived classes provides only simple member functions:

class SingleIM : public IntegrationMethod
{
  public :
   bool isSingleIM() const {return true;}
   bool isDoubleIM() const {return false;}
   SingleIM(IntegrationMethodType imt=_undefIM);
   virtual void print(std::ostream& os) const;
}

class DoubleIM : public IntegrationMethod
{
 public :
  bool isSingleIM() const {return false;}
  bool isDoubleIM() const {return true;}
  DoubleIM(IntegrationMethodType imt=_undefIM);
  virtual void print(std::ostream& os) const;
}

The QuadratureIM class

The QuadratureIM class handles a list of QuadratureIM pointers indexed by ShapeType to deal with multiple quadrature rule in case of domain having more than one shape type. As the quadrature rule are often involved in context of FE computation, the QuadratureIM class may store shape values (pointer):

class QuadratureIM : public SingleIM
{
 protected :
  std::map<ShapeType, Quadrature *> quadratures_;
  std::map<ShapeType, std::vector<ShapeValues>* > shapevalues_;
 ...
}

It provides two basic constructors, accessors and print facilities (all are public):

QuadratureIM( ShapeType, QuadRule =_defaultRule, Number =0);
QuadratureIM(const std::set<ShapeType>&, QuadRule =_defaultRule, Number =0);
virtual bool useQuadraturePoints() const;
Quadrature* getQuadrature(ShapeType) const;
void setShapeValues(ShapeType, std::vector<ShapeValues>*);
std::vector<ShapeValues>* getShapeValues(ShapeType);
virtual void print(std::ostream& os) const;
std::ostream& operator<<(std::ostream&, const QuadratureIM&);

The ProductIM class

In order to represent a product of integration method over a product of geometric domain, the ProductIM class carries two SingleIM pointers:

class ProductIM : public DoubleIM
{
 protected :
  SingleIM* im_x;     //!< integration method along x
  SingleIM* im_y;     //!< integration method along y
}

and provides simple member functions:

ProductIM(SingleIM* imx=0, SingleIM* imy=0);
SingleIM* getxIM();
SingleIM* getyIM();
virtual bool useQuadraturePoints() const;
virtual void print(std::ostream& os) const;  //!< print on stream

The Quadrature, QuadratureRule classes

The quadrature formulae have the following form:

\[\int_{\widehat{E}} f(\widehat{x})d\widehat{x}\approx \sum_{i=1,q} \omega_i\,f(\widehat{x}_i)\]

where \((\widehat{x}_i)_{i=1,q}\) are quadrature points belonging to reference element \(\hat{E}\) and \((w_i)_{i=1,q}\) are quadrature weights.

Up to now, there exist quadrature formulae for unit segment \(]0,1[\), for unit triangle, for unit quadrangle (square), for unit tetrahedron, for unit hexahedron (cube) and for unit prism. The following tables give the list of quadrature formulae :

General rules

	Gauss-Legendre	Gauss-Lobatto	Grundmann-Muller	symmetrical Gauss
segment	any odd degree	any odd degree
quadrangle	any odd degree	any odd degree		any odd degree up to 21
triangle	any odd degree		any odd degree	degree up to 10
hexahedron	any odd degree	any odd degree		odd degree up to 11
tetrahedron	any odd degree -		any odd degree	degree up to 10
prism				degree up to 10
pyramid	any odd degree	any odd degree		degree up to 10

Particular rules

	nodal	miscellaneous
segment	P1 to P4
quadrangle	Q1 to Q4
triangle	P1 to P3	Hammer-Stroud 1 to 6
hexahedron	Q1 to Q4
tetrahedron	P1, P3	Stroud 1 to 5
prism	P1	centroid 1, tensor product 1,3,5
pyramid	P1	centroid 1, Stroud 7

More precisely, the rules that are implemented :

• 1D rules

  • Gauss-Legendre rule, degree \(2n+1\), \(n\) points

  • Gauss-Lobatto rule, degree \(2n-1\), \(n\) points

  • Gauss-Lobatto rule, degree \(2n-1\), \(n\) points

  • trapezoidal rule, degree \(1\), \(2\) points

  • Simpson rule, degree \(3\), \(3\) points

  • Simpson \(3\) \(8\) rule, degree \(3\), \(4\) points

  • Booleb rule, degree \(5\), \(5\) points

  • Wedge rule, degree \(5\), \(6\) points

• Rules over the unit simplex (less quadrature points BUT with negative weights)

  • TN Grundmann-Moller rule, degree \(2n+1\), \(C^n_{n+d+1}\) points

• Rules over the unit triangle \(\{ x > 0, y >0, x+y < 1 \}\)

  • T2P2 MidEdge rule

  • T2P2 Hammer-Stroud rule

  • T2P3 Albrecht-Collatz rule, degree \(3\), \(6\) points, Stroud (p.314)

  • T2P3 Stroud rule, degree \(3\), \(7\) points, Stroud (p.308)

  • T2P5 Radon-Hammer-Marlowe-Stroud rule, degree \(6\), \(7\) points, Stroud (p.314)

  • T2P6 Hammer rule, degree \(6\), \(12\) points, G.Dhatt, G.Touzot (p.298)

  • symmetrical Gauss up to degree 10

• Rules over the unit tetrahedron \(\{x > 0, y >0, z>0, x+y+z < 1\}\)

  • T3P2 Hammer-Stroud rule, degree 2, 4 points, Stroud (p.307)

  • T3P3 Stroud rule, degree \(3\), \(8\) points, Stroud (p.308)

  • T3P5 Stroud rule, degree \(5\), \(15\) points, Stroud (p.315)

  • symmetrical Gauss up to degree 10

• Rules over the unit prism \(\{ 0 < x, y < 1, x+y<1, 0 < z < 1\}\)

  • P1 nodal rule

  • tensor product of 2-points Gauss-Legendre rule on [0,1] and Stroud 7 points formula on the unit triangle (degree 3)

  • tensor product of 3-points Gauss-Legendre rule on [0,1] and Radon-Hammmer-Marlowe-Stroud 7 points formula on the unit triangle (degree 5)

  • symmetrical Gauss up to degree 10

• Rules over the unit pyramid \(\{ 0 < x, y < 1-z, 0 < z < 1\}\)

  • P1 nodal rule

  • symmetrical Gauss up to degree 10

• Rules built from lower dimension rules

  • tensor product of quadrature rules (quadrangle and hexahedron)

  • conical product of quadrature rules (triangle and pyramid)

  • rule built from 1D nodal rule which are quadrangle nodal rule

  • rule built from 1D nodal rule which are hexahedron nodal rule

References :

• A.H. Stroud, Approximate calculation of multiple integrals, Prentice Hall, 1971.

• G.Dhatt & G Touzot, The finite element method displayed, John Wiley & Sons, Chichester, 1984

• 1. Grundmann & H.M. Moller, Invariant Integration Formulas for the N-Simplex by Combinatorial Methods, SIAM Journal on Numerical Analysis, Vol 15, No 2, Apr. 1978, pages 282-290.

• F.D. Witherden & P.E. Vincent, On the identification of symmetric quadrature rules for finite element methods, Computers & Mathematics with Applications, Volume 69, Issue 10, May 2015, Pages 1232-1241.

These quadrature formulae are provided by using two classes : the Quadrature class which handles the main quadrature objects and the QuadratureRule class carrying the quadrature points and weights.

The Quadrature class

The Quadrature class manages all information related to quadrature :

class Quadrature
{
  public:
   GeomRefElement* geomRefElt_p;  // geometric element bearing quadrature rule
   QuadratureRule quadratureRule; // embedded QuadratureRule object
  protected:
   QuadRule rule_;           // type of quadrature rule
   Number degree_;           // degree of quadrature rule
   bool hasPointsOnBoundary_;// as it reads
   String name_;             // name of quadrature formula
  public:
   static std::vector<Quadrature*> theQuadratureRules;
  ...
}

QuadRule is an enumeration of all types of quadrature supported:

enum QuadRule
{
  _defaultRule = 0, _GaussLegendreRule, _symmetricalGaussRule, _GaussLobattoRule, Gauss_Lobatto = _GaussLobattoRule, _nodalRule, _miscRule,_GrundmannMollerRule, _doubleQuadrature
};

degree_ is either the degree of the quadrature rule or the degree of the polynomial interpolation in case of nodal quadrature rules. A Quadrature object has to be instanced only once and the static vector theQuadratureRules collects all the pointers to Quadrature object. Note that the shape bearing the quadrature is given through the GeomRefElement pointer.

The class proposes only a default constructor, a full basic constructor and a destructor :

Quadrature();
Quadrature(ShapeType , QuadRule, Number, const String&, bool pob = false);

Note that the basic constructor initializes member attributes but does not load the quadrature points and quadrature weights (QuadratureRule)! The practical “construction” is done by the external function findQuadrature(). The copy constructor and assignment operator are defined as private member function to avoid duplicated Quadrature object

The following accessors (read only) are provided:

String name() const;
QuadRule rule() const;
Number degree() const;
bool hasPointsOnBoundary() const;
Dimen dim() const;
Number numberOfPoints() const;
ShapeType shapeType() const;
std::vector<Real>::const_iterator point(Number i = 0) const;
std::vector<Real>::const_iterator weight(Number i = 0) const;

Warning

To construct a Quadrature object, the external function findQuadrature has to be used friend Quadrature* findQuadrature(ShapeType shape, QuadRule rule, Number degree, bool hpob= false);

It searches in the static vector theQuadratureRules the quadrature rule defined by a shape, a quadrature rule and a degree. If it is not found, a new Quadrature object is created on the heap. In any case a pointer to the Quadrature object is returned.

There are few functions depending on the shape that are called during the creation process :

Quadrature* segmentQuadrature(QuadRule, Number);
Quadrature* triangleQuadrature(QuadRule, Number);
Quadrature* quadrangleQuadrature(QuadRule, Number);
Quadrature* tetrahedronQuadrature(QuadRule, Number);
Quadrature* hexahedronQuadrature(QuadRule, Number);
Quadrature* prismQuadrature(QuadRule, Number);
Quadrature* pyramidQuadrature(QuadRule, Number);
static void clear();

These functions load quadrature points and quadrature weights in a QuadratureRule object using member functions of QuadratureRule class. The clear() static function deletes quadrature objects stored in the theQuadratureRules static vector.

To choose quadrature rules, the class provides the static function bestQuadRule() returning the ‘best’ quadrature rule for a given shape \(S\) and a given polynomial degree \(d\).

static QuadRule bestQuadRule(ShapeType, Number);

Best rule has to be understood as the rule on shape \(S\) with the minimum of quadrature points integrating exactly polynomials of degree \(d\).

The following table gives the current best rules :

shape	degree	quadrature rule	number of points
segment	\(1\)	Gauss-Legendre \(1\)	\(1\)
segment	\(2,3\)	Gauss-Legendre \(3\)	\(2\)
segment	\(4,5\)	Gauss-Legendre \(5\)	\(3\)
segment	\(2n-1\)	Gauss-Legendre \(2n-1\)	\(n\)

shape	degree	quadrature rule	number of points
quadrangle	\(1\)	Gauss-Legendre \(1\)	\(1\)
quadrangle	\(2,3\)	Gauss-Legendre \(3\)	\(4\)
quadrangle	\(4,5\)	symmetrical Gauss \(5\)	\(8\)
quadrangle	\(6,7\)	symmetrical Gauss \(7\)	\(12\)
quadrangle	\(7,9\)	symmetrical Gauss \(9\)	\(20\)
quadrangle	\(10,11\)	symmetrical Gauss \(11\)	\(28\)
quadrangle	\(12,13\)	symmetrical Gauss \(13\)	\(37\)
quadrangle	\(14,15\)	symmetrical Gauss \(15\)	\(48\)
quadrangle	\(16,17\)	symmetrical Gauss \(17\)	\(60\)
quadrangle	\(18,19\)	symmetrical Gauss \(19\)	\(72\)
quadrangle	\(20,21\)	symmetrical Gauss \(21\)	\(85\)
quadrangle	\(2n-1>21\)	Gauss-Legendre \(2n-1\)	\(n^2\)

shape	degree	quadrature rule	number of points
triangle	\(1\)	centroid rule (misc, \(1\))	\(1\)
triangle	\(2\)	P2 Hammer-Stroud (misc, \(2\))	\(2\)
triangle	\(3\)	Grundmann-Moller \(3\)	\(3\)
triangle	\(4\)	symmetrical Gauss \(4\)	\(6\)
triangle	\(5\)	symmetrical Gauss \(5\)	\(7\)
triangle	\(6\)	symmetrical Gauss \(6\)	\(12\)
triangle	\(7\)	symmetrical Gauss \(7\)	\(15\)
triangle	\(8\)	symmetrical Gauss \(8\)	\(16\)
triangle	\(9\)	symmetrical Gauss \(9\)	\(19\)
triangle	\(10\)	symmetrical Gauss \(10\)	\(25\)
triangle	\(11\)	symmetrical Gauss \(11\)	\(28\)
triangle	\(12\)	symmetrical Gauss \(12\)	\(33\)
triangle	\(13\)	symmetrical Gauss \(13\)	\(37\)
triangle	\(14\)	symmetrical Gauss \(14\)	\(42\)
triangle	\(15\)	symmetrical Gauss \(15\)	\(49\)
triangle	\(16\)	symmetrical Gauss \(16\)	\(55\)
triangle	\(17\)	symmetrical Gauss \(17\)	\(60\)
triangle	\(18\)	symmetrical Gauss \(18\)	\(67\)
triangle	\(19\)	symmetrical Gauss \(19\)	\(73\)
triangle	\(20\)	symmetrical Gauss \(20\)	\(79\)
triangle	\(2n-1>20\)	Gauss-Legendre \(2n-1\)	\(?\)

shape	degree	quadrature rule	number of points
hexahedron	\(1\)	Gauss-Legendre \(1\)	\(1\)
hexahedron	\(2,3\)	symmetrical Gauss \(3\)	\(6\)
hexahedron	\(4,5\)	symmetrical Gauss \(5\)	\(14\)
hexahedron	\(6,7\)	symmetrical Gauss \(7\)	\(34\)
hexahedron	\(8,9\)	symmetrical Gauss \(9\)	\(58\)
hexahedron	\(10,11\)	symmetrical Gauss \(11\)	\(90\)
hexahedron	\(2n-1\ge 12\)	Gauss-Legendre \(2n-1\)	\(n^3\)

shape	degree	quadrature rule	number of points
tetrahedron	\(1\)	centroid rule (misc, \(1\))	\(1\)
tetrahedron	\(2\)	P2 Hammer-Stroud (misc, \(2\))	\(4\)
tetrahedron	\(3\)	Grundmann-Moller \(3\)	\(5\)
tetrahedron	\(4,5\)	symmetrical Gauss \(5\)	\(14\)
tetrahedron	\(6\)	symmetrical Gauss \(6\)	\(24\)
tetrahedron	\(7\)	symmetrical Gauss \(7\)	\(35\)
tetrahedron	\(8\)	symmetrical Gauss \(8\)	\(46\)
tetrahedron	\(9\)	symmetrical Gauss \(9\)	\(59\)
tetrahedron	\(10\)	symmetrical Gauss \(10\)	\(81\)
tetrahedron	\(2n-1\ge 11\)	Grundmann-Moller \(2n-1\)	\(?\)

shape	degree	quadrature rule	number of points
prism	\(1\)	centroid rule (misc, \(1\))	\(1\)
prism	\(2\)	symmetrical Gauss \(2\)	\(5\)
prism	\(3\)	symmetrical Gauss \(3\)	\(8\)
prism	\(4\)	symmetrical Gauss \(4\)	\(11\)
prism	\(5\)	symmetrical Gauss \(5\)	\(16\)
prism	\(6\)	symmetrical Gauss \(6\)	\(28\)
prism	\(7\)	symmetrical Gauss \(7\)	\(35\)
prism	\(8\)	symmetrical Gauss \(8\)	\(46\)
prism	\(9\)	symmetrical Gauss \(9\)	\(60\)
prism	\(10\)	symmetrical Gauss \(10\)	\(85\)

shape	degree	quadrature rule	number of points
pyramid	\(1\)	centroid rule (misc, \(1\))	\(1\)
pyramid	\(2\)	symmetrical Gauss \(2\)	\(5\)
pyramid	\(3\)	symmetrical Gauss \(3\)	\(6\)
pyramid	\(4\)	symmetrical Gauss \(4\)	\(10\)
pyramid	\(5\)	symmetrical Gauss \(5\)	\(15\)
pyramid	\(6\)	symmetrical Gauss \(6\)	\(24\)
pyramid	\(7\)	symmetrical Gauss \(7\)	\(31\)
pyramid	\(8\)	symmetrical Gauss \(8\)	\(47\)
pyramid	\(9\)	symmetrical Gauss \(9\)	\(62\)
pyramid	\(10\)	symmetrical Gauss \(10\)	\(83\)

Finally, Quadrature class provides some print and error utilities:

void badNodeRule(int n_nodes);
void badDegreeRule();
void noEvenDegreeRule();
void upperDegreeRule();
friend std::ostream& operator<<(std::ostream&, const Quadrature&);
friend void alternateRule(QuadRule, ShapeType, const String&);

The QuadratureRule class

The QuadratureRule class stores quadrature points and quadrature weights and provides member functions to load them :

class QuadratureRule
{
 private:
  std::vector<Real> coords_;    // point coordinates of quadrature rule
  std::vector<Real> weights_;   // weights of quadrature rule
  Dimen dim_;                   // dimension of points
  ...
}

The following constructors are proposed:

QuadratureRule(const std::vector<Real>&, Real);
QuadratureRule(const std::vector<Real>&, const std::vector<Real>&);
QuadratureRule(Dimen , Number s);

and the following accessors are provided:

Dimen dim() const;
Number size() const;
std::vector<Real> coords() const;
std::vector<Real>::const_iterator point(Number i = 0) const;
std::vector<Real>::iterator point(Number i = 0);
std::vector<Real>::const_iterator weight(Number i = 0) const;
std::vector<Real>::iterator weight(const Number i = 0);

Some utilities are also helpful:

void point(std::vector<Real>::iterator&, Real, std::vector<Real>::iterator&, Real);
void point(std::vector<Real>::iterator&, Real, Real, std::vector<Real>::iterator&, Real);
void point(std::vector<Real>::iterator&, Real, Real, Real,std::vector<Real>::iterator&, Real);
void coords(std::vector<Real>::iterator&, Dimen, std::vector<Real>::const_iterator&);
void coords(std::vector<Real>::const_iterator);
void coords(const std::vector<Real>&);
void coords(Real);
void weights(const std::vector<Real>&);
void weights(Real);
friend std::ostream& operator<<(std::ostream&, const QuadratureRule&);

The main task of this class is to load quadrature points and weights. Thus, there are a lot of member functions to load these values. Some of them have a general purpose (for instance to build quadrature rules from tensor or conical products of 1D rules) and others are specific to one quadrature rule:

//tensor and conical product
 QuadratureRule& tensorRule(const QuadratureRule&, const QuadratureRule&);
 QuadratureRule& tensorRule(const QuadratureRule&, const QuadratureRule&, const QuadratureRule&);
 QuadratureRule& conicalRule(const QuadratureRule&, const QuadratureRule&);
 QuadratureRule& conicalRule(const QuadratureRule&, const QuadratureRule&, const QuadratureRule&);
//nodal tensor product rules for quadrangle and hexahedron
 QuadratureRule& quadrangleNodalRule(const QuadratureRule&);
 QuadratureRule& hexahedronNodalRule(const QuadratureRule&);
//1D Gauss rules
 QuadratureRule& gaussLegendreRules(Number);
 QuadratureRule& gaussLobattoRules(Number);
 QuadratureRule& ruleOn01(Number);
//1D Newton-Cotes closed rules (P_k node rules on [0,1])
 QuadratureRule& trapezoidalRule();
 QuadratureRule& simpsonRule();
 QuadratureRule& simpson38Rule();
 QuadratureRule& booleRule();
//Grundmann-Moller rules for n dimensional simplex
 QuadratureRule& tNGrundmannMollerRule(int, Dimen);
 void TNGrundmannMollerSet(int);    //utility for Grundmann-Moller
//special (Misc) rules for unit Triangle { x + y < 1 , 0 < x,y < 1}
 QuadratureRule& t2P2MidEdgeRule();
 QuadratureRule& t2P2HammerStroudRule();
 QuadratureRule& t2P3AlbrechtCollatzRule();
 QuadratureRule& t2P3StroudRule();
 QuadratureRule& t2P5RadonHammerMarloweStroudRule();
 QuadratureRule& t2P6HammerRule();
//special (Misc) rules for unit Tetrahedron { x + y + z < 1 , 0 < x,y,z < 1 }
 QuadratureRule& t3P2HammerStroudRule();
 QuadratureRule& t3P3StroudRule();
 QuadratureRule& t3P5StroudRule();
//special (Misc) rules for unit pyramid { 0 < x,y < 1 - z , 0 < z < 1 }
 QuadratureRule& pyramidRule(const QuadratureRule&);//conical product
 QuadratureRule& pyramidStroudRule();               //degree 7, 48 points
//symmetrical Gauss rules for any shape
 QuadratureRule& symmetricalGaussTriangleRule(number_t);    //degree up <= 20
 QuadratureRule& symmetricalGaussQuadrangleRule(number_t);  //odd degree <= 21
 QuadratureRule& symmetricalGaussTetrahedronRule(number_t); //degree <= 10
 QuadratureRule& symmetricalGaussHexahedronRule(number_t);  //odd degree <= 11
 QuadratureRule& symmetricalGaussPrismRule(number_t);       //degre <= 10
 QuadratureRule& symmetricalGaussPyramidRule(number_t);     //degree <= 10

Here is a simple example of how to use quadrature objects :

Quadrature* q=findQuadrature(_triangle,_GaussLegendreRule, 3);
Real S=0, x,y;
Number k=0;
for(Number i =0;i<q->quadratureRule.size();i++,k+=2)
{
  x=q->quadratureRule.coords()[k];
  y=q->quadratureRule.coords()[k+1];
  S+=q->quadratureRule.weights()[i]*std::sin(x*y);
}
std::cout<<"S="<<S;

How to add a new quadrature formula ?

When you want to add a new quadrature formula on a unit element, say emph{geomelt},

• you have to modify the member function :

  Quadrature::geomelt Quadrature (QuadRule rule, Number deg)

by adding in the right case statement the properties of your quadrature formula

• to add in QuadratureRule class the member function defining the quadrature points and weights.

Be care with the type and the number (degree) defining your quadrature formula, because they can be used by another formula. To avoid this difficulty, it is possible to create new type of quadrature. In that case the enumeration QuadRule has to be modified, implying to update the dictionaries and may be, the Quadrature::bestQuadRule() static function.

Integration methods for singular kernels

To deal with single or double integral involving a kernel with a singularity, XLiFE++ provides some additional classes.

DuffyIM class

This class deals with double integral over segments involving 2d kernel \(K\) with \(\log|x-y|\) singularity, in other words:

\[\int_{S_1}\int_{S_2} p_1(x)\,K(x,y)\,p_2(y)\,dy\,dx\]

For adjacent elements (sharing only one vertex) a Duffy transform is used. For self influence (\(S_1=S_2\)) two methods are proposed:

• one based on a splitting of the unit square \(]0,1[\times]0,1[\) in \(n_y\) decreasing rectangles \(]0,1[\times]\frac{1}{2^k},\frac{1}{2^{k-1}}[,\ k=1,n-1\) and \(]0,1[\times]0,\frac{1}{2^{n-1}}[\).

• one based on diagonal splitting and using a \(t^p\) transform in order to concentrate the quadrature points in the neighbor of the diagonal \(x=y\) (for more details see the separate documentation emph{quadrature_self_influence_2D}). This is the default one with \(p=5\).

class DuffyIM : public DoubleIM
{
 public:
 Quadrature* quadSelfX, *quadSelfY;//quadratures on segment [0,1] for self influence
 Quadrature* quadAdjtX, *quadAdjtY;//quadratures on segment [0,1] for adjacent elements
 number_t ordSelfX, ordSelfY;    //order of quadratures for self influence
 number_t ordAdjtX, ordAdjtY;    //order of quadratures for adjacent elements
 number_t selfOrd;               //order of t^p transform or number of number of layers(def=5)
 bool useSelfTrans;              //use improved method for self influence (def=true)

 DuffyIM(IntegrationMethodType);         //default constructor
 DuffyIM(number_t oX=6, number_t oY=6);  //full constructor
 DuffyIM(QuadRule qsX, number_t osX, QuadRule qsY, number_t osY, QuadRule qaX, number_t oaX, QuadRule qaY, number_t oaY, number_t so=5, bool ust=true);
 DuffyIM(const Quadrature&);             //full constructor from a quadrature
}

The class provides computation functions regarding geometrical configuration of segments:

template<typename K>
 void computeIE(...) const;
 void SelfInfluences(...) const;
 void SelfInfluencesTP(...) const;
 void AdjacentSegments(...) const;
 void k4(...) const;
 void k5(...) const;

The computation precision is very sensitive to the order of quadrature used in y direction. When using the first method to compute the self influence coefficients, to get very accurate results in the case of P2 mesh, very high order quadrature must be involved, for instance

DuffyIM dufim(_GaussLegendreRule,60,_GaussLegendreRule,60, _GaussLegendreRule,60,_GaussLegendreRule,60,5,false);

which lead to 300 quadrature points in each direction for self-influence coefficients!

Using the second method improves drastically the efficiency of the computation :

DuffyIM dufim(_GaussLegendreRule,30,_GaussLegendreRule,30, _GaussLegendreRule,50,_GaussLegendreRule,50,5,true);

The following picture shows the convergence rates of th BEM (single layer) in case of the diffraction (\(k=1.35\)) of a plane wave on a disk (radius 1). It can be noted that with a very low quadrature (red curve) the convergence rate is not captured at all, with a very high quadrature using the first method (cyan curve) the convergence rate becomes hard to capture when refining the mesh and finally the second method (black curve) allows to well capture the convergence rate with a reasonable number of quadrature points.

SauterSchwabIM class

This class deals with double integral over triangles involving 3d kernel \(K\) with \(|x-y|^{-1}\) singularity, say

\[\int_{T_1}\int_{T_2} p_1(x)\,K(x,y)\,p_2(y)\,dy\,dx\]

class SauterSchwabIM : public DoubleIM
{
 public:
  Quadrature* quadSelf;    //quadrature for self influence
  Quadrature* quadEdge;    //quadrature for elements adjacent by edge
  Quadrature* quadVertex;  //quadrature for elements adjacent by vertex
  number_t ordSelf;        //order of quadrature for self influence
  number_t ordEdge;        //order of quadrature for edge adjacence
  number_t ordVertex;      //order of quadrature for vertex adjacence

  SauterSchwabIM(number_t =3); //basic constructor (order 3)
  SauterSchwabIM(Quadrature&); //full constructor from quadrature object
}

LenoirSalles2dIM and LenoirSalles3dIM classes

The LenoirSalles2dIM class can compute analytically

\[-\frac{1}{2\pi}\int_{S_1}\int_{S_2} p_1(x)\,\log|x-y|\,p_2(y)\,dy\,dx \ \ \ \mathrm{(single\;layer)}\]

\[\frac{1}{2\pi}\int_{S_1}\int_{S_2} p_1(x)\,\frac{(x-y).n_y}{|x-y|^2}\,p_2(y)\,dy\,dx \ \ \ \text{(double layer)}\]

where \(p_1\) and \(p_2\) are either the P0 or P1 shape functions on segments \(S_1\) and \(S_2\) and the LenoirSalles3dIM class computes analytically

\[\frac{1}{4\pi}\int_{T_1}\int_{T_2} p_1(x)\,\frac{1}{|x-y|}\,p_2(y)\,dy\,dx \ \ \ \text{(single layer)}\]

\[-\frac{1}{4\pi}\int_{T_1}\int_{T_2} p_1(x)\,\frac{(x-y).n_y}{|x-y|^3}\,p_2(y)\,dy\,dx \ \ \ \text{(double layer)}\]

where \(p_1\) and \(p_2\) are the P0 shape functions (say \(1\)) on triangles \(T_1\) and \(T_2\).

These two classes have no member data and only one basic constructor.

LenoirSalles2dIR and LenoirSalles3dIR classes

The LenoirSalles2dIR class can compute analytically

\[-\frac{1}{2\pi}\int_{S} \log|x-y|\,p(y)\,dy\,dx \ \ \ \text{(single layer)}\]

\[\frac{1}{2\pi}\int_{S}\frac{(x-y).n_y}{|x-y|^2}\,p(y)\,dy\,dx \ \ \ \text{(double layer)}\]

where \(p\) is either the P0 or P1 shape functions on segment \(S\) and the LenoirSalles3dIR class that computes analytically

\[\frac{1}{4\pi}\int_{T}\frac{1}{|x-y|}\,p(y)\,dy\,dx \ \ \ \text{(single layer)}\]

\[-\frac{1}{4\pi}\int_{TS}\frac{(x-y).n_y}{|x-y|^3}\,p(y)\,dy\,dx \ \ \ \text{(double layer)}\]

where \(p\) is the P0 or P1 shape functions on triangle \(T\).

These two classes have no member data and only one basic constructor.

All the previous classes provides, a clone method, access to the list of quadratures if there are ones and an interface to the computation functions:

virtual LenoirSalles2dIM* clone() const;  //covariant
virtual std::list<Quadrature*> quadratures() const; //list of quadratures
template<typename K>
void computeIE(const Element* S, const Element* T, AdjacenceInfo& adj, const KernelOperatorOnUnknowns& kuv, Matrix<K>& res, IEcomputationParameters& iep, Vector<K>& vopu, Vector<K>& vopv, Vector<K>& vopk) const;

CollinoIM class

The CollinoIM class is a wrapper to an integration method developed by F. Collino to deal with some integrals involved in 3D Maxwell BEM when using Raviart-Thomas elements of order 1 (triangle). More precisely, it can compute the integrals:

\[\int_{\Gamma}\int_{\Gamma}\left(k\,H(k;x,y)\,wj(y).w_i(x)-\frac{1}{k}H(k;x,y)\,\text{div}\,w_j(y)\,\text{div}\,w_i(x)\right)\]

and

\[\int_{\Gamma}\int_{\Gamma}\left(\nabla_y H(k;x,y)\times w_j(y)\right).w_i(x),\]

where \((w_i)_{i=1,n}\) denote the Raviart-Thomas shape functions and \(H(k;x,y)\) the Green function of the Helmholtz equation (wave number k) in the 3D free space. The class looks like

enum ComputeIntgFlag {_computeI1=1, _computeI2,_computeBoth};

class CollinoIM : public DoubleIM
{
 private:
  myquad_t* quads_;
public:
 number_t ordSNear; //default 12
 number_t ordTNear; //default 64
 number_t ordTFar;  //default 3
 real_t eta;        //default 3.
 ComputeIntgFlag computeFlag;

 CollinoIM();
 CollinoIM(ComputeIntgFlag cf, number_t otf, number_t otn, number_t osn, real_t e);
 CollinoIM(const CollinoIM&);
 ~CollinoIM();
 virtual CollinoIM* clone() const;
 CollinoIM& operator=(const CollinoIM&);

 virtual std::list<Quadrature*> quadratures() const;
 virtual void print(std::ostream& os) const;
 template<typename K>
 void computeIE(const Element*, const Element*, AdjacenceInfo&, const KernelOperatorOnUnknowns&, Matrix<K>&, IEcomputationParameters&,Vector<K>&,Vector<K>&,Vector<K>&) const;
}

quads_ is a pointer to a quadrature structure provided by the Collino package, ordSNear, ordTNear, ordTFar are order of quadrature rules used, eta is a relative distance used to separate far and near triangle interactions and computeFlag allows to choose the integrals that have to be computed. computeIE() is the main computation function calling the real computation function (ElemTools_weakstrong_c()) provided by the Collino package (see files {itshape Collino.hpp} and Collino.cpp).

The IntegrationMethods class

The IntegrationMethods class collects some IntgMeth objects that handle an integration method pointer with some additional parameters:

class IntgMeth
{
  public:
   const IntegrationMethod* intgMeth;
   FunctionPart functionPart; //_allFunction, _regularPart, _singularPart
  real_t bound;              //bound value
}

The bound value may be used to select an integration method regarding a criterion. For instance, in BEM computation the criterion is the relative distance between the centroids of elements. Besides, using functionPart member, different integration method may be called on different part of the kernel, see Kernel class. This class has only some constructors and print stuff:

IntgMeth(const IntegrationMethod& im, FunctionPart fp=_allFunction, real_t b=0);
IntgMeth(const IntgMeth&);
~IntgMeth();
IntgMeth& operator=(const IntgMeth&);
void print(std::ostream&) const;
void print(PrintStream&) const;
friend ostream& operator<<(ostream&, const IntgMeth&)

For safe multi-threading behavior, copy constructor and assign operator do a full copy of IntegrationMethod object.

The IntegrationMethods class is simply a vector of IntMeth object:

class IntegrationMethods
{
 public:
  vector<IntgMeth> intgMethods;
  typedef vector<IntgMeth>::const_iterator const_iterator;
  typedef vector<IntgMeth>::iterator iterator;
 ...
}

This class provides many constructors depending on the number of integration methods given (up to 3 at this time), to some shortcuts :

IntegrationMethods() {};
IntegrationMethods(const IntegrationMethod&,FunctionPart=_allFunction,real_t=0);
...
void add(const IntegrationMethod&, FunctionPart=_allFunction,real_t=0);
void push_back(const IntgMeth&);

To handle more than 3 integration methods, use the add member function.

Besides, the class offers some accessors and print stuff:

const IntgMeth& operator[](number_t) cons ;
vector<IntgMeth>::const_iterator begin() const;
vector<IntgMeth>::const_iterator end() const;
vector<IntgMeth>::iterator begin();
vector<IntgMeth>::iterator end();
bool empty() const {return intgMethods.empty();}
void print(ostream&) const;
void print(PrintStream&) const;
friend ostream& operator<<(ostream&, const IntegrationMethods&);

In BEM computation, the following relative distance between two elements (\(E_i\), \(E_j\)) is used:

\[\DeclareMathOperator{\diam}{diam} dr(E_i,E_j)=\frac{||C_i-C_j||}{\max(\diam E_i,\diam E_j)}\mathrm{\;where\;}C_i\mathrm{\;is\;the\;centroid\;of\;} E_i.\]

So defining for instance

IntegrationMethods ims(SauterSchwabIM(4), 0., symmetrical_Gauss, 5, 1., symmetrical_Gauss, 3);

the BEM computation algorithm will perform on function (all part):

Sauter-Schwab with quadrature of order 4 on segment	if \(dr(E_i,E_j)=0\)
Symmetrical gauss quadrature of order 5 on triangle	if \(0 < dr(E_i,E_j) \leq 1\)
Symmetrical gauss quadrature of order 3 on triangle	if \(dr(E_i,E_j) > 1\)

Standard integration method

For various integration purpose, XLiFE++ provides the uniform rectangle, trapeze and Simpson method (\(h=\frac{b-a}{n}\)):

\[\begin{split}\begin{array}{lll} \textrm{rectangle}&: &\displaystyle \int_a^b f(t)\,dt\sim h\sum_{i=0,n-1} f(a+ih) \\[5mm] \textrm{trapeze} & : &\displaystyle \int_a^b f(t)\,dt\sim \frac{h}{2} \left(f(a)+4\sum_{i=1,n/2} f(a+(2i-1)h) +2\sum_{i=1,n/2-1} f(a+2ih)+f(b)\right)\ \ (n \textrm{ even})\\[5mm] \textrm{Simpson} & : &\displaystyle \int_a^b f(t)\,dt\sim \frac{h}{3} \left(f(a)+2\sum_{i=1,n-1} f(a+ih) + f(b)\right) \end{array}\end{split}\]

with \(f\) a real or complex function. The rectangle method integrates exactly P0 polynomials, the trapeze method integrates exactly P1 polynomials whereas the Simpson method integrates exactly P3 polynomials. They respectively approximate the integral with order 1, 2 and 3.

For each of them, four template functions (with names rectangle, trapz, simpson) are provided according to the user gives the values of \(f\) on a uniform set of points, the function \(f\) or a function \(f\) with some parameters (e.g. trapeze method):

template<typename T, typename Iterator>
T trapz(number_t n, real_t h, Iterator itb, T& intg);
template<typename T>
T trapz(const std::vector<T>& f, real_t h);
T trapz(T(*f)(real_t), real_t a, real_t b, number_t n);
T trapz(T(*f)(real_t, Parameters&), Parameters& pars, real_t a, real_t b, number_t n);

XLiFE++ provides also an adaptive trapeze method using the improved trapeze method (\(0.25(b-a)(f(a)+2f((a+b)/2)+f(b))\)) to get an estimator. The user has to give its function (with parameters if required), the integral bounds a tolerance factor (\(10^{-6}\) if not given):

template<typename T>
T adaptiveTrapz(T(*f)(real_t), real_t a, real_t b, real_t eps=1E-6);
T adaptiveTrapz(T(*f)(real_t,Parameters&), Parameters& pars, real_t a, real_t b, real_t eps=1E-6)

Hint

Because at the first step, the adaptive method uses 5 points uniformly distributed, the result may be surprising if unfortunately the function takes the same value at these points!

Discrete Fast Fourier Transform (FFT)

XLiFE++ provides the standard 1D discrete fast Fourier transform with \(2^n\) values:

\[g_k = \sum_{j=0,n-1} f_je^{-2i\pi\frac{jk}{n}}.\]

FFT and inverse FFT are computed using some functions addressing iterators on a collection of real or complex values:

template<class Iterator>
void  fft(Iterator ita, Iterator itb, number_t log2n);
void ifft(Iterator ita, Iterator itb, number_t log2n);

Computations from real or complex vectors is also available, but the results are always some complex vectors:

template<typename T>
vector<complex_t>&  fft(vector<T>& f, vector<complex_t>& g);
vector<complex_t>& ifft(vector<T>& f, vector<complex_t>& g);

The following example shows how to use them:

Number ln=6, n=64;
Vector<Complex> f(n),g(n),f2(n);
for(number_t i=0; i<n; i++) f[i]=std::sin(2*i*pi_/(n-1));
fft(f.begin(),g.begin(),ln);
ifft(g.begin(),f2.begin(),ln); //f2 should be very close to f
// or
fft(f,g); ifft(g,f2);

Integration methods for oscillatory integrals

Up to now, only one method is provided to compute 1D oscillatory integrals based on the Filon’s method.

FilonIM class

The Filon’s method XLiFE++ proposes is designed to compute an approximation of the following oscillatory integral:

\[I(x)=\int_0^T f(t)\,e^{-ixt}\]

where \(f\) is a slowly varying function with real or complex values (template type).

Considering a uniform grid \(t_n=n\,dt,\ \forall n=0,N\) and an interpolation space on this grid, say \(\text{vect }(w_n)_{n=0,N}\) where \(w_n\) are polynomial on each segment \([t_{n-1},t_n]\); the Filon’s method consists in interpolating the function \(f\):

\[f(t)=\sum_{n=1,N} a_n w_n(t).\]

Substituting \(f\) by its interpolation in integral and collecting terms in a particular way, the integral \(I(x)\) is thus approximated as follows:

\[I_N(x)=dt\sum_{j=1,p} C_j(x) \sum_{n=1,N} f(t_n^j)exp^{-ixt_{n-1}} + dt^2\sum_{j=1,p} C'_j(x) \sum_{n=1,N} f'(t_n^j)exp^{-ixt_{n-1}}\]

where \(Cj(x)\) and \(C'_j\) are coefficients of the form :math:` int_0^1 tau_j(s) exp^{-i,x,s,dt}ds` with \(\tau_j\) some shape functions on the segment \([0,1]\) ant \(t_n^j\) the support of the \(j^{th}\) DoF on the segment \([t_{n-1},t_n]\). The second sum in \(I_N\) appears only when using Hermite interpolation; for Lagrange interpolation, \(C'_j(x)=0\).

If the values \(f(t_n^j)_{n=0,N}\) and \(f'(t_n^j)_{n=0,N}\) (if Hermite interpolation) are pre-computed, the computation of \(I_N(x)\) for many values of \(x\) may be very efficient!

The FilonIMT is a template class, FilonIM being its complex specialization:

template <typename T = complex_t>
class FilonIMT : public SingleIM
{
 typedef T(*funT)(real_t);  // T function pointer alias
 protected:
  number_t ord_;        //interpolation order of (0,1 or 2)
  real_t tf_;           //upper bound of integral (lower bound is 0)
  real_t dt_;           //grid step if uniform
  number_t ns_;         //number of segments
  Vector<T> fn_;        //values of f (and df if required) on the grid
  Vector<real_t> tn_;   //grid points
  ...
}

ord_ corresponds to the interpolation used : 0 for Lagrange P0, 1 for Lagrange P1 and 2 for Hermite P3. The values \(f(t_n^j)_{n=0,N}\) and \(f'(t_n^j)_{n=0,N}\) are collected (only one copy) in the same vector with the following storage:

• Lagrange P0 : \(f(dt/2) f(3\,dt/2),\,\dots\, f((N-1/2)dt)\) (\(N\) values)

• Lagrange P1 : \(f(0) f(dt),\,\dots\, f(N\,dt)\) (\(N+1\) values)

• Hermite P3 : \(f(0) f'(0) f(dt) f'(dt),\,\dots\, f(N\,dt) f'(N\,dt)`\) (\(2(N+1)\) values)

The FilonIMT class provides several constructors that pre-compute the values of \(f\) and \(f'\) from user’s functions of the real variable with real or complex values. More precisely, the following types are available:

Real f(Real t) {...}
Complex f(Real t) {...}
Real f(Real t, Parameters& pars) {...}
Complex f(Real t, Parameters& pars) {...}

It is also possible to construct a FilonIMT object passing the values of \(f\) at uniformly distributed points on:math:` [0,t_f]`. In that case, Filon of order 1 is selected. So the class looks like

typedef T(*funT)(real_t);
typedef T(*funTP)(real_t, Parameters& pa);
FilonIMT();
FilonIMT(number_t o,number_t N,real_t tf,funT f);
FilonIMT(number_t o,number_t N,real_t tf,funT f,funT df);
FilonIMT(number_t o,number_t N,real_t tf,funTP f,Parameters& pars);
FilonIMT(number_t o,number_t N,real_t tf,funTP f,funTP df, Parameters& pars);
FilonIMT(const Vector<T>& vf, real_t tf);
void init(number_t N,funT f,funT df);
void init(number_t N,funTP f,funTP df,Parameters& pars);
number_t size() const {return tn_.size();}     //grid size
template <typename S>
complex_t coef(const S& x, number_t j) const;  //Filon's coefficients
complex_t compute(const S& x) const;           //compute I(x)
complex_t operator()(const S& x) const         //compute I(x)
virtual void print(std::ostream& os) const;    //print FilonIMT on stream

This class is quite easy to use, for instance to compute the integral

\[I(x)=\int_0^{10} e^{-t}\,e^{-ixt}\]

by the P1-Filon method:

Complex f(Real x) {return Complex(std::exp(-x),0.);}
...
FilonIM fim(1, 100, 10, f);
Complex y = fim(0.5);
// or in a more compact syntax if only one evaluation
Complex y = FilonIM(1, 100, 10, f) (0.5);

To invoke Hermite P3 Filon method, define also the derivative of \(f\):

Complex f(Real x)  {return Complex(exp(-x),0.);}
Complex df(Real x) {return Complex(-exp(-x),0.);}
...
FilonIM fim(2, 100, 10, f, df);
Complex y = fim(0.5);

Hint

The P0 Filon method is of order 1, the P1 Filon is of order 2 while the Hermite P3 Filon is of order 4. The Hermite P3 Filon is less stable than the Lagrange Filon methods!

See also

• library= finiteElements,

• header= Quadrature.hpp QuadratureRule.hpp IntegrationMethod.hpp,

• implementation= Quadrature.cpp QuadratureRule.cpp IntegrationMethod.cpp,

• test= test_Quadtrature.cpp,

• header dep= config.h, utils.h, LagrangeQuadrangle.hpp, LagrangeHexahedron.hpp, mathsResources.h.