Reference documentation for deal.II version Git ede8f93e86 2020-12-03 14:59:20 -0700
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2020 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
19
20 #include <algorithm>
21 #include <cmath>
22 #include <functional>
23 #include <limits>
24
25
27
28
29 // please note: for a given dimension, we need the quadrature formulae
30 // for all lower dimensions as well. That is why in this file the check
31 // is for deal_II_dimension >= any_number and not for ==
32
33
34
35 template <>
36 QGauss<0>::QGauss(const unsigned int)
37  : // there are n_q^dim == 1
38  // points
40 {
41  // the single quadrature point gets unit
42  // weight
43  this->weights[0] = 1;
44 }
45
46
47
48 template <>
50  : // there are n_q^dim == 1
51  // points
53 {
54  // the single quadrature point gets unit
55  // weight
56  this->weights[0] = 1;
57 }
58
59
60
61 template <>
62 QGauss<1>::QGauss(const unsigned int n)
64 {
65  if (n == 0)
66  return;
67
68  std::vector<long double> points =
69  Polynomials::jacobi_polynomial_roots<long double>(n, 0, 0);
70
71  for (unsigned int i = 0; i < (points.size() + 1) / 2; ++i)
72  {
74  this->quadrature_points[n - i - 1][0] = 1. - points[i];
75
76  // derivative of Jacobi polynomial
77  const long double pp =
78  0.5 * (n + 1) *
79  Polynomials::jacobi_polynomial_value(n - 1, 1, 1, points[i]);
80  const long double x = -1. + 2. * points[i];
81  const double w = 1. / ((1. - x * x) * pp * pp);
82  this->weights[i] = w;
83  this->weights[n - i - 1] = w;
84  }
85 }
86
87 namespace internal
88 {
89  namespace QGaussLobatto
90  {
95  long double
96  gamma(const unsigned int n)
97  {
98  long double result = n - 1;
99  for (int i = n - 2; i > 1; --i)
100  result *= i;
101  return result;
102  }
103
104
105
113  std::vector<long double>
115  const int alpha,
116  const int beta)
117  {
118  const unsigned int q = x.size();
119  std::vector<long double> w(q);
120
121  const long double factor =
122  std::pow(2., alpha + beta + 1) * gamma(alpha + q) * gamma(beta + q) /
123  ((q - 1) * gamma(q) * gamma(alpha + beta + q + 1));
124  for (unsigned int i = 0; i < q; ++i)
125  {
126  const long double s =
127  Polynomials::jacobi_polynomial_value(q - 1, alpha, beta, x[i]);
128  w[i] = factor / (s * s);
129  }
130  w[0] *= (beta + 1);
131  w[q - 1] *= (alpha + 1);
132
133  return w;
134  }
135  } // namespace QGaussLobatto
136 } // namespace internal
137
138
139 #ifndef DOXYGEN
140 template <>
141 QGaussLobatto<1>::QGaussLobatto(const unsigned int n)
143 {
144  Assert(n >= 2, ExcNotImplemented());
145
146  std::vector<long double> points =
147  Polynomials::jacobi_polynomial_roots<long double>(n - 2, 1, 1);
148  points.insert(points.begin(), 0);
149  points.push_back(1.);
150  std::vector<long double> w =
152
153  // scale weights to the interval [0.0, 1.0]:
154  for (unsigned int i = 0; i < points.size(); ++i)
155  {
157  this->weights[i] = 0.5 * w[i];
158  }
159 }
160 #endif
161
162
163 template <>
166 {
168  this->weights[0] = 1.0;
169 }
170
171
172
173 template <>
176 {
177  static const double xpts[] = {0.0, 1.0};
178  static const double wts[] = {0.5, 0.5};
179
180  for (unsigned int i = 0; i < this->size(); ++i)
181  {
183  this->weights[i] = wts[i];
184  }
185 }
186
187
188
189 template <>
192 {
193  static const double xpts[] = {0.0, 0.5, 1.0};
194  static const double wts[] = {1. / 6., 2. / 3., 1. / 6.};
195
196  for (unsigned int i = 0; i < this->size(); ++i)
197  {
199  this->weights[i] = wts[i];
200  }
201 }
202
203
204
205 template <>
208 {
209  static const double xpts[] = {0.0, .25, .5, .75, 1.0};
210  static const double wts[] = {
211  7. / 90., 32. / 90., 12. / 90., 32. / 90., 7. / 90.};
212
213  for (unsigned int i = 0; i < this->size(); ++i)
214  {
216  this->weights[i] = wts[i];
217  }
218 }
219
220
221
222 template <>
225 {
226  static const double xpts[] = {
227  0.0, 1. / 6., 1. / 3., .5, 2. / 3., 5. / 6., 1.0};
228  static const double wts[] = {41. / 840.,
229  216. / 840.,
230  27. / 840.,
231  272. / 840.,
232  27. / 840.,
233  216. / 840.,
234  41. / 840.};
235
236  for (unsigned int i = 0; i < this->size(); ++i)
237  {
239  this->weights[i] = wts[i];
240  }
241 }
242
243
244 template <>
245 QGaussLog<1>::QGaussLog(const unsigned int n, const bool revert)
247 {
250
251  for (unsigned int i = 0; i < this->size(); ++i)
252  {
253  // Using the change of variables x=1-t, it's possible to show
254  // that int f(x)ln|1-x| = int f(1-t) ln|t|, which implies that
255  // we can use this quadrature formula also with weight ln|1-x|.
257  revert ? Point<1>(1 - p[n - 1 - i]) : Point<1>(p[i]);
258  this->weights[i] = revert ? w[n - 1 - i] : w[i];
259  }
260 }
261
262 template <>
263 std::vector<double>
265 {
266  std::vector<double> q_points(n);
267
268  switch (n)
269  {
270  case 1:
271  q_points[0] = 0.3333333333333333;
272  break;
273
274  case 2:
275  q_points[0] = 0.1120088061669761;
276  q_points[1] = 0.6022769081187381;
277  break;
278
279  case 3:
280  q_points[0] = 0.06389079308732544;
281  q_points[1] = 0.3689970637156184;
282  q_points[2] = 0.766880303938942;
283  break;
284
285  case 4:
286  q_points[0] = 0.04144848019938324;
287  q_points[1] = 0.2452749143206022;
288  q_points[2] = 0.5561654535602751;
289  q_points[3] = 0.848982394532986;
290  break;
291
292  case 5:
293  q_points[0] = 0.02913447215197205;
294  q_points[1] = 0.1739772133208974;
295  q_points[2] = 0.4117025202849029;
296  q_points[3] = 0.6773141745828183;
297  q_points[4] = 0.89477136103101;
298  break;
299
300  case 6:
301  q_points[0] = 0.02163400584411693;
302  q_points[1] = 0.1295833911549506;
303  q_points[2] = 0.3140204499147661;
304  q_points[3] = 0.5386572173517997;
305  q_points[4] = 0.7569153373774084;
306  q_points[5] = 0.922668851372116;
307  break;
308
309
310  case 7:
311  q_points[0] = 0.0167193554082585;
312  q_points[1] = 0.100185677915675;
313  q_points[2] = 0.2462942462079286;
314  q_points[3] = 0.4334634932570557;
315  q_points[4] = 0.6323509880476823;
316  q_points[5] = 0.81111862674023;
317  q_points[6] = 0.940848166743287;
318  break;
319
320  case 8:
321  q_points[0] = 0.01332024416089244;
322  q_points[1] = 0.07975042901389491;
323  q_points[2] = 0.1978710293261864;
324  q_points[3] = 0.354153994351925;
325  q_points[4] = 0.5294585752348643;
326  q_points[5] = 0.7018145299391673;
327  q_points[6] = 0.849379320441094;
328  q_points[7] = 0.953326450056343;
329  break;
330
331  case 9:
332  q_points[0] = 0.01086933608417545;
333  q_points[1] = 0.06498366633800794;
334  q_points[2] = 0.1622293980238825;
335  q_points[3] = 0.2937499039716641;
336  q_points[4] = 0.4466318819056009;
337  q_points[5] = 0.6054816627755208;
338  q_points[6] = 0.7541101371585467;
339  q_points[7] = 0.877265828834263;
340  q_points[8] = 0.96225055941096;
341  break;
342
343  case 10:
344  q_points[0] = 0.00904263096219963;
345  q_points[1] = 0.05397126622250072;
346  q_points[2] = 0.1353118246392511;
347  q_points[3] = 0.2470524162871565;
348  q_points[4] = 0.3802125396092744;
349  q_points[5] = 0.5237923179723384;
350  q_points[6] = 0.6657752055148032;
351  q_points[7] = 0.7941904160147613;
352  q_points[8] = 0.898161091216429;
353  q_points[9] = 0.9688479887196;
354  break;
355
356
357  case 11:
358  q_points[0] = 0.007643941174637681;
359  q_points[1] = 0.04554182825657903;
360  q_points[2] = 0.1145222974551244;
361  q_points[3] = 0.2103785812270227;
362  q_points[4] = 0.3266955532217897;
363  q_points[5] = 0.4554532469286375;
364  q_points[6] = 0.5876483563573721;
365  q_points[7] = 0.7139638500230458;
366  q_points[8] = 0.825453217777127;
367  q_points[9] = 0.914193921640008;
368  q_points[10] = 0.973860256264123;
369  break;
370
371  case 12:
372  q_points[0] = 0.006548722279080035;
373  q_points[1] = 0.03894680956045022;
374  q_points[2] = 0.0981502631060046;
375  q_points[3] = 0.1811385815906331;
376  q_points[4] = 0.2832200676673157;
377  q_points[5] = 0.398434435164983;
378  q_points[6] = 0.5199526267791299;
379  q_points[7] = 0.6405109167754819;
380  q_points[8] = 0.7528650118926111;
381  q_points[9] = 0.850240024421055;
382  q_points[10] = 0.926749682988251;
383  q_points[11] = 0.977756129778486;
384  break;
385
386  default:
387  Assert(false, ExcNotImplemented());
388  break;
389  }
390
391  return q_points;
392 }
393
394
395 template <>
396 std::vector<double>
398 {
400
401  switch (n)
402  {
403  case 1:
405  break;
406  case 2:
409  break;
410
411  case 3:
415  break;
416
417  case 4:
422  break;
423
424  case 5:
430  break;
431
432  case 6:
439  break;
440
441
442  case 7:
450  break;
451
452  case 8:
461  break;
462
463  case 9:
473  break;
474
475  case 10:
486  break;
487
488
489  case 11:
501  break;
502
503  case 12:
516  break;
517
518  default:
519  Assert(false, ExcNotImplemented());
520  break;
521  }
522
524 }
525
526
527 template <>
528 QGaussLogR<1>::QGaussLogR(const unsigned int n,
529  const Point<1> origin,
530  const double alpha,
531  const bool factor_out_singularity)
533  ((origin[0] == 0) || (origin[0] == 1)) ? (alpha == 1 ? n : 2 * n) : 4 * n)
534  , fraction(((origin[0] == 0) || (origin[0] == 1.)) ? 1. : origin[0])
535 {
536  // The three quadrature formulas that make this one up. There are
537  // at most two when the origin is one of the extremes, and there is
538  // only one if the origin is one of the extremes and alpha is
539  // equal to one.
540  //
541  // If alpha is different from one, then we need a correction which
542  // is performed with a standard Gauss quadrature rule on each
543  // segment. This is not needed in the standard case where alpha is
544  // equal to one and the origin is on one of the extremes. We
545  // integrate with weight ln|(x-o)/alpha|. In the easy cases, we
546  // only need n quadrature points. In the most difficult one, we
547  // need 2*n points for the first segment, and 2*n points for the
548  // second segment.
549  QGaussLog<1> quad1(n, origin[0] != 0);
552
553  // Check that the origin is inside 0,1
554  Assert((fraction >= 0) && (fraction <= 1),
555  ExcMessage("Origin is outside [0,1]."));
556
557  // Non singular offset. This is the start of non singular quad
558  // points.
559  unsigned int ns_offset = (fraction == 1) ? n : 2 * n;
560
561  for (unsigned int i = 0, j = ns_offset; i < n; ++i, ++j)
562  {
563  // The first i quadrature points are the same as quad1, and
564  // are by default singular.
566  this->weights[i] = quad1.weight(i) * fraction;
567
568  // We need to scale with -log|fraction*alpha|
569  if ((alpha != 1) || (fraction != 1))
570  {
572  this->weights[j] =
573  -std::log(alpha / fraction) * quad.weight(i) * fraction;
574  }
575  // In case we need the second quadrature as well, do it now.
576  if (fraction != 1)
577  {
579  quad2.point(i) * (1 - fraction) + Point<1>(fraction);
580  this->weights[i + n] = quad2.weight(i) * (1 - fraction);
581
582  // We need to scale with -log|fraction*alpha|
584  quad.point(i) * (1 - fraction) + Point<1>(fraction);
585  this->weights[j + n] =
586  -std::log(alpha / (1 - fraction)) * quad.weight(i) * (1 - fraction);
587  }
588  }
589  if (factor_out_singularity == true)
590  for (unsigned int i = 0; i < size(); ++i)
591  {
592  Assert(
594  ExcMessage(
595  "The singularity cannot be on a Gauss point of the same order!"));
596  double denominator =
597  std::log(std::abs((this->quadrature_points[i] - origin)[0]) / alpha);
598  Assert(denominator != 0.0,
599  ExcMessage(
600  "The quadrature formula you are using does not allow to "
601  "factor out the singularity, which is zero at one point."));
602  this->weights[i] /= denominator;
603  }
604 }
605
606
607 template <>
608 unsigned int
609 QGaussOneOverR<2>::quad_size(const Point<2> singularity, const unsigned int n)
610 {
611  const double eps = 1e-8;
612  const bool on_edge =
613  std::any_of(singularity.begin_raw(),
614  singularity.end_raw(),
615  [eps](double coord) {
616  return std::abs(coord) < eps || std::abs(coord - 1.) < eps;
617  });
618  const bool on_vertex =
619  on_edge &&
620  std::abs((singularity - Point<2>(.5, .5)).norm_square() - .5) < eps;
621  if (on_vertex)
622  return 2 * n * n;
623  else if (on_edge)
624  return 4 * n * n;
625  else
626  return 8 * n * n;
627 }
628
629 template <>
631  const Point<2> singularity,
632  const bool factor_out_singularity)
634 {
635  // We treat all the cases in the
636  // same way. Split the element in 4
637  // pieces, measure the area, if
638  // it's relevant, add the
639  // quadrature connected to that
640  // singularity.
642  std::vector<Point<2>> origins;
643  // Id of the corner with a
644  // singularity
649
650  origins.emplace_back(0., 0.);
651  origins.emplace_back(singularity[0], 0.);
652  origins.emplace_back(0., singularity[1]);
653  origins.push_back(singularity);
654
655  // Lexicographical ordering.
656
657  double eps = 1e-8;
658  unsigned int q_id = 0; // Current quad point index.
659  Tensor<1, 2> dist;
660
661  for (unsigned int box = 0; box < 4; ++box)
662  {
663  dist = (singularity - GeometryInfo<2>::unit_cell_vertex(box));
664  dist = Point<2>(std::abs(dist[0]), std::abs(dist[1]));
665  double area = dist[0] * dist[1];
666  if (area > eps)
667  for (unsigned int q = 0; q < quads[box].size(); ++q, ++q_id)
668  {
669  const Point<2> &qp = quads[box].point(q);
671  origins[box] + Point<2>(dist[0] * qp[0], dist[1] * qp[1]);
672  this->weights[q_id] = quads[box].weight(q) * area;
673  }
674  }
675 }
676
677
678 template <>
680  const unsigned int vertex_index,
681  const bool factor_out_singularity)
682  : Quadrature<2>(2 * n * n)
683 {
684  // This version of the constructor
685  // works only for the 4
686  // vertices. If you need a more
687  // general one, you should use the
688  // one with the Point<2> in the
689  // constructor.
690  AssertIndexRange(vertex_index, 4);
691
693  // element.
694  QGauss<2> gauss(n);
695
696  Assert(gauss.size() == n * n, ExcInternalError());
697
698  // For the moment we only implemented this for the vertices of a
699  // quadrilateral. We are planning to do this also for the support
700  // points of arbitrary FE_Q elements, to allow the use of this
701  // class in boundary element programs with higher order mappings.
702  AssertIndexRange(vertex_index, 4);
703
704  // We create only the first one. All other pieces are rotation of
705  // this one.
706  // In this case the transformation is
707  //
708  // (x,y) = (u, u tan(pi/4 v))
709  //
710  // with Jacobian
711  //
712  // J = pi/4 R / cos(pi/4 v)
713  //
714  // And we get rid of R to take into account the singularity,
715  // unless specified differently in the constructor.
717  std::vector<double> & ws = this->weights;
718  double pi4 = numbers::PI / 4;
719
720  for (unsigned int q = 0; q < gauss.size(); ++q)
721  {
722  const Point<2> &gp = gauss.point(q);
723  ps[q][0] = gp[0];
724  ps[q][1] = gp[0] * std::tan(pi4 * gp[1]);
725  ws[q] = gauss.weight(q) * pi4 / std::cos(pi4 * gp[1]);
726  if (factor_out_singularity)
727  ws[q] *= (ps[q] - GeometryInfo<2>::unit_cell_vertex(0)).norm();
728  // The other half of the quadrilateral is symmetric with
729  // respect to xy plane.
730  ws[gauss.size() + q] = ws[q];
731  ps[gauss.size() + q][0] = ps[q][1];
732  ps[gauss.size() + q][1] = ps[q][0];
733  }
734
735  // Now we distribute these vertices in the correct manner
736  double theta = 0;
737  switch (vertex_index)
738  {
739  case 0:
740  theta = 0;
741  break;
742  case 1:
743  //
744  theta = numbers::PI / 2;
745  break;
746  case 2:
747  theta = -numbers::PI / 2;
748  break;
749  case 3:
750  theta = numbers::PI;
751  break;
752  }
753
754  double R00 = std::cos(theta), R01 = -std::sin(theta);
755  double R10 = std::sin(theta), R11 = std::cos(theta);
756
757  if (vertex_index != 0)
758  for (unsigned int q = 0; q < size(); ++q)
759  {
760  double x = ps[q][0] - .5, y = ps[q][1] - .5;
761
762  ps[q][0] = R00 * x + R01 * y + .5;
763  ps[q][1] = R10 * x + R11 * y + .5;
764  }
765 }
766
767
768 template <int dim>
771 {
773  for (unsigned int i = 0; i < quad.size(); ++i)
774  permutation[i] = i;
775
776  std::sort(permutation.begin(),
777  permutation.end(),
778  [this](const unsigned int x, const unsigned int y) {
779  return this->compare_weights(x, y);
780  });
781
782  // At this point, the variable is_tensor_product_flag is set
783  // to the respective value of the given Quadrature in the base
784  // class copy constructor.
785  // We only call a quadrature formula 'tensor product'
786  // if the quadrature points are also sorted lexicographically.
787  // In particular, any reordering destroys that property
788  // and we might need to modify the variable accordingly.
789  for (unsigned int i = 0; i < quad.size(); ++i)
790  {
793  if (permutation[i] != i)
794  this->is_tensor_product_flag = false;
795  }
796 }
797
798
799 template <int dim>
800 bool
801 QSorted<dim>::compare_weights(const unsigned int a, const unsigned int b) const
802 {
803  return (this->weights[a] < this->weights[b]);
804 }
805
806
807 // construct the quadrature formulae in higher dimensions by
808 // tensor product of lower dimensions
809
810 template <int dim>
811 QGauss<dim>::QGauss(const unsigned int n)
812  : Quadrature<dim>(QGauss<dim - 1>(n), QGauss<1>(n))
813 {}
814
815
816
817 template <int dim>
819  : Quadrature<dim>(QGaussLobatto<dim - 1>(n), QGaussLobatto<1>(n))
820 {}
821
822
823
824 template <int dim>
826  : Quadrature<dim>(QMidpoint<dim - 1>(), QMidpoint<1>())
827 {}
828
829
830
831 template <int dim>
833  : Quadrature<dim>(QTrapezoid<dim - 1>(), QTrapezoid<1>())
834 {}
835
836
837
838 template <int dim>
840  : Quadrature<dim>(QSimpson<dim - 1>(), QSimpson<1>())
841 {}
842
843
844
845 template <int dim>
847  : Quadrature<dim>(QMilne<dim - 1>(), QMilne<1>())
848 {}
849
850
851 template <int dim>
853  : Quadrature<dim>(QWeddle<dim - 1>(), QWeddle<1>())
854 {}
855
856 template <int dim>
858  const Point<dim> & singularity)
859  : // We need the explicit implementation if dim == 1. If dim > 1 we use the
860  // former implementation and apply a tensorial product to obtain the higher
861  // dimensions.
863  dim == 2 ?
866  dim == 3 ?
871 {}
872
873 template <int dim>
874 QTelles<dim>::QTelles(const unsigned int n, const Point<dim> &singularity)
875  : // In this case we map the standard Gauss Legendre formula using the given
876  // singularity point coordinates.
878 {}
879
880
881
882 template <>
884  : // We explicitly implement the Telles' variable change if dim == 1.
886 {
887  // We define all the constants to be used in the implementation of
888  // Telles' rule
889  const double eta_bar = singularity[0] * 2. - 1.;
890  const double eta_star = eta_bar * eta_bar - 1.;
891  double gamma_bar;
892
894  std::vector<double> weights_dummy(weights.size());
895  unsigned int cont = 0;
896  const double tol = 1e-10;
897  for (unsigned int d = 0; d < quadrature_points.size(); ++d)
898  {
899  if (std::abs(quadrature_points[d][0] - singularity[0]) > tol)
900  {
902  weights_dummy[d - cont] = weights[d];
903  }
904  else
905  {
906  // We need to remove the singularity point from the quadrature point
907  // list. To do so we use the variable cont.
908  cont = 1;
909  }
910  }
911  if (cont == 1)
912  {
914  weights.resize(weights_dummy.size() - 1);
915  for (unsigned int d = 0; d < quadrature_points.size(); ++d)
916  {
918  weights[d] = weights_dummy[d];
919  }
920  }
921  // We need to check if the singularity is at the boundary of the interval.
922  if (std::abs(eta_star) <= tol)
923  {
924  gamma_bar =
925  std::pow((eta_bar * eta_star + std::abs(eta_star)), 1.0 / 3.0) +
926  std::pow((eta_bar * eta_star - std::abs(eta_star)), 1.0 / 3.0) +
927  eta_bar;
928  }
929  else
930  {
931  gamma_bar = (eta_bar * eta_star + std::abs(eta_star)) /
932  std::abs(eta_bar * eta_star + std::abs(eta_star)) *
933  std::pow(std::abs(eta_bar * eta_star + std::abs(eta_star)),
934  1.0 / 3.0) +
935  (eta_bar * eta_star - std::abs(eta_star)) /
936  std::abs(eta_bar * eta_star - std::abs(eta_star)) *
937  std::pow(std::abs(eta_bar * eta_star - std::abs(eta_star)),
938  1.0 / 3.0) +
939  eta_bar;
940  }
941  for (unsigned int q = 0; q < quadrature_points.size(); ++q)
942  {
943  double gamma = quadrature_points[q][0] * 2 - 1;
944  double eta = (std::pow(gamma - gamma_bar, 3.0) +
945  gamma_bar * (gamma_bar * gamma_bar + 3)) /
946  (1 + 3 * gamma_bar * gamma_bar);
947
948  double J = 3 * ((gamma - gamma_bar) * (gamma - gamma_bar)) /
949  (1 + 3 * gamma_bar * gamma_bar);
950
951  quadrature_points[q][0] = (eta + 1) / 2.0;
952  weights[q] = J * weights[q];
953  }
954 }
955
956 namespace internal
957 {
958  namespace QGaussChebyshev
959  {
963  std::vector<double>
965  {
966  std::vector<double> points(n);
967  // n point quadrature: index from 0 to n-1
968  for (unsigned short i = 0; i < n; ++i)
969  // would be cos((2i+1)Pi/(2N+2))
970  // put + Pi so we start from the smallest point
971  // then map from [-1,1] to [0,1]
972  points[i] =
973  1. / 2. *
974  (1. + std::cos(numbers::PI *
975  (1. + double(2 * i + 1) / double(2 * (n - 1) + 2))));
976
977  return points;
978  }
979
980
981
985  std::vector<double>
987  {
988  std::vector<double> weights(n);
989
990  for (unsigned short i = 0; i < n; ++i)
991  {
992  // same weights as on [-1,1]
993  weights[i] = numbers::PI / double(n);
994  }
995
996  return weights;
997  }
998  } // namespace QGaussChebyshev
999 } // namespace internal
1000
1001
1002 template <>
1005 {
1006  Assert(n > 0, ExcMessage("Need at least one point for the quadrature rule"));
1009
1010  for (unsigned int i = 0; i < this->size(); ++i)
1011  {
1013  this->weights[i] = w[i];
1014  }
1015 }
1016
1017
1018 template <int dim>
1021 {}
1022
1023
1024 namespace internal
1025 {
1027  {
1028  // Computes the points of the quadrature formula.
1029  std::vector<double>
1032  {
1033  std::vector<double> points(n);
1034  // n point quadrature: index from 0 to n-1
1035  for (unsigned short i = 0; i < n; ++i)
1036  // would be -cos(2i Pi/(2N+1))
1037  // put + Pi so we start from the smallest point
1038  // then map from [-1,1] to [0,1]
1039  switch (ep)
1040  {
1042  {
1043  points[i] =
1044  1. / 2. *
1045  (1. -
1046  std::cos(numbers::PI *
1047  (1 + 2 * double(i) / (2 * double(n - 1) + 1.))));
1048  break;
1049  }
1050
1052  {
1053  points[i] =
1054  1. / 2. *
1055  (1. - std::cos(numbers::PI * (2 * double(n - 1 - i) /
1056  (2 * double(n - 1) + 1.))));
1057  break;
1058  }
1059
1060  default:
1061  Assert(
1062  false,
1063  ExcMessage(
1064  "This constructor can only be called with either "
1066  "second argument."));
1067  }
1068
1069  return points;
1070  }
1071
1072
1073
1074  // Computes the weights of the quadrature formula.
1075  std::vector<double>
1078  {
1079  std::vector<double> weights(n);
1080
1081  for (unsigned short i = 0; i < n; ++i)
1082  {
1083  // same weights as on [-1,1]
1084  weights[i] = 2. * numbers::PI / double(2 * (n - 1) + 1.);
1085  if (ep == ::QGaussRadauChebyshev<1>::left && i == 0)
1086  weights[i] /= 2.;
1087  else if (ep == ::QGaussRadauChebyshev<1>::right &&
1088  i == (n - 1))
1089  weights[i] /= 2.;
1090  }
1091
1092  return weights;
1093  }
1095 } // namespace internal
1096
1097
1098 template <>
1101  , ep(ep)
1102 {
1103  Assert(n > 0, ExcMessage("Need at least one point for quadrature rules"));
1104  std::vector<double> p =
1106  std::vector<double> w =
1108
1109  for (unsigned int i = 0; i < this->size(); ++i)
1110  {
1112  this->weights[i] = w[i];
1113  }
1114 }
1115
1116
1117 template <int dim>
1119  EndPoint ep)
1121  n,
1123  , ep(ep)
1124 {}
1125
1126
1127
1128 namespace internal
1129 {
1131  {
1132  // Computes the points of the quadrature formula.
1133  std::vector<double>
1135  {
1136  std::vector<double> points(n);
1137  // n point quadrature: index from 0 to n-1
1138  for (unsigned short i = 0; i < n; ++i)
1139  // would be cos(i Pi/N)
1140  // put + Pi so we start from the smallest point
1141  // then map from [-1,1] to [0,1]
1142  points[i] =
1143  1. / 2. *
1144  (1. + std::cos(numbers::PI * (1 + double(i) / double(n - 1))));
1145
1146  return points;
1147  }
1148
1149  // Computes the weights of the quadrature formula.
1150  std::vector<double>
1152  {
1153  std::vector<double> weights(n);
1154
1155  for (unsigned short i = 0; i < n; ++i)
1156  {
1157  // same weights as on [-1,1]
1158  weights[i] = numbers::PI / double((n - 1));
1159  if (i == 0 || i == (n - 1))
1160  weights[i] /= 2.;
1161  }
1162
1163  return weights;
1164  }
1165  } // namespace QGaussLobattoChebyshev
1166 } // namespace internal
1167
1168
1169
1170 template <>
1173 {
1174  Assert(n > 1,
1175  ExcMessage(
1176  "Need at least two points for Gauss-Lobatto quadrature rule"));
1177  std::vector<double> p =
1179  std::vector<double> w =
1181
1182  for (unsigned int i = 0; i < this->size(); ++i)
1183  {
1185  this->weights[i] = w[i];
1186  }
1187 }
1188
1189
1190 template <int dim>
1193 {}
1194
1195
1196
1197 template <int dim>
1199 {
1200  std::vector<Point<dim>> qpoints;
1201  std::vector<double> weights;
1202
1203  for (unsigned int i = 0; i < quad.size(); ++i)
1204  {
1205  double r = 0;
1206  for (unsigned int d = 0; d < dim; ++d)
1208  if (r <= 1 + 1e-10)
1209  {
1212  }
1213  }
1214 }
1215
1216
1217
1218 template <int dim>
1221  const std::array<Point<dim>, dim + 1> &vertices) const
1222 {
1223  Tensor<2, dim> B;
1224  for (unsigned int d = 0; d < dim; ++d)
1225  B[d] = vertices[d + 1] - vertices[0];
1226
1227  B = transpose(B);
1228  const double J = std::abs(determinant(B));
1229
1230  // if the determinant is zero, we return an empty quadrature
1231  if (J < 1e-12)
1233
1234  std::vector<Point<dim>> qp(this->size());
1235  std::vector<double> w(this->size());
1236
1237  for (unsigned int i = 0; i < this->size(); ++i)
1238  {
1239  qp[i] = Point<dim>(vertices[0] + B * this->point(i));
1240  w[i] = J * this->weight(i);
1241  }
1242
1244 }
1245
1246
1247
1251 {
1254  this->weights.resize(base.size());
1255  for (unsigned int i = 0; i < base.size(); ++i)
1256  {
1257  const auto q = base.point(i);
1258  const auto w = base.weight(i);
1259
1260  const auto xhat = q[0];
1261  const auto yhat = q[1];
1262
1263  const double t = numbers::PI_2 * yhat;
1264  const double pi = numbers::PI;
1265  const double st = std::sin(t);
1266  const double ct = std::cos(t);
1267  const double r = xhat / (st + ct);
1268
1269  const double J = pi * xhat / (2 * (std::sin(pi * yhat) + 1));
1270
1271  this->quadrature_points[i] = Point<2>(r * ct, r * st);
1272  this->weights[i] = w * J;
1273  }
1274 }
1275
1276
1277
1278 QTrianglePolar::QTrianglePolar(const unsigned int n)
1279  : QTrianglePolar(QGauss<1>(n), QGauss<1>(n))
1280 {}
1281
1282
1283
1286  const double beta)
1288 {
1291  this->weights.resize(base.size());
1292  for (unsigned int i = 0; i < base.size(); ++i)
1293  {
1294  const auto q = base.point(i);
1295  const auto w = base.weight(i);
1296
1297  const auto xhat = q[0];
1298  const auto yhat = q[1];
1299
1300  const double x = std::pow(xhat, beta) * (1 - yhat);
1301  const double y = std::pow(xhat, beta) * yhat;
1302
1303  const double J = beta * std::pow(xhat, 2. * beta - 1.);
1304
1306  this->weights[i] = w * J;
1307  }
1308 }
1309
1310
1311
1312 QDuffy::QDuffy(const unsigned int n, const double beta)
1313  : QDuffy(QGauss<1>(n), QGauss<1>(n), beta)
1314 {}
1315
1316
1317
1318 template <int dim>
1320 {
1322  ExcMessage(
1323  "The split point should be inside the unit reference cell."));
1324
1325  std::array<Point<dim>, dim + 1> vertices;
1326  vertices[0] = split_point;
1327
1328  // Make a simplex from the split_point and the first dim vertices of each
1329  // face. In dimension three, we need to split the face in two triangles, so
1330  // we use once the first dim vertices of each face, and the second time the
1331  // the dim vertices of each face starting from 1.
1332  for (auto f : GeometryInfo<dim>::face_indices())
1333  for (unsigned int start = 0; start < (dim > 2 ? 2 : 1); ++start)
1334  {
1335  for (unsigned int i = 0; i < dim; ++i)
1336  vertices[i + 1] = GeometryInfo<dim>::unit_cell_vertex(
1338  const auto quad = base.compute_affine_transformation(vertices);
1340  {
1344  this->weights.insert(this->weights.end(),
1347  }
1348  }
1349 }
1350
1351
1352
1353 // explicit specialization
1354 // note that 1d formulae are specialized by implementation above
1355 template class QGauss<2>;
1356 template class QGaussLobatto<2>;
1357 template class QMidpoint<2>;
1358 template class QTrapezoid<2>;
1359 template class QSimpson<2>;
1360 template class QMilne<2>;
1361 template class QWeddle<2>;
1362
1363 template class QGauss<3>;
1364 template class QGaussLobatto<3>;
1365 template class QMidpoint<3>;
1366 template class QTrapezoid<3>;
1367 template class QSimpson<3>;
1368 template class QMilne<3>;
1369 template class QWeddle<3>;
1370
1371 template class QSorted<1>;
1372 template class QSorted<2>;
1373 template class QSorted<3>;
1374
1375 template class QTelles<1>;
1376 template class QTelles<2>;
1377 template class QTelles<3>;
1378
1379 template class QGaussChebyshev<1>;
1380 template class QGaussChebyshev<2>;
1381 template class QGaussChebyshev<3>;
1382
1386
1387 template class QGaussLobattoChebyshev<1>;
1388 template class QGaussLobattoChebyshev<2>;
1389 template class QGaussLobattoChebyshev<3>;
1390
1391 template class QSimplex<1>;
1392 template class QSimplex<2>;
1393 template class QSimplex<3>;
1394
1395 template class QSplit<1>;
1396 template class QSplit<2>;
1397 template class QSplit<3>;
1398
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
Number * begin_raw()
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)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1691
const double fraction
QGaussChebyshev(const unsigned int n)
Generate a formula with n quadrature points.
const Point< dim > & point(const unsigned int i) const
std::vector< double > get_quadrature_weights(const unsigned int n)
static std::vector< double > get_quadrature_points(const unsigned int n)
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:472
static Point< dim > unit_cell_vertex(const unsigned int vertex)
std::vector< double > get_quadrature_points(const unsigned int n)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1576
QGauss(const unsigned int n)
void split_point(const Point< dim1+dim2 > &source, Point< dim1 > &p1, Point< dim2 > &p2)
Number jacobi_polynomial_value(const unsigned int degree, const int alpha, const int beta, const Number x)
Definition: polynomial.h:972
Number * end_raw()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1466
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
static std::vector< double > get_quadrature_weights(const unsigned int n)
Generate a formula with n quadrature points.
std::vector< double > get_quadrature_weights(const unsigned int n)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:372
bool compare_weights(const unsigned int a, const unsigned int b) const
Point< 3 > vertices[4]
QGaussLobatto(const unsigned int n)
static constexpr double PI_2
Definition: numbers.h:236
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
bool is_tensor_product_flag
std::vector< Point< dim > > quadrature_points
unsigned int size() const
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
std::vector< long double > compute_quadrature_weights(const std::vector< long double > &x, const int alpha, const int beta)
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Definition: tensor.h:448
static constexpr double PI
Definition: numbers.h:231
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:371
QGaussLobattoChebyshev(const unsigned int n)
Generate a formula with n quadrature points.
static ::ExceptionBase & ExcNotImplemented()