MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
intrules.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
2// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3// LICENSE and NOTICE for details. LLNL-CODE-806117.
4//
5// This file is part of the MFEM library. For more information and source code
6// availability visit https://mfem.org.
7//
8// MFEM is free software; you can redistribute it and/or modify it under the
9// terms of the BSD-3 license. We welcome feedback and contributions, see file
10// CONTRIBUTING.md for details.
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 "../mesh/nurbs.hpp"
20#include <cmath>
21
22#ifdef MFEM_USE_MPFR
23#include <mpfr.h>
24#endif
25
26using namespace std;
27
28namespace mfem
29{
30
32{
33 int i, j, nx, ny;
34
35 nx = irx.GetNPoints();
36 ny = iry.GetNPoints();
37 SetSize(nx * ny);
39 Order = std::min(irx.GetOrder(), iry.GetOrder());
40
41 for (j = 0; j < ny; j++)
42 {
43 IntegrationPoint &ipy = iry.IntPoint(j);
44 for (i = 0; i < nx; i++)
45 {
46 IntegrationPoint &ipx = irx.IntPoint(i);
47 IntegrationPoint &ip = IntPoint(j*nx+i);
48
49 ip.x = ipx.x;
50 ip.y = ipy.x;
51 ip.weight = ipx.weight * ipy.weight;
52 }
53 }
54}
55
57 IntegrationRule &irz)
58{
59 const int nx = irx.GetNPoints();
60 const int ny = iry.GetNPoints();
61 const int nz = irz.GetNPoints();
62 SetSize(nx*ny*nz);
64 Order = std::min({irx.GetOrder(), iry.GetOrder(), irz.GetOrder()});
65
66 for (int iz = 0; iz < nz; ++iz)
67 {
68 IntegrationPoint &ipz = irz.IntPoint(iz);
69 for (int iy = 0; iy < ny; ++iy)
70 {
71 IntegrationPoint &ipy = iry.IntPoint(iy);
72 for (int ix = 0; ix < nx; ++ix)
73 {
74 IntegrationPoint &ipx = irx.IntPoint(ix);
75 IntegrationPoint &ip = IntPoint(iz*nx*ny + iy*nx + ix);
76
77 ip.x = ipx.x;
78 ip.y = ipy.x;
79 ip.z = ipz.x;
80 ip.weight = ipx.weight*ipy.weight*ipz.weight;
81 }
82 }
83 }
84}
85
87{
88 if (weights.Size() != GetNPoints())
89 {
90 weights.SetSize(GetNPoints());
91 for (int i = 0; i < GetNPoints(); i++)
92 {
93 weights[i] = IntPoint(i).weight;
94 }
95 }
96 return weights;
97}
98
100{
101 for (int i = 0; i < Size(); i++)
102 {
103 IntPoint(i).index = i;
104 }
105}
106
107void IntegrationRule::GrundmannMollerSimplexRule(int s, int n)
108{
109 // for pow on older compilers
110 using std::pow;
111 const int d = 2*s + 1;
112 Vector fact(d + n + 1);
113 Array<int> beta(n), sums(n);
114
115 fact(0) = 1.;
116 for (int i = 1; i < fact.Size(); i++)
117 {
118 fact(i) = fact(i - 1)*i;
119 }
120
121 // number of points is \binom{n + s + 1}{n + 1}
122 int np = 1, f = 1;
123 for (int i = 0; i <= n; i++)
124 {
125 np *= (s + i + 1), f *= (i + 1);
126 }
127 np /= f;
128 SetSize(np);
130 Order = 2*s + 1;
131
132 int pt = 0;
133 for (int i = 0; i <= s; i++)
134 {
135 real_t weight;
136
137 weight = pow(2., -2*s)*pow(static_cast<real_t>(d + n - 2*i),
138 d)/fact(i)/fact(d + n - i);
139 if (i%2)
140 {
141 weight = -weight;
142 }
143
144 // loop over all beta : beta_0 + ... + beta_{n-1} <= s - i
145 int k = s - i;
146 beta = 0;
147 sums = 0;
148 while (true)
149 {
150 IntegrationPoint &ip = IntPoint(pt++);
151 ip.weight = weight;
152 ip.x = real_t(2*beta[0] + 1)/(d + n - 2*i);
153 ip.y = real_t(2*beta[1] + 1)/(d + n - 2*i);
154 if (n == 3)
155 {
156 ip.z = real_t(2*beta[2] + 1)/(d + n - 2*i);
157 }
158
159 int j = 0;
160 while (sums[j] == k)
161 {
162 beta[j++] = 0;
163 if (j == n)
164 {
165 goto done_beta;
166 }
167 }
168 beta[j]++;
169 sums[j]++;
170 for (j--; j >= 0; j--)
171 {
172 sums[j] = sums[j+1];
173 }
174 }
175 done_beta:
176 ;
177 }
178}
179
180IntegrationRule*
182{
183 const int np = this->GetNPoints();
184 const int ne = kv.GetNE();
185
186 IntegrationRule *kvir = new IntegrationRule(ne * np);
187 kvir->SetOrder(GetOrder());
188
189 real_t x0 = kv[0];
190 real_t x1 = x0;
191
192 int id = 0;
193 for (int e=0; e<ne; ++e)
194 {
195 x0 = x1;
196
197 if (e == ne-1)
198 {
199 x1 = kv[kv.Size() - 1];
200 }
201 else
202 {
203 // Find the next unique knot
204 while (id < kv.Size() - 1)
205 {
206 id++;
207 if (kv[id] != x0)
208 {
209 x1 = kv[id];
210 break;
211 }
212 }
213 }
214
215 const real_t s = x1 - x0;
216
217 for (int j=0; j<this->GetNPoints(); ++j)
218 {
219 const real_t x = x0 + (s * (*this)[j].x);
220 (*kvir)[(e * np) + j].Set1w(x, (*this)[j].weight);
221 }
222 }
223
224 return kvir;
225}
226
227#ifdef MFEM_USE_MPFR
228
229// Class for computing hi-precision (HP) quadrature in 1D
230class HP_Quadrature1D
231{
232protected:
233 mpfr_t pi, z, pp, p1, p2, p3, dz, w, rtol;
234
235public:
236 static const mpfr_rnd_t rnd = GMP_RNDN;
237 static const int default_prec = 128;
238
239 // prec = MPFR precision in bits
240 HP_Quadrature1D(const int prec = default_prec)
241 {
242 mpfr_inits2(prec, pi, z, pp, p1, p2, p3, dz, w, rtol, (mpfr_ptr) 0);
243 mpfr_const_pi(pi, rnd);
244 mpfr_set_si_2exp(rtol, 1, -32, rnd); // 2^(-32) < 2.33e-10
245 }
246
247 // set rtol = 2^exponent
248 // this is a tolerance for the last correction of x_i in Newton's algorithm;
249 // this gives roughly rtol^2 accuracy for the final x_i.
250 void SetRelTol(const int exponent = -32)
251 {
252 mpfr_set_si_2exp(rtol, 1, exponent, rnd);
253 }
254
255 // n - number of quadrature points
256 // k - index of the point to compute, 0 <= k < n
257 // see also: QuadratureFunctions1D::GaussLegendre
258 void ComputeGaussLegendrePoint(const int n, const int k)
259 {
260 MFEM_ASSERT(n > 0 && 0 <= k && k < n, "invalid n = " << n
261 << " and/or k = " << k);
262
263 int i = (k < (n+1)/2) ? k+1 : n-k;
264
265 // Initial guess for the x-coordinate:
266 // set z = cos(pi * (i - 0.25) / (n + 0.5)) =
267 // = sin(pi * ((n+1-2*i) / (2*n+1)))
268 mpfr_set_si(z, n+1-2*i, rnd);
269 mpfr_div_si(z, z, 2*n+1, rnd);
270 mpfr_mul(z, z, pi, rnd);
271 mpfr_sin(z, z, rnd);
272
273 bool done = false;
274 while (1)
275 {
276 mpfr_set_si(p2, 1, rnd);
277 mpfr_set(p1, z, rnd);
278 for (int j = 2; j <= n; j++)
279 {
280 mpfr_set(p3, p2, rnd);
281 mpfr_set(p2, p1, rnd);
282 // p1 = ((2 * j - 1) * z * p2 - (j - 1) * p3) / j;
283 mpfr_mul_si(p1, z, 2*j-1, rnd);
284 mpfr_mul_si(p3, p3, j-1, rnd);
285 mpfr_fms(p1, p1, p2, p3, rnd);
286 mpfr_div_si(p1, p1, j, rnd);
287 }
288 // p1 is Legendre polynomial
289
290 // derivative of the Legendre polynomial:
291 // pp = n * (z*p1-p2) / (z*z - 1);
292 mpfr_fms(pp, z, p1, p2, rnd);
293 mpfr_mul_si(pp, pp, n, rnd);
294 mpfr_sqr(p2, z, rnd);
295 mpfr_sub_si(p2, p2, 1, rnd);
296 mpfr_div(pp, pp, p2, rnd);
297
298 if (done) { break; }
299
300 // set delta_z: dz = p1/pp;
301 mpfr_div(dz, p1, pp, rnd);
302 // compute absolute tolerance: atol = rtol*(1-z)
303 mpfr_t &atol = w;
304 mpfr_si_sub(atol, 1, z, rnd);
305 mpfr_mul(atol, atol, rtol, rnd);
306 if (mpfr_cmpabs(dz, atol) <= 0)
307 {
308 done = true;
309 // continue the computation: get pp at the new point, then exit
310 }
311 // update z = z - dz
312 mpfr_sub(z, z, dz, rnd);
313 }
314
315 // map z to (0,1): z = (1 - z)/2
316 mpfr_si_sub(z, 1, z, rnd);
317 mpfr_div_2si(z, z, 1, rnd);
318
319 // weight: w = 1/(4*z*(1 - z)*pp*pp)
320 mpfr_sqr(w, pp, rnd);
321 mpfr_mul_2si(w, w, 2, rnd);
322 mpfr_mul(w, w, z, rnd);
323 mpfr_si_sub(p1, 1, z, rnd); // p1 = 1-z
324 mpfr_mul(w, w, p1, rnd);
325 mpfr_si_div(w, 1, w, rnd);
326
327 if (k >= (n+1)/2) { mpfr_swap(z, p1); }
328 }
329
330 // n - number of quadrature points
331 // k - index of the point to compute, 0 <= k < n
332 // see also: QuadratureFunctions1D::GaussLobatto
333 void ComputeGaussLobattoPoint(const int n, const int k)
334 {
335 MFEM_ASSERT(n > 1 && 0 <= k && k < n, "invalid n = " << n
336 << " and/or k = " << k);
337
338 int i = (k < (n+1)/2) ? k : n-1-k;
339
340 if (i == 0)
341 {
342 mpfr_set_si(z, 0, rnd);
343 mpfr_set_si(p1, 1, rnd);
344 mpfr_set_si(w, n*(n-1), rnd);
345 mpfr_si_div(w, 1, w, rnd); // weight = 1/(n*(n-1))
346 return;
347 }
348 // initial guess is the corresponding Chebyshev point, z:
349 // z = -cos(pi * i/(n-1)) = sin(pi * (2*i-n+1)/(2*n-2))
350 mpfr_set_si(z, 2*i-n+1, rnd);
351 mpfr_div_si(z, z, 2*(n-1), rnd);
352 mpfr_mul(z, pi, z, rnd);
353 mpfr_sin(z, z, rnd);
354 bool done = false;
355 for (int iter = 0 ; true ; ++iter)
356 {
357 // build Legendre polynomials, up to P_{n}(z)
358 mpfr_set_si(p1, 1, rnd);
359 mpfr_set(p2, z, rnd);
360
361 for (int l = 1 ; l < (n-1) ; ++l)
362 {
363 // P_{l+1}(x) = [ (2*l+1)*x*P_l(x) - l*P_{l-1}(x) ]/(l+1)
364 mpfr_mul_si(p1, p1, l, rnd);
365 mpfr_mul_si(p3, z, 2*l+1, rnd);
366 mpfr_fms(p3, p3, p2, p1, rnd);
367 mpfr_div_si(p3, p3, l+1, rnd);
368
369 mpfr_set(p1, p2, rnd);
370 mpfr_set(p2, p3, rnd);
371 }
372 if (done) { break; }
373 // compute dz = resid/deriv = (z*p2 - p1) / (n*p2);
374 mpfr_fms(dz, z, p2, p1, rnd);
375 mpfr_mul_si(p3, p2, n, rnd);
376 mpfr_div(dz, dz, p3, rnd);
377 // update: z = z - dz
378 mpfr_sub(z, z, dz, rnd);
379 // compute absolute tolerance: atol = rtol*(1 + z)
380 mpfr_t &atol = w;
381 mpfr_add_si(atol, z, 1, rnd);
382 mpfr_mul(atol, atol, rtol, rnd);
383 // check for convergence
384 if (mpfr_cmpabs(dz, atol) <= 0)
385 {
386 done = true;
387 // continue the computation: get p2 at the new point, then exit
388 }
389 // If the iteration does not converge fast, something is wrong.
390 MFEM_VERIFY(iter < 8, "n = " << n << ", i = " << i
391 << ", dz = " << mpfr_get_d(dz, rnd));
392 }
393 // Map to the interval [0,1] and scale the weights
394 mpfr_add_si(z, z, 1, rnd);
395 mpfr_div_2si(z, z, 1, rnd);
396 // set the symmetric point
397 mpfr_si_sub(p1, 1, z, rnd);
398 // w = 1/[ n*(n-1)*[P_{n-1}(z)]^2 ]
399 mpfr_sqr(w, p2, rnd);
400 mpfr_mul_si(w, w, n*(n-1), rnd);
401 mpfr_si_div(w, 1, w, rnd);
402
403 if (k >= (n+1)/2) { mpfr_swap(z, p1); }
404 }
405
406 real_t GetPoint() const { return mpfr_get_d(z, rnd); }
407 real_t GetSymmPoint() const { return mpfr_get_d(p1, rnd); }
408 real_t GetWeight() const { return mpfr_get_d(w, rnd); }
409
410 const mpfr_t &GetHPPoint() const { return z; }
411 const mpfr_t &GetHPSymmPoint() const { return p1; }
412 const mpfr_t &GetHPWeight() const { return w; }
413
414 ~HP_Quadrature1D()
415 {
416 mpfr_clears(pi, z, pp, p1, p2, p3, dz, w, rtol, (mpfr_ptr) 0);
417 mpfr_free_cache();
418 }
419};
420
421#endif // MFEM_USE_MPFR
422
423
425{
426 ir->SetSize(np);
427 ir->SetPointIndices();
428 ir->SetOrder(2*np - 1);
429
430 switch (np)
431 {
432 case 1:
433 ir->IntPoint(0).Set1w(0.5, 1.0);
434 return;
435 case 2:
436 ir->IntPoint(0).Set1w(0.21132486540518711775, 0.5);
437 ir->IntPoint(1).Set1w(0.78867513459481288225, 0.5);
438 return;
439 case 3:
440 ir->IntPoint(0).Set1w(0.11270166537925831148, 5./18.);
441 ir->IntPoint(1).Set1w(0.5, 4./9.);
442 ir->IntPoint(2).Set1w(0.88729833462074168852, 5./18.);
443 return;
444 }
445
446 const int n = np;
447 const int m = (n+1)/2;
448
449#ifndef MFEM_USE_MPFR
450
451 for (int i = 1; i <= m; i++)
452 {
453 real_t z = cos(M_PI * (i - 0.25) / (n + 0.5));
454 real_t pp, p1, dz, xi = 0.;
455 bool done = false;
456 while (1)
457 {
458 real_t p2 = 1;
459 p1 = z;
460 for (int j = 2; j <= n; j++)
461 {
462 real_t p3 = p2;
463 p2 = p1;
464 p1 = ((2 * j - 1) * z * p2 - (j - 1) * p3) / j;
465 }
466 // p1 is Legendre polynomial
467
468 pp = n * (z*p1-p2) / (z*z - 1);
469 if (done) { break; }
470
471 dz = p1/pp;
472#ifdef MFEM_USE_SINGLE
473 if (std::abs(dz) < 1e-7)
474#elif defined MFEM_USE_DOUBLE
475 if (std::abs(dz) < 1e-16)
476#else
477 MFEM_ABORT("Floating point type undefined");
478 if (std::abs(dz) < 1e-16)
479#endif
480 {
481 done = true;
482 // map the new point (z-dz) to (0,1):
483 xi = ((1 - z) + dz)/2; // (1 - (z - dz))/2 has bad round-off
484 // continue the computation: get pp at the new point, then exit
485 }
486 // update: z = z - dz
487 z -= dz;
488 }
489
490 ir->IntPoint(i-1).x = xi;
491 ir->IntPoint(n-i).x = 1 - xi;
492 ir->IntPoint(i-1).weight =
493 ir->IntPoint(n-i).weight = 1./(4*xi*(1 - xi)*pp*pp);
494 }
495
496#else // MFEM_USE_MPFR is defined
497
498 HP_Quadrature1D hp_quad;
499 for (int i = 1; i <= m; i++)
500 {
501 hp_quad.ComputeGaussLegendrePoint(n, i-1);
502
503 ir->IntPoint(i-1).x = hp_quad.GetPoint();
504 ir->IntPoint(n-i).x = hp_quad.GetSymmPoint();
505 ir->IntPoint(i-1).weight = ir->IntPoint(n-i).weight = hp_quad.GetWeight();
506 }
507
508#endif // MFEM_USE_MPFR
509
510}
511
513{
514 /* An np point Gauss-Lobatto quadrature has (np - 2) free abscissa the other
515 (2) abscissa are the interval endpoints.
516
517 The interior x_i are the zeros of P'_{np-1}(x). The weights of the
518 interior points on the interval [-1,1] are:
519
520 w_i = 2/(np*(np-1)*[P_{np-1}(x_i)]^2)
521
522 The end point weights (on [-1,1]) are: w_{end} = 2/(np*(np-1)).
523
524 The interior abscissa are found via a nonlinear solve, the initial guess
525 for each point is the corresponding Chebyshev point.
526
527 After we find all points on the interval [-1,1], we will map and scale the
528 points and weights to the MFEM natural interval [0,1].
529
530 References:
531 [1] E. E. Lewis and W. F. Millier, "Computational Methods of Neutron
532 Transport", Appendix A
533 [2] the QUADRULE software by John Burkardt,
534 https://people.sc.fsu.edu/~jburkardt/cpp_src/quadrule/quadrule.cpp
535 */
536
537 ir->SetSize(np);
538 ir->SetPointIndices();
539 if ( np == 1 )
540 {
541 ir->IntPoint(0).Set1w(0.5, 1.0);
542 ir->SetOrder(1);
543 }
544 else
545 {
546 ir->SetOrder(2*np - 3);
547
548#ifndef MFEM_USE_MPFR
549
550 // endpoints and respective weights
551 ir->IntPoint(0).x = 0.0;
552 ir->IntPoint(np-1).x = 1.0;
553 ir->IntPoint(0).weight = ir->IntPoint(np-1).weight = 1.0/(np*(np-1));
554
555 // interior points and weights
556 // use symmetry and compute just half of the points
557 for (int i = 1 ; i <= (np-1)/2 ; ++i)
558 {
559 // initial guess is the corresponding Chebyshev point, x_i:
560 // x_i = -cos(\pi * (i / (np-1)))
561 real_t x_i = std::sin(M_PI * ((real_t)(i)/(np-1) - 0.5));
562 real_t z_i = 0., p_l;
563 bool done = false;
564 for (int iter = 0 ; true ; ++iter)
565 {
566 // build Legendre polynomials, up to P_{np}(x_i)
567 real_t p_lm1 = 1.0;
568 p_l = x_i;
569
570 for (int l = 1 ; l < (np-1) ; ++l)
571 {
572 // The Legendre polynomials can be built by recursion:
573 // x * P_l(x) = 1/(2*l+1)*[ (l+1)*P_{l+1}(x) + l*P_{l-1} ], i.e.
574 // P_{l+1}(x) = [ (2*l+1)*x*P_l(x) - l*P_{l-1} ]/(l+1)
575 real_t p_lp1 = ( (2*l + 1)*x_i*p_l - l*p_lm1)/(l + 1);
576
577 p_lm1 = p_l;
578 p_l = p_lp1;
579 }
580 if (done) { break; }
581 // after this loop, p_l holds P_{np-1}(x_i)
582 // resid = (x^2-1)*P'_{np-1}(x_i)
583 // but use the recurrence relationship
584 // (x^2 -1)P'_l(x) = l*[ x*P_l(x) - P_{l-1}(x) ]
585 // thus, resid = (np-1) * (x_i*p_l - p_lm1)
586
587 // The derivative of the residual is:
588 // \frac{d}{d x} \left[ (x^2 -1)P'_l(x) ] \right] =
589 // l * (l+1) * P_l(x), with l = np-1,
590 // therefore, deriv = np * (np-1) * p_l;
591
592 // compute dx = resid/deriv
593 real_t dx = (x_i*p_l - p_lm1) / (np*p_l);
594#ifdef MFEM_USE_SINGLE
595 if (std::abs(dx) < 1e-7)
596#elif defined MFEM_USE_DOUBLE
597 if (std::abs(dx) < 1e-16)
598#else
599 MFEM_ABORT("Floating point type undefined");
600 if (std::abs(dx) < 1e-16)
601#endif
602 {
603 done = true;
604 // Map the point to the interval [0,1]
605 z_i = ((1.0 + x_i) - dx)/2;
606 // continue the computation: get p_l at the new point, then exit
607 }
608 // If the iteration does not converge fast, something is wrong.
609 MFEM_VERIFY(iter < 8, "np = " << np << ", i = " << i
610 << ", dx = " << dx);
611 // update x_i:
612 x_i -= dx;
613 }
614 // Map to the interval [0,1] and scale the weights
615 IntegrationPoint &ip = ir->IntPoint(i);
616 ip.x = z_i;
617 // w_i = (2/[ n*(n-1)*[P_{n-1}(x_i)]^2 ]) / 2
618 ip.weight = (real_t)(1.0 / (np*(np-1)*p_l*p_l));
619
620 // set the symmetric point
621 IntegrationPoint &symm_ip = ir->IntPoint(np-1-i);
622 symm_ip.x = 1.0 - z_i;
623 symm_ip.weight = ip.weight;
624 }
625
626#else // MFEM_USE_MPFR is defined
627
628 HP_Quadrature1D hp_quad;
629 // use symmetry and compute just half of the points
630 for (int i = 0 ; i <= (np-1)/2 ; ++i)
631 {
632 hp_quad.ComputeGaussLobattoPoint(np, i);
633 ir->IntPoint(i).x = hp_quad.GetPoint();
634 ir->IntPoint(np-1-i).x = hp_quad.GetSymmPoint();
635 ir->IntPoint(i).weight =
636 ir->IntPoint(np-1-i).weight = hp_quad.GetWeight();
637 }
638
639#endif // MFEM_USE_MPFR
640
641 }
642}
643
645{
646 ir->SetSize(np);
647 ir->SetPointIndices();
648 ir->SetOrder(np - 1 + np%2);
649
650 // The Newton-Cotes quadrature is based on weights that integrate exactly the
651 // interpolatory polynomial through the equally spaced quadrature points.
652 for (int i = 0; i < np ; ++i)
653 {
654 ir->IntPoint(i).x = real_t(i+1) / real_t(np + 1);
655 }
656
657 CalculateUniformWeights(ir, Quadrature1D::OpenUniform);
658}
659
661 IntegrationRule* ir)
662{
663 ir->SetSize(np);
664 ir->SetPointIndices();
665 ir->SetOrder(np - 1 + np%2);
666 if ( np == 1 ) // allow this case as "closed"
667 {
668 ir->IntPoint(0).Set1w(0.5, 1.0);
669 return;
670 }
671
672 for (int i = 0; i < np ; ++i)
673 {
674 ir->IntPoint(i).x = real_t(i) / (np-1);
675 }
676
677 CalculateUniformWeights(ir, Quadrature1D::ClosedUniform);
678}
679
681{
682 ir->SetSize(np);
683 ir->SetPointIndices();
684 ir->SetOrder(np - 1 + np%2);
685
686 // Open half points: the centers of np uniform intervals
687 for (int i = 0; i < np ; ++i)
688 {
689 ir->IntPoint(i).x = real_t(2*i+1) / (2*np);
690 }
691
692 CalculateUniformWeights(ir, Quadrature1D::OpenHalfUniform);
693}
694
696{
697 ir->SetSize(np);
698 ir->SetPointIndices();
699 ir->IntPoint(0).x = 0.0;
700 ir->IntPoint(np-1).x = 1.0;
701 ir->SetOrder(np - 1 + np%2); // Is this the correct order?
702
703 if ( np > 2 )
704 {
705 IntegrationRule gl_ir;
706 GaussLegendre(np-1, &gl_ir);
707
708 for (int i = 1; i < np-1; ++i)
709 {
710 ir->IntPoint(i).x = (gl_ir.IntPoint(i-1).x + gl_ir.IntPoint(i).x)/2;
711 }
712 }
713
714 CalculateUniformWeights(ir, Quadrature1D::ClosedGL);
715}
716
718 const int type)
719{
720 IntegrationRule ir(np);
721
722 switch (type)
723 {
725 {
726 GaussLegendre(np,&ir);
727 break;
728 }
730 {
731 GaussLobatto(np, &ir);
732 break;
733 }
735 {
736 OpenUniform(np,&ir);
737 break;
738 }
740 {
741 ClosedUniform(np,&ir);
742 break;
743 }
745 {
746 OpenHalfUniform(np, &ir);
747 break;
748 }
750 {
751 ClosedGL(np, &ir);
752 break;
753 }
755 {
756 MFEM_ABORT("Asking for an unknown type of 1D Quadrature points, "
757 "type = " << type);
758 }
759 }
760
761 for (int i = 0 ; i < np ; ++i)
762 {
763 pts[i] = ir.IntPoint(i).x;
764 }
765}
766
767void QuadratureFunctions1D::CalculateUniformWeights(IntegrationRule *ir,
768 const int type)
769{
770 /* The Lagrange polynomials are:
771 p_i = \prod_{j \neq i} {\frac{x - x_j }{x_i - x_j}}
772
773 The weight associated with each abscissa is the integral of p_i over
774 [0,1]. To calculate the integral of p_i, we use a Gauss-Legendre
775 quadrature rule. This approach does not suffer from bad round-off/
776 cancellation errors for large number of points.
777 */
778 const int n = ir->Size();
779 switch (n)
780 {
781 case 1:
782 ir->IntPoint(0).weight = 1.;
783 return;
784 case 2:
785 ir->IntPoint(0).weight = .5;
786 ir->IntPoint(1).weight = .5;
787 return;
788 }
789
790#ifndef MFEM_USE_MPFR
791
792 // This algorithm should work for any set of points, not just uniform
793 const IntegrationRule &glob_ir = IntRules.Get(Geometry::SEGMENT, n-1);
794 const int m = glob_ir.GetNPoints();
795 Vector xv(n);
796 for (int j = 0; j < n; j++)
797 {
798 xv(j) = ir->IntPoint(j).x;
799 }
800 Poly_1D::Basis basis(n-1, xv.GetData()); // nodal basis, with nodes at 'xv'
801 Vector w(n);
802 // Integrate all nodal basis functions using 'glob_ir':
803 w = 0.0;
804 for (int i = 0; i < m; i++)
805 {
806 const IntegrationPoint &ip = glob_ir.IntPoint(i);
807 basis.Eval(ip.x, xv);
808 w.Add(ip.weight, xv); // w += ip.weight * xv
809 }
810 for (int j = 0; j < n; j++)
811 {
812 ir->IntPoint(j).weight = w(j);
813 }
814
815#else // MFEM_USE_MPFR is defined
816
817 static const mpfr_rnd_t rnd = HP_Quadrature1D::rnd;
818 HP_Quadrature1D hp_quad;
819 mpfr_t l, lk, w0, wi, tmp, *weights;
820 mpfr_inits2(hp_quad.default_prec, l, lk, w0, wi, tmp, (mpfr_ptr) 0);
821 weights = new mpfr_t[n];
822 for (int i = 0; i < n; i++)
823 {
824 mpfr_init2(weights[i], hp_quad.default_prec);
825 mpfr_set_si(weights[i], 0, rnd);
826 }
827 hp_quad.SetRelTol(-48); // rtol = 2^(-48) ~ 3.5e-15
828 const int p = n-1;
829 const int m = p/2+1; // number of points for Gauss-Legendre quadrature
830 int hinv = 0, ihoffset = 0; // x_i = (i+ihoffset/2)/hinv
831 switch (type)
832 {
834 // x_i = i/p, i=0,...,p
835 hinv = p;
836 ihoffset = 0;
837 break;
839 // x_i = (i+1)/(p+2), i=0,...,p
840 hinv = p+2;
841 ihoffset = 2;
842 break;
844 // x_i = (i+1/2)/(p+1), i=0,...,p
845 hinv = p+1;
846 ihoffset = 1;
847 break;
852 MFEM_ABORT("invalid Quadrature1D type: " << type);
853 }
854 // set w0 = (-1)^p*(p!)/(hinv^p)
855 mpfr_fac_ui(w0, p, rnd);
856 mpfr_ui_pow_ui(tmp, hinv, p, rnd);
857 mpfr_div(w0, w0, tmp, rnd);
858 if (p%2) { mpfr_neg(w0, w0, rnd); }
859
860 for (int j = 0; j < m; j++)
861 {
862 hp_quad.ComputeGaussLegendrePoint(m, j);
863
864 // Compute l = \prod_{i=0}^p (x-x_i) and lk = l/(x-x_k), where
865 // x = hp_quad.GetHPPoint(), x_i = (i+ihoffset/2)/hinv, and x_k is the
866 // node closest to x, i.e. k = min(max(round(x*hinv-ihoffset/2),0),p)
867 mpfr_mul_si(tmp, hp_quad.GetHPPoint(), hinv, rnd);
868 mpfr_sub_d(tmp, tmp, 0.5*ihoffset, rnd);
869 mpfr_round(tmp, tmp);
870 int k = min(max((int)mpfr_get_si(tmp, rnd), 0), p);
871 mpfr_set_si(lk, 1, rnd);
872 for (int i = 0; i <= p; i++)
873 {
874 mpfr_set_si(tmp, 2*i+ihoffset, rnd);
875 mpfr_div_si(tmp, tmp, 2*hinv, rnd);
876 mpfr_sub(tmp, hp_quad.GetHPPoint(), tmp, rnd);
877 if (i != k)
878 {
879 mpfr_mul(lk, lk, tmp, rnd);
880 }
881 else
882 {
883 mpfr_set(l, tmp, rnd);
884 }
885 }
886 mpfr_mul(l, l, lk, rnd);
887 mpfr_set(wi, w0, rnd);
888 for (int i = 0; true; i++)
889 {
890 if (i != k)
891 {
892 // tmp = l/(wi*(x - x_i))
893 mpfr_set_si(tmp, 2*i+ihoffset, rnd);
894 mpfr_div_si(tmp, tmp, 2*hinv, rnd);
895 mpfr_sub(tmp, hp_quad.GetHPPoint(), tmp, rnd);
896 mpfr_mul(tmp, tmp, wi, rnd);
897 mpfr_div(tmp, l, tmp, rnd);
898 }
899 else
900 {
901 // tmp = lk/wi
902 mpfr_div(tmp, lk, wi, rnd);
903 }
904 // weights[i] += hp_quad.weight*tmp
905 mpfr_mul(tmp, tmp, hp_quad.GetHPWeight(), rnd);
906 mpfr_add(weights[i], weights[i], tmp, rnd);
907
908 if (i == p) { break; }
909
910 // update wi *= (i+1)/(i-p)
911 mpfr_mul_si(wi, wi, i+1, rnd);
912 mpfr_div_si(wi, wi, i-p, rnd);
913 }
914 }
915 for (int i = 0; i < n; i++)
916 {
917 ir->IntPoint(i).weight = mpfr_get_d(weights[i], rnd);
918 mpfr_clear(weights[i]);
919 }
920 delete [] weights;
921 mpfr_clears(l, lk, w0, wi, tmp, (mpfr_ptr) 0);
922
923#endif // MFEM_USE_MPFR
924
925}
926
927
929{
930 switch (type)
931 {
932 case GaussLobatto:
933 case ClosedUniform:
934 case ClosedGL:
935 return type;
936 default:
937 return Invalid;
938 }
939}
940
942{
943 switch (type)
944 {
945 case GaussLegendre:
946 case GaussLobatto:
947 case OpenUniform:
948 case ClosedUniform:
949 case OpenHalfUniform:
950 return type; // all types can work as open
951 default:
952 return Invalid;
953 }
954}
955
956
958
960
962 : quad_type(type)
963{
964 refined = ref;
965
966 if (refined < 0) { own_rules = 0; return; }
967
968 own_rules = 1;
969
970 const MemoryType h_mt = MemoryType::HOST;
971 PointIntRules.SetSize(2, h_mt);
972 PointIntRules = NULL;
973
974 SegmentIntRules.SetSize(32, h_mt);
975 SegmentIntRules = NULL;
976
977 // TriangleIntegrationRule() assumes that this size is >= 26
978 TriangleIntRules.SetSize(32, h_mt);
979 TriangleIntRules = NULL;
980
981 SquareIntRules.SetSize(32, h_mt);
982 SquareIntRules = NULL;
983
984 // TetrahedronIntegrationRule() assumes that this size is >= 10
985 TetrahedronIntRules.SetSize(32, h_mt);
986 TetrahedronIntRules = NULL;
987
988 PyramidIntRules.SetSize(32, h_mt);
989 PyramidIntRules = NULL;
990
991 PrismIntRules.SetSize(32, h_mt);
992 PrismIntRules = NULL;
993
994 CubeIntRules.SetSize(32, h_mt);
995 CubeIntRules = NULL;
996
997#if defined(MFEM_THREAD_SAFE) && defined(MFEM_USE_OPENMP)
998 IntRuleLocks.SetSize(Geometry::NUM_GEOMETRIES, h_mt);
999 for (int i = 0; i < Geometry::NUM_GEOMETRIES; i++)
1000 {
1001 omp_init_lock(&IntRuleLocks[i]);
1002 }
1003#endif
1004}
1005
1006const IntegrationRule &IntegrationRules::Get(int GeomType, int Order)
1007{
1008 Array<IntegrationRule *> *ir_array = NULL;
1009
1010 switch (GeomType)
1011 {
1012 case Geometry::POINT: ir_array = &PointIntRules; Order = 0; break;
1013 case Geometry::SEGMENT: ir_array = &SegmentIntRules; break;
1014 case Geometry::TRIANGLE: ir_array = &TriangleIntRules; break;
1015 case Geometry::SQUARE: ir_array = &SquareIntRules; break;
1016 case Geometry::TETRAHEDRON: ir_array = &TetrahedronIntRules; break;
1017 case Geometry::CUBE: ir_array = &CubeIntRules; break;
1018 case Geometry::PRISM: ir_array = &PrismIntRules; break;
1019 case Geometry::PYRAMID: ir_array = &PyramidIntRules; break;
1020 case Geometry::INVALID:
1022 MFEM_ABORT("Unknown type of reference element!");
1023 }
1024
1025 if (Order < 0)
1026 {
1027 Order = 0;
1028 }
1029
1030#if defined(MFEM_THREAD_SAFE) && defined(MFEM_USE_OPENMP)
1031 omp_set_lock(&IntRuleLocks[GeomType]);
1032#endif
1033
1034 if (!HaveIntRule(*ir_array, Order))
1035 {
1036 IntegrationRule *ir = GenerateIntegrationRule(GeomType, Order);
1037#ifdef MFEM_DEBUG
1038 int RealOrder = Order;
1039 while (RealOrder+1 < ir_array->Size() && (*ir_array)[RealOrder+1] == ir)
1040 {
1041 RealOrder++;
1042 }
1043 MFEM_VERIFY(RealOrder == ir->GetOrder(), "internal error");
1044#else
1045 MFEM_CONTRACT_VAR(ir);
1046#endif
1047 }
1048
1049#if defined(MFEM_THREAD_SAFE) && defined(MFEM_USE_OPENMP)
1050 omp_unset_lock(&IntRuleLocks[GeomType]);
1051#endif
1052
1053 return *(*ir_array)[Order];
1054}
1055
1056void IntegrationRules::Set(int GeomType, int Order, IntegrationRule &IntRule)
1057{
1058 Array<IntegrationRule *> *ir_array = NULL;
1059
1060 switch (GeomType)
1061 {
1062 case Geometry::POINT: ir_array = &PointIntRules; break;
1063 case Geometry::SEGMENT: ir_array = &SegmentIntRules; break;
1064 case Geometry::TRIANGLE: ir_array = &TriangleIntRules; break;
1065 case Geometry::SQUARE: ir_array = &SquareIntRules; break;
1066 case Geometry::TETRAHEDRON: ir_array = &TetrahedronIntRules; break;
1067 case Geometry::CUBE: ir_array = &CubeIntRules; break;
1068 case Geometry::PRISM: ir_array = &PrismIntRules; break;
1069 case Geometry::PYRAMID: ir_array = &PyramidIntRules; break;
1070 case Geometry::INVALID:
1072 MFEM_ABORT("Unknown type of reference element!");
1073 }
1074
1075#if defined(MFEM_THREAD_SAFE) && defined(MFEM_USE_OPENMP)
1076 omp_set_lock(&IntRuleLocks[GeomType]);
1077#endif
1078
1079 if (HaveIntRule(*ir_array, Order))
1080 {
1081 MFEM_ABORT("Overwriting set rules is not supported!");
1082 }
1083
1084 AllocIntRule(*ir_array, Order);
1085
1086 (*ir_array)[Order] = &IntRule;
1087
1088#if defined(MFEM_THREAD_SAFE) && defined(MFEM_USE_OPENMP)
1089 omp_unset_lock(&IntRuleLocks[GeomType]);
1090#endif
1091}
1092
1093void IntegrationRules::DeleteIntRuleArray(
1094 Array<IntegrationRule *> &ir_array) const
1095{
1096 // Many of the intrules have multiple contiguous copies in the ir_array
1097 // so we have to be careful to not delete them twice.
1098 IntegrationRule *ir = NULL;
1099 for (int i = 0; i < ir_array.Size(); i++)
1100 {
1101 if (ir_array[i] != NULL && ir_array[i] != ir)
1102 {
1103 ir = ir_array[i];
1104 delete ir;
1105 }
1106 }
1107}
1108
1110{
1111#if defined(MFEM_THREAD_SAFE) && defined(MFEM_USE_OPENMP)
1112 for (int i = 0; i < Geometry::NUM_GEOMETRIES; i++)
1113 {
1114 omp_destroy_lock(&IntRuleLocks[i]);
1115 }
1116#endif
1117
1118 if (!own_rules) { return; }
1119
1120 DeleteIntRuleArray(PointIntRules);
1121 DeleteIntRuleArray(SegmentIntRules);
1122 DeleteIntRuleArray(TriangleIntRules);
1123 DeleteIntRuleArray(SquareIntRules);
1124 DeleteIntRuleArray(TetrahedronIntRules);
1125 DeleteIntRuleArray(CubeIntRules);
1126 DeleteIntRuleArray(PrismIntRules);
1127 DeleteIntRuleArray(PyramidIntRules);
1128}
1129
1130
1131IntegrationRule *IntegrationRules::GenerateIntegrationRule(int GeomType,
1132 int Order)
1133{
1134 switch (GeomType)
1135 {
1136 case Geometry::POINT:
1137 return PointIntegrationRule(Order);
1138 case Geometry::SEGMENT:
1139 return SegmentIntegrationRule(Order);
1140 case Geometry::TRIANGLE:
1141 return TriangleIntegrationRule(Order);
1142 case Geometry::SQUARE:
1143 return SquareIntegrationRule(Order);
1145 return TetrahedronIntegrationRule(Order);
1146 case Geometry::CUBE:
1147 return CubeIntegrationRule(Order);
1148 case Geometry::PRISM:
1149 return PrismIntegrationRule(Order);
1150 case Geometry::PYRAMID:
1151 return PyramidIntegrationRule(Order);
1152 case Geometry::INVALID:
1154 MFEM_ABORT("Unknown type of reference element!");
1155 }
1156 return NULL;
1157}
1158
1159
1160// Integration rules for a point
1161IntegrationRule *IntegrationRules::PointIntegrationRule(int Order)
1162{
1163 if (Order > 1)
1164 {
1165 MFEM_ABORT("Point Integration Rule of Order > 1 not defined");
1166 return NULL;
1167 }
1168
1169 IntegrationRule *ir = new IntegrationRule(1);
1170 ir->IntPoint(0).x = .0;
1171 ir->IntPoint(0).weight = 1.;
1172 ir->SetOrder(1);
1173
1174 PointIntRules[1] = PointIntRules[0] = ir;
1175
1176 return ir;
1177}
1178
1179// Integration rules for line segment [0,1]
1180IntegrationRule *IntegrationRules::SegmentIntegrationRule(int Order)
1181{
1182 int RealOrder = GetSegmentRealOrder(Order); // RealOrder >= Order
1183 // Order is one of {RealOrder-1,RealOrder}
1184 AllocIntRule(SegmentIntRules, RealOrder);
1185
1186 IntegrationRule *ir = new IntegrationRule;
1187
1188 int n = 0;
1189 // n is the number of points to achieve the exact integral of a
1190 // degree Order polynomial
1191 switch (quad_type)
1192 {
1194 {
1195 // Gauss-Legendre is exact for 2*n-1
1196 n = Order/2 + 1;
1198 break;
1199 }
1201 {
1202 // Gauss-Lobatto is exact for 2*n-3
1203 n = Order/2 + 2;
1205 break;
1206 }
1208 {
1209 // Open Newton Cotes is exact for n-(n+1)%2 = n-1+n%2
1210 n = Order | 1; // n is always odd
1212 break;
1213 }
1215 {
1216 // Closed Newton Cotes is exact for n-(n+1)%2 = n-1+n%2
1217 n = Order | 1; // n is always odd
1219 break;
1220 }
1222 {
1223 // Open half Newton Cotes is exact for n-(n+1)%2 = n-1+n%2
1224 n = Order | 1; // n is always odd
1226 break;
1227 }
1229 {
1230 MFEM_ABORT("unknown Quadrature1D type: " << quad_type);
1231 }
1232 }
1233 if (refined)
1234 {
1235 // Effectively passing memory management to SegmentIntegrationRules
1236 IntegrationRule *refined_ir = new IntegrationRule(2*n);
1237 refined_ir->SetOrder(ir->GetOrder());
1238 for (int j = 0; j < n; j++)
1239 {
1240 refined_ir->IntPoint(j).x = ir->IntPoint(j).x/2.0;
1241 refined_ir->IntPoint(j).weight = ir->IntPoint(j).weight/2.0;
1242 refined_ir->IntPoint(j+n).x = 0.5 + ir->IntPoint(j).x/2.0;
1243 refined_ir->IntPoint(j+n).weight = ir->IntPoint(j).weight/2.0;
1244 }
1245 delete ir;
1246 ir = refined_ir;
1247 }
1248 SegmentIntRules[RealOrder-1] = SegmentIntRules[RealOrder] = ir;
1249 return ir;
1250}
1251
1252// Integration rules for reference triangle {[0,0],[1,0],[0,1]}
1253IntegrationRule *IntegrationRules::TriangleIntegrationRule(int Order)
1254{
1255 IntegrationRule *ir = NULL;
1256 // Note: Set TriangleIntRules[*] to ir only *after* ir is fully constructed.
1257 // This is needed in multithreaded environment.
1258
1259 // assuming that orders <= 25 are pre-allocated
1260 switch (Order)
1261 {
1262 case 0: // 1 point - degree 1
1263 case 1:
1264 ir = new IntegrationRule(1);
1265 ir->AddTriMidPoint(0, 0.5);
1266 ir->SetOrder(1);
1267 TriangleIntRules[0] = TriangleIntRules[1] = ir;
1268 return ir;
1269
1270 case 2: // 3 point - 2 degree
1271 ir = new IntegrationRule(3);
1272 ir->AddTriPoints3(0, 1./6., 1./6.);
1273 ir->SetOrder(2);
1274 TriangleIntRules[2] = ir;
1275 // interior points
1276 return ir;
1277
1278 case 3: // 4 point - 3 degree (has one negative weight)
1279 ir = new IntegrationRule(4);
1280 ir->AddTriMidPoint(0, -0.28125); // -9./32.
1281 ir->AddTriPoints3(1, 0.2, 25./96.);
1282 ir->SetOrder(3);
1283 TriangleIntRules[3] = ir;
1284 return ir;
1285
1286 case 4: // 6 point - 4 degree
1287 ir = new IntegrationRule(6);
1288 ir->AddTriPoints3(0, 0.091576213509770743460, 0.054975871827660933819);
1289 ir->AddTriPoints3(3, 0.44594849091596488632, 0.11169079483900573285);
1290 ir->SetOrder(4);
1291 TriangleIntRules[4] = ir;
1292 return ir;
1293
1294 case 5: // 7 point - 5 degree
1295 ir = new IntegrationRule(7);
1296 ir->AddTriMidPoint(0, 0.1125);
1297 ir->AddTriPoints3(1, 0.10128650732345633880, 0.062969590272413576298);
1298 ir->AddTriPoints3(4, 0.47014206410511508977, 0.066197076394253090369);
1299 ir->SetOrder(5);
1300 TriangleIntRules[5] = ir;
1301 return ir;
1302
1303 case 6: // 12 point - 6 degree
1304 ir = new IntegrationRule(12);
1305 ir->AddTriPoints3(0, 0.063089014491502228340, 0.025422453185103408460);
1306 ir->AddTriPoints3(3, 0.24928674517091042129, 0.058393137863189683013);
1307 ir->AddTriPoints6(6, 0.053145049844816947353, 0.31035245103378440542,
1308 0.041425537809186787597);
1309 ir->SetOrder(6);
1310 TriangleIntRules[6] = ir;
1311 return ir;
1312
1313 case 7: // 12 point - degree 7
1314 ir = new IntegrationRule(12);
1315 ir->AddTriPoints3R(0, 0.062382265094402118174, 0.067517867073916085443,
1316 0.026517028157436251429);
1317 ir->AddTriPoints3R(3, 0.055225456656926611737, 0.32150249385198182267,
1318 0.043881408714446055037);
1319 // slightly better with explicit 3rd area coordinate
1320 ir->AddTriPoints3R(6, 0.034324302945097146470, 0.66094919618673565761,
1321 0.30472650086816719592, 0.028775042784981585738);
1322 ir->AddTriPoints3R(9, 0.51584233435359177926, 0.27771616697639178257,
1323 0.20644149867001643817, 0.067493187009802774463);
1324 ir->SetOrder(7);
1325 TriangleIntRules[7] = ir;
1326 return ir;
1327
1328 case 8: // 16 point - 8 degree
1329 ir = new IntegrationRule(16);
1330 ir->AddTriMidPoint(0, 0.0721578038388935841255455552445323);
1331 ir->AddTriPoints3(1, 0.170569307751760206622293501491464,
1332 0.0516086852673591251408957751460645);
1333 ir->AddTriPoints3(4, 0.0505472283170309754584235505965989,
1334 0.0162292488115990401554629641708902);
1335 ir->AddTriPoints3(7, 0.459292588292723156028815514494169,
1336 0.0475458171336423123969480521942921);
1337 ir->AddTriPoints6(10, 0.008394777409957605337213834539296,
1338 0.263112829634638113421785786284643,
1339 0.0136151570872174971324223450369544);
1340 ir->SetOrder(8);
1341 TriangleIntRules[8] = ir;
1342 return ir;
1343
1344 case 9: // 19 point - 9 degree
1345 ir = new IntegrationRule(19);
1346 ir->AddTriMidPoint(0, 0.0485678981413994169096209912536443);
1347 ir->AddTriPoints3b(1, 0.020634961602524744433,
1348 0.0156673501135695352684274156436046);
1349 ir->AddTriPoints3b(4, 0.12582081701412672546,
1350 0.0389137705023871396583696781497019);
1351 ir->AddTriPoints3(7, 0.188203535619032730240961280467335,
1352 0.0398238694636051265164458871320226);
1353 ir->AddTriPoints3(10, 0.0447295133944527098651065899662763,
1354 0.0127888378293490156308393992794999);
1355 ir->AddTriPoints6(13, 0.0368384120547362836348175987833851,
1356 0.2219629891607656956751025276931919,
1357 0.0216417696886446886446886446886446);
1358 ir->SetOrder(9);
1359 TriangleIntRules[9] = ir;
1360 return ir;
1361
1362 case 10: // 25 point - 10 degree
1363 ir = new IntegrationRule(25);
1364 ir->AddTriMidPoint(0, 0.0454089951913767900476432975500142);
1365 ir->AddTriPoints3b(1, 0.028844733232685245264984935583748,
1366 0.0183629788782333523585030359456832);
1367 ir->AddTriPoints3(4, 0.109481575485037054795458631340522,
1368 0.0226605297177639673913028223692986);
1369 ir->AddTriPoints6(7, 0.141707219414879954756683250476361,
1370 0.307939838764120950165155022930631,
1371 0.0363789584227100543021575883096803);
1372 ir->AddTriPoints6(13, 0.025003534762686386073988481007746,
1373 0.246672560639902693917276465411176,
1374 0.0141636212655287424183685307910495);
1375 ir->AddTriPoints6(19, 0.0095408154002994575801528096228873,
1376 0.0668032510122002657735402127620247,
1377 4.71083348186641172996373548344341E-03);
1378 ir->SetOrder(10);
1379 TriangleIntRules[10] = ir;
1380 return ir;
1381
1382 case 11: // 28 point -- 11 degree
1383 ir = new IntegrationRule(28);
1384 ir->AddTriPoints6(0, 0.0,
1385 0.141129718717363295960826061941652,
1386 3.68119189165027713212944752369032E-03);
1387 ir->AddTriMidPoint(6, 0.0439886505811161193990465846607278);
1388 ir->AddTriPoints3(7, 0.0259891409282873952600324854988407,
1389 4.37215577686801152475821439991262E-03);
1390 ir->AddTriPoints3(10, 0.0942875026479224956305697762754049,
1391 0.0190407859969674687575121697178070);
1392 ir->AddTriPoints3b(13, 0.010726449965572372516734795387128,
1393 9.42772402806564602923839129555767E-03);
1394 ir->AddTriPoints3(16, 0.207343382614511333452934024112966,
1395 0.0360798487723697630620149942932315);
1396 ir->AddTriPoints3b(19, 0.122184388599015809877869236727746,
1397 0.0346645693527679499208828254519072);
1398 ir->AddTriPoints6(22, 0.0448416775891304433090523914688007,
1399 0.2772206675282791551488214673424523,
1400 0.0205281577146442833208261574536469);
1401 ir->SetOrder(11);
1402 TriangleIntRules[11] = ir;
1403 return ir;
1404
1405 case 12: // 33 point - 12 degree
1406 ir = new IntegrationRule(33);
1407 ir->AddTriPoints3b(0, 2.35652204523900E-02, 1.28655332202275E-02);
1408 ir->AddTriPoints3b(3, 1.20551215411079E-01, 2.18462722690190E-02);
1409 ir->AddTriPoints3(6, 2.71210385012116E-01, 3.14291121089425E-02);
1410 ir->AddTriPoints3(9, 1.27576145541586E-01, 1.73980564653545E-02);
1411 ir->AddTriPoints3(12, 2.13173504532100E-02, 3.08313052577950E-03);
1412 ir->AddTriPoints6(15, 1.15343494534698E-01, 2.75713269685514E-01,
1413 2.01857788831905E-02);
1414 ir->AddTriPoints6(21, 2.28383322222570E-02, 2.81325580989940E-01,
1415 1.11783866011515E-02);
1416 ir->AddTriPoints6(27, 2.57340505483300E-02, 1.16251915907597E-01,
1417 8.65811555432950E-03);
1418 ir->SetOrder(12);
1419 TriangleIntRules[12] = ir;
1420 return ir;
1421
1422 case 13: // 37 point - 13 degree
1423 ir = new IntegrationRule(37);
1424 ir->AddTriPoints3b(0, 0.0,
1425 2.67845189554543044455908674650066E-03);
1426 ir->AddTriMidPoint(3, 0.0293480398063595158995969648597808);
1427 ir->AddTriPoints3(4, 0.0246071886432302181878499494124643,
1428 3.92538414805004016372590903990464E-03);
1429 ir->AddTriPoints3b(7, 0.159382493797610632566158925635800,
1430 0.0253344765879434817105476355306468);
1431 ir->AddTriPoints3(10, 0.227900255506160619646298948153592,
1432 0.0250401630452545330803738542916538);
1433 ir->AddTriPoints3(13, 0.116213058883517905247155321839271,
1434 0.0158235572961491595176634480481793);
1435 ir->AddTriPoints3b(16, 0.046794039901841694097491569577008,
1436 0.0157462815379843978450278590138683);
1437 ir->AddTriPoints6(19, 0.0227978945382486125477207592747430,
1438 0.1254265183163409177176192369310890,
1439 7.90126610763037567956187298486575E-03);
1440 ir->AddTriPoints6(25, 0.0162757709910885409437036075960413,
1441 0.2909269114422506044621801030055257,
1442 7.99081889046420266145965132482933E-03);
1443 ir->AddTriPoints6(31, 0.0897330604516053590796290561145196,
1444 0.2723110556841851025078181617634414,
1445 0.0182757511120486476280967518782978);
1446 ir->SetOrder(13);
1447 TriangleIntRules[13] = ir;
1448 return ir;
1449
1450 case 14: // 42 point - 14 degree
1451 ir = new IntegrationRule(42);
1452 ir->AddTriPoints3b(0, 2.20721792756430E-02, 1.09417906847145E-02);
1453 ir->AddTriPoints3b(3, 1.64710561319092E-01, 1.63941767720625E-02);
1454 ir->AddTriPoints3(6, 2.73477528308839E-01, 2.58870522536460E-02);
1455 ir->AddTriPoints3(9, 1.77205532412543E-01, 2.10812943684965E-02);
1456 ir->AddTriPoints3(12, 6.17998830908730E-02, 7.21684983488850E-03);
1457 ir->AddTriPoints3(15, 1.93909612487010E-02, 2.46170180120000E-03);
1458 ir->AddTriPoints6(18, 5.71247574036480E-02, 1.72266687821356E-01,
1459 1.23328766062820E-02);
1460 ir->AddTriPoints6(24, 9.29162493569720E-02, 3.36861459796345E-01,
1461 1.92857553935305E-02);
1462 ir->AddTriPoints6(30, 1.46469500556540E-02, 2.98372882136258E-01,
1463 7.21815405676700E-03);
1464 ir->AddTriPoints6(36, 1.26833093287200E-03, 1.18974497696957E-01,
1465 2.50511441925050E-03);
1466 ir->SetOrder(14);
1467 TriangleIntRules[14] = ir;
1468 return ir;
1469
1470 case 15: // 54 point - 15 degree
1471 ir = new IntegrationRule(54);
1472 ir->AddTriPoints3b(0, 0.0834384072617499333, 0.016330909424402645);
1473 ir->AddTriPoints3b(3, 0.192779070841738867, 0.01370640901568218);
1474 ir->AddTriPoints3(6, 0.293197167913025367, 0.01325501829935165);
1475 ir->AddTriPoints3(9, 0.146467786942772933, 0.014607981068243055);
1476 ir->AddTriPoints3(12, 0.0563628676656034333, 0.005292304033121995);
1477 ir->AddTriPoints3(15, 0.0165751268583703333, 0.0018073215320460175);
1478 ir->AddTriPoints6(18, 0.0099122033092248, 0.239534554154794445,
1479 0.004263874050854718);
1480 ir->AddTriPoints6(24, 0.015803770630228, 0.404878807318339958,
1481 0.006958088258345965);
1482 ir->AddTriPoints6(30, 0.00514360881697066667, 0.0950021131130448885,
1483 0.0021459664703674175);
1484 ir->AddTriPoints6(36, 0.0489223257529888, 0.149753107322273969,
1485 0.008117664640887445);
1486 ir->AddTriPoints6(42, 0.0687687486325192, 0.286919612441334979,
1487 0.012803670460631195);
1488 ir->AddTriPoints6(48, 0.1684044181246992, 0.281835668099084562,
1489 0.016544097765822835);
1490 ir->SetOrder(15);
1491 TriangleIntRules[15] = ir;
1492 return ir;
1493
1494 case 16: // 61 point - 17 degree (used for 16 as well)
1495 case 17:
1496 ir = new IntegrationRule(61);
1497 ir->AddTriMidPoint(0, 1.67185996454015E-02);
1498 ir->AddTriPoints3b(1, 5.65891888645200E-03, 2.54670772025350E-03);
1499 ir->AddTriPoints3b(4, 3.56473547507510E-02, 7.33543226381900E-03);
1500 ir->AddTriPoints3b(7, 9.95200619584370E-02, 1.21754391768360E-02);
1501 ir->AddTriPoints3b(10, 1.99467521245206E-01, 1.55537754344845E-02);
1502 ir->AddTriPoints3 (13, 2.52141267970953E-01, 1.56285556093100E-02);
1503 ir->AddTriPoints3 (16, 1.62047004658461E-01, 1.24078271698325E-02);
1504 ir->AddTriPoints3 (19, 7.58758822607460E-02, 7.02803653527850E-03);
1505 ir->AddTriPoints3 (22, 1.56547269678220E-02, 1.59733808688950E-03);
1506 ir->AddTriPoints6 (25, 1.01869288269190E-02, 3.34319867363658E-01,
1507 4.05982765949650E-03);
1508 ir->AddTriPoints6 (31, 1.35440871671036E-01, 2.92221537796944E-01,
1509 1.34028711415815E-02);
1510 ir->AddTriPoints6 (37, 5.44239242905830E-02, 3.19574885423190E-01,
1511 9.22999660541100E-03);
1512 ir->AddTriPoints6 (43, 1.28685608336370E-02, 1.90704224192292E-01,
1513 4.23843426716400E-03);
1514 ir->AddTriPoints6 (49, 6.71657824135240E-02, 1.80483211648746E-01,
1515 9.14639838501250E-03);
1516 ir->AddTriPoints6 (55, 1.46631822248280E-02, 8.07113136795640E-02,
1517 3.33281600208250E-03);
1518 ir->SetOrder(17);
1519 TriangleIntRules[16] = TriangleIntRules[17] = ir;
1520 return ir;
1521
1522 case 18: // 73 point - 19 degree (used for 18 as well)
1523 case 19:
1524 ir = new IntegrationRule(73);
1525 ir->AddTriMidPoint(0, 0.0164531656944595);
1526 ir->AddTriPoints3b(1, 0.020780025853987, 0.005165365945636);
1527 ir->AddTriPoints3b(4, 0.090926214604215, 0.011193623631508);
1528 ir->AddTriPoints3b(7, 0.197166638701138, 0.015133062934734);
1529 ir->AddTriPoints3 (10, 0.255551654403098, 0.015245483901099);
1530 ir->AddTriPoints3 (13, 0.17707794215213, 0.0120796063708205);
1531 ir->AddTriPoints3 (16, 0.110061053227952, 0.0080254017934005);
1532 ir->AddTriPoints3 (19, 0.05552862425184, 0.004042290130892);
1533 ir->AddTriPoints3 (22, 0.012621863777229, 0.0010396810137425);
1534 ir->AddTriPoints6 (25, 0.003611417848412, 0.395754787356943,
1535 0.0019424384524905);
1536 ir->AddTriPoints6 (31, 0.13446675453078, 0.307929983880436,
1537 0.012787080306011);
1538 ir->AddTriPoints6 (37, 0.014446025776115, 0.26456694840652,
1539 0.004440451786669);
1540 ir->AddTriPoints6 (43, 0.046933578838178, 0.358539352205951,
1541 0.0080622733808655);
1542 ir->AddTriPoints6 (49, 0.002861120350567, 0.157807405968595,
1543 0.0012459709087455);
1544 ir->AddTriPoints6 (55, 0.075050596975911, 0.223861424097916,
1545 0.0091214200594755);
1546 ir->AddTriPoints6 (61, 0.03464707481676, 0.142421601113383,
1547 0.0051292818680995);
1548 ir->AddTriPoints6 (67, 0.065494628082938, 0.010161119296278,
1549 0.001899964427651);
1550 ir->SetOrder(19);
1551 TriangleIntRules[18] = TriangleIntRules[19] = ir;
1552 return ir;
1553
1554 case 20: // 85 point - 20 degree
1555 ir = new IntegrationRule(85);
1556 ir->AddTriMidPoint(0, 0.01380521349884976);
1557 ir->AddTriPoints3b(1, 0.001500649324429, 0.00088951477366337);
1558 ir->AddTriPoints3b(4, 0.0941397519389508667, 0.010056199056980585);
1559 ir->AddTriPoints3b(7, 0.2044721240895264, 0.013408923629665785);
1560 ir->AddTriPoints3(10, 0.264500202532787333, 0.012261566900751005);
1561 ir->AddTriPoints3(13, 0.211018964092076767, 0.008197289205347695);
1562 ir->AddTriPoints3(16, 0.107735607171271333, 0.0073979536993248);
1563 ir->AddTriPoints3(19, 0.0390690878378026667, 0.0022896411388521255);
1564 ir->AddTriPoints3(22, 0.0111743797293296333, 0.0008259132577881085);
1565 ir->AddTriPoints6(25, 0.00534961818733726667, 0.0635496659083522206,
1566 0.001174585454287792);
1567 ir->AddTriPoints6(31, 0.00795481706619893333, 0.157106918940706982,
1568 0.0022329628770908965);
1569 ir->AddTriPoints6(37, 0.0104223982812638, 0.395642114364374018,
1570 0.003049783403953986);
1571 ir->AddTriPoints6(43, 0.0109644147961233333, 0.273167570712910522,
1572 0.0034455406635941015);
1573 ir->AddTriPoints6(49, 0.0385667120854623333, 0.101785382485017108,
1574 0.0039987375362390815);
1575 ir->AddTriPoints6(55, 0.0355805078172182, 0.446658549176413815,
1576 0.003693067142668012);
1577 ir->AddTriPoints6(61, 0.0496708163627641333, 0.199010794149503095,
1578 0.00639966593932413);
1579 ir->AddTriPoints6(67, 0.0585197250843317333, 0.3242611836922827,
1580 0.008629035587848275);
1581 ir->AddTriPoints6(73, 0.121497787004394267, 0.208531363210132855,
1582 0.009336472951467735);
1583 ir->AddTriPoints6(79, 0.140710844943938733, 0.323170566536257485,
1584 0.01140911202919763);
1585 ir->SetOrder(20);
1586 TriangleIntRules[20] = ir;
1587 return ir;
1588
1589 case 21: // 126 point - 25 degree (used also for degrees from 21 to 24)
1590 case 22:
1591 case 23:
1592 case 24:
1593 case 25:
1594 ir = new IntegrationRule(126);
1595 ir->AddTriPoints3b(0, 0.0279464830731742, 0.0040027909400102085);
1596 ir->AddTriPoints3b(3, 0.131178601327651467, 0.00797353841619525);
1597 ir->AddTriPoints3b(6, 0.220221729512072267, 0.006554570615397765);
1598 ir->AddTriPoints3 (9, 0.298443234019804467, 0.00979150048281781);
1599 ir->AddTriPoints3(12, 0.2340441723373718, 0.008235442720768635);
1600 ir->AddTriPoints3(15, 0.151468334609017567, 0.00427363953704605);
1601 ir->AddTriPoints3(18, 0.112733893545993667, 0.004080942928613246);
1602 ir->AddTriPoints3(21, 0.0777156920915263, 0.0030605732699918895);
1603 ir->AddTriPoints3(24, 0.034893093614297, 0.0014542491324683325);
1604 ir->AddTriPoints3(27, 0.00725818462093236667, 0.00034613762283099815);
1605 ir->AddTriPoints6(30, 0.0012923527044422, 0.227214452153364077,
1606 0.0006241445996386985);
1607 ir->AddTriPoints6(36, 0.0053997012721162, 0.435010554853571706,
1608 0.001702376454401511);
1609 ir->AddTriPoints6(42, 0.006384003033975, 0.320309599272204437,
1610 0.0016798271630320255);
1611 ir->AddTriPoints6(48, 0.00502821150199306667, 0.0917503222800051889,
1612 0.000858078269748377);
1613 ir->AddTriPoints6(54, 0.00682675862178186667, 0.0380108358587243835,
1614 0.000740428158357803);
1615 ir->AddTriPoints6(60, 0.0100161996399295333, 0.157425218485311668,
1616 0.0017556563053643425);
1617 ir->AddTriPoints6(66, 0.02575781317339, 0.239889659778533193,
1618 0.003696775074853242);
1619 ir->AddTriPoints6(72, 0.0302278981199158, 0.361943118126060531,
1620 0.003991543738688279);
1621 ir->AddTriPoints6(78, 0.0305049901071620667, 0.0835519609548285602,
1622 0.0021779813065790205);
1623 ir->AddTriPoints6(84, 0.0459565473625693333, 0.148443220732418205,
1624 0.003682528350708916);
1625 ir->AddTriPoints6(90, 0.0674428005402775333, 0.283739708727534955,
1626 0.005481786423209775);
1627 ir->AddTriPoints6(96, 0.0700450914159106, 0.406899375118787573,
1628 0.00587498087177056);
1629 ir->AddTriPoints6(102, 0.0839115246401166, 0.194113987024892542,
1630 0.005007800356899285);
1631 ir->AddTriPoints6(108, 0.120375535677152667, 0.32413434700070316,
1632 0.00665482039381434);
1633 ir->AddTriPoints6(114, 0.148066899157366667, 0.229277483555980969,
1634 0.00707722325261307);
1635 ir->AddTriPoints6(120, 0.191771865867325067, 0.325618122595983752,
1636 0.007440689780584005);
1637 ir->SetOrder(25);
1638 TriangleIntRules[21] =
1639 TriangleIntRules[22] =
1640 TriangleIntRules[23] =
1641 TriangleIntRules[24] =
1642 TriangleIntRules[25] = ir;
1643 return ir;
1644
1645 default:
1646 // Grundmann-Moller rules
1647 int i = (Order / 2) * 2 + 1; // Get closest odd # >= Order
1648 AllocIntRule(TriangleIntRules, i);
1649 ir = new IntegrationRule;
1650 ir->GrundmannMollerSimplexRule(i/2,2);
1651 TriangleIntRules[i-1] = TriangleIntRules[i] = ir;
1652 return ir;
1653 }
1654}
1655
1656// Integration rules for unit square
1657IntegrationRule *IntegrationRules::SquareIntegrationRule(int Order)
1658{
1659 int RealOrder = GetSegmentRealOrder(Order);
1660 // Order is one of {RealOrder-1,RealOrder}
1661 if (!HaveIntRule(SegmentIntRules, RealOrder))
1662 {
1663 SegmentIntegrationRule(RealOrder);
1664 }
1665 AllocIntRule(SquareIntRules, RealOrder); // RealOrder >= Order
1666 SquareIntRules[RealOrder-1] =
1667 SquareIntRules[RealOrder] =
1668 new IntegrationRule(*SegmentIntRules[RealOrder],
1669 *SegmentIntRules[RealOrder]);
1670 return SquareIntRules[Order];
1671}
1672
1673/** Integration rules for reference tetrahedron
1674 {[0,0,0],[1,0,0],[0,1,0],[0,0,1]} */
1675IntegrationRule *IntegrationRules::TetrahedronIntegrationRule(int Order)
1676{
1677 IntegrationRule *ir;
1678 // Note: Set TetrahedronIntRules[*] to ir only *after* ir is fully
1679 // constructed. This is needed in multithreaded environment.
1680
1681 // assuming that orders <= 9 are pre-allocated
1682 switch (Order)
1683 {
1684 case 0: // 1 point - degree 1
1685 case 1:
1686 ir = new IntegrationRule(1);
1687 ir->AddTetMidPoint(0, 1./6.);
1688 ir->SetOrder(1);
1689 TetrahedronIntRules[0] = TetrahedronIntRules[1] = ir;
1690 return ir;
1691
1692 case 2: // 4 points - degree 2
1693 ir = new IntegrationRule(4);
1694 // ir->AddTetPoints4(0, 0.13819660112501051518, 1./24.);
1695 ir->AddTetPoints4b(0, 0.58541019662496845446, 1./24.);
1696 ir->SetOrder(2);
1697 TetrahedronIntRules[2] = ir;
1698 return ir;
1699
1700 case 3: // 5 points - degree 3 (negative weight)
1701 ir = new IntegrationRule(5);
1702 ir->AddTetMidPoint(0, -2./15.);
1703 ir->AddTetPoints4b(1, 0.5, 0.075);
1704 ir->SetOrder(3);
1705 TetrahedronIntRules[3] = ir;
1706 return ir;
1707
1708 case 4: // 11 points - degree 4 (negative weight)
1709 ir = new IntegrationRule(11);
1710 ir->AddTetPoints4(0, 1./14., 343./45000.);
1711 ir->AddTetMidPoint(4, -74./5625.);
1712 ir->AddTetPoints6(5, 0.10059642383320079500, 28./1125.);
1713 ir->SetOrder(4);
1714 TetrahedronIntRules[4] = ir;
1715 return ir;
1716
1717 case 5: // 14 points - degree 5
1718 ir = new IntegrationRule(14);
1719 ir->AddTetPoints6(0, 0.045503704125649649492,
1720 7.0910034628469110730E-03);
1721 ir->AddTetPoints4(6, 0.092735250310891226402, 0.012248840519393658257);
1722 ir->AddTetPoints4b(10, 0.067342242210098170608,
1723 0.018781320953002641800);
1724 ir->SetOrder(5);
1725 TetrahedronIntRules[5] = ir;
1726 return ir;
1727
1728 case 6: // 24 points - degree 6
1729 ir = new IntegrationRule(24);
1730 ir->AddTetPoints4(0, 0.21460287125915202929,
1731 6.6537917096945820166E-03);
1732 ir->AddTetPoints4(4, 0.040673958534611353116,
1733 1.6795351758867738247E-03);
1734 ir->AddTetPoints4b(8, 0.032986329573173468968,
1735 9.2261969239424536825E-03);
1736 ir->AddTetPoints12(12, 0.063661001875017525299, 0.26967233145831580803,
1737 8.0357142857142857143E-03);
1738 ir->SetOrder(6);
1739 TetrahedronIntRules[6] = ir;
1740 return ir;
1741
1742 case 7: // 31 points - degree 7 (negative weight)
1743 ir = new IntegrationRule(31);
1744 ir->AddTetPoints6(0, 0.0, 9.7001763668430335097E-04);
1745 ir->AddTetMidPoint(6, 0.018264223466108820291);
1746 ir->AddTetPoints4(7, 0.078213192330318064374, 0.010599941524413686916);
1747 ir->AddTetPoints4(11, 0.12184321666390517465,
1748 -0.062517740114331851691);
1749 ir->AddTetPoints4b(15, 2.3825066607381275412E-03,
1750 4.8914252630734993858E-03);
1751 ir->AddTetPoints12(19, 0.1, 0.2, 0.027557319223985890653);
1752 ir->SetOrder(7);
1753 TetrahedronIntRules[7] = ir;
1754 return ir;
1755
1756 case 8: // 43 points - degree 8 (negative weight)
1757 ir = new IntegrationRule(43);
1758 ir->AddTetPoints4(0, 5.7819505051979972532E-03,
1759 1.6983410909288737984E-04);
1760 ir->AddTetPoints4(4, 0.082103588310546723091,
1761 1.9670333131339009876E-03);
1762 ir->AddTetPoints12(8, 0.036607749553197423679, 0.19048604193463345570,
1763 2.1405191411620925965E-03);
1764 ir->AddTetPoints6(20, 0.050532740018894224426,
1765 4.5796838244672818007E-03);
1766 ir->AddTetPoints12(26, 0.22906653611681113960, 0.035639582788534043717,
1767 5.7044858086819185068E-03);
1768 ir->AddTetPoints4(38, 0.20682993161067320408, 0.014250305822866901248);
1769 ir->AddTetMidPoint(42, -0.020500188658639915841);
1770 ir->SetOrder(8);
1771 TetrahedronIntRules[8] = ir;
1772 return ir;
1773
1774 case 9: // orders 9 and higher -- Grundmann-Moller rules
1775 ir = new IntegrationRule;
1776 ir->GrundmannMollerSimplexRule(4,3);
1777 TetrahedronIntRules[9] = ir;
1778 return ir;
1779
1780 default: // Grundmann-Moller rules
1781 int i = (Order / 2) * 2 + 1; // Get closest odd # >= Order
1782 AllocIntRule(TetrahedronIntRules, i);
1783 ir = new IntegrationRule;
1784 ir->GrundmannMollerSimplexRule(i/2,3);
1785 TetrahedronIntRules[i-1] = TetrahedronIntRules[i] = ir;
1786 return ir;
1787 }
1788}
1789
1790// Integration rules for reference pyramid
1791IntegrationRule *IntegrationRules::PyramidIntegrationRule(int Order)
1792{
1793 // This is a simple integration rule adapted from an integration
1794 // rule for a cube which seems to be adequate for now. When we
1795 // implement high order finite elements for pyramids we should
1796 // revisit this and see if we can improve upon it.
1797 const IntegrationRule &irc = Get(Geometry::CUBE, Order);
1798 int npts = irc.GetNPoints();
1799 AllocIntRule(PyramidIntRules, Order);
1800 PyramidIntRules[Order] = new IntegrationRule(npts);
1801 PyramidIntRules[Order]->SetOrder(Order); // FIXME: see comment above
1802
1803 for (int k=0; k<npts; k++)
1804 {
1805 const IntegrationPoint &ipc = irc.IntPoint(k);
1806 IntegrationPoint &ipp = PyramidIntRules[Order]->IntPoint(k);
1807 ipp.x = ipc.x * (1.0 - ipc.z);
1808 ipp.y = ipc.y * (1.0 - ipc.z);
1809 ipp.z = ipc.z;
1810 ipp.weight = ipc.weight / 3.0;
1811 }
1812 return PyramidIntRules[Order];
1813}
1814
1815// Integration rules for reference prism
1816IntegrationRule *IntegrationRules::PrismIntegrationRule(int Order)
1817{
1818 const IntegrationRule &irt = Get(Geometry::TRIANGLE, Order);
1819 const IntegrationRule &irs = Get(Geometry::SEGMENT, Order);
1820 int nt = irt.GetNPoints();
1821 int ns = irs.GetNPoints();
1822 AllocIntRule(PrismIntRules, Order);
1823 PrismIntRules[Order] = new IntegrationRule(nt * ns);
1824 PrismIntRules[Order]->SetOrder(std::min(irt.GetOrder(), irs.GetOrder()));
1825 while (Order < std::min(irt.GetOrder(), irs.GetOrder()))
1826 {
1827 AllocIntRule(PrismIntRules, ++Order);
1828 PrismIntRules[Order] = PrismIntRules[Order-1];
1829 }
1830
1831 for (int ks=0; ks<ns; ks++)
1832 {
1833 const IntegrationPoint &ips = irs.IntPoint(ks);
1834 for (int kt=0; kt<nt; kt++)
1835 {
1836 int kp = ks * nt + kt;
1837 const IntegrationPoint &ipt = irt.IntPoint(kt);
1838 IntegrationPoint &ipp = PrismIntRules[Order]->IntPoint(kp);
1839 ipp.x = ipt.x;
1840 ipp.y = ipt.y;
1841 ipp.z = ips.x;
1842 ipp.weight = ipt.weight * ips.weight;
1843 }
1844 }
1845 return PrismIntRules[Order];
1846}
1847
1848// Integration rules for reference cube
1849IntegrationRule *IntegrationRules::CubeIntegrationRule(int Order)
1850{
1851 int RealOrder = GetSegmentRealOrder(Order);
1852 if (!HaveIntRule(SegmentIntRules, RealOrder))
1853 {
1854 SegmentIntegrationRule(RealOrder);
1855 }
1856 AllocIntRule(CubeIntRules, RealOrder);
1857 CubeIntRules[RealOrder-1] =
1858 CubeIntRules[RealOrder] =
1859 new IntegrationRule(*SegmentIntRules[RealOrder],
1860 *SegmentIntRules[RealOrder],
1861 *SegmentIntRules[RealOrder]);
1862 return CubeIntRules[Order];
1863}
1864
1866 const int patch, const int *ijk,
1867 Array<const KnotVector*> const& kv,
1868 bool & deleteRule) const
1869{
1870 deleteRule = false;
1871
1872 // First check whether a rule has been assigned to element index elem.
1873 auto search = elementToRule.find(elem);
1874 if (search != elementToRule.end())
1875 {
1876 return *elementRule[search->second];
1877 }
1878
1879 MFEM_VERIFY(patchRules1D.NumRows(),
1880 "Undefined rule in NURBSMeshRules::GetElementRule");
1881
1882 // Use a tensor product of rules on the patch.
1883 MFEM_VERIFY(kv.Size() == dim, "");
1884
1885 int np = 1;
1886 std::vector<std::vector<real_t>> el(dim);
1887
1888 std::vector<int> npd;
1889 npd.assign(3, 0);
1890
1891 for (int d=0; d<dim; ++d)
1892 {
1893 const int order = kv[d]->GetOrder();
1894
1895 const real_t kv0 = (*kv[d])[order + ijk[d]];
1896 const real_t kv1 = (*kv[d])[order + ijk[d] + 1];
1897
1898 const bool rightEnd = (order + ijk[d] + 1) == (kv[d]->Size() - 1);
1899
1900 for (int i=0; i<patchRules1D(patch,d)->Size(); ++i)
1901 {
1902 const IntegrationPoint& ip = (*patchRules1D(patch,d))[i];
1903 if (kv0 <= ip.x && (ip.x < kv1 || rightEnd))
1904 {
1905 const real_t x = (ip.x - kv0) / (kv1 - kv0);
1906 el[d].push_back(x);
1907 el[d].push_back(ip.weight);
1908 }
1909 }
1910
1911 npd[d] = el[d].size() / 2;
1912 np *= npd[d];
1913 }
1914
1915 IntegrationRule *irp = new IntegrationRule(np);
1916 deleteRule = true;
1917
1918 // Set (*irp)[i + j*npd[0] + k*npd[0]*npd[1]] =
1919 // (el[0][2*i], el[1][2*j], el[2][2*k])
1920
1921 MFEM_VERIFY(npd[0] > 0 && npd[1] > 0, "Assuming 2D or 3D");
1922
1923 for (int i = 0; i < npd[0]; ++i)
1924 {
1925 for (int j = 0; j < npd[1]; ++j)
1926 {
1927 for (int k = 0; k < std::max(npd[2], 1); ++k)
1928 {
1929 const int id = i + j*npd[0] + k*npd[0]*npd[1];
1930 (*irp)[id].x = el[0][2*i];
1931 (*irp)[id].y = el[1][2*j];
1932
1933 (*irp)[id].weight = el[0][(2*i)+1];
1934 (*irp)[id].weight *= el[1][(2*j)+1];
1935
1936 if (npd[2] > 0)
1937 {
1938 (*irp)[id].z = el[2][2*k];
1939 (*irp)[id].weight *= el[2][(2*k)+1];
1940 }
1941 }
1942 }
1943 }
1944
1945 return *irp;
1946}
1947
1948void NURBSMeshRules::GetIntegrationPointFrom1D(const int patch, int i, int j,
1949 int k, IntegrationPoint & ip)
1950{
1951 MFEM_VERIFY(patchRules1D.NumRows() > 0,
1952 "Assuming patchRules1D is set.");
1953
1954 ip.weight = (*patchRules1D(patch,0))[i].weight;
1955 ip.x = (*patchRules1D(patch,0))[i].x;
1956
1957 if (dim > 1)
1958 {
1959 ip.weight *= (*patchRules1D(patch,1))[j].weight;
1960 ip.y = (*patchRules1D(patch,1))[j].x; // 1D rule only has x
1961 }
1962
1963 if (dim > 2)
1964 {
1965 ip.weight *= (*patchRules1D(patch,2))[k].weight;
1966 ip.z = (*patchRules1D(patch,2))[k].x; // 1D rule only has x
1967 }
1968}
1969
1971{
1972 if ((int) pointToElem.size() == npatches) { return; } // Already set
1973
1974 MFEM_VERIFY(elementToRule.empty() && patchRules1D.NumRows() > 0
1975 && npatches > 0, "Assuming patchRules1D is set.");
1976 MFEM_VERIFY(mesh.NURBSext, "");
1977 MFEM_VERIFY(mesh.Dimension() == dim, "");
1978
1979 pointToElem.resize(npatches);
1980 patchRules1D_KnotSpan.resize(npatches);
1981
1982 // First, find all the elements in each patch.
1983 std::vector<std::vector<int>> patchElements(npatches);
1984
1985 for (int e=0; e<mesh.GetNE(); ++e)
1986 {
1987 patchElements[mesh.NURBSext->GetElementPatch(e)].push_back(e);
1988 }
1989
1990 Array<int> ijk(3);
1991 Array<int> maxijk(3);
1992 Array<int> np(3); // Number of points in each dimension
1993 ijk = 0;
1994
1996
1997 for (int p=0; p<npatches; ++p)
1998 {
1999 patchRules1D_KnotSpan[p].resize(dim);
2000
2001 // For each patch, get the range of ijk.
2002 mesh.NURBSext->GetPatchKnotVectors(p, pkv);
2003 MFEM_VERIFY((int) pkv.Size() == dim, "");
2004
2005 maxijk = 1;
2006 np = 1;
2007 for (int d=0; d<dim; ++d)
2008 {
2009 maxijk[d] = pkv[d]->GetNKS();
2010 np[d] = patchRules1D(p,d)->Size();
2011 }
2012
2013 // For each patch, set a map from ijk to element index.
2014 Array3D<int> ijk2elem(maxijk[0], maxijk[1], maxijk[2]);
2015 ijk2elem = -1;
2016
2017 for (auto elem : patchElements[p])
2018 {
2019 mesh.NURBSext->GetElementIJK(elem, ijk);
2020 MFEM_VERIFY(ijk2elem(ijk[0], ijk[1], ijk[2]) == -1, "");
2021 ijk2elem(ijk[0], ijk[1], ijk[2]) = elem;
2022 }
2023
2024 // For each point, find its ijk and from that its element index.
2025 // It is assumed here that the NURBSFiniteElement kv the same as the
2026 // patch kv.
2027
2028 for (int d=0; d<dim; ++d)
2029 {
2030 patchRules1D_KnotSpan[p][d].SetSize(patchRules1D(p,d)->Size());
2031
2032 for (int r=0; r<patchRules1D(p,d)->Size(); ++r)
2033 {
2034 const IntegrationPoint& ip = (*patchRules1D(p,d))[r];
2035
2036 const int order = pkv[d]->GetOrder();
2037
2038 // Find ijk_d such that ip.x is in the corresponding knot-span.
2039 int ijk_d = 0;
2040 bool found = false;
2041 while (!found)
2042 {
2043 const real_t kv0 = (*pkv[d])[order + ijk_d];
2044 const real_t kv1 = (*pkv[d])[order + ijk_d + 1];
2045
2046 const bool rightEnd = (order + ijk_d + 1) == (pkv[d]->Size() - 1);
2047
2048 if (kv0 <= ip.x && (ip.x < kv1 || rightEnd))
2049 {
2050 found = true;
2051 }
2052 else
2053 {
2054 ijk_d++;
2055 }
2056 }
2057
2058 patchRules1D_KnotSpan[p][d][r] = ijk_d;
2059 }
2060 }
2061
2062 pointToElem[p].SetSize(np[0], np[1], np[2]);
2063 for (int i=0; i<np[0]; ++i)
2064 for (int j=0; j<np[1]; ++j)
2065 for (int k=0; k<np[2]; ++k)
2066 {
2067 const int elem = ijk2elem(patchRules1D_KnotSpan[p][0][i],
2068 patchRules1D_KnotSpan[p][1][j],
2069 patchRules1D_KnotSpan[p][2][k]);
2070 MFEM_VERIFY(elem >= 0, "");
2071 pointToElem[p](i,j,k) = elem;
2072 }
2073 } // Loop (p) over patches
2074}
2075
2077 std::vector<const IntegrationRule*> & ir1D)
2078{
2079 MFEM_VERIFY((int) ir1D.size() == dim, "Wrong dimension");
2080
2081 for (int i=0; i<dim; ++i)
2082 {
2083 patchRules1D(patch,i) = ir1D[i];
2084 }
2085}
2086
2088{
2089 for (int i=0; i<patchRules1D.NumRows(); ++i)
2090 for (int j=0; j<patchRules1D.NumCols(); ++j)
2091 {
2092 delete patchRules1D(i, j);
2093 }
2094}
2095
2096}
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:697
int Size() const
Return the logical size of the array.
Definition: array.hpp:144
@ NUM_GEOMETRIES
Definition: geom.hpp:39
Class for integration point with weight.
Definition: intrules.hpp:35
void Set1w(const real_t x1, const real_t w)
Definition: intrules.hpp:93
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:100
int GetOrder() const
Returns the order of the integration rule.
Definition: intrules.hpp:249
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:256
IntegrationRule * ApplyToKnotIntervals(KnotVector const &kv) const
Return an integration rule for KnotVector kv, defined by applying this rule on each knot interval.
Definition: intrules.cpp:181
const Array< real_t > & GetWeights() const
Return the quadrature weights in a contiguous array.
Definition: intrules.cpp:86
void SetOrder(const int order)
Sets the order of the integration rule. This is only for keeping order information,...
Definition: intrules.hpp:253
void SetPointIndices()
Sets the indices of each quadrature point on initialization.
Definition: intrules.cpp:99
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:259
Container class for integration rules.
Definition: intrules.hpp:416
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:1006
IntegrationRules(int ref=0, int type=Quadrature1D::GaussLegendre)
Definition: intrules.cpp:961
void Set(int GeomType, int Order, IntegrationRule &IntRule)
Definition: intrules.cpp:1056
~IntegrationRules()
Destroys an IntegrationRules object.
Definition: intrules.cpp:1109
int Size() const
Definition: nurbs.hpp:54
int GetNE() const
Definition: nurbs.hpp:50
Mesh data type.
Definition: mesh.hpp:56
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:290
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:1226
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1160
void GetElementIJK(int elem, Array< int > &ijk)
Definition: nurbs.cpp:4676
int GetElementPatch(int elem) const
Returns the index of the patch containing element elem.
Definition: nurbs.hpp:592
void GetPatchKnotVectors(int p, Array< KnotVector * > &kv)
Definition: nurbs.cpp:3136
void Finalize(Mesh const &mesh)
Finalize() must be called before this class can be used for assembly. In particular,...
Definition: intrules.cpp:1970
IntegrationRule & GetElementRule(const int elem, const int patch, const int *ijk, Array< const KnotVector * > const &kv, bool &deleteRule) const
Returns a rule for the element.
Definition: intrules.cpp:1865
void SetPatchRules1D(const int patch, std::vector< const IntegrationRule * > &ir1D)
Set 1D integration rules to be used as a tensor product rule on the patch with index patch....
Definition: intrules.cpp:2076
void GetIntegrationPointFrom1D(const int patch, int i, int j, int k, IntegrationPoint &ip)
For tensor product rules defined on each patch by SetPatchRules1D(), return the integration point wit...
Definition: intrules.cpp:1948
@ ClosedUniform
aka closed Newton-Cotes
Definition: intrules.hpp:402
@ ClosedGL
aka closed Gauss Legendre
Definition: intrules.hpp:404
@ OpenHalfUniform
aka "open half" Newton-Cotes
Definition: intrules.hpp:403
@ OpenUniform
aka open Newton-Cotes
Definition: intrules.hpp:401
static int CheckOpen(int type)
If the Quadrature1D type is not open return Invalid; otherwise return type.
Definition: intrules.cpp:941
static int CheckClosed(int type)
If the Quadrature1D type is not closed return Invalid; otherwise return type.
Definition: intrules.cpp:928
static void GaussLegendre(const int np, IntegrationRule *ir)
Definition: intrules.cpp:424
static void ClosedUniform(const int np, IntegrationRule *ir)
Definition: intrules.cpp:660
static void OpenUniform(const int np, IntegrationRule *ir)
Definition: intrules.cpp:644
static void ClosedGL(const int np, IntegrationRule *ir)
Definition: intrules.cpp:695
static void GivePolyPoints(const int np, real_t *pts, const int type)
Definition: intrules.cpp:717
static void OpenHalfUniform(const int np, IntegrationRule *ir)
Definition: intrules.cpp:680
static void GaussLobatto(const int np, IntegrationRule *ir)
Definition: intrules.cpp:512
Vector data type.
Definition: vector.hpp:80
Vector beta
int dim
Definition: ex24.cpp:53
float real_t
Definition: config.hpp:43
MemoryType
Memory types supported by MFEM.
Definition: mem_manager.hpp:39
@ HOST
Host memory; using new[] and delete[].
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition: lor_mms.hpp:30
IntegrationRules RefinedIntRules(1, Quadrature1D::GaussLegendre)
A global object with all refined integration rules.
Definition: intrules.hpp:489
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:486
real_t p(const Vector &x, real_t t)
Definition: navier_mms.cpp:53
RefCoord s[3]
void pts(int iphi, int t, real_t x[])