MFEM  v3.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
intrules.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 // Implementation of IntegrationRule(s) classes
13 
14 // Acknowledgment: Some of the high-precision triangular and tetrahedral
15 // quadrature rules below were obtained from the Encyclopaedia of Cubature
16 // Formulas at http://nines.cs.kuleuven.be/research/ecf/ecf.html
17 
18 #include "fem.hpp"
19 #include <cmath>
20 
21 using namespace std;
22 
23 namespace mfem
24 {
25 
26 IntegrationRule::IntegrationRule(IntegrationRule &irx, IntegrationRule &iry)
27 {
28  int i, j, nx, ny;
29 
30  nx = irx.GetNPoints();
31  ny = iry.GetNPoints();
32  SetSize(nx * ny);
33 
34  for (j = 0; j < ny; j++)
35  {
36  IntegrationPoint &ipy = iry.IntPoint(j);
37  for (i = 0; i < nx; i++)
38  {
39  IntegrationPoint &ipx = irx.IntPoint(i);
40  IntegrationPoint &ip = IntPoint(j*nx+i);
41 
42  ip.x = ipx.x;
43  ip.y = ipy.x;
44  ip.weight = ipx.weight * ipy.weight;
45  }
46  }
47 }
48 
49 IntegrationRule::IntegrationRule(IntegrationRule &irx, IntegrationRule &iry,
50  IntegrationRule &irz)
51 {
52  const int nx = irx.GetNPoints();
53  const int ny = iry.GetNPoints();
54  const int nz = irz.GetNPoints();
55  SetSize(nx*ny*nz);
56 
57  for (int iz = 0; iz < nz; ++iz)
58  {
59  IntegrationPoint &ipz = irz.IntPoint(iz);
60  for (int iy = 0; iy < ny; ++iy)
61  {
62  IntegrationPoint &ipy = iry.IntPoint(iy);
63  for (int ix = 0; ix < nx; ++ix)
64  {
65  IntegrationPoint &ipx = irx.IntPoint(ix);
66  IntegrationPoint &ip = IntPoint(iz*nx*ny + iy*nx + ix);
67 
68  ip.x = ipx.x;
69  ip.y = ipy.x;
70  ip.z = ipz.x;
71  ip.weight = ipx.weight*ipy.weight*ipz.weight;
72  }
73  }
74  }
75 }
76 
77 void IntegrationRule::GaussianRule()
78 {
79  int n = Size();
80  int m = (n+1)/2;
81  int i, j;
82  double p1, p2, p3;
83  double pp, z;
84  for (i = 1; i <= m; i++)
85  {
86  z = cos(M_PI * (i - 0.25) / (n + 0.5));
87 
88  while (1)
89  {
90  p1 = 1;
91  p2 = 0;
92  for (j = 1; j <= n; j++)
93  {
94  p3 = p2;
95  p2 = p1;
96  p1 = ((2 * j - 1) * z * p2 - (j - 1) * p3) / j;
97  }
98  // p1 is Legendre polynomial
99 
100  pp = n * (z*p1-p2) / (z*z - 1);
101 
102  if (fabs(p1/pp) < 2e-16) { break; }
103 
104  z = z - p1/pp;
105  }
106 
107  z = ((1 - z) + p1/pp)/2;
108 
109  IntPoint(i-1).x = z;
110  IntPoint(n-i).x = 1 - z;
111  IntPoint(i-1).weight =
112  IntPoint(n-i).weight = 1./(4*z*(1 - z)*pp*pp);
113  }
114 }
115 
116 void IntegrationRule::UniformRule()
117 {
118  int i;
119  double h;
120 
121  h = 1.0 / (Size() - 1);
122  for (i = 0; i < Size(); i++)
123  {
124  IntPoint(i).x = double(i) / (Size() - 1);
125  IntPoint(i).weight = h;
126  }
127  IntPoint(0).weight = 0.5 * h;
128  IntPoint(Size()-1).weight = 0.5 * h;
129 }
130 
131 void IntegrationRule::GrundmannMollerSimplexRule(int s, int n)
132 {
133  // for pow on older compilers
134  using std::pow;
135  const int d = 2*s + 1;
136  Vector fact(d + n + 1);
137  Array<int> beta(n), sums(n);
138 
139  fact(0) = 1.;
140  for (int i = 1; i < fact.Size(); i++)
141  {
142  fact(i) = fact(i - 1)*i;
143  }
144 
145  // number of points is \binom{n + s + 1}{n + 1}
146  int np = 1, f = 1;
147  for (int i = 0; i <= n; i++)
148  {
149  np *= (s + i + 1), f *= (i + 1);
150  }
151  np /= f;
152  SetSize(np);
153 
154  int pt = 0;
155  for (int i = 0; i <= s; i++)
156  {
157  double weight;
158 
159  weight = pow(2., -2*s)*pow(static_cast<double>(d + n - 2*i),
160  d)/fact(i)/fact(d + n - i);
161  if (i%2)
162  {
163  weight = -weight;
164  }
165 
166  // loop over all beta : beta_0 + ... + beta_{n-1} <= s - i
167  int k = s - i;
168  beta = 0;
169  sums = 0;
170  while (true)
171  {
172  IntegrationPoint &ip = IntPoint(pt++);
173  ip.weight = weight;
174  ip.x = double(2*beta[0] + 1)/(d + n - 2*i);
175  ip.y = double(2*beta[1] + 1)/(d + n - 2*i);
176  if (n == 3)
177  {
178  ip.z = double(2*beta[2] + 1)/(d + n - 2*i);
179  }
180 
181  int j = 0;
182  while (sums[j] == k)
183  {
184  beta[j++] = 0;
185  if (j == n)
186  {
187  goto done_beta;
188  }
189  }
190  beta[j]++;
191  sums[j]++;
192  for (j--; j >= 0; j--)
193  {
194  sums[j] = sums[j+1];
195  }
196  }
197  done_beta:
198  ;
199  }
200 }
201 
202 
203 IntegrationRules IntRules(0);
204 
205 IntegrationRules RefinedIntRules(1);
206 
207 IntegrationRules::IntegrationRules(int Ref)
208 {
209  refined = Ref;
210  if (refined < 0) { own_rules = 0; return; }
211 
212  own_rules = 1;
213 
214  PointIntRules.SetSize(2);
215  PointIntRules = NULL;
216 
217  SegmentIntRules.SetSize(32);
218  SegmentIntRules = NULL;
219 
220  // TriangleIntegrationRule() assumes that this size is >= 26
221  TriangleIntRules.SetSize(32);
222  TriangleIntRules = NULL;
223 
224  SquareIntRules.SetSize(32);
225  SquareIntRules = NULL;
226 
227  // TetrahedronIntegrationRule() assumes that this size is >= 10
228  TetrahedronIntRules.SetSize(32);
229  TetrahedronIntRules = NULL;
230 
231  CubeIntRules.SetSize(32);
232  CubeIntRules = NULL;
233 }
234 
235 const IntegrationRule &IntegrationRules::Get(int GeomType, int Order)
236 {
237  Array<IntegrationRule *> *ir_array;
238 
239  switch (GeomType)
240  {
241  case Geometry::POINT: ir_array = &PointIntRules; Order = 0; break;
242  case Geometry::SEGMENT: ir_array = &SegmentIntRules; break;
243  case Geometry::TRIANGLE: ir_array = &TriangleIntRules; break;
244  case Geometry::SQUARE: ir_array = &SquareIntRules; break;
245  case Geometry::TETRAHEDRON: ir_array = &TetrahedronIntRules; break;
246  case Geometry::CUBE: ir_array = &CubeIntRules; break;
247  default:
248  mfem_error("IntegrationRules::Get(...) : Unknown geometry type!");
249  ir_array = NULL;
250  }
251 
252  if (Order < 0)
253  {
254  Order = 0;
255  }
256 
257  if (!HaveIntRule(*ir_array, Order))
258  {
259  GenerateIntegrationRule(GeomType, Order);
260  }
261 
262  return *(*ir_array)[Order];
263 }
264 
265 void IntegrationRules::Set(int GeomType, int Order, IntegrationRule &IntRule)
266 {
267  Array<IntegrationRule *> *ir_array;
268 
269  switch (GeomType)
270  {
271  case Geometry::POINT: ir_array = &PointIntRules; break;
272  case Geometry::SEGMENT: ir_array = &SegmentIntRules; break;
273  case Geometry::TRIANGLE: ir_array = &TriangleIntRules; break;
274  case Geometry::SQUARE: ir_array = &SquareIntRules; break;
275  case Geometry::TETRAHEDRON: ir_array = &TetrahedronIntRules; break;
276  case Geometry::CUBE: ir_array = &CubeIntRules; break;
277  default:
278  mfem_error("IntegrationRules::Set(...) : Unknown geometry type!");
279  ir_array = NULL;
280  }
281 
282  if (HaveIntRule(*ir_array, Order))
283  {
284  MFEM_ABORT("Overwriting set rules is not supported!");
285  }
286 
287  AllocIntRule(*ir_array, Order);
288 
289  (*ir_array)[Order] = &IntRule;
290 }
291 
292 void IntegrationRules::DeleteIntRuleArray(Array<IntegrationRule *> &ir_array)
293 {
294  int i;
295  IntegrationRule *ir = NULL;
296 
297  // Many of the intrules have multiple contiguous copies in the ir_array
298  // so we have to be careful to not delete them twice.
299  for (i = 0; i < ir_array.Size(); i++)
300  {
301  if (ir_array[i] != NULL && ir_array[i] != ir)
302  {
303  ir = ir_array[i];
304  delete (ir_array[i]);
305  }
306  }
307 }
308 
309 IntegrationRules::~IntegrationRules()
310 {
311  if (!own_rules) { return; }
312 
313  DeleteIntRuleArray(PointIntRules);
314  DeleteIntRuleArray(SegmentIntRules);
315  DeleteIntRuleArray(TriangleIntRules);
316  DeleteIntRuleArray(SquareIntRules);
317  DeleteIntRuleArray(TetrahedronIntRules);
318  DeleteIntRuleArray(CubeIntRules);
319 }
320 
321 
322 IntegrationRule *IntegrationRules::GenerateIntegrationRule(int GeomType,
323  int Order)
324 {
325  switch (GeomType)
326  {
327  case Geometry::POINT:
328  return PointIntegrationRule(Order);
329  case Geometry::SEGMENT:
330  return SegmentIntegrationRule(Order);
331  case Geometry::TRIANGLE:
332  return TriangleIntegrationRule(Order);
333  case Geometry::SQUARE:
334  return SquareIntegrationRule(Order);
335  case Geometry::TETRAHEDRON:
336  return TetrahedronIntegrationRule(Order);
337  case Geometry::CUBE:
338  return CubeIntegrationRule(Order);
339  default:
340  mfem_error("IntegrationRules::Set(...) : Unknown geometry type!");
341  return NULL;
342  }
343 }
344 
345 
346 // Integration rules for a point
347 IntegrationRule *IntegrationRules::PointIntegrationRule(int Order)
348 {
349  if (Order > 1)
350  {
351  mfem_error("Point Integration Rule of Order > 1 not defined");
352  return NULL;
353  }
354 
355  PointIntRules[0] = new IntegrationRule(1);
356  PointIntRules[0] -> IntPoint(0).x = .0;
357  PointIntRules[0] -> IntPoint(0).weight = 1.;
358 
359  PointIntRules[1] = PointIntRules[0];
360 
361  return PointIntRules[0];
362 }
363 
364 // Integration rules for line segment [0,1]
365 IntegrationRule *IntegrationRules::SegmentIntegrationRule(int Order)
366 {
367  int i = (Order / 2) * 2 + 1; // Get closest odd # >= Order
368 
369  AllocIntRule(SegmentIntRules, i);
370 
371  if (refined)
372  {
373  int n = i/2 + 1;
374 
375  IntegrationRule *tmp = new IntegrationRule(n);
376  tmp->GaussianRule();
377 
378  IntegrationRule *ir = new IntegrationRule(2*n);
379  SegmentIntRules[i-1] = SegmentIntRules[i] = ir;
380  for (int j = 0; j < n; j++)
381  {
382  ir->IntPoint(j).x = tmp->IntPoint(j).x/2.0;
383  ir->IntPoint(j).weight = tmp->IntPoint(j).weight/2.0;
384  ir->IntPoint(j+n).x = 0.5 + tmp->IntPoint(j).x/2.0;
385  ir->IntPoint(j+n).weight = tmp->IntPoint(j).weight/2.0;
386  }
387  delete tmp;
388 
389  return ir;
390  }
391 
392  switch (Order / 2)
393  {
394  case 0: // 1 point - 1 degree
395  SegmentIntRules[0] = SegmentIntRules[1] = new IntegrationRule(1);
396  SegmentIntRules[0] -> IntPoint(0).x = .5;
397  SegmentIntRules[0] -> IntPoint(0).weight = 1.;
398  return SegmentIntRules[0];
399  case 1: // 2 point - 3 degree
400  SegmentIntRules[2] = SegmentIntRules[3] = new IntegrationRule(2);
401  SegmentIntRules[2] -> IntPoint(0).x = 0.21132486540518711775;
402  SegmentIntRules[2] -> IntPoint(0).weight = .5;
403  SegmentIntRules[2] -> IntPoint(1).x = 0.78867513459481288225;
404  SegmentIntRules[2] -> IntPoint(1).weight = .5;
405  return SegmentIntRules[2];
406  case 2: // 3 point - 5 degree
407  SegmentIntRules[4] = SegmentIntRules[5] = new IntegrationRule(3);
408  SegmentIntRules[4] -> IntPoint(0).x = 0.11270166537925831148;
409  SegmentIntRules[4] -> IntPoint(0).weight = 5./18.;
410  SegmentIntRules[4] -> IntPoint(1).x = 0.5;
411  SegmentIntRules[4] -> IntPoint(1).weight = 4./9.;
412  SegmentIntRules[4] -> IntPoint(2).x = 0.88729833462074168852;
413  SegmentIntRules[4] -> IntPoint(2).weight = 5./18.;
414  return SegmentIntRules[4];
415  default:
416  SegmentIntRules[i-1] = SegmentIntRules[i] = new IntegrationRule(i/2+1);
417  SegmentIntRules[i]->GaussianRule();
418  return SegmentIntRules[i];
419  }
420 }
421 
422 // Integration rules for reference triangle {[0,0],[1,0],[0,1]}
423 IntegrationRule *IntegrationRules::TriangleIntegrationRule(int Order)
424 {
425  IntegrationRule *ir = NULL;
426 
427  // assuming that orders <= 25 are pre-allocated
428  switch (Order)
429  {
430  case 0: // 1 point - 0 degree
431  case 1:
432  TriangleIntRules[0] =
433  TriangleIntRules[1] = ir = new IntegrationRule(1);
434  TriangleIntRules[0]->AddTriMidPoint(0, 0.5);
435  return ir;
436 
437  case 2: // 3 point - 2 degree
438  TriangleIntRules[2] = ir = new IntegrationRule(3);
439  // interior points
440  TriangleIntRules[2]->AddTriPoints3(0, 1./6., 1./6.);
441  return ir;
442 
443  case 3: // 4 point - 3 degree (has one negative weight)
444  TriangleIntRules[3] = ir = new IntegrationRule(4);
445  ir->AddTriMidPoint(0, -0.28125); // -9./32.
446  ir->AddTriPoints3(1, 0.2, 25./96.);
447  return ir;
448 
449  case 4: // 6 point - 4 degree
450  TriangleIntRules[4] = ir = new IntegrationRule(6);
451  ir->AddTriPoints3(0, 0.091576213509770743460, 0.054975871827660933819);
452  ir->AddTriPoints3(3, 0.44594849091596488632, 0.11169079483900573285);
453  return ir;
454 
455  case 5: // 7 point - 5 degree
456  TriangleIntRules[5] = ir = new IntegrationRule(7);
457  ir->AddTriMidPoint(0, 0.1125);
458  ir->AddTriPoints3(1, 0.10128650732345633880, 0.062969590272413576298);
459  ir->AddTriPoints3(4, 0.47014206410511508977, 0.066197076394253090369);
460  return ir;
461 
462  case 6: // 12 point - 6 degree
463  TriangleIntRules[6] = ir = new IntegrationRule(12);
464  ir->AddTriPoints3(0, 0.063089014491502228340, 0.025422453185103408460);
465  ir->AddTriPoints3(3, 0.24928674517091042129, 0.058393137863189683013);
466  ir->AddTriPoints6(6, 0.053145049844816947353, 0.31035245103378440542,
467  0.041425537809186787597);
468  return ir;
469 
470  case 7: // 12 point - degree 7
471  TriangleIntRules[7] = ir = new IntegrationRule(12);
472  ir->AddTriPoints3R(0, 0.062382265094402118174, 0.067517867073916085443,
473  0.026517028157436251429);
474  ir->AddTriPoints3R(3, 0.055225456656926611737, 0.32150249385198182267,
475  0.043881408714446055037);
476  // slightly better with explicit 3rd area coordinate
477  ir->AddTriPoints3R(6, 0.034324302945097146470, 0.66094919618673565761,
478  0.30472650086816719592, 0.028775042784981585738);
479  ir->AddTriPoints3R(9, 0.51584233435359177926, 0.27771616697639178257,
480  0.20644149867001643817, 0.067493187009802774463);
481  return ir;
482 
483  case 8: // 16 point - 8 degree
484  TriangleIntRules[8] = ir = new IntegrationRule(16);
485  ir->AddTriMidPoint(0, 0.0721578038388935841255455552445323);
486  ir->AddTriPoints3(1, 0.170569307751760206622293501491464,
487  0.0516086852673591251408957751460645);
488  ir->AddTriPoints3(4, 0.0505472283170309754584235505965989,
489  0.0162292488115990401554629641708902);
490  ir->AddTriPoints3(7, 0.459292588292723156028815514494169,
491  0.0475458171336423123969480521942921);
492  ir->AddTriPoints6(10, 0.008394777409957605337213834539296,
493  0.263112829634638113421785786284643,
494  0.0136151570872174971324223450369544);
495  return ir;
496 
497  case 9: // 19 point - 9 degree
498  TriangleIntRules[9] = ir = new IntegrationRule(19);
499  ir->AddTriMidPoint(0, 0.0485678981413994169096209912536443);
500  ir->AddTriPoints3b(1, 0.020634961602524744433,
501  0.0156673501135695352684274156436046);
502  ir->AddTriPoints3b(4, 0.12582081701412672546,
503  0.0389137705023871396583696781497019);
504  ir->AddTriPoints3(7, 0.188203535619032730240961280467335,
505  0.0398238694636051265164458871320226);
506  ir->AddTriPoints3(10, 0.0447295133944527098651065899662763,
507  0.0127888378293490156308393992794999);
508  ir->AddTriPoints6(13, 0.0368384120547362836348175987833851,
509  0.2219629891607656956751025276931919,
510  0.0216417696886446886446886446886446);
511  return ir;
512 
513  case 10: // 25 point - 10 degree
514  TriangleIntRules[10] = ir = new IntegrationRule(25);
515  ir->AddTriMidPoint(0, 0.0454089951913767900476432975500142);
516  ir->AddTriPoints3b(1, 0.028844733232685245264984935583748,
517  0.0183629788782333523585030359456832);
518  ir->AddTriPoints3(4, 0.109481575485037054795458631340522,
519  0.0226605297177639673913028223692986);
520  ir->AddTriPoints6(7, 0.141707219414879954756683250476361,
521  0.307939838764120950165155022930631,
522  0.0363789584227100543021575883096803);
523  ir->AddTriPoints6(13, 0.025003534762686386073988481007746,
524  0.246672560639902693917276465411176,
525  0.0141636212655287424183685307910495);
526  ir->AddTriPoints6(19, 0.0095408154002994575801528096228873,
527  0.0668032510122002657735402127620247,
528  4.71083348186641172996373548344341E-03);
529  return ir;
530 
531  case 11: // 28 point -- 11 degree
532  TriangleIntRules[11] = ir = new IntegrationRule(28);
533  ir->AddTriPoints6(0, 0.0,
534  0.141129718717363295960826061941652,
535  3.68119189165027713212944752369032E-03);
536  ir->AddTriMidPoint(6, 0.0439886505811161193990465846607278);
537  ir->AddTriPoints3(7, 0.0259891409282873952600324854988407,
538  4.37215577686801152475821439991262E-03);
539  ir->AddTriPoints3(10, 0.0942875026479224956305697762754049,
540  0.0190407859969674687575121697178070);
541  ir->AddTriPoints3b(13, 0.010726449965572372516734795387128,
542  9.42772402806564602923839129555767E-03);
543  ir->AddTriPoints3(16, 0.207343382614511333452934024112966,
544  0.0360798487723697630620149942932315);
545  ir->AddTriPoints3b(19, 0.122184388599015809877869236727746,
546  0.0346645693527679499208828254519072);
547  ir->AddTriPoints6(22, 0.0448416775891304433090523914688007,
548  0.2772206675282791551488214673424523,
549  0.0205281577146442833208261574536469);
550  return ir;
551 
552  case 12: // 33 point - 12 degree
553  TriangleIntRules[12] = ir = new IntegrationRule(33);
554  ir->AddTriPoints3b(0, 2.35652204523900E-02, 1.28655332202275E-02);
555  ir->AddTriPoints3b(3, 1.20551215411079E-01, 2.18462722690190E-02);
556  ir->AddTriPoints3(6, 2.71210385012116E-01, 3.14291121089425E-02);
557  ir->AddTriPoints3(9, 1.27576145541586E-01, 1.73980564653545E-02);
558  ir->AddTriPoints3(12, 2.13173504532100E-02, 3.08313052577950E-03);
559  ir->AddTriPoints6(15, 1.15343494534698E-01, 2.75713269685514E-01,
560  2.01857788831905E-02);
561  ir->AddTriPoints6(21, 2.28383322222570E-02, 2.81325580989940E-01,
562  1.11783866011515E-02);
563  ir->AddTriPoints6(27, 2.57340505483300E-02, 1.16251915907597E-01,
564  8.65811555432950E-03);
565  return ir;
566 
567  case 13: // 37 point - 13 degree
568  TriangleIntRules[13] = ir = new IntegrationRule(37);
569  ir->AddTriPoints3b(0, 0.0,
570  2.67845189554543044455908674650066E-03);
571  ir->AddTriMidPoint(3, 0.0293480398063595158995969648597808);
572  ir->AddTriPoints3(4, 0.0246071886432302181878499494124643,
573  3.92538414805004016372590903990464E-03);
574  ir->AddTriPoints3b(7, 0.159382493797610632566158925635800,
575  0.0253344765879434817105476355306468);
576  ir->AddTriPoints3(10, 0.227900255506160619646298948153592,
577  0.0250401630452545330803738542916538);
578  ir->AddTriPoints3(13, 0.116213058883517905247155321839271,
579  0.0158235572961491595176634480481793);
580  ir->AddTriPoints3b(16, 0.046794039901841694097491569577008,
581  0.0157462815379843978450278590138683);
582  ir->AddTriPoints6(19, 0.0227978945382486125477207592747430,
583  0.1254265183163409177176192369310890,
584  7.90126610763037567956187298486575E-03);
585  ir->AddTriPoints6(25, 0.0162757709910885409437036075960413,
586  0.2909269114422506044621801030055257,
587  7.99081889046420266145965132482933E-03);
588  ir->AddTriPoints6(31, 0.0897330604516053590796290561145196,
589  0.2723110556841851025078181617634414,
590  0.0182757511120486476280967518782978);
591  return ir;
592 
593  case 14: // 42 point - 14 degree
594  TriangleIntRules[14] = ir = new IntegrationRule(42);
595  ir->AddTriPoints3b(0, 2.20721792756430E-02, 1.09417906847145E-02);
596  ir->AddTriPoints3b(3, 1.64710561319092E-01, 1.63941767720625E-02);
597  ir->AddTriPoints3(6, 2.73477528308839E-01, 2.58870522536460E-02);
598  ir->AddTriPoints3(9, 1.77205532412543E-01, 2.10812943684965E-02);
599  ir->AddTriPoints3(12, 6.17998830908730E-02, 7.21684983488850E-03);
600  ir->AddTriPoints3(15, 1.93909612487010E-02, 2.46170180120000E-03);
601  ir->AddTriPoints6(18, 5.71247574036480E-02, 1.72266687821356E-01,
602  1.23328766062820E-02);
603  ir->AddTriPoints6(24, 9.29162493569720E-02, 3.36861459796345E-01,
604  1.92857553935305E-02);
605  ir->AddTriPoints6(30, 1.46469500556540E-02, 2.98372882136258E-01,
606  7.21815405676700E-03);
607  ir->AddTriPoints6(36, 1.26833093287200E-03, 1.18974497696957E-01,
608  2.50511441925050E-03);
609  return ir;
610 
611  case 15: // 54 point - 15 degree
612  TriangleIntRules[15] = ir = new IntegrationRule(54);
613  ir->AddTriPoints3b(0, 0.0834384072617499333, 0.016330909424402645);
614  ir->AddTriPoints3b(3, 0.192779070841738867, 0.01370640901568218);
615  ir->AddTriPoints3(6, 0.293197167913025367, 0.01325501829935165);
616  ir->AddTriPoints3(9, 0.146467786942772933, 0.014607981068243055);
617  ir->AddTriPoints3(12, 0.0563628676656034333, 0.005292304033121995);
618  ir->AddTriPoints3(15, 0.0165751268583703333, 0.0018073215320460175);
619  ir->AddTriPoints6(18, 0.0099122033092248, 0.239534554154794445,
620  0.004263874050854718);
621  ir->AddTriPoints6(24, 0.015803770630228, 0.404878807318339958,
622  0.006958088258345965);
623  ir->AddTriPoints6(30, 0.00514360881697066667, 0.0950021131130448885,
624  0.0021459664703674175);
625  ir->AddTriPoints6(36, 0.0489223257529888, 0.149753107322273969,
626  0.008117664640887445);
627  ir->AddTriPoints6(42, 0.0687687486325192, 0.286919612441334979,
628  0.012803670460631195);
629  ir->AddTriPoints6(48, 0.1684044181246992, 0.281835668099084562,
630  0.016544097765822835);
631  return ir;
632 
633  case 16: // 61 point - 17 degree (used for 16 as well)
634  case 17:
635  TriangleIntRules[16] =
636  TriangleIntRules[17] = ir = new IntegrationRule(61);
637  ir->AddTriMidPoint(0, 1.67185996454015E-02);
638  ir->AddTriPoints3b(1, 5.65891888645200E-03, 2.54670772025350E-03);
639  ir->AddTriPoints3b(4, 3.56473547507510E-02, 7.33543226381900E-03);
640  ir->AddTriPoints3b(7, 9.95200619584370E-02, 1.21754391768360E-02);
641  ir->AddTriPoints3b(10, 1.99467521245206E-01, 1.55537754344845E-02);
642  ir->AddTriPoints3 (13, 2.52141267970953E-01, 1.56285556093100E-02);
643  ir->AddTriPoints3 (16, 1.62047004658461E-01, 1.24078271698325E-02);
644  ir->AddTriPoints3 (19, 7.58758822607460E-02, 7.02803653527850E-03);
645  ir->AddTriPoints3 (22, 1.56547269678220E-02, 1.59733808688950E-03);
646  ir->AddTriPoints6 (25, 1.01869288269190E-02, 3.34319867363658E-01,
647  4.05982765949650E-03);
648  ir->AddTriPoints6 (31, 1.35440871671036E-01, 2.92221537796944E-01,
649  1.34028711415815E-02);
650  ir->AddTriPoints6 (37, 5.44239242905830E-02, 3.19574885423190E-01,
651  9.22999660541100E-03);
652  ir->AddTriPoints6 (43, 1.28685608336370E-02, 1.90704224192292E-01,
653  4.23843426716400E-03);
654  ir->AddTriPoints6 (49, 6.71657824135240E-02, 1.80483211648746E-01,
655  9.14639838501250E-03);
656  ir->AddTriPoints6 (55, 1.46631822248280E-02, 8.07113136795640E-02,
657  3.33281600208250E-03);
658  return ir;
659 
660  case 18: // 73 point - 19 degree (used for 18 as well)
661  case 19:
662  TriangleIntRules[18] =
663  TriangleIntRules[19] = ir = new IntegrationRule(73);
664  ir->AddTriMidPoint(0, 0.0164531656944595);
665  ir->AddTriPoints3b(1, 0.020780025853987, 0.005165365945636);
666  ir->AddTriPoints3b(4, 0.090926214604215, 0.011193623631508);
667  ir->AddTriPoints3b(7, 0.197166638701138, 0.015133062934734);
668  ir->AddTriPoints3 (10, 0.255551654403098, 0.015245483901099);
669  ir->AddTriPoints3 (13, 0.17707794215213, 0.0120796063708205);
670  ir->AddTriPoints3 (16, 0.110061053227952, 0.0080254017934005);
671  ir->AddTriPoints3 (19, 0.05552862425184, 0.004042290130892);
672  ir->AddTriPoints3 (22, 0.012621863777229, 0.0010396810137425);
673  ir->AddTriPoints6 (25, 0.003611417848412, 0.395754787356943,
674  0.0019424384524905);
675  ir->AddTriPoints6 (31, 0.13446675453078, 0.307929983880436,
676  0.012787080306011);
677  ir->AddTriPoints6 (37, 0.014446025776115, 0.26456694840652,
678  0.004440451786669);
679  ir->AddTriPoints6 (43, 0.046933578838178, 0.358539352205951,
680  0.0080622733808655);
681  ir->AddTriPoints6 (49, 0.002861120350567, 0.157807405968595,
682  0.0012459709087455);
683  ir->AddTriPoints6 (55, 0.075050596975911, 0.223861424097916,
684  0.0091214200594755);
685  ir->AddTriPoints6 (61, 0.03464707481676, 0.142421601113383,
686  0.0051292818680995);
687  ir->AddTriPoints6 (67, 0.065494628082938, 0.010161119296278,
688  0.001899964427651);
689  return ir;
690 
691  case 20: // 85 point - 20 degree
692  TriangleIntRules[20] = ir = new IntegrationRule(85);
693  ir->AddTriMidPoint(0, 0.01380521349884976);
694  ir->AddTriPoints3b(1, 0.001500649324429, 0.00088951477366337);
695  ir->AddTriPoints3b(4, 0.0941397519389508667, 0.010056199056980585);
696  ir->AddTriPoints3b(7, 0.2044721240895264, 0.013408923629665785);
697  ir->AddTriPoints3(10, 0.264500202532787333, 0.012261566900751005);
698  ir->AddTriPoints3(13, 0.211018964092076767, 0.008197289205347695);
699  ir->AddTriPoints3(16, 0.107735607171271333, 0.0073979536993248);
700  ir->AddTriPoints3(19, 0.0390690878378026667, 0.0022896411388521255);
701  ir->AddTriPoints3(22, 0.0111743797293296333, 0.0008259132577881085);
702  ir->AddTriPoints6(25, 0.00534961818733726667, 0.0635496659083522206,
703  0.001174585454287792);
704  ir->AddTriPoints6(31, 0.00795481706619893333, 0.157106918940706982,
705  0.0022329628770908965);
706  ir->AddTriPoints6(37, 0.0104223982812638, 0.395642114364374018,
707  0.003049783403953986);
708  ir->AddTriPoints6(43, 0.0109644147961233333, 0.273167570712910522,
709  0.0034455406635941015);
710  ir->AddTriPoints6(49, 0.0385667120854623333, 0.101785382485017108,
711  0.0039987375362390815);
712  ir->AddTriPoints6(55, 0.0355805078172182, 0.446658549176413815,
713  0.003693067142668012);
714  ir->AddTriPoints6(61, 0.0496708163627641333, 0.199010794149503095,
715  0.00639966593932413);
716  ir->AddTriPoints6(67, 0.0585197250843317333, 0.3242611836922827,
717  0.008629035587848275);
718  ir->AddTriPoints6(73, 0.121497787004394267, 0.208531363210132855,
719  0.009336472951467735);
720  ir->AddTriPoints6(79, 0.140710844943938733, 0.323170566536257485,
721  0.01140911202919763);
722  return ir;
723 
724  case 21: // 126 point - 25 degree (used also for degrees from 21 to 24)
725  case 22:
726  case 23:
727  case 24:
728  case 25:
729  TriangleIntRules[21] =
730  TriangleIntRules[22] =
731  TriangleIntRules[23] =
732  TriangleIntRules[24] =
733  TriangleIntRules[25] = ir = new IntegrationRule(126);
734  ir->AddTriPoints3b(0, 0.0279464830731742, 0.0040027909400102085);
735  ir->AddTriPoints3b(3, 0.131178601327651467, 0.00797353841619525);
736  ir->AddTriPoints3b(6, 0.220221729512072267, 0.006554570615397765);
737  ir->AddTriPoints3 (9, 0.298443234019804467, 0.00979150048281781);
738  ir->AddTriPoints3(12, 0.2340441723373718, 0.008235442720768635);
739  ir->AddTriPoints3(15, 0.151468334609017567, 0.00427363953704605);
740  ir->AddTriPoints3(18, 0.112733893545993667, 0.004080942928613246);
741  ir->AddTriPoints3(21, 0.0777156920915263, 0.0030605732699918895);
742  ir->AddTriPoints3(24, 0.034893093614297, 0.0014542491324683325);
743  ir->AddTriPoints3(27, 0.00725818462093236667, 0.00034613762283099815);
744  ir->AddTriPoints6(30, 0.0012923527044422, 0.227214452153364077,
745  0.0006241445996386985);
746  ir->AddTriPoints6(36, 0.0053997012721162, 0.435010554853571706,
747  0.001702376454401511);
748  ir->AddTriPoints6(42, 0.006384003033975, 0.320309599272204437,
749  0.0016798271630320255);
750  ir->AddTriPoints6(48, 0.00502821150199306667, 0.0917503222800051889,
751  0.000858078269748377);
752  ir->AddTriPoints6(54, 0.00682675862178186667, 0.0380108358587243835,
753  0.000740428158357803);
754  ir->AddTriPoints6(60, 0.0100161996399295333, 0.157425218485311668,
755  0.0017556563053643425);
756  ir->AddTriPoints6(66, 0.02575781317339, 0.239889659778533193,
757  0.003696775074853242);
758  ir->AddTriPoints6(72, 0.0302278981199158, 0.361943118126060531,
759  0.003991543738688279);
760  ir->AddTriPoints6(78, 0.0305049901071620667, 0.0835519609548285602,
761  0.0021779813065790205);
762  ir->AddTriPoints6(84, 0.0459565473625693333, 0.148443220732418205,
763  0.003682528350708916);
764  ir->AddTriPoints6(90, 0.0674428005402775333, 0.283739708727534955,
765  0.005481786423209775);
766  ir->AddTriPoints6(96, 0.0700450914159106, 0.406899375118787573,
767  0.00587498087177056);
768  ir->AddTriPoints6(102, 0.0839115246401166, 0.194113987024892542,
769  0.005007800356899285);
770  ir->AddTriPoints6(108, 0.120375535677152667, 0.32413434700070316,
771  0.00665482039381434);
772  ir->AddTriPoints6(114, 0.148066899157366667, 0.229277483555980969,
773  0.00707722325261307);
774  ir->AddTriPoints6(120, 0.191771865867325067, 0.325618122595983752,
775  0.007440689780584005);
776  return ir;
777 
778  default:
779  // Grundmann-Moller rules
780  int i = (Order / 2) * 2 + 1; // Get closest odd # >= Order
781  AllocIntRule(TriangleIntRules, i);
782  TriangleIntRules[i-1] = TriangleIntRules[i] = ir = new IntegrationRule;
783  ir->GrundmannMollerSimplexRule(i/2,2);
784  return ir;
785  }
786 }
787 
788 // Integration rules for unit square
789 IntegrationRule *IntegrationRules::SquareIntegrationRule(int Order)
790 {
791  int i = (Order / 2) * 2 + 1; // Get closest odd # >= Order
792 
793  if (!HaveIntRule(SegmentIntRules, i))
794  {
795  SegmentIntegrationRule(i);
796  }
797  AllocIntRule(SquareIntRules, i);
798  SquareIntRules[i-1] =
799  SquareIntRules[i] =
800  new IntegrationRule(*SegmentIntRules[i], *SegmentIntRules[i]);
801  return SquareIntRules[i];
802 }
803 
806 IntegrationRule *IntegrationRules::TetrahedronIntegrationRule(int Order)
807 {
808  IntegrationRule *ir;
809 
810  // assuming that orders <= 9 are pre-allocated
811  switch (Order)
812  {
813  case 0: // 1 point - degree 1
814  case 1:
815  TetrahedronIntRules[0] =
816  TetrahedronIntRules[1] = ir = new IntegrationRule(1);
817  ir->AddTetMidPoint(0, 1./6.);
818  return ir;
819 
820  case 2: // 4 points - degree 2
821  TetrahedronIntRules[2] = ir = new IntegrationRule(4);
822  // ir->AddTetPoints4(0, 0.13819660112501051518, 1./24.);
823  ir->AddTetPoints4b(0, 0.58541019662496845446, 1./24.);
824  return ir;
825 
826  case 3: // 5 points - degree 3 (negative weight)
827  TetrahedronIntRules[3] = ir = new IntegrationRule(5);
828  ir->AddTetMidPoint(0, -2./15.);
829  ir->AddTetPoints4b(1, 0.5, 0.075);
830  return ir;
831 
832  case 4: // 11 points - degree 4 (negative weight)
833  TetrahedronIntRules[4] = ir = new IntegrationRule(11);
834  ir->AddTetPoints4(0, 1./14., 343./45000.);
835  ir->AddTetMidPoint(4, -74./5625.);
836  ir->AddTetPoints6(5, 0.10059642383320079500, 28./1125.);
837  return ir;
838 
839  case 5: // 14 points - degree 5
840  TetrahedronIntRules[5] = ir = new IntegrationRule(14);
841  ir->AddTetPoints6(0, 0.045503704125649649492,
842  7.0910034628469110730E-03);
843  ir->AddTetPoints4(6, 0.092735250310891226402, 0.012248840519393658257);
844  ir->AddTetPoints4b(10, 0.067342242210098170608,
845  0.018781320953002641800);
846  return ir;
847 
848  case 6: // 24 points - degree 6
849  TetrahedronIntRules[6] = ir = new IntegrationRule(24);
850  ir->AddTetPoints4(0, 0.21460287125915202929,
851  6.6537917096945820166E-03);
852  ir->AddTetPoints4(4, 0.040673958534611353116,
853  1.6795351758867738247E-03);
854  ir->AddTetPoints4b(8, 0.032986329573173468968,
855  9.2261969239424536825E-03);
856  ir->AddTetPoints12(12, 0.063661001875017525299, 0.26967233145831580803,
857  8.0357142857142857143E-03);
858  return ir;
859 
860  case 7: // 31 points - degree 7 (negative weight)
861  TetrahedronIntRules[7] = ir = new IntegrationRule(31);
862  ir->AddTetPoints6(0, 0.0, 9.7001763668430335097E-04);
863  ir->AddTetMidPoint(6, 0.018264223466108820291);
864  ir->AddTetPoints4(7, 0.078213192330318064374, 0.010599941524413686916);
865  ir->AddTetPoints4(11, 0.12184321666390517465,
866  -0.062517740114331851691);
867  ir->AddTetPoints4b(15, 2.3825066607381275412E-03,
868  4.8914252630734993858E-03);
869  ir->AddTetPoints12(19, 0.1, 0.2, 0.027557319223985890653);
870  return ir;
871 
872  case 8: // 43 points - degree 8 (negative weight)
873  TetrahedronIntRules[8] = ir = new IntegrationRule(43);
874  ir->AddTetPoints4(0, 5.7819505051979972532E-03,
875  1.6983410909288737984E-04);
876  ir->AddTetPoints4(4, 0.082103588310546723091,
877  1.9670333131339009876E-03);
878  ir->AddTetPoints12(8, 0.036607749553197423679, 0.19048604193463345570,
879  2.1405191411620925965E-03);
880  ir->AddTetPoints6(20, 0.050532740018894224426,
881  4.5796838244672818007E-03);
882  ir->AddTetPoints12(26, 0.22906653611681113960, 0.035639582788534043717,
883  5.7044858086819185068E-03);
884  ir->AddTetPoints4(38, 0.20682993161067320408, 0.014250305822866901248);
885  ir->AddTetMidPoint(42, -0.020500188658639915841);
886  return ir;
887 
888  case 9: // orders 9 and higher -- Grundmann-Moller rules
889  TetrahedronIntRules[9] = ir = new IntegrationRule;
890  ir->GrundmannMollerSimplexRule(4,3);
891  return ir;
892 
893  default: // Grundmann-Moller rules
894  int i = (Order / 2) * 2 + 1; // Get closest odd # >= Order
895  AllocIntRule(TetrahedronIntRules, i);
896  TetrahedronIntRules[i-1] =
897  TetrahedronIntRules[i] = ir = new IntegrationRule;
898  ir->GrundmannMollerSimplexRule(i/2,3);
899  return ir;
900  }
901 }
902 
903 // Integration rules for reference cube
904 IntegrationRule *IntegrationRules::CubeIntegrationRule(int Order)
905 {
906  int i = (Order / 2) * 2 + 1; // Get closest odd # >= Order
907 
908  if (!HaveIntRule(SegmentIntRules, i))
909  {
910  SegmentIntegrationRule(i);
911  }
912  AllocIntRule(CubeIntRules, i);
913  CubeIntRules[i-1] =
914  CubeIntRules[i] =
915  new IntegrationRule(*SegmentIntRules[i], *SegmentIntRules[i],
916  *SegmentIntRules[i]);
917  return CubeIntRules[i];
918 }
919 
920 }
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:228
int Size() const
Logical size of the array.
Definition: array.hpp:109
Class for integration rule.
Definition: intrules.hpp:83
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:231
IntegrationRules IntRules(0)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:292
void mfem_error(const char *msg)
Definition: error.cpp:23
Class for integration point with weight.
Definition: intrules.hpp:25
IntegrationRules RefinedIntRules(1)
A global object with all refined integration rules.
Definition: intrules.hpp:295