16 #include <deal.II/base/quadrature_lib.h> 17 #include <deal.II/base/geometry_info.h> 19 #include <deal.II/base/std_cxx11/bind.h> 26 DEAL_II_NAMESPACE_OPEN
71 const unsigned int m = (n+1)/2;
108 long_double_eps =
static_cast<long double>(std::numeric_limits<long double>::epsilon()),
109 double_eps = static_cast<long double>(std::numeric_limits<double>::epsilon());
121 volatile long double runtime_one = 1.0;
122 const long double tolerance
123 = (runtime_one + long_double_eps != runtime_one
125 std::max (double_eps / 100, long_double_eps * 5)
131 for (
unsigned int i=1; i<=m; ++i)
133 long double z = std::cos(
numbers::PI * (i-.25)/(n+.5));
144 for (
unsigned int j=0; j<n; ++j)
146 const long double p3 = p2;
148 p1 = ((2.*j+1.)*z*p2-j*p3)/(j+1);
150 pp = n*(z*p1-p2)/(z*z-1);
153 while (std::abs(p1/pp) > tolerance);
159 double w = 1./((1.-z*z)*pp*pp);
178 for (
unsigned int i=0; i<points.size(); ++i)
191 const int beta)
const 193 const unsigned int m = q-2;
194 std::vector<long double> x(m);
202 long_double_eps =
static_cast<long double>(std::numeric_limits<long double>::epsilon()),
203 double_eps = static_cast<long double>(std::numeric_limits<double>::epsilon());
208 volatile long double runtime_one = 1.0;
209 const long double tolerance
210 = (runtime_one + long_double_eps != runtime_one
212 std::max (double_eps / 100, long_double_eps * 5)
228 for (
unsigned int i=0; i<m; ++i)
229 x[i] = - std::cos( (
long double) (2*i+1)/(2*m) *
numbers::PI );
231 long double s, J_x, f, delta;
233 for (
unsigned int k=0; k<m; ++k)
235 long double r = x[k];
242 for (
unsigned int i=0; i<k; ++i)
245 J_x = 0.5*(alpha + beta + m + 1)*
JacobiP(r, alpha+1, beta+1, m-1);
246 f =
JacobiP(r, alpha, beta, m);
247 delta = f/(f*s- J_x);
250 while (std::fabs(delta) >= tolerance);
256 x.insert(x.begin(), -1.L);
268 const int beta)
const 270 const unsigned int q = x.size();
271 std::vector<long double> w(q);
274 const long double factor = std::pow(2., alpha+beta+1) *
278 for (
unsigned int i=0; i<q; ++i)
280 s =
JacobiP(x[i], alpha, beta, q-1);
284 w[q-1] *= (alpha + 1);
295 const unsigned int n)
const 299 std::vector<long double> p(n+1);
303 if (n==0)
return p[0];
304 p[1] = ((alpha+beta+2)*x + (alpha-beta))/2;
305 if (n==1)
return p[1];
307 for (
unsigned int i=1; i<=(n-1); ++i)
309 const int v = 2*i + alpha + beta;
310 const int a1 = 2*(i+1)*(i + alpha + beta + 1)*v;
311 const int a2 = (v + 1)*(alpha*alpha - beta*beta);
312 const int a3 = v*(v + 1)*(v + 2);
313 const int a4 = 2*(i+alpha)*(i+beta)*(v + 2);
315 p[i+1] =
static_cast<long double>( (a2 + a3*x)*p[i] - a4*p[i-1])/a1;
325 long double result = n - 1;
326 for (
int i=n-2; i>1; --i)
349 static const double xpts[] = { 0.0, 1.0 };
350 static const double wts[] = { 0.5, 0.5 };
352 for (
unsigned int i=0; i<this->
size(); ++i)
366 static const double xpts[] = { 0.0, 0.5, 1.0 };
367 static const double wts[] = { 1./6., 2./3., 1./6. };
369 for (
unsigned int i=0; i<this->
size(); ++i)
383 static const double xpts[] = { 0.0, .25, .5, .75, 1.0 };
384 static const double wts[] = { 7./90., 32./90., 12./90., 32./90., 7./90. };
386 for (
unsigned int i=0; i<this->
size(); ++i)
400 static const double xpts[] = { 0.0, 1./6., 1./3., .5, 2./3., 5./6., 1.0 };
401 static const double wts[] = { 41./840., 216./840., 27./840., 272./840.,
402 27./840., 216./840., 41./840.
405 for (
unsigned int i=0; i<this->
size(); ++i)
419 std::vector<double> p = get_quadrature_points(n);
420 std::vector<double> w = get_quadrature_weights(n);
422 for (
unsigned int i=0; i<this->
size(); ++i)
428 this->
weights[i] = revert ? w[n-1-i] : w[i];
436 std::vector<double> q_points(n);
441 q_points[0] = 0.3333333333333333;
445 q_points[0] = 0.1120088061669761;
446 q_points[1] = 0.6022769081187381;
450 q_points[0] = 0.06389079308732544;
451 q_points[1] = 0.3689970637156184;
452 q_points[2] = 0.766880303938942;
456 q_points[0] = 0.04144848019938324;
457 q_points[1] = 0.2452749143206022;
458 q_points[2] = 0.5561654535602751;
459 q_points[3] = 0.848982394532986;
463 q_points[0] = 0.02913447215197205;
464 q_points[1] = 0.1739772133208974;
465 q_points[2] = 0.4117025202849029;
466 q_points[3] = 0.6773141745828183;
467 q_points[4] = 0.89477136103101;
471 q_points[0] = 0.02163400584411693;
472 q_points[1] = 0.1295833911549506;
473 q_points[2] = 0.3140204499147661;
474 q_points[3] = 0.5386572173517997;
475 q_points[4] = 0.7569153373774084;
476 q_points[5] = 0.922668851372116;
481 q_points[0] = 0.0167193554082585;
482 q_points[1] = 0.100185677915675;
483 q_points[2] = 0.2462942462079286;
484 q_points[3] = 0.4334634932570557;
485 q_points[4] = 0.6323509880476823;
486 q_points[5] = 0.81111862674023;
487 q_points[6] = 0.940848166743287;
491 q_points[0] = 0.01332024416089244;
492 q_points[1] = 0.07975042901389491;
493 q_points[2] = 0.1978710293261864;
494 q_points[3] = 0.354153994351925;
495 q_points[4] = 0.5294585752348643;
496 q_points[5] = 0.7018145299391673;
497 q_points[6] = 0.849379320441094;
498 q_points[7] = 0.953326450056343;
502 q_points[0] = 0.01086933608417545;
503 q_points[1] = 0.06498366633800794;
504 q_points[2] = 0.1622293980238825;
505 q_points[3] = 0.2937499039716641;
506 q_points[4] = 0.4466318819056009;
507 q_points[5] = 0.6054816627755208;
508 q_points[6] = 0.7541101371585467;
509 q_points[7] = 0.877265828834263;
510 q_points[8] = 0.96225055941096;
514 q_points[0] = 0.00904263096219963;
515 q_points[1] = 0.05397126622250072;
516 q_points[2] = 0.1353118246392511;
517 q_points[3] = 0.2470524162871565;
518 q_points[4] = 0.3802125396092744;
519 q_points[5] = 0.5237923179723384;
520 q_points[6] = 0.6657752055148032;
521 q_points[7] = 0.7941904160147613;
522 q_points[8] = 0.898161091216429;
523 q_points[9] = 0.9688479887196;
528 q_points[0] = 0.007643941174637681;
529 q_points[1] = 0.04554182825657903;
530 q_points[2] = 0.1145222974551244;
531 q_points[3] = 0.2103785812270227;
532 q_points[4] = 0.3266955532217897;
533 q_points[5] = 0.4554532469286375;
534 q_points[6] = 0.5876483563573721;
535 q_points[7] = 0.7139638500230458;
536 q_points[8] = 0.825453217777127;
537 q_points[9] = 0.914193921640008;
538 q_points[10] = 0.973860256264123;
542 q_points[0] = 0.006548722279080035;
543 q_points[1] = 0.03894680956045022;
544 q_points[2] = 0.0981502631060046;
545 q_points[3] = 0.1811385815906331;
546 q_points[4] = 0.2832200676673157;
547 q_points[5] = 0.398434435164983;
548 q_points[6] = 0.5199526267791299;
549 q_points[7] = 0.6405109167754819;
550 q_points[8] = 0.7528650118926111;
551 q_points[9] = 0.850240024421055;
552 q_points[10] = 0.926749682988251;
553 q_points[11] = 0.977756129778486;
569 std::vector<double> quadrature_weights(n);
574 quadrature_weights[0] = -1.0;
577 quadrature_weights[0] = -0.7185393190303845;
578 quadrature_weights[1] = -0.2814606809696154;
582 quadrature_weights[0] = -0.5134045522323634;
583 quadrature_weights[1] = -0.3919800412014877;
584 quadrature_weights[2] = -0.0946154065661483;
588 quadrature_weights[0] =-0.3834640681451353;
589 quadrature_weights[1] =-0.3868753177747627;
590 quadrature_weights[2] =-0.1904351269501432;
591 quadrature_weights[3] =-0.03922548712995894;
595 quadrature_weights[0] =-0.2978934717828955;
596 quadrature_weights[1] =-0.3497762265132236;
597 quadrature_weights[2] =-0.234488290044052;
598 quadrature_weights[3] =-0.0989304595166356;
599 quadrature_weights[4] =-0.01891155214319462;
603 quadrature_weights[0] = -0.2387636625785478;
604 quadrature_weights[1] = -0.3082865732739458;
605 quadrature_weights[2] = -0.2453174265632108;
606 quadrature_weights[3] = -0.1420087565664786;
607 quadrature_weights[4] = -0.05545462232488041;
608 quadrature_weights[5] = -0.01016895869293513;
613 quadrature_weights[0] = -0.1961693894252476;
614 quadrature_weights[1] = -0.2703026442472726;
615 quadrature_weights[2] = -0.239681873007687;
616 quadrature_weights[3] = -0.1657757748104267;
617 quadrature_weights[4] = -0.0889432271377365;
618 quadrature_weights[5] = -0.03319430435645653;
619 quadrature_weights[6] = -0.005932787015162054;
623 quadrature_weights[0] = -0.164416604728002;
624 quadrature_weights[1] = -0.2375256100233057;
625 quadrature_weights[2] = -0.2268419844319134;
626 quadrature_weights[3] = -0.1757540790060772;
627 quadrature_weights[4] = -0.1129240302467932;
628 quadrature_weights[5] = -0.05787221071771947;
629 quadrature_weights[6] = -0.02097907374214317;
630 quadrature_weights[7] = -0.003686407104036044;
634 quadrature_weights[0] = -0.1400684387481339;
635 quadrature_weights[1] = -0.2097722052010308;
636 quadrature_weights[2] = -0.211427149896601;
637 quadrature_weights[3] = -0.1771562339380667;
638 quadrature_weights[4] = -0.1277992280331758;
639 quadrature_weights[5] = -0.07847890261203835;
640 quadrature_weights[6] = -0.0390225049841783;
641 quadrature_weights[7] = -0.01386729555074604;
642 quadrature_weights[8] = -0.002408041036090773;
646 quadrature_weights[0] = -0.12095513195457;
647 quadrature_weights[1] = -0.1863635425640733;
648 quadrature_weights[2] = -0.1956608732777627;
649 quadrature_weights[3] = -0.1735771421828997;
650 quadrature_weights[4] = -0.135695672995467;
651 quadrature_weights[5] = -0.0936467585378491;
652 quadrature_weights[6] = -0.05578772735275126;
653 quadrature_weights[7] = -0.02715981089692378;
654 quadrature_weights[8] = -0.00951518260454442;
655 quadrature_weights[9] = -0.001638157633217673;
660 quadrature_weights[0] = -0.1056522560990997;
661 quadrature_weights[1] = -0.1665716806006314;
662 quadrature_weights[2] = -0.1805632182877528;
663 quadrature_weights[3] = -0.1672787367737502;
664 quadrature_weights[4] = -0.1386970574017174;
665 quadrature_weights[5] = -0.1038334333650771;
666 quadrature_weights[6] = -0.06953669788988512;
667 quadrature_weights[7] = -0.04054160079499477;
668 quadrature_weights[8] = -0.01943540249522013;
669 quadrature_weights[9] = -0.006737429326043388;
670 quadrature_weights[10] = -0.001152486965101561;
674 quadrature_weights[0] = -0.09319269144393;
675 quadrature_weights[1] = -0.1497518275763289;
676 quadrature_weights[2] = -0.166557454364573;
677 quadrature_weights[3] = -0.1596335594369941;
678 quadrature_weights[4] = -0.1384248318647479;
679 quadrature_weights[5] = -0.1100165706360573;
680 quadrature_weights[6] = -0.07996182177673273;
681 quadrature_weights[7] = -0.0524069547809709;
682 quadrature_weights[8] = -0.03007108900074863;
683 quadrature_weights[9] = -0.01424924540252916;
684 quadrature_weights[10] = -0.004899924710875609;
685 quadrature_weights[11] = -0.000834029009809656;
693 return quadrature_weights;
701 const bool factor_out_singularity) :
703 (alpha == 1 ? n : 2*n ) : 4*n ),
704 fraction( ( (origin[0] == 0) || (origin[0] == 1.) ) ? 1. : origin[0] )
724 Assert( (fraction >= 0) && (fraction <= 1),
729 unsigned int ns_offset = (fraction == 1) ? n : 2*n;
731 for (
unsigned int i=0, j=ns_offset; i<n; ++i, ++j)
739 if ( (alpha != 1) || (fraction != 1) )
742 this->
weights[j] = -std::log(alpha/fraction)*quad.
weight(i)*fraction;
752 this->
weights[j+n] = -std::log(alpha/(1-fraction))*quad.
weight(i)*(1-fraction);
755 if (factor_out_singularity ==
true)
756 for (
unsigned int i=0; i<
size(); ++i)
759 ExcMessage(
"The singularity cannot be on a Gauss point of the same order!") );
760 double denominator = std::log(std::abs( (this->
quadrature_points[i]-origin)[0] )/alpha);
761 Assert( denominator != 0.0,
762 ExcMessage(
"The quadrature formula you are using does not allow to " 763 "factor out the singularity, which is zero at one point."));
764 this->
weights[i] /= denominator;
771 const unsigned int n)
775 bool on_vertex=
false;
776 for (
unsigned int i=0; i<2; ++i)
777 if ( ( std::abs(singularity[i] ) < eps ) ||
778 ( std::abs(singularity[i]-1) < eps ) )
780 if (on_edge && (std::abs( (singularity-
Point<2>(.5, .5)).norm_square()-.5)
783 if (on_vertex)
return (2*n*n);
784 if (on_edge)
return (4*n*n);
791 const bool factor_out_singularity) :
800 std::vector<QGaussOneOverR<2> > quads;
801 std::vector<Point<2> > origins;
810 origins.push_back(
Point<2>(singularity[0],0.));
811 origins.push_back(
Point<2>(0.,singularity[1]));
812 origins.push_back(singularity);
817 unsigned int q_id = 0;
820 for (
unsigned int box=0; box<4; ++box)
823 dist =
Point<2>(std::abs(dist[0]), std::abs(dist[1]));
824 double area = dist[0]*dist[1];
826 for (
unsigned int q=0; q<quads[box].size(); ++q, ++q_id)
828 const Point<2> &qp = quads[box].point(q);
831 Point<2>(dist[0]*qp[0], dist[1]*qp[1]);
832 this->
weights[q_id] = quads[box].weight(q)*area;
840 const unsigned int vertex_index,
841 const bool factor_out_singularity) :
877 std::vector<double> &ws = this->
weights;
880 for (
unsigned int q=0; q<gauss.
size(); ++q)
884 ps[q][1] = gp[0] * std::tan(pi4*gp[1]);
885 ws[q] = gauss.
weight(q)*pi4/std::cos(pi4 *gp[1]);
886 if (factor_out_singularity)
890 ws[gauss.
size()+q] = ws[q];
891 ps[gauss.
size()+q][0] = ps[q][1];
892 ps[gauss.
size()+q][1] = ps[q][0];
897 switch (vertex_index)
914 double R00 = std::cos(theta), R01 = -std::sin(theta);
915 double R10 = std::sin(theta), R11 = std::cos(theta);
917 if (vertex_index != 0)
918 for (
unsigned int q=0; q<
size(); ++q)
920 double x = ps[q][0]-.5, y = ps[q][1]-.5;
922 ps[q][0] = R00*x + R01*y + .5;
923 ps[q][1] = R10*x + R11*y + .5;
932 std::vector<unsigned int> permutation(quad.
size());
933 for (
unsigned int i=0; i<quad.
size(); ++i)
936 std::sort(permutation.begin(),
939 std_cxx11::ref(*
this),
943 for (
unsigned int i=0; i<quad.
size(); ++i)
953 const unsigned int b)
const 1038 const unsigned int n,
const Point<dim> &singularity)
1062 const double eta_bar = singularity[0] * 2. - 1.;
1063 const double eta_star = eta_bar * eta_bar - 1.;
1067 std::vector<double> weights_dummy(
weights.size());
1068 unsigned int cont = 0;
1069 const double tol = 1e-10;
1075 weights_dummy[d-cont] =
weights[d];
1088 weights.resize(weights_dummy.size()-1);
1092 weights[d] = weights_dummy[d];
1096 if (std::abs(eta_star) <= tol)
1098 gamma_bar = std::pow((eta_bar * eta_star + std::abs(eta_star)),1.0 / 3.0)
1099 + std::pow((eta_bar * eta_star - std::abs(eta_star)), 1.0 / 3.0)
1104 gamma_bar = (eta_bar * eta_star + std::abs(eta_star))/std::abs(eta_bar * eta_star + std::abs(eta_star))*
1105 std::pow(std::abs(eta_bar * eta_star + std::abs(eta_star)),1.0 / 3.0)
1106 + (eta_bar * eta_star - std::abs(eta_star))/std::abs(eta_bar * eta_star - std::abs(eta_star))*
1107 std::pow(std::abs(eta_bar * eta_star - std::abs(eta_star)), 1.0 / 3.0)
1113 double eta = (std::pow(gamma - gamma_bar, 3.0)
1114 + gamma_bar * (gamma_bar * gamma_bar + 3))
1115 / (1 + 3 * gamma_bar * gamma_bar);
1117 double J = 3 * ((gamma - gamma_bar) *(gamma - gamma_bar))
1118 / (1 + 3 * gamma_bar * gamma_bar);
1133 std::vector<double> points(n);
1135 for (
unsigned short i=0; i<n; ++i)
1139 points[i] = 1./2.*(1.+std::cos(
numbers::PI*(1.+
double(2*i+1)/
double(2*(n-1)+2))));
1150 std::vector<double>
weights(n);
1152 for (
unsigned short i=0; i<n; ++i)
1170 std::vector<double> p=get_quadrature_points(n);
1171 std::vector<double> w=get_quadrature_weights(n);
1173 for (
unsigned int i=0; i<this->
size(); ++i)
1198 std::vector<double> points(n);
1200 for (
unsigned short i=0; i<n; ++i)
1208 points[i] = 1./2.*(1.-std::cos(
numbers::PI*(1+2*
double(i)/(2*
double(n-1)+1.))));
1214 points[i] = 1./2.*(1.-std::cos(
numbers::PI*(2*
double(n-1-i)/(2*
double(n-1)+1.))));
1219 Assert (
false,
ExcMessage (
"This constructor can only be called with either " 1220 "QGaussRadauChebyshev::left or QGaussRadauChebyshev::right as " 1221 "second argument."));
1234 std::vector<double>
weights(n);
1236 for (
unsigned short i=0; i<n; ++i)
1240 if (ep==left && i==0)
1242 else if (ep==right && i==(n-1))
1263 for (
unsigned int i=0; i<this->
size(); ++i)
1296 std::vector<double> points(n);
1298 for (
unsigned short i=0; i<n; ++i)
1302 points[i] = 1./2.*(1.+std::cos(
numbers::PI*(1+
double(i)/
double(n-1))));
1313 std::vector<double>
weights(n);
1315 for (
unsigned short i=0; i<n; ++i)
1319 if (i==0 || i==(n-1))
1334 Assert(n>1,
ExcMessage(
"Need at least two points for Gauss-Lobatto quadrature rule"));
1338 for (
unsigned int i=0; i<this->
size(); ++i)
1391 DEAL_II_NAMESPACE_CLOSE
static long double gamma(const unsigned int n)
QGaussLog(const unsigned int n, const bool revert=false)
std::vector< double > weights
QGaussOneOverR(const unsigned int n, const Point< dim > singularity, const bool factor_out_singular_weight=false)
static std::vector< double > get_quadrature_weights(const unsigned int n)
Computes the weights of the quadrature formula.
static std::vector< double > get_quadrature_weights(const unsigned int n, EndPoint ep)
Computes the weights of the quadrature formula.
QGaussChebyshev(const unsigned int n)
Generate a formula with n quadrature points.
static std::vector< double > get_quadrature_points(const unsigned int n)
static Point< dim > unit_cell_vertex(const unsigned int vertex)
QGauss(const unsigned int n)
static std::vector< double > get_quadrature_weights(const unsigned int n)
Computes the weights of the quadrature formula.
static::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
const Point< dim > & point(const unsigned int i) const
long double JacobiP(const long double x, const int alpha, const int beta, const unsigned int n) const
unsigned int size() const
static::ExceptionBase & ExcMessage(std::string arg1)
std::vector< long double > compute_quadrature_points(const unsigned int q, const int alpha, const int beta) const
#define Assert(cond, exc)
double weight(const unsigned int i) const
static std::vector< double > get_quadrature_weights(const unsigned int n)
QGaussRadauChebyshev(const unsigned int n, EndPoint ep=QGaussRadauChebyshev::left)
Generate a formula with n quadrature points.
static std::vector< double > get_quadrature_points(const unsigned int n)
Computes the points of the quadrature formula.
QGaussLobatto(const unsigned int n)
std::vector< long double > compute_quadrature_weights(const std::vector< long double > &x, const int alpha, const int beta) const
std::vector< Point< dim > > quadrature_points
QTelles(const Quadrature< 1 > &base_quad, const Point< dim > &singularity)
static std::vector< double > get_quadrature_points(const unsigned int n, EndPoint ep)
Computes the points of the quadrature formula.
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
bool compare_weights(const unsigned int a, const unsigned int b) const
static std::vector< double > get_quadrature_points(const unsigned int n)
Computes the points of the quadrature formula.
QGaussLobattoChebyshev(const unsigned int n)
Generate a formula with n quadrature points.
static::ExceptionBase & ExcNotImplemented()
QSorted(const Quadrature< dim > &quad)
QGaussLogR(const unsigned int n, const Point< dim > x0=Point< dim >(), const double alpha=1, const bool factor_out_singular_weight=false)
static unsigned int quad_size(const Point< dim > singularity, const unsigned int n)
static::ExceptionBase & ExcInternalError()