Integration methods#

The IntegrationMethod class#

The IntegrationMethod class, basis class of all integration methods, manages general information about them:

class IntegrationMethod
{public :
  String name;                  //name for doc purpose
  IntegrationMethodType imType; //method type (see below)
  SingularityType singularType; //singularity supported (_notsingular,_r,_logr,_loglogr)
  real_t singularOrder;         //singularity order
  string_t kerName;             //kernel shortname if adapted to (default is empty)
  bool requireRefElement;       //true if method require reference element (false)
  bool requirePhyElement;       //true if method require physical element (false)
  bool requireNormal;           //true if method involves normal (false)
};

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

enum IntegrationMethodType
{ _undefIM, _quadratureIM,_polynomialIM,_productIM,
  _LenoirSalles2dIM,_LenoirSalles3dIM,_LenoirSalles2dIR,_LenoirSalles3dIR,
  _SauterSchwabIM,_SauterSchwabSymIM,_DuffyIM,_DuffySymIM,
  _HMatrixIM, _CollinoIM, _FilonIM
};

The _HMatrixIM method type is a very particular case because it is used to choose the hierarchical matrix representation (HMatrix class) and compression methods which can be seen in some way as an integration method! See H-matrices.

The SingleIM, DoubleIM inherited from the IntegrationMethod are purely interface classes to real integration methods for single or double integral.

The QuadratureIM class#

The QuadratureIM class is dedicated to integration methods adapted to integrals of regular functions, based on quadrature points and quadrature weights. It is derived from SingleIM class and mainly handles a list of Quadrature objects 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 Quadrature class may store shape values:

class QuadratureIM : public SingleIM
{public:
  QuadRule quadRule;         //QuadRule to be used (may be _defaultRule)
  Number degree;             //degree of quadrature (default 3)
  map<ShapeType,Quadrature*> quadratures_; //to store Quadrature obj for different shapes
  map<ShapeType,vector<ShapeValues>*> shapevalues_; //to store shape values
  ...
}

This class offers the setQuadrature and setQuadratures functions that set the Quadrature objects according to the shapes, the quadrature rule and the order provided.

The Quadrature and the QuadratureRule classes#

The quadrature formulae have the following form:

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

The Quadrature class manages all information related to quadrature :

class Quadrature
{public:
  GeomRefElement* geomRefElt_p;  //pointer to geometric reference element
  QuadratureRule quadratureRule; //embedded QuadratureRule object
  QuadRule rule;                 //type of quadrature rule
  number_t degree;               //degree of quadrature rule
  bool hasPointsOnBoundary;      //as it reads
  string_t name;                 //name of quadrature formula (for doc purposes)
  static std::vector<Quadrature*> theQuadratureRules; //global list of Quadrature
  ...
}

where

• geomRefElt_p refers to the geometry (segment, triangle, …) bearing the quadrature.

• quadratureRule refers to a QuadratureRule object which handles the points and weights of quadrature and the functions to build them.

• rule refers to a QuadRule which is an enumeration of all types of quadrature supported:

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

• degree_ is either the degree of the quadrature rule or the degree of the polynomial interpolation in case of nodal quadrature rules.

• the static vector theQuadratures collects all the pointers to Quadrature object in order to have only one instance of a Quadrature object.

The practical “construction” of a Quadrature object is done by the external function findQuadrature which first, searches in the static vector theQuadratures the quadrature 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.

This class offers also the bestQuadrule function that returns the “best” quadrature rule regarding the shape and the degree.

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

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

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.

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;

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) for unit prism and for unit pyramid. The following tables give the list of available quadrature formulae :

General rules

Particular rules

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+d+1}^n$ 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)

  • Gauss-Legendre (conical product), any order

  • symmetrical Gauss up to degree 10

• Rules over the unit quadrangle $[0,1]\times [0,1]$

  • nodal degree 1 to 4 (1D nodal, tensor product)

  • Gauss-Legendre (tensor product), any order

  • Gauss-Lobatto (tensor product), any order

• 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)

  • Gauss-Legendre (conical product), any order

  • symmetrical Gauss up to degree 10

• Rules over the unit hexahedron $[0,1]\times [0,1]\times [0,1]$

  • nodal degree 1 to 4 (1D nodal, tensor product)

  • Gauss-Legendre (tensor product), any order

  • Gauss-Lobatto (tensor product), any order

  • symmetrical Gauss up to degree 11

• 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)

  • Gauss-Legendre (tensor product), any order

  • Gauss-Lobatto (tensor product), any order

  • 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.

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 :

How to add a new quadrature formula ?#

To add a new quadrature formula, say myRule on a unit element, say geomelt:

• first, add myRule to QuadRule enumeration type

• next, add in QuadratureRule class the member function defining the quadrature points and weights, e.g QuadratureRule::myQuadRule.

• finally, modify the geomeltQuadrature function by adding in the right case statement the properties of the new quadrature formula. For instance, in the case of a new quadrature formula on the triangle:

Quadrature* triangleQuadrature(QuadRule rule, number_t numb)
{...
 switch(appliedRule)
 {case _nodalRule           : ...; break;
  case _GaussLegendreRule   : ...; break;
  case _GrundmannMollerRule : ...; break;
  case _symmetricalGaussRule: ...; break;
  case _miscRule            : ...; break;
  case _myRule :
    q=findQuadrature(_triangle, _myRule, numb, false); //find if already loaded
    if(q==nullptr) // create new Quadrature
    {
      q=new Quadrature(_triangle,_myRule,numb,"my rule",false);
      q->quadratureRule.myQuadRule(numb); //set points and weights
    }
    break;
  default: badDegreeRule(numb,"miscellaneous rule",_triangle); break;
 }
 ...
 return q;
}

The ProductIM class#

In order to represent a product of integration method over a product of geometric domain (double integrals), 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
}

ProductIM objects are used to calculate double integrals involving kernels which are regular on the integration domain.

Integration methods for singular kernels#

To deal with single or double integral involving a kernel with a singularity on the integration domain, say

XLiFE++ provides some additional classes:

• DuffyIM (segments) and SauterSchwabIM (triangle) well adapted to Laplace/Helmholtz kernel and any regular functions $p1$, $p2$

• LenoirSalles2dIM (segments) and LenoirSalles3dIM (triangles) well adapted to Laplace kernel with $p1$, $p2$ constant or linear functions

• CollinoIM (triangles) designed to Maxwell 3d BEM integrals using Raviart-Thomas elements (affine).

• LenoirSalles2dIR (segments) and LenoirSalles3dIR (triangles) well adapted to calculate integral representation invoving Laplace kernel with $p1$, $p2$ constant or linear functions

DuffyIM class#

This class deals with double integral over segments ($S_1\times S_2$) involving 2d kernel $K$ with $\log |x-y|$ singularity.

For adjacent elements (sharing only one vertex) a Duffy transform is used (segments $[S_1,S_2]$ and $[T_1,T_2]$ with $T_1=S_1$):

$$\phi :(\eta ;\nu )\in [0,1]\times [0,1]\to (S_1+\eta (S_2-S_1);T_1+\eta \nu (T_2-T_1))\in [S_1,S_2]\times [T_1,T_2]$$

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\text{ 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 ordSelfX, ordSelfY;    //order of quadratures for self influence
 Number ordAdjtX, ordAdjtY;    //order of quadratures for adjacent elements
 Number 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 oX=6, Number oY=6);  //full constructor
 DuffyIM(QuadRule qsX, Number osX, QuadRule qsY, Number osY, QuadRule qaX, Number oaX, QuadRule qaY, Number oaY, Number 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.

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 ordSelf;        //order of quadrature for self influence
  Number ordEdge;        //order of quadrature for edge adjacence
  Number ordVertex;      //order of quadrature for vertex adjacence

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

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;

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

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

 CollinoIM();
 CollinoIM(ComputeIntgFlag cf, Number otf, Number otn, Number osn, Real 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;
}

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 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 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=0);
...
void add(const IntegrationMethod&, FunctionPart=_allFunction,Real=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) 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:

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):

Standard integration method#

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

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

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

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

Number ln=6, n=64;
Vector<Complex> f(n),g(n),f2(n);
for(Number 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);

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

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

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