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