1%%%%%%%%%%%%%%%%%%% 2% XLiFE++ is an extended library of finite elements written in C++ 3% Copyright (C) 2014 Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin 4% 5% This program is free software: you can redistribute it and/or modify 6% it under the terms of the GNU General Public License as published by 7% the Free Software Foundation, either version 3 of the License, or 8% (at your option) any later version. 9% This program is distributed in the hope that it will be useful, 10% but WITHOUT ANY WARRANTY; without even the implied warranty of 11% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12% GNU General Public License for more details. 13% You should have received a copy of the GNU General Public License 14% along with this program. If not, see <http://www.gnu.org/licenses/>. 15%%%%%%%%%%%%%%%%%%% 16 17\section{Integration methods} 18 19To 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. \\ 20Among 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).\\ 21 22The inheritance diagram looks like\\ 23\footnotesize 24\begin{verbatim} 25IntegrationMethod --> | SingleIM (single intg) --> | QuadratureIM (based on Quadrature class) 26 (abstract) | (abstract) | PolynomialIM (exact intg. of polynoms) 27 | | LenoirSalles2dIR (for 2D Laplace kernel) 28 | | LenoirSalles3dIR (for 3D Laplace kernel) 29 | 30 | DoubleIM (double intg) --> | ProductIM (product of single intg meth.) 31 (abstract) | LenoirSalles2dIM (for 2D Laplace kernel) 32 | LenoirSalles3dIM (for 3D Laplace kernel) 33 | SauterSchwabIM (for 3D IE singular integral) 34 | DuffyIM (for 2D IE singular integral) 35\end{verbatim} 36\normalsize 37 38As 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). 39 40 41\subsection{The {\classtitle IntegrationMethod, SingleIM, DoubleIM} classes} 42 43\subsubsection*{The {\classtitle IntegrationMethod} class} 44 45The {\class IntegrationMethod} class, basis class of all integration methods, manages a name, a type of integral and two booleans:\\ 46\vspace{.1cm} 47\begin{lstlisting}[deletekeywords={[3] name}] 48class IntegrationMethod 49{public : 50 String name; 51 IntegrationMethodType imType; 52 bool requireRefElement; 53 bool requirePhyElement; 54}; 55\end{lstlisting} 56\vspace{.2cm} 57Up to now, the types of integral defined in enumeration are: 58\vspace{.1cm} 59\begin{lstlisting}[]{} 60enum IntegrationMethodType 61 {_undefIM, _quadratureIM, _polynomialIM, _productIM, 62 _LenoirSalles2dIM, Lenoir_Salles_2d =_LenoirSalles2dIM, 63 _LenoirSalles3dIM, Lenoir_Salles_3d= _LenoirSalles3dIM, 64 _LenoirSalles2dIR, _LenoirSalles3dIR, 65 _SauterSchwabIM, Sauter_Schwab=_SauterSchwabIM, 66 _DuffyIM, Duffy = _DuffyIM, 67 _HMatrixIM, H_Matrix= _HMatrixIM 68 }; 69\end{lstlisting} 70\vspace{.2cm} 71The 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!\\ 72 73This class provides few general functions, most of them being virtual ones 74\vspace{.1cm} 75\begin{lstlisting}[]{} 76IntegrationMethod(IntegrationMethodType imt=_undefIM, bool ref=false, bool phy=false, const String& na="") 77virtual ~IntegrationMethod(){}; 78virtual bool isSingleIM() const =0; 79virtual bool isDoubleIM() const =0; 80virtual void print(std::ostream& os) const; 81IntegrationMethodType type() const; 82virtual bool useQuadraturePoints() const; 83 84std::ostream& operator<<(std::ostream&, const IntegrationMethod&); 85\end{lstlisting} 86 87\subsubsection*{The {\classtitle SingleIM, DoubleIM} classes} 88 89The {\class SingleIM, DoubleIM} derived classes provides only simple member functions : 90\vspace{.1cm} 91\begin{lstlisting}[]{} 92class SingleIM : public IntegrationMethod 93{public : 94 bool isSingleIM() const {return true;} 95 bool isDoubleIM() const {return false;} 96 SingleIM(IntegrationMethodType imt=_undefIM); 97 virtual void print(std::ostream& os) const; 98} 99\end{lstlisting} 100\vspace{.2cm} 101 102\begin{lstlisting}[]{} 103class DoubleIM : public IntegrationMethod 104{public : 105 bool isSingleIM() const {return false;} 106 bool isDoubleIM() const {return true;} 107 DoubleIM(IntegrationMethodType imt=_undefIM); 108 virtual void print(std::ostream& os) const; 109} 110\end{lstlisting} 111\vspace{.2cm} 112 113\subsection{The {\classtitle QuadratureIM} class}\label{quadrature_section} 114 115The {\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):\\ 116\vspace{.1cm} 117\begin{lstlisting}[deletekeywords={[3] map}] 118class QuadratureIM : public SingleIM 119{protected : 120 std::map<ShapeType, Quadrature *> quadratures_; 121 std::map<ShapeType, std::vector<ShapeValues>* > shapevalues_; 122\end{lstlisting} 123\vspace{.2cm} 124It provides two basic constructors, accessors and print facilities (all are public): 125\vspace{.1cm} 126\begin{lstlisting}[deletekeywords={[3] set}] 127QuadratureIM( ShapeType, QuadRule =_defaultRule, Number =0); 128QuadratureIM(const std::set<ShapeType>&, QuadRule =_defaultRule, Number =0); 129virtual bool useQuadraturePoints() const; 130Quadrature* getQuadrature(ShapeType) const; 131void setShapeValues(ShapeType, std::vector<ShapeValues>*); 132std::vector<ShapeValues>* getShapeValues(ShapeType); 133virtual void print(std::ostream& os) const; 134std::ostream& operator<<(std::ostream&, const QuadratureIM&); 135\end{lstlisting} 136\vspace{.2cm} 137 138\subsection{The {\classtitle ProductIM} class} 139 140In order to represent a product of integration method over a product of geometric domain, the {\class ProductIM} class carries two {\class SingleIM} pointers: 141\vspace{.1cm} 142\begin{lstlisting}[]{} 143class ProductIM : public DoubleIM 144{protected : 145 SingleIM* im_x; //!< integration method along x 146 SingleIM* im_y; //!< integration method along y 147} 148\end{lstlisting} 149\vspace{.2cm} 150and provides simple member functions: 151\vspace{.1cm} 152\begin{lstlisting}[]{} 153 154 ProductIM(SingleIM* imx=0, SingleIM* imy=0); 155 SingleIM* getxIM(); 156 SingleIM* getyIM(); 157 virtual bool useQuadraturePoints() const; 158 virtual void print(std::ostream& os) const; //!< print on stream 159}; 160\end{lstlisting} 161\vspace{.2cm} 162 163\subsection{The {\classtitle Quadrature, QuadratureRule} classes} 164 165The quadrature formulae have the following form : 166$$\int_{\widehat{E}} f(\widehat{x})d\widehat{x}\approx \sum_{i=1,q} \omega_i\,f(\widehat{x}_i)$$ 167where $(\widehat{x}_i)_{i=1,q}$ are quadrature points belonging to reference element $\hat{E}$ and $(w_i)_{i=1,q}$ are quadrature weights. \\ 168Up 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 :\\ 169\goodbreak 170 171\begin{center} 172\textbf{General rules }\\ 173\vspace{2mm} 174\small 175\begin{tabular}{|c|c|c|c|c|} 176 \hline 177 & Gauss-Legendre & Gauss-Lobatto & Grundmann-Muller & symmetrical Gauss\\ 178 \hline 179 segment & any odd degree & any odd degree & & \\ 180 \hline 181 quadrangle & any odd degree & any odd degree & & odd degree up to 21 \\ 182 \hline 183 triangle & any odd degree & & any odd degree & degree up to 10\\ 184 \hline 185 hexahedron & any odd degree & any odd degree & & odd degree up to 11\\ 186 \hline 187 tetrahedron & any odd degree & & any odd degree & degree up to 10\\ 188 \hline 189 prism & & & & degree up to 10\\ 190 \hline 191 pyramid & any odd degree & any odd degree & & degree up to 10 \\ 192 \hline 193\end{tabular} 194\end{center} 195 196 197\begin{center} 198\textbf{Particular rules}\\ 199\vspace{2mm} 200\small 201\begin{tabular}{|c|c|c|} 202 \hline 203 & nodal & miscellaneous\\ 204 \hline 205 segment & P1 to P4 & \\ 206 \hline 207 quadrangle & Q1 to Q4 & \\ 208 \hline 209 triangle & P1 to P3 & Hammer-Stroud 1 to 6\\ 210 \hline 211 hexahedron & Q1 to Q4 & \\ 212 \hline 213 tetrahedron & P1, P3 & Stroud 1 to 5\\ 214 \hline 215 prism & P1 & centroid 1, tensor product 1,3,5\\ 216 \hline 217 pyramid & P1 & centroid 1, Stroud 7 \\ 218 \hline 219\end{tabular} 220\end{center} 221More precisely, the rules that are implemented : 222\begin{itemize} 223\item 1D rules \\ 224 Gauss-Legendre rule, degree $2n+1$, $n$ points\\ 225 Gauss-Lobatto rule, degree $2n-1$, $n$ points\\ 226 trapezoidal rule, degree $1$, $2$ points\\ 227 Simpson rule, degree $3$, $3$ points\\ 228 Simpson $3$ $8$ rule, degree $3$, $4$ points\\ 229 Booleb rule, degree $5$, $5$ points\\ 230 Wedge rule, degree $5$, $6$ points 231 232\item Rules over the unit simplex (less quadrature points BUT with negative weights)\\ 233 TN Grundmann-Moller rule, degree $2n+1$, $C^n_{n+d+1}$ points 234 235\item Rules over the unit triangle $\{ x > 0, y >0, x+y < 1 \}$\\ 236 T2P2 MidEdge rule\\ 237 T2P2 Hammer-Stroud rule\\ 238 T2P3 Albrecht-Collatz rule, degree $3$, $6$ points, Stroud (p.314)\\ 239 T2P3 Stroud rule, degree $3$, $7$ points, Stroud (p.308)\\ 240 T2P5 Radon-Hammer-Marlowe-Stroud rule, degree $6$, $7$ points, Stroud (p.314)\\ 241 T2P6 Hammer rule, degree $6$, $12$ points, G.Dhatt, G.Touzot (p.298)\\ 242 symmetrical Gauss up to degree 10 243 244\item Rules over the unit tetrahedron $\{x > 0, y >0, z>0, x+y+z < 1\}$\\ 245 T3P2 Hammer-Stroud rule, degree 2, 4 points, Stroud (p.307)\\ 246 T3P3 Stroud rule, degree $3$, $8$ points, Stroud (p.308)\\ 247 T3P5 Stroud rule, degree $5$, $15$ points, Stroud (p.315)\\ 248 symmetrical Gauss up to degree 10 249 250 \item Rules over the unit prism $\{ 0 < x, y < 1, x+y<1, 0 < z < 1\}$ \\ 251 P1 nodal rule \\ 252 tensor product of 2-points Gauss-Legendre rule on [0,1] and Stroud 7 points formula on the unit triangle (degree 3)\\ 253 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)\\ 254 symmetrical Gauss up to degree 10 255 256\item Rules over the unit pyramid $\{ 0 < x, y < 1-z, 0 < z < 1\}$ \\ 257P1 nodal rule \\ 258symmetrical Gauss up to degree 10 259 260 261\item Rules built from lower dimension rules \\ 262 tensor product of quadrature rules (quadrangle and hexahedron)\\ 263 conical product of quadrature rules (triangle and pyramid)\\ 264 rule built from 1D nodal rule which are quadrangle nodal rule\\ 265 rule built from 1D nodal rule which are hexahedron nodal rule 266\end{itemize} 267 268References : 269\begin{itemize} 270\item A.H.Stroud, Approximate calculation of multiple integrals, Prentice Hall, 1971. 271\item G.Dhatt \& G Touzot, The finite element method displayed, John Wiley \& Sons, Chichester, 1984 272\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. 273\item F.D. Witherden \& P.E. Vincent, On the identification of symmetric quadrature rules for finite element methods, 274 Computers \& Mathematics with Applications, Volume 69, Issue 10, May 2015, Pages 1232-1241. 275\end{itemize} 276 277These 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. 278 279\subsection*{The {\classtitle Quadrature} class} 280 281The {\class Quadrature} class manages all informations related to quadrature : 282\vspace{.1cm} 283\begin{lstlisting} 284class Quadrature 285{public: 286 GeomRefElement* geomRefElt_p; // geometric element bearing quadrature rule 287 QuadratureRule quadratureRule; // embedded QuadratureRule object 288 protected: 289 QuadRule rule_; // type of quadrature rule 290 Number degree_; // degree of quadrature rule 291 bool hasPointsOnBoundary_;// as it reads 292 String name_; // name of quadrature formula 293 public: 294 static std::vector<Quadrature*> theQuadratureRules; 295 ... 296\end{lstlisting} 297\vspace{.2cm} 298\verb?QuadRule? is a enumeration of all types of quadrature supported: 299\vspace{.1cm} 300\begin{lstlisting}[] 301enum QuadRule 302{_defaultRule = 0, _GaussLegendreRule, _symmetricalGaussRule, _GaussLobattoRule, Gauss_Lobatto = _GaussLobattoRule, 303 _nodalRule, _miscRule,_GrundmannMollerRule, _doubleQuadrature }; 304\end{lstlisting} 305\vspace{.1cm} 306\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. \\ 307 308The class proposes only a default constructor, a full basic constructor and a destructor : 309\vspace{.1cm} 310\begin{lstlisting}[] 311Quadrature(); 312Quadrature(ShapeType , QuadRule, Number, const String&, bool pob = false); 313\end{lstlisting} 314\vspace{.1cm} 315Note that the basic constructor initializes member attributes but does not load the quadrature points and quadrature weights ({\class QuadratureRule})! 316The 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\\ 317 318 319The following accessors (read only) are provided: 320\vspace{.1cm} 321\begin{lstlisting}[] 322String name() const; 323QuadRule rule() const; 324Number degree() const; 325bool hasPointsOnBoundary() const; 326Dimen dim() const; 327Number numberOfPoints() const; 328ShapeType shapeType() const; 329std::vector<Real>::const_iterator point(Number i = 0) const; 330std::vector<Real>::const_iterator weight(Number i = 0) const; 331\end{lstlisting} 332\vspace{.2cm} 333\textdbend \ \ To construct a {\class Quadrature} object, the external function \verb?findQuadrature? has to be used : 334\vspace{.1cm} 335\begin{lstlisting}[] 336friend Quadrature* findQuadrature(ShapeType shape, QuadRule rule, Number degree, bool hpob= false); 337\end{lstlisting} 338\vspace{.1cm} 339It 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} 340object 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 : 341\vspace{.1cm} 342\begin{lstlisting}[] 343Quadrature* segmentQuadrature(QuadRule, Number); 344Quadrature* triangleQuadrature(QuadRule, Number); 345Quadrature* quadrangleQuadrature(QuadRule, Number); 346Quadrature* tetrahedronQuadrature(QuadRule, Number); 347Quadrature* hexahedronQuadrature(QuadRule, Number); 348Quadrature* prismQuadrature(QuadRule, Number); 349Quadrature* pyramidQuadrature(QuadRule, Number); 350static void clear(); 351\end{lstlisting} 352\vspace{.1cm} 353These functions loads quadrature points and quadrature weights in a {\class QuadratureRule} object using member functions of {\class QuadratureRule} class. 354The \verb?clear? static function deletes quadrature objects stored in the \verb?theQuadratureRules? static vector.\\ 355 356To 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$. 357\vspace{.1cm} 358\begin{lstlisting}[] 359static QuadRule bestQuadRule(ShapeType, Number); 360\end{lstlisting} 361\vspace{.1cm} 362Best rule has to be understood as the rule on shape $S$ with the minimum of quadrature points integrating exactly polynomials of degree $d$. 363The following table gives the current best rules : 364\begin{center} 365\begin{tabular}{|c|c|c|c|} 366\hline 367shape& degree & quadrature rule & number of points\\ 368\hline 369segment& $1$& Gauss-Legendre $1$& $1$\\ 370segment& $2,3$& Gauss-Legendre $3$& $2$\\ 371segment& $4,5$& Gauss-Legendre $5$& $3$\\ 372segment& $2n-1$& Gauss-Legendre $2n-1$& $n$\\ 373\hline 374\end{tabular} 375\end{center} 376 377\begin{center} 378\begin{tabular}{|c|c|c|c|} 379\hline 380shape& degree & quadrature rule & number of points\\ 381\hline 382quadrangle& $1$& Gauss-Legendre $1$& $1$\\ 383quadrangle& $2,3$& Gauss-Legendre $3$& $4$\\ 384quadrangle& $4,5$& symmetrical Gauss $5$& $8$\\ 385quadrangle& $6,7$& symmetrical Gauss $7$& $12$\\ 386quadrangle& $7,9$& symmetrical Gauss $9$& $20$\\ 387quadrangle& $10,11$& symmetrical Gauss $11$& $28$\\ 388quadrangle& $12,13$& symmetrical Gauss $13$& $37$\\ 389quadrangle& $14,15$& symmetrical Gauss $15$& $48$\\ 390quadrangle& $16,17$& symmetrical Gauss $17$& $60$\\ 391quadrangle& $18,19$& symmetrical Gauss $19$& $72$\\ 392quadrangle& $20,21$& symmetrical Gauss $21$& $85$\\ 393quadrangle& $2n-1>21$& Gauss-Legendre $2n-1$& $n^2$\\ 394\hline 395\end{tabular} 396\end{center} 397\begin{center} 398\begin{tabular}{|c|c|c|c|} 399\hline 400shape& degree & quadrature rule & number of points\\ 401\hline 402triangle& $1$& centroid rule (misc,$1$)& $1$\\ 403triangle& $2$& P2 Hammer-Stroud (misc,$2$)& $2$\\ 404triangle& $3$& Grundmann-Moller $3$ & $3$\\ 405triangle& $4$& symmetrical Gauss $4$ & $6$\\ 406triangle& $5$& symmetrical Gauss $5$ & $7$\\ 407triangle& $6$& symmetrical Gauss $6$ & $12$\\ 408triangle& $7$& symmetrical Gauss $7$ & $15$\\ 409triangle& $8$& symmetrical Gauss $8$ & $16$\\ 410triangle& $9$& symmetrical Gauss $9$ & $19$\\ 411triangle& $10$& symmetrical Gauss $10$ & $25$\\ 412triangle& $11$& symmetrical Gauss $11$ & $28$\\ 413triangle& $12$& symmetrical Gauss $12$ & $33$\\ 414triangle& $13$& symmetrical Gauss $13$ & $37$\\ 415triangle& $14$& symmetrical Gauss $14$ & $42$\\ 416triangle& $15$& symmetrical Gauss $15$ & $49$\\ 417triangle& $16$& symmetrical Gauss $16$ & $55$\\ 418triangle& $17$& symmetrical Gauss $17$ & $60$\\ 419triangle& $18$& symmetrical Gauss $18$ & $67$\\ 420triangle& $19$& symmetrical Gauss $19$ & $73$\\ 421triangle& $20$& symmetrical Gauss $20$ & $79$\\ 422triangle& $2n-1>20$& Gauss-Legendre $2n-1$& $?$\\ 423\hline 424\end{tabular} 425\end{center} 426\begin{center} 427\begin{tabular}{|c|c|c|c|} 428\hline 429shape& degree & quadrature rule & number of points\\ 430\hline 431hexahedron& $1$& Gauss-Legendre $1$& $1$\\ 432hexahedron& $2,3$& symmetrical Gauss $3$& $6$\\ 433hexahedron& $4,5$& symmetrical Gauss $5$& $14$\\ 434hexahedron& $6,7$& symmetrical Gauss $7$& $34$\\ 435hexahedron& $8,9$& symmetrical Gauss $9$& $58$\\ 436hexahedron& $10,11$& symmetrical Gauss $11$& $90$\\ 437hexahedron& $2n-1\ge 12$& Gauss-Legendre $2n-1$& $n^3$\\ 438\hline 439\end{tabular} 440\end{center} 441\begin{center} 442\begin{tabular}{|c|c|c|c|} 443\hline 444shape& degree & quadrature rule & number of points\\ 445\hline 446tetrahedron& $1$& centroid rule (misc,$1$)& $1$\\ 447tetrahedron& $2$& P2 Hammer-Stroud (misc,$2$)& $4$\\ 448tetrahedron& $3$& Grundmann-Moller $3$ & $5$\\ 449tetrahedron& $4,5$& symmetrical Gauss $5$ & $14$\\ 450tetrahedron& $6$& symmetrical Gauss $6$ & $24$\\ 451tetrahedron& $7$& symmetrical Gauss $7$ & $35$\\ 452tetrahedron& $8$& symmetrical Gauss $8$ & $46$\\ 453tetrahedron& $9$& symmetrical Gauss $9$ & $59$\\ 454tetrahedron& $10$& symmetrical Gauss $10$ & $81$\\ 455tetrahedron& $2n-1\ge 11$& Grundmann-Moller $2n-1$ & $?$\\ 456\hline 457\end{tabular} 458\end{center} 459\begin{center} 460\begin{tabular}{|c|c|c|c|} 461\hline 462shape& degree & quadrature rule & number of points\\ 463\hline 464prism& $1$& centroid rule (misc,$1$)& $1$\\ 465prism& $2$& symmetrical Gauss $2$ & $5$\\ 466prism& $3$& symmetrical Gauss $3$ & $8$\\ 467prism& $4$& symmetrical Gauss $4$ & $11$\\ 468prism& $5$& symmetrical Gauss $5$ & $16$\\ 469prism& $6$& symmetrical Gauss $6$ & $28$\\ 470prism& $7$& symmetrical Gauss $7$ & $35$\\ 471prism& $8$& symmetrical Gauss $8$ & $46$\\ 472prism& $9$& symmetrical Gauss $9$ & $60$\\ 473prism& $10$& symmetrical Gauss $10$ & $85$\\ 474\hline 475\end{tabular} 476\end{center} 477\begin{center} 478\begin{tabular}{|c|c|c|c|} 479\hline 480shape& degree & quadrature rule & number of points\\ 481\hline 482pyramid& $1$ & centroid rule (misc,$1$)& $1$\\ 483pyramid& $2$& symmetrical Gauss $2$ & $5$\\ 484pyramid& $3$& symmetrical Gauss $3$ & $6$\\ 485pyramid& $4$& symmetrical Gauss $4$ & $10$\\ 486pyramid& $5$& symmetrical Gauss $5$ & $15$\\ 487pyramid& $6$& symmetrical Gauss $6$ & $24$\\ 488pyramid& $7$& symmetrical Gauss $7$ & $31$\\ 489pyramid& $8$& symmetrical Gauss $8$ & $47$\\ 490pyramid& $9$& symmetrical Gauss $9$ & $62$\\ 491pyramid& $10$& symmetrical Gauss $10$ & $83$\\ 492\hline 493\end{tabular} 494\end{center} 495\mbox{} \\ 496Finally, {\class Quadrature} class provides some print and error utilities: 497\vspace{.1cm} 498\begin{lstlisting}[] 499void badNodeRule(int n_nodes); 500void badDegreeRule(); 501void noEvenDegreeRule(); 502void upperDegreeRule(); 503friend std::ostream& operator<<(std::ostream&, const Quadrature&); 504friend void alternateRule(QuadRule, ShapeType, const String&); 505\end{lstlisting} 506\vspace{.1cm} 507 508 509\subsection*{The {\classtitle QuadratureRule} class} 510 511The {\class QuadratureRule} class stores quadrature points and quadrature weights and provides member functions to load them : 512\vspace{.1cm} 513\begin{lstlisting} 514class QuadratureRule 515{private: 516 std::vector<Real> coords_; // point coordinates of quadrature rule 517 std::vector<Real> weights_; // weights of quadrature rule 518 Dimen dim_; // dimension of points 519... 520\end{lstlisting} 521\vspace{.2cm} 522the following constructors are proposed: 523\vspace{.1cm} 524\begin{lstlisting}[] 525QuadratureRule(const std::vector<Real>&, Real); 526QuadratureRule(const std::vector<Real>&, const std::vector<Real>&); 527QuadratureRule(Dimen , Number s); 528\end{lstlisting} 529\vspace{.2cm} 530and the following accessors are provided : 531\vspace{.1cm} 532\begin{lstlisting}[] 533Dimen dim() const; 534Number size() const; 535std::vector<Real> coords() const; 536std::vector<Real>::const_iterator point(Number i = 0) const; 537std::vector<Real>::iterator point(Number i = 0); 538std::vector<Real>::const_iterator weight(Number i = 0) const; 539std::vector<Real>::iterator weight(const Number i = 0); 540\end{lstlisting} 541\vspace{.2cm} 542Some utilities are also helpful: 543\vspace{.1cm} 544\begin{lstlisting}[] 545void point(std::vector<Real>::iterator&, Real, std::vector<Real>::iterator&, Real); 546void point(std::vector<Real>::iterator&, Real, Real, std::vector<Real>::iterator&, Real); 547void point(std::vector<Real>::iterator&, Real, Real, Real,std::vector<Real>::iterator&, Real); 548void coords(std::vector<Real>::iterator&, Dimen, std::vector<Real>::const_iterator&); 549void coords(std::vector<Real>::const_iterator); 550void coords(const std::vector<Real>&); 551void coords(Real); 552void weights(const std::vector<Real>&); 553void weights(Real); 554friend std::ostream& operator<<(std::ostream&, const QuadratureRule&); 555\end{lstlisting} 556\vspace{.2cm} 557The 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: 558\vspace{.1cm} 559\begin{lstlisting}[] 560//tensor and conical product 561 QuadratureRule& tensorRule(const QuadratureRule&, const QuadratureRule&); 562 QuadratureRule& tensorRule(const QuadratureRule&, const QuadratureRule&, const QuadratureRule&); 563 QuadratureRule& conicalRule(const QuadratureRule&, const QuadratureRule&); 564 QuadratureRule& conicalRule(const QuadratureRule&, const QuadratureRule&, const QuadratureRule&); 565//nodal tensor product rules for quadrangle and hexahedron 566 QuadratureRule& quadrangleNodalRule(const QuadratureRule&); 567 QuadratureRule& hexahedronNodalRule(const QuadratureRule&); 568//1D Gauss rules 569 QuadratureRule& gaussLegendreRules(Number); 570 QuadratureRule& gaussLobattoRules(Number); 571 QuadratureRule& ruleOn01(Number); 572//1D Newton-Cotes closed rules (P_k node rules on [0,1]) 573 QuadratureRule& trapezoidalRule(); 574 QuadratureRule& simpsonRule(); 575 QuadratureRule& simpson38Rule(); 576 QuadratureRule& booleRule(); 577//Grundmann-Moller rules for n dimensional simplex 578 QuadratureRule& tNGrundmannMollerRule(int, Dimen); 579 void TNGrundmannMollerSet(int); //utility for Grundmann-Moller 580//special (Misc) rules for unit Triangle { x + y < 1 , 0 < x,y < 1} 581 QuadratureRule& t2P2MidEdgeRule(); 582 QuadratureRule& t2P2HammerStroudRule(); 583 QuadratureRule& t2P3AlbrechtCollatzRule(); 584 QuadratureRule& t2P3StroudRule(); 585 QuadratureRule& t2P5RadonHammerMarloweStroudRule(); 586 QuadratureRule& t2P6HammerRule(); 587//special (Misc) rules for unit Tetrahedron { x + y + z < 1 , 0 < x,y,z < 1 } 588 QuadratureRule& t3P2HammerStroudRule(); 589 QuadratureRule& t3P3StroudRule(); 590 QuadratureRule& t3P5StroudRule(); 591//special (Misc) rules for unit pyramid { 0 < x,y < 1 - z , 0 < z < 1 } 592 QuadratureRule& pyramidRule(const QuadratureRule&);//conical product 593 QuadratureRule& pyramidStroudRule(); //degree 7, 48 points 594//symmetrical Gauss rules for any shape 595 QuadratureRule& symmetricalGaussTriangleRule(number_t); //degree up <= 20 596 QuadratureRule& symmetricalGaussQuadrangleRule(number_t); //odd degree <= 21 597 QuadratureRule& symmetricalGaussTetrahedronRule(number_t); //degree <= 10 598 QuadratureRule& symmetricalGaussHexahedronRule(number_t); //odd degree <= 11 599 QuadratureRule& symmetricalGaussPrismRule(number_t); //degre <= 10 600 QuadratureRule& symmetricalGaussPyramidRule(number_t); //degree <= 10 601 \end{lstlisting} 602\vspace{.2cm} 603Here is a simple example of how to use quadrature objects : 604\vspace{.1cm} 605\begin{lstlisting}[deletekeywords={[3] x, y}] 606Quadrature* q=findQuadrature(_triangle,_GaussLegendreRule, 3); 607Real S=0, x,y; 608Number k=0; 609for(Number i =0;i<q->quadratureRule.size();i++,k+=2) 610 { 611 x=q->quadratureRule.coords()[k]; 612 y=q->quadratureRule.coords()[k+1]; 613 S+=q->quadratureRule.weights()[i]*std::sin(x*y); 614 } 615std::cout<<"S="<<S; 616\end{lstlisting} 617\vspace{.2cm} 618 619 620\subsection{How to add a new quadrature formula ?} 621 622When you want to add a new quadrature formula on a unit element, say \emph{geomelt}, 623\begin{itemize} 624\item you have to modify the member function :\\ 625\ \ \ \verb?Quadrature::?\emph{geomelt}\verb?Quadrature(QuadRule rule, Number deg)? \\ 626by adding in the right \emph{case} statement the properties of your quadrature formula 627\item to add in {\class QuadratureRule} class the member function defining the quadrature points and weights. 628\end{itemize} 629Be care with the type and the number (degree) defining your quadrature formula, because they can be used by an other formula. 630To avoid this difficulty, it is possible to create new type of quadrature. In that case the enumeration \verb?QuadRule? has to be modified, 631implying to update the dictionnaries and may be, the \verb?Quadrature::bestQuadRule? static function.\\ 632 633\subsection{Integration methods for singular kernels} 634To deal with single or double integral involving a kernel with a singularity, \xlifepp provides some additional classes. 635 636\subsection*{{\class DuffyIM} class} 637This class deals with double integral over segments involving 2d kernel $K$ with $\log|x-y|$ singularity, in other words: 638 $$\int_{S_1}\int_{S_2} p_1(x)\,K(x,y)\,p_2(y)\,dy\,dx$$ 639\vspace{.1cm} 640\begin{lstlisting}[] 641class DuffyIM : public DoubleIM 642{ 643 public: 644 Quadrature* quadSelf; //quadrature for self influence 645 Quadrature* quadAdjt; //quadrature for adjacent elements 646 number_t ordSelf; //order of quadrature for self influence 647 number_t ordAdjt; //order of quadrature for adjacent elements 648 649 DuffyIM(IntegrationMethodType); 650 DuffyIM(number_t =6); //full constructor from quadrature order 651 DuffyIM(const Quadrature&); //full constructor from a quadrature 652\end{lstlisting} 653 654\subsection*{{\class SauterSchwabIM} class} 655This class deals with double integral over triangles involving 3d kernel $K$ with $|x-y|^{-1}$ singularity, say 656$$\int_{T_1}\int_{T_2} p_1(x)\,K(x,y)\,p_2(y)\,dy\,dx$$ 657\vspace{.1cm} 658\begin{lstlisting}[] 659class SauterSchwabIM : public DoubleIM 660{ 661 public: 662 Quadrature* quadSelf; //quadrature for self influence 663 Quadrature* quadEdge; //quadrature for elements adjacent by edge 664 Quadrature* quadVertex; //quadrature for elements adjacent by vertex 665 number_t ordSelf; //order of quadrature for self influence 666 number_t ordEdge; //order of quadrature for edge adjacence 667 number_t ordVertex; //order of quadrature for vertex adjacence 668 669 SauterSchwabIM(number_t =3); //basic constructor (order 3) 670 SauterSchwabIM(Quadrature&); //full constructor from quadrature object 671 \end{lstlisting} 672 673\subsection*{{\class LenoirSalles2dIM} and {\class LenoirSalles3dIM} classes} 674The {\class LenoirSalles2dIM} class can compute analytically 675 $$-\frac{1}{2\pi}\int_{S_1}\int_{S_2} p_1(x)\,\log|x-y|\,p_2(y)\,dy\,dx \ \ \ \text{(single layer)}$$ 676 $$\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)}$$ 677 where $p_1$ and $p_2$ are either the P0 or P1 shape functions on segments $S_1$ and $S_2$ 678 and the {\class LenoirSalles3dIM} class computes analytically 679 $$\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)}$$ 680 $$-\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)}$$ 681 where $p_1$ and $p_2$ are the P0 shape functions (say $1$) on triangles $T_1$ and $T_2$.\\ 682 These two classes have no member data and only one basic constructor. 683 684 \subsection*{{\class LenoirSalles2dIR} and {\class LenoirSalles3dIR} classes} 685 The {\class LenoirSalles2dIR} class can compute analytically 686 $$-\frac{1}{2\pi}\int_{S} \log|x-y|\,p(y)\,dy\,dx \ \ \ \text{(single layer)}$$ 687 $$\frac{1}{2\pi}\int_{S}\frac{(x-y).n_y}{|x-y|^2}\,p(y)\,dy\,dx \ \ \ \text{(double layer)}$$ 688 where $p$ is either the P0 or P1 shape functions on segment $S$ 689 and the {\class LenoirSalles3dIR} classe that computes analytically 690 $$\frac{1}{4\pi}\int_{T}\frac{1}{|x-y|}\,p(y)\,dy\,dx \ \ \ \text{(single layer)}$$ 691 $$-\frac{1}{4\pi}\int_{TS}\frac{(x-y).n_y}{|x-y|^3}\,p(y)\,dy\,dx \ \ \ \text{(double layer)}$$ 692 where $p$ is the P0 or P1 shape functions on triangle $T$.\\ 693 These two classes have no member data and only one basic constructor. 694\\ 695\\ 696All 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: 697\vspace{.1cm} 698\begin{lstlisting}[] 699virtual LenoirSalles2dIM* clone() const; //covariant 700virtual std::list<Quadrature*> quadratures() const; //list of quadratures 701template<typename K> 702void computeIE(const Element* S, const Element* T, AdjacenceInfo& adj, 703 const KernelOperatorOnUnknowns& kuv, Matrix<K>& res, 704 IEcomputationParameters& iep, 705 Vector<K>& vopu, Vector<K>& vopv, Vector<K>& vopk) const; 706\end{lstlisting} 707 708\subsection{The {\classtitle IntegrationMethods} class} 709The {\class IntegrationMethods} class collects some {\class IntgMeth} objects that handle an integration method pointer with some additional parameters: 710\vspace{.1cm} 711\begin{lstlisting}[deletekeywords={[3] x, y}] 712class IntgMeth 713{ public: 714 const IntegrationMethod* intgMeth; 715 FunctionPart functionPart; //_allFunction, _regularPart, _singularPart 716 real_t bound; //bound value 717} 718\end{lstlisting} 719\vspace{.2cm} 720The {\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 721only some constructors and print stuff: 722\vspace{.1cm} 723\begin{lstlisting}[deletekeywords={[3] x, y}] 724IntgMeth(const IntegrationMethod& im, FunctionPart fp=_allFunction, real_t b=0); 725IntgMeth(const IntgMeth&); 726~IntgMeth(); 727IntgMeth& operator=(const IntgMeth&); 728void print(std::ostream&) const; 729void print(PrintStream&) const; 730friend ostream& operator<<(ostream&, const IntgMeth&) 731\end{lstlisting} 732\vspace{.2cm} 733For safe multithreading behaviour, copy constructor and assign operator do a full copy of {\class IntegrationMethod} object.\\ 734 735The {\class IntegrationMethods} class is simply a vector of {\class IntMeth} object: 736\vspace{.1cm} 737\begin{lstlisting}[deletekeywords={[3] x, y}] 738class IntegrationMethods 739{public: 740 vector<IntgMeth> intgMethods; 741 typedef vector<IntgMeth>::const_iterator const_iterator; 742 typedef vector<IntgMeth>::iterator iterator; 743 ...} 744\end{lstlisting} 745\vspace{.2cm} 746This class provides many constructors depending to the number of integration methods given (up to 3 at this time), to some shortcuts : 747\vspace{.1cm} 748\begin{lstlisting}[deletekeywords={[3] x, y}] 749IntegrationMethods() {}; 750IntegrationMethods(const IntegrationMethod&,FunctionPart=_allFunction,real_t=0); 751... 752void add(const IntegrationMethod&, FunctionPart=_allFunction,real_t=0); 753void push_back(const IntgMeth&); 754\end{lstlisting} 755\vspace{.2cm} 756To handle more than 3 integration methods, use the {\cmd add} member function.\\ 757 758Besides, the class offers some accessors and print stuff: 759\vspace{.1cm} 760\begin{lstlisting}[deletekeywords={[3] x, y}] 761const IntgMeth& operator[](number_t) cons ; 762vector<IntgMeth>::const_iterator begin() const; 763vector<IntgMeth>::const_iterator end() const; 764vector<IntgMeth>::iterator begin(); 765vector<IntgMeth>::iterator end(); 766bool empty() const {return intgMethods.empty();} 767void print(ostream&) const; 768void print(PrintStream&) const; 769friend ostream& operator<<(ostream&, const IntegrationMethods&); 770\end{lstlisting} 771\vspace{.2cm} 772In BEM computation, the following relative distance between two elements ($E_i$, $E_j$) is used: 773 $$ 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.$$ 774 So defining for instance 775\vspace{.1cm} 776\begin{lstlisting}[deletekeywords={[3] x, y}] 777IntegrationMethods ims(SauterSchwabIM(4), 0., symmetrical_Gauss, 5, 1., 778 symmetrical_Gauss, 3); 779\end{lstlisting} 780\vspace{.2cm} 781the BEM computation algorithm will perform on function (all part): 782$$ 783\begin{array}{ll} 784 \text{Sauter-Schwab with quadrature of order 4 on segment}& \text{if } dr(E_i,E_j)=0\\ 785 \text{Symmetrical gauss quadrature of order 5 on triangle}& \text{if } 0< dr(E_i,E_j)<=1\\ 786 \text{Symmetrical gauss quadrature of order 3 on triangle}&\text{if } dr(E_i,E_j) >1 787\end{array} 788$$ 789\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}} 790