Reference documentation for deal.II version GIT 9042b9283b 2023-12-02 14:50:02+00:00
\(\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 - 2023 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
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 #include <deal.II/base/utilities.h>
20 
21 #include <deal.II/fe/fe_nothing.h>
22 #include <deal.II/fe/fe_values.h>
23 
26 #include <deal.II/grid/tria.h>
27 
28 #include <algorithm>
29 #include <cmath>
30 #include <functional>
31 #include <limits>
32 
33 
35 
36 
37 // please note: for a given dimension, we need the quadrature formulae
38 // for all lower dimensions as well. That is why in this file the check
39 // is for deal_II_dimension >= any_number and not for ==
40 
41 
42 
43 template <>
44 QGauss<0>::QGauss(const unsigned int)
45  : // there are n_q^dim == 1
46  // points
47  Quadrature<0>(1)
48 {
49  // the single quadrature point gets unit
50  // weight
51  this->weights[0] = 1;
52 }
53 
54 
55 
56 template <>
58  : // there are n_q^dim == 1
59  // points
60  Quadrature<0>(1)
61 {
62  // the single quadrature point gets unit
63  // weight
64  this->weights[0] = 1;
65 }
66 
67 
68 
69 template <>
70 QGauss<1>::QGauss(const unsigned int n)
71  : Quadrature<1>(n)
72 {
73  if (n == 0)
74  return;
75 
76  std::vector<long double> points =
77  Polynomials::jacobi_polynomial_roots<long double>(n, 0, 0);
78 
79  for (unsigned int i = 0; i < (points.size() + 1) / 2; ++i)
80  {
81  this->quadrature_points[i][0] = points[i];
82  this->quadrature_points[n - i - 1][0] = 1. - points[i];
83 
84  // derivative of Jacobi polynomial
85  const long double pp =
86  0.5 * (n + 1) *
87  Polynomials::jacobi_polynomial_value(n - 1, 1, 1, points[i]);
88  const long double x = -1. + 2. * points[i];
89  const double w = 1. / ((1. - x * x) * pp * pp);
90  this->weights[i] = w;
91  this->weights[n - i - 1] = w;
92  }
93 }
94 
95 namespace internal
96 {
97  namespace QGaussLobatto
98  {
103  long double
104  gamma(const unsigned int n)
105  {
106  long double result = n - 1;
107  for (int i = n - 2; i > 1; --i)
108  result *= i;
109  return result;
110  }
111 
112 
113 
121  std::vector<long double>
122  compute_quadrature_weights(const std::vector<long double> &x,
123  const int alpha,
124  const int beta)
125  {
126  const unsigned int q = x.size();
127  std::vector<long double> w(q);
128 
129  const long double factor =
130  std::pow(2., alpha + beta + 1) * gamma(alpha + q) * gamma(beta + q) /
131  ((q - 1) * gamma(q) * gamma(alpha + beta + q + 1));
132  for (unsigned int i = 0; i < q; ++i)
133  {
134  const long double s =
135  Polynomials::jacobi_polynomial_value(q - 1, alpha, beta, x[i]);
136  w[i] = factor / (s * s);
137  }
138  w[0] *= (beta + 1);
139  w[q - 1] *= (alpha + 1);
140 
141  return w;
142  }
143  } // namespace QGaussLobatto
144 } // namespace internal
145 
146 
147 namespace internal
148 {
149  namespace QGaussRadau
150  {
151 
152  // Implements lookup table after affine transformation to [0,1].
153  //
154  // Analytical values for [-1,1] and n < 4 listed on
155  // https://mathworld.wolfram.com/RadauQuadrature.html
156  // Values for n > 3 calculated with the Julia Package
157  // FastGaussQuadrature.jl
158  // https://github.com/JuliaApproximation/FastGaussQuadrature.jl
159  //
160  std::vector<double>
161  get_left_quadrature_points(const unsigned int n)
162  {
163  std::vector<double> q_points(n);
164  switch (n)
165  {
166  case 1:
167  q_points[0] = 0.;
168  break;
169  case 2:
170  q_points[0] = 0.;
171  q_points[1] = 2. / 3.;
172  break;
173  case 3:
174  q_points[0] = 0.;
175  q_points[1] = (6. - std::sqrt(6)) * 0.1;
176  q_points[2] = (6. + std::sqrt(6)) * 0.1;
177  break;
178 
179  case 4:
180  q_points[0] = 0.000000000000000000;
181  q_points[1] = 0.212340538239152943;
182  q_points[2] = 0.590533135559265343;
183  q_points[3] = 0.911412040487296071;
184  break;
185  case 5:
186  q_points[0] = 0.000000000000000000;
187  q_points[1] = 0.139759864343780571;
188  q_points[2] = 0.416409567631083166;
189  q_points[3] = 0.723156986361876197;
190  q_points[4] = 0.942895803885482331;
191  break;
192  case 6:
193  q_points[0] = 0.000000000000000000;
194  q_points[1] = 0.098535085798826416;
195  q_points[2] = 0.304535726646363913;
196  q_points[3] = 0.562025189752613841;
197  q_points[4] = 0.801986582126391845;
198  q_points[5] = 0.960190142948531222;
199  break;
200  case 7:
201  q_points[0] = 0.000000000000000000;
202  q_points[1] = 0.073054328680258851;
203  q_points[2] = 0.230766137969945495;
204  q_points[3] = 0.441328481228449865;
205  q_points[4] = 0.663015309718845702;
206  q_points[5] = 0.851921400331515644;
207  q_points[6] = 0.970683572840215114;
208  break;
209  case 8:
210  q_points[0] = 0.000000000000000000;
211  q_points[1] = 0.056262560536922135;
212  q_points[2] = 0.180240691736892389;
213  q_points[3] = 0.352624717113169672;
214  q_points[4] = 0.547153626330555420;
215  q_points[5] = 0.734210177215410598;
216  q_points[6] = 0.885320946839095790;
217  q_points[7] = 0.977520613561287499;
218  break;
219  default:
220  Assert(false, ExcNotImplemented());
221  break;
222  }
223  return q_points;
224  }
225 
226  std::vector<double>
227  get_quadrature_points(const unsigned int n,
228  const ::QGaussRadau<1>::EndPoint end_point)
229  {
230  std::vector<double> left_points = get_left_quadrature_points(n);
231  switch (end_point)
232  {
233  case ::QGaussRadau<1>::EndPoint::left:
234  return left_points;
235  case ::QGaussRadau<1>::EndPoint::right:
236  {
237  std::vector<double> points(n);
238  for (unsigned int i = 0; i < n; ++i)
239  {
240  points[n - i - 1] = 1. - left_points[i];
241  }
242  return points;
243  }
244  default:
245  Assert(
246  false,
247  ExcMessage(
248  "This constructor can only be called with either "
249  "QGaussRadau::left or QGaussRadau::right as second argument."));
250  return {};
251  }
252  }
253 
254  // Implements lookup table after affine transformation to [0,1].
255  //
256  // Analytical values for [-1,1] and n < 4 listed on
257  // https://mathworld.wolfram.com/RadauQuadrature.html
258  // Values for n > 3 calculated with the Julia Package
259  // FastGaussQuadrature.jl
260  // https://github.com/JuliaApproximation/FastGaussQuadrature.jl
261  //
262  std::vector<double>
263  get_left_quadrature_weights(const unsigned int n)
264  {
265  std::vector<double> weights(n);
266  switch (n)
267  {
268  case 1:
269  weights[0] = 1.;
270  break;
271  case 2:
272  weights[0] = 0.25;
273  weights[1] = 0.75;
274  break;
275  case 3:
276  weights[0] = 1. / 9.;
277  weights[1] = (16. + std::sqrt(6)) / 36.;
278  weights[2] = (16. - std::sqrt(6)) / 36.;
279  break;
280  case 4:
281  weights[0] = 0.062500000000000000;
282  weights[1] = 0.328844319980059696;
283  weights[2] = 0.388193468843171852;
284  weights[3] = 0.220462211176768369;
285  break;
286  case 5:
287  weights[0] = 0.040000000000000001;
288  weights[1] = 0.223103901083570894;
289  weights[2] = 0.311826522975741427;
290  weights[3] = 0.281356015149462124;
291  weights[4] = 0.143713560791225797;
292  break;
293  case 6:
294  weights[0] = 0.027777777777777776;
295  weights[1] = 0.159820376610255471;
296  weights[2] = 0.242693594234484888;
297  weights[3] = 0.260463391594787597;
298  weights[4] = 0.208450667155953895;
299  weights[5] = 0.100794192626740456;
300  break;
301  case 7:
302  weights[0] = 0.020408163265306121;
303  weights[1] = 0.119613744612656100;
304  weights[2] = 0.190474936822115581;
305  weights[3] = 0.223554914507283209;
306  weights[4] = 0.212351889502977870;
307  weights[5] = 0.159102115733650767;
308  weights[6] = 0.074494235556010341;
309  break;
310  case 8:
311  weights[0] = 0.015625000000000000;
312  weights[1] = 0.092679077401489660;
313  weights[2] = 0.152065310323392683;
314  weights[3] = 0.188258772694559262;
315  weights[4] = 0.195786083726246729;
316  weights[5] = 0.173507397817250691;
317  weights[6] = 0.124823950664932445;
318  weights[7] = 0.057254407372128648;
319  break;
320 
321  default:
323  break;
324  }
325  return weights;
326  }
327 
328  std::vector<double>
329  get_quadrature_weights(const unsigned int n,
330  const ::QGaussRadau<1>::EndPoint end_point)
331  {
332  std::vector<double> left_weights = get_left_quadrature_weights(n);
333  switch (end_point)
334  {
335  case ::QGaussRadau<1>::EndPoint::left:
336  return left_weights;
337  case ::QGaussRadau<1>::EndPoint::right:
338  {
339  std::vector<double> weights(n);
340  for (unsigned int i = 0; i < n; ++i)
341  {
342  weights[n - i - 1] = left_weights[i];
343  }
344  return weights;
345  }
346  default:
347  Assert(false,
348  ExcMessage(
349  "This constructor can only be called with either "
350  "QGaussRadau::EndPoint::left or "
351  "QGaussRadau::EndPoint::right as second argument."));
352  return {};
353  }
354  }
355  } // namespace QGaussRadau
356 } // namespace internal
357 
358 #ifndef DOXYGEN
359 template <>
360 QGaussRadau<1>::QGaussRadau(const unsigned int n, const EndPoint end_point)
361  : Quadrature<1>(n)
362  , end_point(end_point)
363 {
364  Assert(n > 0, ExcMessage("Need at least one point for quadrature rules."));
365  std::vector<double> p =
367  std::vector<double> w =
369 
370  for (unsigned int i = 0; i < this->size(); ++i)
371  {
372  this->quadrature_points[i] = ::Point<1>(p[i]);
373  this->weights[i] = w[i];
374  }
375 }
376 #endif
377 
378 
379 #ifndef DOXYGEN
380 template <>
381 QGaussLobatto<1>::QGaussLobatto(const unsigned int n)
382  : Quadrature<1>(n)
383 {
384  Assert(n >= 2, ExcNotImplemented());
385 
386  std::vector<long double> points =
387  Polynomials::jacobi_polynomial_roots<long double>(n - 2, 1, 1);
388  points.insert(points.begin(), 0);
389  points.push_back(1.);
390  std::vector<long double> w =
392 
393  // scale weights to the interval [0.0, 1.0]:
394  for (unsigned int i = 0; i < points.size(); ++i)
395  {
396  this->quadrature_points[i][0] = points[i];
397  this->weights[i] = 0.5 * w[i];
398  }
399 }
400 #endif
401 
402 
403 template <>
405  : Quadrature<1>(1)
406 {
407  this->quadrature_points[0] = Point<1>(0.5);
408  this->weights[0] = 1.0;
409 }
410 
411 
412 
413 template <>
415  : Quadrature<1>(2)
416 {
417  static const double xpts[] = {0.0, 1.0};
418  static const double wts[] = {0.5, 0.5};
419 
420  for (unsigned int i = 0; i < this->size(); ++i)
421  {
422  this->quadrature_points[i] = Point<1>(xpts[i]);
423  this->weights[i] = wts[i];
424  }
425 }
426 
427 
428 
429 template <>
431  : Quadrature<1>(3)
432 {
433  static const double xpts[] = {0.0, 0.5, 1.0};
434  static const double wts[] = {1. / 6., 2. / 3., 1. / 6.};
435 
436  for (unsigned int i = 0; i < this->size(); ++i)
437  {
438  this->quadrature_points[i] = Point<1>(xpts[i]);
439  this->weights[i] = wts[i];
440  }
441 }
442 
443 
444 
445 template <>
447  : Quadrature<1>(5)
448 {
449  static const double xpts[] = {0.0, .25, .5, .75, 1.0};
450  static const double wts[] = {
451  7. / 90., 32. / 90., 12. / 90., 32. / 90., 7. / 90.};
452 
453  for (unsigned int i = 0; i < this->size(); ++i)
454  {
455  this->quadrature_points[i] = Point<1>(xpts[i]);
456  this->weights[i] = wts[i];
457  }
458 }
459 
460 
461 
462 template <>
464  : Quadrature<1>(7)
465 {
466  static const double xpts[] = {
467  0.0, 1. / 6., 1. / 3., .5, 2. / 3., 5. / 6., 1.0};
468  static const double wts[] = {41. / 840.,
469  216. / 840.,
470  27. / 840.,
471  272. / 840.,
472  27. / 840.,
473  216. / 840.,
474  41. / 840.};
475 
476  for (unsigned int i = 0; i < this->size(); ++i)
477  {
478  this->quadrature_points[i] = Point<1>(xpts[i]);
479  this->weights[i] = wts[i];
480  }
481 }
482 
483 
484 template <>
485 QGaussLog<1>::QGaussLog(const unsigned int n, const bool revert)
486  : Quadrature<1>(n)
487 {
488  std::vector<double> p = get_quadrature_points(n);
489  std::vector<double> w = get_quadrature_weights(n);
490 
491  for (unsigned int i = 0; i < this->size(); ++i)
492  {
493  // Using the change of variables x=1-t, it's possible to show
494  // that int f(x)ln|1-x| = int f(1-t) ln|t|, which implies that
495  // we can use this quadrature formula also with weight ln|1-x|.
496  this->quadrature_points[i] =
497  revert ? Point<1>(1 - p[n - 1 - i]) : Point<1>(p[i]);
498  this->weights[i] = revert ? w[n - 1 - i] : w[i];
499  }
500 }
501 
502 template <>
503 std::vector<double>
505 {
506  std::vector<double> q_points(n);
507 
508  switch (n)
509  {
510  case 1:
511  q_points[0] = 0.3333333333333333;
512  break;
513 
514  case 2:
515  q_points[0] = 0.1120088061669761;
516  q_points[1] = 0.6022769081187381;
517  break;
518 
519  case 3:
520  q_points[0] = 0.06389079308732544;
521  q_points[1] = 0.3689970637156184;
522  q_points[2] = 0.766880303938942;
523  break;
524 
525  case 4:
526  q_points[0] = 0.04144848019938324;
527  q_points[1] = 0.2452749143206022;
528  q_points[2] = 0.5561654535602751;
529  q_points[3] = 0.848982394532986;
530  break;
531 
532  case 5:
533  q_points[0] = 0.02913447215197205;
534  q_points[1] = 0.1739772133208974;
535  q_points[2] = 0.4117025202849029;
536  q_points[3] = 0.6773141745828183;
537  q_points[4] = 0.89477136103101;
538  break;
539 
540  case 6:
541  q_points[0] = 0.02163400584411693;
542  q_points[1] = 0.1295833911549506;
543  q_points[2] = 0.3140204499147661;
544  q_points[3] = 0.5386572173517997;
545  q_points[4] = 0.7569153373774084;
546  q_points[5] = 0.922668851372116;
547  break;
548 
549 
550  case 7:
551  q_points[0] = 0.0167193554082585;
552  q_points[1] = 0.100185677915675;
553  q_points[2] = 0.2462942462079286;
554  q_points[3] = 0.4334634932570557;
555  q_points[4] = 0.6323509880476823;
556  q_points[5] = 0.81111862674023;
557  q_points[6] = 0.940848166743287;
558  break;
559 
560  case 8:
561  q_points[0] = 0.01332024416089244;
562  q_points[1] = 0.07975042901389491;
563  q_points[2] = 0.1978710293261864;
564  q_points[3] = 0.354153994351925;
565  q_points[4] = 0.5294585752348643;
566  q_points[5] = 0.7018145299391673;
567  q_points[6] = 0.849379320441094;
568  q_points[7] = 0.953326450056343;
569  break;
570 
571  case 9:
572  q_points[0] = 0.01086933608417545;
573  q_points[1] = 0.06498366633800794;
574  q_points[2] = 0.1622293980238825;
575  q_points[3] = 0.2937499039716641;
576  q_points[4] = 0.4466318819056009;
577  q_points[5] = 0.6054816627755208;
578  q_points[6] = 0.7541101371585467;
579  q_points[7] = 0.877265828834263;
580  q_points[8] = 0.96225055941096;
581  break;
582 
583  case 10:
584  q_points[0] = 0.00904263096219963;
585  q_points[1] = 0.05397126622250072;
586  q_points[2] = 0.1353118246392511;
587  q_points[3] = 0.2470524162871565;
588  q_points[4] = 0.3802125396092744;
589  q_points[5] = 0.5237923179723384;
590  q_points[6] = 0.6657752055148032;
591  q_points[7] = 0.7941904160147613;
592  q_points[8] = 0.898161091216429;
593  q_points[9] = 0.9688479887196;
594  break;
595 
596 
597  case 11:
598  q_points[0] = 0.007643941174637681;
599  q_points[1] = 0.04554182825657903;
600  q_points[2] = 0.1145222974551244;
601  q_points[3] = 0.2103785812270227;
602  q_points[4] = 0.3266955532217897;
603  q_points[5] = 0.4554532469286375;
604  q_points[6] = 0.5876483563573721;
605  q_points[7] = 0.7139638500230458;
606  q_points[8] = 0.825453217777127;
607  q_points[9] = 0.914193921640008;
608  q_points[10] = 0.973860256264123;
609  break;
610 
611  case 12:
612  q_points[0] = 0.006548722279080035;
613  q_points[1] = 0.03894680956045022;
614  q_points[2] = 0.0981502631060046;
615  q_points[3] = 0.1811385815906331;
616  q_points[4] = 0.2832200676673157;
617  q_points[5] = 0.398434435164983;
618  q_points[6] = 0.5199526267791299;
619  q_points[7] = 0.6405109167754819;
620  q_points[8] = 0.7528650118926111;
621  q_points[9] = 0.850240024421055;
622  q_points[10] = 0.926749682988251;
623  q_points[11] = 0.977756129778486;
624  break;
625 
626  default:
627  Assert(false, ExcNotImplemented());
628  break;
629  }
630 
631  return q_points;
632 }
633 
634 
635 template <>
636 std::vector<double>
638 {
639  std::vector<double> quadrature_weights(n);
640 
641  switch (n)
642  {
643  case 1:
644  quadrature_weights[0] = -1.0;
645  break;
646  case 2:
647  quadrature_weights[0] = -0.7185393190303845;
648  quadrature_weights[1] = -0.2814606809696154;
649  break;
650 
651  case 3:
652  quadrature_weights[0] = -0.5134045522323634;
653  quadrature_weights[1] = -0.3919800412014877;
654  quadrature_weights[2] = -0.0946154065661483;
655  break;
656 
657  case 4:
658  quadrature_weights[0] = -0.3834640681451353;
659  quadrature_weights[1] = -0.3868753177747627;
660  quadrature_weights[2] = -0.1904351269501432;
661  quadrature_weights[3] = -0.03922548712995894;
662  break;
663 
664  case 5:
665  quadrature_weights[0] = -0.2978934717828955;
666  quadrature_weights[1] = -0.3497762265132236;
667  quadrature_weights[2] = -0.234488290044052;
668  quadrature_weights[3] = -0.0989304595166356;
669  quadrature_weights[4] = -0.01891155214319462;
670  break;
671 
672  case 6:
673  quadrature_weights[0] = -0.2387636625785478;
674  quadrature_weights[1] = -0.3082865732739458;
675  quadrature_weights[2] = -0.2453174265632108;
676  quadrature_weights[3] = -0.1420087565664786;
677  quadrature_weights[4] = -0.05545462232488041;
678  quadrature_weights[5] = -0.01016895869293513;
679  break;
680 
681 
682  case 7:
683  quadrature_weights[0] = -0.1961693894252476;
684  quadrature_weights[1] = -0.2703026442472726;
685  quadrature_weights[2] = -0.239681873007687;
686  quadrature_weights[3] = -0.1657757748104267;
687  quadrature_weights[4] = -0.0889432271377365;
688  quadrature_weights[5] = -0.03319430435645653;
689  quadrature_weights[6] = -0.005932787015162054;
690  break;
691 
692  case 8:
693  quadrature_weights[0] = -0.164416604728002;
694  quadrature_weights[1] = -0.2375256100233057;
695  quadrature_weights[2] = -0.2268419844319134;
696  quadrature_weights[3] = -0.1757540790060772;
697  quadrature_weights[4] = -0.1129240302467932;
698  quadrature_weights[5] = -0.05787221071771947;
699  quadrature_weights[6] = -0.02097907374214317;
700  quadrature_weights[7] = -0.003686407104036044;
701  break;
702 
703  case 9:
704  quadrature_weights[0] = -0.1400684387481339;
705  quadrature_weights[1] = -0.2097722052010308;
706  quadrature_weights[2] = -0.211427149896601;
707  quadrature_weights[3] = -0.1771562339380667;
708  quadrature_weights[4] = -0.1277992280331758;
709  quadrature_weights[5] = -0.07847890261203835;
710  quadrature_weights[6] = -0.0390225049841783;
711  quadrature_weights[7] = -0.01386729555074604;
712  quadrature_weights[8] = -0.002408041036090773;
713  break;
714 
715  case 10:
716  quadrature_weights[0] = -0.12095513195457;
717  quadrature_weights[1] = -0.1863635425640733;
718  quadrature_weights[2] = -0.1956608732777627;
719  quadrature_weights[3] = -0.1735771421828997;
720  quadrature_weights[4] = -0.135695672995467;
721  quadrature_weights[5] = -0.0936467585378491;
722  quadrature_weights[6] = -0.05578772735275126;
723  quadrature_weights[7] = -0.02715981089692378;
724  quadrature_weights[8] = -0.00951518260454442;
725  quadrature_weights[9] = -0.001638157633217673;
726  break;
727 
728 
729  case 11:
730  quadrature_weights[0] = -0.1056522560990997;
731  quadrature_weights[1] = -0.1665716806006314;
732  quadrature_weights[2] = -0.1805632182877528;
733  quadrature_weights[3] = -0.1672787367737502;
734  quadrature_weights[4] = -0.1386970574017174;
735  quadrature_weights[5] = -0.1038334333650771;
736  quadrature_weights[6] = -0.06953669788988512;
737  quadrature_weights[7] = -0.04054160079499477;
738  quadrature_weights[8] = -0.01943540249522013;
739  quadrature_weights[9] = -0.006737429326043388;
740  quadrature_weights[10] = -0.001152486965101561;
741  break;
742 
743  case 12:
744  quadrature_weights[0] = -0.09319269144393;
745  quadrature_weights[1] = -0.1497518275763289;
746  quadrature_weights[2] = -0.166557454364573;
747  quadrature_weights[3] = -0.1596335594369941;
748  quadrature_weights[4] = -0.1384248318647479;
749  quadrature_weights[5] = -0.1100165706360573;
750  quadrature_weights[6] = -0.07996182177673273;
751  quadrature_weights[7] = -0.0524069547809709;
752  quadrature_weights[8] = -0.03007108900074863;
753  quadrature_weights[9] = -0.01424924540252916;
754  quadrature_weights[10] = -0.004899924710875609;
755  quadrature_weights[11] = -0.000834029009809656;
756  break;
757 
758  default:
759  Assert(false, ExcNotImplemented());
760  break;
761  }
762 
763  return quadrature_weights;
764 }
765 
766 
767 template <>
768 QGaussLogR<1>::QGaussLogR(const unsigned int n,
769  const Point<1> &origin,
770  const double alpha,
771  const bool factor_out_singularity)
772  : Quadrature<1>(
773  ((origin[0] == 0) || (origin[0] == 1)) ? (alpha == 1 ? n : 2 * n) : 4 * n)
774  , fraction(((origin[0] == 0) || (origin[0] == 1.)) ? 1. : origin[0])
775 {
776  // The three quadrature formulas that make this one up. There are
777  // at most two when the origin is one of the extremes, and there is
778  // only one if the origin is one of the extremes and alpha is
779  // equal to one.
780  //
781  // If alpha is different from one, then we need a correction which
782  // is performed with a standard Gauss quadrature rule on each
783  // segment. This is not needed in the standard case where alpha is
784  // equal to one and the origin is on one of the extremes. We
785  // integrate with weight ln|(x-o)/alpha|. In the easy cases, we
786  // only need n quadrature points. In the most difficult one, we
787  // need 2*n points for the first segment, and 2*n points for the
788  // second segment.
789  QGaussLog<1> quad1(n, origin[0] != 0);
790  QGaussLog<1> quad2(n);
791  QGauss<1> quad(n);
792 
793  // Check that the origin is inside 0,1
794  Assert((fraction >= 0) && (fraction <= 1),
795  ExcMessage("Origin is outside [0,1]."));
796 
797  // Non singular offset. This is the start of non singular quad
798  // points.
799  unsigned int ns_offset = (fraction == 1) ? n : 2 * n;
800 
801  for (unsigned int i = 0, j = ns_offset; i < n; ++i, ++j)
802  {
803  // The first i quadrature points are the same as quad1, and
804  // are by default singular.
805  this->quadrature_points[i] = quad1.point(i) * fraction;
806  this->weights[i] = quad1.weight(i) * fraction;
807 
808  // We need to scale with -log|fraction*alpha|
809  if ((alpha != 1) || (fraction != 1))
810  {
811  this->quadrature_points[j] = quad.point(i) * fraction;
812  this->weights[j] =
813  -std::log(alpha / fraction) * quad.weight(i) * fraction;
814  }
815  // In case we need the second quadrature as well, do it now.
816  if (fraction != 1)
817  {
818  this->quadrature_points[i + n] =
819  quad2.point(i) * (1 - fraction) + Point<1>(fraction);
820  this->weights[i + n] = quad2.weight(i) * (1 - fraction);
821 
822  // We need to scale with -log|fraction*alpha|
823  this->quadrature_points[j + n] =
824  quad.point(i) * (1 - fraction) + Point<1>(fraction);
825  this->weights[j + n] =
826  -std::log(alpha / (1 - fraction)) * quad.weight(i) * (1 - fraction);
827  }
828  }
829  if (factor_out_singularity == true)
830  for (unsigned int i = 0; i < size(); ++i)
831  {
832  Assert(
833  this->quadrature_points[i] != origin,
834  ExcMessage(
835  "The singularity cannot be on a Gauss point of the same order!"));
836  double denominator =
837  std::log(std::abs((this->quadrature_points[i] - origin)[0]) / alpha);
838  Assert(denominator != 0.0,
839  ExcMessage(
840  "The quadrature formula you are using does not allow to "
841  "factor out the singularity, which is zero at one point."));
842  this->weights[i] /= denominator;
843  }
844 }
845 
846 
847 template <>
848 unsigned int
849 QGaussOneOverR<2>::quad_size(const Point<2> &singularity, const unsigned int n)
850 {
851  const double eps = 1e-8;
852  bool on_edge = false;
853  for (unsigned int d = 0; d < 2; ++d)
854  on_edge = on_edge || (std::abs(singularity[d]) < eps ||
855  std::abs(singularity[d] - 1.0) < eps);
856  const bool on_vertex =
857  on_edge &&
858  std::abs((singularity - Point<2>(.5, .5)).norm_square() - .5) < eps;
859  if (on_vertex)
860  return 2 * n * n;
861  else if (on_edge)
862  return 4 * n * n;
863  else
864  return 8 * n * n;
865 }
866 
867 template <>
869  const Point<2> &singularity,
870  const bool factor_out_singularity)
871  : Quadrature<2>(quad_size(singularity, n))
872 {
873  // We treat all the cases in the
874  // same way. Split the element in 4
875  // pieces, measure the area, if
876  // it's relevant, add the
877  // quadrature connected to that
878  // singularity.
879  std::vector<QGaussOneOverR<2>> quads;
880  std::vector<Point<2>> origins;
881  // Id of the corner with a
882  // singularity
883  quads.emplace_back(n, 3, factor_out_singularity);
884  quads.emplace_back(n, 2, factor_out_singularity);
885  quads.emplace_back(n, 1, factor_out_singularity);
886  quads.emplace_back(n, 0, factor_out_singularity);
887 
888  origins.emplace_back(0., 0.);
889  origins.emplace_back(singularity[0], 0.);
890  origins.emplace_back(0., singularity[1]);
891  origins.push_back(singularity);
892 
893  // Lexicographical ordering.
894 
895  double eps = 1e-8;
896  unsigned int q_id = 0; // Current quad point index.
897  Tensor<1, 2> dist;
898 
899  for (unsigned int box = 0; box < 4; ++box)
900  {
901  dist = (singularity - GeometryInfo<2>::unit_cell_vertex(box));
902  dist = Point<2>(std::abs(dist[0]), std::abs(dist[1]));
903  double area = dist[0] * dist[1];
904  if (area > eps)
905  for (unsigned int q = 0; q < quads[box].size(); ++q, ++q_id)
906  {
907  const Point<2> &qp = quads[box].point(q);
908  this->quadrature_points[q_id] =
909  origins[box] + Point<2>(dist[0] * qp[0], dist[1] * qp[1]);
910  this->weights[q_id] = quads[box].weight(q) * area;
911  }
912  }
913 }
914 
915 
916 template <>
918  const unsigned int vertex_index,
919  const bool factor_out_singularity)
920  : Quadrature<2>(2 * n * n)
921 {
922  // This version of the constructor
923  // works only for the 4
924  // vertices. If you need a more
925  // general one, you should use the
926  // one with the Point<2> in the
927  // constructor.
928  AssertIndexRange(vertex_index, 4);
929 
930  // Start with the gauss quadrature formula on the (u,v) reference
931  // element.
932  QGauss<2> gauss(n);
933 
934  Assert(gauss.size() == n * n, ExcInternalError());
935 
936  // For the moment we only implemented this for the vertices of a
937  // quadrilateral. We are planning to do this also for the support
938  // points of arbitrary FE_Q elements, to allow the use of this
939  // class in boundary element programs with higher order mappings.
940  AssertIndexRange(vertex_index, 4);
941 
942  // We create only the first one. All other pieces are rotation of
943  // this one.
944  // In this case the transformation is
945  //
946  // (x,y) = (u, u tan(pi/4 v))
947  //
948  // with Jacobian
949  //
950  // J = pi/4 R / cos(pi/4 v)
951  //
952  // And we get rid of R to take into account the singularity,
953  // unless specified differently in the constructor.
954  std::vector<Point<2>> &ps = this->quadrature_points;
955  std::vector<double> &ws = this->weights;
956  double pi4 = numbers::PI / 4;
957 
958  for (unsigned int q = 0; q < gauss.size(); ++q)
959  {
960  const Point<2> &gp = gauss.point(q);
961  ps[q][0] = gp[0];
962  ps[q][1] = gp[0] * std::tan(pi4 * gp[1]);
963  ws[q] = gauss.weight(q) * pi4 / std::cos(pi4 * gp[1]);
964  if (factor_out_singularity)
965  ws[q] *= (ps[q] - GeometryInfo<2>::unit_cell_vertex(0)).norm();
966  // The other half of the quadrilateral is symmetric with
967  // respect to xy plane.
968  ws[gauss.size() + q] = ws[q];
969  ps[gauss.size() + q][0] = ps[q][1];
970  ps[gauss.size() + q][1] = ps[q][0];
971  }
972 
973  // Now we distribute these vertices in the correct manner
974  double theta = 0;
975  switch (vertex_index)
976  {
977  case 0:
978  theta = 0;
979  break;
980  case 1:
981  //
982  theta = numbers::PI / 2;
983  break;
984  case 2:
985  theta = -numbers::PI / 2;
986  break;
987  case 3:
988  theta = numbers::PI;
989  break;
990  }
991 
992  double R00 = std::cos(theta), R01 = -std::sin(theta);
993  double R10 = std::sin(theta), R11 = std::cos(theta);
994 
995  if (vertex_index != 0)
996  for (unsigned int q = 0; q < size(); ++q)
997  {
998  double x = ps[q][0] - .5, y = ps[q][1] - .5;
999 
1000  ps[q][0] = R00 * x + R01 * y + .5;
1001  ps[q][1] = R10 * x + R11 * y + .5;
1002  }
1003 }
1004 
1005 
1006 template <int dim>
1008  : Quadrature<dim>(quad)
1009 {
1010  std::vector<unsigned int> permutation(quad.size());
1011  for (unsigned int i = 0; i < quad.size(); ++i)
1012  permutation[i] = i;
1013 
1014  std::sort(permutation.begin(),
1015  permutation.end(),
1016  [this](const unsigned int x, const unsigned int y) {
1017  return this->compare_weights(x, y);
1018  });
1019 
1020  // At this point, the variable is_tensor_product_flag is set
1021  // to the respective value of the given Quadrature in the base
1022  // class copy constructor.
1023  // We only call a quadrature formula 'tensor product'
1024  // if the quadrature points are also sorted lexicographically.
1025  // In particular, any reordering destroys that property
1026  // and we might need to modify the variable accordingly.
1027  for (unsigned int i = 0; i < quad.size(); ++i)
1028  {
1029  this->weights[i] = quad.weight(permutation[i]);
1030  this->quadrature_points[i] = quad.point(permutation[i]);
1031  if (permutation[i] != i)
1032  this->is_tensor_product_flag = false;
1033  }
1034 }
1035 
1036 
1037 template <int dim>
1038 bool
1039 QSorted<dim>::compare_weights(const unsigned int a, const unsigned int b) const
1040 {
1041  return (this->weights[a] < this->weights[b]);
1042 }
1043 
1044 
1045 // construct the quadrature formulae in higher dimensions by
1046 // tensor product of lower dimensions
1047 
1048 template <int dim>
1049 QGauss<dim>::QGauss(const unsigned int n)
1050  : Quadrature<dim>(QGauss<dim - 1>(n), QGauss<1>(n))
1051 {}
1052 
1053 
1054 
1055 template <int dim>
1056 QGaussRadau<dim>::QGaussRadau(const unsigned int n, EndPoint end_point)
1057  : Quadrature<dim>(
1058  QGaussRadau<1>(n, static_cast<QGaussRadau<1>::EndPoint>(end_point)))
1059  , end_point(end_point)
1060 {}
1061 
1062 
1063 
1064 template <int dim>
1066  : Quadrature<dim>(QGaussLobatto<dim - 1>(n), QGaussLobatto<1>(n))
1067 {}
1068 
1069 
1070 
1071 template <int dim>
1073  : Quadrature<dim>(QMidpoint<dim - 1>(), QMidpoint<1>())
1074 {}
1075 
1076 
1077 
1078 template <int dim>
1080  : Quadrature<dim>(QTrapezoid<dim - 1>(), QTrapezoid<1>())
1081 {}
1082 
1083 
1084 
1085 template <int dim>
1087  : Quadrature<dim>(QSimpson<dim - 1>(), QSimpson<1>())
1088 {}
1089 
1090 
1091 
1092 template <int dim>
1094  : Quadrature<dim>(QMilne<dim - 1>(), QMilne<1>())
1095 {}
1096 
1097 
1098 template <int dim>
1100  : Quadrature<dim>(QWeddle<dim - 1>(), QWeddle<1>())
1101 {}
1102 
1103 template <int dim>
1105  const Point<dim> &singularity)
1106  : // We need the explicit implementation if dim == 1. If dim > 1 we use the
1107  // former implementation and apply a tensorial product to obtain the higher
1108  // dimensions.
1109  Quadrature<dim>(
1110  dim == 2 ?
1111  QAnisotropic<dim>(QTelles<1>(base_quad, Point<1>(singularity[0])),
1112  QTelles<1>(base_quad, Point<1>(singularity[1]))) :
1113  dim == 3 ?
1114  QAnisotropic<dim>(QTelles<1>(base_quad, Point<1>(singularity[0])),
1115  QTelles<1>(base_quad, Point<1>(singularity[1])),
1116  QTelles<1>(base_quad, Point<1>(singularity[2]))) :
1117  Quadrature<dim>())
1118 {}
1119 
1120 template <int dim>
1121 QTelles<dim>::QTelles(const unsigned int n, const Point<dim> &singularity)
1122  : // In this case we map the standard Gauss Legendre formula using the given
1123  // singularity point coordinates.
1124  Quadrature<dim>(QTelles<dim>(QGauss<1>(n), singularity))
1125 {}
1126 
1127 
1128 
1129 template <>
1130 QTelles<1>::QTelles(const Quadrature<1> &base_quad, const Point<1> &singularity)
1131  : // We explicitly implement the Telles' variable change if dim == 1.
1132  Quadrature<1>(base_quad)
1133 {
1134  // We define all the constants to be used in the implementation of
1135  // Telles' rule
1136  const double eta_bar = singularity[0] * 2. - 1.;
1137  const double eta_star = eta_bar * eta_bar - 1.;
1138  double gamma_bar;
1139 
1140  std::vector<Point<1>> quadrature_points_dummy(quadrature_points.size());
1141  std::vector<double> weights_dummy(weights.size());
1142  unsigned int cont = 0;
1143  const double tol = 1e-10;
1144  for (unsigned int d = 0; d < quadrature_points.size(); ++d)
1145  {
1146  if (std::abs(quadrature_points[d][0] - singularity[0]) > tol)
1147  {
1148  quadrature_points_dummy[d - cont] = quadrature_points[d];
1149  weights_dummy[d - cont] = weights[d];
1150  }
1151  else
1152  {
1153  // We need to remove the singularity point from the quadrature point
1154  // list. To do so we use the variable cont.
1155  cont = 1;
1156  }
1157  }
1158  if (cont == 1)
1159  {
1160  quadrature_points.resize(quadrature_points_dummy.size() - 1);
1161  weights.resize(weights_dummy.size() - 1);
1162  for (unsigned int d = 0; d < quadrature_points.size(); ++d)
1163  {
1164  quadrature_points[d] = quadrature_points_dummy[d];
1165  weights[d] = weights_dummy[d];
1166  }
1167  }
1168  // We need to check if the singularity is at the boundary of the interval.
1169  if (std::abs(eta_star) <= tol)
1170  {
1171  gamma_bar =
1172  std::pow((eta_bar * eta_star + std::abs(eta_star)), 1.0 / 3.0) +
1173  std::pow((eta_bar * eta_star - std::abs(eta_star)), 1.0 / 3.0) +
1174  eta_bar;
1175  }
1176  else
1177  {
1178  gamma_bar = (eta_bar * eta_star + std::abs(eta_star)) /
1179  std::abs(eta_bar * eta_star + std::abs(eta_star)) *
1180  std::pow(std::abs(eta_bar * eta_star + std::abs(eta_star)),
1181  1.0 / 3.0) +
1182  (eta_bar * eta_star - std::abs(eta_star)) /
1183  std::abs(eta_bar * eta_star - std::abs(eta_star)) *
1184  std::pow(std::abs(eta_bar * eta_star - std::abs(eta_star)),
1185  1.0 / 3.0) +
1186  eta_bar;
1187  }
1188  for (unsigned int q = 0; q < quadrature_points.size(); ++q)
1189  {
1190  double gamma = quadrature_points[q][0] * 2 - 1;
1191  double eta = (Utilities::fixed_power<3>(gamma - gamma_bar) +
1192  gamma_bar * (gamma_bar * gamma_bar + 3)) /
1193  (1 + 3 * gamma_bar * gamma_bar);
1194 
1195  double J = 3 * ((gamma - gamma_bar) * (gamma - gamma_bar)) /
1196  (1 + 3 * gamma_bar * gamma_bar);
1197 
1198  quadrature_points[q][0] = (eta + 1) / 2.0;
1199  weights[q] = J * weights[q];
1200  }
1201 }
1202 
1203 namespace internal
1204 {
1206  {
1210  std::vector<double>
1211  get_quadrature_points(const unsigned int n)
1212  {
1213  std::vector<double> points(n);
1214  // n point quadrature: index from 0 to n-1
1215  for (unsigned short i = 0; i < n; ++i)
1216  // would be cos((2i+1)Pi/(2N+2))
1217  // put + Pi so we start from the smallest point
1218  // then map from [-1,1] to [0,1]
1219  points[i] =
1220  1. / 2. *
1221  (1. + std::cos(numbers::PI *
1222  (1. + double(2 * i + 1) / double(2 * (n - 1) + 2))));
1223 
1224  return points;
1225  }
1226 
1227 
1228 
1232  std::vector<double>
1233  get_quadrature_weights(const unsigned int n)
1234  {
1235  std::vector<double> weights(n);
1236 
1237  for (unsigned short i = 0; i < n; ++i)
1238  {
1239  // same weights as on [-1,1]
1240  weights[i] = numbers::PI / double(n);
1241  }
1242 
1243  return weights;
1244  }
1245  } // namespace QGaussChebyshev
1246 } // namespace internal
1247 
1248 
1249 template <>
1251  : Quadrature<1>(n)
1252 {
1253  Assert(n > 0, ExcMessage("Need at least one point for the quadrature rule"));
1254  std::vector<double> p = internal::QGaussChebyshev::get_quadrature_points(n);
1255  std::vector<double> w = internal::QGaussChebyshev::get_quadrature_weights(n);
1256 
1257  for (unsigned int i = 0; i < this->size(); ++i)
1258  {
1259  this->quadrature_points[i] = Point<1>(p[i]);
1260  this->weights[i] = w[i];
1261  }
1262 }
1263 
1264 
1265 template <int dim>
1267  : Quadrature<dim>(QGaussChebyshev<1>(n))
1268 {}
1269 
1270 
1271 namespace internal
1272 {
1274  {
1275  // Computes the points of the quadrature formula.
1276  std::vector<double>
1278  const unsigned int n,
1279  const ::QGaussRadauChebyshev<1>::EndPoint end_point)
1280  {
1281  std::vector<double> points(n);
1282  // n point quadrature: index from 0 to n-1
1283  for (unsigned short i = 0; i < n; ++i)
1284  // would be -cos(2i Pi/(2N+1))
1285  // put + Pi so we start from the smallest point
1286  // then map from [-1,1] to [0,1]
1287  switch (end_point)
1288  {
1289  case ::QGaussRadauChebyshev<1>::EndPoint::left:
1290  {
1291  points[i] =
1292  1. / 2. *
1293  (1. -
1294  std::cos(numbers::PI *
1295  (1 + 2 * double(i) / (2 * double(n - 1) + 1.))));
1296  break;
1297  }
1298 
1299  case ::QGaussRadauChebyshev<1>::EndPoint::right:
1300  {
1301  points[i] =
1302  1. / 2. *
1303  (1. - std::cos(numbers::PI * (2 * double(n - 1 - i) /
1304  (2 * double(n - 1) + 1.))));
1305  break;
1306  }
1307 
1308  default:
1309  Assert(
1310  false,
1311  ExcMessage(
1312  "This constructor can only be called with either "
1313  "QGaussRadauChebyshev::EndPoint::left or "
1314  "QGaussRadauChebyshev::EndPoint:right as second argument."));
1315  }
1316 
1317  return points;
1318  }
1319 
1320 
1321 
1322  // Computes the weights of the quadrature formula.
1323  std::vector<double>
1325  const unsigned int n,
1326  const ::QGaussRadauChebyshev<1>::EndPoint end_point)
1327  {
1328  std::vector<double> weights(n);
1329 
1330  for (unsigned short i = 0; i < n; ++i)
1331  {
1332  // same weights as on [-1,1]
1333  weights[i] = 2. * numbers::PI / double(2 * (n - 1) + 1.);
1334  if (end_point == ::QGaussRadauChebyshev<1>::EndPoint::left &&
1335  i == 0)
1336  weights[i] /= 2.;
1337  else if (end_point ==
1339  i == (n - 1))
1340  weights[i] /= 2.;
1341  }
1342 
1343  return weights;
1344  }
1345  } // namespace QGaussRadauChebyshev
1346 } // namespace internal
1347 
1348 
1349 template <>
1351  const EndPoint end_point)
1352  : Quadrature<1>(n)
1353  , end_point(end_point)
1354 {
1355  Assert(n > 0, ExcMessage("Need at least one point for quadrature rules."));
1356  std::vector<double> points =
1358  std::vector<double> weights =
1360 
1361  for (unsigned int i = 0; i < this->size(); ++i)
1362  {
1363  this->quadrature_points[i] = Point<1>(points[i]);
1364  this->weights[i] = weights[i];
1365  }
1366 }
1367 
1368 
1369 template <int dim>
1371  EndPoint end_point)
1372  : Quadrature<dim>(QGaussRadauChebyshev<1>(
1373  n,
1374  static_cast<QGaussRadauChebyshev<1>::EndPoint>(end_point)))
1375  , end_point(end_point)
1376 {}
1377 
1378 
1379 
1380 namespace internal
1381 {
1383  {
1384  // Computes the points of the quadrature formula.
1385  std::vector<double>
1386  get_quadrature_points(const unsigned int n)
1387  {
1388  std::vector<double> points(n);
1389  // n point quadrature: index from 0 to n-1
1390  for (unsigned short i = 0; i < n; ++i)
1391  // would be cos(i Pi/N)
1392  // put + Pi so we start from the smallest point
1393  // then map from [-1,1] to [0,1]
1394  points[i] =
1395  1. / 2. *
1396  (1. + std::cos(numbers::PI * (1 + double(i) / double(n - 1))));
1397 
1398  return points;
1399  }
1400 
1401  // Computes the weights of the quadrature formula.
1402  std::vector<double>
1403  get_quadrature_weights(const unsigned int n)
1404  {
1405  std::vector<double> weights(n);
1406 
1407  for (unsigned short i = 0; i < n; ++i)
1408  {
1409  // same weights as on [-1,1]
1410  weights[i] = numbers::PI / double((n - 1));
1411  if (i == 0 || i == (n - 1))
1412  weights[i] /= 2.;
1413  }
1414 
1415  return weights;
1416  }
1417  } // namespace QGaussLobattoChebyshev
1418 } // namespace internal
1419 
1420 
1421 
1422 template <>
1424  : Quadrature<1>(n)
1425 {
1426  Assert(n > 1,
1427  ExcMessage(
1428  "Need at least two points for Gauss-Lobatto quadrature rule"));
1429  std::vector<double> p =
1431  std::vector<double> w =
1433 
1434  for (unsigned int i = 0; i < this->size(); ++i)
1435  {
1436  this->quadrature_points[i] = Point<1>(p[i]);
1437  this->weights[i] = w[i];
1438  }
1439 }
1440 
1441 
1442 template <int dim>
1444  : Quadrature<dim>(QGaussLobattoChebyshev<1>(n))
1445 {}
1446 
1447 
1448 
1449 template <int dim>
1451 {
1452  std::vector<Point<dim>> qpoints;
1453  std::vector<double> weights;
1454 
1455  for (unsigned int i = 0; i < quad.size(); ++i)
1456  {
1457  double r = 0;
1458  /* Use "int d" instead of the more natural "unsigned int d" to work
1459  * around a wrong diagnostic in gcc-10.3.0 that warns about that the
1460  * comparison "d < dim" is always false in case of "dim == 0".
1461  * MM 2021 */
1462  for (int d = 0; d < dim; ++d)
1463  r += quad.point(i)[d];
1464  if (r <= 1 + 1e-10)
1465  {
1466  this->quadrature_points.push_back(quad.point(i));
1467  this->weights.push_back(quad.weight(i));
1468  }
1469  }
1470 }
1471 
1472 
1473 
1474 template <int dim>
1475 template <int spacedim>
1478  const std::array<Point<spacedim>, dim + 1> &vertices) const
1479 {
1480  Assert(dim <= spacedim,
1481  ExcMessage("Invalid combination of dim and spacedim ."));
1483  for (unsigned int d = 0; d < dim; ++d)
1484  Bt[d] = vertices[d + 1] - vertices[0];
1485 
1486  const auto B = Bt.transpose();
1487  const double J = std::abs(B.determinant());
1488 
1489  // if the determinant is zero, we return an empty quadrature
1490  if (J < 1e-12)
1491  return Quadrature<spacedim>();
1492 
1493  std::vector<Point<spacedim>> qp(this->size());
1494  std::vector<double> w(this->size());
1495 
1496  for (unsigned int i = 0; i < this->size(); ++i)
1497  {
1498  qp[i] =
1500  w[i] = J * this->weight(i);
1501  }
1502 
1503  return Quadrature<spacedim>(qp, w);
1504 }
1505 
1506 
1507 
1508 template <int dim>
1509 template <int spacedim>
1512  const std::vector<std::array<Point<spacedim>, dim + 1>> &simplices) const
1513 {
1514  Assert(!(dim == 1 && spacedim == 1),
1515  ExcMessage("This function is not supposed to work in 1D-1d case."));
1516  Assert(dim <= spacedim,
1517  ExcMessage("Invalid combination of dim and spacedim ."));
1518 
1519  std::vector<Point<spacedim>> qp;
1520  std::vector<double> ws;
1521  for (const auto &simplex : simplices)
1522  {
1523  const auto rule = this->compute_affine_transformation(simplex);
1524  std::transform(rule.get_points().begin(),
1525  rule.get_points().end(),
1526  std::back_inserter(qp),
1527  [&](const Point<spacedim> &p) { return p; });
1528  std::transform(rule.get_weights().begin(),
1529  rule.get_weights().end(),
1530  std::back_inserter(ws),
1531  [&](const double w) { return w; });
1532  }
1533  return Quadrature<spacedim>(qp, ws);
1534 }
1535 
1536 
1537 
1539  const Quadrature<1> &angular_quadrature)
1540  : QSimplex<2>(Quadrature<2>())
1541 {
1542  const QAnisotropic<2> base(radial_quadrature, angular_quadrature);
1543  this->quadrature_points.resize(base.size());
1544  this->weights.resize(base.size());
1545  for (unsigned int i = 0; i < base.size(); ++i)
1546  {
1547  const auto &q = base.point(i);
1548  const auto w = base.weight(i);
1549 
1550  const auto xhat = q[0];
1551  const auto yhat = q[1];
1552 
1553  const double t = numbers::PI_2 * yhat;
1554  const double pi = numbers::PI;
1555  const double st = std::sin(t);
1556  const double ct = std::cos(t);
1557  const double r = xhat / (st + ct);
1558 
1559  const double J = pi * xhat / (2 * (std::sin(pi * yhat) + 1));
1560 
1561  this->quadrature_points[i] = Point<2>(r * ct, r * st);
1562  this->weights[i] = w * J;
1563  }
1564 }
1565 
1566 
1567 
1568 QTrianglePolar::QTrianglePolar(const unsigned int n)
1569  : QTrianglePolar(QGauss<1>(n), QGauss<1>(n))
1570 {}
1571 
1572 
1573 
1574 QDuffy::QDuffy(const Quadrature<1> &radial_quadrature,
1575  const Quadrature<1> &angular_quadrature,
1576  const double beta)
1577  : QSimplex<2>(Quadrature<2>())
1578 {
1579  const QAnisotropic<2> base(radial_quadrature, angular_quadrature);
1580  this->quadrature_points.resize(base.size());
1581  this->weights.resize(base.size());
1582  for (unsigned int i = 0; i < base.size(); ++i)
1583  {
1584  const auto &q = base.point(i);
1585  const auto w = base.weight(i);
1586 
1587  const auto xhat = q[0];
1588  const auto yhat = q[1];
1589 
1590  const double x = std::pow(xhat, beta) * (1 - yhat);
1591  const double y = std::pow(xhat, beta) * yhat;
1592 
1593  const double J = beta * std::pow(xhat, 2. * beta - 1.);
1594 
1595  this->quadrature_points[i] = Point<2>(x, y);
1596  this->weights[i] = w * J;
1597  }
1598 }
1599 
1600 
1601 
1602 QDuffy::QDuffy(const unsigned int n, const double beta)
1603  : QDuffy(QGauss<1>(n), QGauss<1>(n), beta)
1604 {}
1605 
1606 
1607 
1608 template <int dim>
1610 {
1612  ExcMessage(
1613  "The split point should be inside the unit reference cell."));
1614 
1615  std::array<Point<dim>, dim + 1> vertices;
1616  vertices[0] = split_point;
1617 
1618  // Make a simplex from the split_point and the first dim vertices of each
1619  // face. In dimension three, we need to split the face in two triangles, so
1620  // we use once the first dim vertices of each face, and the second time the
1621  // the dim vertices of each face starting from 1.
1622  for (auto f : GeometryInfo<dim>::face_indices())
1623  for (unsigned int start = 0; start < (dim > 2 ? 2 : 1); ++start)
1624  {
1625  for (unsigned int i = 0; i < dim; ++i)
1628  const auto quad = base.compute_affine_transformation(vertices);
1629  if (quad.size())
1630  {
1631  this->quadrature_points.insert(this->quadrature_points.end(),
1632  quad.get_points().begin(),
1633  quad.get_points().end());
1634  this->weights.insert(this->weights.end(),
1635  quad.get_weights().begin(),
1636  quad.get_weights().end());
1637  }
1638  }
1639 }
1640 
1641 
1642 
1643 template <int dim>
1644 QGaussSimplex<dim>::QGaussSimplex(const unsigned int n_points_1D)
1645  : QSimplex<dim>(Quadrature<dim>())
1646 {
1647  // fill quadrature points and quadrature weights
1648  if (dim == 0 || dim == 1)
1649  {
1650  const ::QGauss<dim> quad(n_points_1D);
1651 
1652  this->quadrature_points = quad.get_points();
1653  this->weights = quad.get_weights();
1654  }
1655  else if (dim == 2)
1656  {
1657  if (n_points_1D == 1)
1658  {
1659  const double p = 1.0 / 3.0;
1660  this->quadrature_points.emplace_back(p, p);
1661  this->weights.emplace_back(0.5);
1662  }
1663  else if (n_points_1D == 2)
1664  {
1665  // The Hillion 7 scheme, as communicated by quadpy
1666  //
1667  // See: Numerical Integration on a Triangle, International Journal for
1668  // Numerical Methods in Engineering, 1977
1669  const double Q12 = 1.0 / 2.0;
1670  this->quadrature_points.emplace_back(0.17855872826361643,
1671  0.1550510257216822);
1672  this->quadrature_points.emplace_back(0.07503111022260812,
1673  0.6449489742783178);
1674  this->quadrature_points.emplace_back(0.6663902460147014,
1675  0.1550510257216822);
1676  this->quadrature_points.emplace_back(0.28001991549907407,
1677  0.6449489742783178);
1678 
1679  this->weights.emplace_back(0.31804138174397717 * Q12);
1680  this->weights.emplace_back(0.18195861825602283 * Q12);
1681  this->weights.emplace_back(0.31804138174397717 * Q12);
1682  this->weights.emplace_back(0.18195861825602283 * Q12);
1683  }
1684  else if (n_points_1D == 3)
1685  {
1686  // The Hammer-Marlowe-Stroud 5 Scheme, as communicated by quadpy
1687  const double p0 = 2.0 / 7.0 - std::sqrt(15.0) / 21.0;
1688  const double p1 = 2.0 / 7.0 + std::sqrt(15.0) / 21.0;
1689  const double p2 = 3.0 / 7.0 - 2.0 * std::sqrt(15.0) / 21.0;
1690  const double p3 = 3.0 / 7.0 + 2.0 * std::sqrt(15.0) / 21.0;
1691  this->quadrature_points.emplace_back(1.0 / 3.0, 1.0 / 3.0);
1692  this->quadrature_points.emplace_back(p3, p0);
1693  this->quadrature_points.emplace_back(p0, p3);
1694  this->quadrature_points.emplace_back(p0, p0);
1695  this->quadrature_points.emplace_back(p2, p1);
1696  this->quadrature_points.emplace_back(p1, p2);
1697  this->quadrature_points.emplace_back(p1, p1);
1698 
1699  const double q12 = 0.5;
1700  const double w0 = 9.0 / 40.0;
1701  const double w1 = 31.0 / 240.0 - std::sqrt(15.0) / 1200.0;
1702  const double w2 = 31.0 / 240.0 + std::sqrt(15.0) / 1200.0;
1703  this->weights.emplace_back(q12 * w0);
1704  this->weights.emplace_back(q12 * w1);
1705  this->weights.emplace_back(q12 * w1);
1706  this->weights.emplace_back(q12 * w1);
1707  this->weights.emplace_back(q12 * w2);
1708  this->weights.emplace_back(q12 * w2);
1709  this->weights.emplace_back(q12 * w2);
1710  }
1711  else if (n_points_1D == 4)
1712  {
1714  QWitherdenVincentSimplex<dim>(n_points_1D));
1715  }
1716  }
1717  else if (dim == 3)
1718  {
1719  if (n_points_1D == 1)
1720  {
1721  const double Q14 = 1.0 / 4.0;
1722  const double Q16 = 1.0 / 6.0;
1723 
1724  this->quadrature_points.emplace_back(Q14, Q14, Q14);
1725  this->weights.emplace_back(Q16);
1726  }
1727  // The Xiao Gimbutas 03 scheme, as communicated by quadpy
1728  //
1729  // See: A numerical algorithm for the construction of efficient quadrature
1730  // rules in two and higher dimensions, Computers & Mathematics with
1731  // Applications, 2010
1732  else if (n_points_1D == 2)
1733  {
1734  const double Q16 = 1.0 / 6.0;
1735  this->weights.emplace_back(0.1223220027573451 * Q16);
1736  this->weights.emplace_back(0.1280664127107469 * Q16);
1737  this->weights.emplace_back(0.1325680271444452 * Q16);
1738  this->weights.emplace_back(0.1406244096604032 * Q16);
1739  this->weights.emplace_back(0.2244151669175574 * Q16);
1740  this->weights.emplace_back(0.2520039808095023 * Q16);
1741 
1742  this->quadrature_points.emplace_back(0.1620014916985245,
1743  0.1838503504920977,
1744  0.01271836631368145);
1745  this->quadrature_points.emplace_back(0.01090521221118924,
1746  0.2815238021235462,
1747  0.3621268299455338);
1748  this->quadrature_points.emplace_back(0.1901170024392839,
1749  0.01140332944455717,
1750  0.3586207204668839);
1751  this->quadrature_points.emplace_back(0.170816925164989,
1752  0.1528181430909273,
1753  0.6384932999617267);
1754  this->quadrature_points.emplace_back(0.1586851632274406,
1755  0.5856628056552158,
1756  0.1308471689520965);
1757  this->quadrature_points.emplace_back(0.5712260521491151,
1758  0.1469183900871696,
1759  0.1403728057942107);
1760  }
1761  // Past this point the best rules (positive weights, minimal number of
1762  // points) we have right now are the Witherden-Vincent ones
1763  else if (n_points_1D == 3)
1764  {
1766  QWitherdenVincentSimplex<dim>(n_points_1D));
1767  }
1768  else if (n_points_1D == 4)
1769  {
1771  QWitherdenVincentSimplex<dim>(n_points_1D));
1772  }
1773  }
1774 
1775  AssertDimension(this->quadrature_points.size(), this->weights.size());
1776  Assert(this->quadrature_points.size() > 0,
1778  "QGaussSimplex is currently only implemented for "
1779  "n_points_1D = 1, 2, 3, and 4 while you are asking for "
1780  "n_points_1D = " +
1781  Utilities::to_string(n_points_1D)));
1782 }
1783 
1784 namespace
1785 {
1786  template <std::size_t b_dim>
1787  std::vector<std::array<double, b_dim>>
1788  all_permutations(const std::array<double, b_dim> &b_point)
1789  {
1790  std::vector<std::array<double, b_dim>> output;
1791 
1792  // We want all possible permutations of the barycentric coordinates.
1793  // The easiest way to get all of them is to sort the input first and
1794  // then use next_permutation to cycle through them all.
1795  std::array<double, b_dim> temp = b_point;
1796  std::sort(temp.begin(), temp.end());
1797  do
1798  {
1799  output.push_back(temp);
1800  }
1801  while (std::next_permutation(temp.begin(), temp.end()));
1802 
1803  return output;
1804  }
1805 } // namespace
1806 
1807 
1808 
1809 template <int dim>
1811  const unsigned int n_points_1D,
1812  const bool use_odd_order)
1813  : QSimplex<dim>(Quadrature<dim>())
1814 {
1815  Assert(1 <= dim && dim <= 3, ExcNotImplemented());
1816  // Just use Gauss in 1d: this is a high-order open rule so this is a
1817  // reasonable equivalent for generic programming.
1818  if (dim == 1)
1819  {
1821  return;
1822  }
1823 
1824  std::array<double, dim + 1> centroid;
1825  std::fill(centroid.begin(), centroid.end(), 1.0 / (dim + 1.0));
1826  std::vector<std::vector<std::array<double, dim + 1>>> b_point_permutations;
1827  std::vector<double> b_weights;
1828 
1829  // We can simplify the implementation of these quadrature rules
1830  // by quite a bit by exploiting symmetry - we do essentially the
1831  // same thing for each barycentric coordinate, so we can express
1832  // our quadrature rule as permutations of barycentric points
1833  // instead of writing things out explicitly.
1834 
1835  // Apply a Barycentric permutation where one point is different.
1836  // Equivalent to d3_aa and s31 in quadpy.
1837  auto process_point_1 = [&](const double a, const double w) {
1838  const double b = 1.0 - dim * a;
1839  std::array<double, dim + 1> b_point;
1840  std::fill(b_point.begin(), b_point.begin() + dim, a);
1841  b_point[dim] = b;
1842 
1843  b_weights.push_back(w);
1844  b_point_permutations.push_back(all_permutations(b_point));
1845  };
1846 
1847  // Apply a Barycentric permutation where two points (in 3d) are different.
1848  // Equivalent to s22 in quadpy.
1849  auto process_point_2 = [&](const double a, const double w) {
1850  Assert(dim == 3, ExcInternalError());
1851  const double b = (1.0 - 2.0 * a) / 2.0;
1852  std::array<double, dim + 1> b_point;
1853  std::fill(b_point.begin(), b_point.begin() + dim - 1, a);
1854  b_point[dim - 1] = b;
1855  b_point[dim] = b;
1856 
1857  b_weights.push_back(w);
1858  b_point_permutations.push_back(all_permutations(b_point));
1859  };
1860 
1861  // Apply a Barycentric permutation where three (or four) points
1862  // are different (since there are two inputs).
1863  // Equivalent to d3_ab and s211 in quadpy.
1864  auto process_point_3 = [&](const double a, const double b, const double w) {
1865  const double c = 1.0 - (dim - 1.0) * a - b;
1866  std::array<double, dim + 1> b_point;
1867  std::fill(b_point.begin(), b_point.begin() + dim - 1, a);
1868  b_point[dim - 1] = b;
1869  b_point[dim] = c;
1870 
1871  b_weights.push_back(w);
1872  b_point_permutations.push_back(all_permutations(b_point));
1873  };
1874 
1875  switch (n_points_1D)
1876  {
1877  case 1:
1878  switch (dim)
1879  {
1880  case 2:
1881  if (use_odd_order)
1882  {
1883  // WV-1, 2d
1884  b_point_permutations.push_back({centroid});
1885  b_weights.push_back(1.0000000000000000e+00);
1886  }
1887  else
1888  {
1889  // WV-2, 2d
1890  process_point_1(1.6666666666666669e-01,
1891  3.3333333333333331e-01);
1892  }
1893  break;
1894  case 3:
1895  if (use_odd_order)
1896  {
1897  // WV-1, 3d
1898  b_point_permutations.push_back({centroid});
1899  b_weights.push_back(1.0000000000000000e+00);
1900  }
1901  else
1902  {
1903  // WV-2, 3d
1904  process_point_1(1.3819660112501050e-01,
1905  2.5000000000000000e-01);
1906  }
1907  break;
1908  default:
1909  Assert(false, ExcNotImplemented());
1910  }
1911  break;
1912  case 2:
1913  switch (dim)
1914  {
1915  case 2:
1916  // WV-4 in both cases (no WV-3 in 2d)
1917  process_point_1(9.1576213509770743e-02, 1.0995174365532187e-01);
1918  process_point_1(4.4594849091596489e-01, 2.2338158967801147e-01);
1919  break;
1920  case 3:
1921  if (use_odd_order)
1922  {
1923  // WV-3, 3d
1924  process_point_1(3.2816330251638171e-01,
1925  1.3621784253708741e-01);
1926  process_point_1(1.0804724989842859e-01,
1927  1.1378215746291261e-01);
1928  }
1929  else
1930  {
1931  // WV-5 (no WV-4 in 3d)
1933  }
1934  break;
1935  default:
1936  Assert(false, ExcInternalError());
1937  }
1938  break;
1939  case 3:
1940  switch (dim)
1941  {
1942  case 2:
1943  if (use_odd_order)
1944  {
1945  // WV-5, 2d
1946  b_point_permutations.push_back({centroid});
1947  b_weights.push_back(2.2500000000000001e-01);
1948  process_point_1(1.0128650732345634e-01,
1949  1.2593918054482714e-01);
1950  process_point_1(4.7014206410511511e-01,
1951  1.3239415278850619e-01);
1952  }
1953  else
1954  {
1955  // WV-6, 2d
1956  process_point_1(6.3089014491502227e-02,
1957  5.0844906370206819e-02);
1958  process_point_1(2.4928674517091043e-01,
1959  1.1678627572637937e-01);
1960  process_point_3(5.3145049844816938e-02,
1961  3.1035245103378439e-01,
1962  8.2851075618373571e-02);
1963  }
1964  break;
1965  case 3:
1966  if (use_odd_order)
1967  {
1968  // WV-5, 3d
1969  process_point_1(3.1088591926330061e-01,
1970  1.1268792571801590e-01);
1971  process_point_1(9.2735250310891248e-02,
1972  7.3493043116361956e-02);
1973  process_point_2(4.5503704125649642e-02,
1974  4.2546020777081472e-02);
1975  }
1976  else
1977  {
1978  // WV-6, 3d
1979  process_point_1(4.0673958534611372e-02,
1980  1.0077211055320640e-02);
1981  process_point_1(3.2233789014227548e-01,
1982  5.5357181543654717e-02);
1983  process_point_1(2.1460287125915201e-01,
1984  3.9922750258167487e-02);
1985  process_point_3(6.3661001875017442e-02,
1986  6.0300566479164919e-01,
1987  4.8214285714285710e-02);
1988  }
1989  break;
1990  default:
1991  Assert(false, ExcInternalError());
1992  }
1993  break;
1994  case 4:
1995  switch (dim)
1996  {
1997  case 2:
1998  if (use_odd_order)
1999  {
2000  // WV-7, 2d
2001  process_point_1(3.3730648554587850e-02,
2002  1.6545050110792131e-02);
2003  process_point_1(4.7430969250471822e-01,
2004  7.7086646185986069e-02);
2005  process_point_1(2.4157738259540357e-01,
2006  1.2794417123015558e-01);
2007  process_point_3(4.7036644652595216e-02,
2008  1.9868331479735168e-01,
2009  5.5878732903199779e-02);
2010  }
2011  else
2012  {
2013  // WV-8, 2d
2014  b_point_permutations.push_back({centroid});
2015  b_weights.push_back(1.4431560767778717e-01);
2016  process_point_1(5.0547228317030957e-02,
2017  3.2458497623198079e-02);
2018  process_point_1(4.5929258829272313e-01,
2019  9.5091634267284619e-02);
2020  process_point_1(1.7056930775176021e-01,
2021  1.0321737053471824e-01);
2022  process_point_3(8.3947774099575878e-03,
2023  2.6311282963463811e-01,
2024  2.7230314174434993e-02);
2025  }
2026  break;
2027  case 3:
2028  if (use_odd_order)
2029  {
2030  // WV-7, 3d
2031  b_point_permutations.push_back({centroid});
2032  b_weights.push_back(9.5485289464130846e-02);
2033  process_point_1(3.1570114977820279e-01,
2034  4.2329581209967028e-02);
2035  process_point_2(5.0489822598396350e-02,
2036  3.1896927832857580e-02);
2037  process_point_3(1.8883383102600099e-01,
2038  5.7517163758699996e-01,
2039  3.7207130728334620e-02);
2040  process_point_3(2.1265472541483140e-02,
2041  8.1083024109854862e-01,
2042  8.1107708299033420e-03);
2043  }
2044  else
2045  {
2046  // WV-8, 3d
2047  process_point_1(1.0795272496221089e-01,
2048  2.6426650908408830e-02);
2049  process_point_1(1.8510948778258660e-01,
2050  5.2031747563738531e-02);
2051  process_point_1(4.2316543684767283e-02,
2052  7.5252561535401989e-03);
2053  process_point_1(3.1418170912403898e-01,
2054  4.1763782856934897e-02);
2055  process_point_2(4.3559132858383021e-01,
2056  3.6280930261308818e-02);
2057  process_point_3(2.1433930127130570e-02,
2058  7.1746406342630831e-01,
2059  7.1569028908444327e-03);
2060  process_point_3(2.0413933387602909e-01,
2061  5.8379737830214440e-01,
2062  1.5453486150960340e-02);
2063  }
2064  break;
2065  default:
2066  Assert(false, ExcInternalError());
2067  }
2068  break;
2069  case 5:
2070  switch (dim)
2071  {
2072  case 2:
2073  if (use_odd_order)
2074  {
2075  // WV-9, 2d
2076  b_point_permutations.push_back({centroid});
2077  b_weights.push_back(9.7135796282798836e-02);
2078  process_point_1(4.4729513394452691e-02,
2079  2.5577675658698031e-02);
2080  process_point_1(4.8968251919873762e-01,
2081  3.1334700227139071e-02);
2082  process_point_1(4.3708959149293664e-01,
2083  7.7827541004774278e-02);
2084  process_point_1(1.8820353561903275e-01,
2085  7.9647738927210249e-02);
2086  process_point_3(3.6838412054736258e-02,
2087  2.2196298916076568e-01,
2088  4.3283539377289376e-02);
2089  }
2090  else
2091  {
2092  // WV-10, 2d
2093  b_point_permutations.push_back({centroid});
2094  b_weights.push_back(8.1743329146285973e-02);
2095  process_point_1(3.2055373216943517e-02,
2096  1.3352968813149567e-02);
2097  process_point_1(1.4216110105656438e-01,
2098  4.5957963604744731e-02);
2099  process_point_3(2.8367665339938453e-02,
2100  1.6370173373718250e-01,
2101  2.5297757707288385e-02);
2102  process_point_3(2.9619889488729734e-02,
2103  3.6914678182781102e-01,
2104  3.4184648162959429e-02);
2105  process_point_3(1.4813288578382056e-01,
2106  3.2181299528883545e-01,
2107  6.3904906396424044e-02);
2108  }
2109  break;
2110  case 3:
2111  if (use_odd_order)
2112  {
2113  // WV-9, 3d
2114  b_point_permutations.push_back({centroid});
2115  b_weights.push_back(5.8010548912480253e-02);
2116  process_point_1(6.1981697552226933e-10,
2117  6.4319281759256394e-05);
2118  process_point_1(1.6077453539526160e-01,
2119  2.3173338462425461e-02);
2120  process_point_1(3.2227652182142102e-01,
2121  2.9562912335429289e-02);
2122  process_point_1(4.5108918345413578e-02,
2123  8.0639799796161822e-03);
2124  process_point_2(1.1229654600437609e-01,
2125  3.8134080103702457e-02);
2126  process_point_3(4.5887144875245922e-01,
2127  2.5545792330413102e-03,
2128  8.3844221982985519e-03);
2129  process_point_3(3.3775870685338598e-02,
2130  7.1835032644207453e-01,
2131  1.0234559352745330e-02);
2132  process_point_3(1.8364136980992790e-01,
2133  3.4415910578175279e-02,
2134  2.0524915967988139e-02);
2135  }
2136  else
2137  {
2138  // WV-10, 3d
2139  b_point_permutations.push_back({centroid});
2140  b_weights.push_back(4.7399773556020743e-02);
2141  process_point_1(3.1225006869518868e-01,
2142  2.6937059992268701e-02);
2143  process_point_1(1.1430965385734609e-01,
2144  9.8691597167933822e-03);
2145  process_point_3(4.1043073921896539e-01,
2146  1.6548602561961109e-01,
2147  1.1393881220195230e-02);
2148  process_point_3(6.1380088247906528e-03,
2149  9.4298876734520487e-01,
2150  3.6194434433925362e-04);
2151  process_point_3(1.2105018114558939e-01,
2152  4.7719037990428043e-01,
2153  2.5739731980456069e-02);
2154  process_point_3(3.2779468216442620e-02,
2155  5.9425626948000698e-01,
2156  1.0135871679755789e-02);
2157  process_point_3(3.2485281564823047e-02,
2158  8.0117728465834437e-01,
2159  6.5761472770359038e-03);
2160  process_point_3(1.7497934218393901e-01,
2161  6.2807184547536599e-01,
2162  1.2907035798861989e-02);
2163  }
2164  break;
2165  default:
2166  Assert(false, ExcNotImplemented());
2167  }
2168  break;
2169  case 6:
2170  // There is no WV-11 rule in 3d yet
2171  Assert(dim == 2, ExcNotImplemented());
2172  if (use_odd_order)
2173  {
2174  // WV-11, 2d
2175  b_point_permutations.push_back({centroid});
2176  b_weights.push_back(8.5761179732224219e-02);
2177  process_point_1(2.8485417614371900e-02, 1.0431870512894697e-02);
2178  process_point_1(4.9589190096589092e-01, 1.6606273054585369e-02);
2179  process_point_1(1.0263548271224643e-01, 3.8630759237019321e-02);
2180  process_point_1(4.3846592676435220e-01, 6.7316154079468296e-02);
2181  process_point_1(2.1021995670317828e-01, 7.0515684111716576e-02);
2182  process_point_3(7.3254276860644785e-03,
2183  1.4932478865208237e-01,
2184  1.0290289572953278e-02);
2185  process_point_3(4.6010500165429957e-02,
2186  2.8958112563770588e-01,
2187  4.0332476640500554e-02);
2188  }
2189  else
2190  {
2191  // WV-12, 2d
2192  process_point_1(2.4646363436335583e-02, 7.9316425099736389e-03);
2193  process_point_1(4.8820375094554153e-01, 2.4266838081452032e-02);
2194  process_point_1(1.0925782765935427e-01, 2.8486052068877544e-02);
2195  process_point_1(4.4011164865859309e-01, 4.9918334928060942e-02);
2196  process_point_1(2.7146250701492608e-01, 6.2541213195902765e-02);
2197  process_point_3(2.1382490256170616e-02,
2198  1.2727971723358933e-01,
2199  1.5083677576511438e-02);
2200  process_point_3(2.3034156355267121e-02,
2201  2.9165567973834094e-01,
2202  2.1783585038607559e-02);
2203  process_point_3(1.1629601967792658e-01,
2204  2.5545422863851736e-01,
2205  4.3227363659414209e-02);
2206  }
2207  break;
2208  case 7:
2209  // There is no WV-13 rule in 3d yet
2210  Assert(dim == 2, ExcNotImplemented());
2211  if (use_odd_order)
2212  {
2213  // WV-13, 2d
2214  b_point_permutations.push_back({centroid});
2215  b_weights.push_back(6.7960036586831640e-02);
2216  process_point_1(2.1509681108843159e-02, 6.0523371035391717e-03);
2217  process_point_1(4.8907694645253935e-01, 2.3994401928894731e-02);
2218  process_point_1(4.2694141425980042e-01, 5.5601967530453329e-02);
2219  process_point_1(2.2137228629183292e-01, 5.8278485119199981e-02);
2220  process_point_3(5.1263891023823893e-03,
2221  2.7251581777342970e-01,
2222  9.5906810035432631e-03);
2223  process_point_3(2.4370186901093827e-02,
2224  1.1092204280346341e-01,
2225  1.4965401105165668e-02);
2226  process_point_3(8.7895483032197297e-02,
2227  1.6359740106785048e-01,
2228  2.4179039811593819e-02);
2229  process_point_3(6.8012243554206653e-02,
2230  3.0844176089211778e-01,
2231  3.4641276140848373e-02);
2232  }
2233  else
2234  {
2235  // WV-14, 2d
2236  process_point_1(1.9390961248701044e-02, 4.9234036024000819e-03);
2237  process_point_1(6.1799883090872587e-02, 1.4433699669776668e-02);
2238  process_point_1(4.8896391036217862e-01, 2.1883581369428889e-02);
2239  process_point_1(4.1764471934045394e-01, 3.2788353544125348e-02);
2240  process_point_1(1.7720553241254344e-01, 4.2162588736993016e-02);
2241  process_point_1(2.7347752830883865e-01, 5.1774104507291585e-02);
2242  process_point_3(1.2683309328720416e-03,
2243  1.1897449769695684e-01,
2244  5.0102288385006719e-03);
2245  process_point_3(1.4646950055654417e-02,
2246  2.9837288213625779e-01,
2247  1.4436308113533840e-02);
2248  process_point_3(5.7124757403647919e-02,
2249  1.7226668782135557e-01,
2250  2.4665753212563674e-02);
2251  process_point_3(9.2916249356971847e-02,
2252  3.3686145979634496e-01,
2253  3.8571510787060684e-02);
2254  }
2255  break;
2256  default:
2257  Assert(false, ExcNotImplemented());
2258  }
2259 
2260  Assert(b_point_permutations.size() == b_weights.size(), ExcInternalError());
2261  for (unsigned int permutation_n = 0; permutation_n < b_weights.size();
2262  ++permutation_n)
2263  {
2264  for (const std::array<double, dim + 1> &b_point :
2265  b_point_permutations[permutation_n])
2266  {
2267  const double volume = (dim == 2 ? 1.0 / 2.0 : 1.0 / 6.0);
2268  this->weights.emplace_back(volume * b_weights[permutation_n]);
2269  Point<dim> c_point;
2270  for (int d = 0; d < dim; ++d)
2271  c_point[d] = b_point[d];
2272  this->quadrature_points.emplace_back(c_point);
2273  }
2274  }
2275 }
2276 
2277 
2278 
2279 namespace
2280 {
2281  template <int dim>
2283  setup_qiterated_1D(const Quadrature<dim> &, const unsigned int)
2284  {
2285  Assert(false, ExcInternalError());
2286  return Quadrature<dim>();
2287  }
2288 
2289 
2290 
2292  setup_qiterated_1D(const Quadrature<1> &base_quad,
2293  const unsigned int n_copies)
2294  {
2295  return QIterated<1>(base_quad, n_copies);
2296  }
2297 } // namespace
2298 
2299 
2300 
2301 template <int dim>
2303  const unsigned int n_copies)
2304 {
2305  switch (dim)
2306  {
2307  case 1:
2308  static_cast<Quadrature<dim> &>(*this) =
2309  setup_qiterated_1D(base_quad, n_copies);
2310  break;
2311  case 2:
2312  case 3:
2313  {
2314  const auto n_refinements =
2315  static_cast<unsigned int>(std::round(std::log2(n_copies)));
2316  Assert((1u << n_refinements) == n_copies,
2317  ExcMessage("The number of copies must be a power of 2."));
2319  const auto reference_cell = ReferenceCells::get_simplex<dim>();
2321  tria.refine_global(n_refinements);
2322  const Mapping<dim> &mapping =
2323  reference_cell.template get_default_linear_mapping<dim>();
2325 
2326  FEValues<dim> fe_values(mapping,
2327  fe,
2328  base_quad,
2330  std::vector<Point<dim>> points;
2331  std::vector<double> weights;
2332  for (const auto &cell : tria.active_cell_iterators())
2333  {
2334  fe_values.reinit(cell);
2335  for (unsigned int qp = 0; qp < base_quad.size(); ++qp)
2336  {
2337  points.push_back(fe_values.quadrature_point(qp));
2338  weights.push_back(fe_values.JxW(qp));
2339  }
2340  }
2341 
2342  static_cast<Quadrature<dim> &>(*this) =
2343  Quadrature<dim>(points, weights);
2344 
2345  break;
2346  }
2347  default:
2348  Assert(false, ExcNotImplemented());
2349  }
2350 }
2351 
2352 
2353 
2354 template <int dim>
2355 QGaussWedge<dim>::QGaussWedge(const unsigned int n_points)
2356  : Quadrature<dim>()
2357 {
2358  AssertDimension(dim, 3);
2359 
2360  QGaussSimplex<2> quad_tri(n_points);
2361  QGauss<1> quad_line(n_points);
2362 
2363  for (unsigned int i = 0; i < quad_line.size(); ++i)
2364  for (unsigned int j = 0; j < quad_tri.size(); ++j)
2365  {
2366  this->quadrature_points.emplace_back(quad_tri.point(j)[0],
2367  quad_tri.point(j)[1],
2368  quad_line.point(i)[0]);
2369  this->weights.emplace_back(quad_tri.weight(j) * quad_line.weight(i));
2370  }
2371 
2372  AssertDimension(this->quadrature_points.size(), this->weights.size());
2373  Assert(this->quadrature_points.size() > 0,
2374  ExcMessage("No valid quadrature points!"));
2375 }
2376 
2377 
2378 
2379 template <int dim>
2380 QGaussPyramid<dim>::QGaussPyramid(const unsigned int n_points_1D)
2381  : Quadrature<dim>()
2382 {
2383  AssertDimension(dim, 3);
2384 
2385  if (n_points_1D == 1)
2386  {
2387  const double Q14 = 1.0 / 4.0;
2388  const double Q43 = 4.0 / 3.0;
2389 
2390  this->quadrature_points.emplace_back(0, 0, Q14);
2391  this->weights.emplace_back(Q43);
2392  }
2393  else if (n_points_1D == 2)
2394  {
2395  // clang-format off
2396  this->quadrature_points.emplace_back(-0.26318405556971, -0.26318405556971, 0.54415184401122);
2397  this->quadrature_points.emplace_back(-0.50661630334979, -0.50661630334979, 0.12251482265544);
2398  this->quadrature_points.emplace_back(-0.26318405556971, +0.26318405556971, 0.54415184401122);
2399  this->quadrature_points.emplace_back(-0.50661630334979, +0.50661630334979, 0.12251482265544);
2400  this->quadrature_points.emplace_back(+0.26318405556971, -0.26318405556971, 0.54415184401122);
2401  this->quadrature_points.emplace_back(+0.50661630334979, -0.50661630334979, 0.12251482265544);
2402  this->quadrature_points.emplace_back(+0.26318405556971, +0.26318405556971, 0.54415184401122);
2403  this->quadrature_points.emplace_back(+0.50661630334979, +0.50661630334979, 0.12251482265544);
2404  // clang-format on
2405 
2406  this->weights.emplace_back(0.10078588207983);
2407  this->weights.emplace_back(0.23254745125351);
2408  this->weights.emplace_back(0.10078588207983);
2409  this->weights.emplace_back(0.23254745125351);
2410  this->weights.emplace_back(0.10078588207983);
2411  this->weights.emplace_back(0.23254745125351);
2412  this->weights.emplace_back(0.10078588207983);
2413  this->weights.emplace_back(0.23254745125351);
2414  }
2415 
2416  AssertDimension(this->quadrature_points.size(), this->weights.size());
2417  Assert(this->quadrature_points.size() > 0,
2418  ExcMessage("No valid quadrature points!"));
2419 }
2420 
2421 
2422 
2423 // explicit specialization
2424 // note that 1d formulae are specialized by implementation above
2425 template class QGauss<2>;
2426 template class QGaussRadau<2>;
2427 template class QGaussLobatto<2>;
2428 template class QMidpoint<2>;
2429 template class QTrapezoid<2>;
2430 template class QSimpson<2>;
2431 template class QMilne<2>;
2432 template class QWeddle<2>;
2433 
2434 template class QGauss<3>;
2435 template class QGaussRadau<3>;
2436 template class QGaussLobatto<3>;
2437 template class QMidpoint<3>;
2438 template class QTrapezoid<3>;
2439 template class QSimpson<3>;
2440 template class QMilne<3>;
2441 template class QWeddle<3>;
2442 
2443 template class QSorted<1>;
2444 template class QSorted<2>;
2445 template class QSorted<3>;
2446 
2447 template class QTelles<1>;
2448 template class QTelles<2>;
2449 template class QTelles<3>;
2450 
2451 template class QGaussChebyshev<1>;
2452 template class QGaussChebyshev<2>;
2453 template class QGaussChebyshev<3>;
2454 
2455 template class QGaussRadauChebyshev<1>;
2456 template class QGaussRadauChebyshev<2>;
2457 template class QGaussRadauChebyshev<3>;
2458 
2459 template class QGaussLobattoChebyshev<1>;
2460 template class QGaussLobattoChebyshev<2>;
2461 template class QGaussLobattoChebyshev<3>;
2462 
2463 template class QSimplex<1>;
2464 template class QSimplex<2>;
2465 template class QSimplex<3>;
2466 
2467 template class QIteratedSimplex<1>;
2468 template class QIteratedSimplex<2>;
2469 template class QIteratedSimplex<3>;
2470 
2471 template class QSplit<1>;
2472 template class QSplit<2>;
2473 template class QSplit<3>;
2474 
2475 template class QGaussSimplex<0>;
2476 template class QGaussSimplex<1>;
2477 template class QGaussSimplex<2>;
2478 template class QGaussSimplex<3>;
2479 template class QGaussWedge<0>;
2480 template class QGaussWedge<1>;
2481 template class QGaussWedge<2>;
2482 template class QGaussWedge<3>;
2483 template class QGaussPyramid<0>;
2484 template class QGaussPyramid<1>;
2485 template class QGaussPyramid<2>;
2486 template class QGaussPyramid<3>;
2487 
2488 template class QWitherdenVincentSimplex<1>;
2489 template class QWitherdenVincentSimplex<2>;
2490 template class QWitherdenVincentSimplex<3>;
2491 
2492 #ifndef DOXYGEN
2493 template Quadrature<1>
2495  const std::array<Point<1>, 1 + 1> &vertices) const;
2496 
2497 template Quadrature<2>
2499  const std::array<Point<2>, 1 + 1> &vertices) const;
2500 
2501 template Quadrature<2>
2503  const std::array<Point<2>, 2 + 1> &vertices) const;
2504 
2505 template Quadrature<3>
2507  const std::array<Point<3>, 1 + 1> &vertices) const;
2508 
2509 template Quadrature<3>
2511  const std::array<Point<3>, 2 + 1> &vertices) const;
2512 
2513 template Quadrature<3>
2515  const std::array<Point<3>, 3 + 1> &vertices) const;
2516 
2517 template Quadrature<2>
2519  const std::vector<std::array<Point<2>, 1 + 1>> &simplices) const;
2520 
2521 template Quadrature<3>
2523  const std::vector<std::array<Point<3>, 1 + 1>> &simplices) const;
2524 
2525 template Quadrature<2>
2527  const std::vector<std::array<Point<2>, 2 + 1>> &simplices) const;
2528 
2529 template Quadrature<3>
2531  const std::vector<std::array<Point<3>, 2 + 1>> &simplices) const;
2532 
2533 template Quadrature<3>
2535  const std::vector<std::array<Point<3>, 3 + 1>> &simplices) const;
2536 #endif
2537 
Tensor< 1, spacedim, typename ProductType< Number1, Number2 >::type > apply_transformation(const DerivativeForm< 1, dim, spacedim, Number1 > &grad_F, const Tensor< 1, dim, Number2 > &d_x)
DerivativeForm< 1, spacedim, dim, Number > transpose() const
const Point< spacedim > & quadrature_point(const unsigned int q_point) const
double JxW(const unsigned int q_point) const
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access >> &cell)
Abstract base class for mapping classes.
Definition: mapping.h:317
QDuffy(const Quadrature< 1 > &radial_quadrature, const Quadrature< 1 > &angular_quadrature, const double beta=1.0)
QGaussChebyshev(const unsigned int n)
Generate a formula with n quadrature points.
QGaussLobattoChebyshev(const unsigned int n)
Generate a formula with n quadrature points.
QGaussLobatto(const unsigned int n)
const double fraction
QGaussLogR(const unsigned int n, const Point< dim > &x0=Point< dim >(), const double alpha=1, const bool factor_out_singular_weight=false)
QGaussLog(const unsigned int n, const bool revert=false)
static std::vector< double > get_quadrature_points(const unsigned int n)
static std::vector< double > get_quadrature_weights(const unsigned int n)
QGaussOneOverR(const unsigned int n, const Point< dim > &singularity, const bool factor_out_singular_weight=false)
static unsigned int quad_size(const Point< dim > &singularity, const unsigned int n)
QGaussPyramid(const unsigned int n_points_1D)
QGaussRadauChebyshev(const unsigned int n, const EndPoint end_point=QGaussRadauChebyshev::EndPoint::left)
Generate a formula with n quadrature points.
const EndPoint end_point
QGaussRadau(const unsigned int n, const EndPoint end_point=QGaussRadau::EndPoint::left)
QGaussSimplex(const unsigned int n_points_1D)
QGaussWedge(const unsigned int n_points_1D)
QGauss(const unsigned int n)
QIteratedSimplex(const Quadrature< dim > &base_quadrature, const unsigned int n_copies)
Quadrature< spacedim > compute_affine_transformation(const std::array< Point< spacedim >, dim+1 > &vertices) const
QSimplex(const Quadrature< dim > &quad)
Quadrature< spacedim > mapped_quadrature(const std::vector< std::array< Point< spacedim >, dim+1 >> &simplices) const
QSorted(const Quadrature< dim > &quad)
bool compare_weights(const unsigned int a, const unsigned int b) const
QSplit(const QSimplex< dim > &base, const Point< dim > &split_point)
QTelles(const Quadrature< 1 > &base_quad, const Point< dim > &singularity)
QTrianglePolar(const Quadrature< 1 > &radial_quadrature, const Quadrature< 1 > &angular_quadrature)
QWitherdenVincentSimplex(const unsigned int n_points_1D, const bool use_odd_order=true)
std::vector< Point< dim > > quadrature_points
Definition: quadrature.h:339
Quadrature & operator=(const Quadrature< dim > &)
Definition: quadrature.cc:297
bool is_tensor_product_flag
Definition: quadrature.h:354
double weight(const unsigned int i) const
const Point< dim > & point(const unsigned int i) const
std::vector< double > weights
Definition: quadrature.h:345
unsigned int size() const
Definition: tensor.h:516
void refine_global(const unsigned int times=1)
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
Point< 3 > vertices[4]
IteratorRange< active_cell_iterator > active_cell_iterators() const
static ::ExceptionBase & ExcInternalError()
#define Assert(cond, exc)
Definition: exceptions.h:1631
static ::ExceptionBase & ExcNotImplemented()
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1820
#define AssertIndexRange(index, range)
Definition: exceptions.h:1888
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1732
@ update_JxW_values
Transformed quadrature weights.
@ update_quadrature_points
Transformed quadrature points.
void simplex(Triangulation< dim, dim > &tria, const std::vector< Point< dim >> &vertices)
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
double volume(const Triangulation< dim, spacedim > &tria)
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:472
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:192
void quadrature_points(const Triangulation< dim, spacedim > &triangulation, const Quadrature< dim > &quadrature, const std::vector< std::vector< BoundingBox< spacedim >>> &global_bounding_boxes, ParticleHandler< dim, spacedim > &particle_handler, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()), const std::vector< std::vector< double >> &properties={})
Definition: generators.cc:490
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
Number jacobi_polynomial_value(const unsigned int degree, const int alpha, const int beta, const Number x)
Definition: polynomial.h:1039
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:480
std::vector< double > get_quadrature_weights(const unsigned int n)
std::vector< double > get_quadrature_points(const unsigned int n)
std::vector< double > get_quadrature_weights(const unsigned int n)
std::vector< double > get_quadrature_points(const unsigned int n)
std::vector< long double > compute_quadrature_weights(const std::vector< long double > &x, const int alpha, const int beta)
long double gamma(const unsigned int n)
std::vector< double > get_quadrature_points(const unsigned int n, const ::QGaussRadauChebyshev< 1 >::EndPoint end_point)
std::vector< double > get_quadrature_weights(const unsigned int n, const ::QGaussRadauChebyshev< 1 >::EndPoint end_point)
std::vector< double > get_quadrature_points(const unsigned int n, const ::QGaussRadau< 1 >::EndPoint end_point)
std::vector< double > get_quadrature_weights(const unsigned int n, const ::QGaussRadau< 1 >::EndPoint end_point)
std::vector< double > get_left_quadrature_weights(const unsigned int n)
std::vector< double > get_left_quadrature_points(const unsigned int n)
void split_point(const Point< dim1+dim2 > &source, Point< dim1 > &p1, Point< dim2 > &p2)
static constexpr double PI_2
Definition: numbers.h:265
static constexpr double PI
Definition: numbers.h:260
void transform(const InputIterator &begin_in, const InputIterator &end_in, OutputIterator out, const Function &function, const unsigned int grainsize)
Definition: parallel.h:166
static Point< dim > unit_cell_vertex(const unsigned int vertex)
const ::Triangulation< dim, spacedim > & tria