Reference documentation for deal.II version 9.3.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
quadrature_lib.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2021 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
19 
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
39  Quadrature<0>(1)
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
52  Quadrature<0>(1)
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)
63  : Quadrature<1>(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  {
73  this->quadrature_points[i][0] = points[i];
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>
114  compute_quadrature_weights(const std::vector<long double> &x,
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)
142  : Quadrature<1>(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  {
156  this->quadrature_points[i][0] = points[i];
157  this->weights[i] = 0.5 * w[i];
158  }
159 }
160 #endif
161 
162 
163 template <>
165  : Quadrature<1>(1)
166 {
167  this->quadrature_points[0] = Point<1>(0.5);
168  this->weights[0] = 1.0;
169 }
170 
171 
172 
173 template <>
175  : Quadrature<1>(2)
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  {
182  this->quadrature_points[i] = Point<1>(xpts[i]);
183  this->weights[i] = wts[i];
184  }
185 }
186 
187 
188 
189 template <>
191  : Quadrature<1>(3)
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  {
198  this->quadrature_points[i] = Point<1>(xpts[i]);
199  this->weights[i] = wts[i];
200  }
201 }
202 
203 
204 
205 template <>
207  : Quadrature<1>(5)
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  {
215  this->quadrature_points[i] = Point<1>(xpts[i]);
216  this->weights[i] = wts[i];
217  }
218 }
219 
220 
221 
222 template <>
224  : Quadrature<1>(7)
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  {
238  this->quadrature_points[i] = Point<1>(xpts[i]);
239  this->weights[i] = wts[i];
240  }
241 }
242 
243 
244 template <>
245 QGaussLog<1>::QGaussLog(const unsigned int n, const bool revert)
246  : Quadrature<1>(n)
247 {
248  std::vector<double> p = get_quadrature_points(n);
249  std::vector<double> w = get_quadrature_weights(n);
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|.
256  this->quadrature_points[i] =
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 {
399  std::vector<double> quadrature_weights(n);
400 
401  switch (n)
402  {
403  case 1:
404  quadrature_weights[0] = -1.0;
405  break;
406  case 2:
407  quadrature_weights[0] = -0.7185393190303845;
408  quadrature_weights[1] = -0.2814606809696154;
409  break;
410 
411  case 3:
412  quadrature_weights[0] = -0.5134045522323634;
413  quadrature_weights[1] = -0.3919800412014877;
414  quadrature_weights[2] = -0.0946154065661483;
415  break;
416 
417  case 4:
418  quadrature_weights[0] = -0.3834640681451353;
419  quadrature_weights[1] = -0.3868753177747627;
420  quadrature_weights[2] = -0.1904351269501432;
421  quadrature_weights[3] = -0.03922548712995894;
422  break;
423 
424  case 5:
425  quadrature_weights[0] = -0.2978934717828955;
426  quadrature_weights[1] = -0.3497762265132236;
427  quadrature_weights[2] = -0.234488290044052;
428  quadrature_weights[3] = -0.0989304595166356;
429  quadrature_weights[4] = -0.01891155214319462;
430  break;
431 
432  case 6:
433  quadrature_weights[0] = -0.2387636625785478;
434  quadrature_weights[1] = -0.3082865732739458;
435  quadrature_weights[2] = -0.2453174265632108;
436  quadrature_weights[3] = -0.1420087565664786;
437  quadrature_weights[4] = -0.05545462232488041;
438  quadrature_weights[5] = -0.01016895869293513;
439  break;
440 
441 
442  case 7:
443  quadrature_weights[0] = -0.1961693894252476;
444  quadrature_weights[1] = -0.2703026442472726;
445  quadrature_weights[2] = -0.239681873007687;
446  quadrature_weights[3] = -0.1657757748104267;
447  quadrature_weights[4] = -0.0889432271377365;
448  quadrature_weights[5] = -0.03319430435645653;
449  quadrature_weights[6] = -0.005932787015162054;
450  break;
451 
452  case 8:
453  quadrature_weights[0] = -0.164416604728002;
454  quadrature_weights[1] = -0.2375256100233057;
455  quadrature_weights[2] = -0.2268419844319134;
456  quadrature_weights[3] = -0.1757540790060772;
457  quadrature_weights[4] = -0.1129240302467932;
458  quadrature_weights[5] = -0.05787221071771947;
459  quadrature_weights[6] = -0.02097907374214317;
460  quadrature_weights[7] = -0.003686407104036044;
461  break;
462 
463  case 9:
464  quadrature_weights[0] = -0.1400684387481339;
465  quadrature_weights[1] = -0.2097722052010308;
466  quadrature_weights[2] = -0.211427149896601;
467  quadrature_weights[3] = -0.1771562339380667;
468  quadrature_weights[4] = -0.1277992280331758;
469  quadrature_weights[5] = -0.07847890261203835;
470  quadrature_weights[6] = -0.0390225049841783;
471  quadrature_weights[7] = -0.01386729555074604;
472  quadrature_weights[8] = -0.002408041036090773;
473  break;
474 
475  case 10:
476  quadrature_weights[0] = -0.12095513195457;
477  quadrature_weights[1] = -0.1863635425640733;
478  quadrature_weights[2] = -0.1956608732777627;
479  quadrature_weights[3] = -0.1735771421828997;
480  quadrature_weights[4] = -0.135695672995467;
481  quadrature_weights[5] = -0.0936467585378491;
482  quadrature_weights[6] = -0.05578772735275126;
483  quadrature_weights[7] = -0.02715981089692378;
484  quadrature_weights[8] = -0.00951518260454442;
485  quadrature_weights[9] = -0.001638157633217673;
486  break;
487 
488 
489  case 11:
490  quadrature_weights[0] = -0.1056522560990997;
491  quadrature_weights[1] = -0.1665716806006314;
492  quadrature_weights[2] = -0.1805632182877528;
493  quadrature_weights[3] = -0.1672787367737502;
494  quadrature_weights[4] = -0.1386970574017174;
495  quadrature_weights[5] = -0.1038334333650771;
496  quadrature_weights[6] = -0.06953669788988512;
497  quadrature_weights[7] = -0.04054160079499477;
498  quadrature_weights[8] = -0.01943540249522013;
499  quadrature_weights[9] = -0.006737429326043388;
500  quadrature_weights[10] = -0.001152486965101561;
501  break;
502 
503  case 12:
504  quadrature_weights[0] = -0.09319269144393;
505  quadrature_weights[1] = -0.1497518275763289;
506  quadrature_weights[2] = -0.166557454364573;
507  quadrature_weights[3] = -0.1596335594369941;
508  quadrature_weights[4] = -0.1384248318647479;
509  quadrature_weights[5] = -0.1100165706360573;
510  quadrature_weights[6] = -0.07996182177673273;
511  quadrature_weights[7] = -0.0524069547809709;
512  quadrature_weights[8] = -0.03007108900074863;
513  quadrature_weights[9] = -0.01424924540252916;
514  quadrature_weights[10] = -0.004899924710875609;
515  quadrature_weights[11] = -0.000834029009809656;
516  break;
517 
518  default:
519  Assert(false, ExcNotImplemented());
520  break;
521  }
522 
523  return quadrature_weights;
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)
532  : Quadrature<1>(
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);
550  QGaussLog<1> quad2(n);
551  QGauss<1> quad(n);
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.
565  this->quadrature_points[i] = quad1.point(i) * fraction;
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  {
571  this->quadrature_points[j] = quad.point(i) * fraction;
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  {
578  this->quadrature_points[i + n] =
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|
583  this->quadrature_points[j + n] =
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(
593  this->quadrature_points[i] != origin,
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)
633  : Quadrature<2>(quad_size(singularity, n))
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.
641  std::vector<QGaussOneOverR<2>> quads;
642  std::vector<Point<2>> origins;
643  // Id of the corner with a
644  // singularity
645  quads.emplace_back(n, 3, factor_out_singularity);
646  quads.emplace_back(n, 2, factor_out_singularity);
647  quads.emplace_back(n, 1, factor_out_singularity);
648  quads.emplace_back(n, 0, factor_out_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);
670  this->quadrature_points[q_id] =
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 
692  // Start with the gauss quadrature formula on the (u,v) reference
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.
716  std::vector<Point<2>> &ps = this->quadrature_points;
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>
770  : Quadrature<dim>(quad)
771 {
772  std::vector<unsigned int> permutation(quad.size());
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  {
791  this->weights[i] = quad.weight(permutation[i]);
792  this->quadrature_points[i] = quad.point(permutation[i]);
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.
862  Quadrature<dim>(
863  dim == 2 ?
864  QAnisotropic<dim>(QTelles<1>(base_quad, Point<1>(singularity[0])),
865  QTelles<1>(base_quad, Point<1>(singularity[1]))) :
866  dim == 3 ?
867  QAnisotropic<dim>(QTelles<1>(base_quad, Point<1>(singularity[0])),
868  QTelles<1>(base_quad, Point<1>(singularity[1])),
869  QTelles<1>(base_quad, Point<1>(singularity[2]))) :
870  Quadrature<dim>())
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.
877  Quadrature<dim>(QTelles<dim>(QGauss<1>(n), singularity))
878 {}
879 
880 
881 
882 template <>
883 QTelles<1>::QTelles(const Quadrature<1> &base_quad, const Point<1> &singularity)
884  : // We explicitly implement the Telles' variable change if dim == 1.
885  Quadrature<1>(base_quad)
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 
893  std::vector<Point<1>> quadrature_points_dummy(quadrature_points.size());
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  {
901  quadrature_points_dummy[d - cont] = quadrature_points[d];
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  {
913  quadrature_points.resize(quadrature_points_dummy.size() - 1);
914  weights.resize(weights_dummy.size() - 1);
915  for (unsigned int d = 0; d < quadrature_points.size(); ++d)
916  {
917  quadrature_points[d] = quadrature_points_dummy[d];
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>
964  get_quadrature_points(const unsigned int n)
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>
986  get_quadrature_weights(const unsigned int n)
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 <>
1004  : Quadrature<1>(n)
1005 {
1006  Assert(n > 0, ExcMessage("Need at least one point for the quadrature rule"));
1007  std::vector<double> p = internal::QGaussChebyshev::get_quadrature_points(n);
1008  std::vector<double> w = internal::QGaussChebyshev::get_quadrature_weights(n);
1009 
1010  for (unsigned int i = 0; i < this->size(); ++i)
1011  {
1012  this->quadrature_points[i] = Point<1>(p[i]);
1013  this->weights[i] = w[i];
1014  }
1015 }
1016 
1017 
1018 template <int dim>
1020  : Quadrature<dim>(QGaussChebyshev<1>(n))
1021 {}
1022 
1023 
1024 namespace internal
1025 {
1027  {
1028  // Computes the points of the quadrature formula.
1029  std::vector<double>
1030  get_quadrature_points(const unsigned int n,
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  {
1041  case ::QGaussRadauChebyshev<1>::left:
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 
1051  case ::QGaussRadauChebyshev<1>::right:
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 "
1065  "QGaussRadauChebyshev::left or QGaussRadauChebyshev::right as "
1066  "second argument."));
1067  }
1068 
1069  return points;
1070  }
1071 
1072 
1073 
1074  // Computes the weights of the quadrature formula.
1075  std::vector<double>
1076  get_quadrature_weights(const unsigned int n,
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  }
1094  } // namespace QGaussRadauChebyshev
1095 } // namespace internal
1096 
1097 
1098 template <>
1100  : Quadrature<1>(n)
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  {
1111  this->quadrature_points[i] = Point<1>(p[i]);
1112  this->weights[i] = w[i];
1113  }
1114 }
1115 
1116 
1117 template <int dim>
1119  EndPoint ep)
1120  : Quadrature<dim>(QGaussRadauChebyshev<1>(
1121  n,
1122  static_cast<QGaussRadauChebyshev<1>::EndPoint>(ep)))
1123  , ep(ep)
1124 {}
1125 
1126 
1127 
1128 namespace internal
1129 {
1131  {
1132  // Computes the points of the quadrature formula.
1133  std::vector<double>
1134  get_quadrature_points(const unsigned int n)
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>
1151  get_quadrature_weights(const unsigned int n)
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 <>
1172  : Quadrature<1>(n)
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  {
1184  this->quadrature_points[i] = Point<1>(p[i]);
1185  this->weights[i] = w[i];
1186  }
1187 }
1188 
1189 
1190 template <int dim>
1192  : Quadrature<dim>(QGaussLobattoChebyshev<1>(n))
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  /* Use "int d" instead of the more natural "unsigned int d" to work
1207  * around a wrong diagnostic in gcc-10.3.0 that warns about that the
1208  * comparison "d < dim" is always false in case of "dim == 0".
1209  * MM 2021 */
1210  for (int d = 0; d < dim; ++d)
1211  r += quad.point(i)[d];
1212  if (r <= 1 + 1e-10)
1213  {
1214  this->quadrature_points.push_back(quad.point(i));
1215  this->weights.push_back(quad.weight(i));
1216  }
1217  }
1218 }
1219 
1220 
1221 
1222 template <int dim>
1225  const std::array<Point<dim>, dim + 1> &vertices) const
1226 {
1227  Tensor<2, dim> B;
1228  for (unsigned int d = 0; d < dim; ++d)
1229  B[d] = vertices[d + 1] - vertices[0];
1230 
1231  B = transpose(B);
1232  const double J = std::abs(determinant(B));
1233 
1234  // if the determinant is zero, we return an empty quadrature
1235  if (J < 1e-12)
1236  return Quadrature<dim>();
1237 
1238  std::vector<Point<dim>> qp(this->size());
1239  std::vector<double> w(this->size());
1240 
1241  for (unsigned int i = 0; i < this->size(); ++i)
1242  {
1243  qp[i] = Point<dim>(vertices[0] + B * this->point(i));
1244  w[i] = J * this->weight(i);
1245  }
1246 
1247  return Quadrature<dim>(qp, w);
1248 }
1249 
1250 
1251 
1253  const Quadrature<1> &angular_quadrature)
1254  : QSimplex<2>(Quadrature<2>())
1255 {
1256  const QAnisotropic<2> base(radial_quadrature, angular_quadrature);
1257  this->quadrature_points.resize(base.size());
1258  this->weights.resize(base.size());
1259  for (unsigned int i = 0; i < base.size(); ++i)
1260  {
1261  const auto &q = base.point(i);
1262  const auto w = base.weight(i);
1263 
1264  const auto xhat = q[0];
1265  const auto yhat = q[1];
1266 
1267  const double t = numbers::PI_2 * yhat;
1268  const double pi = numbers::PI;
1269  const double st = std::sin(t);
1270  const double ct = std::cos(t);
1271  const double r = xhat / (st + ct);
1272 
1273  const double J = pi * xhat / (2 * (std::sin(pi * yhat) + 1));
1274 
1275  this->quadrature_points[i] = Point<2>(r * ct, r * st);
1276  this->weights[i] = w * J;
1277  }
1278 }
1279 
1280 
1281 
1282 QTrianglePolar::QTrianglePolar(const unsigned int n)
1283  : QTrianglePolar(QGauss<1>(n), QGauss<1>(n))
1284 {}
1285 
1286 
1287 
1288 QDuffy::QDuffy(const Quadrature<1> &radial_quadrature,
1289  const Quadrature<1> &angular_quadrature,
1290  const double beta)
1291  : QSimplex<2>(Quadrature<2>())
1292 {
1293  const QAnisotropic<2> base(radial_quadrature, angular_quadrature);
1294  this->quadrature_points.resize(base.size());
1295  this->weights.resize(base.size());
1296  for (unsigned int i = 0; i < base.size(); ++i)
1297  {
1298  const auto &q = base.point(i);
1299  const auto w = base.weight(i);
1300 
1301  const auto xhat = q[0];
1302  const auto yhat = q[1];
1303 
1304  const double x = std::pow(xhat, beta) * (1 - yhat);
1305  const double y = std::pow(xhat, beta) * yhat;
1306 
1307  const double J = beta * std::pow(xhat, 2. * beta - 1.);
1308 
1309  this->quadrature_points[i] = Point<2>(x, y);
1310  this->weights[i] = w * J;
1311  }
1312 }
1313 
1314 
1315 
1316 QDuffy::QDuffy(const unsigned int n, const double beta)
1317  : QDuffy(QGauss<1>(n), QGauss<1>(n), beta)
1318 {}
1319 
1320 
1321 
1322 template <int dim>
1324 {
1326  ExcMessage(
1327  "The split point should be inside the unit reference cell."));
1328 
1329  std::array<Point<dim>, dim + 1> vertices;
1330  vertices[0] = split_point;
1331 
1332  // Make a simplex from the split_point and the first dim vertices of each
1333  // face. In dimension three, we need to split the face in two triangles, so
1334  // we use once the first dim vertices of each face, and the second time the
1335  // the dim vertices of each face starting from 1.
1336  for (auto f : GeometryInfo<dim>::face_indices())
1337  for (unsigned int start = 0; start < (dim > 2 ? 2 : 1); ++start)
1338  {
1339  for (unsigned int i = 0; i < dim; ++i)
1340  vertices[i + 1] = GeometryInfo<dim>::unit_cell_vertex(
1342  const auto quad = base.compute_affine_transformation(vertices);
1343  if (quad.size())
1344  {
1345  this->quadrature_points.insert(this->quadrature_points.end(),
1346  quad.get_points().begin(),
1347  quad.get_points().end());
1348  this->weights.insert(this->weights.end(),
1349  quad.get_weights().begin(),
1350  quad.get_weights().end());
1351  }
1352  }
1353 }
1354 
1355 
1356 
1357 template <int dim>
1358 QGaussSimplex<dim>::QGaussSimplex(const unsigned int n_points_1D)
1359  : QSimplex<dim>(Quadrature<dim>())
1360 {
1361  // fill quadrature points and quadrature weights
1362  if (dim == 0 || dim == 1)
1363  {
1364  const ::QGauss<dim> quad(n_points_1D);
1365 
1366  this->quadrature_points = quad.get_points();
1367  this->weights = quad.get_weights();
1368  }
1369  else if (dim == 2)
1370  {
1371  if (n_points_1D == 1)
1372  {
1373  const double p = 1.0 / 3.0;
1374  this->quadrature_points.emplace_back(p, p);
1375  this->weights.emplace_back(0.5);
1376  }
1377  else if (n_points_1D == 2)
1378  {
1379  const double Q23 = 2.0 / 3.0;
1380  const double Q16 = 1.0 / 6.0;
1381 
1382  this->quadrature_points.emplace_back(Q23, Q16);
1383  this->quadrature_points.emplace_back(Q16, Q23);
1384  this->quadrature_points.emplace_back(Q16, Q16);
1385  this->weights.emplace_back(Q16);
1386  this->weights.emplace_back(Q16);
1387  this->weights.emplace_back(Q16);
1388  }
1389  else if (n_points_1D == 3)
1390  {
1391  const double q12 = 0.5;
1392 
1393  // clang-format off
1394  this->quadrature_points.emplace_back(0.3333333333330, 0.3333333333330);
1395  this->quadrature_points.emplace_back(0.7974269853530, 0.1012865073230);
1396  this->quadrature_points.emplace_back(0.1012865073230, 0.7974269853530);
1397  this->quadrature_points.emplace_back(0.1012865073230, 0.1012865073230);
1398  this->quadrature_points.emplace_back(0.0597158717898, 0.4701420641050);
1399  this->quadrature_points.emplace_back(0.4701420641050, 0.0597158717898);
1400  this->quadrature_points.emplace_back(0.4701420641050, 0.4701420641050);
1401  // clang-format on
1402 
1403  this->weights.emplace_back(q12 * 0.225);
1404  this->weights.emplace_back(q12 * 0.125939180545);
1405  this->weights.emplace_back(q12 * 0.125939180545);
1406  this->weights.emplace_back(q12 * 0.125939180545);
1407  this->weights.emplace_back(q12 * 0.132394152789);
1408  this->weights.emplace_back(q12 * 0.132394152789);
1409  this->weights.emplace_back(q12 * 0.132394152789);
1410  }
1411  else if (n_points_1D == 4)
1412  {
1414  QWitherdenVincentSimplex<dim>(n_points_1D));
1415  }
1416  }
1417  else if (dim == 3)
1418  {
1419  if (n_points_1D == 1)
1420  {
1421  const double Q14 = 1.0 / 4.0;
1422  const double Q16 = 1.0 / 6.0;
1423 
1424  this->quadrature_points.emplace_back(Q14, Q14, Q14);
1425  this->weights.emplace_back(Q16);
1426  }
1427  else if (n_points_1D == 2)
1428  {
1429  const double Q124 = 1.0 / 6.0 / 4.0;
1430 
1431  const double palpha = (5.0 + 3.0 * sqrt(5.0)) / 20.0;
1432  const double pbeta = (5.0 - sqrt(5.0)) / 20.0;
1433  this->quadrature_points.emplace_back(pbeta, pbeta, pbeta);
1434  this->quadrature_points.emplace_back(palpha, pbeta, pbeta);
1435  this->quadrature_points.emplace_back(pbeta, palpha, pbeta);
1436  this->quadrature_points.emplace_back(pbeta, pbeta, palpha);
1437  this->weights.emplace_back(Q124);
1438  this->weights.emplace_back(Q124);
1439  this->weights.emplace_back(Q124);
1440  this->weights.emplace_back(Q124);
1441  }
1442  else if (n_points_1D == 3)
1443  {
1444  const double Q16 = 1.0 / 6.0;
1445 
1446  // clang-format off
1447  this->quadrature_points.emplace_back(0.5684305841968444, 0.1438564719343852, 0.1438564719343852);
1448  this->quadrature_points.emplace_back(0.1438564719343852, 0.1438564719343852, 0.1438564719343852);
1449  this->quadrature_points.emplace_back(0.1438564719343852, 0.1438564719343852, 0.5684305841968444);
1450  this->quadrature_points.emplace_back(0.1438564719343852, 0.5684305841968444, 0.1438564719343852);
1451  this->quadrature_points.emplace_back(0.0000000000000000, 0.5000000000000000, 0.5000000000000000);
1452  this->quadrature_points.emplace_back(0.5000000000000000, 0.0000000000000000, 0.5000000000000000);
1453  this->quadrature_points.emplace_back(0.5000000000000000, 0.5000000000000000, 0.0000000000000000);
1454  this->quadrature_points.emplace_back(0.5000000000000000, 0.0000000000000000, 0.0000000000000000);
1455  this->quadrature_points.emplace_back(0.0000000000000000, 0.5000000000000000, 0.0000000000000000);
1456  this->quadrature_points.emplace_back(0.0000000000000000, 0.0000000000000000, 0.5000000000000000);
1457  // clang-format on
1458 
1459  this->weights.emplace_back(0.2177650698804054 * Q16);
1460  this->weights.emplace_back(0.2177650698804054 * Q16);
1461  this->weights.emplace_back(0.2177650698804054 * Q16);
1462  this->weights.emplace_back(0.2177650698804054 * Q16);
1463  this->weights.emplace_back(0.0214899534130631 * Q16);
1464  this->weights.emplace_back(0.0214899534130631 * Q16);
1465  this->weights.emplace_back(0.0214899534130631 * Q16);
1466  this->weights.emplace_back(0.0214899534130631 * Q16);
1467  this->weights.emplace_back(0.0214899534130631 * Q16);
1468  this->weights.emplace_back(0.0214899534130631 * Q16);
1469  }
1470  else if (n_points_1D == 4)
1471  {
1473  QWitherdenVincentSimplex<dim>(n_points_1D));
1474  }
1475  }
1476 
1477  AssertDimension(this->quadrature_points.size(), this->weights.size());
1478  Assert(this->quadrature_points.size() > 0,
1480  "QGaussSimplex is currently only implemented for "
1481  "n_points_1D = 1, 2, 3, and 4 while you are asking for "
1482  "n_points_1D = " +
1483  Utilities::to_string(n_points_1D)));
1484 }
1485 
1486 namespace
1487 {
1488  template <std::size_t b_dim>
1489  std::vector<std::array<double, b_dim>>
1490  all_permutations(const std::array<double, b_dim> &b_point)
1491  {
1492  std::vector<std::array<double, b_dim>> output;
1493 
1494  // We want all possible permutations of the barycentric coordinates.
1495  // The easiest way to get all of them is to sort the input first and
1496  // then use next_permutation to cycle through them all.
1497  std::array<double, b_dim> temp = b_point;
1498  std::sort(temp.begin(), temp.end());
1499  do
1500  {
1501  output.push_back(temp);
1502  }
1503  while (std::next_permutation(temp.begin(), temp.end()));
1504 
1505  return output;
1506  }
1507 } // namespace
1508 
1509 
1510 
1511 template <int dim>
1513  const unsigned int n_points_1D)
1514  : QSimplex<dim>(Quadrature<dim>())
1515 {
1516  Assert(1 <= dim && dim <= 3, ExcNotImplemented());
1517  // Just use Gauss in 1D: this is a high-order open rule so this is a
1518  // reasonable equivalent for generic programming.
1519  if (dim == 1)
1520  {
1521  Quadrature<dim>::operator=(::QGauss<dim>(n_points_1D));
1522  return;
1523  }
1524 
1525  std::array<double, dim + 1> centroid;
1526  std::fill(centroid.begin(), centroid.end(), 1.0 / (dim + 1.0));
1527  std::vector<std::vector<std::array<double, dim + 1>>> b_point_permutations;
1528  std::vector<double> b_weights;
1529 
1530  // We can simplify the implementation of these quadrature rules
1531  // by quite a bit by exploiting symmetry - we do essentially the
1532  // same thing for each barycentric coordinate, so we can express
1533  // our quadrature rule as permutations of barycentric points
1534  // instead of writing things out explicitly.
1535 
1536  // Apply a Barycentric permutation where one point is different.
1537  auto process_point_1 = [&](const double a, const double w) {
1538  const double b = 1.0 - dim * a;
1539  std::array<double, dim + 1> b_point;
1540  std::fill(b_point.begin(), b_point.begin() + dim, a);
1541  b_point[dim] = b;
1542 
1543  b_weights.push_back(w);
1544  b_point_permutations.push_back(all_permutations(b_point));
1545  };
1546 
1547  // Apply a Barycentric permutation where two points (in 3D) are different.
1548  auto process_point_2 = [&](const double a, const double w) {
1549  Assert(dim == 3, ExcInternalError());
1550  const double b = (1.0 - 2.0 * a) / 2.0;
1551  std::array<double, dim + 1> b_point;
1552  std::fill(b_point.begin(), b_point.begin() + dim - 1, a);
1553  b_point[dim - 1] = b;
1554  b_point[dim] = b;
1555 
1556  b_weights.push_back(w);
1557  b_point_permutations.push_back(all_permutations(b_point));
1558  };
1559 
1560  // Apply a Barycentric permutation where three (or four) points
1561  // are different (since there are two inputs).
1562  auto process_point_3 = [&](const double a, const double b, const double w) {
1563  const double c = 1.0 - (dim - 1.0) * a - b;
1564  std::array<double, dim + 1> b_point;
1565  std::fill(b_point.begin(), b_point.begin() + dim - 1, a);
1566  b_point[dim - 1] = b;
1567  b_point[dim] = c;
1568 
1569  b_weights.push_back(w);
1570  b_point_permutations.push_back(all_permutations(b_point));
1571  };
1572 
1573  if (n_points_1D == 1)
1574  {
1575  b_point_permutations.push_back({centroid});
1576  b_weights.push_back(1.0);
1577  }
1578  else if (n_points_1D == 2)
1579  {
1580  // This is WV-4 in 2D and WV-3 in 3D
1581  if (dim == 2)
1582  {
1583  process_point_1(9.1576213509770743e-02, 1.0995174365532187e-01);
1584  process_point_1(4.4594849091596489e-01, 2.2338158967801147e-01);
1585  }
1586  else if (dim == 3)
1587  {
1588  process_point_1(3.281633025163817e-01, 1.362178425370874e-01);
1589  process_point_1(1.080472498984286e-01, 1.137821574629126e-01);
1590  }
1591  }
1592  else if (n_points_1D == 3)
1593  {
1594  // This is the WV-5 rule in both 2D and 3D
1595  if (dim == 2)
1596  {
1597  b_weights.push_back(0.225);
1598  b_point_permutations.push_back({centroid});
1599 
1600  process_point_1(1.0128650732345634e-01, 1.2593918054482714e-01);
1601  process_point_1(4.7014206410511511e-01, 1.3239415278850619e-01);
1602  }
1603  else if (dim == 3)
1604  {
1605  process_point_1(3.108859192633006e-01, 1.126879257180159e-01);
1606  process_point_1(9.273525031089125e-02, 7.349304311636196e-02);
1607 
1608  process_point_2(4.550370412564964e-02, 4.254602077708147e-02);
1609  }
1610  }
1611  else if (n_points_1D == 4)
1612  {
1613  // This is the WV-7 rule in both 2D and 3D
1614  if (dim == 2)
1615  {
1616  process_point_1(3.3730648554587850e-02, 1.6545050110792131e-02);
1617  process_point_1(4.7430969250471822e-01, 7.7086646185986069e-02);
1618  process_point_1(2.4157738259540357e-01, 1.2794417123015558e-01);
1619  process_point_3(4.7036644652595216e-02,
1620  1.9868331479735168e-01,
1621  5.5878732903199779e-02);
1622  }
1623  else if (dim == 3)
1624  {
1625  b_point_permutations.push_back({centroid});
1626  b_weights.push_back(9.548528946413085e-02);
1627 
1628  process_point_1(3.157011497782028e-01, 4.232958120996703e-02);
1629  process_point_2(5.048982259839635e-02, 3.189692783285758e-02);
1630 
1631  process_point_3(1.888338310260010e-01,
1632  5.751716375870000e-01,
1633  3.720713072833462e-02);
1634  process_point_3(2.126547254148314e-02,
1635  8.108302410985486e-01,
1636  8.110770829903342e-03);
1637  }
1638  }
1639  else if (n_points_1D == 5)
1640  {
1641  // This is the WV-9 rule in both 2D and 3D
1642  if (dim == 2)
1643  {
1644  b_point_permutations.push_back({centroid});
1645  b_weights.push_back(9.7135796282798836e-02);
1646 
1647  process_point_1(4.4729513394452691e-02, 2.5577675658698031e-02);
1648  process_point_1(4.8968251919873762e-01, 3.1334700227139071e-02);
1649  process_point_1(4.3708959149293664e-01, 7.7827541004774278e-02);
1650  process_point_1(1.8820353561903275e-01, 7.9647738927210249e-02);
1651 
1652  process_point_3(3.6838412054736258e-02,
1653  2.2196298916076568e-01,
1654  4.3283539377289376e-02);
1655  }
1656  else if (dim == 3)
1657  {
1658  b_point_permutations.push_back({centroid});
1659  b_weights.push_back(5.801054891248025e-02);
1660 
1661  process_point_1(6.198169755222693e-10, 6.431928175925639e-05);
1662  process_point_1(1.607745353952616e-01, 2.317333846242546e-02);
1663  process_point_1(3.222765218214210e-01, 2.956291233542929e-02);
1664  process_point_1(4.510891834541358e-02, 8.063979979616182e-03);
1665 
1666  process_point_2(1.122965460043761e-01, 3.813408010370246e-02);
1667 
1668  process_point_3(4.588714487524592e-01,
1669  2.554579233041310e-03,
1670  8.384422198298552e-03);
1671  process_point_3(3.377587068533860e-02,
1672  7.183503264420745e-01,
1673  1.023455935274533e-02);
1674  process_point_3(1.836413698099279e-01,
1675  3.441591057817528e-02,
1676  2.052491596798814e-02);
1677  }
1678  }
1679  else if (n_points_1D == 6)
1680  {
1681  // There is no WV-11 rule in 3D yet
1682  if (dim == 2)
1683  {
1684  b_point_permutations.push_back({centroid});
1685  b_weights.push_back(8.5761179732224219e-02);
1686 
1687  process_point_1(2.8485417614371900e-02, 1.0431870512894697e-02);
1688  process_point_1(4.9589190096589092e-01, 1.6606273054585369e-02);
1689  process_point_1(1.0263548271224643e-01, 3.8630759237019321e-02);
1690  process_point_1(4.3846592676435220e-01, 6.7316154079468296e-02);
1691  process_point_1(2.1021995670317828e-01, 7.0515684111716576e-02);
1692 
1693  process_point_3(7.3254276860644785e-03,
1694  1.4932478865208237e-01,
1695  1.0290289572953278e-02);
1696  process_point_3(4.6010500165429957e-02,
1697  2.8958112563770588e-01,
1698  4.0332476640500554e-02);
1699  }
1700  else if (dim == 3)
1701  {
1702  Assert(false, ExcNotImplemented());
1703  }
1704  }
1705  else
1706  {
1707  Assert(false, ExcNotImplemented());
1708  }
1709 
1710  Assert(b_point_permutations.size() == b_weights.size(), ExcInternalError());
1711  for (unsigned int permutation_n = 0; permutation_n < b_weights.size();
1712  ++permutation_n)
1713  {
1714  for (const std::array<double, dim + 1> &b_point :
1715  b_point_permutations[permutation_n])
1716  {
1717  const double volume = (dim == 2 ? 1.0 / 2.0 : 1.0 / 6.0);
1718  this->weights.emplace_back(volume * b_weights[permutation_n]);
1719  Point<dim> c_point;
1720  std::copy(b_point.begin(),
1721  b_point.begin() + dim,
1722  c_point.begin_raw());
1723  this->quadrature_points.emplace_back(c_point);
1724  }
1725  }
1726 }
1727 
1728 
1729 
1730 template <int dim>
1731 QGaussWedge<dim>::QGaussWedge(const unsigned int n_points)
1732  : Quadrature<dim>()
1733 {
1734  AssertDimension(dim, 3);
1735 
1736  QGaussSimplex<2> quad_tri(n_points);
1737  QGauss<1> quad_line(n_points);
1738 
1739  for (unsigned int i = 0; i < quad_line.size(); ++i)
1740  for (unsigned int j = 0; j < quad_tri.size(); ++j)
1741  {
1742  this->quadrature_points.emplace_back(quad_tri.point(j)[0],
1743  quad_tri.point(j)[1],
1744  quad_line.point(i)[0]);
1745  this->weights.emplace_back(quad_tri.weight(j) * quad_line.weight(i));
1746  }
1747 
1748  AssertDimension(this->quadrature_points.size(), this->weights.size());
1749  Assert(this->quadrature_points.size() > 0,
1750  ExcMessage("No valid quadrature points!"));
1751 }
1752 
1753 
1754 
1755 template <int dim>
1756 QGaussPyramid<dim>::QGaussPyramid(const unsigned int n_points_1D)
1757  : Quadrature<dim>()
1758 {
1759  AssertDimension(dim, 3);
1760 
1761  if (n_points_1D == 1)
1762  {
1763  const double Q14 = 1.0 / 4.0;
1764  const double Q43 = 4.0 / 3.0;
1765 
1766  this->quadrature_points.emplace_back(0, 0, Q14);
1767  this->weights.emplace_back(Q43);
1768  }
1769  else if (n_points_1D == 2)
1770  {
1771  // clang-format off
1772  this->quadrature_points.emplace_back(-0.26318405556971, -0.26318405556971, 0.54415184401122);
1773  this->quadrature_points.emplace_back(-0.50661630334979, -0.50661630334979, 0.12251482265544);
1774  this->quadrature_points.emplace_back(-0.26318405556971, +0.26318405556971, 0.54415184401122);
1775  this->quadrature_points.emplace_back(-0.50661630334979, +0.50661630334979, 0.12251482265544);
1776  this->quadrature_points.emplace_back(+0.26318405556971, -0.26318405556971, 0.54415184401122);
1777  this->quadrature_points.emplace_back(+0.50661630334979, -0.50661630334979, 0.12251482265544);
1778  this->quadrature_points.emplace_back(+0.26318405556971, +0.26318405556971, 0.54415184401122);
1779  this->quadrature_points.emplace_back(+0.50661630334979, +0.50661630334979, 0.12251482265544);
1780  // clang-format on
1781 
1782  this->weights.emplace_back(0.10078588207983);
1783  this->weights.emplace_back(0.23254745125351);
1784  this->weights.emplace_back(0.10078588207983);
1785  this->weights.emplace_back(0.23254745125351);
1786  this->weights.emplace_back(0.10078588207983);
1787  this->weights.emplace_back(0.23254745125351);
1788  this->weights.emplace_back(0.10078588207983);
1789  this->weights.emplace_back(0.23254745125351);
1790  }
1791 
1792  AssertDimension(this->quadrature_points.size(), this->weights.size());
1793  Assert(this->quadrature_points.size() > 0,
1794  ExcMessage("No valid quadrature points!"));
1795 }
1796 
1797 
1798 
1799 // explicit specialization
1800 // note that 1d formulae are specialized by implementation above
1801 template class QGauss<2>;
1802 template class QGaussLobatto<2>;
1803 template class QMidpoint<2>;
1804 template class QTrapezoid<2>;
1805 template class QSimpson<2>;
1806 template class QMilne<2>;
1807 template class QWeddle<2>;
1808 
1809 template class QGauss<3>;
1810 template class QGaussLobatto<3>;
1811 template class QMidpoint<3>;
1812 template class QTrapezoid<3>;
1813 template class QSimpson<3>;
1814 template class QMilne<3>;
1815 template class QWeddle<3>;
1816 
1817 template class QSorted<1>;
1818 template class QSorted<2>;
1819 template class QSorted<3>;
1820 
1821 template class QTelles<1>;
1822 template class QTelles<2>;
1823 template class QTelles<3>;
1824 
1825 template class QGaussChebyshev<1>;
1826 template class QGaussChebyshev<2>;
1827 template class QGaussChebyshev<3>;
1828 
1829 template class QGaussRadauChebyshev<1>;
1830 template class QGaussRadauChebyshev<2>;
1831 template class QGaussRadauChebyshev<3>;
1832 
1833 template class QGaussLobattoChebyshev<1>;
1834 template class QGaussLobattoChebyshev<2>;
1835 template class QGaussLobattoChebyshev<3>;
1836 
1837 template class QSimplex<1>;
1838 template class QSimplex<2>;
1839 template class QSimplex<3>;
1840 
1841 template class QSplit<1>;
1842 template class QSplit<2>;
1843 template class QSplit<3>;
1844 
1845 template class QGaussSimplex<0>;
1846 template class QGaussSimplex<1>;
1847 template class QGaussSimplex<2>;
1848 template class QGaussSimplex<3>;
1849 template class QGaussWedge<1>;
1850 template class QGaussWedge<2>;
1851 template class QGaussWedge<3>;
1852 template class QGaussPyramid<1>;
1853 template class QGaussPyramid<2>;
1854 template class QGaussPyramid<3>;
1855 
1856 template class QWitherdenVincentSimplex<1>;
1857 template class QWitherdenVincentSimplex<2>;
1858 template class QWitherdenVincentSimplex<3>;
1859 
QGaussLogR(const unsigned int n, const Point< dim > &x0=Point< dim >(), const double alpha=1, const bool factor_out_singular_weight=false)
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
Number * begin_raw()
QGaussSimplex(const unsigned int n_points_1D)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
QGaussLog(const unsigned int n, const bool revert=false)
std::vector< double > weights
Definition: quadrature.h:289
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)
QGaussWedge(const unsigned int n_points_1D)
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
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)
QSimplex(const Quadrature< dim > &quad)
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:1575
QGauss(const unsigned int n)
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:482
Quadrature & operator=(const Quadrature< dim > &)
Definition: quadrature.cc:282
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:976
Number * end_raw()
QGaussPyramid(const unsigned int n_points_1D)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1465
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
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.
QDuffy(const Quadrature< 1 > &radial_quadrature, const Quadrature< 1 > &angular_quadrature, const double beta=1.0)
std::vector< double > get_quadrature_weights(const unsigned int n)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:395
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
QWitherdenVincentSimplex(const unsigned int n_points_1D)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
bool is_tensor_product_flag
Definition: quadrature.h:298
std::vector< Point< dim > > quadrature_points
Definition: quadrature.h:283
unsigned int size() const
static unsigned int quad_size(const Point< dim > &singularity, const unsigned int n)
QTelles(const Quadrature< 1 > &base_quad, const Point< dim > &singularity)
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)
std::vector< double > get_quadrature_points(const unsigned int n, ::QGaussRadauChebyshev< 1 >::EndPoint ep)
Definition: tensor.h:462
static constexpr double PI
Definition: numbers.h:231
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:394
double volume(const Triangulation< dim, spacedim > &tria, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
Definition: grid_tools.cc:137
QGaussLobattoChebyshev(const unsigned int n)
Generate a formula with n quadrature points.
static ::ExceptionBase & ExcNotImplemented()
QSorted(const Quadrature< dim > &quad)
long double gamma(const unsigned int n)
QSplit(const QSimplex< dim > &base, const Point< dim > &split_point)
std::vector< double > get_quadrature_weights(const unsigned int n, ::QGaussRadauChebyshev< 1 >::EndPoint ep)
void copy(const T *begin, const T *end, U *dest)
Quadrature< dim > compute_affine_transformation(const std::array< Point< dim >, dim+1 > &vertices) const
std::vector< double > get_quadrature_points(const unsigned int n)
double weight(const unsigned int i) const
static ::ExceptionBase & ExcInternalError()
QTrianglePolar(const Quadrature< 1 > &radial_quadrature, const Quadrature< 1 > &angular_quadrature)