75 std::vector<long double> points =
76 Polynomials::jacobi_polynomial_roots<long double>(n, 0, 0);
78 for (
unsigned int i = 0; i < (points.size() + 1) / 2; ++i)
80 this->quadrature_points[i][0] = points[i];
81 this->quadrature_points[n - i - 1][0] = 1. - points[i];
84 const long double pp =
87 const long double x = -1. + 2. * points[i];
88 const double w = 1. / ((1. -
x *
x) *
pp *
pp);
90 this->weights[n - i - 1] = w;
105 long double result = n - 1;
106 for (
int i = n - 2; i > 1; --i)
120 std::vector<long double>
125 const unsigned int q =
x.size();
126 std::vector<long double> w(
q);
128 const long double factor =
131 for (
unsigned int i = 0; i <
q; ++i)
133 const long double s =
135 w[i] = factor / (s * s);
162 std::vector<double> q_points(n);
170 q_points[1] = 2. / 3.;
179 q_points[0] = 0.000000000000000000;
180 q_points[1] = 0.212340538239152943;
181 q_points[2] = 0.590533135559265343;
182 q_points[3] = 0.911412040487296071;
185 q_points[0] = 0.000000000000000000;
186 q_points[1] = 0.139759864343780571;
187 q_points[2] = 0.416409567631083166;
188 q_points[3] = 0.723156986361876197;
189 q_points[4] = 0.942895803885482331;
192 q_points[0] = 0.000000000000000000;
193 q_points[1] = 0.098535085798826416;
194 q_points[2] = 0.304535726646363913;
195 q_points[3] = 0.562025189752613841;
196 q_points[4] = 0.801986582126391845;
197 q_points[5] = 0.960190142948531222;
200 q_points[0] = 0.000000000000000000;
201 q_points[1] = 0.073054328680258851;
202 q_points[2] = 0.230766137969945495;
203 q_points[3] = 0.441328481228449865;
204 q_points[4] = 0.663015309718845702;
205 q_points[5] = 0.851921400331515644;
206 q_points[6] = 0.970683572840215114;
209 q_points[0] = 0.000000000000000000;
210 q_points[1] = 0.056262560536922135;
211 q_points[2] = 0.180240691736892389;
212 q_points[3] = 0.352624717113169672;
213 q_points[4] = 0.547153626330555420;
214 q_points[5] = 0.734210177215410598;
215 q_points[6] = 0.885320946839095790;
216 q_points[7] = 0.977520613561287499;
227 const ::QGaussRadau<1>::EndPoint end_point)
232 case ::QGaussRadau<1>::EndPoint::left:
234 case ::QGaussRadau<1>::EndPoint::right:
236 std::vector<double> points(n);
237 for (
unsigned int i = 0; i < n; ++i)
247 "This constructor can only be called with either "
248 "QGaussRadau::left or QGaussRadau::right as second argument."));
264 std::vector<double> weights(n);
275 weights[0] = 1. / 9.;
280 weights[0] = 0.062500000000000000;
281 weights[1] = 0.328844319980059696;
282 weights[2] = 0.388193468843171852;
283 weights[3] = 0.220462211176768369;
286 weights[0] = 0.040000000000000001;
287 weights[1] = 0.223103901083570894;
288 weights[2] = 0.311826522975741427;
289 weights[3] = 0.281356015149462124;
290 weights[4] = 0.143713560791225797;
293 weights[0] = 0.027777777777777776;
294 weights[1] = 0.159820376610255471;
295 weights[2] = 0.242693594234484888;
296 weights[3] = 0.260463391594787597;
297 weights[4] = 0.208450667155953895;
298 weights[5] = 0.100794192626740456;
301 weights[0] = 0.020408163265306121;
302 weights[1] = 0.119613744612656100;
303 weights[2] = 0.190474936822115581;
304 weights[3] = 0.223554914507283209;
305 weights[4] = 0.212351889502977870;
306 weights[5] = 0.159102115733650767;
307 weights[6] = 0.074494235556010341;
310 weights[0] = 0.015625000000000000;
311 weights[1] = 0.092679077401489660;
312 weights[2] = 0.152065310323392683;
313 weights[3] = 0.188258772694559262;
314 weights[4] = 0.195786083726246729;
315 weights[5] = 0.173507397817250691;
316 weights[6] = 0.124823950664932445;
317 weights[7] = 0.057254407372128648;
329 const ::QGaussRadau<1>::EndPoint end_point)
334 case ::QGaussRadau<1>::EndPoint::left:
336 case ::QGaussRadau<1>::EndPoint::right:
338 std::vector<double> weights(n);
339 for (
unsigned int i = 0; i < n; ++i)
348 "This constructor can only be called with either "
349 "QGaussRadau::EndPoint::left or "
350 "QGaussRadau::EndPoint::right as second argument."));
361 , end_point(end_point)
364 std::vector<double> p =
366 std::vector<double>
w =
369 for (
unsigned int i = 0; i < this->
size(); ++i)
372 this->weights[i] =
w[i];
385 std::vector<long double> points =
386 Polynomials::jacobi_polynomial_roots<long double>(n - 2, 1, 1);
387 points.insert(points.begin(), 0);
388 points.push_back(1.);
389 std::vector<long double>
w =
393 for (
unsigned int i = 0; i < points.size(); ++i)
396 this->weights[i] = 0.5 *
w[i];
406 this->quadrature_points[0] =
Point<1>(0.5);
407 this->weights[0] = 1.0;
416 static const double xpts[] = {0.0, 1.0};
417 static const double wts[] = {0.5, 0.5};
419 for (
unsigned int i = 0; i < this->
size(); ++i)
422 this->weights[i] =
wts[i];
432 static const double xpts[] = {0.0, 0.5, 1.0};
433 static const double wts[] = {1. / 6., 2. / 3., 1. / 6.};
435 for (
unsigned int i = 0; i < this->
size(); ++i)
438 this->weights[i] =
wts[i];
448 static const double xpts[] = {0.0, .25, .5, .75, 1.0};
449 static const double wts[] = {
450 7. / 90., 32. / 90., 12. / 90., 32. / 90., 7. / 90.};
452 for (
unsigned int i = 0; i < this->
size(); ++i)
455 this->weights[i] =
wts[i];
465 static const double xpts[] = {
466 0.0, 1. / 6., 1. / 3., .5, 2. / 3., 5. / 6., 1.0};
467 static const double wts[] = {41. / 840.,
475 for (
unsigned int i = 0; i < this->
size(); ++i)
478 this->weights[i] =
wts[i];
490 for (
unsigned int i = 0; i < this->
size(); ++i)
495 this->quadrature_points[i] =
497 this->weights[i] =
revert ? w[n - 1 - i] : w[i];
505 std::vector<double> q_points(n);
510 q_points[0] = 0.3333333333333333;
514 q_points[0] = 0.1120088061669761;
515 q_points[1] = 0.6022769081187381;
519 q_points[0] = 0.06389079308732544;
520 q_points[1] = 0.3689970637156184;
521 q_points[2] = 0.766880303938942;
525 q_points[0] = 0.04144848019938324;
526 q_points[1] = 0.2452749143206022;
527 q_points[2] = 0.5561654535602751;
528 q_points[3] = 0.848982394532986;
532 q_points[0] = 0.02913447215197205;
533 q_points[1] = 0.1739772133208974;
534 q_points[2] = 0.4117025202849029;
535 q_points[3] = 0.6773141745828183;
536 q_points[4] = 0.89477136103101;
540 q_points[0] = 0.02163400584411693;
541 q_points[1] = 0.1295833911549506;
542 q_points[2] = 0.3140204499147661;
543 q_points[3] = 0.5386572173517997;
544 q_points[4] = 0.7569153373774084;
545 q_points[5] = 0.922668851372116;
550 q_points[0] = 0.0167193554082585;
551 q_points[1] = 0.100185677915675;
552 q_points[2] = 0.2462942462079286;
553 q_points[3] = 0.4334634932570557;
554 q_points[4] = 0.6323509880476823;
555 q_points[5] = 0.81111862674023;
556 q_points[6] = 0.940848166743287;
560 q_points[0] = 0.01332024416089244;
561 q_points[1] = 0.07975042901389491;
562 q_points[2] = 0.1978710293261864;
563 q_points[3] = 0.354153994351925;
564 q_points[4] = 0.5294585752348643;
565 q_points[5] = 0.7018145299391673;
566 q_points[6] = 0.849379320441094;
567 q_points[7] = 0.953326450056343;
571 q_points[0] = 0.01086933608417545;
572 q_points[1] = 0.06498366633800794;
573 q_points[2] = 0.1622293980238825;
574 q_points[3] = 0.2937499039716641;
575 q_points[4] = 0.4466318819056009;
576 q_points[5] = 0.6054816627755208;
577 q_points[6] = 0.7541101371585467;
578 q_points[7] = 0.877265828834263;
579 q_points[8] = 0.96225055941096;
583 q_points[0] = 0.00904263096219963;
584 q_points[1] = 0.05397126622250072;
585 q_points[2] = 0.1353118246392511;
586 q_points[3] = 0.2470524162871565;
587 q_points[4] = 0.3802125396092744;
588 q_points[5] = 0.5237923179723384;
589 q_points[6] = 0.6657752055148032;
590 q_points[7] = 0.7941904160147613;
591 q_points[8] = 0.898161091216429;
592 q_points[9] = 0.9688479887196;
597 q_points[0] = 0.007643941174637681;
598 q_points[1] = 0.04554182825657903;
599 q_points[2] = 0.1145222974551244;
600 q_points[3] = 0.2103785812270227;
601 q_points[4] = 0.3266955532217897;
602 q_points[5] = 0.4554532469286375;
603 q_points[6] = 0.5876483563573721;
604 q_points[7] = 0.7139638500230458;
605 q_points[8] = 0.825453217777127;
606 q_points[9] = 0.914193921640008;
607 q_points[10] = 0.973860256264123;
611 q_points[0] = 0.006548722279080035;
612 q_points[1] = 0.03894680956045022;
613 q_points[2] = 0.0981502631060046;
614 q_points[3] = 0.1811385815906331;
615 q_points[4] = 0.2832200676673157;
616 q_points[5] = 0.398434435164983;
617 q_points[6] = 0.5199526267791299;
618 q_points[7] = 0.6405109167754819;
619 q_points[8] = 0.7528650118926111;
620 q_points[9] = 0.850240024421055;
621 q_points[10] = 0.926749682988251;
622 q_points[11] = 0.977756129778486;
638 std::vector<double> quadrature_weights(n);
643 quadrature_weights[0] = -1.0;
646 quadrature_weights[0] = -0.7185393190303845;
647 quadrature_weights[1] = -0.2814606809696154;
651 quadrature_weights[0] = -0.5134045522323634;
652 quadrature_weights[1] = -0.3919800412014877;
653 quadrature_weights[2] = -0.0946154065661483;
657 quadrature_weights[0] = -0.3834640681451353;
658 quadrature_weights[1] = -0.3868753177747627;
659 quadrature_weights[2] = -0.1904351269501432;
660 quadrature_weights[3] = -0.03922548712995894;
664 quadrature_weights[0] = -0.2978934717828955;
665 quadrature_weights[1] = -0.3497762265132236;
666 quadrature_weights[2] = -0.234488290044052;
667 quadrature_weights[3] = -0.0989304595166356;
668 quadrature_weights[4] = -0.01891155214319462;
672 quadrature_weights[0] = -0.2387636625785478;
673 quadrature_weights[1] = -0.3082865732739458;
674 quadrature_weights[2] = -0.2453174265632108;
675 quadrature_weights[3] = -0.1420087565664786;
676 quadrature_weights[4] = -0.05545462232488041;
677 quadrature_weights[5] = -0.01016895869293513;
682 quadrature_weights[0] = -0.1961693894252476;
683 quadrature_weights[1] = -0.2703026442472726;
684 quadrature_weights[2] = -0.239681873007687;
685 quadrature_weights[3] = -0.1657757748104267;
686 quadrature_weights[4] = -0.0889432271377365;
687 quadrature_weights[5] = -0.03319430435645653;
688 quadrature_weights[6] = -0.005932787015162054;
692 quadrature_weights[0] = -0.164416604728002;
693 quadrature_weights[1] = -0.2375256100233057;
694 quadrature_weights[2] = -0.2268419844319134;
695 quadrature_weights[3] = -0.1757540790060772;
696 quadrature_weights[4] = -0.1129240302467932;
697 quadrature_weights[5] = -0.05787221071771947;
698 quadrature_weights[6] = -0.02097907374214317;
699 quadrature_weights[7] = -0.003686407104036044;
703 quadrature_weights[0] = -0.1400684387481339;
704 quadrature_weights[1] = -0.2097722052010308;
705 quadrature_weights[2] = -0.211427149896601;
706 quadrature_weights[3] = -0.1771562339380667;
707 quadrature_weights[4] = -0.1277992280331758;
708 quadrature_weights[5] = -0.07847890261203835;
709 quadrature_weights[6] = -0.0390225049841783;
710 quadrature_weights[7] = -0.01386729555074604;
711 quadrature_weights[8] = -0.002408041036090773;
715 quadrature_weights[0] = -0.12095513195457;
716 quadrature_weights[1] = -0.1863635425640733;
717 quadrature_weights[2] = -0.1956608732777627;
718 quadrature_weights[3] = -0.1735771421828997;
719 quadrature_weights[4] = -0.135695672995467;
720 quadrature_weights[5] = -0.0936467585378491;
721 quadrature_weights[6] = -0.05578772735275126;
722 quadrature_weights[7] = -0.02715981089692378;
723 quadrature_weights[8] = -0.00951518260454442;
724 quadrature_weights[9] = -0.001638157633217673;
729 quadrature_weights[0] = -0.1056522560990997;
730 quadrature_weights[1] = -0.1665716806006314;
731 quadrature_weights[2] = -0.1805632182877528;
732 quadrature_weights[3] = -0.1672787367737502;
733 quadrature_weights[4] = -0.1386970574017174;
734 quadrature_weights[5] = -0.1038334333650771;
735 quadrature_weights[6] = -0.06953669788988512;
736 quadrature_weights[7] = -0.04054160079499477;
737 quadrature_weights[8] = -0.01943540249522013;
738 quadrature_weights[9] = -0.006737429326043388;
739 quadrature_weights[10] = -0.001152486965101561;
743 quadrature_weights[0] = -0.09319269144393;
744 quadrature_weights[1] = -0.1497518275763289;
745 quadrature_weights[2] = -0.166557454364573;
746 quadrature_weights[3] = -0.1596335594369941;
747 quadrature_weights[4] = -0.1384248318647479;
748 quadrature_weights[5] = -0.1100165706360573;
749 quadrature_weights[6] = -0.07996182177673273;
750 quadrature_weights[7] = -0.0524069547809709;
751 quadrature_weights[8] = -0.03007108900074863;
752 quadrature_weights[9] = -0.01424924540252916;
753 quadrature_weights[10] = -0.004899924710875609;
754 quadrature_weights[11] = -0.000834029009809656;
762 return quadrature_weights;
800 for (
unsigned int i = 0,
j =
ns_offset; i < n; ++i, ++
j)
810 this->quadrature_points[
j] = quad.point(i) *
fraction;
817 this->quadrature_points[i + n] =
822 this->quadrature_points[
j + n] =
824 this->weights[
j + n] =
829 for (
unsigned int i = 0; i <
size(); ++i)
832 this->quadrature_points[i] !=
origin,
834 "The singularity cannot be on a Gauss point of the same order!"));
839 "The quadrature formula you are using does not allow to "
840 "factor out the singularity, which is zero at one point."));
850 const double eps = 1e-8;
852 for (
unsigned int d = 0; d < 2; ++d)
878 std::vector<QGaussOneOverR<2>> quads;
895 unsigned int q_id = 0;
898 for (
unsigned int box = 0;
box < 4; ++
box)
904 for (
unsigned int q = 0;
q < quads[
box].size(); ++
q, ++
q_id)
907 this->quadrature_points[
q_id] =
917 const unsigned int vertex_index,
953 std::vector<Point<2>> &
ps = this->quadrature_points;
954 std::vector<double> &
ws = this->weights;
957 for (
unsigned int q = 0;
q <
gauss.size(); ++
q)
974 switch (vertex_index)
994 if (vertex_index != 0)
995 for (
unsigned int q = 0;
q <
size(); ++
q)
997 double x =
ps[
q][0] - .5,
y =
ps[
q][1] - .5;
1009 std::vector<unsigned int> permutation(quad.size());
1010 for (
unsigned int i = 0; i < quad.size(); ++i)
1013 std::sort(permutation.begin(),
1015 [
this](
const unsigned int x,
const unsigned int y) {
1016 return this->compare_weights(x, y);
1026 for (
unsigned int i = 0; i < quad.size(); ++i)
1028 this->weights[i] = quad.weight(permutation[i]);
1029 this->quadrature_points[i] = quad.point(permutation[i]);
1030 if (permutation[i] != i)
1031 this->is_tensor_product_flag =
false;
1040 return (this->weights[a] < this->weights[b]);
1058 , end_point(end_point)
1141 unsigned int cont = 0;
1142 const double tol = 1e-10;
1143 for (
unsigned int d = 0; d < quadrature_points.size(); ++d)
1161 for (
unsigned int d = 0; d < quadrature_points.size(); ++d)
1184 for (
unsigned int q = 0;
q < quadrature_points.size(); ++
q)
1186 double gamma = quadrature_points[
q][0] * 2 - 1;
1187 double eta = (Utilities::fixed_power<3>(gamma -
gamma_bar) +
1194 quadrature_points[
q][0] = (
eta + 1) / 2.0;
1195 weights[
q] = J * weights[
q];
1209 std::vector<double> points(n);
1211 for (
unsigned short i = 0; i < n; ++i)
1218 (1. +
double(2 * i + 1) /
double(2 * (n - 1) + 2))));
1231 std::vector<double> weights(n);
1233 for (
unsigned short i = 0; i < n; ++i)
1249 Assert(n > 0,
ExcMessage(
"Need at least one point for the quadrature rule"));
1253 for (
unsigned int i = 0; i < this->
size(); ++i)
1255 this->quadrature_points[i] =
Point<1>(p[i]);
1256 this->weights[i] = w[i];
1274 const unsigned int n,
1275 const ::QGaussRadauChebyshev<1>::EndPoint end_point)
1277 std::vector<double> points(n);
1279 for (
unsigned short i = 0; i < n; ++i)
1285 case ::QGaussRadauChebyshev<1>::EndPoint::left:
1291 (1 + 2 *
double(i) / (2 *
double(n - 1) + 1.))));
1295 case ::QGaussRadauChebyshev<1>::EndPoint::right:
1300 (2 *
double(n - 1) + 1.))));
1308 "This constructor can only be called with either "
1309 "QGaussRadauChebyshev::EndPoint::left or "
1310 "QGaussRadauChebyshev::EndPoint:right as second argument."));
1321 const unsigned int n,
1322 const ::QGaussRadauChebyshev<1>::EndPoint end_point)
1324 std::vector<double> weights(n);
1326 for (
unsigned short i = 0; i < n; ++i)
1329 weights[i] = 2. *
numbers::PI / double(2 * (n - 1) + 1.);
1333 else if (end_point ==
1349 , end_point(end_point)
1351 Assert(n > 0,
ExcMessage(
"Need at least one point for quadrature rules."));
1352 std::vector<double> points =
1357 for (
unsigned int i = 0; i < this->
size(); ++i)
1359 this->quadrature_points[i] =
Point<1>(points[i]);
1371 , end_point(end_point)
1384 std::vector<double> points(n);
1386 for (
unsigned short i = 0; i < n; ++i)
1401 std::vector<double> weights(n);
1403 for (
unsigned short i = 0; i < n; ++i)
1407 if (i == 0 || i == (n - 1))
1424 "Need at least two points for Gauss-Lobatto quadrature rule"));
1425 std::vector<double> p =
1427 std::vector<double> w =
1430 for (
unsigned int i = 0; i < this->
size(); ++i)
1432 this->quadrature_points[i] =
Point<1>(p[i]);
1433 this->weights[i] = w[i];
1448 std::vector<Point<dim>>
qpoints;
1449 std::vector<double> weights;
1451 for (
unsigned int i = 0; i < quad.size(); ++i)
1458 for (
int d = 0; d < dim; ++d)
1459 r += quad.point(i)[d];
1462 this->quadrature_points.push_back(quad.point(i));
1463 this->weights.push_back(quad.weight(i));
1471template <
int spacedim>
1477 ExcMessage(
"Invalid combination of dim and spacedim ."));
1479 for (
unsigned int d = 0; d < dim; ++d)
1480 Bt[d] = vertices[d + 1] - vertices[0];
1482 const auto B =
Bt.transpose();
1483 const double J =
std::abs(
B.determinant());
1489 std::vector<Point<spacedim>>
qp(this->
size());
1490 std::vector<double> w(this->
size());
1492 for (
unsigned int i = 0; i < this->
size(); ++i)
1496 w[i] = J * this->weight(i);
1505template <
int spacedim>
1510 Assert(!(dim == 1 && spacedim == 1),
1511 ExcMessage(
"This function is not supposed to work in 1D-1d case."));
1513 ExcMessage(
"Invalid combination of dim and spacedim ."));
1515 std::vector<Point<spacedim>>
qp;
1516 std::vector<double>
ws;
1519 const auto rule = this->compute_affine_transformation(simplex);
1520 std::transform(
rule.get_points().begin(),
1521 rule.get_points().end(),
1522 std::back_inserter(
qp),
1524 std::transform(
rule.get_weights().begin(),
1525 rule.get_weights().end(),
1526 std::back_inserter(
ws),
1527 [&](
const double w) { return w; });
1539 this->quadrature_points.resize(base.size());
1540 this->weights.resize(base.size());
1541 for (
unsigned int i = 0; i < base.size(); ++i)
1543 const auto &
q = base.point(i);
1544 const auto w = base.weight(i);
1546 const auto xhat =
q[0];
1547 const auto yhat =
q[1];
1553 const double r =
xhat / (st +
ct);
1558 this->weights[i] = w * J;
1576 this->quadrature_points.resize(base.size());
1577 this->weights.resize(base.size());
1578 for (
unsigned int i = 0; i < base.size(); ++i)
1580 const auto &
q = base.point(i);
1581 const auto w = base.weight(i);
1583 const auto xhat =
q[0];
1584 const auto yhat =
q[1];
1592 this->weights[i] = w * J;
1609 "The split point should be inside the unit reference cell."));
1611 std::array<Point<dim>, dim + 1> vertices;
1612 vertices[0] = split_point;
1619 for (
unsigned int start = 0; start < (dim > 2 ? 2 : 1); ++start)
1621 for (
unsigned int i = 0; i < dim; ++i)
1624 const auto quad = base.compute_affine_transformation(vertices);
1627 this->quadrature_points.insert(this->quadrature_points.end(),
1628 quad.get_points().begin(),
1629 quad.get_points().end());
1630 this->weights.insert(this->weights.end(),
1631 quad.get_weights().begin(),
1632 quad.get_weights().end());
1644 if (dim == 0 || dim == 1)
1648 this->quadrature_points = quad.get_points();
1649 this->weights = quad.get_weights();
1655 const double p = 1.0 / 3.0;
1656 this->quadrature_points.emplace_back(p, p);
1657 this->weights.emplace_back(0.5);
1665 const double Q12 = 1.0 / 2.0;
1666 this->quadrature_points.emplace_back(0.17855872826361643,
1667 0.1550510257216822);
1668 this->quadrature_points.emplace_back(0.07503111022260812,
1669 0.6449489742783178);
1670 this->quadrature_points.emplace_back(0.6663902460147014,
1671 0.1550510257216822);
1672 this->quadrature_points.emplace_back(0.28001991549907407,
1673 0.6449489742783178);
1675 this->weights.emplace_back(0.31804138174397717 *
Q12);
1676 this->weights.emplace_back(0.18195861825602283 *
Q12);
1677 this->weights.emplace_back(0.31804138174397717 *
Q12);
1678 this->weights.emplace_back(0.18195861825602283 *
Q12);
1683 const double p0 = 2.0 / 7.0 -
std::sqrt(15.0) / 21.0;
1684 const double p1 = 2.0 / 7.0 +
std::sqrt(15.0) / 21.0;
1685 const double p2 = 3.0 / 7.0 - 2.0 *
std::sqrt(15.0) / 21.0;
1686 const double p3 = 3.0 / 7.0 + 2.0 *
std::sqrt(15.0) / 21.0;
1687 this->quadrature_points.emplace_back(1.0 / 3.0, 1.0 / 3.0);
1688 this->quadrature_points.emplace_back(
p3,
p0);
1689 this->quadrature_points.emplace_back(
p0,
p3);
1690 this->quadrature_points.emplace_back(
p0,
p0);
1691 this->quadrature_points.emplace_back(
p2,
p1);
1692 this->quadrature_points.emplace_back(
p1,
p2);
1693 this->quadrature_points.emplace_back(
p1,
p1);
1695 const double q12 = 0.5;
1696 const double w0 = 9.0 / 40.0;
1697 const double w1 = 31.0 / 240.0 -
std::sqrt(15.0) / 1200.0;
1698 const double w2 = 31.0 / 240.0 +
std::sqrt(15.0) / 1200.0;
1699 this->weights.emplace_back(
q12 *
w0);
1700 this->weights.emplace_back(
q12 *
w1);
1701 this->weights.emplace_back(
q12 *
w1);
1702 this->weights.emplace_back(
q12 *
w1);
1703 this->weights.emplace_back(
q12 *
w2);
1704 this->weights.emplace_back(
q12 *
w2);
1705 this->weights.emplace_back(
q12 *
w2);
1717 const double Q14 = 1.0 / 4.0;
1718 const double Q16 = 1.0 / 6.0;
1720 this->quadrature_points.emplace_back(
Q14,
Q14,
Q14);
1721 this->weights.emplace_back(
Q16);
1730 const double Q16 = 1.0 / 6.0;
1731 this->weights.emplace_back(0.1223220027573451 *
Q16);
1732 this->weights.emplace_back(0.1280664127107469 *
Q16);
1733 this->weights.emplace_back(0.1325680271444452 *
Q16);
1734 this->weights.emplace_back(0.1406244096604032 *
Q16);
1735 this->weights.emplace_back(0.2244151669175574 *
Q16);
1736 this->weights.emplace_back(0.2520039808095023 *
Q16);
1738 this->quadrature_points.emplace_back(0.1620014916985245,
1740 0.01271836631368145);
1741 this->quadrature_points.emplace_back(0.01090521221118924,
1743 0.3621268299455338);
1744 this->quadrature_points.emplace_back(0.1901170024392839,
1745 0.01140332944455717,
1746 0.3586207204668839);
1747 this->quadrature_points.emplace_back(0.170816925164989,
1749 0.6384932999617267);
1750 this->quadrature_points.emplace_back(0.1586851632274406,
1752 0.1308471689520965);
1753 this->quadrature_points.emplace_back(0.5712260521491151,
1755 0.1403728057942107);
1772 Assert(this->quadrature_points.size() > 0,
1774 "QGaussSimplex is currently only implemented for "
1775 "n_points_1D = 1, 2, 3, and 4 while you are asking for "
1782 template <std::
size_t b_dim>
1783 std::vector<std::array<double, b_dim>>
1786 std::vector<std::array<double, b_dim>> output;
1792 std::sort(
temp.begin(),
temp.end());
1795 output.push_back(
temp);
1797 while (std::next_permutation(
temp.begin(),
temp.end()));
1820 std::array<double, dim + 1>
centroid;
1834 const double b = 1.0 - dim * a;
1835 std::array<double, dim + 1>
b_point;
1847 const double b = (1.0 - 2.0 * a) / 2.0;
1848 std::array<double, dim + 1>
b_point;
1860 auto process_point_3 = [&](
const double a,
const double b,
const double w) {
1861 const double c = 1.0 - (dim - 1.0) * a - b;
1862 std::array<double, dim + 1>
b_point;
1881 b_weights.push_back(1.0000000000000000e+00);
1887 3.3333333333333331e-01);
1895 b_weights.push_back(1.0000000000000000e+00);
1901 2.5000000000000000e-01);
1921 1.3621784253708741e-01);
1923 1.1378215746291261e-01);
1943 b_weights.push_back(2.2500000000000001e-01);
1945 1.2593918054482714e-01);
1947 1.3239415278850619e-01);
1953 5.0844906370206819e-02);
1955 1.1678627572637937e-01);
1957 3.1035245103378439e-01,
1958 8.2851075618373571e-02);
1966 1.1268792571801590e-01);
1968 7.3493043116361956e-02);
1970 4.2546020777081472e-02);
1976 1.0077211055320640e-02);
1978 5.5357181543654717e-02);
1980 3.9922750258167487e-02);
1982 6.0300566479164919e-01,
1983 4.8214285714285710e-02);
1998 1.6545050110792131e-02);
2000 7.7086646185986069e-02);
2002 1.2794417123015558e-01);
2004 1.9868331479735168e-01,
2005 5.5878732903199779e-02);
2011 b_weights.push_back(1.4431560767778717e-01);
2013 3.2458497623198079e-02);
2015 9.5091634267284619e-02);
2017 1.0321737053471824e-01);
2019 2.6311282963463811e-01,
2020 2.7230314174434993e-02);
2028 b_weights.push_back(9.5485289464130846e-02);
2030 4.2329581209967028e-02);
2032 3.1896927832857580e-02);
2034 5.7517163758699996e-01,
2035 3.7207130728334620e-02);
2037 8.1083024109854862e-01,
2038 8.1107708299033420e-03);
2044 2.6426650908408830e-02);
2046 5.2031747563738531e-02);
2048 7.5252561535401989e-03);
2050 4.1763782856934897e-02);
2052 3.6280930261308818e-02);
2054 7.1746406342630831e-01,
2055 7.1569028908444327e-03);
2057 5.8379737830214440e-01,
2058 1.5453486150960340e-02);
2073 b_weights.push_back(9.7135796282798836e-02);
2075 2.5577675658698031e-02);
2077 3.1334700227139071e-02);
2079 7.7827541004774278e-02);
2081 7.9647738927210249e-02);
2083 2.2196298916076568e-01,
2084 4.3283539377289376e-02);
2090 b_weights.push_back(8.1743329146285973e-02);
2092 1.3352968813149567e-02);
2094 4.5957963604744731e-02);
2096 1.6370173373718250e-01,
2097 2.5297757707288385e-02);
2099 3.6914678182781102e-01,
2100 3.4184648162959429e-02);
2102 3.2181299528883545e-01,
2103 6.3904906396424044e-02);
2111 b_weights.push_back(5.8010548912480253e-02);
2113 6.4319281759256394e-05);
2115 2.3173338462425461e-02);
2117 2.9562912335429289e-02);
2119 8.0639799796161822e-03);
2121 3.8134080103702457e-02);
2123 2.5545792330413102e-03,
2124 8.3844221982985519e-03);
2126 7.1835032644207453e-01,
2127 1.0234559352745330e-02);
2129 3.4415910578175279e-02,
2130 2.0524915967988139e-02);
2136 b_weights.push_back(4.7399773556020743e-02);
2138 2.6937059992268701e-02);
2140 9.8691597167933822e-03);
2142 1.6548602561961109e-01,
2143 1.1393881220195230e-02);
2145 9.4298876734520487e-01,
2146 3.6194434433925362e-04);
2148 4.7719037990428043e-01,
2149 2.5739731980456069e-02);
2151 5.9425626948000698e-01,
2152 1.0135871679755789e-02);
2154 8.0117728465834437e-01,
2155 6.5761472770359038e-03);
2157 6.2807184547536599e-01,
2158 1.2907035798861989e-02);
2172 b_weights.push_back(8.5761179732224219e-02);
2179 1.4932478865208237e-01,
2180 1.0290289572953278e-02);
2182 2.8958112563770588e-01,
2183 4.0332476640500554e-02);
2194 1.2727971723358933e-01,
2195 1.5083677576511438e-02);
2197 2.9165567973834094e-01,
2198 2.1783585038607559e-02);
2200 2.5545422863851736e-01,
2201 4.3227363659414209e-02);
2211 b_weights.push_back(6.7960036586831640e-02);
2217 2.7251581777342970e-01,
2218 9.5906810035432631e-03);
2220 1.1092204280346341e-01,
2221 1.4965401105165668e-02);
2223 1.6359740106785048e-01,
2224 2.4179039811593819e-02);
2226 3.0844176089211778e-01,
2227 3.4641276140848373e-02);
2239 1.1897449769695684e-01,
2240 5.0102288385006719e-03);
2242 2.9837288213625779e-01,
2243 1.4436308113533840e-02);
2245 1.7226668782135557e-01,
2246 2.4665753212563674e-02);
2248 3.3686145979634496e-01,
2249 3.8571510787060684e-02);
2260 for (
const std::array<double, dim + 1> &
b_point :
2263 const double volume = (dim == 2 ? 1.0 / 2.0 : 1.0 / 6.0);
2266 for (
int d = 0; d < dim; ++d)
2268 this->quadrature_points.emplace_back(
c_point);
2311 static_cast<unsigned int>(std::round(std::log2(
n_copies)));
2313 ExcMessage(
"The number of copies must be a power of 2."));
2315 const auto reference_cell = ReferenceCells::get_simplex<dim>();
2326 std::vector<Point<dim>> points;
2327 std::vector<double> weights;
2328 for (
const auto &cell : tria.active_cell_iterators())
2330 fe_values.reinit(cell);
2333 points.push_back(fe_values.quadrature_point(
qp));
2334 weights.push_back(fe_values.JxW(
qp));
2359 for (
unsigned int i = 0; i <
quad_line.size(); ++i)
2360 for (
unsigned int j = 0;
j <
quad_tri.size(); ++
j)
2362 this->quadrature_points.emplace_back(
quad_tri.point(
j)[0],
2369 Assert(this->quadrature_points.size() > 0,
2383 const double Q14 = 1.0 / 4.0;
2384 const double Q43 = 4.0 / 3.0;
2386 this->quadrature_points.emplace_back(0, 0,
Q14);
2387 this->weights.emplace_back(
Q43);
2392 this->quadrature_points.emplace_back(-0.26318405556971, -0.26318405556971, 0.54415184401122);
2393 this->quadrature_points.emplace_back(-0.50661630334979, -0.50661630334979, 0.12251482265544);
2394 this->quadrature_points.emplace_back(-0.26318405556971, +0.26318405556971, 0.54415184401122);
2395 this->quadrature_points.emplace_back(-0.50661630334979, +0.50661630334979, 0.12251482265544);
2396 this->quadrature_points.emplace_back(+0.26318405556971, -0.26318405556971, 0.54415184401122);
2397 this->quadrature_points.emplace_back(+0.50661630334979, -0.50661630334979, 0.12251482265544);
2398 this->quadrature_points.emplace_back(+0.26318405556971, +0.26318405556971, 0.54415184401122);
2399 this->quadrature_points.emplace_back(+0.50661630334979, +0.50661630334979, 0.12251482265544);
2402 this->weights.emplace_back(0.10078588207983);
2403 this->weights.emplace_back(0.23254745125351);
2404 this->weights.emplace_back(0.10078588207983);
2405 this->weights.emplace_back(0.23254745125351);
2406 this->weights.emplace_back(0.10078588207983);
2407 this->weights.emplace_back(0.23254745125351);
2408 this->weights.emplace_back(0.10078588207983);
2409 this->weights.emplace_back(0.23254745125351);
2413 Assert(this->quadrature_points.size() > 0,
2491 const std::array<
Point<1>, 1 + 1> &vertices)
const;
2495 const std::array<
Point<2>, 1 + 1> &vertices)
const;
2499 const std::array<
Point<2>, 2 + 1> &vertices)
const;
2503 const std::array<
Point<3>, 1 + 1> &vertices)
const;
2507 const std::array<
Point<3>, 2 + 1> &vertices)
const;
2511 const std::array<
Point<3>, 3 + 1> &vertices)
const;
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)
QGaussLogR(const unsigned int n, const Point< dim > &x0=Point< dim >(), const double alpha=1, const bool factor_out_singular_weight=false)
QGaussLog(const unsigned int n, const bool revert=false)
static std::vector< double > get_quadrature_points(const unsigned int n)
static std::vector< double > get_quadrature_weights(const unsigned int n)
QGaussOneOverR(const unsigned int n, const Point< dim > &singularity, const bool factor_out_singular_weight=false)
static unsigned int quad_size(const Point< dim > &singularity, const unsigned int n)
QGaussPyramid(const unsigned int n_points_1D)
QGaussRadauChebyshev(const unsigned int n, const EndPoint end_point=QGaussRadauChebyshev::EndPoint::left)
Generate a formula with n quadrature points.
QGaussRadau(const unsigned int n, const EndPoint end_point=QGaussRadau::EndPoint::left)
QGaussSimplex(const unsigned int n_points_1D)
QGaussWedge(const unsigned int n_points_1D)
QGauss(const unsigned int n)
QIteratedSimplex(const Quadrature< dim > &base_quadrature, const unsigned int n_copies)
Quadrature< spacedim > compute_affine_transformation(const std::array< Point< spacedim >, dim+1 > &vertices) const
QSimplex(const Quadrature< dim > &quad)
Quadrature< spacedim > mapped_quadrature(const std::vector< std::array< Point< spacedim >, dim+1 > > &simplices) const
QSorted(const Quadrature< dim > &quad)
bool compare_weights(const unsigned int a, const unsigned int b) const
QSplit(const QSimplex< dim > &base, const Point< dim > &split_point)
QTelles(const Quadrature< 1 > &base_quad, const Point< dim > &singularity)
QTrianglePolar(const Quadrature< 1 > &radial_quadrature, const Quadrature< 1 > &angular_quadrature)
QWitherdenVincentSimplex(const unsigned int n_points_1D, const bool use_odd_order=true)
Quadrature & operator=(const Quadrature< dim > &)
constexpr numbers::NumberTraits< Number >::real_type norm_square() const
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_CLOSE
static ::ExceptionBase & ExcNotImplemented()
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
#define AssertIndexRange(index, range)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
@ update_JxW_values
Transformed quadrature weights.
@ update_quadrature_points
Transformed quadrature points.
#define DEAL_II_ASSERT_UNREACHABLE()
#define DEAL_II_NOT_IMPLEMENTED()
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)
constexpr T pow(const T base, const int iexp)
std::vector< double > get_quadrature_weights(const unsigned int n)
std::vector< double > get_quadrature_points(const unsigned int n)
std::vector< double > get_quadrature_points(const unsigned int n)
std::vector< double > get_quadrature_weights(const unsigned int n)
std::vector< long double > compute_quadrature_weights(const std::vector< long double > &x, const int alpha, const int beta)
long double gamma(const unsigned int n)
std::vector< double > get_quadrature_weights(const unsigned int n, const ::QGaussRadauChebyshev< 1 >::EndPoint end_point)
std::vector< double > get_quadrature_points(const unsigned int n, const ::QGaussRadauChebyshev< 1 >::EndPoint end_point)
std::vector< double > get_left_quadrature_weights(const unsigned int n)
std::vector< double > get_left_quadrature_points(const unsigned int n)
std::vector< double > get_quadrature_weights(const unsigned int n, const ::QGaussRadau< 1 >::EndPoint end_point)
std::vector< double > get_quadrature_points(const unsigned int n, const ::QGaussRadau< 1 >::EndPoint end_point)
::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)