MFEM  v3.3
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 #ifdef MFEM_USE_MPFR
22 #include <mpfr.h>
23 #endif
24 
25 using namespace std;
26 
27 namespace mfem
28 {
29 
30 IntegrationRule::IntegrationRule(IntegrationRule &irx, IntegrationRule &iry)
31 {
32  int i, j, nx, ny;
33 
34  nx = irx.GetNPoints();
35  ny = iry.GetNPoints();
36  SetSize(nx * ny);
37 
38  for (j = 0; j < ny; j++)
39  {
40  IntegrationPoint &ipy = iry.IntPoint(j);
41  for (i = 0; i < nx; i++)
42  {
43  IntegrationPoint &ipx = irx.IntPoint(i);
44  IntegrationPoint &ip = IntPoint(j*nx+i);
45 
46  ip.x = ipx.x;
47  ip.y = ipy.x;
48  ip.weight = ipx.weight * ipy.weight;
49  }
50  }
51 }
52 
53 IntegrationRule::IntegrationRule(IntegrationRule &irx, IntegrationRule &iry,
54  IntegrationRule &irz)
55 {
56  const int nx = irx.GetNPoints();
57  const int ny = iry.GetNPoints();
58  const int nz = irz.GetNPoints();
59  SetSize(nx*ny*nz);
60 
61  for (int iz = 0; iz < nz; ++iz)
62  {
63  IntegrationPoint &ipz = irz.IntPoint(iz);
64  for (int iy = 0; iy < ny; ++iy)
65  {
66  IntegrationPoint &ipy = iry.IntPoint(iy);
67  for (int ix = 0; ix < nx; ++ix)
68  {
69  IntegrationPoint &ipx = irx.IntPoint(ix);
70  IntegrationPoint &ip = IntPoint(iz*nx*ny + iy*nx + ix);
71 
72  ip.x = ipx.x;
73  ip.y = ipy.x;
74  ip.z = ipz.x;
75  ip.weight = ipx.weight*ipy.weight*ipz.weight;
76  }
77  }
78  }
79 }
80 
81 void IntegrationRule::GrundmannMollerSimplexRule(int s, int n)
82 {
83  // for pow on older compilers
84  using std::pow;
85  const int d = 2*s + 1;
86  Vector fact(d + n + 1);
87  Array<int> beta(n), sums(n);
88 
89  fact(0) = 1.;
90  for (int i = 1; i < fact.Size(); i++)
91  {
92  fact(i) = fact(i - 1)*i;
93  }
94 
95  // number of points is \binom{n + s + 1}{n + 1}
96  int np = 1, f = 1;
97  for (int i = 0; i <= n; i++)
98  {
99  np *= (s + i + 1), f *= (i + 1);
100  }
101  np /= f;
102  SetSize(np);
103 
104  int pt = 0;
105  for (int i = 0; i <= s; i++)
106  {
107  double weight;
108 
109  weight = pow(2., -2*s)*pow(static_cast<double>(d + n - 2*i),
110  d)/fact(i)/fact(d + n - i);
111  if (i%2)
112  {
113  weight = -weight;
114  }
115 
116  // loop over all beta : beta_0 + ... + beta_{n-1} <= s - i
117  int k = s - i;
118  beta = 0;
119  sums = 0;
120  while (true)
121  {
122  IntegrationPoint &ip = IntPoint(pt++);
123  ip.weight = weight;
124  ip.x = double(2*beta[0] + 1)/(d + n - 2*i);
125  ip.y = double(2*beta[1] + 1)/(d + n - 2*i);
126  if (n == 3)
127  {
128  ip.z = double(2*beta[2] + 1)/(d + n - 2*i);
129  }
130 
131  int j = 0;
132  while (sums[j] == k)
133  {
134  beta[j++] = 0;
135  if (j == n)
136  {
137  goto done_beta;
138  }
139  }
140  beta[j]++;
141  sums[j]++;
142  for (j--; j >= 0; j--)
143  {
144  sums[j] = sums[j+1];
145  }
146  }
147  done_beta:
148  ;
149  }
150 }
151 
152 
153 #ifdef MFEM_USE_MPFR
154 
155 // Class for computing hi-precision (HP) quadrature in 1D
156 class HP_Quadrature1D
157 {
158 protected:
159  mpfr_t pi, z, pp, p1, p2, p3, dz, w, rtol;
160 
161 public:
162  static const mpfr_rnd_t rnd = GMP_RNDN;
163  static const int default_prec = 128;
164 
165  // prec = MPFR precision in bits
166  HP_Quadrature1D(const int prec = default_prec)
167  {
168  mpfr_inits2(prec, pi, z, pp, p1, p2, p3, dz, w, rtol, (mpfr_ptr) 0);
169  mpfr_const_pi(pi, rnd);
170  mpfr_set_si_2exp(rtol, 1, -32, rnd); // 2^(-32) < 2.33e-10
171  }
172 
173  // set rtol = 2^exponent
174  // this is a tolerance for the last correction of x_i in Newton's algorithm;
175  // this gives roughly rtol^2 accuracy for the final x_i.
176  void SetRelTol(const int exponent = -32)
177  {
178  mpfr_set_si_2exp(rtol, 1, exponent, rnd);
179  }
180 
181  // n - number of quadrature points
182  // k - index of the point to compute, 0 <= k < n
183  // see also: QuadratureFunctions1D::GaussLegendre
184  void ComputeGaussLegendrePoint(const int n, const int k)
185  {
186  MFEM_ASSERT(n > 0 && 0 <= k && k < n, "invalid n = " << n
187  << " and/or k = " << k);
188 
189  int i = (k < (n+1)/2) ? k+1 : n-k;
190 
191  // Initial guess for the x-coordinate:
192  // set z = cos(pi * (i - 0.25) / (n + 0.5)) =
193  // = sin(pi * ((n+1-2*i) / (2*n+1)))
194  mpfr_set_si(z, n+1-2*i, rnd);
195  mpfr_div_si(z, z, 2*n+1, rnd);
196  mpfr_mul(z, z, pi, rnd);
197  mpfr_sin(z, z, rnd);
198 
199  bool done = false;
200  while (1)
201  {
202  mpfr_set_si(p2, 1, rnd);
203  mpfr_set(p1, z, rnd);
204  for (int j = 2; j <= n; j++)
205  {
206  mpfr_set(p3, p2, rnd);
207  mpfr_set(p2, p1, rnd);
208  // p1 = ((2 * j - 1) * z * p2 - (j - 1) * p3) / j;
209  mpfr_mul_si(p1, z, 2*j-1, rnd);
210  mpfr_mul_si(p3, p3, j-1, rnd);
211  mpfr_fms(p1, p1, p2, p3, rnd);
212  mpfr_div_si(p1, p1, j, rnd);
213  }
214  // p1 is Legendre polynomial
215 
216  // derivative of the Legendre polynomial:
217  // pp = n * (z*p1-p2) / (z*z - 1);
218  mpfr_fms(pp, z, p1, p2, rnd);
219  mpfr_mul_si(pp, pp, n, rnd);
220  mpfr_sqr(p2, z, rnd);
221  mpfr_sub_si(p2, p2, 1, rnd);
222  mpfr_div(pp, pp, p2, rnd);
223 
224  if (done) { break; }
225 
226  // set delta_z: dz = p1/pp;
227  mpfr_div(dz, p1, pp, rnd);
228  // compute absolute tolerance: atol = rtol*(1-z)
229  mpfr_t &atol = w;
230  mpfr_si_sub(atol, 1, z, rnd);
231  mpfr_mul(atol, atol, rtol, rnd);
232  if (mpfr_cmpabs(dz, atol) <= 0)
233  {
234  done = true;
235  // continue the computation: get pp at the new point, then exit
236  }
237  // update z = z - dz
238  mpfr_sub(z, z, dz, rnd);
239  }
240 
241  // map z to (0,1): z = (1 - z)/2
242  mpfr_si_sub(z, 1, z, rnd);
243  mpfr_div_2si(z, z, 1, rnd);
244 
245  // weight: w = 1/(4*z*(1 - z)*pp*pp)
246  mpfr_sqr(w, pp, rnd);
247  mpfr_mul_2si(w, w, 2, rnd);
248  mpfr_mul(w, w, z, rnd);
249  mpfr_si_sub(p1, 1, z, rnd); // p1 = 1-z
250  mpfr_mul(w, w, p1, rnd);
251  mpfr_si_div(w, 1, w, rnd);
252 
253  if (k >= (n+1)/2) { mpfr_swap(z, p1); }
254  }
255 
256  // n - number of quadrature points
257  // k - index of the point to compute, 0 <= k < n
258  // see also: QuadratureFunctions1D::GaussLobatto
259  void ComputeGaussLobattoPoint(const int n, const int k)
260  {
261  MFEM_ASSERT(n > 1 && 0 <= k && k < n, "invalid n = " << n
262  << " and/or k = " << k);
263 
264  int i = (k < (n+1)/2) ? k : n-1-k;
265 
266  if (i == 0)
267  {
268  mpfr_set_si(z, 0, rnd);
269  mpfr_set_si(p1, 1, rnd);
270  mpfr_set_si(w, n*(n-1), rnd);
271  mpfr_si_div(w, 1, w, rnd); // weight = 1/(n*(n-1))
272  return;
273  }
274  // initial guess is the corresponding Chebyshev point, z:
275  // z = -cos(pi * i/(n-1)) = sin(pi * (2*i-n+1)/(2*n-2))
276  mpfr_set_si(z, 2*i-n+1, rnd);
277  mpfr_div_si(z, z, 2*(n-1), rnd);
278  mpfr_mul(z, pi, z, rnd);
279  mpfr_sin(z, z, rnd);
280  bool done = false;
281  for (int iter = 0 ; true ; ++iter)
282  {
283  // build Legendre polynomials, up to P_{n}(z)
284  mpfr_set_si(p1, 1, rnd);
285  mpfr_set(p2, z, rnd);
286 
287  for (int l = 1 ; l < (n-1) ; ++l)
288  {
289  // P_{l+1}(x) = [ (2*l+1)*x*P_l(x) - l*P_{l-1}(x) ]/(l+1)
290  mpfr_mul_si(p1, p1, l, rnd);
291  mpfr_mul_si(p3, z, 2*l+1, rnd);
292  mpfr_fms(p3, p3, p2, p1, rnd);
293  mpfr_div_si(p3, p3, l+1, rnd);
294 
295  mpfr_set(p1, p2, rnd);
296  mpfr_set(p2, p3, rnd);
297  }
298  if (done) { break; }
299  // compute dz = resid/deriv = (z*p2 - p1) / (n*p2);
300  mpfr_fms(dz, z, p2, p1, rnd);
301  mpfr_mul_si(p3, p2, n, rnd);
302  mpfr_div(dz, dz, p3, rnd);
303  // update: z = z - dz
304  mpfr_sub(z, z, dz, rnd);
305  // compute absolute tolerance: atol = rtol*(1 + z)
306  mpfr_t &atol = w;
307  mpfr_add_si(atol, z, 1, rnd);
308  mpfr_mul(atol, atol, rtol, rnd);
309  // check for convergence
310  if (mpfr_cmpabs(dz, atol) <= 0)
311  {
312  done = true;
313  // continue the computation: get p2 at the new point, then exit
314  }
315  // If the iteration does not converge fast, something is wrong.
316  MFEM_VERIFY(iter < 8, "n = " << n << ", i = " << i
317  << ", dz = " << mpfr_get_d(dz, rnd));
318  }
319  // Map to the interval [0,1] and scale the weights
320  mpfr_add_si(z, z, 1, rnd);
321  mpfr_div_2si(z, z, 1, rnd);
322  // set the symmetric point
323  mpfr_si_sub(p1, 1, z, rnd);
324  // w = 1/[ n*(n-1)*[P_{n-1}(z)]^2 ]
325  mpfr_sqr(w, p2, rnd);
326  mpfr_mul_si(w, w, n*(n-1), rnd);
327  mpfr_si_div(w, 1, w, rnd);
328 
329  if (k >= (n+1)/2) { mpfr_swap(z, p1); }
330  }
331 
332  double GetPoint() const { return mpfr_get_d(z, rnd); }
333  double GetSymmPoint() const { return mpfr_get_d(p1, rnd); }
334  double GetWeight() const { return mpfr_get_d(w, rnd); }
335 
336  const mpfr_t &GetHPPoint() const { return z; }
337  const mpfr_t &GetHPSymmPoint() const { return p1; }
338  const mpfr_t &GetHPWeight() const { return w; }
339 
340  ~HP_Quadrature1D()
341  {
342  mpfr_clears(pi, z, pp, p1, p2, p3, dz, w, rtol, (mpfr_ptr) 0);
343  mpfr_free_cache();
344  }
345 };
346 
347 #endif // MFEM_USE_MPFR
348 
349 
350 void QuadratureFunctions1D::GaussLegendre(const int np, IntegrationRule* ir)
351 {
352  ir->SetSize(np);
353 
354  switch (np)
355  {
356  case 1:
357  ir->IntPoint(0).Set1w(0.5, 1.0);
358  return;
359  case 2:
360  ir->IntPoint(0).Set1w(0.21132486540518711775, 0.5);
361  ir->IntPoint(1).Set1w(0.78867513459481288225, 0.5);
362  return;
363  case 3:
364  ir->IntPoint(0).Set1w(0.11270166537925831148, 5./18.);
365  ir->IntPoint(1).Set1w(0.5, 4./9.);
366  ir->IntPoint(2).Set1w(0.88729833462074168852, 5./18.);
367  return;
368  }
369 
370  const int n = np;
371  const int m = (n+1)/2;
372 
373 #ifndef MFEM_USE_MPFR
374 
375  for (int i = 1; i <= m; i++)
376  {
377  double z = cos(M_PI * (i - 0.25) / (n + 0.5));
378  double pp, p1, dz, xi = 0.;
379  bool done = false;
380  while (1)
381  {
382  double p2 = 1;
383  p1 = z;
384  for (int j = 2; j <= n; j++)
385  {
386  double p3 = p2;
387  p2 = p1;
388  p1 = ((2 * j - 1) * z * p2 - (j - 1) * p3) / j;
389  }
390  // p1 is Legendre polynomial
391 
392  pp = n * (z*p1-p2) / (z*z - 1);
393  if (done) { break; }
394 
395  dz = p1/pp;
396  if (fabs(dz) < 1e-16)
397  {
398  done = true;
399  // map the new point (z-dz) to (0,1):
400  xi = ((1 - z) + dz)/2; // (1 - (z - dz))/2 has bad round-off
401  // continue the computation: get pp at the new point, then exit
402  }
403  // update: z = z - dz
404  z -= dz;
405  }
406 
407  ir->IntPoint(i-1).x = xi;
408  ir->IntPoint(n-i).x = 1 - xi;
409  ir->IntPoint(i-1).weight =
410  ir->IntPoint(n-i).weight = 1./(4*xi*(1 - xi)*pp*pp);
411  }
412 
413 #else // MFEM_USE_MPFR is defined
414 
415  HP_Quadrature1D hp_quad;
416  for (int i = 1; i <= m; i++)
417  {
418  hp_quad.ComputeGaussLegendrePoint(n, i-1);
419 
420  ir->IntPoint(i-1).x = hp_quad.GetPoint();
421  ir->IntPoint(n-i).x = hp_quad.GetSymmPoint();
422  ir->IntPoint(i-1).weight = ir->IntPoint(n-i).weight = hp_quad.GetWeight();
423  }
424 
425 #endif // MFEM_USE_MPFR
426 
427 }
428 
429 void QuadratureFunctions1D::GaussLobatto(const int np, IntegrationRule* ir)
430 {
431  /* An np point Gauss-Lobatto quadrature has (np - 2) free abscissa the other
432  (2) abscissa are the interval endpoints.
433 
434  The interior x_i are the zeros of P'_{np-1}(x). The weights of the
435  interior points on the interval [-1,1] are:
436 
437  w_i = 2/(np*(np-1)*[P_{np-1}(x_i)]^2)
438 
439  The end point weights (on [-1,1]) are: w_{end} = 2/(np*(np-1)).
440 
441  The interior abscissa are found via a nonlinear solve, the initial guess
442  for each point is the corresponding Chebyshev point.
443 
444  After we find all points on the interval [-1,1], we will map and scale the
445  points and weights to the MFEM natural interval [0,1].
446 
447  References:
448  [1] E. E. Lewis and W. F. Millier, "Computational Methods of Neutron
449  Transport", Appendix A
450  [2] the QUADRULE software by John Burkardt,
451  https://people.sc.fsu.edu/~jburkardt/cpp_src/quadrule/quadrule.cpp
452  */
453 
454  ir->SetSize(np);
455  if ( np == 1 )
456  {
457  ir->IntPoint(0).Set1w(0.5, 1.0);
458  }
459  else
460  {
461 
462 #ifndef MFEM_USE_MPFR
463 
464  // endpoints and respective weights
465  ir->IntPoint(0).x = 0.0;
466  ir->IntPoint(np-1).x = 1.0;
467  ir->IntPoint(0).weight = ir->IntPoint(np-1).weight = 1.0/(np*(np-1));
468 
469  // interior points and weights
470  // use symmetry and compute just half of the points
471  for (int i = 1 ; i <= (np-1)/2 ; ++i)
472  {
473  // initial guess is the corresponding Chebyshev point, x_i:
474  // x_i = -cos(\pi * (i / (np-1)))
475  double x_i = std::sin(M_PI * ((double)(i)/(np-1) - 0.5));
476  double z_i = 0., p_l;
477  bool done = false;
478  for (int iter = 0 ; true ; ++iter)
479  {
480  // build Legendre polynomials, up to P_{np}(x_i)
481  double p_lm1 = 1.0;
482  p_l = x_i;
483 
484  for (int l = 1 ; l < (np-1) ; ++l)
485  {
486  // The Legendre polynomials can be built by recursion:
487  // x * P_l(x) = 1/(2*l+1)*[ (l+1)*P_{l+1}(x) + l*P_{l-1} ], i.e.
488  // P_{l+1}(x) = [ (2*l+1)*x*P_l(x) - l*P_{l-1} ]/(l+1)
489  double p_lp1 = ( (2*l + 1)*x_i*p_l - l*p_lm1)/(l + 1);
490 
491  p_lm1 = p_l;
492  p_l = p_lp1;
493  }
494  if (done) { break; }
495  // after this loop, p_l holds P_{np-1}(x_i)
496  // resid = (x^2-1)*P'_{np-1}(x_i)
497  // but use the recurrence relationship
498  // (x^2 -1)P'_l(x) = l*[ x*P_l(x) - P_{l-1}(x) ]
499  // thus, resid = (np-1) * (x_i*p_l - p_lm1)
500 
501  // The derivative of the residual is:
502  // \frac{d}{d x} \left[ (x^2 -1)P'_l(x) ] \right] =
503  // l * (l+1) * P_l(x), with l = np-1,
504  // therefore, deriv = np * (np-1) * p_l;
505 
506  // compute dx = resid/deriv
507  double dx = (x_i*p_l - p_lm1) / (np*p_l);
508  if (std::abs(dx) < 1e-16)
509  {
510  done = true;
511  // Map the point to the interval [0,1]
512  z_i = ((1.0 + x_i) - dx)/2;
513  // continue the computation: get p_l at the new point, then exit
514  }
515  // If the iteration does not converge fast, something is wrong.
516  MFEM_VERIFY(iter < 8, "np = " << np << ", i = " << i
517  << ", dx = " << dx);
518  // update x_i:
519  x_i -= dx;
520  }
521  // Map to the interval [0,1] and scale the weights
522  IntegrationPoint &ip = ir->IntPoint(i);
523  ip.x = z_i;
524  // w_i = (2/[ n*(n-1)*[P_{n-1}(x_i)]^2 ]) / 2
525  ip.weight = (double)(1.0 / (np*(np-1)*p_l*p_l));
526 
527  // set the symmetric point
528  IntegrationPoint &symm_ip = ir->IntPoint(np-1-i);
529  symm_ip.x = 1.0 - z_i;
530  symm_ip.weight = ip.weight;
531  }
532 
533 #else // MFEM_USE_MPFR is defined
534 
535  HP_Quadrature1D hp_quad;
536  // use symmetry and compute just half of the points
537  for (int i = 0 ; i <= (np-1)/2 ; ++i)
538  {
539  hp_quad.ComputeGaussLobattoPoint(np, i);
540  ir->IntPoint(i).x = hp_quad.GetPoint();
541  ir->IntPoint(np-1-i).x = hp_quad.GetSymmPoint();
542  ir->IntPoint(i).weight =
543  ir->IntPoint(np-1-i).weight = hp_quad.GetWeight();
544  }
545 
546 #endif // MFEM_USE_MPFR
547 
548  }
549 }
550 
551 void QuadratureFunctions1D::OpenUniform(const int np, IntegrationRule* ir)
552 {
553  ir->SetSize(np);
554 
555  // The Newton-Cotes quadrature is based on weights that integrate exactly the
556  // interpolatory polynomial through the equally spaced quadrature points.
557  for (int i = 0; i < np ; ++i)
558  {
559  ir->IntPoint(i).x = double(i+1) / double(np + 1);
560  }
561 
562  CalculateUniformWeights(ir, Quadrature1D::OpenUniform);
563 }
564 
565 void QuadratureFunctions1D::ClosedUniform(const int np,
566  IntegrationRule* ir)
567 {
568  ir->SetSize(np);
569  if ( np == 1 ) // allow this case as "closed"
570  {
571  ir->IntPoint(0).Set1w(0.5, 1.0);
572  return;
573  }
574 
575  for (int i = 0; i < np ; ++i)
576  {
577  ir->IntPoint(i).x = double(i) / (np-1);
578  }
579 
580  CalculateUniformWeights(ir, Quadrature1D::ClosedUniform);
581 }
582 
583 void QuadratureFunctions1D::OpenHalfUniform(const int np, IntegrationRule* ir)
584 {
585  ir->SetSize(np);
586 
587  // Open half points: the centers of np uniform intervals
588  for (int i = 0; i < np ; ++i)
589  {
590  ir->IntPoint(i).x = double(2*i+1) / (2*np);
591  }
592 
593  CalculateUniformWeights(ir, Quadrature1D::OpenHalfUniform);
594 }
595 
596 void QuadratureFunctions1D::GivePolyPoints(const int np, double *pts,
597  const int type)
598 {
599  IntegrationRule ir(np);
600 
601  switch (type)
602  {
603  case Quadrature1D::GaussLegendre:
604  {
605  GaussLegendre(np,&ir);
606  break;
607  }
608  case Quadrature1D::GaussLobatto:
609  {
610  GaussLobatto(np, &ir);
611  break;
612  }
613  case Quadrature1D::OpenUniform:
614  {
615  OpenUniform(np,&ir);
616  break;
617  }
618  case Quadrature1D::ClosedUniform:
619  {
620  ClosedUniform(np,&ir);
621  break;
622  }
623  case Quadrature1D::OpenHalfUniform:
624  {
625  OpenHalfUniform(np, &ir);
626  break;
627  }
628  default:
629  {
630  MFEM_ABORT("Asking for an unknown type of 1D Quadrature points, "
631  "type = " << type);
632  }
633  }
634 
635  for (int i = 0 ; i < np ; ++i)
636  {
637  pts[i] = ir.IntPoint(i).x;
638  }
639 }
640 
641 void QuadratureFunctions1D::CalculateUniformWeights(IntegrationRule *ir,
642  const int type)
643 {
644  /* The Lagrange polynomials are:
645  p_i = \prod_{j \neq i} {\frac{x - x_j }{x_i - x_j}}
646 
647  The weight associated with each abscissa is the integral of p_i over
648  [0,1]. To calculate the integral of p_i, we use a Gauss-Legendre
649  quadrature rule. This approach does not suffer from bad round-off/
650  cancellation errors for large number of points.
651  */
652  const int n = ir->Size();
653  switch (n)
654  {
655  case 1:
656  ir->IntPoint(0).weight = 1.;
657  return;
658  case 2:
659  ir->IntPoint(0).weight = .5;
660  ir->IntPoint(1).weight = .5;
661  return;
662  }
663 
664 #ifndef MFEM_USE_MPFR
665 
666  // This algorithm should work for any set of points, not just uniform
667  const IntegrationRule &glob_ir = IntRules.Get(Geometry::SEGMENT, n-1);
668  const int m = glob_ir.GetNPoints();
669  Vector xv(n);
670  for (int j = 0; j < n; j++)
671  {
672  xv(j) = ir->IntPoint(j).x;
673  }
674  Poly_1D::Basis basis(n-1, xv.GetData()); // nodal basis, with nodes at 'xv'
675  Vector w(n);
676  // Integrate all nodal basis functions using 'glob_ir':
677  w = 0.0;
678  for (int i = 0; i < m; i++)
679  {
680  const IntegrationPoint &ip = glob_ir.IntPoint(i);
681  basis.Eval(ip.x, xv);
682  w.Add(ip.weight, xv); // w += ip.weight * xv
683  }
684  for (int j = 0; j < n; j++)
685  {
686  ir->IntPoint(j).weight = w(j);
687  }
688 
689 #else // MFEM_USE_MPFR is defined
690 
691  static const mpfr_rnd_t rnd = HP_Quadrature1D::rnd;
692  HP_Quadrature1D hp_quad;
693  mpfr_t l, lk, w0, wi, tmp, *weights;
694  mpfr_inits2(hp_quad.default_prec, l, lk, w0, wi, tmp, (mpfr_ptr) 0);
695  weights = new mpfr_t[n];
696  for (int i = 0; i < n; i++)
697  {
698  mpfr_init2(weights[i], hp_quad.default_prec);
699  mpfr_set_si(weights[i], 0, rnd);
700  }
701  hp_quad.SetRelTol(-48); // rtol = 2^(-48) ~ 3.5e-15
702  const int p = n-1;
703  const int m = p/2+1; // number of points for Gauss-Legendre quadrature
704  int hinv = 0, ihoffset = 0; // x_i = (i+ihoffset/2)/hinv
705  switch (type)
706  {
707  case Quadrature1D::ClosedUniform:
708  // x_i = i/p, i=0,...,p
709  hinv = p;
710  ihoffset = 0;
711  break;
712  case Quadrature1D::OpenUniform:
713  // x_i = (i+1)/(p+2), i=0,...,p
714  hinv = p+2;
715  ihoffset = 2;
716  break;
717  case Quadrature1D::OpenHalfUniform:
718  // x_i = (i+1/2)/(p+1), i=0,...,p
719  hinv = p+1;
720  ihoffset = 1;
721  break;
722  default:
723  MFEM_ABORT("invalid Quadrature1D type: " << type);
724  }
725  // set w0 = (-1)^p*(p!)/(hinv^p)
726  mpfr_fac_ui(w0, p, rnd);
727  mpfr_ui_pow_ui(tmp, hinv, p, rnd);
728  mpfr_div(w0, w0, tmp, rnd);
729  if (p%2) { mpfr_neg(w0, w0, rnd); }
730 
731  for (int j = 0; j < m; j++)
732  {
733  hp_quad.ComputeGaussLegendrePoint(m, j);
734 
735  // Compute l = \prod_{i=0}^p (x-x_i) and lk = l/(x-x_k), where
736  // x = hp_quad.GetHPPoint(), x_i = (i+ihoffset/2)/hinv, and x_k is the
737  // node closest to x, i.e. k = min(max(round(x*hinv-ihoffset/2),0),p)
738  mpfr_mul_si(tmp, hp_quad.GetHPPoint(), hinv, rnd);
739  mpfr_sub_d(tmp, tmp, 0.5*ihoffset, rnd);
740  mpfr_round(tmp, tmp);
741  int k = min(max((int)mpfr_get_si(tmp, rnd), 0), p);
742  mpfr_set_si(lk, 1, rnd);
743  for (int i = 0; i <= p; i++)
744  {
745  mpfr_set_si(tmp, 2*i+ihoffset, rnd);
746  mpfr_div_si(tmp, tmp, 2*hinv, rnd);
747  mpfr_sub(tmp, hp_quad.GetHPPoint(), tmp, rnd);
748  if (i != k)
749  {
750  mpfr_mul(lk, lk, tmp, rnd);
751  }
752  else
753  {
754  mpfr_set(l, tmp, rnd);
755  }
756  }
757  mpfr_mul(l, l, lk, rnd);
758  mpfr_set(wi, w0, rnd);
759  for (int i = 0; true; i++)
760  {
761  if (i != k)
762  {
763  // tmp = l/(wi*(x - x_i))
764  mpfr_set_si(tmp, 2*i+ihoffset, rnd);
765  mpfr_div_si(tmp, tmp, 2*hinv, rnd);
766  mpfr_sub(tmp, hp_quad.GetHPPoint(), tmp, rnd);
767  mpfr_mul(tmp, tmp, wi, rnd);
768  mpfr_div(tmp, l, tmp, rnd);
769  }
770  else
771  {
772  // tmp = lk/wi
773  mpfr_div(tmp, lk, wi, rnd);
774  }
775  // weights[i] += hp_quad.weight*tmp
776  mpfr_mul(tmp, tmp, hp_quad.GetHPWeight(), rnd);
777  mpfr_add(weights[i], weights[i], tmp, rnd);
778 
779  if (i == p) { break; }
780 
781  // update wi *= (i+1)/(i-p)
782  mpfr_mul_si(wi, wi, i+1, rnd);
783  mpfr_div_si(wi, wi, i-p, rnd);
784  }
785  }
786  for (int i = 0; i < n; i++)
787  {
788  ir->IntPoint(i).weight = mpfr_get_d(weights[i], rnd);
789  mpfr_clear(weights[i]);
790  }
791  delete [] weights;
792  mpfr_clears(l, lk, w0, wi, tmp, (mpfr_ptr) 0);
793 
794 #endif // MFEM_USE_MPFR
795 
796 }
797 
798 
799 int Quadrature1D::CheckClosed(int type)
800 {
801  switch (type)
802  {
803  case GaussLobatto:
804  case ClosedUniform:
805  return type;
806  default:
807  return Invalid;
808  }
809 }
810 
811 int Quadrature1D::CheckOpen(int type)
812 {
813  switch (type)
814  {
815  case GaussLegendre:
816  case GaussLobatto:
817  case OpenUniform:
818  case ClosedUniform:
819  case OpenHalfUniform:
820  return type; // all types can work as open
821  default:
822  return Invalid;
823  }
824 }
825 
826 
827 IntegrationRules IntRules(0, Quadrature1D::GaussLegendre);
828 
829 IntegrationRules RefinedIntRules(1, Quadrature1D::GaussLegendre);
830 
831 IntegrationRules::IntegrationRules(int Ref, int _type):
832  quad_type(_type)
833 {
834  refined = Ref;
835 
836  if (refined < 0) { own_rules = 0; return; }
837 
838  own_rules = 1;
839 
840  PointIntRules.SetSize(2);
841  PointIntRules = NULL;
842 
843  SegmentIntRules.SetSize(32);
844  SegmentIntRules = NULL;
845 
846  // TriangleIntegrationRule() assumes that this size is >= 26
847  TriangleIntRules.SetSize(32);
848  TriangleIntRules = NULL;
849 
850  SquareIntRules.SetSize(32);
851  SquareIntRules = NULL;
852 
853  // TetrahedronIntegrationRule() assumes that this size is >= 10
854  TetrahedronIntRules.SetSize(32);
855  TetrahedronIntRules = NULL;
856 
857  CubeIntRules.SetSize(32);
858  CubeIntRules = NULL;
859 }
860 
861 const IntegrationRule &IntegrationRules::Get(int GeomType, int Order)
862 {
863  Array<IntegrationRule *> *ir_array;
864 
865  switch (GeomType)
866  {
867  case Geometry::POINT: ir_array = &PointIntRules; Order = 0; break;
868  case Geometry::SEGMENT: ir_array = &SegmentIntRules; break;
869  case Geometry::TRIANGLE: ir_array = &TriangleIntRules; break;
870  case Geometry::SQUARE: ir_array = &SquareIntRules; break;
871  case Geometry::TETRAHEDRON: ir_array = &TetrahedronIntRules; break;
872  case Geometry::CUBE: ir_array = &CubeIntRules; break;
873  default:
874  mfem_error("IntegrationRules::Get(...) : Unknown geometry type!");
875  ir_array = NULL;
876  }
877 
878  if (Order < 0)
879  {
880  Order = 0;
881  }
882 
883  if (!HaveIntRule(*ir_array, Order))
884  {
885 #ifdef MFEM_USE_OPENMP
886  #pragma omp critical
887 #endif
888  {
889  if (!HaveIntRule(*ir_array, Order))
890  {
891  GenerateIntegrationRule(GeomType, Order);
892  }
893  }
894  }
895 
896  return *(*ir_array)[Order];
897 }
898 
899 void IntegrationRules::Set(int GeomType, int Order, IntegrationRule &IntRule)
900 {
901  Array<IntegrationRule *> *ir_array;
902 
903  switch (GeomType)
904  {
905  case Geometry::POINT: ir_array = &PointIntRules; break;
906  case Geometry::SEGMENT: ir_array = &SegmentIntRules; break;
907  case Geometry::TRIANGLE: ir_array = &TriangleIntRules; break;
908  case Geometry::SQUARE: ir_array = &SquareIntRules; break;
909  case Geometry::TETRAHEDRON: ir_array = &TetrahedronIntRules; break;
910  case Geometry::CUBE: ir_array = &CubeIntRules; break;
911  default:
912  mfem_error("IntegrationRules::Set(...) : Unknown geometry type!");
913  ir_array = NULL;
914  }
915 
916  if (HaveIntRule(*ir_array, Order))
917  {
918  MFEM_ABORT("Overwriting set rules is not supported!");
919  }
920 
921  AllocIntRule(*ir_array, Order);
922 
923  (*ir_array)[Order] = &IntRule;
924 }
925 
926 void IntegrationRules::DeleteIntRuleArray(Array<IntegrationRule *> &ir_array)
927 {
928  int i;
929  IntegrationRule *ir = NULL;
930 
931  // Many of the intrules have multiple contiguous copies in the ir_array
932  // so we have to be careful to not delete them twice.
933  for (i = 0; i < ir_array.Size(); i++)
934  {
935  if (ir_array[i] != NULL && ir_array[i] != ir)
936  {
937  ir = ir_array[i];
938  delete ir;
939  }
940  }
941 }
942 
944 {
945  if (!own_rules) { return; }
946 
947  DeleteIntRuleArray(PointIntRules);
948  DeleteIntRuleArray(SegmentIntRules);
949  DeleteIntRuleArray(TriangleIntRules);
950  DeleteIntRuleArray(SquareIntRules);
951  DeleteIntRuleArray(TetrahedronIntRules);
952  DeleteIntRuleArray(CubeIntRules);
953 }
954 
955 
956 IntegrationRule *IntegrationRules::GenerateIntegrationRule(int GeomType,
957  int Order)
958 {
959  switch (GeomType)
960  {
961  case Geometry::POINT:
962  return PointIntegrationRule(Order);
963  case Geometry::SEGMENT:
964  return SegmentIntegrationRule(Order);
965  case Geometry::TRIANGLE:
966  return TriangleIntegrationRule(Order);
967  case Geometry::SQUARE:
968  return SquareIntegrationRule(Order);
970  return TetrahedronIntegrationRule(Order);
971  case Geometry::CUBE:
972  return CubeIntegrationRule(Order);
973  default:
974  mfem_error("IntegrationRules::Set(...) : Unknown geometry type!");
975  return NULL;
976  }
977 }
978 
979 
980 // Integration rules for a point
981 IntegrationRule *IntegrationRules::PointIntegrationRule(int Order)
982 {
983  if (Order > 1)
984  {
985  mfem_error("Point Integration Rule of Order > 1 not defined");
986  return NULL;
987  }
988 
989  IntegrationRule *ir = new IntegrationRule(1);
990  ir->IntPoint(0).x = .0;
991  ir->IntPoint(0).weight = 1.;
992 
993  PointIntRules[1] = PointIntRules[0] = ir;
994 
995  return ir;
996 }
997 
998 // Integration rules for line segment [0,1]
999 IntegrationRule *IntegrationRules::SegmentIntegrationRule(int Order)
1000 {
1001  int RealOrder = GetSegmentRealOrder(Order); // RealOrder >= Order
1002  // Order is one of {RealOrder-1,RealOrder}
1003  AllocIntRule(SegmentIntRules, RealOrder);
1004 
1005  IntegrationRule tmp, *ir;
1006  ir = refined ? &tmp : new IntegrationRule;
1007 
1008  int n = 0;
1009  // n is the number of points to achieve the exact integral of a
1010  // degree Order polynomial
1011  switch (quad_type)
1012  {
1014  {
1015  // Gauss-Legendre is exact for 2*n-1
1016  n = Order/2 + 1;
1017  quad_func.GaussLegendre(n, ir);
1018  break;
1019  }
1021  {
1022  // Gauss-Lobatto is exact for 2*n-3
1023  n = Order/2 + 2;
1024  quad_func.GaussLobatto(n, ir);
1025  break;
1026  }
1028  {
1029  // Open Newton Cotes is exact for n-(n+1)%2 = n-1+n%2
1030  n = Order | 1; // n is always odd
1031  quad_func.OpenUniform(n, ir);
1032  break;
1033  }
1035  {
1036  // Closed Newton Cotes is exact for n-(n+1)%2 = n-1+n%2
1037  n = Order | 1; // n is always odd
1038  quad_func.ClosedUniform(n, ir);
1039  break;
1040  }
1042  {
1043  // Open half Newton Cotes is exact for n-(n+1)%2 = n-1+n%2
1044  n = Order | 1; // n is always odd
1045  quad_func.OpenHalfUniform(n, ir);
1046  break;
1047  }
1048  default:
1049  {
1050  MFEM_ABORT("unknown Quadrature1D type: " << quad_type);
1051  }
1052  }
1053  if (refined)
1054  {
1055  // Effectively passing memory management to SegmentIntegrationRules
1056  ir = new IntegrationRule(2*n);
1057  for (int j = 0; j < n; j++)
1058  {
1059  ir->IntPoint(j).x = tmp.IntPoint(j).x/2.0;
1060  ir->IntPoint(j).weight = tmp.IntPoint(j).weight/2.0;
1061  ir->IntPoint(j+n).x = 0.5 + tmp.IntPoint(j).x/2.0;
1062  ir->IntPoint(j+n).weight = tmp.IntPoint(j).weight/2.0;
1063  }
1064  }
1065  SegmentIntRules[RealOrder-1] = SegmentIntRules[RealOrder] = ir;
1066  return ir;
1067 }
1068 
1069 // Integration rules for reference triangle {[0,0],[1,0],[0,1]}
1070 IntegrationRule *IntegrationRules::TriangleIntegrationRule(int Order)
1071 {
1072  IntegrationRule *ir = NULL;
1073  // Note: Set TriangleIntRules[*] to ir only *after* ir is fully constructed.
1074  // This is needed in multithreaded environment.
1075 
1076  // assuming that orders <= 25 are pre-allocated
1077  switch (Order)
1078  {
1079  case 0: // 1 point - 0 degree
1080  case 1:
1081  ir = new IntegrationRule(1);
1082  ir->AddTriMidPoint(0, 0.5);
1083  TriangleIntRules[0] = TriangleIntRules[1] = ir;
1084  return ir;
1085 
1086  case 2: // 3 point - 2 degree
1087  ir = new IntegrationRule(3);
1088  ir->AddTriPoints3(0, 1./6., 1./6.);
1089  TriangleIntRules[2] = ir;
1090  // interior points
1091  return ir;
1092 
1093  case 3: // 4 point - 3 degree (has one negative weight)
1094  ir = new IntegrationRule(4);
1095  ir->AddTriMidPoint(0, -0.28125); // -9./32.
1096  ir->AddTriPoints3(1, 0.2, 25./96.);
1097  TriangleIntRules[3] = ir;
1098  return ir;
1099 
1100  case 4: // 6 point - 4 degree
1101  ir = new IntegrationRule(6);
1102  ir->AddTriPoints3(0, 0.091576213509770743460, 0.054975871827660933819);
1103  ir->AddTriPoints3(3, 0.44594849091596488632, 0.11169079483900573285);
1104  TriangleIntRules[4] = ir;
1105  return ir;
1106 
1107  case 5: // 7 point - 5 degree
1108  ir = new IntegrationRule(7);
1109  ir->AddTriMidPoint(0, 0.1125);
1110  ir->AddTriPoints3(1, 0.10128650732345633880, 0.062969590272413576298);
1111  ir->AddTriPoints3(4, 0.47014206410511508977, 0.066197076394253090369);
1112  TriangleIntRules[5] = ir;
1113  return ir;
1114 
1115  case 6: // 12 point - 6 degree
1116  ir = new IntegrationRule(12);
1117  ir->AddTriPoints3(0, 0.063089014491502228340, 0.025422453185103408460);
1118  ir->AddTriPoints3(3, 0.24928674517091042129, 0.058393137863189683013);
1119  ir->AddTriPoints6(6, 0.053145049844816947353, 0.31035245103378440542,
1120  0.041425537809186787597);
1121  TriangleIntRules[6] = ir;
1122  return ir;
1123 
1124  case 7: // 12 point - degree 7
1125  ir = new IntegrationRule(12);
1126  ir->AddTriPoints3R(0, 0.062382265094402118174, 0.067517867073916085443,
1127  0.026517028157436251429);
1128  ir->AddTriPoints3R(3, 0.055225456656926611737, 0.32150249385198182267,
1129  0.043881408714446055037);
1130  // slightly better with explicit 3rd area coordinate
1131  ir->AddTriPoints3R(6, 0.034324302945097146470, 0.66094919618673565761,
1132  0.30472650086816719592, 0.028775042784981585738);
1133  ir->AddTriPoints3R(9, 0.51584233435359177926, 0.27771616697639178257,
1134  0.20644149867001643817, 0.067493187009802774463);
1135  TriangleIntRules[7] = ir;
1136  return ir;
1137 
1138  case 8: // 16 point - 8 degree
1139  ir = new IntegrationRule(16);
1140  ir->AddTriMidPoint(0, 0.0721578038388935841255455552445323);
1141  ir->AddTriPoints3(1, 0.170569307751760206622293501491464,
1142  0.0516086852673591251408957751460645);
1143  ir->AddTriPoints3(4, 0.0505472283170309754584235505965989,
1144  0.0162292488115990401554629641708902);
1145  ir->AddTriPoints3(7, 0.459292588292723156028815514494169,
1146  0.0475458171336423123969480521942921);
1147  ir->AddTriPoints6(10, 0.008394777409957605337213834539296,
1148  0.263112829634638113421785786284643,
1149  0.0136151570872174971324223450369544);
1150  TriangleIntRules[8] = ir;
1151  return ir;
1152 
1153  case 9: // 19 point - 9 degree
1154  ir = new IntegrationRule(19);
1155  ir->AddTriMidPoint(0, 0.0485678981413994169096209912536443);
1156  ir->AddTriPoints3b(1, 0.020634961602524744433,
1157  0.0156673501135695352684274156436046);
1158  ir->AddTriPoints3b(4, 0.12582081701412672546,
1159  0.0389137705023871396583696781497019);
1160  ir->AddTriPoints3(7, 0.188203535619032730240961280467335,
1161  0.0398238694636051265164458871320226);
1162  ir->AddTriPoints3(10, 0.0447295133944527098651065899662763,
1163  0.0127888378293490156308393992794999);
1164  ir->AddTriPoints6(13, 0.0368384120547362836348175987833851,
1165  0.2219629891607656956751025276931919,
1166  0.0216417696886446886446886446886446);
1167  TriangleIntRules[9] = ir;
1168  return ir;
1169 
1170  case 10: // 25 point - 10 degree
1171  ir = new IntegrationRule(25);
1172  ir->AddTriMidPoint(0, 0.0454089951913767900476432975500142);
1173  ir->AddTriPoints3b(1, 0.028844733232685245264984935583748,
1174  0.0183629788782333523585030359456832);
1175  ir->AddTriPoints3(4, 0.109481575485037054795458631340522,
1176  0.0226605297177639673913028223692986);
1177  ir->AddTriPoints6(7, 0.141707219414879954756683250476361,
1178  0.307939838764120950165155022930631,
1179  0.0363789584227100543021575883096803);
1180  ir->AddTriPoints6(13, 0.025003534762686386073988481007746,
1181  0.246672560639902693917276465411176,
1182  0.0141636212655287424183685307910495);
1183  ir->AddTriPoints6(19, 0.0095408154002994575801528096228873,
1184  0.0668032510122002657735402127620247,
1185  4.71083348186641172996373548344341E-03);
1186  TriangleIntRules[10] = ir;
1187  return ir;
1188 
1189  case 11: // 28 point -- 11 degree
1190  ir = new IntegrationRule(28);
1191  ir->AddTriPoints6(0, 0.0,
1192  0.141129718717363295960826061941652,
1193  3.68119189165027713212944752369032E-03);
1194  ir->AddTriMidPoint(6, 0.0439886505811161193990465846607278);
1195  ir->AddTriPoints3(7, 0.0259891409282873952600324854988407,
1196  4.37215577686801152475821439991262E-03);
1197  ir->AddTriPoints3(10, 0.0942875026479224956305697762754049,
1198  0.0190407859969674687575121697178070);
1199  ir->AddTriPoints3b(13, 0.010726449965572372516734795387128,
1200  9.42772402806564602923839129555767E-03);
1201  ir->AddTriPoints3(16, 0.207343382614511333452934024112966,
1202  0.0360798487723697630620149942932315);
1203  ir->AddTriPoints3b(19, 0.122184388599015809877869236727746,
1204  0.0346645693527679499208828254519072);
1205  ir->AddTriPoints6(22, 0.0448416775891304433090523914688007,
1206  0.2772206675282791551488214673424523,
1207  0.0205281577146442833208261574536469);
1208  TriangleIntRules[11] = ir;
1209  return ir;
1210 
1211  case 12: // 33 point - 12 degree
1212  ir = new IntegrationRule(33);
1213  ir->AddTriPoints3b(0, 2.35652204523900E-02, 1.28655332202275E-02);
1214  ir->AddTriPoints3b(3, 1.20551215411079E-01, 2.18462722690190E-02);
1215  ir->AddTriPoints3(6, 2.71210385012116E-01, 3.14291121089425E-02);
1216  ir->AddTriPoints3(9, 1.27576145541586E-01, 1.73980564653545E-02);
1217  ir->AddTriPoints3(12, 2.13173504532100E-02, 3.08313052577950E-03);
1218  ir->AddTriPoints6(15, 1.15343494534698E-01, 2.75713269685514E-01,
1219  2.01857788831905E-02);
1220  ir->AddTriPoints6(21, 2.28383322222570E-02, 2.81325580989940E-01,
1221  1.11783866011515E-02);
1222  ir->AddTriPoints6(27, 2.57340505483300E-02, 1.16251915907597E-01,
1223  8.65811555432950E-03);
1224  TriangleIntRules[12] = ir;
1225  return ir;
1226 
1227  case 13: // 37 point - 13 degree
1228  ir = new IntegrationRule(37);
1229  ir->AddTriPoints3b(0, 0.0,
1230  2.67845189554543044455908674650066E-03);
1231  ir->AddTriMidPoint(3, 0.0293480398063595158995969648597808);
1232  ir->AddTriPoints3(4, 0.0246071886432302181878499494124643,
1233  3.92538414805004016372590903990464E-03);
1234  ir->AddTriPoints3b(7, 0.159382493797610632566158925635800,
1235  0.0253344765879434817105476355306468);
1236  ir->AddTriPoints3(10, 0.227900255506160619646298948153592,
1237  0.0250401630452545330803738542916538);
1238  ir->AddTriPoints3(13, 0.116213058883517905247155321839271,
1239  0.0158235572961491595176634480481793);
1240  ir->AddTriPoints3b(16, 0.046794039901841694097491569577008,
1241  0.0157462815379843978450278590138683);
1242  ir->AddTriPoints6(19, 0.0227978945382486125477207592747430,
1243  0.1254265183163409177176192369310890,
1244  7.90126610763037567956187298486575E-03);
1245  ir->AddTriPoints6(25, 0.0162757709910885409437036075960413,
1246  0.2909269114422506044621801030055257,
1247  7.99081889046420266145965132482933E-03);
1248  ir->AddTriPoints6(31, 0.0897330604516053590796290561145196,
1249  0.2723110556841851025078181617634414,
1250  0.0182757511120486476280967518782978);
1251  TriangleIntRules[13] = ir;
1252  return ir;
1253 
1254  case 14: // 42 point - 14 degree
1255  ir = new IntegrationRule(42);
1256  ir->AddTriPoints3b(0, 2.20721792756430E-02, 1.09417906847145E-02);
1257  ir->AddTriPoints3b(3, 1.64710561319092E-01, 1.63941767720625E-02);
1258  ir->AddTriPoints3(6, 2.73477528308839E-01, 2.58870522536460E-02);
1259  ir->AddTriPoints3(9, 1.77205532412543E-01, 2.10812943684965E-02);
1260  ir->AddTriPoints3(12, 6.17998830908730E-02, 7.21684983488850E-03);
1261  ir->AddTriPoints3(15, 1.93909612487010E-02, 2.46170180120000E-03);
1262  ir->AddTriPoints6(18, 5.71247574036480E-02, 1.72266687821356E-01,
1263  1.23328766062820E-02);
1264  ir->AddTriPoints6(24, 9.29162493569720E-02, 3.36861459796345E-01,
1265  1.92857553935305E-02);
1266  ir->AddTriPoints6(30, 1.46469500556540E-02, 2.98372882136258E-01,
1267  7.21815405676700E-03);
1268  ir->AddTriPoints6(36, 1.26833093287200E-03, 1.18974497696957E-01,
1269  2.50511441925050E-03);
1270  TriangleIntRules[14] = ir;
1271  return ir;
1272 
1273  case 15: // 54 point - 15 degree
1274  ir = new IntegrationRule(54);
1275  ir->AddTriPoints3b(0, 0.0834384072617499333, 0.016330909424402645);
1276  ir->AddTriPoints3b(3, 0.192779070841738867, 0.01370640901568218);
1277  ir->AddTriPoints3(6, 0.293197167913025367, 0.01325501829935165);
1278  ir->AddTriPoints3(9, 0.146467786942772933, 0.014607981068243055);
1279  ir->AddTriPoints3(12, 0.0563628676656034333, 0.005292304033121995);
1280  ir->AddTriPoints3(15, 0.0165751268583703333, 0.0018073215320460175);
1281  ir->AddTriPoints6(18, 0.0099122033092248, 0.239534554154794445,
1282  0.004263874050854718);
1283  ir->AddTriPoints6(24, 0.015803770630228, 0.404878807318339958,
1284  0.006958088258345965);
1285  ir->AddTriPoints6(30, 0.00514360881697066667, 0.0950021131130448885,
1286  0.0021459664703674175);
1287  ir->AddTriPoints6(36, 0.0489223257529888, 0.149753107322273969,
1288  0.008117664640887445);
1289  ir->AddTriPoints6(42, 0.0687687486325192, 0.286919612441334979,
1290  0.012803670460631195);
1291  ir->AddTriPoints6(48, 0.1684044181246992, 0.281835668099084562,
1292  0.016544097765822835);
1293  TriangleIntRules[15] = ir;
1294  return ir;
1295 
1296  case 16: // 61 point - 17 degree (used for 16 as well)
1297  case 17:
1298  ir = new IntegrationRule(61);
1299  ir->AddTriMidPoint(0, 1.67185996454015E-02);
1300  ir->AddTriPoints3b(1, 5.65891888645200E-03, 2.54670772025350E-03);
1301  ir->AddTriPoints3b(4, 3.56473547507510E-02, 7.33543226381900E-03);
1302  ir->AddTriPoints3b(7, 9.95200619584370E-02, 1.21754391768360E-02);
1303  ir->AddTriPoints3b(10, 1.99467521245206E-01, 1.55537754344845E-02);
1304  ir->AddTriPoints3 (13, 2.52141267970953E-01, 1.56285556093100E-02);
1305  ir->AddTriPoints3 (16, 1.62047004658461E-01, 1.24078271698325E-02);
1306  ir->AddTriPoints3 (19, 7.58758822607460E-02, 7.02803653527850E-03);
1307  ir->AddTriPoints3 (22, 1.56547269678220E-02, 1.59733808688950E-03);
1308  ir->AddTriPoints6 (25, 1.01869288269190E-02, 3.34319867363658E-01,
1309  4.05982765949650E-03);
1310  ir->AddTriPoints6 (31, 1.35440871671036E-01, 2.92221537796944E-01,
1311  1.34028711415815E-02);
1312  ir->AddTriPoints6 (37, 5.44239242905830E-02, 3.19574885423190E-01,
1313  9.22999660541100E-03);
1314  ir->AddTriPoints6 (43, 1.28685608336370E-02, 1.90704224192292E-01,
1315  4.23843426716400E-03);
1316  ir->AddTriPoints6 (49, 6.71657824135240E-02, 1.80483211648746E-01,
1317  9.14639838501250E-03);
1318  ir->AddTriPoints6 (55, 1.46631822248280E-02, 8.07113136795640E-02,
1319  3.33281600208250E-03);
1320  TriangleIntRules[16] = TriangleIntRules[17] = ir;
1321  return ir;
1322 
1323  case 18: // 73 point - 19 degree (used for 18 as well)
1324  case 19:
1325  ir = new IntegrationRule(73);
1326  ir->AddTriMidPoint(0, 0.0164531656944595);
1327  ir->AddTriPoints3b(1, 0.020780025853987, 0.005165365945636);
1328  ir->AddTriPoints3b(4, 0.090926214604215, 0.011193623631508);
1329  ir->AddTriPoints3b(7, 0.197166638701138, 0.015133062934734);
1330  ir->AddTriPoints3 (10, 0.255551654403098, 0.015245483901099);
1331  ir->AddTriPoints3 (13, 0.17707794215213, 0.0120796063708205);
1332  ir->AddTriPoints3 (16, 0.110061053227952, 0.0080254017934005);
1333  ir->AddTriPoints3 (19, 0.05552862425184, 0.004042290130892);
1334  ir->AddTriPoints3 (22, 0.012621863777229, 0.0010396810137425);
1335  ir->AddTriPoints6 (25, 0.003611417848412, 0.395754787356943,
1336  0.0019424384524905);
1337  ir->AddTriPoints6 (31, 0.13446675453078, 0.307929983880436,
1338  0.012787080306011);
1339  ir->AddTriPoints6 (37, 0.014446025776115, 0.26456694840652,
1340  0.004440451786669);
1341  ir->AddTriPoints6 (43, 0.046933578838178, 0.358539352205951,
1342  0.0080622733808655);
1343  ir->AddTriPoints6 (49, 0.002861120350567, 0.157807405968595,
1344  0.0012459709087455);
1345  ir->AddTriPoints6 (55, 0.075050596975911, 0.223861424097916,
1346  0.0091214200594755);
1347  ir->AddTriPoints6 (61, 0.03464707481676, 0.142421601113383,
1348  0.0051292818680995);
1349  ir->AddTriPoints6 (67, 0.065494628082938, 0.010161119296278,
1350  0.001899964427651);
1351  TriangleIntRules[18] = TriangleIntRules[19] = ir;
1352  return ir;
1353 
1354  case 20: // 85 point - 20 degree
1355  ir = new IntegrationRule(85);
1356  ir->AddTriMidPoint(0, 0.01380521349884976);
1357  ir->AddTriPoints3b(1, 0.001500649324429, 0.00088951477366337);
1358  ir->AddTriPoints3b(4, 0.0941397519389508667, 0.010056199056980585);
1359  ir->AddTriPoints3b(7, 0.2044721240895264, 0.013408923629665785);
1360  ir->AddTriPoints3(10, 0.264500202532787333, 0.012261566900751005);
1361  ir->AddTriPoints3(13, 0.211018964092076767, 0.008197289205347695);
1362  ir->AddTriPoints3(16, 0.107735607171271333, 0.0073979536993248);
1363  ir->AddTriPoints3(19, 0.0390690878378026667, 0.0022896411388521255);
1364  ir->AddTriPoints3(22, 0.0111743797293296333, 0.0008259132577881085);
1365  ir->AddTriPoints6(25, 0.00534961818733726667, 0.0635496659083522206,
1366  0.001174585454287792);
1367  ir->AddTriPoints6(31, 0.00795481706619893333, 0.157106918940706982,
1368  0.0022329628770908965);
1369  ir->AddTriPoints6(37, 0.0104223982812638, 0.395642114364374018,
1370  0.003049783403953986);
1371  ir->AddTriPoints6(43, 0.0109644147961233333, 0.273167570712910522,
1372  0.0034455406635941015);
1373  ir->AddTriPoints6(49, 0.0385667120854623333, 0.101785382485017108,
1374  0.0039987375362390815);
1375  ir->AddTriPoints6(55, 0.0355805078172182, 0.446658549176413815,
1376  0.003693067142668012);
1377  ir->AddTriPoints6(61, 0.0496708163627641333, 0.199010794149503095,
1378  0.00639966593932413);
1379  ir->AddTriPoints6(67, 0.0585197250843317333, 0.3242611836922827,
1380  0.008629035587848275);
1381  ir->AddTriPoints6(73, 0.121497787004394267, 0.208531363210132855,
1382  0.009336472951467735);
1383  ir->AddTriPoints6(79, 0.140710844943938733, 0.323170566536257485,
1384  0.01140911202919763);
1385  TriangleIntRules[20] = ir;
1386  return ir;
1387 
1388  case 21: // 126 point - 25 degree (used also for degrees from 21 to 24)
1389  case 22:
1390  case 23:
1391  case 24:
1392  case 25:
1393  ir = new IntegrationRule(126);
1394  ir->AddTriPoints3b(0, 0.0279464830731742, 0.0040027909400102085);
1395  ir->AddTriPoints3b(3, 0.131178601327651467, 0.00797353841619525);
1396  ir->AddTriPoints3b(6, 0.220221729512072267, 0.006554570615397765);
1397  ir->AddTriPoints3 (9, 0.298443234019804467, 0.00979150048281781);
1398  ir->AddTriPoints3(12, 0.2340441723373718, 0.008235442720768635);
1399  ir->AddTriPoints3(15, 0.151468334609017567, 0.00427363953704605);
1400  ir->AddTriPoints3(18, 0.112733893545993667, 0.004080942928613246);
1401  ir->AddTriPoints3(21, 0.0777156920915263, 0.0030605732699918895);
1402  ir->AddTriPoints3(24, 0.034893093614297, 0.0014542491324683325);
1403  ir->AddTriPoints3(27, 0.00725818462093236667, 0.00034613762283099815);
1404  ir->AddTriPoints6(30, 0.0012923527044422, 0.227214452153364077,
1405  0.0006241445996386985);
1406  ir->AddTriPoints6(36, 0.0053997012721162, 0.435010554853571706,
1407  0.001702376454401511);
1408  ir->AddTriPoints6(42, 0.006384003033975, 0.320309599272204437,
1409  0.0016798271630320255);
1410  ir->AddTriPoints6(48, 0.00502821150199306667, 0.0917503222800051889,
1411  0.000858078269748377);
1412  ir->AddTriPoints6(54, 0.00682675862178186667, 0.0380108358587243835,
1413  0.000740428158357803);
1414  ir->AddTriPoints6(60, 0.0100161996399295333, 0.157425218485311668,
1415  0.0017556563053643425);
1416  ir->AddTriPoints6(66, 0.02575781317339, 0.239889659778533193,
1417  0.003696775074853242);
1418  ir->AddTriPoints6(72, 0.0302278981199158, 0.361943118126060531,
1419  0.003991543738688279);
1420  ir->AddTriPoints6(78, 0.0305049901071620667, 0.0835519609548285602,
1421  0.0021779813065790205);
1422  ir->AddTriPoints6(84, 0.0459565473625693333, 0.148443220732418205,
1423  0.003682528350708916);
1424  ir->AddTriPoints6(90, 0.0674428005402775333, 0.283739708727534955,
1425  0.005481786423209775);
1426  ir->AddTriPoints6(96, 0.0700450914159106, 0.406899375118787573,
1427  0.00587498087177056);
1428  ir->AddTriPoints6(102, 0.0839115246401166, 0.194113987024892542,
1429  0.005007800356899285);
1430  ir->AddTriPoints6(108, 0.120375535677152667, 0.32413434700070316,
1431  0.00665482039381434);
1432  ir->AddTriPoints6(114, 0.148066899157366667, 0.229277483555980969,
1433  0.00707722325261307);
1434  ir->AddTriPoints6(120, 0.191771865867325067, 0.325618122595983752,
1435  0.007440689780584005);
1436  TriangleIntRules[21] =
1437  TriangleIntRules[22] =
1438  TriangleIntRules[23] =
1439  TriangleIntRules[24] =
1440  TriangleIntRules[25] = ir;
1441  return ir;
1442 
1443  default:
1444  // Grundmann-Moller rules
1445  int i = (Order / 2) * 2 + 1; // Get closest odd # >= Order
1446  AllocIntRule(TriangleIntRules, i);
1447  ir = new IntegrationRule;
1448  ir->GrundmannMollerSimplexRule(i/2,2);
1449  TriangleIntRules[i-1] = TriangleIntRules[i] = ir;
1450  return ir;
1451  }
1452 }
1453 
1454 // Integration rules for unit square
1455 IntegrationRule *IntegrationRules::SquareIntegrationRule(int Order)
1456 {
1457  int RealOrder = GetSegmentRealOrder(Order);
1458  // Order is one of {RealOrder-1,RealOrder}
1459  if (!HaveIntRule(SegmentIntRules, RealOrder))
1460  {
1461  SegmentIntegrationRule(RealOrder);
1462  }
1463  AllocIntRule(SquareIntRules, RealOrder); // RealOrder >= Order
1464  SquareIntRules[RealOrder-1] =
1465  SquareIntRules[RealOrder] =
1466  new IntegrationRule(*SegmentIntRules[RealOrder],
1467  *SegmentIntRules[RealOrder]);
1468  return SquareIntRules[Order];
1469 }
1470 
1471 /** Integration rules for reference tetrahedron
1472  {[0,0,0],[1,0,0],[0,1,0],[0,0,1]} */
1473 IntegrationRule *IntegrationRules::TetrahedronIntegrationRule(int Order)
1474 {
1475  IntegrationRule *ir;
1476  // Note: Set TetrahedronIntRules[*] to ir only *after* ir is fully
1477  // constructed. This is needed in multithreaded environment.
1478 
1479  // assuming that orders <= 9 are pre-allocated
1480  switch (Order)
1481  {
1482  case 0: // 1 point - degree 1
1483  case 1:
1484  ir = new IntegrationRule(1);
1485  ir->AddTetMidPoint(0, 1./6.);
1486  TetrahedronIntRules[0] = TetrahedronIntRules[1] = ir;
1487  return ir;
1488 
1489  case 2: // 4 points - degree 2
1490  ir = new IntegrationRule(4);
1491  // ir->AddTetPoints4(0, 0.13819660112501051518, 1./24.);
1492  ir->AddTetPoints4b(0, 0.58541019662496845446, 1./24.);
1493  TetrahedronIntRules[2] = ir;
1494  return ir;
1495 
1496  case 3: // 5 points - degree 3 (negative weight)
1497  ir = new IntegrationRule(5);
1498  ir->AddTetMidPoint(0, -2./15.);
1499  ir->AddTetPoints4b(1, 0.5, 0.075);
1500  TetrahedronIntRules[3] = ir;
1501  return ir;
1502 
1503  case 4: // 11 points - degree 4 (negative weight)
1504  ir = new IntegrationRule(11);
1505  ir->AddTetPoints4(0, 1./14., 343./45000.);
1506  ir->AddTetMidPoint(4, -74./5625.);
1507  ir->AddTetPoints6(5, 0.10059642383320079500, 28./1125.);
1508  TetrahedronIntRules[4] = ir;
1509  return ir;
1510 
1511  case 5: // 14 points - degree 5
1512  ir = new IntegrationRule(14);
1513  ir->AddTetPoints6(0, 0.045503704125649649492,
1514  7.0910034628469110730E-03);
1515  ir->AddTetPoints4(6, 0.092735250310891226402, 0.012248840519393658257);
1516  ir->AddTetPoints4b(10, 0.067342242210098170608,
1517  0.018781320953002641800);
1518  TetrahedronIntRules[5] = ir;
1519  return ir;
1520 
1521  case 6: // 24 points - degree 6
1522  ir = new IntegrationRule(24);
1523  ir->AddTetPoints4(0, 0.21460287125915202929,
1524  6.6537917096945820166E-03);
1525  ir->AddTetPoints4(4, 0.040673958534611353116,
1526  1.6795351758867738247E-03);
1527  ir->AddTetPoints4b(8, 0.032986329573173468968,
1528  9.2261969239424536825E-03);
1529  ir->AddTetPoints12(12, 0.063661001875017525299, 0.26967233145831580803,
1530  8.0357142857142857143E-03);
1531  TetrahedronIntRules[6] = ir;
1532  return ir;
1533 
1534  case 7: // 31 points - degree 7 (negative weight)
1535  ir = new IntegrationRule(31);
1536  ir->AddTetPoints6(0, 0.0, 9.7001763668430335097E-04);
1537  ir->AddTetMidPoint(6, 0.018264223466108820291);
1538  ir->AddTetPoints4(7, 0.078213192330318064374, 0.010599941524413686916);
1539  ir->AddTetPoints4(11, 0.12184321666390517465,
1540  -0.062517740114331851691);
1541  ir->AddTetPoints4b(15, 2.3825066607381275412E-03,
1542  4.8914252630734993858E-03);
1543  ir->AddTetPoints12(19, 0.1, 0.2, 0.027557319223985890653);
1544  TetrahedronIntRules[7] = ir;
1545  return ir;
1546 
1547  case 8: // 43 points - degree 8 (negative weight)
1548  ir = new IntegrationRule(43);
1549  ir->AddTetPoints4(0, 5.7819505051979972532E-03,
1550  1.6983410909288737984E-04);
1551  ir->AddTetPoints4(4, 0.082103588310546723091,
1552  1.9670333131339009876E-03);
1553  ir->AddTetPoints12(8, 0.036607749553197423679, 0.19048604193463345570,
1554  2.1405191411620925965E-03);
1555  ir->AddTetPoints6(20, 0.050532740018894224426,
1556  4.5796838244672818007E-03);
1557  ir->AddTetPoints12(26, 0.22906653611681113960, 0.035639582788534043717,
1558  5.7044858086819185068E-03);
1559  ir->AddTetPoints4(38, 0.20682993161067320408, 0.014250305822866901248);
1560  ir->AddTetMidPoint(42, -0.020500188658639915841);
1561  TetrahedronIntRules[8] = ir;
1562  return ir;
1563 
1564  case 9: // orders 9 and higher -- Grundmann-Moller rules
1565  ir = new IntegrationRule;
1566  ir->GrundmannMollerSimplexRule(4,3);
1567  TetrahedronIntRules[9] = ir;
1568  return ir;
1569 
1570  default: // Grundmann-Moller rules
1571  int i = (Order / 2) * 2 + 1; // Get closest odd # >= Order
1572  AllocIntRule(TetrahedronIntRules, i);
1573  ir = new IntegrationRule;
1574  ir->GrundmannMollerSimplexRule(i/2,3);
1575  TetrahedronIntRules[i-1] = TetrahedronIntRules[i] = ir;
1576  return ir;
1577  }
1578 }
1579 
1580 // Integration rules for reference cube
1581 IntegrationRule *IntegrationRules::CubeIntegrationRule(int Order)
1582 {
1583  int RealOrder = GetSegmentRealOrder(Order);
1584  if (!HaveIntRule(SegmentIntRules, RealOrder))
1585  {
1586  SegmentIntegrationRule(RealOrder);
1587  }
1588  AllocIntRule(CubeIntRules, RealOrder);
1589  CubeIntRules[RealOrder-1] =
1590  CubeIntRules[RealOrder] =
1591  new IntegrationRule(*SegmentIntRules[RealOrder],
1592  *SegmentIntRules[RealOrder],
1593  *SegmentIntRules[RealOrder]);
1594  return CubeIntRules[Order];
1595 }
1596 
1597 }
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:222
int Size() const
Logical size of the array.
Definition: array.hpp:109
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:83
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:861
void ClosedUniform(const int np, IntegrationRule *ir)
Definition: intrules.cpp:565
~IntegrationRules()
Destroys an IntegrationRules object.
Definition: intrules.cpp:943
Container class for integration rules.
Definition: intrules.hpp:279
void OpenUniform(const int np, IntegrationRule *ir)
Definition: intrules.cpp:551
void GaussLobatto(const int np, IntegrationRule *ir)
Definition: intrules.cpp:429
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:225
aka open Newton-Cotes
Definition: intrules.hpp:266
void OpenHalfUniform(const int np, IntegrationRule *ir)
Definition: intrules.cpp:583
void mfem_error(const char *msg)
Definition: error.cpp:106
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:349
Class for integration point with weight.
Definition: intrules.hpp:25
aka closed Newton-Cotes
Definition: intrules.hpp:267
aka &quot;open half&quot; Newton-Cotes
Definition: intrules.hpp:268
Vector data type.
Definition: vector.hpp:36
IntegrationRules RefinedIntRules(1, Quadrature1D::GaussLegendre)
A global object with all refined integration rules.
Definition: intrules.hpp:346
void Set(int GeomType, int Order, IntegrationRule &IntRule)
Definition: intrules.cpp:899
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:343
void Set1w(const double x1, const double w)
Definition: intrules.hpp:77
void GaussLegendre(const int np, IntegrationRule *ir)
Definition: intrules.cpp:350