inputs/dev/IntegrationMethod.tex

%%%%%%%%%%%%%%%%%%%
% XLiFE++ is an extended library of finite elements written in C++
%     Copyright (C) 2014  Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin
%
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.
%%%%%%%%%%%%%%%%%%%

\section{Integration methods}

To perform computation of integrals over elements, the library provides {\class IntegrationMethod} class which is an abstract class. Two abstract classes inherit from it : the {\class SingleIM} class to deal with single integral and the {\class 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. \xlifepp provides the {\class QuadratureIM} class, inherited fom {\class SingleIM} class, which encapsulates the usual quadrature formulae defined in {\class QuadratureIM} class. Besides, the  {\class ProductIM} class, inherited from  {\class doubleIM} class, manages a pair of {\class SingleIM} object and may be used to take into double quadrature method (adapted to compute double integral with regular kernel).\\

The inheritance diagram looks like\\
\footnotesize
\begin{verbatim}
IntegrationMethod --> | SingleIM (single intg) --> | QuadratureIM (based on Quadrature class)
    (abstract)        |       (abstract)           | PolynomialIM (exact intg. of polynoms)
                      |                            | LenoirSalles2dIR (for 2D Laplace kernel)
                      |                            | LenoirSalles3dIR (for 3D Laplace kernel)
                      |
                      | DoubleIM (double intg) --> | ProductIM (product of single intg meth.)
                               (abstract)          | LenoirSalles2dIM (for 2D Laplace kernel)
                                                   | LenoirSalles3dIM (for 3D Laplace kernel)
                                                   | SauterSchwabIM (for 3D IE singular integral)
                                                   | DuffyIM  (for 2D IE singular integral)
\end{verbatim}
\normalsize

As some computations, in particular BEM computation, require more than one integration method, \xlifepp provides the {\class IntegrationMethods} class that collects {\class IntegrationMethod}  with additional parameters ({\class IntgMeth} class).


\subsection{The {\classtitle IntegrationMethod, SingleIM, DoubleIM} classes}

\subsubsection*{The {\classtitle IntegrationMethod} class}

The {\class IntegrationMethod} class, basis class of all integration methods, manages a name, a type of integral and two booleans:\\
\vspace{.1cm}
\begin{lstlisting}[deletekeywords={[3] name}]
class IntegrationMethod
{public :
  String name;
  IntegrationMethodType imType;
  bool requireRefElement;
  bool requirePhyElement;
};
\end{lstlisting}
\vspace{.2cm}
Up to now, the types of integral defined in enumeration are:
\vspace{.1cm}
\begin{lstlisting}[]{}
enum IntegrationMethodType
  {_undefIM, _quadratureIM, _polynomialIM, _productIM,
   _LenoirSalles2dIM, Lenoir_Salles_2d =_LenoirSalles2dIM,
   _LenoirSalles3dIM, Lenoir_Salles_3d= _LenoirSalles3dIM,
   _LenoirSalles2dIR, _LenoirSalles3dIR,
   _SauterSchwabIM, Sauter_Schwab=_SauterSchwabIM,
   _DuffyIM, Duffy = _DuffyIM,
   _HMatrixIM, H_Matrix= _HMatrixIM
   };
\end{lstlisting}
\vspace{.2cm}
The type {\tt \_HMatrixIM} 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
\vspace{.1cm}
\begin{lstlisting}[]{}
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&);
\end{lstlisting}

\subsubsection*{The {\classtitle SingleIM, DoubleIM} classes}

The {\class SingleIM, DoubleIM} derived classes provides only simple member functions :
\vspace{.1cm}
\begin{lstlisting}[]{}
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;
}
\end{lstlisting}
\vspace{.2cm}

\begin{lstlisting}[]{}
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;
}
\end{lstlisting}
\vspace{.2cm}

\subsection{The {\classtitle QuadratureIM} class}\label{quadrature_section}

The {\class QuadratureIM} class handles a list of {\class QuadratureIM} pointers indexed by {\class ShapeType} to deal with multiple quadrature rule in case of domain having more than one shape type. As theq uadrature rule are often involved in context of FE computation, the {\class QuadratureIM} class may store shape values (pointer):\\
\vspace{.1cm}
\begin{lstlisting}[deletekeywords={[3] map}]
class QuadratureIM : public SingleIM
{protected :
 std::map<ShapeType, Quadrature *> quadratures_;
 std::map<ShapeType, std::vector<ShapeValues>* > shapevalues_;
\end{lstlisting}
\vspace{.2cm}
It provides two basic constructors, accessors and print facilities (all are public):
\vspace{.1cm}
\begin{lstlisting}[deletekeywords={[3] set}]
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&);
\end{lstlisting}
\vspace{.2cm}

\subsection{The {\classtitle ProductIM} class}

In order to represent a product of integration method over a product of geometric domain, the {\class ProductIM} class carries two {\class SingleIM} pointers:
\vspace{.1cm}
\begin{lstlisting}[]{}
class ProductIM : public DoubleIM
{protected :
 SingleIM* im_x;     //!< integration method along x
 SingleIM* im_y;     //!< integration method along y
}
\end{lstlisting}
\vspace{.2cm}
and provides simple member functions:
\vspace{.1cm}
\begin{lstlisting}[]{}

    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
};
\end{lstlisting}
\vspace{.2cm}

\subsection{The {\classtitle 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 :\\
\goodbreak

\begin{center}
\textbf{General rules }\\
\vspace{2mm}
\small
\begin{tabular}{|c|c|c|c|c|}
  \hline
  & Gauss-Legendre & Gauss-Lobatto & Grundmann-Muller & symmetrical Gauss\\
  \hline
  segment & any odd degree & any odd degree &   &  \\
	\hline
  quadrangle & any odd degree & any odd degree & & odd degree up to 21 \\
	\hline
	triangle & any odd degree &   &  any odd degree & degree up to 10\\
  \hline
	hexahedron & any odd degree & any odd degree  & & odd degree up to 11\\
	\hline
	tetrahedron & any odd degree &   & any odd degree & degree up to 10\\
  \hline
	prism &  &    &  & degree up to 10\\
  \hline
	pyramid &  any odd degree & any odd degree  &  &   degree up to 10 \\
  \hline
\end{tabular}
\end{center}


\begin{center}
\textbf{Particular rules}\\
\vspace{2mm}
\small
\begin{tabular}{|c|c|c|}
  \hline
  &  nodal & miscellaneous\\
  \hline
  segment & P1 to P4 &  \\
	\hline
  quadrangle & Q1 to Q4 &  \\
	\hline
	triangle & P1 to P3 & Hammer-Stroud 1 to 6\\
  \hline
	hexahedron & Q1 to Q4 &  \\
	\hline
	tetrahedron &  P1, P3 & Stroud 1 to 5\\
  \hline
	prism & P1  &  centroid 1, tensor product 1,3,5\\
  \hline
	pyramid & P1 & centroid 1, Stroud 7 \\
  \hline
\end{tabular}
\end{center}
More precisely, the rules that are implemented :
\begin{itemize}
\item 1D rules \\
 Gauss-Legendre 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

\item 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

\item 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

\item 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

 \item 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

\item Rules over the unit pyramid $\{ 0 < x, y < 1-z, 0 < z < 1\}$ \\
P1 nodal rule \\
symmetrical Gauss up to degree 10


\item 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
\end{itemize}

References :
\begin{itemize}
\item A.H.Stroud, Approximate calculation of multiple integrals, Prentice Hall, 1971.
\item G.Dhatt \& G Touzot, The finite element method displayed,  John Wiley \& Sons, Chichester, 1984
\item A. 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.
\item 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.
\end{itemize}

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

\subsection*{The {\classtitle Quadrature} class}

The {\class Quadrature} class manages all informations related to quadrature :
\vspace{.1cm}
\begin{lstlisting}
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;
 ...
\end{lstlisting}
\vspace{.2cm}
\verb?QuadRule? is a enumeration of all types of quadrature supported:
\vspace{.1cm}
\begin{lstlisting}[]
enum QuadRule
{_defaultRule = 0, _GaussLegendreRule, _symmetricalGaussRule, _GaussLobattoRule, Gauss_Lobatto = _GaussLobattoRule,
  _nodalRule, _miscRule,_GrundmannMollerRule, _doubleQuadrature };
\end{lstlisting}
\vspace{.1cm}
\verb?degree_? is either the degree of the quadrature rule or the degree of the polynomial interpolation in case of nodal quadrature rules. A {\class Quadrature} object has to be instanced only once and the static vector \verb?theQuadratureRules? collects all the pointers to \verb?Quadrature? object. Note that the shape bearing the quadrature is given through the {\class GeomRefElement} pointer. \\

The class proposes only a default constructor, a full basic constructor and a destructor :
\vspace{.1cm}
\begin{lstlisting}[]
Quadrature();
Quadrature(ShapeType , QuadRule, Number, const String&, bool pob = false);
\end{lstlisting}
\vspace{.1cm}
Note that the basic constructor initializes member attributes but does not load the quadrature points and quadrature weights ({\class QuadratureRule})!
The practical "construction" is done by the external function  \verb?findQuadrature?. The copy constructor and assignment operator are defined as private member function to avoid duplicated {\class Quadrature} object\\


The following accessors (read only) are provided:
\vspace{.1cm}
\begin{lstlisting}[]
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;
\end{lstlisting}
\vspace{.2cm}
\textdbend \ \  To construct a {\class Quadrature} object, the external function \verb?findQuadrature? has to be used :
\vspace{.1cm}
\begin{lstlisting}[]
friend Quadrature* findQuadrature(ShapeType shape, QuadRule rule, Number degree, bool hpob= false);
\end{lstlisting}
\vspace{.1cm}
It searches in the static vector \verb?theQuadratureRules? the quadrature rule defined by a shape, a quadrature rule and a degree. If it is not found, a new {\class Quadrature}
object is created on the heap. In any case a pointer to the {\class Quadrature} object is returned. There are few functions depending on the shape that are called during the creation process :
\vspace{.1cm}
\begin{lstlisting}[]
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();
\end{lstlisting}
\vspace{.1cm}
These functions loads quadrature points and quadrature weights in a {\class QuadratureRule} object using member functions of {\class QuadratureRule} class.
The \verb?clear? static function deletes quadrature objects stored in the  \verb?theQuadratureRules? static vector.\\

To choose quadrature rules, the class provides the static function \verb?bestQuadRule? returning the 'best' quadrature rule for a given shape $S$ and a given polynomial degree $d$.
\vspace{.1cm}
\begin{lstlisting}[]
static QuadRule bestQuadRule(ShapeType, Number);
\end{lstlisting}
\vspace{.1cm}
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 :
\begin{center}
\begin{tabular}{|c|c|c|c|}
\hline
shape& degree & quadrature rule & number of points\\
\hline
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$\\
\hline
\end{tabular}
\end{center}

\begin{center}
\begin{tabular}{|c|c|c|c|}
\hline
shape& degree & quadrature rule & number of points\\
\hline
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$\\
\hline
\end{tabular}
\end{center}
\begin{center}
\begin{tabular}{|c|c|c|c|}
\hline
shape& degree & quadrature rule & number of points\\
\hline
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$& $?$\\
\hline
\end{tabular}
\end{center}
\begin{center}
\begin{tabular}{|c|c|c|c|}
\hline
shape& degree & quadrature rule & number of points\\
\hline
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$\\
\hline
\end{tabular}
\end{center}
\begin{center}
\begin{tabular}{|c|c|c|c|}
\hline
shape& degree & quadrature rule & number of points\\
\hline
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$ & $?$\\
\hline
\end{tabular}
\end{center}
\begin{center}
\begin{tabular}{|c|c|c|c|}
\hline
shape& degree & quadrature rule & number of points\\
\hline
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$\\
\hline
\end{tabular}
\end{center}
\begin{center}
\begin{tabular}{|c|c|c|c|}
\hline
shape& degree & quadrature rule & number of points\\
\hline
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$\\
\hline
\end{tabular}
\end{center}
\mbox{} \\
Finally, {\class Quadrature} class provides some print and error utilities:
\vspace{.1cm}
\begin{lstlisting}[]
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&);
\end{lstlisting}
\vspace{.1cm}


\subsection*{The {\classtitle QuadratureRule} class}

The {\class QuadratureRule} class stores quadrature points and quadrature weights and provides member functions to load them :
\vspace{.1cm}
\begin{lstlisting}
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
...
\end{lstlisting}
\vspace{.2cm}
the following constructors are proposed:
\vspace{.1cm}
\begin{lstlisting}[]
QuadratureRule(const std::vector<Real>&, Real);
QuadratureRule(const std::vector<Real>&, const std::vector<Real>&);
QuadratureRule(Dimen , Number s);
\end{lstlisting}
\vspace{.2cm}
and the following accessors are provided :
\vspace{.1cm}
\begin{lstlisting}[]
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);
\end{lstlisting}
\vspace{.2cm}
Some utilities are also helpful:
\vspace{.1cm}
\begin{lstlisting}[]
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&);
\end{lstlisting}
\vspace{.2cm}
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:
\vspace{.1cm}
\begin{lstlisting}[]
//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
 \end{lstlisting}
\vspace{.2cm}
Here is a simple example of how to use quadrature objects :
\vspace{.1cm}
\begin{lstlisting}[deletekeywords={[3] x, y}]
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;
\end{lstlisting}
\vspace{.2cm}


\subsection{How to add a new quadrature formula ?}

When you want to add a new quadrature formula on a unit element, say \emph{geomelt},
\begin{itemize}
\item you have to modify the member function :\\
\ \ \ \verb?Quadrature::?\emph{geomelt}\verb?Quadrature(QuadRule rule, Number deg)? \\
by adding in the right \emph{case} statement the properties of your quadrature formula
\item to add in {\class QuadratureRule} class the member function defining the quadrature points and weights.
\end{itemize}
Be care with the type and the number (degree) defining your quadrature formula, because they can be used by an other formula.
To avoid this difficulty, it is possible to create new type of quadrature. In that case the enumeration \verb?QuadRule? has to be modified,
implying to update the dictionnaries and may be, the \verb?Quadrature::bestQuadRule? static function.\\

\subsection{Integration methods for singular kernels}
To deal with single or double integral involving a kernel with a singularity, \xlifepp provides some additional classes.

\subsection*{{\class 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$$
\vspace{.1cm}
\begin{lstlisting}[]
class DuffyIM : public DoubleIM
{
 public:
 Quadrature* quadSelf; //quadrature for self influence
 Quadrature* quadAdjt; //quadrature for adjacent elements
 number_t ordSelf;     //order of quadrature for self influence
 number_t ordAdjt;     //order of quadrature for adjacent elements

 DuffyIM(IntegrationMethodType);
 DuffyIM(number_t =6);       //full constructor from quadrature order
 DuffyIM(const Quadrature&); //full constructor from a quadrature
\end{lstlisting}

\subsection*{{\class 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$$
\vspace{.1cm}
\begin{lstlisting}[]
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
 \end{lstlisting}

\subsection*{{\class LenoirSalles2dIM} and {\class LenoirSalles3dIM} classes}
The {\class 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 \ \ \ \text{(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 {\class 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.

	\subsection*{{\class LenoirSalles2dIR} and {\class LenoirSalles3dIR} classes}
	The {\class 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 {\class LenoirSalles3dIR} classe 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, an access to the list of quadratures if there are ones and an interface to the computation functions:
\vspace{.1cm}
\begin{lstlisting}[]
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;
\end{lstlisting}

\subsection{The {\classtitle IntegrationMethods} class}
The {\class IntegrationMethods} class collects some {\class IntgMeth} objects that handle an integration method pointer with some additional parameters:
\vspace{.1cm}
\begin{lstlisting}[deletekeywords={[3] x, y}]
class IntgMeth
{ public:
  const IntegrationMethod* intgMeth;
  FunctionPart functionPart; //_allFunction, _regularPart, _singularPart
  real_t bound;              //bound value
}
\end{lstlisting}
\vspace{.2cm}
The {\var bound} value may be used to select an integration method regarding a criteria. For instance, in BEM computation the criteria is the relative distance between the centroids of elements. Besides, using {\var functionPart} member, different integration method may be called on different part of the kernel, see {\class Kernel} class. This class has
only some constructors and print stuff:
\vspace{.1cm}
\begin{lstlisting}[deletekeywords={[3] x, y}]
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&)
\end{lstlisting}
\vspace{.2cm}
For safe multithreading behaviour, copy constructor and assign operator do a full copy of {\class IntegrationMethod} object.\\

The {\class IntegrationMethods} class is simply a vector of {\class IntMeth} object:
\vspace{.1cm}
\begin{lstlisting}[deletekeywords={[3] x, y}]
class IntegrationMethods
{public:
 vector<IntgMeth> intgMethods;
 typedef vector<IntgMeth>::const_iterator const_iterator;
 typedef vector<IntgMeth>::iterator iterator;
 ...}
\end{lstlisting}
\vspace{.2cm}
This class provides many constructors depending to the number of integration methods given (up to 3 at this time), to some shortcuts :
\vspace{.1cm}
\begin{lstlisting}[deletekeywords={[3] x, y}]
IntegrationMethods() {};
IntegrationMethods(const IntegrationMethod&,FunctionPart=_allFunction,real_t=0);
...
void add(const IntegrationMethod&, FunctionPart=_allFunction,real_t=0);
void push_back(const IntgMeth&);
\end{lstlisting}
\vspace{.2cm}
To handle more than 3 integration methods, use the {\cmd add} member function.\\

Besides, the class offers some accessors and print stuff:
\vspace{.1cm}
\begin{lstlisting}[deletekeywords={[3] x, y}]
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&);
\end{lstlisting}
\vspace{.2cm}
In BEM computation, the following relative distance between two elements ($E_i$, $E_j$) is used:
 $$ dr(E_i,E_j)=\frac{||C_i-C_j||}{\max(diam(E_i),diam(E_j))}\text{ where }C_i\text{ is the centroid of } E_i.$$
 So defining for instance
\vspace{.1cm}
\begin{lstlisting}[deletekeywords={[3] x, y}]
IntegrationMethods ims(SauterSchwabIM(4), 0., symmetrical_Gauss, 5, 1.,
                       symmetrical_Gauss, 3);
\end{lstlisting}
\vspace{.2cm}
the BEM computation algorithm will perform on function (all part):
$$
\begin{array}{ll}
	 \text{Sauter-Schwab with quadrature of order 4 on segment}& \text{if } dr(E_i,E_j)=0\\
	 \text{Symmetrical gauss quadrature of order 5 on triangle}& \text{if } 0< dr(E_i,E_j)<=1\\
	 \text{Symmetrical gauss quadrature of order 3 on triangle}&\text{if }  dr(E_i,E_j) >1
\end{array}
$$
\displayInfos{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}}