Reference documentation for deal.II version 9.5.0
Searching...
No Matches
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 1998 - 2023 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
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
20
23
26#include <deal.II/grid/tria.h>
27
28#include <algorithm>
29#include <cmath>
30#include <functional>
31#include <limits>
32
33
35
36
37// please note: for a given dimension, we need the quadrature formulae
38// for all lower dimensions as well. That is why in this file the check
39// is for deal_II_dimension >= any_number and not for ==
40
41
42
43template <>
44QGauss<0>::QGauss(const unsigned int)
45 : // there are n_q^dim == 1
46 // points
48{
49 // the single quadrature point gets unit
50 // weight
51 this->weights[0] = 1;
52}
53
54
55
56template <>
58 : // there are n_q^dim == 1
59 // points
61{
62 // the single quadrature point gets unit
63 // weight
64 this->weights[0] = 1;
65}
66
67
68
69template <>
70QGauss<1>::QGauss(const unsigned int n)
72{
73 if (n == 0)
74 return;
75
76 std::vector<long double> points =
77 Polynomials::jacobi_polynomial_roots<long double>(n, 0, 0);
78
79 for (unsigned int i = 0; i < (points.size() + 1) / 2; ++i)
80 {
82 this->quadrature_points[n - i - 1][0] = 1. - points[i];
83
84 // derivative of Jacobi polynomial
85 const long double pp =
86 0.5 * (n + 1) *
87 Polynomials::jacobi_polynomial_value(n - 1, 1, 1, points[i]);
88 const long double x = -1. + 2. * points[i];
89 const double w = 1. / ((1. - x * x) * pp * pp);
90 this->weights[i] = w;
91 this->weights[n - i - 1] = w;
92 }
93}
94
95namespace internal
96{
97 namespace QGaussLobatto
98 {
103 long double
104 gamma(const unsigned int n)
105 {
106 long double result = n - 1;
107 for (int i = n - 2; i > 1; --i)
108 result *= i;
109 return result;
110 }
111
112
113
121 std::vector<long double>
123 const int alpha,
124 const int beta)
125 {
126 const unsigned int q = x.size();
127 std::vector<long double> w(q);
128
129 const long double factor =
130 std::pow(2., alpha + beta + 1) * gamma(alpha + q) * gamma(beta + q) /
131 ((q - 1) * gamma(q) * gamma(alpha + beta + q + 1));
132 for (unsigned int i = 0; i < q; ++i)
133 {
134 const long double s =
135 Polynomials::jacobi_polynomial_value(q - 1, alpha, beta, x[i]);
136 w[i] = factor / (s * s);
137 }
138 w[0] *= (beta + 1);
139 w[q - 1] *= (alpha + 1);
140
141 return w;
142 }
143 } // namespace QGaussLobatto
144} // namespace internal
145
146
147#ifndef DOXYGEN
148template <>
149QGaussLobatto<1>::QGaussLobatto(const unsigned int n)
151{
152 Assert(n >= 2, ExcNotImplemented());
153
154 std::vector<long double> points =
155 Polynomials::jacobi_polynomial_roots<long double>(n - 2, 1, 1);
156 points.insert(points.begin(), 0);
157 points.push_back(1.);
158 std::vector<long double> w =
160
161 // scale weights to the interval [0.0, 1.0]:
162 for (unsigned int i = 0; i < points.size(); ++i)
163 {
165 this->weights[i] = 0.5 * w[i];
166 }
167}
168#endif
169
170
171template <>
174{
176 this->weights[0] = 1.0;
177}
178
179
180
181template <>
184{
185 static const double xpts[] = {0.0, 1.0};
186 static const double wts[] = {0.5, 0.5};
187
188 for (unsigned int i = 0; i < this->size(); ++i)
189 {
191 this->weights[i] = wts[i];
192 }
193}
194
195
196
197template <>
200{
201 static const double xpts[] = {0.0, 0.5, 1.0};
202 static const double wts[] = {1. / 6., 2. / 3., 1. / 6.};
203
204 for (unsigned int i = 0; i < this->size(); ++i)
205 {
207 this->weights[i] = wts[i];
208 }
209}
210
211
212
213template <>
216{
217 static const double xpts[] = {0.0, .25, .5, .75, 1.0};
218 static const double wts[] = {
219 7. / 90., 32. / 90., 12. / 90., 32. / 90., 7. / 90.};
220
221 for (unsigned int i = 0; i < this->size(); ++i)
222 {
224 this->weights[i] = wts[i];
225 }
226}
227
228
229
230template <>
233{
234 static const double xpts[] = {
235 0.0, 1. / 6., 1. / 3., .5, 2. / 3., 5. / 6., 1.0};
236 static const double wts[] = {41. / 840.,
237 216. / 840.,
238 27. / 840.,
239 272. / 840.,
240 27. / 840.,
241 216. / 840.,
242 41. / 840.};
243
244 for (unsigned int i = 0; i < this->size(); ++i)
245 {
247 this->weights[i] = wts[i];
248 }
249}
250
251
252template <>
253QGaussLog<1>::QGaussLog(const unsigned int n, const bool revert)
255{
258
259 for (unsigned int i = 0; i < this->size(); ++i)
260 {
261 // Using the change of variables x=1-t, it's possible to show
262 // that int f(x)ln|1-x| = int f(1-t) ln|t|, which implies that
263 // we can use this quadrature formula also with weight ln|1-x|.
265 revert ? Point<1>(1 - p[n - 1 - i]) : Point<1>(p[i]);
266 this->weights[i] = revert ? w[n - 1 - i] : w[i];
267 }
268}
269
270template <>
271std::vector<double>
273{
274 std::vector<double> q_points(n);
275
276 switch (n)
277 {
278 case 1:
279 q_points[0] = 0.3333333333333333;
280 break;
281
282 case 2:
283 q_points[0] = 0.1120088061669761;
284 q_points[1] = 0.6022769081187381;
285 break;
286
287 case 3:
288 q_points[0] = 0.06389079308732544;
289 q_points[1] = 0.3689970637156184;
290 q_points[2] = 0.766880303938942;
291 break;
292
293 case 4:
294 q_points[0] = 0.04144848019938324;
295 q_points[1] = 0.2452749143206022;
296 q_points[2] = 0.5561654535602751;
297 q_points[3] = 0.848982394532986;
298 break;
299
300 case 5:
301 q_points[0] = 0.02913447215197205;
302 q_points[1] = 0.1739772133208974;
303 q_points[2] = 0.4117025202849029;
304 q_points[3] = 0.6773141745828183;
305 q_points[4] = 0.89477136103101;
306 break;
307
308 case 6:
309 q_points[0] = 0.02163400584411693;
310 q_points[1] = 0.1295833911549506;
311 q_points[2] = 0.3140204499147661;
312 q_points[3] = 0.5386572173517997;
313 q_points[4] = 0.7569153373774084;
314 q_points[5] = 0.922668851372116;
315 break;
316
317
318 case 7:
319 q_points[0] = 0.0167193554082585;
320 q_points[1] = 0.100185677915675;
321 q_points[2] = 0.2462942462079286;
322 q_points[3] = 0.4334634932570557;
323 q_points[4] = 0.6323509880476823;
324 q_points[5] = 0.81111862674023;
325 q_points[6] = 0.940848166743287;
326 break;
327
328 case 8:
329 q_points[0] = 0.01332024416089244;
330 q_points[1] = 0.07975042901389491;
331 q_points[2] = 0.1978710293261864;
332 q_points[3] = 0.354153994351925;
333 q_points[4] = 0.5294585752348643;
334 q_points[5] = 0.7018145299391673;
335 q_points[6] = 0.849379320441094;
336 q_points[7] = 0.953326450056343;
337 break;
338
339 case 9:
340 q_points[0] = 0.01086933608417545;
341 q_points[1] = 0.06498366633800794;
342 q_points[2] = 0.1622293980238825;
343 q_points[3] = 0.2937499039716641;
344 q_points[4] = 0.4466318819056009;
345 q_points[5] = 0.6054816627755208;
346 q_points[6] = 0.7541101371585467;
347 q_points[7] = 0.877265828834263;
348 q_points[8] = 0.96225055941096;
349 break;
350
351 case 10:
352 q_points[0] = 0.00904263096219963;
353 q_points[1] = 0.05397126622250072;
354 q_points[2] = 0.1353118246392511;
355 q_points[3] = 0.2470524162871565;
356 q_points[4] = 0.3802125396092744;
357 q_points[5] = 0.5237923179723384;
358 q_points[6] = 0.6657752055148032;
359 q_points[7] = 0.7941904160147613;
360 q_points[8] = 0.898161091216429;
361 q_points[9] = 0.9688479887196;
362 break;
363
364
365 case 11:
366 q_points[0] = 0.007643941174637681;
367 q_points[1] = 0.04554182825657903;
368 q_points[2] = 0.1145222974551244;
369 q_points[3] = 0.2103785812270227;
370 q_points[4] = 0.3266955532217897;
371 q_points[5] = 0.4554532469286375;
372 q_points[6] = 0.5876483563573721;
373 q_points[7] = 0.7139638500230458;
374 q_points[8] = 0.825453217777127;
375 q_points[9] = 0.914193921640008;
376 q_points[10] = 0.973860256264123;
377 break;
378
379 case 12:
380 q_points[0] = 0.006548722279080035;
381 q_points[1] = 0.03894680956045022;
382 q_points[2] = 0.0981502631060046;
383 q_points[3] = 0.1811385815906331;
384 q_points[4] = 0.2832200676673157;
385 q_points[5] = 0.398434435164983;
386 q_points[6] = 0.5199526267791299;
387 q_points[7] = 0.6405109167754819;
388 q_points[8] = 0.7528650118926111;
389 q_points[9] = 0.850240024421055;
390 q_points[10] = 0.926749682988251;
391 q_points[11] = 0.977756129778486;
392 break;
393
394 default:
395 Assert(false, ExcNotImplemented());
396 break;
397 }
398
399 return q_points;
400}
401
402
403template <>
404std::vector<double>
406{
408
409 switch (n)
410 {
411 case 1:
413 break;
414 case 2:
417 break;
418
419 case 3:
423 break;
424
425 case 4:
430 break;
431
432 case 5:
438 break;
439
440 case 6:
447 break;
448
449
450 case 7:
458 break;
459
460 case 8:
469 break;
470
471 case 9:
481 break;
482
483 case 10:
494 break;
495
496
497 case 11:
509 break;
510
511 case 12:
524 break;
525
526 default:
527 Assert(false, ExcNotImplemented());
528 break;
529 }
530
532}
533
534
535template <>
536QGaussLogR<1>::QGaussLogR(const unsigned int n,
537 const Point<1> & origin,
538 const double alpha,
539 const bool factor_out_singularity)
541 ((origin[0] == 0) || (origin[0] == 1)) ? (alpha == 1 ? n : 2 * n) : 4 * n)
542 , fraction(((origin[0] == 0) || (origin[0] == 1.)) ? 1. : origin[0])
543{
544 // The three quadrature formulas that make this one up. There are
545 // at most two when the origin is one of the extremes, and there is
546 // only one if the origin is one of the extremes and alpha is
547 // equal to one.
548 //
549 // If alpha is different from one, then we need a correction which
550 // is performed with a standard Gauss quadrature rule on each
551 // segment. This is not needed in the standard case where alpha is
552 // equal to one and the origin is on one of the extremes. We
553 // integrate with weight ln|(x-o)/alpha|. In the easy cases, we
554 // only need n quadrature points. In the most difficult one, we
555 // need 2*n points for the first segment, and 2*n points for the
556 // second segment.
557 QGaussLog<1> quad1(n, origin[0] != 0);
560
561 // Check that the origin is inside 0,1
562 Assert((fraction >= 0) && (fraction <= 1),
563 ExcMessage("Origin is outside [0,1]."));
564
565 // Non singular offset. This is the start of non singular quad
566 // points.
567 unsigned int ns_offset = (fraction == 1) ? n : 2 * n;
568
569 for (unsigned int i = 0, j = ns_offset; i < n; ++i, ++j)
570 {
571 // The first i quadrature points are the same as quad1, and
572 // are by default singular.
574 this->weights[i] = quad1.weight(i) * fraction;
575
576 // We need to scale with -log|fraction*alpha|
577 if ((alpha != 1) || (fraction != 1))
578 {
580 this->weights[j] =
581 -std::log(alpha / fraction) * quad.weight(i) * fraction;
582 }
583 // In case we need the second quadrature as well, do it now.
584 if (fraction != 1)
585 {
587 quad2.point(i) * (1 - fraction) + Point<1>(fraction);
588 this->weights[i + n] = quad2.weight(i) * (1 - fraction);
589
590 // We need to scale with -log|fraction*alpha|
592 quad.point(i) * (1 - fraction) + Point<1>(fraction);
593 this->weights[j + n] =
594 -std::log(alpha / (1 - fraction)) * quad.weight(i) * (1 - fraction);
595 }
596 }
597 if (factor_out_singularity == true)
598 for (unsigned int i = 0; i < size(); ++i)
599 {
600 Assert(
603 "The singularity cannot be on a Gauss point of the same order!"));
604 double denominator =
605 std::log(std::abs((this->quadrature_points[i] - origin)[0]) / alpha);
606 Assert(denominator != 0.0,
608 "The quadrature formula you are using does not allow to "
609 "factor out the singularity, which is zero at one point."));
610 this->weights[i] /= denominator;
612}
613
614
615template <>
616unsigned int
617QGaussOneOverR<2>::quad_size(const Point<2> &singularity, const unsigned int n)
618{
619 const double eps = 1e-8;
620 bool on_edge = false;
621 for (unsigned int d = 0; d < 2; ++d)
622 on_edge = on_edge || (std::abs(singularity[d]) < eps ||
623 std::abs(singularity[d] - 1.0) < eps);
624 const bool on_vertex =
625 on_edge &&
626 std::abs((singularity - Point<2>(.5, .5)).norm_square() - .5) < eps;
627 if (on_vertex)
628 return 2 * n * n;
629 else if (on_edge)
630 return 4 * n * n;
631 else
632 return 8 * n * n;
633}
634
635template <>
637 const Point<2> & singularity,
638 const bool factor_out_singularity)
640{
641 // We treat all the cases in the
642 // same way. Split the element in 4
643 // pieces, measure the area, if
644 // it's relevant, add the
645 // quadrature connected to that
646 // singularity.
648 std::vector<Point<2>> origins;
649 // Id of the corner with a
650 // singularity
655
656 origins.emplace_back(0., 0.);
657 origins.emplace_back(singularity[0], 0.);
658 origins.emplace_back(0., singularity[1]);
659 origins.push_back(singularity);
661 // Lexicographical ordering.
662
663 double eps = 1e-8;
664 unsigned int q_id = 0; // Current quad point index.
665 Tensor<1, 2> dist;
666
667 for (unsigned int box = 0; box < 4; ++box)
668 {
669 dist = (singularity - GeometryInfo<2>::unit_cell_vertex(box));
670 dist = Point<2>(std::abs(dist[0]), std::abs(dist[1]));
671 double area = dist[0] * dist[1];
672 if (area > eps)
673 for (unsigned int q = 0; q < quads[box].size(); ++q, ++q_id)
674 {
675 const Point<2> &qp = quads[box].point(q);
677 origins[box] + Point<2>(dist[0] * qp[0], dist[1] * qp[1]);
678 this->weights[q_id] = quads[box].weight(q) * area;
679 }
680 }
681}
682
683
684template <>
686 const unsigned int vertex_index,
687 const bool factor_out_singularity)
688 : Quadrature<2>(2 * n * n)
689{
690 // This version of the constructor
691 // works only for the 4
692 // vertices. If you need a more
693 // general one, you should use the
694 // one with the Point<2> in the
695 // constructor.
696 AssertIndexRange(vertex_index, 4);
697
699 // element.
700 QGauss<2> gauss(n);
701
702 Assert(gauss.size() == n * n, ExcInternalError());
703
704 // For the moment we only implemented this for the vertices of a
705 // quadrilateral. We are planning to do this also for the support
706 // points of arbitrary FE_Q elements, to allow the use of this
707 // class in boundary element programs with higher order mappings.
708 AssertIndexRange(vertex_index, 4);
709
710 // We create only the first one. All other pieces are rotation of
711 // this one.
712 // In this case the transformation is
713 //
714 // (x,y) = (u, u tan(pi/4 v))
715 //
716 // with Jacobian
717 //
718 // J = pi/4 R / cos(pi/4 v)
719 //
720 // And we get rid of R to take into account the singularity,
721 // unless specified differently in the constructor.
723 std::vector<double> & ws = this->weights;
724 double pi4 = numbers::PI / 4;
725
726 for (unsigned int q = 0; q < gauss.size(); ++q)
727 {
728 const Point<2> &gp = gauss.point(q);
729 ps[q][0] = gp[0];
730 ps[q][1] = gp[0] * std::tan(pi4 * gp[1]);
731 ws[q] = gauss.weight(q) * pi4 / std::cos(pi4 * gp[1]);
732 if (factor_out_singularity)
733 ws[q] *= (ps[q] - GeometryInfo<2>::unit_cell_vertex(0)).norm();
734 // The other half of the quadrilateral is symmetric with
735 // respect to xy plane.
736 ws[gauss.size() + q] = ws[q];
737 ps[gauss.size() + q][0] = ps[q][1];
738 ps[gauss.size() + q][1] = ps[q][0];
739 }
740
741 // Now we distribute these vertices in the correct manner
742 double theta = 0;
743 switch (vertex_index)
744 {
745 case 0:
746 theta = 0;
747 break;
748 case 1:
749 //
750 theta = numbers::PI / 2;
751 break;
752 case 2:
753 theta = -numbers::PI / 2;
754 break;
755 case 3:
756 theta = numbers::PI;
757 break;
758 }
759
760 double R00 = std::cos(theta), R01 = -std::sin(theta);
761 double R10 = std::sin(theta), R11 = std::cos(theta);
762
763 if (vertex_index != 0)
764 for (unsigned int q = 0; q < size(); ++q)
765 {
766 double x = ps[q][0] - .5, y = ps[q][1] - .5;
767
768 ps[q][0] = R00 * x + R01 * y + .5;
769 ps[q][1] = R10 * x + R11 * y + .5;
770 }
771}
772
773
774template <int dim>
777{
779 for (unsigned int i = 0; i < quad.size(); ++i)
780 permutation[i] = i;
781
782 std::sort(permutation.begin(),
783 permutation.end(),
784 [this](const unsigned int x, const unsigned int y) {
785 return this->compare_weights(x, y);
786 });
787
788 // At this point, the variable is_tensor_product_flag is set
789 // to the respective value of the given Quadrature in the base
790 // class copy constructor.
791 // We only call a quadrature formula 'tensor product'
792 // if the quadrature points are also sorted lexicographically.
793 // In particular, any reordering destroys that property
794 // and we might need to modify the variable accordingly.
795 for (unsigned int i = 0; i < quad.size(); ++i)
796 {
799 if (permutation[i] != i)
800 this->is_tensor_product_flag = false;
801 }
802}
803
804
805template <int dim>
806bool
807QSorted<dim>::compare_weights(const unsigned int a, const unsigned int b) const
808{
809 return (this->weights[a] < this->weights[b]);
810}
811
812
813// construct the quadrature formulae in higher dimensions by
814// tensor product of lower dimensions
815
816template <int dim>
817QGauss<dim>::QGauss(const unsigned int n)
818 : Quadrature<dim>(QGauss<dim - 1>(n), QGauss<1>(n))
819{}
820
821
822
823template <int dim>
825 : Quadrature<dim>(QGaussLobatto<dim - 1>(n), QGaussLobatto<1>(n))
826{}
827
828
829
830template <int dim>
832 : Quadrature<dim>(QMidpoint<dim - 1>(), QMidpoint<1>())
833{}
834
835
836
837template <int dim>
839 : Quadrature<dim>(QTrapezoid<dim - 1>(), QTrapezoid<1>())
840{}
841
842
843
844template <int dim>
846 : Quadrature<dim>(QSimpson<dim - 1>(), QSimpson<1>())
847{}
848
849
850
851template <int dim>
853 : Quadrature<dim>(QMilne<dim - 1>(), QMilne<1>())
854{}
855
856
857template <int dim>
859 : Quadrature<dim>(QWeddle<dim - 1>(), QWeddle<1>())
860{}
861
862template <int dim>
864 const Point<dim> & singularity)
865 : // We need the explicit implementation if dim == 1. If dim > 1 we use the
866 // former implementation and apply a tensorial product to obtain the higher
867 // dimensions.
869 dim == 2 ?
872 dim == 3 ?
877{}
878
879template <int dim>
880QTelles<dim>::QTelles(const unsigned int n, const Point<dim> &singularity)
881 : // In this case we map the standard Gauss Legendre formula using the given
882 // singularity point coordinates.
884{}
885
886
887
888template <>
890 : // We explicitly implement the Telles' variable change if dim == 1.
892{
893 // We define all the constants to be used in the implementation of
894 // Telles' rule
895 const double eta_bar = singularity[0] * 2. - 1.;
896 const double eta_star = eta_bar * eta_bar - 1.;
897 double gamma_bar;
898
900 std::vector<double> weights_dummy(weights.size());
901 unsigned int cont = 0;
902 const double tol = 1e-10;
903 for (unsigned int d = 0; d < quadrature_points.size(); ++d)
904 {
905 if (std::abs(quadrature_points[d][0] - singularity[0]) > tol)
906 {
908 weights_dummy[d - cont] = weights[d];
909 }
910 else
911 {
912 // We need to remove the singularity point from the quadrature point
913 // list. To do so we use the variable cont.
914 cont = 1;
915 }
916 }
917 if (cont == 1)
918 {
920 weights.resize(weights_dummy.size() - 1);
921 for (unsigned int d = 0; d < quadrature_points.size(); ++d)
922 {
924 weights[d] = weights_dummy[d];
925 }
926 }
927 // We need to check if the singularity is at the boundary of the interval.
928 if (std::abs(eta_star) <= tol)
929 {
930 gamma_bar =
931 std::pow((eta_bar * eta_star + std::abs(eta_star)), 1.0 / 3.0) +
932 std::pow((eta_bar * eta_star - std::abs(eta_star)), 1.0 / 3.0) +
933 eta_bar;
934 }
935 else
936 {
937 gamma_bar = (eta_bar * eta_star + std::abs(eta_star)) /
938 std::abs(eta_bar * eta_star + std::abs(eta_star)) *
939 std::pow(std::abs(eta_bar * eta_star + std::abs(eta_star)),
940 1.0 / 3.0) +
941 (eta_bar * eta_star - std::abs(eta_star)) /
942 std::abs(eta_bar * eta_star - std::abs(eta_star)) *
943 std::pow(std::abs(eta_bar * eta_star - std::abs(eta_star)),
944 1.0 / 3.0) +
945 eta_bar;
946 }
947 for (unsigned int q = 0; q < quadrature_points.size(); ++q)
948 {
949 double gamma = quadrature_points[q][0] * 2 - 1;
950 double eta = (Utilities::fixed_power<3>(gamma - gamma_bar) +
951 gamma_bar * (gamma_bar * gamma_bar + 3)) /
952 (1 + 3 * gamma_bar * gamma_bar);
953
954 double J = 3 * ((gamma - gamma_bar) * (gamma - gamma_bar)) /
955 (1 + 3 * gamma_bar * gamma_bar);
956
957 quadrature_points[q][0] = (eta + 1) / 2.0;
958 weights[q] = J * weights[q];
959 }
960}
961
962namespace internal
963{
965 {
969 std::vector<double>
971 {
972 std::vector<double> points(n);
973 // n point quadrature: index from 0 to n-1
974 for (unsigned short i = 0; i < n; ++i)
975 // would be cos((2i+1)Pi/(2N+2))
976 // put + Pi so we start from the smallest point
977 // then map from [-1,1] to [0,1]
978 points[i] =
979 1. / 2. *
980 (1. + std::cos(numbers::PI *
981 (1. + double(2 * i + 1) / double(2 * (n - 1) + 2))));
982
983 return points;
984 }
985
986
987
991 std::vector<double>
993 {
994 std::vector<double> weights(n);
995
996 for (unsigned short i = 0; i < n; ++i)
997 {
998 // same weights as on [-1,1]
999 weights[i] = numbers::PI / double(n);
1000 }
1001
1002 return weights;
1003 }
1004 } // namespace QGaussChebyshev
1005} // namespace internal
1006
1007
1008template <>
1011{
1012 Assert(n > 0, ExcMessage("Need at least one point for the quadrature rule"));
1015
1016 for (unsigned int i = 0; i < this->size(); ++i)
1017 {
1019 this->weights[i] = w[i];
1020 }
1021}
1022
1023
1024template <int dim>
1027{}
1028
1029
1030namespace internal
1031{
1033 {
1034 // Computes the points of the quadrature formula.
1035 std::vector<double>
1038 {
1039 std::vector<double> points(n);
1040 // n point quadrature: index from 0 to n-1
1041 for (unsigned short i = 0; i < n; ++i)
1042 // would be -cos(2i Pi/(2N+1))
1043 // put + Pi so we start from the smallest point
1044 // then map from [-1,1] to [0,1]
1045 switch (ep)
1046 {
1048 {
1049 points[i] =
1050 1. / 2. *
1051 (1. -
1053 (1 + 2 * double(i) / (2 * double(n - 1) + 1.))));
1054 break;
1055 }
1056
1058 {
1059 points[i] =
1060 1. / 2. *
1061 (1. - std::cos(numbers::PI * (2 * double(n - 1 - i) /
1062 (2 * double(n - 1) + 1.))));
1063 break;
1064 }
1065
1066 default:
1067 Assert(
1068 false,
1069 ExcMessage(
1070 "This constructor can only be called with either "
1072 "second argument."));
1073 }
1074
1075 return points;
1076 }
1077
1078
1079
1080 // Computes the weights of the quadrature formula.
1081 std::vector<double>
1084 {
1085 std::vector<double> weights(n);
1086
1087 for (unsigned short i = 0; i < n; ++i)
1088 {
1089 // same weights as on [-1,1]
1090 weights[i] = 2. * numbers::PI / double(2 * (n - 1) + 1.);
1091 if (ep == ::QGaussRadauChebyshev<1>::left && i == 0)
1092 weights[i] /= 2.;
1093 else if (ep == ::QGaussRadauChebyshev<1>::right &&
1094 i == (n - 1))
1095 weights[i] /= 2.;
1096 }
1097
1098 return weights;
1099 }
1101} // namespace internal
1102
1103
1104template <>
1107 , ep(ep)
1108{
1109 Assert(n > 0, ExcMessage("Need at least one point for quadrature rules"));
1110 std::vector<double> p =
1112 std::vector<double> w =
1114
1115 for (unsigned int i = 0; i < this->size(); ++i)
1116 {
1118 this->weights[i] = w[i];
1119 }
1120}
1121
1122
1123template <int dim>
1125 EndPoint ep)
1127 n,
1129 , ep(ep)
1130{}
1131
1132
1133
1134namespace internal
1135{
1137 {
1138 // Computes the points of the quadrature formula.
1139 std::vector<double>
1141 {
1142 std::vector<double> points(n);
1143 // n point quadrature: index from 0 to n-1
1144 for (unsigned short i = 0; i < n; ++i)
1145 // would be cos(i Pi/N)
1146 // put + Pi so we start from the smallest point
1147 // then map from [-1,1] to [0,1]
1148 points[i] =
1149 1. / 2. *
1150 (1. + std::cos(numbers::PI * (1 + double(i) / double(n - 1))));
1151
1152 return points;
1153 }
1154
1155 // Computes the weights of the quadrature formula.
1156 std::vector<double>
1158 {
1159 std::vector<double> weights(n);
1160
1161 for (unsigned short i = 0; i < n; ++i)
1162 {
1163 // same weights as on [-1,1]
1164 weights[i] = numbers::PI / double((n - 1));
1165 if (i == 0 || i == (n - 1))
1166 weights[i] /= 2.;
1167 }
1168
1169 return weights;
1170 }
1171 } // namespace QGaussLobattoChebyshev
1172} // namespace internal
1173
1174
1175
1176template <>
1179{
1180 Assert(n > 1,
1181 ExcMessage(
1182 "Need at least two points for Gauss-Lobatto quadrature rule"));
1183 std::vector<double> p =
1185 std::vector<double> w =
1187
1188 for (unsigned int i = 0; i < this->size(); ++i)
1189 {
1191 this->weights[i] = w[i];
1192 }
1193}
1194
1195
1196template <int dim>
1199{}
1200
1201
1202
1203template <int dim>
1205{
1206 std::vector<Point<dim>> qpoints;
1207 std::vector<double> weights;
1208
1209 for (unsigned int i = 0; i < quad.size(); ++i)
1210 {
1211 double r = 0;
1212 /* Use "int d" instead of the more natural "unsigned int d" to work
1213 * around a wrong diagnostic in gcc-10.3.0 that warns about that the
1214 * comparison "d < dim" is always false in case of "dim == 0".
1215 * MM 2021 */
1216 for (int d = 0; d < dim; ++d)
1218 if (r <= 1 + 1e-10)
1219 {
1222 }
1223 }
1224}
1225
1226
1227
1228template <int dim>
1229template <int spacedim>
1232 const std::array<Point<spacedim>, dim + 1> &vertices) const
1233{
1234 Assert(dim <= spacedim,
1235 ExcMessage("Invalid combination of dim and spacedim ."));
1237 for (unsigned int d = 0; d < dim; ++d)
1238 Bt[d] = vertices[d + 1] - vertices[0];
1239
1240 const auto B = Bt.transpose();
1241 const double J = std::abs(B.determinant());
1242
1243 // if the determinant is zero, we return an empty quadrature
1244 if (J < 1e-12)
1246
1247 std::vector<Point<spacedim>> qp(this->size());
1248 std::vector<double> w(this->size());
1249
1250 for (unsigned int i = 0; i < this->size(); ++i)
1251 {
1252 qp[i] =
1253 Point<spacedim>(vertices[0] + apply_transformation(B, this->point(i)));
1254 w[i] = J * this->weight(i);
1255 }
1256
1258}
1259
1260
1261
1262template <int dim>
1263template <int spacedim>
1266 const std::vector<std::array<Point<spacedim>, dim + 1>> &simplices) const
1267{
1268 Assert(!(dim == 1 && spacedim == 1),
1269 ExcMessage("This function is not supposed to work in 1D-1d case."));
1270 Assert(dim <= spacedim,
1271 ExcMessage("Invalid combination of dim and spacedim ."));
1272
1273 std::vector<Point<spacedim>> qp;
1274 std::vector<double> ws;
1275 for (const auto &simplex : simplices)
1276 {
1277 const auto rule = this->compute_affine_transformation(simplex);
1278 std::transform(rule.get_points().begin(),
1279 rule.get_points().end(),
1280 std::back_inserter(qp),
1281 [&](const Point<spacedim> &p) { return p; });
1282 std::transform(rule.get_weights().begin(),
1283 rule.get_weights().end(),
1284 std::back_inserter(ws),
1285 [&](const double w) { return w; });
1286 }
1288}
1289
1290
1291
1295{
1298 this->weights.resize(base.size());
1299 for (unsigned int i = 0; i < base.size(); ++i)
1300 {
1301 const auto &q = base.point(i);
1302 const auto w = base.weight(i);
1303
1304 const auto xhat = q[0];
1305 const auto yhat = q[1];
1306
1307 const double t = numbers::PI_2 * yhat;
1308 const double pi = numbers::PI;
1309 const double st = std::sin(t);
1310 const double ct = std::cos(t);
1311 const double r = xhat / (st + ct);
1312
1313 const double J = pi * xhat / (2 * (std::sin(pi * yhat) + 1));
1314
1315 this->quadrature_points[i] = Point<2>(r * ct, r * st);
1316 this->weights[i] = w * J;
1317 }
1318}
1319
1320
1321
1323 : QTrianglePolar(QGauss<1>(n), QGauss<1>(n))
1324{}
1325
1326
1327
1330 const double beta)
1332{
1335 this->weights.resize(base.size());
1336 for (unsigned int i = 0; i < base.size(); ++i)
1337 {
1338 const auto &q = base.point(i);
1339 const auto w = base.weight(i);
1340
1341 const auto xhat = q[0];
1342 const auto yhat = q[1];
1343
1344 const double x = std::pow(xhat, beta) * (1 - yhat);
1345 const double y = std::pow(xhat, beta) * yhat;
1346
1347 const double J = beta * std::pow(xhat, 2. * beta - 1.);
1348
1350 this->weights[i] = w * J;
1351 }
1352}
1353
1354
1355
1356QDuffy::QDuffy(const unsigned int n, const double beta)
1357 : QDuffy(QGauss<1>(n), QGauss<1>(n), beta)
1358{}
1359
1360
1361
1362template <int dim>
1363QSplit<dim>::QSplit(const QSimplex<dim> &base, const Point<dim> &split_point)
1364{
1366 ExcMessage(
1367 "The split point should be inside the unit reference cell."));
1368
1369 std::array<Point<dim>, dim + 1> vertices;
1370 vertices[0] = split_point;
1371
1372 // Make a simplex from the split_point and the first dim vertices of each
1373 // face. In dimension three, we need to split the face in two triangles, so
1374 // we use once the first dim vertices of each face, and the second time the
1375 // the dim vertices of each face starting from 1.
1376 for (auto f : GeometryInfo<dim>::face_indices())
1377 for (unsigned int start = 0; start < (dim > 2 ? 2 : 1); ++start)
1378 {
1379 for (unsigned int i = 0; i < dim; ++i)
1382 const auto quad = base.compute_affine_transformation(vertices);
1384 {
1388 this->weights.insert(this->weights.end(),
1391 }
1392 }
1393}
1394
1395
1396
1397template <int dim>
1398QGaussSimplex<dim>::QGaussSimplex(const unsigned int n_points_1D)
1400{
1402 if (dim == 0 || dim == 1)
1403 {
1405
1408 }
1409 else if (dim == 2)
1410 {
1411 if (n_points_1D == 1)
1412 {
1413 const double p = 1.0 / 3.0;
1415 this->weights.emplace_back(0.5);
1416 }
1417 else if (n_points_1D == 2)
1418 {
1419 // The Hillion 7 scheme, as communicated by quadpy
1420 //
1421 // See: Numerical Integration on a Triangle, International Journal for
1422 // Numerical Methods in Engineering, 1977
1423 const double Q12 = 1.0 / 2.0;
1425 0.1550510257216822);
1427 0.6449489742783178);
1429 0.1550510257216822);
1431 0.6449489742783178);
1432
1433 this->weights.emplace_back(0.31804138174397717 * Q12);
1434 this->weights.emplace_back(0.18195861825602283 * Q12);
1435 this->weights.emplace_back(0.31804138174397717 * Q12);
1436 this->weights.emplace_back(0.18195861825602283 * Q12);
1437 }
1438 else if (n_points_1D == 3)
1439 {
1440 // The Hammer-Marlowe-Stroud 5 Scheme, as communicated by quadpy
1441 const double p0 = 2.0 / 7.0 - std::sqrt(15.0) / 21.0;
1442 const double p1 = 2.0 / 7.0 + std::sqrt(15.0) / 21.0;
1443 const double p2 = 3.0 / 7.0 - 2.0 * std::sqrt(15.0) / 21.0;
1444 const double p3 = 3.0 / 7.0 + 2.0 * std::sqrt(15.0) / 21.0;
1445 this->quadrature_points.emplace_back(1.0 / 3.0, 1.0 / 3.0);
1452
1453 const double q12 = 0.5;
1454 const double w0 = 9.0 / 40.0;
1455 const double w1 = 31.0 / 240.0 - std::sqrt(15.0) / 1200.0;
1456 const double w2 = 31.0 / 240.0 + std::sqrt(15.0) / 1200.0;
1457 this->weights.emplace_back(q12 * w0);
1458 this->weights.emplace_back(q12 * w1);
1459 this->weights.emplace_back(q12 * w1);
1460 this->weights.emplace_back(q12 * w1);
1461 this->weights.emplace_back(q12 * w2);
1462 this->weights.emplace_back(q12 * w2);
1463 this->weights.emplace_back(q12 * w2);
1464 }
1465 else if (n_points_1D == 4)
1466 {
1468 QWitherdenVincentSimplex<dim>(n_points_1D));
1469 }
1470 }
1471 else if (dim == 3)
1472 {
1473 if (n_points_1D == 1)
1474 {
1475 const double Q14 = 1.0 / 4.0;
1476 const double Q16 = 1.0 / 6.0;
1477
1479 this->weights.emplace_back(Q16);
1480 }
1481 // The Xiao Gimbutas 03 scheme, as communicated by quadpy
1482 //
1483 // See: A numerical algorithm for the construction of efficient quadrature
1484 // rules in two and higher dimensions, Computers & Mathematics with
1485 // Applications, 2010
1486 else if (n_points_1D == 2)
1487 {
1488 const double Q16 = 1.0 / 6.0;
1489 this->weights.emplace_back(0.1223220027573451 * Q16);
1490 this->weights.emplace_back(0.1280664127107469 * Q16);
1491 this->weights.emplace_back(0.1325680271444452 * Q16);
1492 this->weights.emplace_back(0.1406244096604032 * Q16);
1493 this->weights.emplace_back(0.2244151669175574 * Q16);
1494 this->weights.emplace_back(0.2520039808095023 * Q16);
1495
1497 0.1838503504920977,
1498 0.01271836631368145);
1500 0.2815238021235462,
1501 0.3621268299455338);
1503 0.01140332944455717,
1504 0.3586207204668839);
1506 0.1528181430909273,
1507 0.6384932999617267);
1509 0.5856628056552158,
1510 0.1308471689520965);
1512 0.1469183900871696,
1513 0.1403728057942107);
1514 }
1515 // Past this point the best rules (positive weights, minimal number of
1516 // points) we have right now are the Witherden-Vincent ones
1517 else if (n_points_1D == 3)
1518 {
1520 QWitherdenVincentSimplex<dim>(n_points_1D));
1521 }
1522 else if (n_points_1D == 4)
1523 {
1525 QWitherdenVincentSimplex<dim>(n_points_1D));
1526 }
1527 }
1528
1532 "QGaussSimplex is currently only implemented for "
1533 "n_points_1D = 1, 2, 3, and 4 while you are asking for "
1534 "n_points_1D = " +
1535 Utilities::to_string(n_points_1D)));
1536}
1537
1538namespace
1539{
1540 template <std::size_t b_dim>
1541 std::vector<std::array<double, b_dim>>
1542 all_permutations(const std::array<double, b_dim> &b_point)
1543 {
1544 std::vector<std::array<double, b_dim>> output;
1545
1546 // We want all possible permutations of the barycentric coordinates.
1547 // The easiest way to get all of them is to sort the input first and
1548 // then use next_permutation to cycle through them all.
1549 std::array<double, b_dim> temp = b_point;
1550 std::sort(temp.begin(), temp.end());
1551 do
1552 {
1553 output.push_back(temp);
1554 }
1555 while (std::next_permutation(temp.begin(), temp.end()));
1556
1557 return output;
1558 }
1559} // namespace
1560
1561
1562
1563template <int dim>
1565 const unsigned int n_points_1D,
1566 const bool use_odd_order)
1568{
1569 Assert(1 <= dim && dim <= 3, ExcNotImplemented());
1570 // Just use Gauss in 1d: this is a high-order open rule so this is a
1571 // reasonable equivalent for generic programming.
1572 if (dim == 1)
1573 {
1575 return;
1576 }
1577
1578 std::array<double, dim + 1> centroid;
1579 std::fill(centroid.begin(), centroid.end(), 1.0 / (dim + 1.0));
1580 std::vector<std::vector<std::array<double, dim + 1>>> b_point_permutations;
1581 std::vector<double> b_weights;
1582
1583 // We can simplify the implementation of these quadrature rules
1584 // by quite a bit by exploiting symmetry - we do essentially the
1585 // same thing for each barycentric coordinate, so we can express
1586 // our quadrature rule as permutations of barycentric points
1587 // instead of writing things out explicitly.
1588
1589 // Apply a Barycentric permutation where one point is different.
1590 // Equivalent to d3_aa and s31 in quadpy.
1591 auto process_point_1 = [&](const double a, const double w) {
1592 const double b = 1.0 - dim * a;
1593 std::array<double, dim + 1> b_point;
1594 std::fill(b_point.begin(), b_point.begin() + dim, a);
1595 b_point[dim] = b;
1596
1597 b_weights.push_back(w);
1598 b_point_permutations.push_back(all_permutations(b_point));
1599 };
1600
1601 // Apply a Barycentric permutation where two points (in 3d) are different.
1602 // Equivalent to s22 in quadpy.
1603 auto process_point_2 = [&](const double a, const double w) {
1604 Assert(dim == 3, ExcInternalError());
1605 const double b = (1.0 - 2.0 * a) / 2.0;
1606 std::array<double, dim + 1> b_point;
1607 std::fill(b_point.begin(), b_point.begin() + dim - 1, a);
1608 b_point[dim - 1] = b;
1609 b_point[dim] = b;
1610
1611 b_weights.push_back(w);
1612 b_point_permutations.push_back(all_permutations(b_point));
1613 };
1614
1615 // Apply a Barycentric permutation where three (or four) points
1616 // are different (since there are two inputs).
1617 // Equivalent to d3_ab and s211 in quadpy.
1618 auto process_point_3 = [&](const double a, const double b, const double w) {
1619 const double c = 1.0 - (dim - 1.0) * a - b;
1620 std::array<double, dim + 1> b_point;
1621 std::fill(b_point.begin(), b_point.begin() + dim - 1, a);
1622 b_point[dim - 1] = b;
1623 b_point[dim] = c;
1624
1625 b_weights.push_back(w);
1626 b_point_permutations.push_back(all_permutations(b_point));
1627 };
1628
1629 switch (n_points_1D)
1630 {
1631 case 1:
1632 switch (dim)
1633 {
1634 case 2:
1635 if (use_odd_order)
1636 {
1637 // WV-1, 2d
1638 b_point_permutations.push_back({centroid});
1639 b_weights.push_back(1.0000000000000000e+00);
1640 }
1641 else
1642 {
1643 // WV-2, 2d
1644 process_point_1(1.6666666666666669e-01,
1645 3.3333333333333331e-01);
1646 }
1647 break;
1648 case 3:
1649 if (use_odd_order)
1650 {
1651 // WV-1, 3d
1652 b_point_permutations.push_back({centroid});
1653 b_weights.push_back(1.0000000000000000e+00);
1654 }
1655 else
1656 {
1657 // WV-2, 3d
1658 process_point_1(1.3819660112501050e-01,
1659 2.5000000000000000e-01);
1660 }
1661 break;
1662 default:
1663 Assert(false, ExcNotImplemented());
1664 }
1665 break;
1666 case 2:
1667 switch (dim)
1668 {
1669 case 2:
1670 // WV-4 in both cases (no WV-3 in 2d)
1671 process_point_1(9.1576213509770743e-02, 1.0995174365532187e-01);
1672 process_point_1(4.4594849091596489e-01, 2.2338158967801147e-01);
1673 break;
1674 case 3:
1675 if (use_odd_order)
1676 {
1677 // WV-3, 3d
1678 process_point_1(3.2816330251638171e-01,
1679 1.3621784253708741e-01);
1680 process_point_1(1.0804724989842859e-01,
1681 1.1378215746291261e-01);
1682 }
1683 else
1684 {
1685 // WV-5 (no WV-4 in 3d)
1687 }
1688 break;
1689 default:
1690 Assert(false, ExcInternalError());
1691 }
1692 break;
1693 case 3:
1694 switch (dim)
1695 {
1696 case 2:
1697 if (use_odd_order)
1698 {
1699 // WV-5, 2d
1700 b_point_permutations.push_back({centroid});
1701 b_weights.push_back(2.2500000000000001e-01);
1702 process_point_1(1.0128650732345634e-01,
1703 1.2593918054482714e-01);
1704 process_point_1(4.7014206410511511e-01,
1705 1.3239415278850619e-01);
1706 }
1707 else
1708 {
1709 // WV-6, 2d
1710 process_point_1(6.3089014491502227e-02,
1711 5.0844906370206819e-02);
1712 process_point_1(2.4928674517091043e-01,
1713 1.1678627572637937e-01);
1714 process_point_3(5.3145049844816938e-02,
1715 3.1035245103378439e-01,
1716 8.2851075618373571e-02);
1717 }
1718 break;
1719 case 3:
1720 if (use_odd_order)
1721 {
1722 // WV-5, 3d
1723 process_point_1(3.1088591926330061e-01,
1724 1.1268792571801590e-01);
1725 process_point_1(9.2735250310891248e-02,
1726 7.3493043116361956e-02);
1727 process_point_2(4.5503704125649642e-02,
1728 4.2546020777081472e-02);
1729 }
1730 else
1731 {
1732 // WV-6, 3d
1733 process_point_1(4.0673958534611372e-02,
1734 1.0077211055320640e-02);
1735 process_point_1(3.2233789014227548e-01,
1736 5.5357181543654717e-02);
1737 process_point_1(2.1460287125915201e-01,
1738 3.9922750258167487e-02);
1739 process_point_3(6.3661001875017442e-02,
1740 6.0300566479164919e-01,
1741 4.8214285714285710e-02);
1742 }
1743 break;
1744 default:
1745 Assert(false, ExcInternalError());
1746 }
1747 break;
1748 case 4:
1749 switch (dim)
1750 {
1751 case 2:
1752 if (use_odd_order)
1753 {
1754 // WV-7, 2d
1755 process_point_1(3.3730648554587850e-02,
1756 1.6545050110792131e-02);
1757 process_point_1(4.7430969250471822e-01,
1758 7.7086646185986069e-02);
1759 process_point_1(2.4157738259540357e-01,
1760 1.2794417123015558e-01);
1761 process_point_3(4.7036644652595216e-02,
1762 1.9868331479735168e-01,
1763 5.5878732903199779e-02);
1764 }
1765 else
1766 {
1767 // WV-8, 2d
1768 b_point_permutations.push_back({centroid});
1769 b_weights.push_back(1.4431560767778717e-01);
1770 process_point_1(5.0547228317030957e-02,
1771 3.2458497623198079e-02);
1772 process_point_1(4.5929258829272313e-01,
1773 9.5091634267284619e-02);
1774 process_point_1(1.7056930775176021e-01,
1775 1.0321737053471824e-01);
1776 process_point_3(8.3947774099575878e-03,
1777 2.6311282963463811e-01,
1778 2.7230314174434993e-02);
1779 }
1780 break;
1781 case 3:
1782 if (use_odd_order)
1783 {
1784 // WV-7, 3d
1785 b_point_permutations.push_back({centroid});
1786 b_weights.push_back(9.5485289464130846e-02);
1787 process_point_1(3.1570114977820279e-01,
1788 4.2329581209967028e-02);
1789 process_point_2(5.0489822598396350e-02,
1790 3.1896927832857580e-02);
1791 process_point_3(1.8883383102600099e-01,
1792 5.7517163758699996e-01,
1793 3.7207130728334620e-02);
1794 process_point_3(2.1265472541483140e-02,
1795 8.1083024109854862e-01,
1796 8.1107708299033420e-03);
1797 }
1798 else
1799 {
1800 // WV-8, 3d
1801 process_point_1(1.0795272496221089e-01,
1802 2.6426650908408830e-02);
1803 process_point_1(1.8510948778258660e-01,
1804 5.2031747563738531e-02);
1805 process_point_1(4.2316543684767283e-02,
1806 7.5252561535401989e-03);
1807 process_point_1(3.1418170912403898e-01,
1808 4.1763782856934897e-02);
1809 process_point_2(4.3559132858383021e-01,
1810 3.6280930261308818e-02);
1811 process_point_3(2.1433930127130570e-02,
1812 7.1746406342630831e-01,
1813 7.1569028908444327e-03);
1814 process_point_3(2.0413933387602909e-01,
1815 5.8379737830214440e-01,
1816 1.5453486150960340e-02);
1817 }
1818 break;
1819 default:
1820 Assert(false, ExcInternalError());
1821 }
1822 break;
1823 case 5:
1824 switch (dim)
1825 {
1826 case 2:
1827 if (use_odd_order)
1828 {
1829 // WV-9, 2d
1830 b_point_permutations.push_back({centroid});
1831 b_weights.push_back(9.7135796282798836e-02);
1832 process_point_1(4.4729513394452691e-02,
1833 2.5577675658698031e-02);
1834 process_point_1(4.8968251919873762e-01,
1835 3.1334700227139071e-02);
1836 process_point_1(4.3708959149293664e-01,
1837 7.7827541004774278e-02);
1838 process_point_1(1.8820353561903275e-01,
1839 7.9647738927210249e-02);
1840 process_point_3(3.6838412054736258e-02,
1841 2.2196298916076568e-01,
1842 4.3283539377289376e-02);
1843 }
1844 else
1845 {
1846 // WV-10, 2d
1847 b_point_permutations.push_back({centroid});
1848 b_weights.push_back(8.1743329146285973e-02);
1849 process_point_1(3.2055373216943517e-02,
1850 1.3352968813149567e-02);
1851 process_point_1(1.4216110105656438e-01,
1852 4.5957963604744731e-02);
1853 process_point_3(2.8367665339938453e-02,
1854 1.6370173373718250e-01,
1855 2.5297757707288385e-02);
1856 process_point_3(2.9619889488729734e-02,
1857 3.6914678182781102e-01,
1858 3.4184648162959429e-02);
1859 process_point_3(1.4813288578382056e-01,
1860 3.2181299528883545e-01,
1861 6.3904906396424044e-02);
1862 }
1863 break;
1864 case 3:
1865 if (use_odd_order)
1866 {
1867 // WV-9, 3d
1868 b_point_permutations.push_back({centroid});
1869 b_weights.push_back(5.8010548912480253e-02);
1870 process_point_1(6.1981697552226933e-10,
1871 6.4319281759256394e-05);
1872 process_point_1(1.6077453539526160e-01,
1873 2.3173338462425461e-02);
1874 process_point_1(3.2227652182142102e-01,
1875 2.9562912335429289e-02);
1876 process_point_1(4.5108918345413578e-02,
1877 8.0639799796161822e-03);
1878 process_point_2(1.1229654600437609e-01,
1879 3.8134080103702457e-02);
1880 process_point_3(4.5887144875245922e-01,
1881 2.5545792330413102e-03,
1882 8.3844221982985519e-03);
1883 process_point_3(3.3775870685338598e-02,
1884 7.1835032644207453e-01,
1885 1.0234559352745330e-02);
1886 process_point_3(1.8364136980992790e-01,
1887 3.4415910578175279e-02,
1888 2.0524915967988139e-02);
1889 }
1890 else
1891 {
1892 // WV-10, 3d
1893 b_point_permutations.push_back({centroid});
1894 b_weights.push_back(4.7399773556020743e-02);
1895 process_point_1(3.1225006869518868e-01,
1896 2.6937059992268701e-02);
1897 process_point_1(1.1430965385734609e-01,
1898 9.8691597167933822e-03);
1899 process_point_3(4.1043073921896539e-01,
1900 1.6548602561961109e-01,
1901 1.1393881220195230e-02);
1902 process_point_3(6.1380088247906528e-03,
1903 9.4298876734520487e-01,
1904 3.6194434433925362e-04);
1905 process_point_3(1.2105018114558939e-01,
1906 4.7719037990428043e-01,
1907 2.5739731980456069e-02);
1908 process_point_3(3.2779468216442620e-02,
1909 5.9425626948000698e-01,
1910 1.0135871679755789e-02);
1911 process_point_3(3.2485281564823047e-02,
1912 8.0117728465834437e-01,
1913 6.5761472770359038e-03);
1914 process_point_3(1.7497934218393901e-01,
1915 6.2807184547536599e-01,
1916 1.2907035798861989e-02);
1917 }
1918 break;
1919 default:
1920 Assert(false, ExcNotImplemented());
1921 }
1922 break;
1923 case 6:
1924 // There is no WV-11 rule in 3d yet
1925 Assert(dim == 2, ExcNotImplemented());
1926 if (use_odd_order)
1927 {
1928 // WV-11, 2d
1929 b_point_permutations.push_back({centroid});
1930 b_weights.push_back(8.5761179732224219e-02);
1931 process_point_1(2.8485417614371900e-02, 1.0431870512894697e-02);
1932 process_point_1(4.9589190096589092e-01, 1.6606273054585369e-02);
1933 process_point_1(1.0263548271224643e-01, 3.8630759237019321e-02);
1934 process_point_1(4.3846592676435220e-01, 6.7316154079468296e-02);
1935 process_point_1(2.1021995670317828e-01, 7.0515684111716576e-02);
1936 process_point_3(7.3254276860644785e-03,
1937 1.4932478865208237e-01,
1938 1.0290289572953278e-02);
1939 process_point_3(4.6010500165429957e-02,
1940 2.8958112563770588e-01,
1941 4.0332476640500554e-02);
1942 }
1943 else
1944 {
1945 // WV-12, 2d
1946 process_point_1(2.4646363436335583e-02, 7.9316425099736389e-03);
1947 process_point_1(4.8820375094554153e-01, 2.4266838081452032e-02);
1948 process_point_1(1.0925782765935427e-01, 2.8486052068877544e-02);
1949 process_point_1(4.4011164865859309e-01, 4.9918334928060942e-02);
1950 process_point_1(2.7146250701492608e-01, 6.2541213195902765e-02);
1951 process_point_3(2.1382490256170616e-02,
1952 1.2727971723358933e-01,
1953 1.5083677576511438e-02);
1954 process_point_3(2.3034156355267121e-02,
1955 2.9165567973834094e-01,
1956 2.1783585038607559e-02);
1957 process_point_3(1.1629601967792658e-01,
1958 2.5545422863851736e-01,
1959 4.3227363659414209e-02);
1960 }
1961 break;
1962 case 7:
1963 // There is no WV-13 rule in 3d yet
1964 Assert(dim == 2, ExcNotImplemented());
1965 if (use_odd_order)
1966 {
1967 // WV-13, 2d
1968 b_point_permutations.push_back({centroid});
1969 b_weights.push_back(6.7960036586831640e-02);
1970 process_point_1(2.1509681108843159e-02, 6.0523371035391717e-03);
1971 process_point_1(4.8907694645253935e-01, 2.3994401928894731e-02);
1972 process_point_1(4.2694141425980042e-01, 5.5601967530453329e-02);
1973 process_point_1(2.2137228629183292e-01, 5.8278485119199981e-02);
1974 process_point_3(5.1263891023823893e-03,
1975 2.7251581777342970e-01,
1976 9.5906810035432631e-03);
1977 process_point_3(2.4370186901093827e-02,
1978 1.1092204280346341e-01,
1979 1.4965401105165668e-02);
1980 process_point_3(8.7895483032197297e-02,
1981 1.6359740106785048e-01,
1982 2.4179039811593819e-02);
1983 process_point_3(6.8012243554206653e-02,
1984 3.0844176089211778e-01,
1985 3.4641276140848373e-02);
1986 }
1987 else
1988 {
1989 // WV-14, 2d
1990 process_point_1(1.9390961248701044e-02, 4.9234036024000819e-03);
1991 process_point_1(6.1799883090872587e-02, 1.4433699669776668e-02);
1992 process_point_1(4.8896391036217862e-01, 2.1883581369428889e-02);
1993 process_point_1(4.1764471934045394e-01, 3.2788353544125348e-02);
1994 process_point_1(1.7720553241254344e-01, 4.2162588736993016e-02);
1995 process_point_1(2.7347752830883865e-01, 5.1774104507291585e-02);
1996 process_point_3(1.2683309328720416e-03,
1997 1.1897449769695684e-01,
1998 5.0102288385006719e-03);
1999 process_point_3(1.4646950055654417e-02,
2000 2.9837288213625779e-01,
2001 1.4436308113533840e-02);
2002 process_point_3(5.7124757403647919e-02,
2003 1.7226668782135557e-01,
2004 2.4665753212563674e-02);
2005 process_point_3(9.2916249356971847e-02,
2006 3.3686145979634496e-01,
2007 3.8571510787060684e-02);
2008 }
2009 break;
2010 default:
2011 Assert(false, ExcNotImplemented());
2012 }
2013
2014 Assert(b_point_permutations.size() == b_weights.size(), ExcInternalError());
2015 for (unsigned int permutation_n = 0; permutation_n < b_weights.size();
2016 ++permutation_n)
2017 {
2018 for (const std::array<double, dim + 1> &b_point :
2019 b_point_permutations[permutation_n])
2020 {
2021 const double volume = (dim == 2 ? 1.0 / 2.0 : 1.0 / 6.0);
2022 this->weights.emplace_back(volume * b_weights[permutation_n]);
2023 Point<dim> c_point;
2024 for (int d = 0; d < dim; ++d)
2025 c_point[d] = b_point[d];
2027 }
2028 }
2029}
2030
2031
2032
2033namespace
2034{
2035 template <int dim>
2037 setup_qiterated_1D(const Quadrature<dim> &, const unsigned int)
2038 {
2039 Assert(false, ExcInternalError());
2041 }
2042
2043
2044
2047 const unsigned int n_copies)
2048 {
2050 }
2051} // namespace
2052
2053
2054
2055template <int dim>
2057 const unsigned int n_copies)
2058{
2059 switch (dim)
2060 {
2061 case 1:
2064 break;
2065 case 2:
2066 case 3:
2067 {
2068 const auto n_refinements =
2069 static_cast<unsigned int>(std::round(std::log2(n_copies)));
2070 Assert((1u << n_refinements) == n_copies,
2071 ExcMessage("The number of copies must be a power of 2."));
2073 const auto reference_cell = ReferenceCells::get_simplex<dim>();
2074 GridGenerator::reference_cell(tria, reference_cell);
2075 tria.refine_global(n_refinements);
2076 const Mapping<dim> &mapping =
2077 reference_cell.template get_default_linear_mapping<dim>();
2078 FE_Nothing<dim> fe(reference_cell);
2079
2080 FEValues<dim> fe_values(mapping,
2081 fe,
2084 std::vector<Point<dim>> points;
2085 std::vector<double> weights;
2086 for (const auto &cell : tria.active_cell_iterators())
2087 {
2088 fe_values.reinit(cell);
2089 for (unsigned int qp = 0; qp < base_quad.size(); ++qp)
2090 {
2092 weights.push_back(fe_values.JxW(qp));
2093 }
2094 }
2095
2098
2099 break;
2100 }
2101 default:
2102 Assert(false, ExcNotImplemented());
2103 }
2104}
2105
2106
2107
2108template <int dim>
2109QGaussWedge<dim>::QGaussWedge(const unsigned int n_points)
2111{
2112 AssertDimension(dim, 3);
2113
2116
2117 for (unsigned int i = 0; i < quad_line.size(); ++i)
2118 for (unsigned int j = 0; j < quad_tri.size(); ++j)
2119 {
2124 }
2125
2129}
2130
2131
2132
2133template <int dim>
2134QGaussPyramid<dim>::QGaussPyramid(const unsigned int n_points_1D)
2136{
2137 AssertDimension(dim, 3);
2138
2139 if (n_points_1D == 1)
2140 {
2141 const double Q14 = 1.0 / 4.0;
2142 const double Q43 = 4.0 / 3.0;
2143
2145 this->weights.emplace_back(Q43);
2146 }
2147 else if (n_points_1D == 2)
2148 {
2149 // clang-format off
2158 // clang-format on
2159
2160 this->weights.emplace_back(0.10078588207983);
2161 this->weights.emplace_back(0.23254745125351);
2162 this->weights.emplace_back(0.10078588207983);
2163 this->weights.emplace_back(0.23254745125351);
2164 this->weights.emplace_back(0.10078588207983);
2165 this->weights.emplace_back(0.23254745125351);
2166 this->weights.emplace_back(0.10078588207983);
2167 this->weights.emplace_back(0.23254745125351);
2168 }
2169
2173}
2174
2175
2176
2177// explicit specialization
2178// note that 1d formulae are specialized by implementation above
2179template class QGauss<2>;
2180template class QGaussLobatto<2>;
2181template class QMidpoint<2>;
2182template class QTrapezoid<2>;
2183template class QSimpson<2>;
2184template class QMilne<2>;
2185template class QWeddle<2>;
2186
2187template class QGauss<3>;
2188template class QGaussLobatto<3>;
2189template class QMidpoint<3>;
2190template class QTrapezoid<3>;
2191template class QSimpson<3>;
2192template class QMilne<3>;
2193template class QWeddle<3>;
2194
2195template class QSorted<1>;
2196template class QSorted<2>;
2197template class QSorted<3>;
2198
2199template class QTelles<1>;
2200template class QTelles<2>;
2201template class QTelles<3>;
2202
2203template class QGaussChebyshev<1>;
2204template class QGaussChebyshev<2>;
2205template class QGaussChebyshev<3>;
2206
2210
2211template class QGaussLobattoChebyshev<1>;
2212template class QGaussLobattoChebyshev<2>;
2213template class QGaussLobattoChebyshev<3>;
2214
2215template class QSimplex<1>;
2216template class QSimplex<2>;
2217template class QSimplex<3>;
2218
2219template class QIteratedSimplex<1>;
2220template class QIteratedSimplex<2>;
2221template class QIteratedSimplex<3>;
2222
2223template class QSplit<1>;
2224template class QSplit<2>;
2225template class QSplit<3>;
2226
2227template class QGaussSimplex<0>;
2228template class QGaussSimplex<1>;
2229template class QGaussSimplex<2>;
2230template class QGaussSimplex<3>;
2231template class QGaussWedge<0>;
2232template class QGaussWedge<1>;
2233template class QGaussWedge<2>;
2234template class QGaussWedge<3>;
2235template class QGaussPyramid<0>;
2236template class QGaussPyramid<1>;
2237template class QGaussPyramid<2>;
2238template class QGaussPyramid<3>;
2239
2240template class QWitherdenVincentSimplex<1>;
2241template class QWitherdenVincentSimplex<2>;
2242template class QWitherdenVincentSimplex<3>;
2243
2244#ifndef DOXYGEN
2247 const std::array<Point<1>, 1 + 1> &vertices) const;
2248
2251 const std::array<Point<2>, 1 + 1> &vertices) const;
2252
2255 const std::array<Point<2>, 2 + 1> &vertices) const;
2256
2259 const std::array<Point<3>, 1 + 1> &vertices) const;
2260
2263 const std::array<Point<3>, 2 + 1> &vertices) const;
2264
2267 const std::array<Point<3>, 3 + 1> &vertices) const;
2268
2271 const std::vector<std::array<Point<2>, 1 + 1>> &simplices) const;
2272
2275 const std::vector<std::array<Point<3>, 1 + 1>> &simplices) const;
2276
2279 const std::vector<std::array<Point<2>, 2 + 1>> &simplices) const;
2280
2283 const std::vector<std::array<Point<3>, 2 + 1>> &simplices) const;
2284
2287 const std::vector<std::array<Point<3>, 3 + 1>> &simplices) const;
2288#endif
2289
DerivativeForm< 1, spacedim, dim, Number > transpose() const
const Point< spacedim > & quadrature_point(const unsigned int q_point) const
double JxW(const unsigned int q_point) const
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access > > &cell)
Abstract base class for mapping classes.
Definition mapping.h:317
Definition point.h:112
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)
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)
Quadrature< spacedim > compute_affine_transformation(const std::array< Point< spacedim >, dim+1 > &vertices) const
Quadrature< spacedim > mapped_quadrature(const std::vector< std::array< Point< spacedim >, dim+1 > > &simplices) const
bool compare_weights(const unsigned int a, const unsigned int b) const
QSplit(const QSimplex< dim > &base, const Point< dim > &split_point)
QWitherdenVincentSimplex(const unsigned int n_points_1D, const bool use_odd_order=true)
std::vector< Point< dim > > quadrature_points
const Point< dim > & point(const unsigned int i) const
bool is_tensor_product_flag
double weight(const unsigned int i) const
std::vector< double > weights
unsigned int size() const
void refine_global(const unsigned int times=1)
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
Point< 3 > vertices[4]
Tensor< 1, spacedim, typename ProductType< Number1, Number2 >::type > apply_transformation(const DerivativeForm< 1, dim, spacedim, Number1 > &grad_F, const Tensor< 1, dim, Number2 > &d_x)
IteratorRange< active_cell_iterator > active_cell_iterators() const
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
@ update_JxW_values
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
void quadrature_points(const Triangulation< dim, spacedim > &triangulation, const Quadrature< dim > &quadrature, const std::vector< std::vector< BoundingBox< spacedim > > > &global_bounding_boxes, ParticleHandler< dim, spacedim > &particle_handler, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()), const std::vector< std::vector< double > > &properties={})
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Number jacobi_polynomial_value(const unsigned int degree, const int alpha, const int beta, const Number x)
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:480
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)
static constexpr double PI_2
Definition numbers.h:264
static constexpr double PI
Definition numbers.h:259
::VectorizedArray< Number, width > log(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > tan(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
static std_cxx20::ranges::iota_view< unsigned int, unsigned int > face_indices()
static Point< dim > unit_cell_vertex(const unsigned int vertex)
const ::Triangulation< dim, spacedim > & tria