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