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