deal.II version GIT relicensing-2220-ge8e5360a96 2024-12-05 09:50:00+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
quadrature_lib.cc
Go to the documentation of this file.
1// ------------------------------------------------------------------------
2//
3// SPDX-License-Identifier: LGPL-2.1-or-later
4// Copyright (C) 1998 - 2024 by the deal.II authors
5//
6// This file is part of the deal.II library.
7//
8// Part of the source code is dual licensed under Apache-2.0 WITH
9// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
10// governing the source code and code contributions can be found in
11// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
12//
13// ------------------------------------------------------------------------
14
19
22
25#include <deal.II/grid/tria.h>
26
27#include <algorithm>
28#include <cmath>
29#include <functional>
30#include <limits>
31
32
34
35
36// please note: for a given dimension, we need the quadrature formulae
37// for all lower dimensions as well. That is why in this file the check
38// is for deal_II_dimension >= any_number and not for ==
39
40
41
42template <>
43QGauss<0>::QGauss(const unsigned int)
44 : // there are n_q^dim == 1
45 // points
46 Quadrature<0>(1)
47{
48 // the single quadrature point gets unit
49 // weight
50 this->weights[0] = 1;
51}
52
53
54
55template <>
57 : // there are n_q^dim == 1
58 // points
59 Quadrature<0>(1)
60{
61 // the single quadrature point gets unit
62 // weight
63 this->weights[0] = 1;
64}
65
66
67
68template <>
69QGauss<1>::QGauss(const unsigned int n)
70 : Quadrature<1>(n)
71{
72 if (n == 0)
73 return;
74
75 std::vector<long double> points =
76 Polynomials::jacobi_polynomial_roots<long double>(n, 0, 0);
77
78 for (unsigned int i = 0; i < (points.size() + 1) / 2; ++i)
79 {
80 this->quadrature_points[i][0] = points[i];
81 this->quadrature_points[n - i - 1][0] = 1. - points[i];
82
83 // derivative of Jacobi polynomial
84 const long double pp =
85 0.5 * (n + 1) *
86 Polynomials::jacobi_polynomial_value(n - 1, 1, 1, points[i]);
87 const long double x = -1. + 2. * points[i];
88 const double w = 1. / ((1. - x * x) * pp * pp);
89 this->weights[i] = w;
90 this->weights[n - i - 1] = w;
91 }
92}
93
94namespace internal
95{
96 namespace QGaussLobatto
97 {
102 long double
103 gamma(const unsigned int n)
104 {
105 long double result = n - 1;
106 for (int i = n - 2; i > 1; --i)
107 result *= i;
108 return result;
109 }
110
111
112
120 std::vector<long double>
121 compute_quadrature_weights(const std::vector<long double> &x,
122 const int alpha,
123 const int beta)
124 {
125 const unsigned int q = x.size();
126 std::vector<long double> w(q);
127
128 const long double factor =
129 Utilities::pow(2, alpha + beta + 1) * gamma(alpha + q) *
130 gamma(beta + q) / ((q - 1) * gamma(q) * gamma(alpha + beta + q + 1));
131 for (unsigned int i = 0; i < q; ++i)
132 {
133 const long double s =
134 Polynomials::jacobi_polynomial_value(q - 1, alpha, beta, x[i]);
135 w[i] = factor / (s * s);
136 }
137 w[0] *= (beta + 1);
138 w[q - 1] *= (alpha + 1);
139
140 return w;
141 }
142 } // namespace QGaussLobatto
143} // namespace internal
144
145
146namespace internal
147{
148 namespace QGaussRadau
149 {
150
151 // Implements lookup table after affine transformation to [0,1].
152 //
153 // Analytical values for [-1,1] and n < 4 listed on
154 // https://mathworld.wolfram.com/RadauQuadrature.html
155 // Values for n > 3 calculated with the Julia Package
156 // FastGaussQuadrature.jl
157 // https://github.com/JuliaApproximation/FastGaussQuadrature.jl
158 //
159 std::vector<double>
160 get_left_quadrature_points(const unsigned int n)
161 {
162 std::vector<double> q_points(n);
163 switch (n)
164 {
165 case 1:
166 q_points[0] = 0.;
167 break;
168 case 2:
169 q_points[0] = 0.;
170 q_points[1] = 2. / 3.;
171 break;
172 case 3:
173 q_points[0] = 0.;
174 q_points[1] = (6. - std::sqrt(6)) * 0.1;
175 q_points[2] = (6. + std::sqrt(6)) * 0.1;
176 break;
177
178 case 4:
179 q_points[0] = 0.000000000000000000;
180 q_points[1] = 0.212340538239152943;
181 q_points[2] = 0.590533135559265343;
182 q_points[3] = 0.911412040487296071;
183 break;
184 case 5:
185 q_points[0] = 0.000000000000000000;
186 q_points[1] = 0.139759864343780571;
187 q_points[2] = 0.416409567631083166;
188 q_points[3] = 0.723156986361876197;
189 q_points[4] = 0.942895803885482331;
190 break;
191 case 6:
192 q_points[0] = 0.000000000000000000;
193 q_points[1] = 0.098535085798826416;
194 q_points[2] = 0.304535726646363913;
195 q_points[3] = 0.562025189752613841;
196 q_points[4] = 0.801986582126391845;
197 q_points[5] = 0.960190142948531222;
198 break;
199 case 7:
200 q_points[0] = 0.000000000000000000;
201 q_points[1] = 0.073054328680258851;
202 q_points[2] = 0.230766137969945495;
203 q_points[3] = 0.441328481228449865;
204 q_points[4] = 0.663015309718845702;
205 q_points[5] = 0.851921400331515644;
206 q_points[6] = 0.970683572840215114;
207 break;
208 case 8:
209 q_points[0] = 0.000000000000000000;
210 q_points[1] = 0.056262560536922135;
211 q_points[2] = 0.180240691736892389;
212 q_points[3] = 0.352624717113169672;
213 q_points[4] = 0.547153626330555420;
214 q_points[5] = 0.734210177215410598;
215 q_points[6] = 0.885320946839095790;
216 q_points[7] = 0.977520613561287499;
217 break;
218 default:
220 break;
221 }
222 return q_points;
223 }
224
225 std::vector<double>
226 get_quadrature_points(const unsigned int n,
227 const ::QGaussRadau<1>::EndPoint end_point)
228 {
229 std::vector<double> left_points = get_left_quadrature_points(n);
230 switch (end_point)
231 {
232 case ::QGaussRadau<1>::EndPoint::left:
233 return left_points;
234 case ::QGaussRadau<1>::EndPoint::right:
235 {
236 std::vector<double> points(n);
237 for (unsigned int i = 0; i < n; ++i)
238 {
239 points[n - i - 1] = 1. - left_points[i];
240 }
241 return points;
242 }
243 default:
244 Assert(
245 false,
247 "This constructor can only be called with either "
248 "QGaussRadau::left or QGaussRadau::right as second argument."));
249 return {};
250 }
251 }
252
253 // Implements lookup table after affine transformation to [0,1].
254 //
255 // Analytical values for [-1,1] and n < 4 listed on
256 // https://mathworld.wolfram.com/RadauQuadrature.html
257 // Values for n > 3 calculated with the Julia Package
258 // FastGaussQuadrature.jl
259 // https://github.com/JuliaApproximation/FastGaussQuadrature.jl
260 //
261 std::vector<double>
262 get_left_quadrature_weights(const unsigned int n)
263 {
264 std::vector<double> weights(n);
265 switch (n)
266 {
267 case 1:
268 weights[0] = 1.;
269 break;
270 case 2:
271 weights[0] = 0.25;
272 weights[1] = 0.75;
273 break;
274 case 3:
275 weights[0] = 1. / 9.;
276 weights[1] = (16. + std::sqrt(6)) / 36.;
277 weights[2] = (16. - std::sqrt(6)) / 36.;
278 break;
279 case 4:
280 weights[0] = 0.062500000000000000;
281 weights[1] = 0.328844319980059696;
282 weights[2] = 0.388193468843171852;
283 weights[3] = 0.220462211176768369;
284 break;
285 case 5:
286 weights[0] = 0.040000000000000001;
287 weights[1] = 0.223103901083570894;
288 weights[2] = 0.311826522975741427;
289 weights[3] = 0.281356015149462124;
290 weights[4] = 0.143713560791225797;
291 break;
292 case 6:
293 weights[0] = 0.027777777777777776;
294 weights[1] = 0.159820376610255471;
295 weights[2] = 0.242693594234484888;
296 weights[3] = 0.260463391594787597;
297 weights[4] = 0.208450667155953895;
298 weights[5] = 0.100794192626740456;
299 break;
300 case 7:
301 weights[0] = 0.020408163265306121;
302 weights[1] = 0.119613744612656100;
303 weights[2] = 0.190474936822115581;
304 weights[3] = 0.223554914507283209;
305 weights[4] = 0.212351889502977870;
306 weights[5] = 0.159102115733650767;
307 weights[6] = 0.074494235556010341;
308 break;
309 case 8:
310 weights[0] = 0.015625000000000000;
311 weights[1] = 0.092679077401489660;
312 weights[2] = 0.152065310323392683;
313 weights[3] = 0.188258772694559262;
314 weights[4] = 0.195786083726246729;
315 weights[5] = 0.173507397817250691;
316 weights[6] = 0.124823950664932445;
317 weights[7] = 0.057254407372128648;
318 break;
319
320 default:
322 break;
323 }
324 return weights;
325 }
326
327 std::vector<double>
328 get_quadrature_weights(const unsigned int n,
329 const ::QGaussRadau<1>::EndPoint end_point)
330 {
331 std::vector<double> left_weights = get_left_quadrature_weights(n);
332 switch (end_point)
333 {
334 case ::QGaussRadau<1>::EndPoint::left:
335 return left_weights;
336 case ::QGaussRadau<1>::EndPoint::right:
337 {
338 std::vector<double> weights(n);
339 for (unsigned int i = 0; i < n; ++i)
340 {
341 weights[n - i - 1] = left_weights[i];
342 }
343 return weights;
344 }
345 default:
346 Assert(false,
348 "This constructor can only be called with either "
349 "QGaussRadau::EndPoint::left or "
350 "QGaussRadau::EndPoint::right as second argument."));
351 return {};
352 }
353 }
354 } // namespace QGaussRadau
355} // namespace internal
356
357#ifndef DOXYGEN
358template <>
359QGaussRadau<1>::QGaussRadau(const unsigned int n, const EndPoint end_point)
360 : Quadrature<1>(n)
361 , end_point(end_point)
362{
363 Assert(n > 0, ExcMessage("Need at least one point for quadrature rules."));
364 std::vector<double> p =
366 std::vector<double> w =
368
369 for (unsigned int i = 0; i < this->size(); ++i)
370 {
371 this->quadrature_points[i] = ::Point<1>(p[i]);
372 this->weights[i] = w[i];
373 }
374}
375#endif
376
377
378#ifndef DOXYGEN
379template <>
380QGaussLobatto<1>::QGaussLobatto(const unsigned int n)
381 : Quadrature<1>(n)
382{
383 Assert(n >= 2, ExcNotImplemented());
384
385 std::vector<long double> points =
386 Polynomials::jacobi_polynomial_roots<long double>(n - 2, 1, 1);
387 points.insert(points.begin(), 0);
388 points.push_back(1.);
389 std::vector<long double> w =
391
392 // scale weights to the interval [0.0, 1.0]:
393 for (unsigned int i = 0; i < points.size(); ++i)
394 {
395 this->quadrature_points[i][0] = points[i];
396 this->weights[i] = 0.5 * w[i];
397 }
398}
399#endif
400
401
402template <>
404 : Quadrature<1>(1)
405{
406 this->quadrature_points[0] = Point<1>(0.5);
407 this->weights[0] = 1.0;
408}
409
410
411
412template <>
414 : Quadrature<1>(2)
415{
416 static const double xpts[] = {0.0, 1.0};
417 static const double wts[] = {0.5, 0.5};
418
419 for (unsigned int i = 0; i < this->size(); ++i)
420 {
421 this->quadrature_points[i] = Point<1>(xpts[i]);
422 this->weights[i] = wts[i];
423 }
424}
425
426
427
428template <>
430 : Quadrature<1>(3)
431{
432 static const double xpts[] = {0.0, 0.5, 1.0};
433 static const double wts[] = {1. / 6., 2. / 3., 1. / 6.};
434
435 for (unsigned int i = 0; i < this->size(); ++i)
436 {
437 this->quadrature_points[i] = Point<1>(xpts[i]);
438 this->weights[i] = wts[i];
439 }
440}
441
442
443
444template <>
446 : Quadrature<1>(5)
447{
448 static const double xpts[] = {0.0, .25, .5, .75, 1.0};
449 static const double wts[] = {
450 7. / 90., 32. / 90., 12. / 90., 32. / 90., 7. / 90.};
451
452 for (unsigned int i = 0; i < this->size(); ++i)
453 {
454 this->quadrature_points[i] = Point<1>(xpts[i]);
455 this->weights[i] = wts[i];
456 }
457}
458
459
460
461template <>
463 : Quadrature<1>(7)
464{
465 static const double xpts[] = {
466 0.0, 1. / 6., 1. / 3., .5, 2. / 3., 5. / 6., 1.0};
467 static const double wts[] = {41. / 840.,
468 216. / 840.,
469 27. / 840.,
470 272. / 840.,
471 27. / 840.,
472 216. / 840.,
473 41. / 840.};
474
475 for (unsigned int i = 0; i < this->size(); ++i)
476 {
477 this->quadrature_points[i] = Point<1>(xpts[i]);
478 this->weights[i] = wts[i];
479 }
480}
481
482
483template <>
484QGaussLog<1>::QGaussLog(const unsigned int n, const bool revert)
485 : Quadrature<1>(n)
486{
487 std::vector<double> p = get_quadrature_points(n);
488 std::vector<double> w = get_quadrature_weights(n);
489
490 for (unsigned int i = 0; i < this->size(); ++i)
491 {
492 // Using the change of variables x=1-t, it's possible to show
493 // that int f(x)ln|1-x| = int f(1-t) ln|t|, which implies that
494 // we can use this quadrature formula also with weight ln|1-x|.
495 this->quadrature_points[i] =
496 revert ? Point<1>(1 - p[n - 1 - i]) : Point<1>(p[i]);
497 this->weights[i] = revert ? w[n - 1 - i] : w[i];
498 }
499}
500
501template <>
502std::vector<double>
504{
505 std::vector<double> q_points(n);
506
507 switch (n)
508 {
509 case 1:
510 q_points[0] = 0.3333333333333333;
511 break;
512
513 case 2:
514 q_points[0] = 0.1120088061669761;
515 q_points[1] = 0.6022769081187381;
516 break;
517
518 case 3:
519 q_points[0] = 0.06389079308732544;
520 q_points[1] = 0.3689970637156184;
521 q_points[2] = 0.766880303938942;
522 break;
523
524 case 4:
525 q_points[0] = 0.04144848019938324;
526 q_points[1] = 0.2452749143206022;
527 q_points[2] = 0.5561654535602751;
528 q_points[3] = 0.848982394532986;
529 break;
530
531 case 5:
532 q_points[0] = 0.02913447215197205;
533 q_points[1] = 0.1739772133208974;
534 q_points[2] = 0.4117025202849029;
535 q_points[3] = 0.6773141745828183;
536 q_points[4] = 0.89477136103101;
537 break;
538
539 case 6:
540 q_points[0] = 0.02163400584411693;
541 q_points[1] = 0.1295833911549506;
542 q_points[2] = 0.3140204499147661;
543 q_points[3] = 0.5386572173517997;
544 q_points[4] = 0.7569153373774084;
545 q_points[5] = 0.922668851372116;
546 break;
547
548
549 case 7:
550 q_points[0] = 0.0167193554082585;
551 q_points[1] = 0.100185677915675;
552 q_points[2] = 0.2462942462079286;
553 q_points[3] = 0.4334634932570557;
554 q_points[4] = 0.6323509880476823;
555 q_points[5] = 0.81111862674023;
556 q_points[6] = 0.940848166743287;
557 break;
558
559 case 8:
560 q_points[0] = 0.01332024416089244;
561 q_points[1] = 0.07975042901389491;
562 q_points[2] = 0.1978710293261864;
563 q_points[3] = 0.354153994351925;
564 q_points[4] = 0.5294585752348643;
565 q_points[5] = 0.7018145299391673;
566 q_points[6] = 0.849379320441094;
567 q_points[7] = 0.953326450056343;
568 break;
569
570 case 9:
571 q_points[0] = 0.01086933608417545;
572 q_points[1] = 0.06498366633800794;
573 q_points[2] = 0.1622293980238825;
574 q_points[3] = 0.2937499039716641;
575 q_points[4] = 0.4466318819056009;
576 q_points[5] = 0.6054816627755208;
577 q_points[6] = 0.7541101371585467;
578 q_points[7] = 0.877265828834263;
579 q_points[8] = 0.96225055941096;
580 break;
581
582 case 10:
583 q_points[0] = 0.00904263096219963;
584 q_points[1] = 0.05397126622250072;
585 q_points[2] = 0.1353118246392511;
586 q_points[3] = 0.2470524162871565;
587 q_points[4] = 0.3802125396092744;
588 q_points[5] = 0.5237923179723384;
589 q_points[6] = 0.6657752055148032;
590 q_points[7] = 0.7941904160147613;
591 q_points[8] = 0.898161091216429;
592 q_points[9] = 0.9688479887196;
593 break;
594
595
596 case 11:
597 q_points[0] = 0.007643941174637681;
598 q_points[1] = 0.04554182825657903;
599 q_points[2] = 0.1145222974551244;
600 q_points[3] = 0.2103785812270227;
601 q_points[4] = 0.3266955532217897;
602 q_points[5] = 0.4554532469286375;
603 q_points[6] = 0.5876483563573721;
604 q_points[7] = 0.7139638500230458;
605 q_points[8] = 0.825453217777127;
606 q_points[9] = 0.914193921640008;
607 q_points[10] = 0.973860256264123;
608 break;
609
610 case 12:
611 q_points[0] = 0.006548722279080035;
612 q_points[1] = 0.03894680956045022;
613 q_points[2] = 0.0981502631060046;
614 q_points[3] = 0.1811385815906331;
615 q_points[4] = 0.2832200676673157;
616 q_points[5] = 0.398434435164983;
617 q_points[6] = 0.5199526267791299;
618 q_points[7] = 0.6405109167754819;
619 q_points[8] = 0.7528650118926111;
620 q_points[9] = 0.850240024421055;
621 q_points[10] = 0.926749682988251;
622 q_points[11] = 0.977756129778486;
623 break;
624
625 default:
627 break;
628 }
629
630 return q_points;
631}
632
633
634template <>
635std::vector<double>
637{
638 std::vector<double> quadrature_weights(n);
639
640 switch (n)
641 {
642 case 1:
643 quadrature_weights[0] = -1.0;
644 break;
645 case 2:
646 quadrature_weights[0] = -0.7185393190303845;
647 quadrature_weights[1] = -0.2814606809696154;
648 break;
649
650 case 3:
651 quadrature_weights[0] = -0.5134045522323634;
652 quadrature_weights[1] = -0.3919800412014877;
653 quadrature_weights[2] = -0.0946154065661483;
654 break;
655
656 case 4:
657 quadrature_weights[0] = -0.3834640681451353;
658 quadrature_weights[1] = -0.3868753177747627;
659 quadrature_weights[2] = -0.1904351269501432;
660 quadrature_weights[3] = -0.03922548712995894;
661 break;
662
663 case 5:
664 quadrature_weights[0] = -0.2978934717828955;
665 quadrature_weights[1] = -0.3497762265132236;
666 quadrature_weights[2] = -0.234488290044052;
667 quadrature_weights[3] = -0.0989304595166356;
668 quadrature_weights[4] = -0.01891155214319462;
669 break;
670
671 case 6:
672 quadrature_weights[0] = -0.2387636625785478;
673 quadrature_weights[1] = -0.3082865732739458;
674 quadrature_weights[2] = -0.2453174265632108;
675 quadrature_weights[3] = -0.1420087565664786;
676 quadrature_weights[4] = -0.05545462232488041;
677 quadrature_weights[5] = -0.01016895869293513;
678 break;
679
680
681 case 7:
682 quadrature_weights[0] = -0.1961693894252476;
683 quadrature_weights[1] = -0.2703026442472726;
684 quadrature_weights[2] = -0.239681873007687;
685 quadrature_weights[3] = -0.1657757748104267;
686 quadrature_weights[4] = -0.0889432271377365;
687 quadrature_weights[5] = -0.03319430435645653;
688 quadrature_weights[6] = -0.005932787015162054;
689 break;
690
691 case 8:
692 quadrature_weights[0] = -0.164416604728002;
693 quadrature_weights[1] = -0.2375256100233057;
694 quadrature_weights[2] = -0.2268419844319134;
695 quadrature_weights[3] = -0.1757540790060772;
696 quadrature_weights[4] = -0.1129240302467932;
697 quadrature_weights[5] = -0.05787221071771947;
698 quadrature_weights[6] = -0.02097907374214317;
699 quadrature_weights[7] = -0.003686407104036044;
700 break;
701
702 case 9:
703 quadrature_weights[0] = -0.1400684387481339;
704 quadrature_weights[1] = -0.2097722052010308;
705 quadrature_weights[2] = -0.211427149896601;
706 quadrature_weights[3] = -0.1771562339380667;
707 quadrature_weights[4] = -0.1277992280331758;
708 quadrature_weights[5] = -0.07847890261203835;
709 quadrature_weights[6] = -0.0390225049841783;
710 quadrature_weights[7] = -0.01386729555074604;
711 quadrature_weights[8] = -0.002408041036090773;
712 break;
713
714 case 10:
715 quadrature_weights[0] = -0.12095513195457;
716 quadrature_weights[1] = -0.1863635425640733;
717 quadrature_weights[2] = -0.1956608732777627;
718 quadrature_weights[3] = -0.1735771421828997;
719 quadrature_weights[4] = -0.135695672995467;
720 quadrature_weights[5] = -0.0936467585378491;
721 quadrature_weights[6] = -0.05578772735275126;
722 quadrature_weights[7] = -0.02715981089692378;
723 quadrature_weights[8] = -0.00951518260454442;
724 quadrature_weights[9] = -0.001638157633217673;
725 break;
727
728 case 11:
729 quadrature_weights[0] = -0.1056522560990997;
730 quadrature_weights[1] = -0.1665716806006314;
731 quadrature_weights[2] = -0.1805632182877528;
732 quadrature_weights[3] = -0.1672787367737502;
733 quadrature_weights[4] = -0.1386970574017174;
734 quadrature_weights[5] = -0.1038334333650771;
735 quadrature_weights[6] = -0.06953669788988512;
736 quadrature_weights[7] = -0.04054160079499477;
737 quadrature_weights[8] = -0.01943540249522013;
738 quadrature_weights[9] = -0.006737429326043388;
739 quadrature_weights[10] = -0.001152486965101561;
740 break;
741
742 case 12:
743 quadrature_weights[0] = -0.09319269144393;
744 quadrature_weights[1] = -0.1497518275763289;
745 quadrature_weights[2] = -0.166557454364573;
746 quadrature_weights[3] = -0.1596335594369941;
747 quadrature_weights[4] = -0.1384248318647479;
748 quadrature_weights[5] = -0.1100165706360573;
749 quadrature_weights[6] = -0.07996182177673273;
750 quadrature_weights[7] = -0.0524069547809709;
751 quadrature_weights[8] = -0.03007108900074863;
752 quadrature_weights[9] = -0.01424924540252916;
753 quadrature_weights[10] = -0.004899924710875609;
754 quadrature_weights[11] = -0.000834029009809656;
755 break;
756
757 default:
759 break;
760 }
761
762 return quadrature_weights;
763}
764
765
766template <>
767QGaussLogR<1>::QGaussLogR(const unsigned int n,
768 const Point<1> &origin,
769 const double alpha,
770 const bool factor_out_singularity)
771 : Quadrature<1>(
772 ((origin[0] == 0) || (origin[0] == 1)) ? (alpha == 1 ? n : 2 * n) : 4 * n)
773 , fraction(((origin[0] == 0) || (origin[0] == 1.)) ? 1. : origin[0])
774{
775 // The three quadrature formulas that make this one up. There are
776 // at most two when the origin is one of the extremes, and there is
777 // only one if the origin is one of the extremes and alpha is
778 // equal to one.
779 //
780 // If alpha is different from one, then we need a correction which
781 // is performed with a standard Gauss quadrature rule on each
782 // segment. This is not needed in the standard case where alpha is
783 // equal to one and the origin is on one of the extremes. We
784 // integrate with weight ln|(x-o)/alpha|. In the easy cases, we
785 // only need n quadrature points. In the most difficult one, we
786 // need 2*n points for the first segment, and 2*n points for the
787 // second segment.
788 const QGaussLog<1> quad1(n, origin[0] != 0);
789 const QGaussLog<1> quad2(n);
790 const QGauss<1> quad(n);
791
792 // Check that the origin is inside 0,1
793 Assert((fraction >= 0) && (fraction <= 1),
794 ExcMessage("Origin is outside [0,1]."));
795
796 // Non singular offset. This is the start of non singular quad
797 // points.
798 unsigned int ns_offset = (fraction == 1) ? n : 2 * n;
799
800 for (unsigned int i = 0, j = ns_offset; i < n; ++i, ++j)
801 {
802 // The first i quadrature points are the same as quad1, and
803 // are by default singular.
804 this->quadrature_points[i] = quad1.point(i) * fraction;
805 this->weights[i] = quad1.weight(i) * fraction;
806
807 // We need to scale with -log|fraction*alpha|
808 if ((alpha != 1) || (fraction != 1))
809 {
810 this->quadrature_points[j] = quad.point(i) * fraction;
811 this->weights[j] =
812 -std::log(alpha / fraction) * quad.weight(i) * fraction;
813 }
814 // In case we need the second quadrature as well, do it now.
815 if (fraction != 1)
816 {
817 this->quadrature_points[i + n] =
818 quad2.point(i) * (1 - fraction) + Point<1>(fraction);
819 this->weights[i + n] = quad2.weight(i) * (1 - fraction);
820
821 // We need to scale with -log|fraction*alpha|
822 this->quadrature_points[j + n] =
823 quad.point(i) * (1 - fraction) + Point<1>(fraction);
824 this->weights[j + n] =
825 -std::log(alpha / (1 - fraction)) * quad.weight(i) * (1 - fraction);
826 }
827 }
828 if (factor_out_singularity == true)
829 for (unsigned int i = 0; i < size(); ++i)
830 {
831 Assert(
832 this->quadrature_points[i] != origin,
834 "The singularity cannot be on a Gauss point of the same order!"));
835 double denominator =
836 std::log(std::abs((this->quadrature_points[i] - origin)[0]) / alpha);
837 Assert(denominator != 0.0,
839 "The quadrature formula you are using does not allow to "
840 "factor out the singularity, which is zero at one point."));
841 this->weights[i] /= denominator;
842 }
843}
844
845
846template <>
847unsigned int
848QGaussOneOverR<2>::quad_size(const Point<2> &singularity, const unsigned int n)
849{
850 const double eps = 1e-8;
851 bool on_edge = false;
852 for (unsigned int d = 0; d < 2; ++d)
853 on_edge = on_edge || (std::abs(singularity[d]) < eps ||
854 std::abs(singularity[d] - 1.0) < eps);
855 const bool on_vertex =
856 on_edge &&
857 std::abs((singularity - Point<2>(.5, .5)).norm_square() - .5) < eps;
858 if (on_vertex)
859 return 2 * n * n;
860 else if (on_edge)
861 return 4 * n * n;
862 else
863 return 8 * n * n;
864}
865
866template <>
868 const Point<2> &singularity,
869 const bool factor_out_singularity)
870 : Quadrature<2>(quad_size(singularity, n))
871{
872 // We treat all the cases in the
873 // same way. Split the element in 4
874 // pieces, measure the area, if
875 // it's relevant, add the
876 // quadrature connected to that
877 // singularity.
878 std::vector<QGaussOneOverR<2>> quads;
879 std::vector<Point<2>> origins;
880 // Id of the corner with a
881 // singularity
882 quads.emplace_back(n, 3, factor_out_singularity);
883 quads.emplace_back(n, 2, factor_out_singularity);
884 quads.emplace_back(n, 1, factor_out_singularity);
885 quads.emplace_back(n, 0, factor_out_singularity);
886
887 origins.emplace_back(0., 0.);
888 origins.emplace_back(singularity[0], 0.);
889 origins.emplace_back(0., singularity[1]);
890 origins.push_back(singularity);
891
892 // Lexicographical ordering.
893
894 double eps = 1e-8;
895 unsigned int q_id = 0; // Current quad point index.
896 Tensor<1, 2> dist;
897
898 for (unsigned int box = 0; box < 4; ++box)
899 {
900 dist = (singularity - GeometryInfo<2>::unit_cell_vertex(box));
901 dist = Point<2>(std::abs(dist[0]), std::abs(dist[1]));
902 double area = dist[0] * dist[1];
903 if (area > eps)
904 for (unsigned int q = 0; q < quads[box].size(); ++q, ++q_id)
905 {
906 const Point<2> &qp = quads[box].point(q);
907 this->quadrature_points[q_id] =
908 origins[box] + Point<2>(dist[0] * qp[0], dist[1] * qp[1]);
909 this->weights[q_id] = quads[box].weight(q) * area;
910 }
911 }
912}
913
914
915template <>
917 const unsigned int vertex_index,
918 const bool factor_out_singularity)
919 : Quadrature<2>(2 * n * n)
920{
921 // This version of the constructor
922 // works only for the 4
923 // vertices. If you need a more
924 // general one, you should use the
925 // one with the Point<2> in the
926 // constructor.
927 AssertIndexRange(vertex_index, 4);
928
929 // Start with the gauss quadrature formula on the (u,v) reference
930 // element.
931 const QGauss<2> gauss(n);
932
933 Assert(gauss.size() == n * n, ExcInternalError());
934
935 // For the moment we only implemented this for the vertices of a
936 // quadrilateral. We are planning to do this also for the support
937 // points of arbitrary FE_Q elements, to allow the use of this
938 // class in boundary element programs with higher order mappings.
939 AssertIndexRange(vertex_index, 4);
940
941 // We create only the first one. All other pieces are rotation of
942 // this one.
943 // In this case the transformation is
944 //
945 // (x,y) = (u, u tan(pi/4 v))
946 //
947 // with Jacobian
948 //
949 // J = pi/4 R / cos(pi/4 v)
950 //
951 // And we get rid of R to take into account the singularity,
952 // unless specified differently in the constructor.
953 std::vector<Point<2>> &ps = this->quadrature_points;
954 std::vector<double> &ws = this->weights;
955 double pi4 = numbers::PI / 4;
956
957 for (unsigned int q = 0; q < gauss.size(); ++q)
958 {
959 const Point<2> &gp = gauss.point(q);
960 ps[q][0] = gp[0];
961 ps[q][1] = gp[0] * std::tan(pi4 * gp[1]);
962 ws[q] = gauss.weight(q) * pi4 / std::cos(pi4 * gp[1]);
963 if (factor_out_singularity)
964 ws[q] *= (ps[q] - GeometryInfo<2>::unit_cell_vertex(0)).norm();
965 // The other half of the quadrilateral is symmetric with
966 // respect to xy plane.
967 ws[gauss.size() + q] = ws[q];
968 ps[gauss.size() + q][0] = ps[q][1];
969 ps[gauss.size() + q][1] = ps[q][0];
970 }
971
972 // Now we distribute these vertices in the correct manner
973 double theta = 0;
974 switch (vertex_index)
975 {
976 case 0:
977 theta = 0;
978 break;
979 case 1:
980 //
981 theta = numbers::PI / 2;
982 break;
983 case 2:
984 theta = -numbers::PI / 2;
985 break;
986 case 3:
987 theta = numbers::PI;
988 break;
989 }
990
991 double R00 = std::cos(theta), R01 = -std::sin(theta);
992 double R10 = std::sin(theta), R11 = std::cos(theta);
993
994 if (vertex_index != 0)
995 for (unsigned int q = 0; q < size(); ++q)
996 {
997 double x = ps[q][0] - .5, y = ps[q][1] - .5;
998
999 ps[q][0] = R00 * x + R01 * y + .5;
1000 ps[q][1] = R10 * x + R11 * y + .5;
1001 }
1002}
1003
1004
1005template <int dim>
1007 : Quadrature<dim>(quad)
1008{
1009 std::vector<unsigned int> permutation(quad.size());
1010 for (unsigned int i = 0; i < quad.size(); ++i)
1011 permutation[i] = i;
1012
1013 std::sort(permutation.begin(),
1014 permutation.end(),
1015 [this](const unsigned int x, const unsigned int y) {
1016 return this->compare_weights(x, y);
1017 });
1018
1019 // At this point, the variable is_tensor_product_flag is set
1020 // to the respective value of the given Quadrature in the base
1021 // class copy constructor.
1022 // We only call a quadrature formula 'tensor product'
1023 // if the quadrature points are also sorted lexicographically.
1024 // In particular, any reordering destroys that property
1025 // and we might need to modify the variable accordingly.
1026 for (unsigned int i = 0; i < quad.size(); ++i)
1027 {
1028 this->weights[i] = quad.weight(permutation[i]);
1029 this->quadrature_points[i] = quad.point(permutation[i]);
1030 if (permutation[i] != i)
1031 this->is_tensor_product_flag = false;
1032 }
1033}
1034
1035
1036template <int dim>
1037bool
1038QSorted<dim>::compare_weights(const unsigned int a, const unsigned int b) const
1039{
1040 return (this->weights[a] < this->weights[b]);
1041}
1042
1043
1044// construct the quadrature formulae in higher dimensions by
1045// tensor product of lower dimensions
1046
1047template <int dim>
1048QGauss<dim>::QGauss(const unsigned int n)
1049 : Quadrature<dim>(QGauss<dim - 1>(n), QGauss<1>(n))
1050{}
1051
1052
1053
1054template <int dim>
1055QGaussRadau<dim>::QGaussRadau(const unsigned int n, EndPoint end_point)
1056 : Quadrature<dim>(
1057 QGaussRadau<1>(n, static_cast<QGaussRadau<1>::EndPoint>(end_point)))
1058 , end_point(end_point)
1059{}
1060
1061
1062
1063template <int dim>
1065 : Quadrature<dim>(QGaussLobatto<dim - 1>(n), QGaussLobatto<1>(n))
1066{}
1067
1068
1069
1070template <int dim>
1072 : Quadrature<dim>(QMidpoint<dim - 1>(), QMidpoint<1>())
1073{}
1074
1075
1076
1077template <int dim>
1079 : Quadrature<dim>(QTrapezoid<dim - 1>(), QTrapezoid<1>())
1080{}
1081
1082
1083
1084template <int dim>
1086 : Quadrature<dim>(QSimpson<dim - 1>(), QSimpson<1>())
1087{}
1088
1089
1090
1091template <int dim>
1093 : Quadrature<dim>(QMilne<dim - 1>(), QMilne<1>())
1094{}
1095
1096
1097template <int dim>
1099 : Quadrature<dim>(QWeddle<dim - 1>(), QWeddle<1>())
1100{}
1101
1102template <int dim>
1104 const Point<dim> &singularity)
1105 : // We need the explicit implementation if dim == 1. If dim > 1 we use the
1106 // former implementation and apply a tensorial product to obtain the higher
1107 // dimensions.
1108 Quadrature<dim>(
1109 dim == 2 ?
1110 QAnisotropic<dim>(QTelles<1>(base_quad, Point<1>(singularity[0])),
1111 QTelles<1>(base_quad, Point<1>(singularity[1]))) :
1112 dim == 3 ?
1113 QAnisotropic<dim>(QTelles<1>(base_quad, Point<1>(singularity[0])),
1114 QTelles<1>(base_quad, Point<1>(singularity[1])),
1115 QTelles<1>(base_quad, Point<1>(singularity[2]))) :
1116 Quadrature<dim>())
1117{}
1118
1119template <int dim>
1120QTelles<dim>::QTelles(const unsigned int n, const Point<dim> &singularity)
1121 : // In this case we map the standard Gauss Legendre formula using the given
1122 // singularity point coordinates.
1123 Quadrature<dim>(QTelles<dim>(QGauss<1>(n), singularity))
1124{}
1125
1126
1127
1128template <>
1129QTelles<1>::QTelles(const Quadrature<1> &base_quad, const Point<1> &singularity)
1130 : // We explicitly implement the Telles' variable change if dim == 1.
1131 Quadrature<1>(base_quad)
1132{
1133 // We define all the constants to be used in the implementation of
1134 // Telles' rule
1135 const double eta_bar = singularity[0] * 2. - 1.;
1136 const double eta_star = eta_bar * eta_bar - 1.;
1137 double gamma_bar;
1138
1139 std::vector<Point<1>> quadrature_points_dummy(quadrature_points.size());
1140 std::vector<double> weights_dummy(weights.size());
1141 unsigned int cont = 0;
1142 const double tol = 1e-10;
1143 for (unsigned int d = 0; d < quadrature_points.size(); ++d)
1144 {
1145 if (std::abs(quadrature_points[d][0] - singularity[0]) > tol)
1146 {
1147 quadrature_points_dummy[d - cont] = quadrature_points[d];
1148 weights_dummy[d - cont] = weights[d];
1149 }
1150 else
1151 {
1152 // We need to remove the singularity point from the quadrature point
1153 // list. To do so we use the variable cont.
1154 cont = 1;
1155 }
1156 }
1157 if (cont == 1)
1158 {
1159 quadrature_points.resize(quadrature_points_dummy.size() - 1);
1160 weights.resize(weights_dummy.size() - 1);
1161 for (unsigned int d = 0; d < quadrature_points.size(); ++d)
1162 {
1163 quadrature_points[d] = quadrature_points_dummy[d];
1164 weights[d] = weights_dummy[d];
1165 }
1166 }
1167 // We need to check if the singularity is at the boundary of the interval.
1168 if (std::abs(eta_star) <= tol)
1169 {
1170 gamma_bar = std::cbrt(eta_bar * eta_star + std::abs(eta_star)) +
1171 std::cbrt(eta_bar * eta_star - std::abs(eta_star)) + eta_bar;
1172 }
1173 else
1174 {
1175 gamma_bar =
1176 (eta_bar * eta_star + std::abs(eta_star)) /
1177 std::abs(eta_bar * eta_star + std::abs(eta_star)) *
1178 std::cbrt(std::abs(eta_bar * eta_star + std::abs(eta_star))) +
1179 (eta_bar * eta_star - std::abs(eta_star)) /
1180 std::abs(eta_bar * eta_star - std::abs(eta_star)) *
1181 std::cbrt(std::abs(eta_bar * eta_star - std::abs(eta_star))) +
1182 eta_bar;
1183 }
1184 for (unsigned int q = 0; q < quadrature_points.size(); ++q)
1185 {
1186 double gamma = quadrature_points[q][0] * 2 - 1;
1187 double eta = (Utilities::fixed_power<3>(gamma - gamma_bar) +
1188 gamma_bar * (gamma_bar * gamma_bar + 3)) /
1189 (1 + 3 * gamma_bar * gamma_bar);
1190
1191 double J = 3 * ((gamma - gamma_bar) * (gamma - gamma_bar)) /
1192 (1 + 3 * gamma_bar * gamma_bar);
1193
1194 quadrature_points[q][0] = (eta + 1) / 2.0;
1195 weights[q] = J * weights[q];
1196 }
1197}
1198
1199namespace internal
1200{
1202 {
1206 std::vector<double>
1207 get_quadrature_points(const unsigned int n)
1208 {
1209 std::vector<double> points(n);
1210 // n point quadrature: index from 0 to n-1
1211 for (unsigned short i = 0; i < n; ++i)
1212 // would be cos((2i+1)Pi/(2N+2))
1213 // put + Pi so we start from the smallest point
1214 // then map from [-1,1] to [0,1]
1215 points[i] =
1216 1. / 2. *
1217 (1. + std::cos(numbers::PI *
1218 (1. + double(2 * i + 1) / double(2 * (n - 1) + 2))));
1219
1220 return points;
1221 }
1222
1223
1224
1228 std::vector<double>
1229 get_quadrature_weights(const unsigned int n)
1230 {
1231 std::vector<double> weights(n);
1232
1233 for (unsigned short i = 0; i < n; ++i)
1234 {
1235 // same weights as on [-1,1]
1236 weights[i] = numbers::PI / double(n);
1237 }
1238
1239 return weights;
1240 }
1241 } // namespace QGaussChebyshev
1242} // namespace internal
1243
1244
1245template <>
1247 : Quadrature<1>(n)
1248{
1249 Assert(n > 0, ExcMessage("Need at least one point for the quadrature rule"));
1250 std::vector<double> p = internal::QGaussChebyshev::get_quadrature_points(n);
1251 std::vector<double> w = internal::QGaussChebyshev::get_quadrature_weights(n);
1252
1253 for (unsigned int i = 0; i < this->size(); ++i)
1254 {
1255 this->quadrature_points[i] = Point<1>(p[i]);
1256 this->weights[i] = w[i];
1257 }
1258}
1259
1260
1261template <int dim>
1263 : Quadrature<dim>(QGaussChebyshev<1>(n))
1264{}
1265
1266
1267namespace internal
1268{
1270 {
1271 // Computes the points of the quadrature formula.
1272 std::vector<double>
1274 const unsigned int n,
1275 const ::QGaussRadauChebyshev<1>::EndPoint end_point)
1276 {
1277 std::vector<double> points(n);
1278 // n point quadrature: index from 0 to n-1
1279 for (unsigned short i = 0; i < n; ++i)
1280 // would be -cos(2i Pi/(2N+1))
1281 // put + Pi so we start from the smallest point
1282 // then map from [-1,1] to [0,1]
1283 switch (end_point)
1284 {
1285 case ::QGaussRadauChebyshev<1>::EndPoint::left:
1286 {
1287 points[i] =
1288 1. / 2. *
1289 (1. -
1291 (1 + 2 * double(i) / (2 * double(n - 1) + 1.))));
1292 break;
1293 }
1294
1295 case ::QGaussRadauChebyshev<1>::EndPoint::right:
1296 {
1297 points[i] =
1298 1. / 2. *
1299 (1. - std::cos(numbers::PI * (2 * double(n - 1 - i) /
1300 (2 * double(n - 1) + 1.))));
1301 break;
1302 }
1303
1304 default:
1305 Assert(
1306 false,
1307 ExcMessage(
1308 "This constructor can only be called with either "
1309 "QGaussRadauChebyshev::EndPoint::left or "
1310 "QGaussRadauChebyshev::EndPoint:right as second argument."));
1311 }
1312
1313 return points;
1314 }
1315
1316
1317
1318 // Computes the weights of the quadrature formula.
1319 std::vector<double>
1321 const unsigned int n,
1322 const ::QGaussRadauChebyshev<1>::EndPoint end_point)
1323 {
1324 std::vector<double> weights(n);
1325
1326 for (unsigned short i = 0; i < n; ++i)
1327 {
1328 // same weights as on [-1,1]
1329 weights[i] = 2. * numbers::PI / double(2 * (n - 1) + 1.);
1331 i == 0)
1332 weights[i] /= 2.;
1333 else if (end_point ==
1335 i == (n - 1))
1336 weights[i] /= 2.;
1337 }
1338
1339 return weights;
1340 }
1341 } // namespace QGaussRadauChebyshev
1342} // namespace internal
1343
1344
1345template <>
1347 const EndPoint end_point)
1348 : Quadrature<1>(n)
1349 , end_point(end_point)
1350{
1351 Assert(n > 0, ExcMessage("Need at least one point for quadrature rules."));
1352 std::vector<double> points =
1354 std::vector<double> new_weights =
1356
1357 for (unsigned int i = 0; i < this->size(); ++i)
1358 {
1359 this->quadrature_points[i] = Point<1>(points[i]);
1360 this->weights[i] = new_weights[i];
1361 }
1362}
1363
1364
1365template <int dim>
1367 EndPoint end_point)
1369 n,
1370 static_cast<QGaussRadauChebyshev<1>::EndPoint>(end_point)))
1371 , end_point(end_point)
1372{}
1373
1374
1375
1376namespace internal
1377{
1379 {
1380 // Computes the points of the quadrature formula.
1381 std::vector<double>
1382 get_quadrature_points(const unsigned int n)
1383 {
1384 std::vector<double> points(n);
1385 // n point quadrature: index from 0 to n-1
1386 for (unsigned short i = 0; i < n; ++i)
1387 // would be cos(i Pi/N)
1388 // put + Pi so we start from the smallest point
1389 // then map from [-1,1] to [0,1]
1390 points[i] =
1391 1. / 2. *
1392 (1. + std::cos(numbers::PI * (1 + double(i) / double(n - 1))));
1393
1394 return points;
1395 }
1396
1397 // Computes the weights of the quadrature formula.
1398 std::vector<double>
1399 get_quadrature_weights(const unsigned int n)
1400 {
1401 std::vector<double> weights(n);
1402
1403 for (unsigned short i = 0; i < n; ++i)
1404 {
1405 // same weights as on [-1,1]
1406 weights[i] = numbers::PI / double((n - 1));
1407 if (i == 0 || i == (n - 1))
1408 weights[i] /= 2.;
1409 }
1410
1411 return weights;
1412 }
1413 } // namespace QGaussLobattoChebyshev
1414} // namespace internal
1415
1416
1417
1418template <>
1420 : Quadrature<1>(n)
1421{
1422 Assert(n > 1,
1423 ExcMessage(
1424 "Need at least two points for Gauss-Lobatto quadrature rule"));
1425 std::vector<double> p =
1427 std::vector<double> w =
1429
1430 for (unsigned int i = 0; i < this->size(); ++i)
1431 {
1432 this->quadrature_points[i] = Point<1>(p[i]);
1433 this->weights[i] = w[i];
1434 }
1435}
1436
1437
1438template <int dim>
1442
1443
1444
1445template <int dim>
1447{
1448 std::vector<Point<dim>> qpoints;
1449 std::vector<double> weights;
1450
1451 for (unsigned int i = 0; i < quad.size(); ++i)
1452 {
1453 double r = 0;
1454 /* Use "int d" instead of the more natural "unsigned int d" to work
1455 * around a wrong diagnostic in gcc-10.3.0 that warns about that the
1456 * comparison "d < dim" is always false in case of "dim == 0".
1457 * MM 2021 */
1458 for (int d = 0; d < dim; ++d)
1459 r += quad.point(i)[d];
1460 if (r <= 1 + 1e-10)
1461 {
1462 this->quadrature_points.push_back(quad.point(i));
1463 this->weights.push_back(quad.weight(i));
1464 }
1465 }
1466}
1467
1468
1469
1470template <int dim>
1471template <int spacedim>
1474 const std::array<Point<spacedim>, dim + 1> &vertices) const
1475{
1476 Assert(dim <= spacedim,
1477 ExcMessage("Invalid combination of dim and spacedim ."));
1479 for (unsigned int d = 0; d < dim; ++d)
1480 Bt[d] = vertices[d + 1] - vertices[0];
1481
1482 const auto B = Bt.transpose();
1483 const double J = std::abs(B.determinant());
1484
1485 // if the determinant is zero, we return an empty quadrature
1486 if (J < 1e-12)
1487 return Quadrature<spacedim>();
1488
1489 std::vector<Point<spacedim>> qp(this->size());
1490 std::vector<double> w(this->size());
1491
1492 for (unsigned int i = 0; i < this->size(); ++i)
1493 {
1494 qp[i] =
1495 Point<spacedim>(vertices[0] + apply_transformation(B, this->point(i)));
1496 w[i] = J * this->weight(i);
1497 }
1498
1499 return Quadrature<spacedim>(qp, w);
1500}
1501
1502
1503
1504template <int dim>
1505template <int spacedim>
1508 const std::vector<std::array<Point<spacedim>, dim + 1>> &simplices) const
1509{
1510 Assert(!(dim == 1 && spacedim == 1),
1511 ExcMessage("This function is not supposed to work in 1D-1d case."));
1512 Assert(dim <= spacedim,
1513 ExcMessage("Invalid combination of dim and spacedim ."));
1514
1515 std::vector<Point<spacedim>> qp;
1516 std::vector<double> ws;
1517 for (const auto &simplex : simplices)
1518 {
1519 const auto rule = this->compute_affine_transformation(simplex);
1520 std::transform(rule.get_points().begin(),
1521 rule.get_points().end(),
1522 std::back_inserter(qp),
1523 [&](const Point<spacedim> &p) { return p; });
1524 std::transform(rule.get_weights().begin(),
1525 rule.get_weights().end(),
1526 std::back_inserter(ws),
1527 [&](const double w) { return w; });
1528 }
1529 return Quadrature<spacedim>(qp, ws);
1530}
1531
1532
1533
1535 const Quadrature<1> &angular_quadrature)
1536 : QSimplex<2>(Quadrature<2>())
1537{
1538 const QAnisotropic<2> base(radial_quadrature, angular_quadrature);
1539 this->quadrature_points.resize(base.size());
1540 this->weights.resize(base.size());
1541 for (unsigned int i = 0; i < base.size(); ++i)
1542 {
1543 const auto &q = base.point(i);
1544 const auto w = base.weight(i);
1545
1546 const auto xhat = q[0];
1547 const auto yhat = q[1];
1548
1549 const double t = numbers::PI_2 * yhat;
1550 const double pi = numbers::PI;
1551 const double st = std::sin(t);
1552 const double ct = std::cos(t);
1553 const double r = xhat / (st + ct);
1554
1555 const double J = pi * xhat / (2 * (std::sin(pi * yhat) + 1));
1556
1557 this->quadrature_points[i] = Point<2>(r * ct, r * st);
1558 this->weights[i] = w * J;
1559 }
1560}
1561
1562
1563
1565 : QTrianglePolar(QGauss<1>(n), QGauss<1>(n))
1566{}
1567
1568
1569
1570QDuffy::QDuffy(const Quadrature<1> &radial_quadrature,
1571 const Quadrature<1> &angular_quadrature,
1572 const double beta)
1573 : QSimplex<2>(Quadrature<2>())
1574{
1575 const QAnisotropic<2> base(radial_quadrature, angular_quadrature);
1576 this->quadrature_points.resize(base.size());
1577 this->weights.resize(base.size());
1578 for (unsigned int i = 0; i < base.size(); ++i)
1579 {
1580 const auto &q = base.point(i);
1581 const auto w = base.weight(i);
1582
1583 const auto xhat = q[0];
1584 const auto yhat = q[1];
1585
1586 const double x = std::pow(xhat, beta) * (1 - yhat);
1587 const double y = std::pow(xhat, beta) * yhat;
1588
1589 const double J = beta * std::pow(xhat, 2. * beta - 1.);
1590
1591 this->quadrature_points[i] = Point<2>(x, y);
1592 this->weights[i] = w * J;
1593 }
1594}
1595
1596
1597
1598QDuffy::QDuffy(const unsigned int n, const double beta)
1599 : QDuffy(QGauss<1>(n), QGauss<1>(n), beta)
1600{}
1601
1602
1603
1604template <int dim>
1605QSplit<dim>::QSplit(const QSimplex<dim> &base, const Point<dim> &split_point)
1606{
1608 ExcMessage(
1609 "The split point should be inside the unit reference cell."));
1610
1611 std::array<Point<dim>, dim + 1> vertices;
1612 vertices[0] = split_point;
1613
1614 // Make a simplex from the split_point and the first dim vertices of each
1615 // face. In dimension three, we need to split the face in two triangles, so
1616 // we use once the first dim vertices of each face, and the second time the
1617 // the dim vertices of each face starting from 1.
1618 for (auto f : GeometryInfo<dim>::face_indices())
1619 for (unsigned int start = 0; start < (dim > 2 ? 2 : 1); ++start)
1620 {
1621 for (unsigned int i = 0; i < dim; ++i)
1622 vertices[i + 1] = GeometryInfo<dim>::unit_cell_vertex(
1624 const auto quad = base.compute_affine_transformation(vertices);
1625 if (quad.size())
1626 {
1627 this->quadrature_points.insert(this->quadrature_points.end(),
1628 quad.get_points().begin(),
1629 quad.get_points().end());
1630 this->weights.insert(this->weights.end(),
1631 quad.get_weights().begin(),
1632 quad.get_weights().end());
1633 }
1634 }
1635}
1636
1637
1638
1639template <int dim>
1640QGaussSimplex<dim>::QGaussSimplex(const unsigned int n_points_1D)
1641 : QSimplex<dim>(Quadrature<dim>())
1642{
1643 // fill quadrature points and quadrature weights
1644 if (dim == 0 || dim == 1)
1645 {
1646 const ::QGauss<dim> quad(n_points_1D);
1647
1648 this->quadrature_points = quad.get_points();
1649 this->weights = quad.get_weights();
1650 }
1651 else if (dim == 2)
1652 {
1653 if (n_points_1D == 1)
1654 {
1655 const double p = 1.0 / 3.0;
1656 this->quadrature_points.emplace_back(p, p);
1657 this->weights.emplace_back(0.5);
1658 }
1659 else if (n_points_1D == 2)
1660 {
1661 // The Hillion 7 scheme, as communicated by quadpy
1662 //
1663 // See: Numerical Integration on a Triangle, International Journal for
1664 // Numerical Methods in Engineering, 1977
1665 const double Q12 = 1.0 / 2.0;
1666 this->quadrature_points.emplace_back(0.17855872826361643,
1667 0.1550510257216822);
1668 this->quadrature_points.emplace_back(0.07503111022260812,
1669 0.6449489742783178);
1670 this->quadrature_points.emplace_back(0.6663902460147014,
1671 0.1550510257216822);
1672 this->quadrature_points.emplace_back(0.28001991549907407,
1673 0.6449489742783178);
1674
1675 this->weights.emplace_back(0.31804138174397717 * Q12);
1676 this->weights.emplace_back(0.18195861825602283 * Q12);
1677 this->weights.emplace_back(0.31804138174397717 * Q12);
1678 this->weights.emplace_back(0.18195861825602283 * Q12);
1679 }
1680 else if (n_points_1D == 3)
1681 {
1682 // The Hammer-Marlowe-Stroud 5 Scheme, as communicated by quadpy
1683 const double p0 = 2.0 / 7.0 - std::sqrt(15.0) / 21.0;
1684 const double p1 = 2.0 / 7.0 + std::sqrt(15.0) / 21.0;
1685 const double p2 = 3.0 / 7.0 - 2.0 * std::sqrt(15.0) / 21.0;
1686 const double p3 = 3.0 / 7.0 + 2.0 * std::sqrt(15.0) / 21.0;
1687 this->quadrature_points.emplace_back(1.0 / 3.0, 1.0 / 3.0);
1688 this->quadrature_points.emplace_back(p3, p0);
1689 this->quadrature_points.emplace_back(p0, p3);
1690 this->quadrature_points.emplace_back(p0, p0);
1691 this->quadrature_points.emplace_back(p2, p1);
1692 this->quadrature_points.emplace_back(p1, p2);
1693 this->quadrature_points.emplace_back(p1, p1);
1694
1695 const double q12 = 0.5;
1696 const double w0 = 9.0 / 40.0;
1697 const double w1 = 31.0 / 240.0 - std::sqrt(15.0) / 1200.0;
1698 const double w2 = 31.0 / 240.0 + std::sqrt(15.0) / 1200.0;
1699 this->weights.emplace_back(q12 * w0);
1700 this->weights.emplace_back(q12 * w1);
1701 this->weights.emplace_back(q12 * w1);
1702 this->weights.emplace_back(q12 * w1);
1703 this->weights.emplace_back(q12 * w2);
1704 this->weights.emplace_back(q12 * w2);
1705 this->weights.emplace_back(q12 * w2);
1706 }
1707 else if (n_points_1D == 4)
1708 {
1710 QWitherdenVincentSimplex<dim>(n_points_1D));
1711 }
1712 }
1713 else if (dim == 3)
1714 {
1715 if (n_points_1D == 1)
1716 {
1717 const double Q14 = 1.0 / 4.0;
1718 const double Q16 = 1.0 / 6.0;
1719
1720 this->quadrature_points.emplace_back(Q14, Q14, Q14);
1721 this->weights.emplace_back(Q16);
1722 }
1723 // The Xiao Gimbutas 03 scheme, as communicated by quadpy
1724 //
1725 // See: A numerical algorithm for the construction of efficient quadrature
1726 // rules in two and higher dimensions, Computers & Mathematics with
1727 // Applications, 2010
1728 else if (n_points_1D == 2)
1729 {
1730 const double Q16 = 1.0 / 6.0;
1731 this->weights.emplace_back(0.1223220027573451 * Q16);
1732 this->weights.emplace_back(0.1280664127107469 * Q16);
1733 this->weights.emplace_back(0.1325680271444452 * Q16);
1734 this->weights.emplace_back(0.1406244096604032 * Q16);
1735 this->weights.emplace_back(0.2244151669175574 * Q16);
1736 this->weights.emplace_back(0.2520039808095023 * Q16);
1737
1738 this->quadrature_points.emplace_back(0.1620014916985245,
1739 0.1838503504920977,
1740 0.01271836631368145);
1741 this->quadrature_points.emplace_back(0.01090521221118924,
1742 0.2815238021235462,
1743 0.3621268299455338);
1744 this->quadrature_points.emplace_back(0.1901170024392839,
1745 0.01140332944455717,
1746 0.3586207204668839);
1747 this->quadrature_points.emplace_back(0.170816925164989,
1748 0.1528181430909273,
1749 0.6384932999617267);
1750 this->quadrature_points.emplace_back(0.1586851632274406,
1751 0.5856628056552158,
1752 0.1308471689520965);
1753 this->quadrature_points.emplace_back(0.5712260521491151,
1754 0.1469183900871696,
1755 0.1403728057942107);
1756 }
1757 // Past this point the best rules (positive weights, minimal number of
1758 // points) we have right now are the Witherden-Vincent ones
1759 else if (n_points_1D == 3)
1760 {
1762 QWitherdenVincentSimplex<dim>(n_points_1D));
1763 }
1764 else if (n_points_1D == 4)
1765 {
1767 QWitherdenVincentSimplex<dim>(n_points_1D));
1768 }
1769 }
1770
1771 AssertDimension(this->quadrature_points.size(), this->weights.size());
1772 Assert(this->quadrature_points.size() > 0,
1774 "QGaussSimplex is currently only implemented for "
1775 "n_points_1D = 1, 2, 3, and 4 while you are asking for "
1776 "n_points_1D = " +
1777 Utilities::to_string(n_points_1D)));
1778}
1779
1780namespace
1781{
1782 template <std::size_t b_dim>
1783 std::vector<std::array<double, b_dim>>
1784 all_permutations(const std::array<double, b_dim> &b_point)
1785 {
1786 std::vector<std::array<double, b_dim>> output;
1787
1788 // We want all possible permutations of the barycentric coordinates.
1789 // The easiest way to get all of them is to sort the input first and
1790 // then use next_permutation to cycle through them all.
1791 std::array<double, b_dim> temp = b_point;
1792 std::sort(temp.begin(), temp.end());
1793 do
1794 {
1795 output.push_back(temp);
1796 }
1797 while (std::next_permutation(temp.begin(), temp.end()));
1798
1799 return output;
1800 }
1801} // namespace
1802
1803
1804
1805template <int dim>
1807 const unsigned int n_points_1D,
1808 const bool use_odd_order)
1809 : QSimplex<dim>(Quadrature<dim>())
1810{
1811 Assert(1 <= dim && dim <= 3, ExcNotImplemented());
1812 // Just use Gauss in 1d: this is a high-order open rule so this is a
1813 // reasonable equivalent for generic programming.
1814 if (dim == 1)
1815 {
1817 return;
1818 }
1819
1820 std::array<double, dim + 1> centroid;
1821 std::fill(centroid.begin(), centroid.end(), 1.0 / (dim + 1.0));
1822 std::vector<std::vector<std::array<double, dim + 1>>> b_point_permutations;
1823 std::vector<double> b_weights;
1824
1825 // We can simplify the implementation of these quadrature rules
1826 // by quite a bit by exploiting symmetry - we do essentially the
1827 // same thing for each barycentric coordinate, so we can express
1828 // our quadrature rule as permutations of barycentric points
1829 // instead of writing things out explicitly.
1830
1831 // Apply a Barycentric permutation where one point is different.
1832 // Equivalent to d3_aa and s31 in quadpy.
1833 auto process_point_1 = [&](const double a, const double w) {
1834 const double b = 1.0 - dim * a;
1835 std::array<double, dim + 1> b_point;
1836 std::fill(b_point.begin(), b_point.begin() + dim, a);
1837 b_point[dim] = b;
1838
1839 b_weights.push_back(w);
1840 b_point_permutations.push_back(all_permutations(b_point));
1841 };
1842
1843 // Apply a Barycentric permutation where two points (in 3d) are different.
1844 // Equivalent to s22 in quadpy.
1845 auto process_point_2 = [&](const double a, const double w) {
1846 Assert(dim == 3, ExcInternalError());
1847 const double b = (1.0 - 2.0 * a) / 2.0;
1848 std::array<double, dim + 1> b_point;
1849 std::fill(b_point.begin(), b_point.begin() + dim - 1, a);
1850 b_point[dim - 1] = b;
1851 b_point[dim] = b;
1852
1853 b_weights.push_back(w);
1854 b_point_permutations.push_back(all_permutations(b_point));
1855 };
1856
1857 // Apply a Barycentric permutation where three (or four) points
1858 // are different (since there are two inputs).
1859 // Equivalent to d3_ab and s211 in quadpy.
1860 auto process_point_3 = [&](const double a, const double b, const double w) {
1861 const double c = 1.0 - (dim - 1.0) * a - b;
1862 std::array<double, dim + 1> b_point;
1863 std::fill(b_point.begin(), b_point.begin() + dim - 1, a);
1864 b_point[dim - 1] = b;
1865 b_point[dim] = c;
1866
1867 b_weights.push_back(w);
1868 b_point_permutations.push_back(all_permutations(b_point));
1869 };
1870
1871 switch (n_points_1D)
1872 {
1873 case 1:
1874 switch (dim)
1875 {
1876 case 2:
1877 if (use_odd_order)
1878 {
1879 // WV-1, 2d
1880 b_point_permutations.push_back({centroid});
1881 b_weights.push_back(1.0000000000000000e+00);
1882 }
1883 else
1884 {
1885 // WV-2, 2d
1886 process_point_1(1.6666666666666669e-01,
1887 3.3333333333333331e-01);
1888 }
1889 break;
1890 case 3:
1891 if (use_odd_order)
1892 {
1893 // WV-1, 3d
1894 b_point_permutations.push_back({centroid});
1895 b_weights.push_back(1.0000000000000000e+00);
1896 }
1897 else
1898 {
1899 // WV-2, 3d
1900 process_point_1(1.3819660112501050e-01,
1901 2.5000000000000000e-01);
1902 }
1903 break;
1904 default:
1906 }
1907 break;
1908 case 2:
1909 switch (dim)
1910 {
1911 case 2:
1912 // WV-4 in both cases (no WV-3 in 2d)
1913 process_point_1(9.1576213509770743e-02, 1.0995174365532187e-01);
1914 process_point_1(4.4594849091596489e-01, 2.2338158967801147e-01);
1915 break;
1916 case 3:
1917 if (use_odd_order)
1918 {
1919 // WV-3, 3d
1920 process_point_1(3.2816330251638171e-01,
1921 1.3621784253708741e-01);
1922 process_point_1(1.0804724989842859e-01,
1923 1.1378215746291261e-01);
1924 }
1925 else
1926 {
1927 // WV-5 (no WV-4 in 3d)
1929 }
1930 break;
1931 default:
1933 }
1934 break;
1935 case 3:
1936 switch (dim)
1937 {
1938 case 2:
1939 if (use_odd_order)
1940 {
1941 // WV-5, 2d
1942 b_point_permutations.push_back({centroid});
1943 b_weights.push_back(2.2500000000000001e-01);
1944 process_point_1(1.0128650732345634e-01,
1945 1.2593918054482714e-01);
1946 process_point_1(4.7014206410511511e-01,
1947 1.3239415278850619e-01);
1948 }
1949 else
1950 {
1951 // WV-6, 2d
1952 process_point_1(6.3089014491502227e-02,
1953 5.0844906370206819e-02);
1954 process_point_1(2.4928674517091043e-01,
1955 1.1678627572637937e-01);
1956 process_point_3(5.3145049844816938e-02,
1957 3.1035245103378439e-01,
1958 8.2851075618373571e-02);
1959 }
1960 break;
1961 case 3:
1962 if (use_odd_order)
1963 {
1964 // WV-5, 3d
1965 process_point_1(3.1088591926330061e-01,
1966 1.1268792571801590e-01);
1967 process_point_1(9.2735250310891248e-02,
1968 7.3493043116361956e-02);
1969 process_point_2(4.5503704125649642e-02,
1970 4.2546020777081472e-02);
1971 }
1972 else
1973 {
1974 // WV-6, 3d
1975 process_point_1(4.0673958534611372e-02,
1976 1.0077211055320640e-02);
1977 process_point_1(3.2233789014227548e-01,
1978 5.5357181543654717e-02);
1979 process_point_1(2.1460287125915201e-01,
1980 3.9922750258167487e-02);
1981 process_point_3(6.3661001875017442e-02,
1982 6.0300566479164919e-01,
1983 4.8214285714285710e-02);
1984 }
1985 break;
1986 default:
1988 }
1989 break;
1990 case 4:
1991 switch (dim)
1992 {
1993 case 2:
1994 if (use_odd_order)
1995 {
1996 // WV-7, 2d
1997 process_point_1(3.3730648554587850e-02,
1998 1.6545050110792131e-02);
1999 process_point_1(4.7430969250471822e-01,
2000 7.7086646185986069e-02);
2001 process_point_1(2.4157738259540357e-01,
2002 1.2794417123015558e-01);
2003 process_point_3(4.7036644652595216e-02,
2004 1.9868331479735168e-01,
2005 5.5878732903199779e-02);
2006 }
2007 else
2008 {
2009 // WV-8, 2d
2010 b_point_permutations.push_back({centroid});
2011 b_weights.push_back(1.4431560767778717e-01);
2012 process_point_1(5.0547228317030957e-02,
2013 3.2458497623198079e-02);
2014 process_point_1(4.5929258829272313e-01,
2015 9.5091634267284619e-02);
2016 process_point_1(1.7056930775176021e-01,
2017 1.0321737053471824e-01);
2018 process_point_3(8.3947774099575878e-03,
2019 2.6311282963463811e-01,
2020 2.7230314174434993e-02);
2021 }
2022 break;
2023 case 3:
2024 if (use_odd_order)
2025 {
2026 // WV-7, 3d
2027 b_point_permutations.push_back({centroid});
2028 b_weights.push_back(9.5485289464130846e-02);
2029 process_point_1(3.1570114977820279e-01,
2030 4.2329581209967028e-02);
2031 process_point_2(5.0489822598396350e-02,
2032 3.1896927832857580e-02);
2033 process_point_3(1.8883383102600099e-01,
2034 5.7517163758699996e-01,
2035 3.7207130728334620e-02);
2036 process_point_3(2.1265472541483140e-02,
2037 8.1083024109854862e-01,
2038 8.1107708299033420e-03);
2039 }
2040 else
2041 {
2042 // WV-8, 3d
2043 process_point_1(1.0795272496221089e-01,
2044 2.6426650908408830e-02);
2045 process_point_1(1.8510948778258660e-01,
2046 5.2031747563738531e-02);
2047 process_point_1(4.2316543684767283e-02,
2048 7.5252561535401989e-03);
2049 process_point_1(3.1418170912403898e-01,
2050 4.1763782856934897e-02);
2051 process_point_2(4.3559132858383021e-01,
2052 3.6280930261308818e-02);
2053 process_point_3(2.1433930127130570e-02,
2054 7.1746406342630831e-01,
2055 7.1569028908444327e-03);
2056 process_point_3(2.0413933387602909e-01,
2057 5.8379737830214440e-01,
2058 1.5453486150960340e-02);
2059 }
2060 break;
2061 default:
2063 }
2064 break;
2065 case 5:
2066 switch (dim)
2067 {
2068 case 2:
2069 if (use_odd_order)
2070 {
2071 // WV-9, 2d
2072 b_point_permutations.push_back({centroid});
2073 b_weights.push_back(9.7135796282798836e-02);
2074 process_point_1(4.4729513394452691e-02,
2075 2.5577675658698031e-02);
2076 process_point_1(4.8968251919873762e-01,
2077 3.1334700227139071e-02);
2078 process_point_1(4.3708959149293664e-01,
2079 7.7827541004774278e-02);
2080 process_point_1(1.8820353561903275e-01,
2081 7.9647738927210249e-02);
2082 process_point_3(3.6838412054736258e-02,
2083 2.2196298916076568e-01,
2084 4.3283539377289376e-02);
2085 }
2086 else
2087 {
2088 // WV-10, 2d
2089 b_point_permutations.push_back({centroid});
2090 b_weights.push_back(8.1743329146285973e-02);
2091 process_point_1(3.2055373216943517e-02,
2092 1.3352968813149567e-02);
2093 process_point_1(1.4216110105656438e-01,
2094 4.5957963604744731e-02);
2095 process_point_3(2.8367665339938453e-02,
2096 1.6370173373718250e-01,
2097 2.5297757707288385e-02);
2098 process_point_3(2.9619889488729734e-02,
2099 3.6914678182781102e-01,
2100 3.4184648162959429e-02);
2101 process_point_3(1.4813288578382056e-01,
2102 3.2181299528883545e-01,
2103 6.3904906396424044e-02);
2104 }
2105 break;
2106 case 3:
2107 if (use_odd_order)
2108 {
2109 // WV-9, 3d
2110 b_point_permutations.push_back({centroid});
2111 b_weights.push_back(5.8010548912480253e-02);
2112 process_point_1(6.1981697552226933e-10,
2113 6.4319281759256394e-05);
2114 process_point_1(1.6077453539526160e-01,
2115 2.3173338462425461e-02);
2116 process_point_1(3.2227652182142102e-01,
2117 2.9562912335429289e-02);
2118 process_point_1(4.5108918345413578e-02,
2119 8.0639799796161822e-03);
2120 process_point_2(1.1229654600437609e-01,
2121 3.8134080103702457e-02);
2122 process_point_3(4.5887144875245922e-01,
2123 2.5545792330413102e-03,
2124 8.3844221982985519e-03);
2125 process_point_3(3.3775870685338598e-02,
2126 7.1835032644207453e-01,
2127 1.0234559352745330e-02);
2128 process_point_3(1.8364136980992790e-01,
2129 3.4415910578175279e-02,
2130 2.0524915967988139e-02);
2131 }
2132 else
2133 {
2134 // WV-10, 3d
2135 b_point_permutations.push_back({centroid});
2136 b_weights.push_back(4.7399773556020743e-02);
2137 process_point_1(3.1225006869518868e-01,
2138 2.6937059992268701e-02);
2139 process_point_1(1.1430965385734609e-01,
2140 9.8691597167933822e-03);
2141 process_point_3(4.1043073921896539e-01,
2142 1.6548602561961109e-01,
2143 1.1393881220195230e-02);
2144 process_point_3(6.1380088247906528e-03,
2145 9.4298876734520487e-01,
2146 3.6194434433925362e-04);
2147 process_point_3(1.2105018114558939e-01,
2148 4.7719037990428043e-01,
2149 2.5739731980456069e-02);
2150 process_point_3(3.2779468216442620e-02,
2151 5.9425626948000698e-01,
2152 1.0135871679755789e-02);
2153 process_point_3(3.2485281564823047e-02,
2154 8.0117728465834437e-01,
2155 6.5761472770359038e-03);
2156 process_point_3(1.7497934218393901e-01,
2157 6.2807184547536599e-01,
2158 1.2907035798861989e-02);
2159 }
2160 break;
2161 default:
2163 }
2164 break;
2165 case 6:
2166 // There is no WV-11 rule in 3d yet
2167 Assert(dim == 2, ExcNotImplemented());
2168 if (use_odd_order)
2169 {
2170 // WV-11, 2d
2171 b_point_permutations.push_back({centroid});
2172 b_weights.push_back(8.5761179732224219e-02);
2173 process_point_1(2.8485417614371900e-02, 1.0431870512894697e-02);
2174 process_point_1(4.9589190096589092e-01, 1.6606273054585369e-02);
2175 process_point_1(1.0263548271224643e-01, 3.8630759237019321e-02);
2176 process_point_1(4.3846592676435220e-01, 6.7316154079468296e-02);
2177 process_point_1(2.1021995670317828e-01, 7.0515684111716576e-02);
2178 process_point_3(7.3254276860644785e-03,
2179 1.4932478865208237e-01,
2180 1.0290289572953278e-02);
2181 process_point_3(4.6010500165429957e-02,
2182 2.8958112563770588e-01,
2183 4.0332476640500554e-02);
2184 }
2185 else
2186 {
2187 // WV-12, 2d
2188 process_point_1(2.4646363436335583e-02, 7.9316425099736389e-03);
2189 process_point_1(4.8820375094554153e-01, 2.4266838081452032e-02);
2190 process_point_1(1.0925782765935427e-01, 2.8486052068877544e-02);
2191 process_point_1(4.4011164865859309e-01, 4.9918334928060942e-02);
2192 process_point_1(2.7146250701492608e-01, 6.2541213195902765e-02);
2193 process_point_3(2.1382490256170616e-02,
2194 1.2727971723358933e-01,
2195 1.5083677576511438e-02);
2196 process_point_3(2.3034156355267121e-02,
2197 2.9165567973834094e-01,
2198 2.1783585038607559e-02);
2199 process_point_3(1.1629601967792658e-01,
2200 2.5545422863851736e-01,
2201 4.3227363659414209e-02);
2202 }
2203 break;
2204 case 7:
2205 // There is no WV-13 rule in 3d yet
2206 Assert(dim == 2, ExcNotImplemented());
2207 if (use_odd_order)
2208 {
2209 // WV-13, 2d
2210 b_point_permutations.push_back({centroid});
2211 b_weights.push_back(6.7960036586831640e-02);
2212 process_point_1(2.1509681108843159e-02, 6.0523371035391717e-03);
2213 process_point_1(4.8907694645253935e-01, 2.3994401928894731e-02);
2214 process_point_1(4.2694141425980042e-01, 5.5601967530453329e-02);
2215 process_point_1(2.2137228629183292e-01, 5.8278485119199981e-02);
2216 process_point_3(5.1263891023823893e-03,
2217 2.7251581777342970e-01,
2218 9.5906810035432631e-03);
2219 process_point_3(2.4370186901093827e-02,
2220 1.1092204280346341e-01,
2221 1.4965401105165668e-02);
2222 process_point_3(8.7895483032197297e-02,
2223 1.6359740106785048e-01,
2224 2.4179039811593819e-02);
2225 process_point_3(6.8012243554206653e-02,
2226 3.0844176089211778e-01,
2227 3.4641276140848373e-02);
2228 }
2229 else
2230 {
2231 // WV-14, 2d
2232 process_point_1(1.9390961248701044e-02, 4.9234036024000819e-03);
2233 process_point_1(6.1799883090872587e-02, 1.4433699669776668e-02);
2234 process_point_1(4.8896391036217862e-01, 2.1883581369428889e-02);
2235 process_point_1(4.1764471934045394e-01, 3.2788353544125348e-02);
2236 process_point_1(1.7720553241254344e-01, 4.2162588736993016e-02);
2237 process_point_1(2.7347752830883865e-01, 5.1774104507291585e-02);
2238 process_point_3(1.2683309328720416e-03,
2239 1.1897449769695684e-01,
2240 5.0102288385006719e-03);
2241 process_point_3(1.4646950055654417e-02,
2242 2.9837288213625779e-01,
2243 1.4436308113533840e-02);
2244 process_point_3(5.7124757403647919e-02,
2245 1.7226668782135557e-01,
2246 2.4665753212563674e-02);
2247 process_point_3(9.2916249356971847e-02,
2248 3.3686145979634496e-01,
2249 3.8571510787060684e-02);
2250 }
2251 break;
2252 default:
2254 }
2255
2256 Assert(b_point_permutations.size() == b_weights.size(), ExcInternalError());
2257 for (unsigned int permutation_n = 0; permutation_n < b_weights.size();
2258 ++permutation_n)
2259 {
2260 for (const std::array<double, dim + 1> &b_point :
2261 b_point_permutations[permutation_n])
2262 {
2263 const double volume = (dim == 2 ? 1.0 / 2.0 : 1.0 / 6.0);
2264 this->weights.emplace_back(volume * b_weights[permutation_n]);
2265 Point<dim> c_point;
2266 for (int d = 0; d < dim; ++d)
2267 c_point[d] = b_point[d];
2268 this->quadrature_points.emplace_back(c_point);
2269 }
2270 }
2271}
2272
2273
2274
2275namespace
2276{
2277 template <int dim>
2279 setup_qiterated_1D(const Quadrature<dim> &, const unsigned int)
2280 {
2282 return Quadrature<dim>();
2283 }
2284
2285
2286
2288 setup_qiterated_1D(const Quadrature<1> &base_quad,
2289 const unsigned int n_copies)
2290 {
2291 return QIterated<1>(base_quad, n_copies);
2292 }
2293} // namespace
2294
2295
2296
2297template <int dim>
2299 const unsigned int n_copies)
2300{
2301 switch (dim)
2302 {
2303 case 1:
2304 static_cast<Quadrature<dim> &>(*this) =
2305 setup_qiterated_1D(base_quad, n_copies);
2306 break;
2307 case 2:
2308 case 3:
2309 {
2310 const auto n_refinements =
2311 static_cast<unsigned int>(std::round(std::log2(n_copies)));
2312 Assert((1u << n_refinements) == n_copies,
2313 ExcMessage("The number of copies must be a power of 2."));
2314 Triangulation<dim> tria;
2315 const auto reference_cell = ReferenceCells::get_simplex<dim>();
2316 GridGenerator::reference_cell(tria, reference_cell);
2317 tria.refine_global(n_refinements);
2318 const Mapping<dim> &mapping =
2319 reference_cell.template get_default_linear_mapping<dim>();
2320 FE_Nothing<dim> fe(reference_cell);
2321
2322 FEValues<dim> fe_values(mapping,
2323 fe,
2324 base_quad,
2326 std::vector<Point<dim>> points;
2327 std::vector<double> weights;
2328 for (const auto &cell : tria.active_cell_iterators())
2329 {
2330 fe_values.reinit(cell);
2331 for (unsigned int qp = 0; qp < base_quad.size(); ++qp)
2332 {
2333 points.push_back(fe_values.quadrature_point(qp));
2334 weights.push_back(fe_values.JxW(qp));
2335 }
2336 }
2337
2338 static_cast<Quadrature<dim> &>(*this) =
2339 Quadrature<dim>(points, weights);
2340
2341 break;
2342 }
2343 default:
2345 }
2346}
2347
2348
2349
2350template <int dim>
2351QGaussWedge<dim>::QGaussWedge(const unsigned int n_points)
2352 : Quadrature<dim>()
2353{
2354 AssertDimension(dim, 3);
2355
2356 const QGaussSimplex<2> quad_tri(n_points);
2357 const QGauss<1> quad_line(n_points);
2358
2359 for (unsigned int i = 0; i < quad_line.size(); ++i)
2360 for (unsigned int j = 0; j < quad_tri.size(); ++j)
2361 {
2362 this->quadrature_points.emplace_back(quad_tri.point(j)[0],
2363 quad_tri.point(j)[1],
2364 quad_line.point(i)[0]);
2365 this->weights.emplace_back(quad_tri.weight(j) * quad_line.weight(i));
2366 }
2367
2368 AssertDimension(this->quadrature_points.size(), this->weights.size());
2369 Assert(this->quadrature_points.size() > 0,
2370 ExcMessage("No valid quadrature points!"));
2371}
2372
2373
2374
2375template <int dim>
2376QGaussPyramid<dim>::QGaussPyramid(const unsigned int n_points_1D)
2377 : Quadrature<dim>()
2378{
2379 AssertDimension(dim, 3);
2380
2381 if (n_points_1D == 1)
2382 {
2383 const double Q14 = 1.0 / 4.0;
2384 const double Q43 = 4.0 / 3.0;
2385
2386 this->quadrature_points.emplace_back(0, 0, Q14);
2387 this->weights.emplace_back(Q43);
2388 }
2389 else if (n_points_1D == 2)
2390 {
2391 // clang-format off
2392 this->quadrature_points.emplace_back(-0.26318405556971, -0.26318405556971, 0.54415184401122);
2393 this->quadrature_points.emplace_back(-0.50661630334979, -0.50661630334979, 0.12251482265544);
2394 this->quadrature_points.emplace_back(-0.26318405556971, +0.26318405556971, 0.54415184401122);
2395 this->quadrature_points.emplace_back(-0.50661630334979, +0.50661630334979, 0.12251482265544);
2396 this->quadrature_points.emplace_back(+0.26318405556971, -0.26318405556971, 0.54415184401122);
2397 this->quadrature_points.emplace_back(+0.50661630334979, -0.50661630334979, 0.12251482265544);
2398 this->quadrature_points.emplace_back(+0.26318405556971, +0.26318405556971, 0.54415184401122);
2399 this->quadrature_points.emplace_back(+0.50661630334979, +0.50661630334979, 0.12251482265544);
2400 // clang-format on
2401
2402 this->weights.emplace_back(0.10078588207983);
2403 this->weights.emplace_back(0.23254745125351);
2404 this->weights.emplace_back(0.10078588207983);
2405 this->weights.emplace_back(0.23254745125351);
2406 this->weights.emplace_back(0.10078588207983);
2407 this->weights.emplace_back(0.23254745125351);
2408 this->weights.emplace_back(0.10078588207983);
2409 this->weights.emplace_back(0.23254745125351);
2410 }
2411
2412 AssertDimension(this->quadrature_points.size(), this->weights.size());
2413 Assert(this->quadrature_points.size() > 0,
2414 ExcMessage("No valid quadrature points!"));
2415}
2416
2417
2418
2419// explicit specialization
2420// note that 1d formulae are specialized by implementation above
2421template class QGauss<2>;
2422template class QGaussRadau<2>;
2423template class QGaussLobatto<2>;
2424template class QMidpoint<2>;
2425template class QTrapezoid<2>;
2426template class QSimpson<2>;
2427template class QMilne<2>;
2428template class QWeddle<2>;
2429
2430template class QGauss<3>;
2431template class QGaussRadau<3>;
2432template class QGaussLobatto<3>;
2433template class QMidpoint<3>;
2434template class QTrapezoid<3>;
2435template class QSimpson<3>;
2436template class QMilne<3>;
2437template class QWeddle<3>;
2438
2439template class QSorted<1>;
2440template class QSorted<2>;
2441template class QSorted<3>;
2442
2443template class QTelles<1>;
2444template class QTelles<2>;
2445template class QTelles<3>;
2446
2447template class QGaussChebyshev<1>;
2448template class QGaussChebyshev<2>;
2449template class QGaussChebyshev<3>;
2450
2451template class QGaussRadauChebyshev<1>;
2452template class QGaussRadauChebyshev<2>;
2453template class QGaussRadauChebyshev<3>;
2454
2455template class QGaussLobattoChebyshev<1>;
2456template class QGaussLobattoChebyshev<2>;
2457template class QGaussLobattoChebyshev<3>;
2458
2459template class QSimplex<1>;
2460template class QSimplex<2>;
2461template class QSimplex<3>;
2462
2463template class QIteratedSimplex<1>;
2464template class QIteratedSimplex<2>;
2465template class QIteratedSimplex<3>;
2466
2467template class QSplit<1>;
2468template class QSplit<2>;
2469template class QSplit<3>;
2470
2471template class QGaussSimplex<0>;
2472template class QGaussSimplex<1>;
2473template class QGaussSimplex<2>;
2474template class QGaussSimplex<3>;
2475template class QGaussWedge<0>;
2476template class QGaussWedge<1>;
2477template class QGaussWedge<2>;
2478template class QGaussWedge<3>;
2479template class QGaussPyramid<0>;
2480template class QGaussPyramid<1>;
2481template class QGaussPyramid<2>;
2482template class QGaussPyramid<3>;
2483
2484template class QWitherdenVincentSimplex<1>;
2485template class QWitherdenVincentSimplex<2>;
2486template class QWitherdenVincentSimplex<3>;
2487
2488#ifndef DOXYGEN
2489template Quadrature<1>
2491 const std::array<Point<1>, 1 + 1> &vertices) const;
2492
2493template Quadrature<2>
2495 const std::array<Point<2>, 1 + 1> &vertices) const;
2496
2497template Quadrature<2>
2499 const std::array<Point<2>, 2 + 1> &vertices) const;
2500
2501template Quadrature<3>
2503 const std::array<Point<3>, 1 + 1> &vertices) const;
2504
2505template Quadrature<3>
2507 const std::array<Point<3>, 2 + 1> &vertices) const;
2508
2509template Quadrature<3>
2511 const std::array<Point<3>, 3 + 1> &vertices) const;
2512
2513template Quadrature<2>
2515 const std::vector<std::array<Point<2>, 1 + 1>> &simplices) const;
2516
2517template Quadrature<3>
2519 const std::vector<std::array<Point<3>, 1 + 1>> &simplices) const;
2520
2521template Quadrature<2>
2523 const std::vector<std::array<Point<2>, 2 + 1>> &simplices) const;
2524
2525template Quadrature<3>
2527 const std::vector<std::array<Point<3>, 2 + 1>> &simplices) const;
2528
2529template Quadrature<3>
2531 const std::vector<std::array<Point<3>, 3 + 1>> &simplices) const;
2532#endif
2533
DerivativeForm< 1, spacedim, dim, Number > transpose() const
const Point< spacedim > & quadrature_point(const unsigned int q_point) const
double JxW(const unsigned int q_point) const
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell)
Abstract base class for mapping classes.
Definition mapping.h:318
Definition point.h:111
QDuffy(const Quadrature< 1 > &radial_quadrature, const Quadrature< 1 > &angular_quadrature, const double beta=1.0)
QGaussChebyshev(const unsigned int n)
Generate a formula with n quadrature points.
QGaussLobattoChebyshev(const unsigned int n)
Generate a formula with n quadrature points.
QGaussLobatto(const unsigned int n)
const double fraction
QGaussLogR(const unsigned int n, const Point< dim > &x0=Point< dim >(), const double alpha=1, const bool factor_out_singular_weight=false)
QGaussLog(const unsigned int n, const bool revert=false)
static std::vector< double > get_quadrature_points(const unsigned int n)
static std::vector< double > get_quadrature_weights(const unsigned int n)
QGaussOneOverR(const unsigned int n, const Point< dim > &singularity, const bool factor_out_singular_weight=false)
static unsigned int quad_size(const Point< dim > &singularity, const unsigned int n)
QGaussPyramid(const unsigned int n_points_1D)
QGaussRadauChebyshev(const unsigned int n, const EndPoint end_point=QGaussRadauChebyshev::EndPoint::left)
Generate a formula with n quadrature points.
const EndPoint end_point
QGaussRadau(const unsigned int n, const EndPoint end_point=QGaussRadau::EndPoint::left)
QGaussSimplex(const unsigned int n_points_1D)
QGaussWedge(const unsigned int n_points_1D)
QGauss(const unsigned int n)
QIteratedSimplex(const Quadrature< dim > &base_quadrature, const unsigned int n_copies)
Quadrature< spacedim > compute_affine_transformation(const std::array< Point< spacedim >, dim+1 > &vertices) const
QSimplex(const Quadrature< dim > &quad)
Quadrature< spacedim > mapped_quadrature(const std::vector< std::array< Point< spacedim >, dim+1 > > &simplices) const
QSorted(const Quadrature< dim > &quad)
bool compare_weights(const unsigned int a, const unsigned int b) const
QSplit(const QSimplex< dim > &base, const Point< dim > &split_point)
QTelles(const Quadrature< 1 > &base_quad, const Point< dim > &singularity)
QTrianglePolar(const Quadrature< 1 > &radial_quadrature, const Quadrature< 1 > &angular_quadrature)
QWitherdenVincentSimplex(const unsigned int n_points_1D, const bool use_odd_order=true)
Quadrature & operator=(const Quadrature< dim > &)
const Point< dim > & point(const unsigned int i) const
double weight(const unsigned int i) const
unsigned int size() const
void refine_global(const unsigned int times=1)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:498
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:499
Tensor< 1, spacedim, typename ProductType< Number1, Number2 >::type > apply_transformation(const DerivativeForm< 1, dim, spacedim, Number1 > &grad_F, const Tensor< 1, dim, Number2 > &d_x)
IteratorRange< active_cell_iterator > active_cell_iterators() const
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
@ update_JxW_values
Transformed quadrature weights.
@ update_quadrature_points
Transformed quadrature points.
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
std::size_t size
Definition mpi.cc:734
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
void quadrature_points(const Triangulation< dim, spacedim > &triangulation, const Quadrature< dim > &quadrature, const std::vector< std::vector< BoundingBox< spacedim > > > &global_bounding_boxes, ParticleHandler< dim, spacedim > &particle_handler, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()), const std::vector< std::vector< double > > &properties={})
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Number jacobi_polynomial_value(const unsigned int degree, const int alpha, const int beta, const Number x)
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:479
constexpr T pow(const T base, const int iexp)
Definition utilities.h:966
std::vector< double > get_quadrature_weights(const unsigned int n)
std::vector< double > get_quadrature_points(const unsigned int n)
std::vector< double > get_quadrature_points(const unsigned int n)
std::vector< double > get_quadrature_weights(const unsigned int n)
std::vector< long double > compute_quadrature_weights(const std::vector< long double > &x, const int alpha, const int beta)
long double gamma(const unsigned int n)
std::vector< double > get_quadrature_weights(const unsigned int n, const ::QGaussRadauChebyshev< 1 >::EndPoint end_point)
std::vector< double > get_quadrature_points(const unsigned int n, const ::QGaussRadauChebyshev< 1 >::EndPoint end_point)
std::vector< double > get_left_quadrature_weights(const unsigned int n)
std::vector< double > get_left_quadrature_points(const unsigned int n)
std::vector< double > get_quadrature_weights(const unsigned int n, const ::QGaussRadau< 1 >::EndPoint end_point)
std::vector< double > get_quadrature_points(const unsigned int n, const ::QGaussRadau< 1 >::EndPoint end_point)
static constexpr double PI_2
Definition numbers.h:259
static constexpr double PI
Definition numbers.h:254
::VectorizedArray< Number, width > log(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > tan(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices()
static Point< dim > unit_cell_vertex(const unsigned int vertex)