1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- 2 // vi: set et ts=4 sw=2 sts=2: 3 #ifndef DUNE_GEOMETRY_QUADRATURE_SIMPLEX_HH 4 #define DUNE_GEOMETRY_QUADRATURE_SIMPLEX_HH 5 6 #ifndef DUNE_INCLUDING_IMPLEMENTATION 7 #error This is a private header that should not be included directly. 8 #error Use #include <dune/geometry/quadraturerules.hh> instead. 9 #endif 10 11 namespace Dune { 12 13 /************************************************ 14 * Quadraturerule for Simplices/Triangle 15 *************************************************/ 16 17 /** \brief Quadrature rules for simplices 18 \ingroup Quadrature 19 */ 20 template<typename ct, int dim> 21 class SimplexQuadratureRule; 22 23 /** \brief Quadrature rules for triangles 24 \ingroup Quadrature 25 */ 26 template<typename ct> 27 class SimplexQuadratureRule<ct,2> : public QuadratureRule<ct,2> 28 { 29 public: 30 /** \brief The highest quadrature order available */ 31 enum { highest_order = 12 }; 32 private: 33 friend class QuadratureRuleFactory<ct,2>; 34 SimplexQuadratureRule (int p); ~SimplexQuadratureRule()35 ~SimplexQuadratureRule(){} 36 }; 37 38 /** \brief Quadrature rules for tetrahedra 39 \ingroup Quadrature 40 */ 41 template<typename ct> 42 class SimplexQuadratureRule<ct,3> : public QuadratureRule<ct,3> 43 { 44 public: 45 /** \brief The highest quadrature order available */ 46 enum { highest_order = 5 }; 47 private: 48 friend class QuadratureRuleFactory<ct,3>; 49 SimplexQuadratureRule (int p); ~SimplexQuadratureRule()50 ~SimplexQuadratureRule(){} 51 }; 52 53 //! 54 template<int dim> 55 class SimplexQuadraturePoints; 56 57 /** Singleton holding the Gauss points on the interval */ 58 template<int dim> 59 struct SimplexQuadraturePointsSingleton {}; 60 61 template<> 62 struct SimplexQuadraturePointsSingleton<2> { 63 static SimplexQuadraturePoints<2> sqp; 64 }; 65 66 template<> 67 struct SimplexQuadraturePointsSingleton<3> { 68 static SimplexQuadraturePoints<3> sqp; 69 }; 70 71 template<> 72 class SimplexQuadraturePoints<2> 73 { 74 public: 75 enum { MAXP=33}; 76 enum { highest_order=12 }; 77 78 //! initialize quadrature points on the interval for all orders SimplexQuadraturePoints()79 SimplexQuadraturePoints () 80 { 81 init(); 82 } 83 init()84 void init() 85 { 86 int m = 0; 87 O[m] = 0; 88 89 // polynom degree 1 90 91 // Rule t2-1-1 of the "Encyclopaedia of Cubature Formulas" at 92 // http://www.cs.kuleuven.ac.be/~nines/research/ecf/ecf.html 93 // maintained by Ronald Cools. 94 95 // For further reference: Rule 1-1, P. 307, A.H. Stroud, Approximate Calculation of Multiple Integrals 96 97 m = 1; 98 G[m][0][0] = 0.333333333333333333333333333333333; 99 G[m][0][1] = 0.333333333333333333333333333333333; 100 W[m][0] = 0.5; 101 O[m] = 1; 102 103 // polynom degree 2 104 // symmetric 105 106 // Rule t2-2-3a of the "Encyclopaedia of Cubature Formulas" at 107 // http://www.cs.kuleuven.ac.be/~nines/research/ecf/ecf.html 108 // maintained by Ronald Cools. 109 110 // For further reference: Rule 2-1, P. 307, A.H. Stroud, Approximate Calculation of Multiple Integrals 111 112 m = 3; 113 G[m][0][0] = 4.0/6.0; 114 G[m][0][1] = 1.0/6.0; 115 G[m][1][0] = 1.0/6.0; 116 G[m][1][1] = 4.0/6.0; 117 G[m][2][0] = 1.0/6.0; 118 G[m][2][1] = 1.0/6.0; 119 W[m][0] = 0.5/3.0; 120 W[m][1] = 0.5/3.0; 121 W[m][2] = 0.5/3.0; 122 O[m] = 2; 123 124 // polynom degree 3 125 // symmetric 126 127 // Rule t2-3-4a of the "Encyclopaedia of Cubature Formulas" at 128 // http://www.cs.kuleuven.ac.be/~nines/research/ecf/ecf.html 129 // maintained by Ronald Cools. 130 131 // For further reference: Rule 3-1, P. 308, A.H. Stroud, Approximate Calculation of Multiple Integrals 132 133 m = 4; 134 G[m][0][0] = 10.0/30.0; 135 G[m][0][1] = 10.0/30.0; 136 G[m][1][0] = 18.0/30.0; 137 G[m][1][1] = 6.0/30.0; 138 G[m][2][0] = 6.0/30.0; 139 G[m][2][1] = 18.0/30.0; 140 G[m][3][0] = 6.0/30.0; 141 G[m][3][1] = 6.0/30.0; 142 143 W[m][0] = 0.5 * -27.0/48.0; 144 W[m][1] = 0.5 * 25.0/48.0; 145 W[m][2] = 0.5 * 25.0/48.0; 146 W[m][3] = 0.5 * 25.0/48.0; 147 O[m] = 3; 148 149 // polynomial degree 4 150 // symmetric points 151 152 // Rule t2-4-6a of the "Encyclopaedia of Cubature Formulas" at 153 // http://www.cs.kuleuven.ac.be/~nines/research/ecf/ecf.html 154 // maintained by Ronald Cools. 155 156 // For further reference: Appendix II, D.A. Dunavant, High degree efficient symmetrical Gaussian quadrature rules for the triangle 157 158 m = 6; 159 G[m][0][0] = 0.81684757298045851308085707319560; 160 G[m][0][1] = 0.091576213509770743459571463402202; 161 G[m][1][0] = 0.091576213509770743459571463402202; 162 G[m][1][1] = 0.81684757298045851308085707319560; 163 G[m][2][0] = 0.091576213509770743459571463402202; 164 G[m][2][1] = 0.091576213509770743459571463402202; 165 G[m][3][0] = 0.10810301816807022736334149223390; 166 G[m][3][1] = 0.44594849091596488631832925388305; 167 G[m][4][0] = 0.44594849091596488631832925388305; 168 G[m][4][1] = 0.10810301816807022736334149223390; 169 G[m][5][0] = 0.44594849091596488631832925388305; 170 G[m][5][1] = 0.44594849091596488631832925388305; 171 172 W[m][0] = 0.5 * 0.10995174365532186763832632490021; 173 W[m][1] = 0.5 * 0.10995174365532186763832632490021; 174 W[m][2] = 0.5 * 0.10995174365532186763832632490021; 175 W[m][3] = 0.5 * 0.22338158967801146569500700843312; 176 W[m][4] = 0.5 * 0.22338158967801146569500700843312; 177 W[m][5] = 0.5 * 0.22338158967801146569500700843312; 178 O[m] = 4; 179 180 // polynomial degree 5 181 // symmetric points 182 183 // Rule t2-5-7 of the "Encyclopaedia of Cubature Formulas" at 184 // http://www.cs.kuleuven.ac.be/~nines/research/ecf/ecf.html 185 // maintained by Ronald Cools. 186 187 // For further reference: Rule 5-1, P. 314, A.H. Stroud, Approximate Calculation of Multiple Integrals 188 189 m = 7; 190 G[m][0][0] = 0.333333333333333333333333333333333; 191 G[m][0][1] = 0.333333333333333333333333333333333; 192 G[m][1][0] = 0.79742698535308732239802527616975; 193 G[m][1][1] = 0.1012865073234563388009873619151; 194 G[m][2][0] = 0.10128650732345633880098736191512; 195 G[m][2][1] = 0.79742698535308732239802527616975; 196 G[m][3][0] = 0.10128650732345633880098736191512; 197 G[m][3][1] = 0.10128650732345633880098736191512; 198 G[m][4][0] = 0.05971587178976982045911758097311; 199 G[m][4][1] = 0.47014206410511508977044120951345; 200 G[m][5][0] = 0.47014206410511508977044120951345; 201 G[m][5][1] = 0.05971587178976982045911758097311; 202 G[m][6][0] = 0.47014206410511508977044120951345; 203 G[m][6][1] = 0.47014206410511508977044120951345; 204 205 W[m][0] = 0.5 * 0.225; 206 W[m][1] = 0.5 * 0.12593918054482715259568394550018; 207 W[m][2] = 0.5 * 0.12593918054482715259568394550018; 208 W[m][3] = 0.5 * 0.12593918054482715259568394550018; 209 W[m][4] = 0.5 * 0.13239415278850618073764938783315; 210 W[m][5] = 0.5 * 0.13239415278850618073764938783315; 211 W[m][6] = 0.5 * 0.13239415278850618073764938783315; 212 O[m] = 5; 213 214 // polynomial degree 6 215 /* 12 inner Gauss points, positive weights */ 216 217 // Rule t2-6-12a of the "Encyclopaedia of Cubature Formulas" at 218 // http://www.cs.kuleuven.ac.be/~nines/research/ecf/ecf.html 219 // maintained by Ronald Cools. 220 221 // For further reference: Appendix II, D.A. Dunavant, High degree efficient symmetrical Gaussian quadrature rules for the triangle 222 223 m=12; 224 G[m][0][0] = 0.063089014491502228340331602870819; 225 G[m][0][1] = 0.063089014491502228340331602870819; 226 G[m][1][0] = 0.063089014491502228340331602870819; 227 G[m][1][1] = 0.87382197101699554331933679425836; 228 G[m][2][0] = 0.87382197101699554331933679425836; 229 G[m][2][1] = 0.063089014491502228340331602870819; 230 G[m][3][0] = 0.24928674517091042129163855310702; 231 G[m][3][1] = 0.24928674517091042129163855310702; 232 G[m][4][0] = 0.24928674517091042129163855310702; 233 G[m][4][1] = 0.50142650965817915741672289378596; 234 G[m][5][0] = 0.50142650965817915741672289378596; 235 G[m][5][1] = 0.24928674517091042129163855310702; 236 G[m][6][0] = 0.053145049844816947353249671631398; 237 G[m][6][1] = 0.31035245103378440541660773395655; 238 G[m][7][0] = 0.053145049844816947353249671631398; 239 G[m][7][1] = 0.63650249912139864723014259441205; 240 G[m][8][0] = 0.31035245103378440541660773395655; 241 G[m][8][1] = 0.053145049844816947353249671631398; 242 G[m][9][0] = 0.31035245103378440541660773395655; 243 G[m][9][1] = 0.63650249912139864723014259441205; 244 G[m][10][0] = 0.63650249912139864723014259441205; 245 G[m][10][1] = 0.053145049844816947353249671631398; 246 G[m][11][0] = 0.63650249912139864723014259441205; 247 G[m][11][1] = 0.31035245103378440541660773395655; 248 249 W[m][0] = 0.5 * 0.050844906370206816920936809106869; 250 W[m][1] = 0.5 * 0.050844906370206816920936809106869; 251 W[m][2] = 0.5 * 0.050844906370206816920936809106869; 252 W[m][3] = 0.5 * 0.11678627572637936602528961138558; 253 W[m][4] = 0.5 * 0.11678627572637936602528961138558; 254 W[m][5] = 0.5 * 0.11678627572637936602528961138558; 255 W[m][6] = 0.5 * 0.082851075618373575193553456420442; 256 W[m][7] = 0.5 * 0.082851075618373575193553456420442; 257 W[m][8] = 0.5 * 0.082851075618373575193553456420442; 258 W[m][9] = 0.5 * 0.082851075618373575193553456420442; 259 W[m][10] = 0.5 * 0.082851075618373575193553456420442; 260 W[m][11] = 0.5 * 0.082851075618373575193553456420442; 261 O[m] = 6; 262 263 // polynomial degree 7 264 /* 12 inner Gauss points, positive weights */ 265 // Rule t2-7-12 of the "Encyclopaedia of Cubature Formulas" at 266 // http://www.cs.kuleuven.ac.be/~nines/research/ecf/ecf.html 267 // maintained by Ronald Cools. 268 269 // For further reference: Table 5, K. Gatermann, The construction of symmetric cubature formulas for the square and the triangle 270 271 m=12; 272 G[m][0][0] = 0.0623822650944021181736830009963499; 273 G[m][0][1] = 0.0675178670739160854425571310508685; 274 G[m][1][0] = 0.0675178670739160854425571310508685; 275 G[m][1][1] = 0.870099867831681796383759867952782; 276 G[m][2][0] = 0.870099867831681796383759867952782; 277 G[m][2][1] = 0.0623822650944021181736830009963499; 278 G[m][3][0] = 0.0552254566569266117374791902756449; 279 G[m][3][1] = 0.321502493851981822666307849199202; 280 G[m][4][0] = 0.321502493851981822666307849199202; 281 G[m][4][1] = 0.623272049491091565596212960525153; 282 G[m][5][0] = 0.623272049491091565596212960525153; 283 G[m][5][1] = 0.0552254566569266117374791902756449; 284 G[m][6][0] = 0.0343243029450971464696306424839376; 285 G[m][6][1] = 0.660949196186735657611980310197799; 286 G[m][7][0] = 0.660949196186735657611980310197799; 287 G[m][7][1] = 0.304726500868167195918389047318263; 288 G[m][8][0] = 0.304726500868167195918389047318263; 289 G[m][8][1] = 0.0343243029450971464696306424839376; 290 G[m][9][0] = 0.515842334353591779257463386826430; 291 G[m][9][1] = 0.277716166976391782569581871393723; 292 G[m][10][0] = 0.277716166976391782569581871393723; 293 G[m][10][1] = 0.20644149867001643817295474177985; 294 G[m][11][0] = 0.20644149867001643817295474177985; 295 G[m][11][1] = 0.515842334353591779257463386826430; 296 297 W[m][0] = 0.5 * 0.053034056314872502857508360921478; 298 W[m][1] = 0.5 * 0.053034056314872502857508360921478; 299 W[m][2] = 0.5 * 0.053034056314872502857508360921478; 300 W[m][3] = 0.5 * 0.087762817428892110073539806278575; 301 W[m][4] = 0.5 * 0.087762817428892110073539806278575; 302 W[m][5] = 0.5 * 0.087762817428892110073539806278575; 303 W[m][6] = 0.5 * 0.057550085569963171476890993800437; 304 W[m][7] = 0.5 * 0.057550085569963171476890993800437; 305 W[m][8] = 0.5 * 0.057550085569963171476890993800437; 306 W[m][9] = 0.5 * 0.13498637401960554892539417233284; 307 W[m][10] = 0.5 * 0.13498637401960554892539417233284; 308 W[m][11] = 0.5 * 0.13498637401960554892539417233284; 309 O[m] = 7; 310 311 // polynomial degree 8 312 /* 16 inner Gauss points, positive weights */ 313 // Rule t2-8-16a of the "Encyclopaedia of Cubature Formulas" at 314 // http://www.cs.kuleuven.ac.be/~nines/research/ecf/ecf.html 315 // maintained by Ronald Cools. 316 317 // For further reference: Appendix II, A. Dunavant, High degree efficient symmetrical Gaussian quadrature rules for the triangle 318 319 m=16; 320 G[m][0][0] = 0.33333333333333333333333333333333; 321 G[m][0][1] = 0.33333333333333333333333333333333; 322 G[m][1][0] = 0.17056930775176020662229350149146; 323 G[m][1][1] = 0.17056930775176020662229350149146; 324 G[m][2][0] = 0.17056930775176020662229350149146; 325 G[m][2][1] = 0.65886138449647958675541299701707; 326 G[m][3][0] = 0.65886138449647958675541299701707; 327 G[m][3][1] = 0.17056930775176020662229350149146; 328 G[m][4][0] = 0.050547228317030975458423550596599; 329 G[m][4][1] = 0.050547228317030975458423550596599; 330 G[m][5][0] = 0.050547228317030975458423550596599; 331 G[m][5][1] = 0.89890554336593804908315289880680; 332 G[m][6][0] = 0.89890554336593804908315289880680; 333 G[m][6][1] = 0.050547228317030975458423550596599; 334 G[m][7][0] = 0.45929258829272315602881551449417; 335 G[m][7][1] = 0.45929258829272315602881551449417; 336 G[m][8][0] = 0.45929258829272315602881551449417; 337 G[m][8][1] = 0.08141482341455368794236897101166; 338 G[m][9][0] = 0.08141482341455368794236897101166; 339 G[m][9][1] = 0.45929258829272315602881551449417; 340 G[m][10][0] = 0.72849239295540428124100037917606; 341 G[m][10][1] = 0.26311282963463811342178578628464; 342 G[m][11][0] = 0.72849239295540428124100037917606; 343 G[m][11][1] = 0.00839477740995760533721383453930; 344 G[m][12][0] = 0.26311282963463811342178578628464; 345 G[m][12][1] = 0.72849239295540428124100037917606; 346 G[m][13][0] = 0.26311282963463811342178578628464; 347 G[m][13][1] = 0.00839477740995760533721383453930; 348 G[m][14][0] = 0.00839477740995760533721383453930; 349 G[m][14][1] = 0.72849239295540428124100037917606; 350 G[m][15][0] = 0.00839477740995760533721383453930; 351 G[m][15][1] = 0.26311282963463811342178578628464; 352 353 W[m][0] = 0.5 * 0.14431560767778716825109111048906; 354 W[m][1] = 0.5 * 0.10321737053471825028179155029213; 355 W[m][2] = 0.5 * 0.10321737053471825028179155029213; 356 W[m][3] = 0.5 * 0.10321737053471825028179155029213; 357 W[m][4] = 0.5 * 0.032458497623198080310925928341780; 358 W[m][5] = 0.5 * 0.032458497623198080310925928341780; 359 W[m][6] = 0.5 * 0.032458497623198080310925928341780; 360 W[m][7] = 0.5 * 0.095091634267284624793896104388584; 361 W[m][8] = 0.5 * 0.095091634267284624793896104388584; 362 W[m][9] = 0.5 * 0.095091634267284624793896104388584; 363 W[m][10] = 0.5 * 0.027230314174434994264844690073909; 364 W[m][11] = 0.5 * 0.027230314174434994264844690073909; 365 W[m][12] = 0.5 * 0.027230314174434994264844690073909; 366 W[m][13] = 0.5 * 0.027230314174434994264844690073909; 367 W[m][14] = 0.5 * 0.027230314174434994264844690073909; 368 W[m][15] = 0.5 * 0.027230314174434994264844690073909; 369 O[m] = 8; 370 371 // polynomial degree 9 372 /* 19 inner Gauss points, positive weights */ 373 374 // Rule t2-9-19 of the "Encyclopaedia of Cubature Formulas" at 375 // http://www.cs.kuleuven.ac.be/~nines/research/ecf/ecf.html 376 // maintained by Ronald Cools. 377 378 // For further reference: Appendix II, A. Dunavant, High degree efficient symmetrical Gaussian quadrature rules for the triangle 379 380 m=19; 381 G[m][0][0] = 0.333333333333333333333333333333333; 382 G[m][0][1] = 0.333333333333333333333333333333333; 383 G[m][1][0] = 0.48968251919873762778370692483619; 384 G[m][1][1] = 0.48968251919873762778370692483619; 385 G[m][2][0] = 0.48968251919873762778370692483619; 386 G[m][2][1] = 0.02063496160252474443258615032762; 387 G[m][3][0] = 0.02063496160252474443258615032762; 388 G[m][3][1] = 0.48968251919873762778370692483619; 389 G[m][4][0] = 0.43708959149293663726993036443535; 390 G[m][4][1] = 0.43708959149293663726993036443535; 391 G[m][5][0] = 0.43708959149293663726993036443535; 392 G[m][5][1] = 0.12582081701412672546013927112929; 393 G[m][6][0] = 0.12582081701412672546013927112929; 394 G[m][6][1] = 0.43708959149293663726993036443535; 395 G[m][7][0] = 0.18820353561903273024096128046733; 396 G[m][7][1] = 0.18820353561903273024096128046733; 397 G[m][8][0] = 0.18820353561903273024096128046733; 398 G[m][8][1] = 0.62359292876193453951807743906533; 399 G[m][9][0] = 0.62359292876193453951807743906533; 400 G[m][9][1] = 0.18820353561903273024096128046733; 401 G[m][10][0] = 0.044729513394452709865106589966276; 402 G[m][10][1] = 0.044729513394452709865106589966276; 403 G[m][11][0] = 0.044729513394452709865106589966276; 404 G[m][11][1] = 0.91054097321109458026978682006745; 405 G[m][12][0] = 0.91054097321109458026978682006745; 406 G[m][12][1] = 0.044729513394452709865106589966276; 407 G[m][13][0] = 0.74119859878449802069007987352342; 408 G[m][13][1] = 0.036838412054736283634817598783385; 409 G[m][14][0] = 0.74119859878449802069007987352342; 410 G[m][14][1] = 0.22196298916076569567510252769319; 411 G[m][15][0] = 0.036838412054736283634817598783385; 412 G[m][15][1] = 0.74119859878449802069007987352342; 413 G[m][16][0] = 0.036838412054736283634817598783385; 414 G[m][16][1] = 0.22196298916076569567510252769319; 415 G[m][17][0] = 0.22196298916076569567510252769319; 416 G[m][17][1] = 0.74119859878449802069007987352342; 417 G[m][18][0] = 0.22196298916076569567510252769319; 418 G[m][18][1] = 0.036838412054736283634817598783385; 419 420 W[m][0] = 0.5 * 0.097135796282798833819241982507289; 421 W[m][1] = 0.5 * 0.031334700227139070536854831287209; 422 W[m][2] = 0.5 * 0.031334700227139070536854831287209; 423 W[m][3] = 0.5 * 0.031334700227139070536854831287209; 424 W[m][4] = 0.5 * 0.077827541004774279316739356299404; 425 W[m][5] = 0.5 * 0.077827541004774279316739356299404; 426 W[m][6] = 0.5 * 0.077827541004774279316739356299404; 427 W[m][7] = 0.5 * 0.079647738927210253032891774264045; 428 W[m][8] = 0.5 * 0.079647738927210253032891774264045; 429 W[m][9] = 0.5 * 0.079647738927210253032891774264045; 430 W[m][10] = 0.5 * 0.025577675658698031261678798559000; 431 W[m][11] = 0.5 * 0.025577675658698031261678798559000; 432 W[m][12] = 0.5 * 0.025577675658698031261678798559000; 433 W[m][13] = 0.5 * 0.043283539377289377289377289377289; 434 W[m][14] = 0.5 * 0.043283539377289377289377289377289; 435 W[m][15] = 0.5 * 0.043283539377289377289377289377289; 436 W[m][16] = 0.5 * 0.043283539377289377289377289377289; 437 W[m][17] = 0.5 * 0.043283539377289377289377289377289; 438 W[m][18] = 0.5 * 0.043283539377289377289377289377289; 439 O[m] = 9; 440 441 // polynomial degree 10 442 /* 25 inner Gauss points, positive weights */ 443 444 // Rule t2-10-25a of the "Encyclopaedia of Cubature Formulas" at 445 // http://www.cs.kuleuven.ac.be/~nines/research/ecf/ecf.html 446 // maintained by Ronald Cools. 447 448 // For further reference: M.E. Laursen and M. Gellert, Some criteria for numerically integrated matrices and quadrature formulas for triangles 449 450 m= 25; 451 G[m][0][0] = 0.333333333333333333333333333333333; 452 G[m][0][1] = 0.333333333333333333333333333333333; 453 G[m][1][0] = 0.42508621060209057296952951163804; 454 G[m][1][1] = 0.42508621060209057296952951163804; 455 G[m][2][0] = 0.42508621060209057296952951163804; 456 G[m][2][1] = 0.14982757879581885406094097672391; 457 G[m][3][0] = 0.14982757879581885406094097672391; 458 G[m][3][1] = 0.42508621060209057296952951163804; 459 G[m][4][0] = 0.023308867510000190714466386895980; 460 G[m][4][1] = 0.023308867510000190714466386895980; 461 G[m][5][0] = 0.023308867510000190714466386895980; 462 G[m][5][1] = 0.95338226497999961857106722620804; 463 G[m][6][0] = 0.95338226497999961857106722620804; 464 G[m][6][1] = 0.023308867510000190714466386895980; 465 G[m][7][0] = 0.62830740021349255642083766607883; 466 G[m][7][1] = 0.22376697357697300622568649026820; 467 G[m][8][0] = 0.62830740021349255642083766607883; 468 G[m][8][1] = 0.14792562620953443735347584365296; 469 G[m][9][0] = 0.22376697357697300622568649026820; 470 G[m][9][1] = 0.62830740021349255642083766607883; 471 G[m][10][0] = 0.22376697357697300622568649026820; 472 G[m][10][1] = 0.14792562620953443735347584365296; 473 G[m][11][0] = 0.14792562620953443735347584365296; 474 G[m][11][1] = 0.62830740021349255642083766607883; 475 G[m][12][0] = 0.14792562620953443735347584365296; 476 G[m][12][1] = 0.22376697357697300622568649026820; 477 G[m][13][0] = 0.61131382618139764891875500225390; 478 G[m][13][1] = 0.35874014186443146457815530072385; 479 G[m][14][0] = 0.61131382618139764891875500225390; 480 G[m][14][1] = 0.02994603195417088650308969702225; 481 G[m][15][0] = 0.35874014186443146457815530072385; 482 G[m][15][1] = 0.61131382618139764891875500225390; 483 G[m][16][0] = 0.35874014186443146457815530072385; 484 G[m][16][1] = 0.02994603195417088650308969702225; 485 G[m][17][0] = 0.02994603195417088650308969702225; 486 G[m][17][1] = 0.61131382618139764891875500225390; 487 G[m][18][0] = 0.02994603195417088650308969702225; 488 G[m][18][1] = 0.35874014186443146457815530072385; 489 G[m][19][0] = 0.82107206998562937337354441347218; 490 G[m][19][1] = 0.14329537042686714530585663061732; 491 G[m][20][0] = 0.82107206998562937337354441347218; 492 G[m][20][1] = 0.03563255958750348132059895591050; 493 G[m][21][0] = 0.14329537042686714530585663061732; 494 G[m][21][1] = 0.82107206998562937337354441347218; 495 G[m][22][0] = 0.14329537042686714530585663061732; 496 G[m][22][1] = 0.03563255958750348132059895591050; 497 G[m][23][0] = 0.03563255958750348132059895591050; 498 G[m][23][1] = 0.82107206998562937337354441347218; 499 G[m][24][0] = 0.03563255958750348132059895591050; 500 G[m][24][1] = 0.14329537042686714530585663061732; 501 502 W[m][0] = 0.5 * 0.079894504741239707831247045213386; 503 W[m][1] = 0.5 * 0.071123802232377334639291287398658; 504 W[m][2] = 0.5 * 0.071123802232377334639291287398658; 505 W[m][3] = 0.5 * 0.071123802232377334639291287398658; 506 W[m][4] = 0.5 * 0.0082238186904641955186466203624719; 507 W[m][5] = 0.5 * 0.0082238186904641955186466203624719; 508 W[m][6] = 0.5 * 0.0082238186904641955186466203624719; 509 W[m][7] = 0.5 * 0.045430592296170018007073629243933; 510 W[m][8] = 0.5 * 0.045430592296170018007073629243933; 511 W[m][9] = 0.5 * 0.045430592296170018007073629243933; 512 W[m][10] = 0.5 * 0.045430592296170018007073629243933; 513 W[m][11] = 0.5 * 0.045430592296170018007073629243933; 514 W[m][12] = 0.5 * 0.045430592296170018007073629243933; 515 W[m][13] = 0.5 * 0.037359856234305276826236499001975; 516 W[m][14] = 0.5 * 0.037359856234305276826236499001975; 517 W[m][15] = 0.5 * 0.037359856234305276826236499001975; 518 W[m][16] = 0.5 * 0.037359856234305276826236499001975; 519 W[m][17] = 0.5 * 0.037359856234305276826236499001975; 520 W[m][18] = 0.5 * 0.037359856234305276826236499001975; 521 W[m][19] = 0.5 * 0.030886656884563988782513077004629; 522 W[m][20] = 0.5 * 0.030886656884563988782513077004629; 523 W[m][21] = 0.5 * 0.030886656884563988782513077004629; 524 W[m][22] = 0.5 * 0.030886656884563988782513077004629; 525 W[m][23] = 0.5 * 0.030886656884563988782513077004629; 526 W[m][24] = 0.5 * 0.030886656884563988782513077004629; 527 O[m] = 10; 528 529 // polynomial degree 11 530 /* 28 inner Gauss points, positive weights */ 531 532 // Rule t2-11-28 of the "Encyclopaedia of Cubature Formulas" at 533 // http://www.cs.kuleuven.ac.be/~nines/research/ecf/ecf.html 534 // maintained by Ronald Cools. 535 536 // For further reference: J.N. Lyness and D. Jespersen, Moderate degree symmetric quadrature rules for the triangle 537 538 m=28; 539 G[m][0][0] = 0.858870281282636704039173938058347; 540 G[m][0][1] = 0.141129718717363295960826061941652; 541 G[m][1][0] = 0.858870281282636704039173938058347; 542 G[m][1][1] = 0.0; 543 G[m][2][0] = 0.141129718717363295960826061941652; 544 G[m][2][1] = 0.858870281282636704039173938058347; 545 G[m][3][0] = 0.141129718717363295960826061941652; 546 G[m][3][1] = 0.0; 547 G[m][4][0] = 0.0; 548 G[m][4][1] = 0.858870281282636704039173938058347; 549 G[m][5][0] = 0.0; 550 G[m][5][1] = 0.141129718717363295960826061941652; 551 G[m][6][0] = 0.333333333333333333333333333333333; 552 G[m][6][1] = 0.333333333333333333333333333333333; 553 G[m][7][0] = 0.025989140928287395260032485498841; 554 G[m][7][1] = 0.025989140928287395260032485498841; 555 G[m][8][0] = 0.025989140928287395260032485498841; 556 G[m][8][1] = 0.94802171814342520947993502900232; 557 G[m][9][0] = 0.94802171814342520947993502900232; 558 G[m][9][1] = 0.025989140928287395260032485498841; 559 G[m][10][0] = 0.094287502647922495630569776275405; 560 G[m][10][1] = 0.094287502647922495630569776275405; 561 G[m][11][0] = 0.094287502647922495630569776275405; 562 G[m][11][1] = 0.81142499470415500873886044744919; 563 G[m][12][0] = 0.81142499470415500873886044744919; 564 G[m][12][1] = 0.094287502647922495630569776275405; 565 G[m][13][0] = 0.49463677501721381374163260230644; 566 G[m][13][1] = 0.49463677501721381374163260230644; 567 G[m][14][0] = 0.49463677501721381374163260230644; 568 G[m][14][1] = 0.01072644996557237251673479538713; 569 G[m][15][0] = 0.01072644996557237251673479538713; 570 G[m][15][1] = 0.49463677501721381374163260230644; 571 G[m][16][0] = 0.20734338261451133345293402411297; 572 G[m][16][1] = 0.20734338261451133345293402411297; 573 G[m][17][0] = 0.20734338261451133345293402411297; 574 G[m][17][1] = 0.58531323477097733309413195177407; 575 G[m][18][0] = 0.58531323477097733309413195177407; 576 G[m][18][1] = 0.20734338261451133345293402411297; 577 G[m][19][0] = 0.43890780570049209506106538163613; 578 G[m][19][1] = 0.43890780570049209506106538163613; 579 G[m][20][0] = 0.43890780570049209506106538163613; 580 G[m][20][1] = 0.12218438859901580987786923672775; 581 G[m][21][0] = 0.12218438859901580987786923672775; 582 G[m][21][1] = 0.43890780570049209506106538163613; 583 G[m][22][0] = 0.67793765488259040154212614118875; 584 G[m][22][1] = 0.044841677589130443309052391468801; 585 G[m][23][0] = 0.67793765488259040154212614118875; 586 G[m][23][1] = 0.27722066752827915514882146734245; 587 G[m][24][0] = 0.044841677589130443309052391468801; 588 G[m][24][1] = 0.67793765488259040154212614118875; 589 G[m][25][0] = 0.044841677589130443309052391468801; 590 G[m][25][1] = 0.27722066752827915514882146734245; 591 G[m][26][0] = 0.27722066752827915514882146734245; 592 G[m][26][1] = 0.67793765488259040154212614118875; 593 G[m][27][0] = 0.27722066752827915514882146734245; 594 G[m][27][1] = 0.044841677589130443309052391468801; 595 596 W[m][0] = 0.5 * 0.0073623837833005542642588950473806; 597 W[m][1] = 0.5 * 0.0073623837833005542642588950473806; 598 W[m][2] = 0.5 * 0.0073623837833005542642588950473806; 599 600 W[m][3] = 0.5 * 0.0073623837833005542642588950473806; 601 W[m][4] = 0.5 * 0.0073623837833005542642588950473806; 602 W[m][5] = 0.5 * 0.0073623837833005542642588950473806; 603 604 W[m][6] = 0.5 * 0.087977301162232238798093169321456; 605 W[m][7] = 0.5 * 0.0087443115537360230495164287998252; 606 W[m][8] = 0.5 * 0.0087443115537360230495164287998252; 607 W[m][9] = 0.5 * 0.0087443115537360230495164287998252; 608 609 W[m][10] = 0.5 * 0.038081571993934937515024339435614; 610 W[m][11] = 0.5 * 0.038081571993934937515024339435614; 611 W[m][12] = 0.5 * 0.038081571993934937515024339435614; 612 613 W[m][13] = 0.5 * 0.018855448056131292058476782591115; 614 W[m][14] = 0.5 * 0.018855448056131292058476782591115; 615 W[m][15] = 0.5 * 0.018855448056131292058476782591115; 616 617 W[m][16] = 0.5 * 0.072159697544739526124029988586463; 618 W[m][17] = 0.5 * 0.072159697544739526124029988586463; 619 W[m][18] = 0.5 * 0.072159697544739526124029988586463; 620 621 W[m][19] = 0.5 * 0.069329138705535899841765650903814; 622 W[m][20] = 0.5 * 0.069329138705535899841765650903814; 623 W[m][21] = 0.5 * 0.069329138705535899841765650903814; 624 625 W[m][22] = 0.5 * 0.041056315429288566641652314907294; 626 W[m][23] = 0.5 * 0.041056315429288566641652314907294; 627 W[m][24] = 0.5 * 0.041056315429288566641652314907294; 628 629 W[m][25] = 0.5 * 0.041056315429288566641652314907294; 630 W[m][26] = 0.5 * 0.041056315429288566641652314907294; 631 W[m][27] = 0.5 * 0.041056315429288566641652314907294; 632 O[m] = 11; 633 634 // polynomial degree 12 635 /* 33 inner Gauss points, positive weights */ 636 637 // Rule t2-12-33 of the "Encyclopaedia of Cubature Formulas" at 638 // http://www.cs.kuleuven.ac.be/~nines/research/ecf/ecf.html 639 // maintained by Ronald Cools. 640 641 // For further reference: Appendix II, A. Dunavant, High degree efficient symmetrical Gaussian quadrature rules for the triangle 642 643 m=33; 644 G[m][0][0] = 0.02356522045239; 645 G[m][0][1] = 0.488217389773805; 646 G[m][1][0] = 0.488217389773805; 647 G[m][1][1] = 0.02356522045239; 648 G[m][2][0] = 0.488217389773805; 649 G[m][2][1] = 0.488217389773805; 650 G[m][3][0] = 0.43972439229446; 651 G[m][3][1] = 0.43972439229446; 652 G[m][4][0] = 0.43972439229446; 653 G[m][4][1] = 0.120551215411079; 654 G[m][5][0] = 0.120551215411079; 655 G[m][5][1] = 0.43972439229446; 656 G[m][6][0] = 0.271210385012116; 657 G[m][6][1] = 0.271210385012116; 658 G[m][7][0] = 0.271210385012116; 659 G[m][7][1] = 0.457579229975768; 660 G[m][8][0] = 0.457579229975768; 661 G[m][8][1] = 0.271210385012116; 662 G[m][9][0] = 0.127576145541586; 663 G[m][9][1] = 0.127576145541586; 664 G[m][10][0] = 0.127576145541586; 665 G[m][10][1] = 0.7448477089168279; 666 G[m][11][0] = 0.7448477089168279; 667 G[m][11][1] = 0.127576145541586; 668 G[m][12][0] = 0.02131735045321; 669 G[m][12][1] = 0.02131735045321; 670 G[m][13][0] = 0.02131735045321; 671 G[m][13][1] = 0.9573652990935799; 672 G[m][14][0] = 0.9573652990935799; 673 G[m][14][1] = 0.02131735045321; 674 G[m][15][0] = 0.115343494534698; 675 G[m][15][1] = 0.275713269685514; 676 G[m][16][0] = 0.115343494534698; 677 G[m][16][1] = 0.6089432357797879; 678 G[m][17][0] = 0.275713269685514; 679 G[m][17][1] = 0.115343494534698; 680 G[m][18][0] = 0.275713269685514; 681 G[m][18][1] = 0.6089432357797879; 682 G[m][19][0] = 0.6089432357797879; 683 G[m][19][1] = 0.115343494534698; 684 G[m][20][0] = 0.6089432357797879; 685 G[m][20][1] = 0.275713269685514; 686 G[m][21][0] = 0.022838332222257; 687 G[m][21][1] = 0.28132558098994; 688 G[m][22][0] = 0.022838332222257; 689 G[m][22][1] = 0.6958360867878031; 690 G[m][23][0] = 0.28132558098994; 691 G[m][23][1] = 0.022838332222257; 692 G[m][24][0] = 0.28132558098994; 693 G[m][24][1] = 0.6958360867878031; 694 G[m][25][0] = 0.6958360867878031; 695 G[m][25][1] = 0.022838332222257; 696 G[m][26][0] = 0.6958360867878031; 697 G[m][26][1] = 0.28132558098994; 698 G[m][27][0] = 0.02573405054833; 699 G[m][27][1] = 0.116251915907597; 700 G[m][28][0] = 0.02573405054833; 701 G[m][28][1] = 0.858014033544073; 702 G[m][29][0] = 0.116251915907597; 703 G[m][29][1] = 0.02573405054833; 704 G[m][30][0] = 0.116251915907597; 705 G[m][30][1] = 0.858014033544073; 706 G[m][31][0]= 0.858014033544073; 707 G[m][31][1] =0.02573405054833; 708 G[m][32][0] =0.858014033544073; 709 G[m][32][1]= 0.116251915907597; 710 711 W[m][0] = 0.5 * 0.025731066440455; 712 W[m][1] = 0.5 * 0.025731066440455; 713 W[m][2] = 0.5 * 0.025731066440455; 714 W[m][3] = 0.5 * 0.043692544538038; 715 W[m][4] = 0.5 * 0.043692544538038; 716 W[m][5] = 0.5 * 0.043692544538038; 717 W[m][6] = 0.5 * 0.062858224217885; 718 W[m][7] = 0.5 * 0.062858224217885; 719 W[m][8] = 0.5 * 0.062858224217885; 720 W[m][9] = 0.5 * 0.034796112930709; 721 W[m][10] = 0.5 * 0.034796112930709; 722 W[m][11] = 0.5 * 0.034796112930709; 723 W[m][12] = 0.5 * 0.006166261051559; 724 W[m][13] = 0.5 * 0.006166261051559; 725 W[m][14] = 0.5 * 0.006166261051559; 726 W[m][15] = 0.5 * 0.040371557766381; 727 W[m][16] = 0.5 * 0.040371557766381; 728 W[m][17] = 0.5 * 0.040371557766381; 729 W[m][18] = 0.5 * 0.040371557766381; 730 W[m][19] = 0.5 * 0.040371557766381; 731 W[m][20] = 0.5 * 0.040371557766381; 732 W[m][21] = 0.5 * 0.022356773202303; 733 W[m][22] = 0.5 * 0.022356773202303; 734 W[m][23] = 0.5 * 0.022356773202303; 735 W[m][24] = 0.5 * 0.022356773202303; 736 W[m][25] = 0.5 * 0.022356773202303; 737 W[m][26] = 0.5 * 0.022356773202303; 738 W[m][27] = 0.5 * 0.017316231108659; 739 W[m][28] = 0.5 * 0.017316231108659; 740 W[m][29] = 0.5 * 0.017316231108659; 741 W[m][30] = 0.5 * 0.017316231108659; 742 W[m][31] = 0.5 * 0.017316231108659; 743 W[m][32] = 0.5 * 0.017316231108659; 744 O[m] = 12; 745 } 746 point(int m,int i)747 FieldVector<double, 2> point(int m, int i) 748 { 749 return G[m][i]; 750 } 751 weight(int m,int i)752 double weight (int m, int i) 753 { 754 return W[m][i]; 755 } 756 order(int m)757 int order (int m) 758 { 759 return O[m]; 760 } 761 762 private: 763 FieldVector<double, 2> G[MAXP+1][MAXP]; 764 765 double W[MAXP+1][MAXP]; // weights associated with points 766 int O[MAXP+1]; // order of the rule 767 }; 768 769 template<typename ct> SimplexQuadratureRule(int p)770 SimplexQuadratureRule<ct,2>::SimplexQuadratureRule(int p) : QuadratureRule<ct,2>(GeometryTypes::triangle) 771 { 772 int m; 773 if (p>SimplexQuadraturePoints<2>::highest_order) 774 DUNE_THROW(QuadratureOrderOutOfRange, 775 "QuadratureRule for order " << p << " and GeometryType " 776 << this->type() << " not available"); 777 778 switch(p) 779 { 780 case 0 : // to be verified 781 m=1; // to be verified 782 break; 783 case 1 : 784 m=1; 785 break; 786 case 2 : 787 m=3; 788 break; 789 case 3 : 790 m=4; 791 break; 792 case 4 : 793 m=6; 794 break; 795 case 5 : 796 m=7; 797 break; 798 case 6 : 799 m=12; 800 break; 801 case 7 : 802 m=12; 803 break; 804 case 8 : 805 m=16; 806 break; 807 case 9 : 808 m=19; 809 break; 810 case 10 : 811 m=25; 812 break; 813 case 11 : 814 m=28; 815 break; 816 case 12 : 817 m=33; 818 break; 819 default : m=33; 820 } 821 822 this->delivered_order = SimplexQuadraturePointsSingleton<2>::sqp.order(m); 823 824 for(int i=0; i<m; ++i) 825 { 826 FieldVector<ct,2> local = SimplexQuadraturePointsSingleton<2>::sqp.point(m,i); 827 ct weight = SimplexQuadraturePointsSingleton<2>::sqp.weight(m,i); 828 // put in container 829 this->push_back(QuadraturePoint<ct,2>(local,weight)); 830 } 831 } 832 833 //! 834 template<> 835 class SimplexQuadraturePoints<3> 836 { 837 public: 838 enum { MAXP=15}; 839 enum { highest_order=5 }; 840 841 //! initialize quadrature points on the interval for all orders SimplexQuadraturePoints()842 SimplexQuadraturePoints() 843 { 844 init(); 845 } 846 init()847 void init() 848 { 849 int m = 0; 850 O[m] = 0; 851 852 // polynom degree 1 853 854 // Rule t3-1-1 of the "Encyclopaedia of Cubature Formulas" at 855 // http://www.cs.kuleuven.ac.be/~nines/research/ecf/ecf.html 856 // maintained by Ronald Cools. 857 858 // For further reference: Rule 1-1, P. 307, A.H. Stroud, Approximate Calculation of Multiple Integrals 859 860 m = 1; 861 G[m][0][0] = 0.25; 862 G[m][0][1] = 0.25; 863 G[m][0][2] = 0.25; 864 865 W[m][0] = 1.0/6.0; 866 O[m] = 1; 867 868 // polynom degree 2 869 // symmetric 870 871 // Rule t3-2-4a of the "Encyclopaedia of Cubature Formulas" at 872 // http://www.cs.kuleuven.ac.be/~nines/research/ecf/ecf.html 873 // maintained by Ronald Cools. 874 875 // For further reference: Rule 2-1, P. 307, A.H. Stroud, Approximate Calculation of Multiple Integrals 876 877 m = 4; 878 static const double m_4_a = 0.585410196624968500; 879 static const double m_4_b = 0.138196601125010500; 880 G[m][0] = m_4_b; 881 G[m][1] = m_4_b; 882 G[m][2] = m_4_b; 883 G[m][3] = m_4_b; 884 G[m][0][0] = m_4_a; 885 G[m][1][1] = m_4_a; 886 G[m][2][2] = m_4_a; 887 888 W[m][0] = 1.0/4.0/6.0; 889 W[m][1] = 1.0/4.0/6.0; 890 W[m][2] = 1.0/4.0/6.0; 891 W[m][3] = 1.0/4.0/6.0; 892 O[m] = 2; 893 894 // polynom degree 3 895 // symmetric 896 897 // Rule t3-3-8b of the "Encyclopaedia of Cubature Formulas" at 898 // http://www.cs.kuleuven.ac.be/~nines/research/ecf/ecf.html 899 // maintained by Ronald Cools. 900 901 // For further reference: Rule 3-7, P. 309, A.H. Stroud, Approximate Calculation of Multiple Integrals 902 903 m = 8; 904 G[m][0][0] = 0.0; 905 G[m][0][1] = 0.0; 906 G[m][0][2] = 0.0; 907 G[m][1][0] = 1.0; 908 G[m][1][1] = 0.0; 909 G[m][1][2] = 0.0; 910 G[m][2][0] = 0.0; 911 G[m][2][1] = 1.0; 912 G[m][2][2] = 0.0; 913 G[m][3][0] = 0.0; 914 G[m][3][1] = 0.0; 915 G[m][3][2] = 1.0; 916 G[m][4][0] = 1.0/3.0; 917 G[m][4][1] = 1.0/3.0; 918 G[m][4][2] = 0.0; 919 G[m][5][0] = 1.0/3.0; 920 G[m][5][1] = 0.0; 921 G[m][5][2] = 1.0/3.0; 922 G[m][6][0] = 0.0; 923 G[m][6][1] = 1.0/3.0; 924 G[m][6][2] = 1.0/3.0; 925 G[m][7][0] = 1.0/3.0; 926 G[m][7][1] = 1.0/3.0; 927 G[m][7][2] = 1.0/3.0; 928 W[m][0] = 0.025/6.0; 929 W[m][1] = 0.025/6.0; 930 W[m][2] = 0.025/6.0; 931 W[m][3] = 0.025/6.0; 932 W[m][4] = 0.225/6.0; 933 W[m][5] = 0.225/6.0; 934 W[m][6] = 0.225/6.0; 935 W[m][7] = 0.225/6.0; 936 O[m] = 3; 937 938 939 // polynomial degree 5 940 // symmetric points 941 942 // Rule t3-5-15a (actual values for weights etc. not available here) of the "Encyclopaedia of Cubature Formulas" at 943 // http://www.cs.kuleuven.ac.be/~nines/research/ecf/ecf.html 944 // maintained by Ronald Cools. 945 946 // For further reference: Rule 5-1, P. 315, A.H. Stroud, Approximate Calculation of Multiple Integrals 947 948 static const double s_1=0.09197107805272303279; /* (7 - sqrt(15) ) / 34 */ 949 static const double s_2=0.31979362782962990839; /* (7 + sqrt(15) ) / 34 */ 950 static const double t_1=0.72408676584183090164; /* (13 + 3*sqrt(15) ) / 34 */ 951 static const double t_2=0.04061911651111027484; /* (13 - 3*sqrt(15) ) / 34 */ 952 static const double u =0.05635083268962915574; /* (10 - 2*sqrt(15) ) / 40 */ 953 static const double v =0.44364916731037084426; /* (10 + 2*sqrt(15) ) / 40 */ 954 static const double A =0.019753086419753086420; /* 16 / 135 / vol */ 955 static const double B_1=0.011989513963169770001; /* (2665 + 14*sqrt(15) ) / 37800 / vol */ 956 static const double B_2=0.011511367871045397547; /* (2665 - 14*sqrt(15) ) / 37800 / vol */ 957 static const double C =0.0088183421516754850088; /* 20 / 378 / vol */ 958 959 m=15; 960 G[m][0][0] = 0.25; 961 G[m][0][1] = 0.25; 962 G[m][0][2] = 0.25; 963 G[m][1][0] = s_1; 964 G[m][1][1] = s_1; 965 G[m][1][2] = s_1; 966 G[m][2][0] = t_1; 967 G[m][2][1] = s_1; 968 G[m][2][2] = s_1; 969 G[m][3][0] = s_1; 970 G[m][3][1] = t_1; 971 G[m][3][2] = s_1; 972 G[m][4][0] = s_1; 973 G[m][4][1] = s_1; 974 G[m][4][2] = t_1; 975 G[m][5][0] = s_2; 976 G[m][5][1] = s_2; 977 G[m][5][2] = s_2; 978 G[m][6][0] = t_2; 979 G[m][6][1] = s_2; 980 G[m][6][2] = s_2; 981 G[m][7][0] = s_2; 982 G[m][7][1] = t_2; 983 G[m][7][2] = s_2; 984 G[m][8][0] = s_2; 985 G[m][8][1] = s_2; 986 G[m][8][2] = t_2; 987 G[m][9][0] = v; 988 G[m][9][1] = u; 989 G[m][9][2] = u; 990 G[m][10][0] = u; 991 G[m][10][1] = v; 992 G[m][10][2] = u; 993 G[m][11][0] = u; 994 G[m][11][1] = u; 995 G[m][11][2] = v; 996 G[m][12][0] = v; 997 G[m][12][1] = v; 998 G[m][12][2] = u; 999 G[m][13][0] = v; 1000 G[m][13][1] = u; 1001 G[m][13][2] = v; 1002 G[m][14][0] = u; 1003 G[m][14][1] = v; 1004 G[m][14][2] = v; 1005 1006 W[m][0] = A; 1007 W[m][1] = B_1; 1008 W[m][2] = B_1; 1009 W[m][3] = B_1; 1010 W[m][4] = B_1; 1011 W[m][5] = B_2; 1012 W[m][6] = B_2; 1013 W[m][7] = B_2; 1014 W[m][8] = B_2; 1015 W[m][9] = C; 1016 W[m][10] = C; 1017 W[m][11] = C; 1018 W[m][12] = C; 1019 W[m][13] = C; 1020 W[m][14] = C; 1021 O[m] = 5; 1022 1023 } 1024 point(int m,int i)1025 FieldVector<double, 3> point(int m, int i) 1026 { 1027 return G[m][i]; 1028 } 1029 weight(int m,int i)1030 double weight (int m, int i) 1031 { 1032 return W[m][i]; 1033 } 1034 order(int m)1035 int order (int m) 1036 { 1037 return O[m]; 1038 } 1039 1040 private: 1041 FieldVector<double, 3> G[MAXP+1][MAXP]; 1042 double W[MAXP+1][MAXP]; // weights associated with points 1043 int O[MAXP+1]; // order of the rule 1044 }; 1045 1046 template<typename ct> SimplexQuadratureRule(int p)1047 SimplexQuadratureRule<ct,3>::SimplexQuadratureRule(int p) : QuadratureRule<ct,3>(GeometryTypes::tetrahedron) 1048 { 1049 int m; 1050 if (p>SimplexQuadraturePoints<3>::highest_order) 1051 DUNE_THROW(QuadratureOrderOutOfRange, 1052 "QuadratureRule for order " << p << " and GeometryType " 1053 << this->type() << " not available"); 1054 switch(p) 1055 { 1056 case 0 : 1057 m=1; 1058 break; 1059 case 1 : 1060 m=1; 1061 break; 1062 case 2 : 1063 m=4; 1064 break; 1065 case 3 : 1066 m=8; 1067 break; 1068 case 4 : 1069 case 5 : 1070 m=15; 1071 break; 1072 default : m=15; 1073 } 1074 this->delivered_order = SimplexQuadraturePointsSingleton<3>::sqp.order(m); 1075 1076 for(int i=0; i<m; ++i) 1077 { 1078 FieldVector<ct,3> local = SimplexQuadraturePointsSingleton<3>::sqp.point(m,i); 1079 ct weight = SimplexQuadraturePointsSingleton<3>::sqp.weight(m,i); 1080 // put in container 1081 this->push_back(QuadraturePoint<ct,3>(local,weight)); 1082 } 1083 1084 } 1085 1086 } // end namespace Dune 1087 1088 #endif // DUNE_GEOMETRY_QUADRATURE_SIMPLEX_HH 1089