Reference documentation for deal.II version 9.4.0
\(\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\}}\)
quadrature_lib.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 1998 - 2021 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
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 std::pow(2., alpha + beta + 1) * gamma(alpha + q) * gamma(beta + q) /
130 ((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
146#ifndef DOXYGEN
147template <>
148QGaussLobatto<1>::QGaussLobatto(const unsigned int n)
149 : Quadrature<1>(n)
150{
151 Assert(n >= 2, ExcNotImplemented());
152
153 std::vector<long double> points =
154 Polynomials::jacobi_polynomial_roots<long double>(n - 2, 1, 1);
155 points.insert(points.begin(), 0);
156 points.push_back(1.);
157 std::vector<long double> w =
159
160 // scale weights to the interval [0.0, 1.0]:
161 for (unsigned int i = 0; i < points.size(); ++i)
162 {
163 this->quadrature_points[i][0] = points[i];
164 this->weights[i] = 0.5 * w[i];
165 }
166}
167#endif
168
169
170template <>
172 : Quadrature<1>(1)
173{
174 this->quadrature_points[0] = Point<1>(0.5);
175 this->weights[0] = 1.0;
176}
177
178
179
180template <>
182 : Quadrature<1>(2)
183{
184 static const double xpts[] = {0.0, 1.0};
185 static const double wts[] = {0.5, 0.5};
186
187 for (unsigned int i = 0; i < this->size(); ++i)
188 {
189 this->quadrature_points[i] = Point<1>(xpts[i]);
190 this->weights[i] = wts[i];
191 }
192}
193
194
195
196template <>
198 : Quadrature<1>(3)
199{
200 static const double xpts[] = {0.0, 0.5, 1.0};
201 static const double wts[] = {1. / 6., 2. / 3., 1. / 6.};
202
203 for (unsigned int i = 0; i < this->size(); ++i)
204 {
205 this->quadrature_points[i] = Point<1>(xpts[i]);
206 this->weights[i] = wts[i];
207 }
208}
209
210
211
212template <>
214 : Quadrature<1>(5)
215{
216 static const double xpts[] = {0.0, .25, .5, .75, 1.0};
217 static const double wts[] = {
218 7. / 90., 32. / 90., 12. / 90., 32. / 90., 7. / 90.};
219
220 for (unsigned int i = 0; i < this->size(); ++i)
221 {
222 this->quadrature_points[i] = Point<1>(xpts[i]);
223 this->weights[i] = wts[i];
224 }
225}
226
227
228
229template <>
231 : Quadrature<1>(7)
232{
233 static const double xpts[] = {
234 0.0, 1. / 6., 1. / 3., .5, 2. / 3., 5. / 6., 1.0};
235 static const double wts[] = {41. / 840.,
236 216. / 840.,
237 27. / 840.,
238 272. / 840.,
239 27. / 840.,
240 216. / 840.,
241 41. / 840.};
242
243 for (unsigned int i = 0; i < this->size(); ++i)
244 {
245 this->quadrature_points[i] = Point<1>(xpts[i]);
246 this->weights[i] = wts[i];
247 }
248}
249
250
251template <>
252QGaussLog<1>::QGaussLog(const unsigned int n, const bool revert)
253 : Quadrature<1>(n)
254{
255 std::vector<double> p = get_quadrature_points(n);
256 std::vector<double> w = get_quadrature_weights(n);
257
258 for (unsigned int i = 0; i < this->size(); ++i)
259 {
260 // Using the change of variables x=1-t, it's possible to show
261 // that int f(x)ln|1-x| = int f(1-t) ln|t|, which implies that
262 // we can use this quadrature formula also with weight ln|1-x|.
263 this->quadrature_points[i] =
264 revert ? Point<1>(1 - p[n - 1 - i]) : Point<1>(p[i]);
265 this->weights[i] = revert ? w[n - 1 - i] : w[i];
266 }
267}
268
269template <>
270std::vector<double>
272{
273 std::vector<double> q_points(n);
274
275 switch (n)
276 {
277 case 1:
278 q_points[0] = 0.3333333333333333;
279 break;
280
281 case 2:
282 q_points[0] = 0.1120088061669761;
283 q_points[1] = 0.6022769081187381;
284 break;
285
286 case 3:
287 q_points[0] = 0.06389079308732544;
288 q_points[1] = 0.3689970637156184;
289 q_points[2] = 0.766880303938942;
290 break;
291
292 case 4:
293 q_points[0] = 0.04144848019938324;
294 q_points[1] = 0.2452749143206022;
295 q_points[2] = 0.5561654535602751;
296 q_points[3] = 0.848982394532986;
297 break;
298
299 case 5:
300 q_points[0] = 0.02913447215197205;
301 q_points[1] = 0.1739772133208974;
302 q_points[2] = 0.4117025202849029;
303 q_points[3] = 0.6773141745828183;
304 q_points[4] = 0.89477136103101;
305 break;
306
307 case 6:
308 q_points[0] = 0.02163400584411693;
309 q_points[1] = 0.1295833911549506;
310 q_points[2] = 0.3140204499147661;
311 q_points[3] = 0.5386572173517997;
312 q_points[4] = 0.7569153373774084;
313 q_points[5] = 0.922668851372116;
314 break;
315
316
317 case 7:
318 q_points[0] = 0.0167193554082585;
319 q_points[1] = 0.100185677915675;
320 q_points[2] = 0.2462942462079286;
321 q_points[3] = 0.4334634932570557;
322 q_points[4] = 0.6323509880476823;
323 q_points[5] = 0.81111862674023;
324 q_points[6] = 0.940848166743287;
325 break;
326
327 case 8:
328 q_points[0] = 0.01332024416089244;
329 q_points[1] = 0.07975042901389491;
330 q_points[2] = 0.1978710293261864;
331 q_points[3] = 0.354153994351925;
332 q_points[4] = 0.5294585752348643;
333 q_points[5] = 0.7018145299391673;
334 q_points[6] = 0.849379320441094;
335 q_points[7] = 0.953326450056343;
336 break;
337
338 case 9:
339 q_points[0] = 0.01086933608417545;
340 q_points[1] = 0.06498366633800794;
341 q_points[2] = 0.1622293980238825;
342 q_points[3] = 0.2937499039716641;
343 q_points[4] = 0.4466318819056009;
344 q_points[5] = 0.6054816627755208;
345 q_points[6] = 0.7541101371585467;
346 q_points[7] = 0.877265828834263;
347 q_points[8] = 0.96225055941096;
348 break;
349
350 case 10:
351 q_points[0] = 0.00904263096219963;
352 q_points[1] = 0.05397126622250072;
353 q_points[2] = 0.1353118246392511;
354 q_points[3] = 0.2470524162871565;
355 q_points[4] = 0.3802125396092744;
356 q_points[5] = 0.5237923179723384;
357 q_points[6] = 0.6657752055148032;
358 q_points[7] = 0.7941904160147613;
359 q_points[8] = 0.898161091216429;
360 q_points[9] = 0.9688479887196;
361 break;
362
363
364 case 11:
365 q_points[0] = 0.007643941174637681;
366 q_points[1] = 0.04554182825657903;
367 q_points[2] = 0.1145222974551244;
368 q_points[3] = 0.2103785812270227;
369 q_points[4] = 0.3266955532217897;
370 q_points[5] = 0.4554532469286375;
371 q_points[6] = 0.5876483563573721;
372 q_points[7] = 0.7139638500230458;
373 q_points[8] = 0.825453217777127;
374 q_points[9] = 0.914193921640008;
375 q_points[10] = 0.973860256264123;
376 break;
377
378 case 12:
379 q_points[0] = 0.006548722279080035;
380 q_points[1] = 0.03894680956045022;
381 q_points[2] = 0.0981502631060046;
382 q_points[3] = 0.1811385815906331;
383 q_points[4] = 0.2832200676673157;
384 q_points[5] = 0.398434435164983;
385 q_points[6] = 0.5199526267791299;
386 q_points[7] = 0.6405109167754819;
387 q_points[8] = 0.7528650118926111;
388 q_points[9] = 0.850240024421055;
389 q_points[10] = 0.926749682988251;
390 q_points[11] = 0.977756129778486;
391 break;
392
393 default:
394 Assert(false, ExcNotImplemented());
395 break;
396 }
397
398 return q_points;
399}
400
401
402template <>
403std::vector<double>
405{
406 std::vector<double> quadrature_weights(n);
407
408 switch (n)
409 {
410 case 1:
411 quadrature_weights[0] = -1.0;
412 break;
413 case 2:
414 quadrature_weights[0] = -0.7185393190303845;
415 quadrature_weights[1] = -0.2814606809696154;
416 break;
417
418 case 3:
419 quadrature_weights[0] = -0.5134045522323634;
420 quadrature_weights[1] = -0.3919800412014877;
421 quadrature_weights[2] = -0.0946154065661483;
422 break;
423
424 case 4:
425 quadrature_weights[0] = -0.3834640681451353;
426 quadrature_weights[1] = -0.3868753177747627;
427 quadrature_weights[2] = -0.1904351269501432;
428 quadrature_weights[3] = -0.03922548712995894;
429 break;
430
431 case 5:
432 quadrature_weights[0] = -0.2978934717828955;
433 quadrature_weights[1] = -0.3497762265132236;
434 quadrature_weights[2] = -0.234488290044052;
435 quadrature_weights[3] = -0.0989304595166356;
436 quadrature_weights[4] = -0.01891155214319462;
437 break;
438
439 case 6:
440 quadrature_weights[0] = -0.2387636625785478;
441 quadrature_weights[1] = -0.3082865732739458;
442 quadrature_weights[2] = -0.2453174265632108;
443 quadrature_weights[3] = -0.1420087565664786;
444 quadrature_weights[4] = -0.05545462232488041;
445 quadrature_weights[5] = -0.01016895869293513;
446 break;
447
448
449 case 7:
450 quadrature_weights[0] = -0.1961693894252476;
451 quadrature_weights[1] = -0.2703026442472726;
452 quadrature_weights[2] = -0.239681873007687;
453 quadrature_weights[3] = -0.1657757748104267;
454 quadrature_weights[4] = -0.0889432271377365;
455 quadrature_weights[5] = -0.03319430435645653;
456 quadrature_weights[6] = -0.005932787015162054;
457 break;
458
459 case 8:
460 quadrature_weights[0] = -0.164416604728002;
461 quadrature_weights[1] = -0.2375256100233057;
462 quadrature_weights[2] = -0.2268419844319134;
463 quadrature_weights[3] = -0.1757540790060772;
464 quadrature_weights[4] = -0.1129240302467932;
465 quadrature_weights[5] = -0.05787221071771947;
466 quadrature_weights[6] = -0.02097907374214317;
467 quadrature_weights[7] = -0.003686407104036044;
468 break;
469
470 case 9:
471 quadrature_weights[0] = -0.1400684387481339;
472 quadrature_weights[1] = -0.2097722052010308;
473 quadrature_weights[2] = -0.211427149896601;
474 quadrature_weights[3] = -0.1771562339380667;
475 quadrature_weights[4] = -0.1277992280331758;
476 quadrature_weights[5] = -0.07847890261203835;
477 quadrature_weights[6] = -0.0390225049841783;
478 quadrature_weights[7] = -0.01386729555074604;
479 quadrature_weights[8] = -0.002408041036090773;
480 break;
481
482 case 10:
483 quadrature_weights[0] = -0.12095513195457;
484 quadrature_weights[1] = -0.1863635425640733;
485 quadrature_weights[2] = -0.1956608732777627;
486 quadrature_weights[3] = -0.1735771421828997;
487 quadrature_weights[4] = -0.135695672995467;
488 quadrature_weights[5] = -0.0936467585378491;
489 quadrature_weights[6] = -0.05578772735275126;
490 quadrature_weights[7] = -0.02715981089692378;
491 quadrature_weights[8] = -0.00951518260454442;
492 quadrature_weights[9] = -0.001638157633217673;
493 break;
494
495
496 case 11:
497 quadrature_weights[0] = -0.1056522560990997;
498 quadrature_weights[1] = -0.1665716806006314;
499 quadrature_weights[2] = -0.1805632182877528;
500 quadrature_weights[3] = -0.1672787367737502;
501 quadrature_weights[4] = -0.1386970574017174;
502 quadrature_weights[5] = -0.1038334333650771;
503 quadrature_weights[6] = -0.06953669788988512;
504 quadrature_weights[7] = -0.04054160079499477;
505 quadrature_weights[8] = -0.01943540249522013;
506 quadrature_weights[9] = -0.006737429326043388;
507 quadrature_weights[10] = -0.001152486965101561;
508 break;
509
510 case 12:
511 quadrature_weights[0] = -0.09319269144393;
512 quadrature_weights[1] = -0.1497518275763289;
513 quadrature_weights[2] = -0.166557454364573;
514 quadrature_weights[3] = -0.1596335594369941;
515 quadrature_weights[4] = -0.1384248318647479;
516 quadrature_weights[5] = -0.1100165706360573;
517 quadrature_weights[6] = -0.07996182177673273;
518 quadrature_weights[7] = -0.0524069547809709;
519 quadrature_weights[8] = -0.03007108900074863;
520 quadrature_weights[9] = -0.01424924540252916;
521 quadrature_weights[10] = -0.004899924710875609;
522 quadrature_weights[11] = -0.000834029009809656;
523 break;
524
525 default:
526 Assert(false, ExcNotImplemented());
527 break;
528 }
529
530 return quadrature_weights;
531}
532
533
534template <>
535QGaussLogR<1>::QGaussLogR(const unsigned int n,
536 const Point<1> & origin,
537 const double alpha,
538 const bool factor_out_singularity)
539 : Quadrature<1>(
540 ((origin[0] == 0) || (origin[0] == 1)) ? (alpha == 1 ? n : 2 * n) : 4 * n)
541 , fraction(((origin[0] == 0) || (origin[0] == 1.)) ? 1. : origin[0])
542{
543 // The three quadrature formulas that make this one up. There are
544 // at most two when the origin is one of the extremes, and there is
545 // only one if the origin is one of the extremes and alpha is
546 // equal to one.
547 //
548 // If alpha is different from one, then we need a correction which
549 // is performed with a standard Gauss quadrature rule on each
550 // segment. This is not needed in the standard case where alpha is
551 // equal to one and the origin is on one of the extremes. We
552 // integrate with weight ln|(x-o)/alpha|. In the easy cases, we
553 // only need n quadrature points. In the most difficult one, we
554 // need 2*n points for the first segment, and 2*n points for the
555 // second segment.
556 QGaussLog<1> quad1(n, origin[0] != 0);
557 QGaussLog<1> quad2(n);
558 QGauss<1> quad(n);
559
560 // Check that the origin is inside 0,1
561 Assert((fraction >= 0) && (fraction <= 1),
562 ExcMessage("Origin is outside [0,1]."));
563
564 // Non singular offset. This is the start of non singular quad
565 // points.
566 unsigned int ns_offset = (fraction == 1) ? n : 2 * n;
567
568 for (unsigned int i = 0, j = ns_offset; i < n; ++i, ++j)
569 {
570 // The first i quadrature points are the same as quad1, and
571 // are by default singular.
572 this->quadrature_points[i] = quad1.point(i) * fraction;
573 this->weights[i] = quad1.weight(i) * fraction;
574
575 // We need to scale with -log|fraction*alpha|
576 if ((alpha != 1) || (fraction != 1))
577 {
578 this->quadrature_points[j] = quad.point(i) * fraction;
579 this->weights[j] =
580 -std::log(alpha / fraction) * quad.weight(i) * fraction;
581 }
582 // In case we need the second quadrature as well, do it now.
583 if (fraction != 1)
584 {
585 this->quadrature_points[i + n] =
586 quad2.point(i) * (1 - fraction) + Point<1>(fraction);
587 this->weights[i + n] = quad2.weight(i) * (1 - fraction);
588
589 // We need to scale with -log|fraction*alpha|
590 this->quadrature_points[j + n] =
591 quad.point(i) * (1 - fraction) + Point<1>(fraction);
592 this->weights[j + n] =
593 -std::log(alpha / (1 - fraction)) * quad.weight(i) * (1 - fraction);
594 }
595 }
596 if (factor_out_singularity == true)
597 for (unsigned int i = 0; i < size(); ++i)
598 {
599 Assert(
600 this->quadrature_points[i] != origin,
602 "The singularity cannot be on a Gauss point of the same order!"));
603 double denominator =
604 std::log(std::abs((this->quadrature_points[i] - origin)[0]) / alpha);
605 Assert(denominator != 0.0,
607 "The quadrature formula you are using does not allow to "
608 "factor out the singularity, which is zero at one point."));
609 this->weights[i] /= denominator;
610 }
611}
612
613
614template <>
615unsigned int
616QGaussOneOverR<2>::quad_size(const Point<2> &singularity, const unsigned int n)
617{
618 const double eps = 1e-8;
619 const bool on_edge =
620 std::any_of(singularity.begin_raw(),
621 singularity.end_raw(),
622 [eps](double coord) {
623 return std::abs(coord) < eps || std::abs(coord - 1.) < eps;
624 });
625 const bool on_vertex =
626 on_edge &&
627 std::abs((singularity - Point<2>(.5, .5)).norm_square() - .5) < eps;
628 if (on_vertex)
629 return 2 * n * n;
630 else if (on_edge)
631 return 4 * n * n;
632 else
633 return 8 * n * n;
634}
635
636template <>
638 const Point<2> & singularity,
639 const bool factor_out_singularity)
640 : Quadrature<2>(quad_size(singularity, n))
641{
642 // We treat all the cases in the
643 // same way. Split the element in 4
644 // pieces, measure the area, if
645 // it's relevant, add the
646 // quadrature connected to that
647 // singularity.
648 std::vector<QGaussOneOverR<2>> quads;
649 std::vector<Point<2>> origins;
650 // Id of the corner with a
651 // singularity
652 quads.emplace_back(n, 3, factor_out_singularity);
653 quads.emplace_back(n, 2, factor_out_singularity);
654 quads.emplace_back(n, 1, factor_out_singularity);
655 quads.emplace_back(n, 0, factor_out_singularity);
656
657 origins.emplace_back(0., 0.);
658 origins.emplace_back(singularity[0], 0.);
659 origins.emplace_back(0., singularity[1]);
660 origins.push_back(singularity);
661
662 // Lexicographical ordering.
663
664 double eps = 1e-8;
665 unsigned int q_id = 0; // Current quad point index.
666 Tensor<1, 2> dist;
667
668 for (unsigned int box = 0; box < 4; ++box)
669 {
670 dist = (singularity - GeometryInfo<2>::unit_cell_vertex(box));
671 dist = Point<2>(std::abs(dist[0]), std::abs(dist[1]));
672 double area = dist[0] * dist[1];
673 if (area > eps)
674 for (unsigned int q = 0; q < quads[box].size(); ++q, ++q_id)
675 {
676 const Point<2> &qp = quads[box].point(q);
677 this->quadrature_points[q_id] =
678 origins[box] + Point<2>(dist[0] * qp[0], dist[1] * qp[1]);
679 this->weights[q_id] = quads[box].weight(q) * area;
680 }
681 }
682}
683
684
685template <>
687 const unsigned int vertex_index,
688 const bool factor_out_singularity)
689 : Quadrature<2>(2 * n * n)
690{
691 // This version of the constructor
692 // works only for the 4
693 // vertices. If you need a more
694 // general one, you should use the
695 // one with the Point<2> in the
696 // constructor.
697 AssertIndexRange(vertex_index, 4);
698
699 // Start with the gauss quadrature formula on the (u,v) reference
700 // element.
701 QGauss<2> gauss(n);
702
703 Assert(gauss.size() == n * n, ExcInternalError());
704
705 // For the moment we only implemented this for the vertices of a
706 // quadrilateral. We are planning to do this also for the support
707 // points of arbitrary FE_Q elements, to allow the use of this
708 // class in boundary element programs with higher order mappings.
709 AssertIndexRange(vertex_index, 4);
710
711 // We create only the first one. All other pieces are rotation of
712 // this one.
713 // In this case the transformation is
714 //
715 // (x,y) = (u, u tan(pi/4 v))
716 //
717 // with Jacobian
718 //
719 // J = pi/4 R / cos(pi/4 v)
720 //
721 // And we get rid of R to take into account the singularity,
722 // unless specified differently in the constructor.
723 std::vector<Point<2>> &ps = this->quadrature_points;
724 std::vector<double> & ws = this->weights;
725 double pi4 = numbers::PI / 4;
726
727 for (unsigned int q = 0; q < gauss.size(); ++q)
728 {
729 const Point<2> &gp = gauss.point(q);
730 ps[q][0] = gp[0];
731 ps[q][1] = gp[0] * std::tan(pi4 * gp[1]);
732 ws[q] = gauss.weight(q) * pi4 / std::cos(pi4 * gp[1]);
733 if (factor_out_singularity)
734 ws[q] *= (ps[q] - GeometryInfo<2>::unit_cell_vertex(0)).norm();
735 // The other half of the quadrilateral is symmetric with
736 // respect to xy plane.
737 ws[gauss.size() + q] = ws[q];
738 ps[gauss.size() + q][0] = ps[q][1];
739 ps[gauss.size() + q][1] = ps[q][0];
740 }
741
742 // Now we distribute these vertices in the correct manner
743 double theta = 0;
744 switch (vertex_index)
745 {
746 case 0:
747 theta = 0;
748 break;
749 case 1:
750 //
751 theta = numbers::PI / 2;
752 break;
753 case 2:
754 theta = -numbers::PI / 2;
755 break;
756 case 3:
757 theta = numbers::PI;
758 break;
759 }
760
761 double R00 = std::cos(theta), R01 = -std::sin(theta);
762 double R10 = std::sin(theta), R11 = std::cos(theta);
763
764 if (vertex_index != 0)
765 for (unsigned int q = 0; q < size(); ++q)
766 {
767 double x = ps[q][0] - .5, y = ps[q][1] - .5;
768
769 ps[q][0] = R00 * x + R01 * y + .5;
770 ps[q][1] = R10 * x + R11 * y + .5;
771 }
772}
773
774
775template <int dim>
777 : Quadrature<dim>(quad)
778{
779 std::vector<unsigned int> permutation(quad.size());
780 for (unsigned int i = 0; i < quad.size(); ++i)
781 permutation[i] = i;
782
783 std::sort(permutation.begin(),
784 permutation.end(),
785 [this](const unsigned int x, const unsigned int y) {
786 return this->compare_weights(x, y);
787 });
788
789 // At this point, the variable is_tensor_product_flag is set
790 // to the respective value of the given Quadrature in the base
791 // class copy constructor.
792 // We only call a quadrature formula 'tensor product'
793 // if the quadrature points are also sorted lexicographically.
794 // In particular, any reordering destroys that property
795 // and we might need to modify the variable accordingly.
796 for (unsigned int i = 0; i < quad.size(); ++i)
797 {
798 this->weights[i] = quad.weight(permutation[i]);
799 this->quadrature_points[i] = quad.point(permutation[i]);
800 if (permutation[i] != i)
801 this->is_tensor_product_flag = false;
802 }
803}
804
805
806template <int dim>
807bool
808QSorted<dim>::compare_weights(const unsigned int a, const unsigned int b) const
809{
810 return (this->weights[a] < this->weights[b]);
811}
812
813
814// construct the quadrature formulae in higher dimensions by
815// tensor product of lower dimensions
816
817template <int dim>
818QGauss<dim>::QGauss(const unsigned int n)
819 : Quadrature<dim>(QGauss<dim - 1>(n), QGauss<1>(n))
820{}
821
822
823
824template <int dim>
826 : Quadrature<dim>(QGaussLobatto<dim - 1>(n), QGaussLobatto<1>(n))
827{}
828
829
830
831template <int dim>
833 : Quadrature<dim>(QMidpoint<dim - 1>(), QMidpoint<1>())
834{}
835
836
837
838template <int dim>
840 : Quadrature<dim>(QTrapezoid<dim - 1>(), QTrapezoid<1>())
841{}
842
843
844
845template <int dim>
847 : Quadrature<dim>(QSimpson<dim - 1>(), QSimpson<1>())
848{}
849
850
851
852template <int dim>
854 : Quadrature<dim>(QMilne<dim - 1>(), QMilne<1>())
855{}
856
857
858template <int dim>
860 : Quadrature<dim>(QWeddle<dim - 1>(), QWeddle<1>())
861{}
862
863template <int dim>
865 const Point<dim> & singularity)
866 : // We need the explicit implementation if dim == 1. If dim > 1 we use the
867 // former implementation and apply a tensorial product to obtain the higher
868 // dimensions.
869 Quadrature<dim>(
870 dim == 2 ?
871 QAnisotropic<dim>(QTelles<1>(base_quad, Point<1>(singularity[0])),
872 QTelles<1>(base_quad, Point<1>(singularity[1]))) :
873 dim == 3 ?
874 QAnisotropic<dim>(QTelles<1>(base_quad, Point<1>(singularity[0])),
875 QTelles<1>(base_quad, Point<1>(singularity[1])),
876 QTelles<1>(base_quad, Point<1>(singularity[2]))) :
877 Quadrature<dim>())
878{}
879
880template <int dim>
881QTelles<dim>::QTelles(const unsigned int n, const Point<dim> &singularity)
882 : // In this case we map the standard Gauss Legendre formula using the given
883 // singularity point coordinates.
884 Quadrature<dim>(QTelles<dim>(QGauss<1>(n), singularity))
885{}
886
887
888
889template <>
890QTelles<1>::QTelles(const Quadrature<1> &base_quad, const Point<1> &singularity)
891 : // We explicitly implement the Telles' variable change if dim == 1.
892 Quadrature<1>(base_quad)
893{
894 // We define all the constants to be used in the implementation of
895 // Telles' rule
896 const double eta_bar = singularity[0] * 2. - 1.;
897 const double eta_star = eta_bar * eta_bar - 1.;
898 double gamma_bar;
899
900 std::vector<Point<1>> quadrature_points_dummy(quadrature_points.size());
901 std::vector<double> weights_dummy(weights.size());
902 unsigned int cont = 0;
903 const double tol = 1e-10;
904 for (unsigned int d = 0; d < quadrature_points.size(); ++d)
905 {
906 if (std::abs(quadrature_points[d][0] - singularity[0]) > tol)
907 {
908 quadrature_points_dummy[d - cont] = quadrature_points[d];
909 weights_dummy[d - cont] = weights[d];
910 }
911 else
912 {
913 // We need to remove the singularity point from the quadrature point
914 // list. To do so we use the variable cont.
915 cont = 1;
916 }
917 }
918 if (cont == 1)
919 {
920 quadrature_points.resize(quadrature_points_dummy.size() - 1);
921 weights.resize(weights_dummy.size() - 1);
922 for (unsigned int d = 0; d < quadrature_points.size(); ++d)
923 {
924 quadrature_points[d] = quadrature_points_dummy[d];
925 weights[d] = weights_dummy[d];
926 }
927 }
928 // We need to check if the singularity is at the boundary of the interval.
929 if (std::abs(eta_star) <= tol)
930 {
931 gamma_bar =
932 std::pow((eta_bar * eta_star + std::abs(eta_star)), 1.0 / 3.0) +
933 std::pow((eta_bar * eta_star - std::abs(eta_star)), 1.0 / 3.0) +
934 eta_bar;
935 }
936 else
937 {
938 gamma_bar = (eta_bar * eta_star + std::abs(eta_star)) /
939 std::abs(eta_bar * eta_star + std::abs(eta_star)) *
940 std::pow(std::abs(eta_bar * eta_star + std::abs(eta_star)),
941 1.0 / 3.0) +
942 (eta_bar * eta_star - std::abs(eta_star)) /
943 std::abs(eta_bar * eta_star - std::abs(eta_star)) *
944 std::pow(std::abs(eta_bar * eta_star - std::abs(eta_star)),
945 1.0 / 3.0) +
946 eta_bar;
947 }
948 for (unsigned int q = 0; q < quadrature_points.size(); ++q)
949 {
950 double gamma = quadrature_points[q][0] * 2 - 1;
951 double eta = (std::pow(gamma - gamma_bar, 3.0) +
952 gamma_bar * (gamma_bar * gamma_bar + 3)) /
953 (1 + 3 * gamma_bar * gamma_bar);
954
955 double J = 3 * ((gamma - gamma_bar) * (gamma - gamma_bar)) /
956 (1 + 3 * gamma_bar * gamma_bar);
957
958 quadrature_points[q][0] = (eta + 1) / 2.0;
959 weights[q] = J * weights[q];
960 }
961}
962
963namespace internal
964{
966 {
970 std::vector<double>
971 get_quadrature_points(const unsigned int n)
972 {
973 std::vector<double> points(n);
974 // n point quadrature: index from 0 to n-1
975 for (unsigned short i = 0; i < n; ++i)
976 // would be cos((2i+1)Pi/(2N+2))
977 // put + Pi so we start from the smallest point
978 // then map from [-1,1] to [0,1]
979 points[i] =
980 1. / 2. *
981 (1. + std::cos(numbers::PI *
982 (1. + double(2 * i + 1) / double(2 * (n - 1) + 2))));
983
984 return points;
985 }
986
987
988
992 std::vector<double>
993 get_quadrature_weights(const unsigned int n)
994 {
995 std::vector<double> weights(n);
996
997 for (unsigned short i = 0; i < n; ++i)
998 {
999 // same weights as on [-1,1]
1000 weights[i] = numbers::PI / double(n);
1001 }
1002
1003 return weights;
1004 }
1005 } // namespace QGaussChebyshev
1006} // namespace internal
1007
1008
1009template <>
1011 : Quadrature<1>(n)
1012{
1013 Assert(n > 0, ExcMessage("Need at least one point for the quadrature rule"));
1014 std::vector<double> p = internal::QGaussChebyshev::get_quadrature_points(n);
1016
1017 for (unsigned int i = 0; i < this->size(); ++i)
1018 {
1019 this->quadrature_points[i] = Point<1>(p[i]);
1020 this->weights[i] = w[i];
1021 }
1022}
1023
1024
1025template <int dim>
1027 : Quadrature<dim>(QGaussChebyshev<1>(n))
1028{}
1029
1030
1031namespace internal
1032{
1034 {
1035 // Computes the points of the quadrature formula.
1036 std::vector<double>
1037 get_quadrature_points(const unsigned int n,
1039 {
1040 std::vector<double> points(n);
1041 // n point quadrature: index from 0 to n-1
1042 for (unsigned short i = 0; i < n; ++i)
1043 // would be -cos(2i Pi/(2N+1))
1044 // put + Pi so we start from the smallest point
1045 // then map from [-1,1] to [0,1]
1046 switch (ep)
1047 {
1048 case ::QGaussRadauChebyshev<1>::left:
1049 {
1050 points[i] =
1051 1. / 2. *
1052 (1. -
1054 (1 + 2 * double(i) / (2 * double(n - 1) + 1.))));
1055 break;
1056 }
1057
1058 case ::QGaussRadauChebyshev<1>::right:
1059 {
1060 points[i] =
1061 1. / 2. *
1062 (1. - std::cos(numbers::PI * (2 * double(n - 1 - i) /
1063 (2 * double(n - 1) + 1.))));
1064 break;
1065 }
1066
1067 default:
1068 Assert(
1069 false,
1070 ExcMessage(
1071 "This constructor can only be called with either "
1072 "QGaussRadauChebyshev::left or QGaussRadauChebyshev::right as "
1073 "second argument."));
1074 }
1075
1076 return points;
1077 }
1078
1079
1080
1081 // Computes the weights of the quadrature formula.
1082 std::vector<double>
1083 get_quadrature_weights(const unsigned int n,
1085 {
1086 std::vector<double> weights(n);
1087
1088 for (unsigned short i = 0; i < n; ++i)
1089 {
1090 // same weights as on [-1,1]
1091 weights[i] = 2. * numbers::PI / double(2 * (n - 1) + 1.);
1092 if (ep == ::QGaussRadauChebyshev<1>::left && i == 0)
1093 weights[i] /= 2.;
1094 else if (ep == ::QGaussRadauChebyshev<1>::right &&
1095 i == (n - 1))
1096 weights[i] /= 2.;
1097 }
1098
1099 return weights;
1100 }
1101 } // namespace QGaussRadauChebyshev
1102} // namespace internal
1103
1104
1105template <>
1107 : Quadrature<1>(n)
1108 , ep(ep)
1109{
1110 Assert(n > 0, ExcMessage("Need at least one point for quadrature rules"));
1111 std::vector<double> p =
1113 std::vector<double> w =
1115
1116 for (unsigned int i = 0; i < this->size(); ++i)
1117 {
1118 this->quadrature_points[i] = Point<1>(p[i]);
1119 this->weights[i] = w[i];
1120 }
1121}
1122
1123
1124template <int dim>
1126 EndPoint ep)
1128 n,
1129 static_cast<QGaussRadauChebyshev<1>::EndPoint>(ep)))
1130 , ep(ep)
1131{}
1132
1133
1134
1135namespace internal
1136{
1138 {
1139 // Computes the points of the quadrature formula.
1140 std::vector<double>
1141 get_quadrature_points(const unsigned int n)
1142 {
1143 std::vector<double> points(n);
1144 // n point quadrature: index from 0 to n-1
1145 for (unsigned short i = 0; i < n; ++i)
1146 // would be cos(i Pi/N)
1147 // put + Pi so we start from the smallest point
1148 // then map from [-1,1] to [0,1]
1149 points[i] =
1150 1. / 2. *
1151 (1. + std::cos(numbers::PI * (1 + double(i) / double(n - 1))));
1152
1153 return points;
1154 }
1155
1156 // Computes the weights of the quadrature formula.
1157 std::vector<double>
1158 get_quadrature_weights(const unsigned int n)
1159 {
1160 std::vector<double> weights(n);
1161
1162 for (unsigned short i = 0; i < n; ++i)
1163 {
1164 // same weights as on [-1,1]
1165 weights[i] = numbers::PI / double((n - 1));
1166 if (i == 0 || i == (n - 1))
1167 weights[i] /= 2.;
1168 }
1169
1170 return weights;
1171 }
1172 } // namespace QGaussLobattoChebyshev
1173} // namespace internal
1174
1175
1176
1177template <>
1179 : Quadrature<1>(n)
1180{
1181 Assert(n > 1,
1182 ExcMessage(
1183 "Need at least two points for Gauss-Lobatto quadrature rule"));
1184 std::vector<double> p =
1186 std::vector<double> w =
1188
1189 for (unsigned int i = 0; i < this->size(); ++i)
1190 {
1191 this->quadrature_points[i] = Point<1>(p[i]);
1192 this->weights[i] = w[i];
1193 }
1194}
1195
1196
1197template <int dim>
1200{}
1201
1202
1203
1204template <int dim>
1206{
1207 std::vector<Point<dim>> qpoints;
1208 std::vector<double> weights;
1209
1210 for (unsigned int i = 0; i < quad.size(); ++i)
1211 {
1212 double r = 0;
1213 /* Use "int d" instead of the more natural "unsigned int d" to work
1214 * around a wrong diagnostic in gcc-10.3.0 that warns about that the
1215 * comparison "d < dim" is always false in case of "dim == 0".
1216 * MM 2021 */
1217 for (int d = 0; d < dim; ++d)
1218 r += quad.point(i)[d];
1219 if (r <= 1 + 1e-10)
1220 {
1221 this->quadrature_points.push_back(quad.point(i));
1222 this->weights.push_back(quad.weight(i));
1223 }
1224 }
1225}
1226
1227
1228
1229template <int dim>
1232 const std::array<Point<dim>, dim + 1> &vertices) const
1233{
1235 for (unsigned int d = 0; d < dim; ++d)
1236 B[d] = vertices[d + 1] - vertices[0];
1237
1238 B = transpose(B);
1239 const double J = std::abs(determinant(B));
1240
1241 // if the determinant is zero, we return an empty quadrature
1242 if (J < 1e-12)
1243 return Quadrature<dim>();
1244
1245 std::vector<Point<dim>> qp(this->size());
1246 std::vector<double> w(this->size());
1247
1248 for (unsigned int i = 0; i < this->size(); ++i)
1249 {
1250 qp[i] = Point<dim>(vertices[0] + B * this->point(i));
1251 w[i] = J * this->weight(i);
1252 }
1253
1254 return Quadrature<dim>(qp, w);
1255}
1256
1257
1258
1260 const Quadrature<1> &angular_quadrature)
1261 : QSimplex<2>(Quadrature<2>())
1262{
1263 const QAnisotropic<2> base(radial_quadrature, angular_quadrature);
1264 this->quadrature_points.resize(base.size());
1265 this->weights.resize(base.size());
1266 for (unsigned int i = 0; i < base.size(); ++i)
1267 {
1268 const auto &q = base.point(i);
1269 const auto w = base.weight(i);
1270
1271 const auto xhat = q[0];
1272 const auto yhat = q[1];
1273
1274 const double t = numbers::PI_2 * yhat;
1275 const double pi = numbers::PI;
1276 const double st = std::sin(t);
1277 const double ct = std::cos(t);
1278 const double r = xhat / (st + ct);
1279
1280 const double J = pi * xhat / (2 * (std::sin(pi * yhat) + 1));
1281
1282 this->quadrature_points[i] = Point<2>(r * ct, r * st);
1283 this->weights[i] = w * J;
1284 }
1285}
1286
1287
1288
1290 : QTrianglePolar(QGauss<1>(n), QGauss<1>(n))
1291{}
1292
1293
1294
1295QDuffy::QDuffy(const Quadrature<1> &radial_quadrature,
1296 const Quadrature<1> &angular_quadrature,
1297 const double beta)
1298 : QSimplex<2>(Quadrature<2>())
1299{
1300 const QAnisotropic<2> base(radial_quadrature, angular_quadrature);
1301 this->quadrature_points.resize(base.size());
1302 this->weights.resize(base.size());
1303 for (unsigned int i = 0; i < base.size(); ++i)
1304 {
1305 const auto &q = base.point(i);
1306 const auto w = base.weight(i);
1307
1308 const auto xhat = q[0];
1309 const auto yhat = q[1];
1310
1311 const double x = std::pow(xhat, beta) * (1 - yhat);
1312 const double y = std::pow(xhat, beta) * yhat;
1313
1314 const double J = beta * std::pow(xhat, 2. * beta - 1.);
1315
1316 this->quadrature_points[i] = Point<2>(x, y);
1317 this->weights[i] = w * J;
1318 }
1319}
1320
1321
1322
1323QDuffy::QDuffy(const unsigned int n, const double beta)
1324 : QDuffy(QGauss<1>(n), QGauss<1>(n), beta)
1325{}
1326
1327
1328
1329template <int dim>
1331{
1333 ExcMessage(
1334 "The split point should be inside the unit reference cell."));
1335
1336 std::array<Point<dim>, dim + 1> vertices;
1337 vertices[0] = split_point;
1338
1339 // Make a simplex from the split_point and the first dim vertices of each
1340 // face. In dimension three, we need to split the face in two triangles, so
1341 // we use once the first dim vertices of each face, and the second time the
1342 // the dim vertices of each face starting from 1.
1343 for (auto f : GeometryInfo<dim>::face_indices())
1344 for (unsigned int start = 0; start < (dim > 2 ? 2 : 1); ++start)
1345 {
1346 for (unsigned int i = 0; i < dim; ++i)
1349 const auto quad = base.compute_affine_transformation(vertices);
1350 if (quad.size())
1351 {
1352 this->quadrature_points.insert(this->quadrature_points.end(),
1353 quad.get_points().begin(),
1354 quad.get_points().end());
1355 this->weights.insert(this->weights.end(),
1356 quad.get_weights().begin(),
1357 quad.get_weights().end());
1358 }
1359 }
1360}
1361
1362
1363
1364template <int dim>
1365QGaussSimplex<dim>::QGaussSimplex(const unsigned int n_points_1D)
1366 : QSimplex<dim>(Quadrature<dim>())
1367{
1368 // fill quadrature points and quadrature weights
1369 if (dim == 0 || dim == 1)
1370 {
1371 const ::QGauss<dim> quad(n_points_1D);
1372
1373 this->quadrature_points = quad.get_points();
1374 this->weights = quad.get_weights();
1375 }
1376 else if (dim == 2)
1377 {
1378 if (n_points_1D == 1)
1379 {
1380 const double p = 1.0 / 3.0;
1381 this->quadrature_points.emplace_back(p, p);
1382 this->weights.emplace_back(0.5);
1383 }
1384 else if (n_points_1D == 2)
1385 {
1386 // The Hillion 7 scheme, as communicated by quadpy
1387 //
1388 // See: Numerical Integration on a Triangle, International Journal for
1389 // Numerical Methods in Engineering, 1977
1390 const double Q12 = 1.0 / 2.0;
1391 this->quadrature_points.emplace_back(0.17855872826361643,
1392 0.1550510257216822);
1393 this->quadrature_points.emplace_back(0.07503111022260812,
1394 0.6449489742783178);
1395 this->quadrature_points.emplace_back(0.6663902460147014,
1396 0.1550510257216822);
1397 this->quadrature_points.emplace_back(0.28001991549907407,
1398 0.6449489742783178);
1399
1400 this->weights.emplace_back(0.31804138174397717 * Q12);
1401 this->weights.emplace_back(0.18195861825602283 * Q12);
1402 this->weights.emplace_back(0.31804138174397717 * Q12);
1403 this->weights.emplace_back(0.18195861825602283 * Q12);
1404 }
1405 else if (n_points_1D == 3)
1406 {
1407 // The Hammer-Marlowe-Stroud 5 Scheme, as communicated by quadpy
1408 const double p0 = 2.0 / 7.0 - std::sqrt(15.0) / 21.0;
1409 const double p1 = 2.0 / 7.0 + std::sqrt(15.0) / 21.0;
1410 const double p2 = 3.0 / 7.0 - 2.0 * std::sqrt(15.0) / 21.0;
1411 const double p3 = 3.0 / 7.0 + 2.0 * std::sqrt(15.0) / 21.0;
1412 this->quadrature_points.emplace_back(1.0 / 3.0, 1.0 / 3.0);
1413 this->quadrature_points.emplace_back(p3, p0);
1414 this->quadrature_points.emplace_back(p0, p3);
1415 this->quadrature_points.emplace_back(p0, p0);
1416 this->quadrature_points.emplace_back(p2, p1);
1417 this->quadrature_points.emplace_back(p1, p2);
1418 this->quadrature_points.emplace_back(p1, p1);
1419
1420 const double q12 = 0.5;
1421 const double w0 = 9.0 / 40.0;
1422 const double w1 = 31.0 / 240.0 - std::sqrt(15.0) / 1200.0;
1423 const double w2 = 31.0 / 240.0 + std::sqrt(15.0) / 1200.0;
1424 this->weights.emplace_back(q12 * w0);
1425 this->weights.emplace_back(q12 * w1);
1426 this->weights.emplace_back(q12 * w1);
1427 this->weights.emplace_back(q12 * w1);
1428 this->weights.emplace_back(q12 * w2);
1429 this->weights.emplace_back(q12 * w2);
1430 this->weights.emplace_back(q12 * w2);
1431 }
1432 else if (n_points_1D == 4)
1433 {
1435 QWitherdenVincentSimplex<dim>(n_points_1D));
1436 }
1437 }
1438 else if (dim == 3)
1439 {
1440 if (n_points_1D == 1)
1441 {
1442 const double Q14 = 1.0 / 4.0;
1443 const double Q16 = 1.0 / 6.0;
1444
1445 this->quadrature_points.emplace_back(Q14, Q14, Q14);
1446 this->weights.emplace_back(Q16);
1447 }
1448 // The Xiao Gimbutas 03 scheme, as communicated by quadpy
1449 //
1450 // See: A numerical algorithm for the construction of efficient quadrature
1451 // rules in two and higher dimensions, Computers & Mathematics with
1452 // Applications, 2010
1453 else if (n_points_1D == 2)
1454 {
1455 const double Q16 = 1.0 / 6.0;
1456 this->weights.emplace_back(0.1223220027573451 * Q16);
1457 this->weights.emplace_back(0.1280664127107469 * Q16);
1458 this->weights.emplace_back(0.1325680271444452 * Q16);
1459 this->weights.emplace_back(0.1406244096604032 * Q16);
1460 this->weights.emplace_back(0.2244151669175574 * Q16);
1461 this->weights.emplace_back(0.2520039808095023 * Q16);
1462
1463 this->quadrature_points.emplace_back(0.1620014916985245,
1464 0.1838503504920977,
1465 0.01271836631368145);
1466 this->quadrature_points.emplace_back(0.01090521221118924,
1467 0.2815238021235462,
1468 0.3621268299455338);
1469 this->quadrature_points.emplace_back(0.1901170024392839,
1470 0.01140332944455717,
1471 0.3586207204668839);
1472 this->quadrature_points.emplace_back(0.170816925164989,
1473 0.1528181430909273,
1474 0.6384932999617267);
1475 this->quadrature_points.emplace_back(0.1586851632274406,
1476 0.5856628056552158,
1477 0.1308471689520965);
1478 this->quadrature_points.emplace_back(0.5712260521491151,
1479 0.1469183900871696,
1480 0.1403728057942107);
1481 }
1482 // Past this point the best rules (positive weights, minimal number of
1483 // points) we have right now are the Witherden-Vincent ones
1484 else if (n_points_1D == 3)
1485 {
1487 QWitherdenVincentSimplex<dim>(n_points_1D));
1488 }
1489 else if (n_points_1D == 4)
1490 {
1492 QWitherdenVincentSimplex<dim>(n_points_1D));
1493 }
1494 }
1495
1496 AssertDimension(this->quadrature_points.size(), this->weights.size());
1497 Assert(this->quadrature_points.size() > 0,
1499 "QGaussSimplex is currently only implemented for "
1500 "n_points_1D = 1, 2, 3, and 4 while you are asking for "
1501 "n_points_1D = " +
1502 Utilities::to_string(n_points_1D)));
1503}
1504
1505namespace
1506{
1507 template <std::size_t b_dim>
1508 std::vector<std::array<double, b_dim>>
1509 all_permutations(const std::array<double, b_dim> &b_point)
1510 {
1511 std::vector<std::array<double, b_dim>> output;
1512
1513 // We want all possible permutations of the barycentric coordinates.
1514 // The easiest way to get all of them is to sort the input first and
1515 // then use next_permutation to cycle through them all.
1516 std::array<double, b_dim> temp = b_point;
1517 std::sort(temp.begin(), temp.end());
1518 do
1519 {
1520 output.push_back(temp);
1521 }
1522 while (std::next_permutation(temp.begin(), temp.end()));
1523
1524 return output;
1525 }
1526} // namespace
1527
1528
1529
1530template <int dim>
1532 const unsigned int n_points_1D,
1533 const bool use_odd_order)
1534 : QSimplex<dim>(Quadrature<dim>())
1535{
1536 Assert(1 <= dim && dim <= 3, ExcNotImplemented());
1537 // Just use Gauss in 1D: this is a high-order open rule so this is a
1538 // reasonable equivalent for generic programming.
1539 if (dim == 1)
1540 {
1542 return;
1543 }
1544
1545 std::array<double, dim + 1> centroid;
1546 std::fill(centroid.begin(), centroid.end(), 1.0 / (dim + 1.0));
1547 std::vector<std::vector<std::array<double, dim + 1>>> b_point_permutations;
1548 std::vector<double> b_weights;
1549
1550 // We can simplify the implementation of these quadrature rules
1551 // by quite a bit by exploiting symmetry - we do essentially the
1552 // same thing for each barycentric coordinate, so we can express
1553 // our quadrature rule as permutations of barycentric points
1554 // instead of writing things out explicitly.
1555
1556 // Apply a Barycentric permutation where one point is different.
1557 // Equivalent to d3_aa and s31 in quadpy.
1558 auto process_point_1 = [&](const double a, const double w) {
1559 const double b = 1.0 - dim * a;
1560 std::array<double, dim + 1> b_point;
1561 std::fill(b_point.begin(), b_point.begin() + dim, a);
1562 b_point[dim] = b;
1563
1564 b_weights.push_back(w);
1565 b_point_permutations.push_back(all_permutations(b_point));
1566 };
1567
1568 // Apply a Barycentric permutation where two points (in 3D) are different.
1569 // Equivalent to s22 in quadpy.
1570 auto process_point_2 = [&](const double a, const double w) {
1571 Assert(dim == 3, ExcInternalError());
1572 const double b = (1.0 - 2.0 * a) / 2.0;
1573 std::array<double, dim + 1> b_point;
1574 std::fill(b_point.begin(), b_point.begin() + dim - 1, a);
1575 b_point[dim - 1] = b;
1576 b_point[dim] = b;
1577
1578 b_weights.push_back(w);
1579 b_point_permutations.push_back(all_permutations(b_point));
1580 };
1581
1582 // Apply a Barycentric permutation where three (or four) points
1583 // are different (since there are two inputs).
1584 // Equivalent to d3_ab and s211 in quadpy.
1585 auto process_point_3 = [&](const double a, const double b, const double w) {
1586 const double c = 1.0 - (dim - 1.0) * a - b;
1587 std::array<double, dim + 1> b_point;
1588 std::fill(b_point.begin(), b_point.begin() + dim - 1, a);
1589 b_point[dim - 1] = b;
1590 b_point[dim] = c;
1591
1592 b_weights.push_back(w);
1593 b_point_permutations.push_back(all_permutations(b_point));
1594 };
1595
1596 switch (n_points_1D)
1597 {
1598 case 1:
1599 switch (dim)
1600 {
1601 case 2:
1602 if (use_odd_order)
1603 {
1604 // WV-1, 2D
1605 b_point_permutations.push_back({centroid});
1606 b_weights.push_back(1.0000000000000000e+00);
1607 }
1608 else
1609 {
1610 // WV-2, 2D
1611 process_point_1(1.6666666666666669e-01,
1612 3.3333333333333331e-01);
1613 }
1614 break;
1615 case 3:
1616 if (use_odd_order)
1617 {
1618 // WV-1, 3D
1619 b_point_permutations.push_back({centroid});
1620 b_weights.push_back(1.0000000000000000e+00);
1621 }
1622 else
1623 {
1624 // WV-2, 3D
1625 process_point_1(1.3819660112501050e-01,
1626 2.5000000000000000e-01);
1627 }
1628 break;
1629 default:
1630 Assert(false, ExcNotImplemented());
1631 }
1632 break;
1633 case 2:
1634 switch (dim)
1635 {
1636 case 2:
1637 // WV-4 in both cases (no WV-3 in 2D)
1638 process_point_1(9.1576213509770743e-02, 1.0995174365532187e-01);
1639 process_point_1(4.4594849091596489e-01, 2.2338158967801147e-01);
1640 break;
1641 case 3:
1642 if (use_odd_order)
1643 {
1644 // WV-3, 3D
1645 process_point_1(3.2816330251638171e-01,
1646 1.3621784253708741e-01);
1647 process_point_1(1.0804724989842859e-01,
1648 1.1378215746291261e-01);
1649 }
1650 else
1651 {
1652 // WV-5 (no WV-4 in 3D)
1654 }
1655 break;
1656 default:
1657 Assert(false, ExcInternalError());
1658 }
1659 break;
1660 case 3:
1661 switch (dim)
1662 {
1663 case 2:
1664 if (use_odd_order)
1665 {
1666 // WV-5, 2D
1667 b_point_permutations.push_back({centroid});
1668 b_weights.push_back(2.2500000000000001e-01);
1669 process_point_1(1.0128650732345634e-01,
1670 1.2593918054482714e-01);
1671 process_point_1(4.7014206410511511e-01,
1672 1.3239415278850619e-01);
1673 }
1674 else
1675 {
1676 // WV-6, 2D
1677 process_point_1(6.3089014491502227e-02,
1678 5.0844906370206819e-02);
1679 process_point_1(2.4928674517091043e-01,
1680 1.1678627572637937e-01);
1681 process_point_3(5.3145049844816938e-02,
1682 3.1035245103378439e-01,
1683 8.2851075618373571e-02);
1684 }
1685 break;
1686 case 3:
1687 if (use_odd_order)
1688 {
1689 // WV-5, 3D
1690 process_point_1(3.1088591926330061e-01,
1691 1.1268792571801590e-01);
1692 process_point_1(9.2735250310891248e-02,
1693 7.3493043116361956e-02);
1694 process_point_2(4.5503704125649642e-02,
1695 4.2546020777081472e-02);
1696 }
1697 else
1698 {
1699 // WV-6, 3D
1700 process_point_1(4.0673958534611372e-02,
1701 1.0077211055320640e-02);
1702 process_point_1(3.2233789014227548e-01,
1703 5.5357181543654717e-02);
1704 process_point_1(2.1460287125915201e-01,
1705 3.9922750258167487e-02);
1706 process_point_3(6.3661001875017442e-02,
1707 6.0300566479164919e-01,
1708 4.8214285714285710e-02);
1709 }
1710 break;
1711 default:
1712 Assert(false, ExcInternalError());
1713 }
1714 break;
1715 case 4:
1716 switch (dim)
1717 {
1718 case 2:
1719 if (use_odd_order)
1720 {
1721 // WV-7, 2D
1722 process_point_1(3.3730648554587850e-02,
1723 1.6545050110792131e-02);
1724 process_point_1(4.7430969250471822e-01,
1725 7.7086646185986069e-02);
1726 process_point_1(2.4157738259540357e-01,
1727 1.2794417123015558e-01);
1728 process_point_3(4.7036644652595216e-02,
1729 1.9868331479735168e-01,
1730 5.5878732903199779e-02);
1731 }
1732 else
1733 {
1734 // WV-8, 2D
1735 b_point_permutations.push_back({centroid});
1736 b_weights.push_back(1.4431560767778717e-01);
1737 process_point_1(5.0547228317030957e-02,
1738 3.2458497623198079e-02);
1739 process_point_1(4.5929258829272313e-01,
1740 9.5091634267284619e-02);
1741 process_point_1(1.7056930775176021e-01,
1742 1.0321737053471824e-01);
1743 process_point_3(8.3947774099575878e-03,
1744 2.6311282963463811e-01,
1745 2.7230314174434993e-02);
1746 }
1747 break;
1748 case 3:
1749 if (use_odd_order)
1750 {
1751 // WV-7, 3D
1752 b_point_permutations.push_back({centroid});
1753 b_weights.push_back(9.5485289464130846e-02);
1754 process_point_1(3.1570114977820279e-01,
1755 4.2329581209967028e-02);
1756 process_point_2(5.0489822598396350e-02,
1757 3.1896927832857580e-02);
1758 process_point_3(1.8883383102600099e-01,
1759 5.7517163758699996e-01,
1760 3.7207130728334620e-02);
1761 process_point_3(2.1265472541483140e-02,
1762 8.1083024109854862e-01,
1763 8.1107708299033420e-03);
1764 }
1765 else
1766 {
1767 // WV-8, 3D
1768 process_point_1(1.0795272496221089e-01,
1769 2.6426650908408830e-02);
1770 process_point_1(1.8510948778258660e-01,
1771 5.2031747563738531e-02);
1772 process_point_1(4.2316543684767283e-02,
1773 7.5252561535401989e-03);
1774 process_point_1(3.1418170912403898e-01,
1775 4.1763782856934897e-02);
1776 process_point_2(4.3559132858383021e-01,
1777 3.6280930261308818e-02);
1778 process_point_3(2.1433930127130570e-02,
1779 7.1746406342630831e-01,
1780 7.1569028908444327e-03);
1781 process_point_3(2.0413933387602909e-01,
1782 5.8379737830214440e-01,
1783 1.5453486150960340e-02);
1784 }
1785 break;
1786 default:
1787 Assert(false, ExcInternalError());
1788 }
1789 break;
1790 case 5:
1791 switch (dim)
1792 {
1793 case 2:
1794 if (use_odd_order)
1795 {
1796 // WV-9, 2D
1797 b_point_permutations.push_back({centroid});
1798 b_weights.push_back(9.7135796282798836e-02);
1799 process_point_1(4.4729513394452691e-02,
1800 2.5577675658698031e-02);
1801 process_point_1(4.8968251919873762e-01,
1802 3.1334700227139071e-02);
1803 process_point_1(4.3708959149293664e-01,
1804 7.7827541004774278e-02);
1805 process_point_1(1.8820353561903275e-01,
1806 7.9647738927210249e-02);
1807 process_point_3(3.6838412054736258e-02,
1808 2.2196298916076568e-01,
1809 4.3283539377289376e-02);
1810 }
1811 else
1812 {
1813 // WV-10, 2D
1814 b_point_permutations.push_back({centroid});
1815 b_weights.push_back(8.1743329146285973e-02);
1816 process_point_1(3.2055373216943517e-02,
1817 1.3352968813149567e-02);
1818 process_point_1(1.4216110105656438e-01,
1819 4.5957963604744731e-02);
1820 process_point_3(2.8367665339938453e-02,
1821 1.6370173373718250e-01,
1822 2.5297757707288385e-02);
1823 process_point_3(2.9619889488729734e-02,
1824 3.6914678182781102e-01,
1825 3.4184648162959429e-02);
1826 process_point_3(1.4813288578382056e-01,
1827 3.2181299528883545e-01,
1828 6.3904906396424044e-02);
1829 }
1830 break;
1831 case 3:
1832 if (use_odd_order)
1833 {
1834 // WV-9, 3D
1835 b_point_permutations.push_back({centroid});
1836 b_weights.push_back(5.8010548912480253e-02);
1837 process_point_1(6.1981697552226933e-10,
1838 6.4319281759256394e-05);
1839 process_point_1(1.6077453539526160e-01,
1840 2.3173338462425461e-02);
1841 process_point_1(3.2227652182142102e-01,
1842 2.9562912335429289e-02);
1843 process_point_1(4.5108918345413578e-02,
1844 8.0639799796161822e-03);
1845 process_point_2(1.1229654600437609e-01,
1846 3.8134080103702457e-02);
1847 process_point_3(4.5887144875245922e-01,
1848 2.5545792330413102e-03,
1849 8.3844221982985519e-03);
1850 process_point_3(3.3775870685338598e-02,
1851 7.1835032644207453e-01,
1852 1.0234559352745330e-02);
1853 process_point_3(1.8364136980992790e-01,
1854 3.4415910578175279e-02,
1855 2.0524915967988139e-02);
1856 }
1857 else
1858 {
1859 // WV-10, 3D
1860 b_point_permutations.push_back({centroid});
1861 b_weights.push_back(4.7399773556020743e-02);
1862 process_point_1(3.1225006869518868e-01,
1863 2.6937059992268701e-02);
1864 process_point_1(1.1430965385734609e-01,
1865 9.8691597167933822e-03);
1866 process_point_3(4.1043073921896539e-01,
1867 1.6548602561961109e-01,
1868 1.1393881220195230e-02);
1869 process_point_3(6.1380088247906528e-03,
1870 9.4298876734520487e-01,
1871 3.6194434433925362e-04);
1872 process_point_3(1.2105018114558939e-01,
1873 4.7719037990428043e-01,
1874 2.5739731980456069e-02);
1875 process_point_3(3.2779468216442620e-02,
1876 5.9425626948000698e-01,
1877 1.0135871679755789e-02);
1878 process_point_3(3.2485281564823047e-02,
1879 8.0117728465834437e-01,
1880 6.5761472770359038e-03);
1881 process_point_3(1.7497934218393901e-01,
1882 6.2807184547536599e-01,
1883 1.2907035798861989e-02);
1884 }
1885 break;
1886 default:
1887 Assert(false, ExcNotImplemented());
1888 }
1889 break;
1890 case 6:
1891 // There is no WV-11 rule in 3D yet
1892 Assert(dim == 2, ExcNotImplemented());
1893 if (use_odd_order)
1894 {
1895 // WV-11, 2D
1896 b_point_permutations.push_back({centroid});
1897 b_weights.push_back(8.5761179732224219e-02);
1898 process_point_1(2.8485417614371900e-02, 1.0431870512894697e-02);
1899 process_point_1(4.9589190096589092e-01, 1.6606273054585369e-02);
1900 process_point_1(1.0263548271224643e-01, 3.8630759237019321e-02);
1901 process_point_1(4.3846592676435220e-01, 6.7316154079468296e-02);
1902 process_point_1(2.1021995670317828e-01, 7.0515684111716576e-02);
1903 process_point_3(7.3254276860644785e-03,
1904 1.4932478865208237e-01,
1905 1.0290289572953278e-02);
1906 process_point_3(4.6010500165429957e-02,
1907 2.8958112563770588e-01,
1908 4.0332476640500554e-02);
1909 }
1910 else
1911 {
1912 // WV-12, 2D
1913 process_point_1(2.4646363436335583e-02, 7.9316425099736389e-03);
1914 process_point_1(4.8820375094554153e-01, 2.4266838081452032e-02);
1915 process_point_1(1.0925782765935427e-01, 2.8486052068877544e-02);
1916 process_point_1(4.4011164865859309e-01, 4.9918334928060942e-02);
1917 process_point_1(2.7146250701492608e-01, 6.2541213195902765e-02);
1918 process_point_3(2.1382490256170616e-02,
1919 1.2727971723358933e-01,
1920 1.5083677576511438e-02);
1921 process_point_3(2.3034156355267121e-02,
1922 2.9165567973834094e-01,
1923 2.1783585038607559e-02);
1924 process_point_3(1.1629601967792658e-01,
1925 2.5545422863851736e-01,
1926 4.3227363659414209e-02);
1927 }
1928 break;
1929 case 7:
1930 // There is no WV-13 rule in 3D yet
1931 Assert(dim == 2, ExcNotImplemented());
1932 if (use_odd_order)
1933 {
1934 // WV-13, 2D
1935 b_point_permutations.push_back({centroid});
1936 b_weights.push_back(6.7960036586831640e-02);
1937 process_point_1(2.1509681108843159e-02, 6.0523371035391717e-03);
1938 process_point_1(4.8907694645253935e-01, 2.3994401928894731e-02);
1939 process_point_1(4.2694141425980042e-01, 5.5601967530453329e-02);
1940 process_point_1(2.2137228629183292e-01, 5.8278485119199981e-02);
1941 process_point_3(5.1263891023823893e-03,
1942 2.7251581777342970e-01,
1943 9.5906810035432631e-03);
1944 process_point_3(2.4370186901093827e-02,
1945 1.1092204280346341e-01,
1946 1.4965401105165668e-02);
1947 process_point_3(8.7895483032197297e-02,
1948 1.6359740106785048e-01,
1949 2.4179039811593819e-02);
1950 process_point_3(6.8012243554206653e-02,
1951 3.0844176089211778e-01,
1952 3.4641276140848373e-02);
1953 }
1954 else
1955 {
1956 // WV-14, 2D
1957 process_point_1(1.9390961248701044e-02, 4.9234036024000819e-03);
1958 process_point_1(6.1799883090872587e-02, 1.4433699669776668e-02);
1959 process_point_1(4.8896391036217862e-01, 2.1883581369428889e-02);
1960 process_point_1(4.1764471934045394e-01, 3.2788353544125348e-02);
1961 process_point_1(1.7720553241254344e-01, 4.2162588736993016e-02);
1962 process_point_1(2.7347752830883865e-01, 5.1774104507291585e-02);
1963 process_point_3(1.2683309328720416e-03,
1964 1.1897449769695684e-01,
1965 5.0102288385006719e-03);
1966 process_point_3(1.4646950055654417e-02,
1967 2.9837288213625779e-01,
1968 1.4436308113533840e-02);
1969 process_point_3(5.7124757403647919e-02,
1970 1.7226668782135557e-01,
1971 2.4665753212563674e-02);
1972 process_point_3(9.2916249356971847e-02,
1973 3.3686145979634496e-01,
1974 3.8571510787060684e-02);
1975 }
1976 break;
1977 default:
1978 Assert(false, ExcNotImplemented());
1979 }
1980
1981 Assert(b_point_permutations.size() == b_weights.size(), ExcInternalError());
1982 for (unsigned int permutation_n = 0; permutation_n < b_weights.size();
1983 ++permutation_n)
1984 {
1985 for (const std::array<double, dim + 1> &b_point :
1986 b_point_permutations[permutation_n])
1987 {
1988 const double volume = (dim == 2 ? 1.0 / 2.0 : 1.0 / 6.0);
1989 this->weights.emplace_back(volume * b_weights[permutation_n]);
1990 Point<dim> c_point;
1991 std::copy(b_point.begin(),
1992 b_point.begin() + dim,
1993 c_point.begin_raw());
1994 this->quadrature_points.emplace_back(c_point);
1995 }
1996 }
1997}
1998
1999
2000
2001namespace
2002{
2003 template <int dim>
2005 setup_qiterated_1D(const Quadrature<dim> &, const unsigned int)
2006 {
2007 Assert(false, ExcInternalError());
2008 return Quadrature<dim>();
2009 }
2010
2011
2012
2014 setup_qiterated_1D(const Quadrature<1> &base_quad,
2015 const unsigned int n_copies)
2016 {
2017 return QIterated<1>(base_quad, n_copies);
2018 }
2019} // namespace
2020
2021
2022
2023template <int dim>
2025 const unsigned int n_copies)
2026{
2027 switch (dim)
2028 {
2029 case 1:
2030 static_cast<Quadrature<dim> &>(*this) =
2031 setup_qiterated_1D(base_quad, n_copies);
2032 break;
2033 case 2:
2034 case 3:
2035 {
2036 const auto n_refinements =
2037 static_cast<unsigned int>(std::round(std::log2(n_copies)));
2038 Assert((1u << n_refinements) == n_copies,
2039 ExcMessage("The number of copies must be a power of 2."));
2041 const auto reference_cell = ReferenceCells::get_simplex<dim>();
2043 tria.refine_global(n_refinements);
2044 const Mapping<dim> &mapping =
2045 reference_cell.template get_default_linear_mapping<dim>();
2047
2048 FEValues<dim> fe_values(mapping,
2049 fe,
2050 base_quad,
2052 std::vector<Point<dim>> points;
2053 std::vector<double> weights;
2054 for (const auto &cell : tria.active_cell_iterators())
2055 {
2056 fe_values.reinit(cell);
2057 for (unsigned int qp = 0; qp < base_quad.size(); ++qp)
2058 {
2059 points.push_back(fe_values.quadrature_point(qp));
2060 weights.push_back(fe_values.JxW(qp));
2061 }
2062 }
2063
2064 static_cast<Quadrature<dim> &>(*this) =
2065 Quadrature<dim>(points, weights);
2066
2067 break;
2068 }
2069 default:
2070 Assert(false, ExcNotImplemented());
2071 }
2072}
2073
2074
2075
2076template <int dim>
2077QGaussWedge<dim>::QGaussWedge(const unsigned int n_points)
2078 : Quadrature<dim>()
2079{
2080 AssertDimension(dim, 3);
2081
2082 QGaussSimplex<2> quad_tri(n_points);
2083 QGauss<1> quad_line(n_points);
2084
2085 for (unsigned int i = 0; i < quad_line.size(); ++i)
2086 for (unsigned int j = 0; j < quad_tri.size(); ++j)
2087 {
2088 this->quadrature_points.emplace_back(quad_tri.point(j)[0],
2089 quad_tri.point(j)[1],
2090 quad_line.point(i)[0]);
2091 this->weights.emplace_back(quad_tri.weight(j) * quad_line.weight(i));
2092 }
2093
2094 AssertDimension(this->quadrature_points.size(), this->weights.size());
2095 Assert(this->quadrature_points.size() > 0,
2096 ExcMessage("No valid quadrature points!"));
2097}
2098
2099
2100
2101template <int dim>
2102QGaussPyramid<dim>::QGaussPyramid(const unsigned int n_points_1D)
2103 : Quadrature<dim>()
2104{
2105 AssertDimension(dim, 3);
2106
2107 if (n_points_1D == 1)
2108 {
2109 const double Q14 = 1.0 / 4.0;
2110 const double Q43 = 4.0 / 3.0;
2111
2112 this->quadrature_points.emplace_back(0, 0, Q14);
2113 this->weights.emplace_back(Q43);
2114 }
2115 else if (n_points_1D == 2)
2116 {
2117 // clang-format off
2118 this->quadrature_points.emplace_back(-0.26318405556971, -0.26318405556971, 0.54415184401122);
2119 this->quadrature_points.emplace_back(-0.50661630334979, -0.50661630334979, 0.12251482265544);
2120 this->quadrature_points.emplace_back(-0.26318405556971, +0.26318405556971, 0.54415184401122);
2121 this->quadrature_points.emplace_back(-0.50661630334979, +0.50661630334979, 0.12251482265544);
2122 this->quadrature_points.emplace_back(+0.26318405556971, -0.26318405556971, 0.54415184401122);
2123 this->quadrature_points.emplace_back(+0.50661630334979, -0.50661630334979, 0.12251482265544);
2124 this->quadrature_points.emplace_back(+0.26318405556971, +0.26318405556971, 0.54415184401122);
2125 this->quadrature_points.emplace_back(+0.50661630334979, +0.50661630334979, 0.12251482265544);
2126 // clang-format on
2127
2128 this->weights.emplace_back(0.10078588207983);
2129 this->weights.emplace_back(0.23254745125351);
2130 this->weights.emplace_back(0.10078588207983);
2131 this->weights.emplace_back(0.23254745125351);
2132 this->weights.emplace_back(0.10078588207983);
2133 this->weights.emplace_back(0.23254745125351);
2134 this->weights.emplace_back(0.10078588207983);
2135 this->weights.emplace_back(0.23254745125351);
2136 }
2137
2138 AssertDimension(this->quadrature_points.size(), this->weights.size());
2139 Assert(this->quadrature_points.size() > 0,
2140 ExcMessage("No valid quadrature points!"));
2141}
2142
2143
2144
2145// explicit specialization
2146// note that 1d formulae are specialized by implementation above
2147template class QGauss<2>;
2148template class QGaussLobatto<2>;
2149template class QMidpoint<2>;
2150template class QTrapezoid<2>;
2151template class QSimpson<2>;
2152template class QMilne<2>;
2153template class QWeddle<2>;
2154
2155template class QGauss<3>;
2156template class QGaussLobatto<3>;
2157template class QMidpoint<3>;
2158template class QTrapezoid<3>;
2159template class QSimpson<3>;
2160template class QMilne<3>;
2161template class QWeddle<3>;
2162
2163template class QSorted<1>;
2164template class QSorted<2>;
2165template class QSorted<3>;
2166
2167template class QTelles<1>;
2168template class QTelles<2>;
2169template class QTelles<3>;
2170
2171template class QGaussChebyshev<1>;
2172template class QGaussChebyshev<2>;
2173template class QGaussChebyshev<3>;
2174
2175template class QGaussRadauChebyshev<1>;
2176template class QGaussRadauChebyshev<2>;
2177template class QGaussRadauChebyshev<3>;
2178
2179template class QGaussLobattoChebyshev<1>;
2180template class QGaussLobattoChebyshev<2>;
2181template class QGaussLobattoChebyshev<3>;
2182
2183template class QSimplex<1>;
2184template class QSimplex<2>;
2185template class QSimplex<3>;
2186
2187template class QIteratedSimplex<1>;
2188template class QIteratedSimplex<2>;
2189template class QIteratedSimplex<3>;
2190
2191template class QSplit<1>;
2192template class QSplit<2>;
2193template class QSplit<3>;
2194
2195template class QGaussSimplex<0>;
2196template class QGaussSimplex<1>;
2197template class QGaussSimplex<2>;
2198template class QGaussSimplex<3>;
2199template class QGaussWedge<1>;
2200template class QGaussWedge<2>;
2201template class QGaussWedge<3>;
2202template class QGaussPyramid<1>;
2203template class QGaussPyramid<2>;
2204template class QGaussPyramid<3>;
2205
2206template class QWitherdenVincentSimplex<1>;
2207template class QWitherdenVincentSimplex<2>;
2208template class QWitherdenVincentSimplex<3>;
2209
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
double JxW(const unsigned int quadrature_point) const
const Point< spacedim > & quadrature_point(const unsigned int q) const
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell)
Abstract base class for mapping classes.
Definition: mapping.h:311
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, EndPoint ep=QGaussRadauChebyshev::left)
Generate a formula with n quadrature points.
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)
QSimplex(const Quadrature< dim > &quad)
Quadrature< dim > compute_affine_transformation(const std::array< Point< dim >, dim+1 > &vertices) 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)
std::vector< Point< dim > > quadrature_points
Definition: quadrature.h:290
Quadrature & operator=(const Quadrature< dim > &)
Definition: quadrature.cc:284
const Point< dim > & point(const unsigned int i) const
bool is_tensor_product_flag
Definition: quadrature.h:305
double weight(const unsigned int i) const
std::vector< double > weights
Definition: quadrature.h:296
unsigned int size() const
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &t)
Definition: tensor.h:503
Number * begin_raw()
Number * end_raw()
void refine_global(const unsigned int times=1)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
Point< 3 > vertices[4]
@ update_JxW_values
Transformed quadrature weights.
@ update_quadrature_points
Transformed quadrature points.
IteratorRange< active_cell_iterator > active_cell_iterators() const
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
Definition: exceptions.h:1473
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1667
#define AssertIndexRange(index, range)
Definition: exceptions.h:1732
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
double volume(const Triangulation< dim, spacedim > &tria, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
Definition: grid_tools.cc:139
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Definition: divergence.h:472
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:190
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={})
Definition: generators.cc:493
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(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)
Definition: polynomial.h:1037
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:482
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_points(const unsigned int n, ::QGaussRadauChebyshev< 1 >::EndPoint ep)
std::vector< double > get_quadrature_weights(const unsigned int n, ::QGaussRadauChebyshev< 1 >::EndPoint ep)
void split_point(const Point< dim1+dim2 > &source, Point< dim1 > &p1, Point< dim2 > &p2)
void copy(const T *begin, const T *end, U *dest)
static constexpr double PI_2
Definition: numbers.h:238
static constexpr double PI
Definition: numbers.h:233
::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 Point< dim > unit_cell_vertex(const unsigned int vertex)
const ::Triangulation< dim, spacedim > & tria