Reference documentation for deal.II version Git 3f1f337db3 2021-10-23 13:19:02 -0600
\(\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>
1232  const std::array<Point<dim>, dim + 1> &vertices) const
1233 {
1234  Tensor<2, dim> B;
1235  for (unsigned int d = 0; d < dim; ++d)
1236  B[d] = vertices[d + 1] - vertices[0];
1237 
1238  B = transpose(B);
1239  const double J = std::abs(determinant(B));
1240 
1241  // if the determinant is zero, we return an empty quadrature
1242  if (J < 1e-12)
1243  return Quadrature<dim>();
1244 
1245  std::vector<Point<dim>> qp(this->size());
1246  std::vector<double> w(this->size());
1247 
1248  for (unsigned int i = 0; i < this->size(); ++i)
1249  {
1250  qp[i] = Point<dim>(vertices[0] + B * this->point(i));
1251  w[i] = J * this->weight(i);
1252  }
1253 
1254  return Quadrature<dim>(qp, w);
1255 }
1256 
1257 
1258 
1260  const Quadrature<1> &angular_quadrature)
1261  : QSimplex<2>(Quadrature<2>())
1262 {
1263  const QAnisotropic<2> base(radial_quadrature, angular_quadrature);
1264  this->quadrature_points.resize(base.size());
1265  this->weights.resize(base.size());
1266  for (unsigned int i = 0; i < base.size(); ++i)
1267  {
1268  const auto &q = base.point(i);
1269  const auto w = base.weight(i);
1270 
1271  const auto xhat = q[0];
1272  const auto yhat = q[1];
1273 
1274  const double t = numbers::PI_2 * yhat;
1275  const double pi = numbers::PI;
1276  const double st = std::sin(t);
1277  const double ct = std::cos(t);
1278  const double r = xhat / (st + ct);
1279 
1280  const double J = pi * xhat / (2 * (std::sin(pi * yhat) + 1));
1281 
1282  this->quadrature_points[i] = Point<2>(r * ct, r * st);
1283  this->weights[i] = w * J;
1284  }
1285 }
1286 
1287 
1288 
1289 QTrianglePolar::QTrianglePolar(const unsigned int n)
1290  : QTrianglePolar(QGauss<1>(n), QGauss<1>(n))
1291 {}
1292 
1293 
1294 
1295 QDuffy::QDuffy(const Quadrature<1> &radial_quadrature,
1296  const Quadrature<1> &angular_quadrature,
1297  const double beta)
1298  : QSimplex<2>(Quadrature<2>())
1299 {
1300  const QAnisotropic<2> base(radial_quadrature, angular_quadrature);
1301  this->quadrature_points.resize(base.size());
1302  this->weights.resize(base.size());
1303  for (unsigned int i = 0; i < base.size(); ++i)
1304  {
1305  const auto &q = base.point(i);
1306  const auto w = base.weight(i);
1307 
1308  const auto xhat = q[0];
1309  const auto yhat = q[1];
1310 
1311  const double x = std::pow(xhat, beta) * (1 - yhat);
1312  const double y = std::pow(xhat, beta) * yhat;
1313 
1314  const double J = beta * std::pow(xhat, 2. * beta - 1.);
1315 
1316  this->quadrature_points[i] = Point<2>(x, y);
1317  this->weights[i] = w * J;
1318  }
1319 }
1320 
1321 
1322 
1323 QDuffy::QDuffy(const unsigned int n, const double beta)
1324  : QDuffy(QGauss<1>(n), QGauss<1>(n), beta)
1325 {}
1326 
1327 
1328 
1329 template <int dim>
1331 {
1333  ExcMessage(
1334  "The split point should be inside the unit reference cell."));
1335 
1336  std::array<Point<dim>, dim + 1> vertices;
1337  vertices[0] = split_point;
1338 
1339  // Make a simplex from the split_point and the first dim vertices of each
1340  // face. In dimension three, we need to split the face in two triangles, so
1341  // we use once the first dim vertices of each face, and the second time the
1342  // the dim vertices of each face starting from 1.
1343  for (auto f : GeometryInfo<dim>::face_indices())
1344  for (unsigned int start = 0; start < (dim > 2 ? 2 : 1); ++start)
1345  {
1346  for (unsigned int i = 0; i < dim; ++i)
1347  vertices[i + 1] = GeometryInfo<dim>::unit_cell_vertex(
1349  const auto quad = base.compute_affine_transformation(vertices);
1350  if (quad.size())
1351  {
1352  this->quadrature_points.insert(this->quadrature_points.end(),
1353  quad.get_points().begin(),
1354  quad.get_points().end());
1355  this->weights.insert(this->weights.end(),
1356  quad.get_weights().begin(),
1357  quad.get_weights().end());
1358  }
1359  }
1360 }
1361 
1362 
1363 
1364 template <int dim>
1365 QGaussSimplex<dim>::QGaussSimplex(const unsigned int n_points_1D)
1366  : QSimplex<dim>(Quadrature<dim>())
1367 {
1368  // fill quadrature points and quadrature weights
1369  if (dim == 0 || dim == 1)
1370  {
1371  const ::QGauss<dim> quad(n_points_1D);
1372 
1373  this->quadrature_points = quad.get_points();
1374  this->weights = quad.get_weights();
1375  }
1376  else if (dim == 2)
1377  {
1378  if (n_points_1D == 1)
1379  {
1380  const double p = 1.0 / 3.0;
1381  this->quadrature_points.emplace_back(p, p);
1382  this->weights.emplace_back(0.5);
1383  }
1384  else if (n_points_1D == 2)
1385  {
1386  // The Hillion 7 scheme, as communicated by quadpy
1387  //
1388  // See: Numerical Integration on a Triangle, International Journal for
1389  // Numerical Methods in Engineering, 1977
1390  const double Q12 = 1.0 / 2.0;
1391  this->quadrature_points.emplace_back(0.17855872826361643,
1392  0.1550510257216822);
1393  this->quadrature_points.emplace_back(0.07503111022260812,
1394  0.6449489742783178);
1395  this->quadrature_points.emplace_back(0.6663902460147014,
1396  0.1550510257216822);
1397  this->quadrature_points.emplace_back(0.28001991549907407,
1398  0.6449489742783178);
1399 
1400  this->weights.emplace_back(0.31804138174397717 * Q12);
1401  this->weights.emplace_back(0.18195861825602283 * Q12);
1402  this->weights.emplace_back(0.31804138174397717 * Q12);
1403  this->weights.emplace_back(0.18195861825602283 * Q12);
1404  }
1405  else if (n_points_1D == 3)
1406  {
1407  // The Hammer-Marlowe-Stroud 5 Scheme, as communicated by quadpy
1408  const double p0 = 2.0 / 7.0 - std::sqrt(15.0) / 21.0;
1409  const double p1 = 2.0 / 7.0 + std::sqrt(15.0) / 21.0;
1410  const double p2 = 3.0 / 7.0 - 2.0 * std::sqrt(15.0) / 21.0;
1411  const double p3 = 3.0 / 7.0 + 2.0 * std::sqrt(15.0) / 21.0;
1412  this->quadrature_points.emplace_back(1.0 / 3.0, 1.0 / 3.0);
1413  this->quadrature_points.emplace_back(p3, p0);
1414  this->quadrature_points.emplace_back(p0, p3);
1415  this->quadrature_points.emplace_back(p0, p0);
1416  this->quadrature_points.emplace_back(p2, p1);
1417  this->quadrature_points.emplace_back(p1, p2);
1418  this->quadrature_points.emplace_back(p1, p1);
1419 
1420  const double q12 = 0.5;
1421  const double w0 = 9.0 / 40.0;
1422  const double w1 = 31.0 / 240.0 - std::sqrt(15.0) / 1200.0;
1423  const double w2 = 31.0 / 240.0 + std::sqrt(15.0) / 1200.0;
1424  this->weights.emplace_back(q12 * w0);
1425  this->weights.emplace_back(q12 * w1);
1426  this->weights.emplace_back(q12 * w1);
1427  this->weights.emplace_back(q12 * w1);
1428  this->weights.emplace_back(q12 * w2);
1429  this->weights.emplace_back(q12 * w2);
1430  this->weights.emplace_back(q12 * w2);
1431  }
1432  else if (n_points_1D == 4)
1433  {
1435  QWitherdenVincentSimplex<dim>(n_points_1D));
1436  }
1437  }
1438  else if (dim == 3)
1439  {
1440  if (n_points_1D == 1)
1441  {
1442  const double Q14 = 1.0 / 4.0;
1443  const double Q16 = 1.0 / 6.0;
1444 
1445  this->quadrature_points.emplace_back(Q14, Q14, Q14);
1446  this->weights.emplace_back(Q16);
1447  }
1448  // The Xiao Gimbutas 03 scheme, as communicated by quadpy
1449  //
1450  // See: A numerical algorithm for the construction of efficient quadrature
1451  // rules in two and higher dimensions, Computers & Mathematics with
1452  // Applications, 2010
1453  else if (n_points_1D == 2)
1454  {
1455  const double Q16 = 1.0 / 6.0;
1456  this->weights.emplace_back(0.1223220027573451 * Q16);
1457  this->weights.emplace_back(0.1280664127107469 * Q16);
1458  this->weights.emplace_back(0.1325680271444452 * Q16);
1459  this->weights.emplace_back(0.1406244096604032 * Q16);
1460  this->weights.emplace_back(0.2244151669175574 * Q16);
1461  this->weights.emplace_back(0.2520039808095023 * Q16);
1462 
1463  this->quadrature_points.emplace_back(0.1620014916985245,
1464  0.1838503504920977,
1465  0.01271836631368145);
1466  this->quadrature_points.emplace_back(0.01090521221118924,
1467  0.2815238021235462,
1468  0.3621268299455338);
1469  this->quadrature_points.emplace_back(0.1901170024392839,
1470  0.01140332944455717,
1471  0.3586207204668839);
1472  this->quadrature_points.emplace_back(0.170816925164989,
1473  0.1528181430909273,
1474  0.6384932999617267);
1475  this->quadrature_points.emplace_back(0.1586851632274406,
1476  0.5856628056552158,
1477  0.1308471689520965);
1478  this->quadrature_points.emplace_back(0.5712260521491151,
1479  0.1469183900871696,
1480  0.1403728057942107);
1481  }
1482  // Past this point the best rules (positive weights, minimal number of
1483  // points) we have right now are the Witherden-Vincent ones
1484  else if (n_points_1D == 3)
1485  {
1487  QWitherdenVincentSimplex<dim>(n_points_1D));
1488  }
1489  else if (n_points_1D == 4)
1490  {
1492  QWitherdenVincentSimplex<dim>(n_points_1D));
1493  }
1494  }
1495 
1496  AssertDimension(this->quadrature_points.size(), this->weights.size());
1497  Assert(this->quadrature_points.size() > 0,
1499  "QGaussSimplex is currently only implemented for "
1500  "n_points_1D = 1, 2, 3, and 4 while you are asking for "
1501  "n_points_1D = " +
1502  Utilities::to_string(n_points_1D)));
1503 }
1504 
1505 namespace
1506 {
1507  template <std::size_t b_dim>
1508  std::vector<std::array<double, b_dim>>
1509  all_permutations(const std::array<double, b_dim> &b_point)
1510  {
1511  std::vector<std::array<double, b_dim>> output;
1512 
1513  // We want all possible permutations of the barycentric coordinates.
1514  // The easiest way to get all of them is to sort the input first and
1515  // then use next_permutation to cycle through them all.
1516  std::array<double, b_dim> temp = b_point;
1517  std::sort(temp.begin(), temp.end());
1518  do
1519  {
1520  output.push_back(temp);
1521  }
1522  while (std::next_permutation(temp.begin(), temp.end()));
1523 
1524  return output;
1525  }
1526 } // namespace
1527 
1528 
1529 
1530 template <int dim>
1532  const unsigned int n_points_1D,
1533  const bool use_odd_order)
1534  : QSimplex<dim>(Quadrature<dim>())
1535 {
1536  Assert(1 <= dim && dim <= 3, ExcNotImplemented());
1537  // Just use Gauss in 1D: this is a high-order open rule so this is a
1538  // reasonable equivalent for generic programming.
1539  if (dim == 1)
1540  {
1542  return;
1543  }
1544 
1545  std::array<double, dim + 1> centroid;
1546  std::fill(centroid.begin(), centroid.end(), 1.0 / (dim + 1.0));
1547  std::vector<std::vector<std::array<double, dim + 1>>> b_point_permutations;
1548  std::vector<double> b_weights;
1549 
1550  // We can simplify the implementation of these quadrature rules
1551  // by quite a bit by exploiting symmetry - we do essentially the
1552  // same thing for each barycentric coordinate, so we can express
1553  // our quadrature rule as permutations of barycentric points
1554  // instead of writing things out explicitly.
1555 
1556  // Apply a Barycentric permutation where one point is different.
1557  // Equivalent to d3_aa and s31 in quadpy.
1558  auto process_point_1 = [&](const double a, const double w) {
1559  const double b = 1.0 - dim * a;
1560  std::array<double, dim + 1> b_point;
1561  std::fill(b_point.begin(), b_point.begin() + dim, a);
1562  b_point[dim] = b;
1563 
1564  b_weights.push_back(w);
1565  b_point_permutations.push_back(all_permutations(b_point));
1566  };
1567 
1568  // Apply a Barycentric permutation where two points (in 3D) are different.
1569  // Equivalent to s22 in quadpy.
1570  auto process_point_2 = [&](const double a, const double w) {
1571  Assert(dim == 3, ExcInternalError());
1572  const double b = (1.0 - 2.0 * a) / 2.0;
1573  std::array<double, dim + 1> b_point;
1574  std::fill(b_point.begin(), b_point.begin() + dim - 1, a);
1575  b_point[dim - 1] = b;
1576  b_point[dim] = b;
1577 
1578  b_weights.push_back(w);
1579  b_point_permutations.push_back(all_permutations(b_point));
1580  };
1581 
1582  // Apply a Barycentric permutation where three (or four) points
1583  // are different (since there are two inputs).
1584  // Equivalent to d3_ab and s211 in quadpy.
1585  auto process_point_3 = [&](const double a, const double b, const double w) {
1586  const double c = 1.0 - (dim - 1.0) * a - b;
1587  std::array<double, dim + 1> b_point;
1588  std::fill(b_point.begin(), b_point.begin() + dim - 1, a);
1589  b_point[dim - 1] = b;
1590  b_point[dim] = c;
1591 
1592  b_weights.push_back(w);
1593  b_point_permutations.push_back(all_permutations(b_point));
1594  };
1595 
1596  switch (n_points_1D)
1597  {
1598  case 1:
1599  switch (dim)
1600  {
1601  case 2:
1602  if (use_odd_order)
1603  {
1604  // WV-1, 2D
1605  b_point_permutations.push_back({centroid});
1606  b_weights.push_back(1.0000000000000000e+00);
1607  }
1608  else
1609  {
1610  // WV-2, 2D
1611  process_point_1(1.6666666666666669e-01,
1612  3.3333333333333331e-01);
1613  }
1614  break;
1615  case 3:
1616  if (use_odd_order)
1617  {
1618  // WV-1, 3D
1619  b_point_permutations.push_back({centroid});
1620  b_weights.push_back(1.0000000000000000e+00);
1621  }
1622  else
1623  {
1624  // WV-2, 3D
1625  process_point_1(1.3819660112501050e-01,
1626  2.5000000000000000e-01);
1627  }
1628  break;
1629  default:
1630  Assert(false, ExcNotImplemented());
1631  }
1632  break;
1633  case 2:
1634  switch (dim)
1635  {
1636  case 2:
1637  // WV-4 in both cases (no WV-3 in 2D)
1638  process_point_1(9.1576213509770743e-02, 1.0995174365532187e-01);
1639  process_point_1(4.4594849091596489e-01, 2.2338158967801147e-01);
1640  break;
1641  case 3:
1642  if (use_odd_order)
1643  {
1644  // WV-3, 3D
1645  process_point_1(3.2816330251638171e-01,
1646  1.3621784253708741e-01);
1647  process_point_1(1.0804724989842859e-01,
1648  1.1378215746291261e-01);
1649  }
1650  else
1651  {
1652  // WV-5 (no WV-4 in 3D)
1654  }
1655  break;
1656  default:
1657  Assert(false, ExcInternalError());
1658  }
1659  break;
1660  case 3:
1661  switch (dim)
1662  {
1663  case 2:
1664  if (use_odd_order)
1665  {
1666  // WV-5, 2D
1667  b_point_permutations.push_back({centroid});
1668  b_weights.push_back(2.2500000000000001e-01);
1669  process_point_1(1.0128650732345634e-01,
1670  1.2593918054482714e-01);
1671  process_point_1(4.7014206410511511e-01,
1672  1.3239415278850619e-01);
1673  }
1674  else
1675  {
1676  // WV-6, 2D
1677  process_point_1(6.3089014491502227e-02,
1678  5.0844906370206819e-02);
1679  process_point_1(2.4928674517091043e-01,
1680  1.1678627572637937e-01);
1681  process_point_3(5.3145049844816938e-02,
1682  3.1035245103378439e-01,
1683  8.2851075618373571e-02);
1684  }
1685  break;
1686  case 3:
1687  if (use_odd_order)
1688  {
1689  // WV-5, 3D
1690  process_point_1(3.1088591926330061e-01,
1691  1.1268792571801590e-01);
1692  process_point_1(9.2735250310891248e-02,
1693  7.3493043116361956e-02);
1694  process_point_2(4.5503704125649642e-02,
1695  4.2546020777081472e-02);
1696  }
1697  else
1698  {
1699  // WV-6, 3D
1700  process_point_1(4.0673958534611372e-02,
1701  1.0077211055320640e-02);
1702  process_point_1(3.2233789014227548e-01,
1703  5.5357181543654717e-02);
1704  process_point_1(2.1460287125915201e-01,
1705  3.9922750258167487e-02);
1706  process_point_3(6.3661001875017442e-02,
1707  6.0300566479164919e-01,
1708  4.8214285714285710e-02);
1709  }
1710  break;
1711  default:
1712  Assert(false, ExcInternalError());
1713  }
1714  break;
1715  case 4:
1716  switch (dim)
1717  {
1718  case 2:
1719  if (use_odd_order)
1720  {
1721  // WV-7, 2D
1722  process_point_1(3.3730648554587850e-02,
1723  1.6545050110792131e-02);
1724  process_point_1(4.7430969250471822e-01,
1725  7.7086646185986069e-02);
1726  process_point_1(2.4157738259540357e-01,
1727  1.2794417123015558e-01);
1728  process_point_3(4.7036644652595216e-02,
1729  1.9868331479735168e-01,
1730  5.5878732903199779e-02);
1731  }
1732  else
1733  {
1734  // WV-8, 2D
1735  b_point_permutations.push_back({centroid});
1736  b_weights.push_back(1.4431560767778717e-01);
1737  process_point_1(5.0547228317030957e-02,
1738  3.2458497623198079e-02);
1739  process_point_1(4.5929258829272313e-01,
1740  9.5091634267284619e-02);
1741  process_point_1(1.7056930775176021e-01,
1742  1.0321737053471824e-01);
1743  process_point_3(8.3947774099575878e-03,
1744  2.6311282963463811e-01,
1745  2.7230314174434993e-02);
1746  }
1747  break;
1748  case 3:
1749  if (use_odd_order)
1750  {
1751  // WV-7, 3D
1752  b_point_permutations.push_back({centroid});
1753  b_weights.push_back(9.5485289464130846e-02);
1754  process_point_1(3.1570114977820279e-01,
1755  4.2329581209967028e-02);
1756  process_point_2(5.0489822598396350e-02,
1757  3.1896927832857580e-02);
1758  process_point_3(1.8883383102600099e-01,
1759  5.7517163758699996e-01,
1760  3.7207130728334620e-02);
1761  process_point_3(2.1265472541483140e-02,
1762  8.1083024109854862e-01,
1763  8.1107708299033420e-03);
1764  }
1765  else
1766  {
1767  // WV-8, 3D
1768  process_point_1(1.0795272496221089e-01,
1769  2.6426650908408830e-02);
1770  process_point_1(1.8510948778258660e-01,
1771  5.2031747563738531e-02);
1772  process_point_1(4.2316543684767283e-02,
1773  7.5252561535401989e-03);
1774  process_point_1(3.1418170912403898e-01,
1775  4.1763782856934897e-02);
1776  process_point_2(4.3559132858383021e-01,
1777  3.6280930261308818e-02);
1778  process_point_3(2.1433930127130570e-02,
1779  7.1746406342630831e-01,
1780  7.1569028908444327e-03);
1781  process_point_3(2.0413933387602909e-01,
1782  5.8379737830214440e-01,
1783  1.5453486150960340e-02);
1784  }
1785  break;
1786  default:
1787  Assert(false, ExcInternalError());
1788  }
1789  break;
1790  case 5:
1791  switch (dim)
1792  {
1793  case 2:
1794  if (use_odd_order)
1795  {
1796  // WV-9, 2D
1797  b_point_permutations.push_back({centroid});
1798  b_weights.push_back(9.7135796282798836e-02);
1799  process_point_1(4.4729513394452691e-02,
1800  2.5577675658698031e-02);
1801  process_point_1(4.8968251919873762e-01,
1802  3.1334700227139071e-02);
1803  process_point_1(4.3708959149293664e-01,
1804  7.7827541004774278e-02);
1805  process_point_1(1.8820353561903275e-01,
1806  7.9647738927210249e-02);
1807  process_point_3(3.6838412054736258e-02,
1808  2.2196298916076568e-01,
1809  4.3283539377289376e-02);
1810  }
1811  else
1812  {
1813  // WV-10, 2D
1814  b_point_permutations.push_back({centroid});
1815  b_weights.push_back(8.1743329146285973e-02);
1816  process_point_1(3.2055373216943517e-02,
1817  1.3352968813149567e-02);
1818  process_point_1(1.4216110105656438e-01,
1819  4.5957963604744731e-02);
1820  process_point_3(2.8367665339938453e-02,
1821  1.6370173373718250e-01,
1822  2.5297757707288385e-02);
1823  process_point_3(2.9619889488729734e-02,
1824  3.6914678182781102e-01,
1825  3.4184648162959429e-02);
1826  process_point_3(1.4813288578382056e-01,
1827  3.2181299528883545e-01,
1828  6.3904906396424044e-02);
1829  }
1830  break;
1831  case 3:
1832  if (use_odd_order)
1833  {
1834  // WV-9, 3D
1835  b_point_permutations.push_back({centroid});
1836  b_weights.push_back(5.8010548912480253e-02);
1837  process_point_1(6.1981697552226933e-10,
1838  6.4319281759256394e-05);
1839  process_point_1(1.6077453539526160e-01,
1840  2.3173338462425461e-02);
1841  process_point_1(3.2227652182142102e-01,
1842  2.9562912335429289e-02);
1843  process_point_1(4.5108918345413578e-02,
1844  8.0639799796161822e-03);
1845  process_point_2(1.1229654600437609e-01,
1846  3.8134080103702457e-02);
1847  process_point_3(4.5887144875245922e-01,
1848  2.5545792330413102e-03,
1849  8.3844221982985519e-03);
1850  process_point_3(3.3775870685338598e-02,
1851  7.1835032644207453e-01,
1852  1.0234559352745330e-02);
1853  process_point_3(1.8364136980992790e-01,
1854  3.4415910578175279e-02,
1855  2.0524915967988139e-02);
1856  }
1857  else
1858  {
1859  // WV-10, 3D
1860  b_point_permutations.push_back({centroid});
1861  b_weights.push_back(4.7399773556020743e-02);
1862  process_point_1(3.1225006869518868e-01,
1863  2.6937059992268701e-02);
1864  process_point_1(1.1430965385734609e-01,
1865  9.8691597167933822e-03);
1866  process_point_3(4.1043073921896539e-01,
1867  1.6548602561961109e-01,
1868  1.1393881220195230e-02);
1869  process_point_3(6.1380088247906528e-03,
1870  9.4298876734520487e-01,
1871  3.6194434433925362e-04);
1872  process_point_3(1.2105018114558939e-01,
1873  4.7719037990428043e-01,
1874  2.5739731980456069e-02);
1875  process_point_3(3.2779468216442620e-02,
1876  5.9425626948000698e-01,
1877  1.0135871679755789e-02);
1878  process_point_3(3.2485281564823047e-02,
1879  8.0117728465834437e-01,
1880  6.5761472770359038e-03);
1881  process_point_3(1.7497934218393901e-01,
1882  6.2807184547536599e-01,
1883  1.2907035798861989e-02);
1884  }
1885  break;
1886  default:
1887  Assert(false, ExcNotImplemented());
1888  }
1889  break;
1890  case 6:
1891  // There is no WV-11 rule in 3D yet
1892  Assert(dim == 2, ExcNotImplemented());
1893  if (use_odd_order)
1894  {
1895  // WV-11, 2D
1896  b_point_permutations.push_back({centroid});
1897  b_weights.push_back(8.5761179732224219e-02);
1898  process_point_1(2.8485417614371900e-02, 1.0431870512894697e-02);
1899  process_point_1(4.9589190096589092e-01, 1.6606273054585369e-02);
1900  process_point_1(1.0263548271224643e-01, 3.8630759237019321e-02);
1901  process_point_1(4.3846592676435220e-01, 6.7316154079468296e-02);
1902  process_point_1(2.1021995670317828e-01, 7.0515684111716576e-02);
1903  process_point_3(7.3254276860644785e-03,
1904  1.4932478865208237e-01,
1905  1.0290289572953278e-02);
1906  process_point_3(4.6010500165429957e-02,
1907  2.8958112563770588e-01,
1908  4.0332476640500554e-02);
1909  }
1910  else
1911  {
1912  // WV-12, 2D
1913  process_point_1(2.4646363436335583e-02, 7.9316425099736389e-03);
1914  process_point_1(4.8820375094554153e-01, 2.4266838081452032e-02);
1915  process_point_1(1.0925782765935427e-01, 2.8486052068877544e-02);
1916  process_point_1(4.4011164865859309e-01, 4.9918334928060942e-02);
1917  process_point_1(2.7146250701492608e-01, 6.2541213195902765e-02);
1918  process_point_3(2.1382490256170616e-02,
1919  1.2727971723358933e-01,
1920  1.5083677576511438e-02);
1921  process_point_3(2.3034156355267121e-02,
1922  2.9165567973834094e-01,
1923  2.1783585038607559e-02);
1924  process_point_3(1.1629601967792658e-01,
1925  2.5545422863851736e-01,
1926  4.3227363659414209e-02);
1927  }
1928  break;
1929  case 7:
1930  // There is no WV-13 rule in 3D yet
1931  Assert(dim == 2, ExcNotImplemented());
1932  if (use_odd_order)
1933  {
1934  // WV-13, 2D
1935  b_point_permutations.push_back({centroid});
1936  b_weights.push_back(6.7960036586831640e-02);
1937  process_point_1(2.1509681108843159e-02, 6.0523371035391717e-03);
1938  process_point_1(4.8907694645253935e-01, 2.3994401928894731e-02);
1939  process_point_1(4.2694141425980042e-01, 5.5601967530453329e-02);
1940  process_point_1(2.2137228629183292e-01, 5.8278485119199981e-02);
1941  process_point_3(5.1263891023823893e-03,
1942  2.7251581777342970e-01,
1943  9.5906810035432631e-03);
1944  process_point_3(2.4370186901093827e-02,
1945  1.1092204280346341e-01,
1946  1.4965401105165668e-02);
1947  process_point_3(8.7895483032197297e-02,
1948  1.6359740106785048e-01,
1949  2.4179039811593819e-02);
1950  process_point_3(6.8012243554206653e-02,
1951  3.0844176089211778e-01,
1952  3.4641276140848373e-02);
1953  }
1954  else
1955  {
1956  // WV-14, 2D
1957  process_point_1(1.9390961248701044e-02, 4.9234036024000819e-03);
1958  process_point_1(6.1799883090872587e-02, 1.4433699669776668e-02);
1959  process_point_1(4.8896391036217862e-01, 2.1883581369428889e-02);
1960  process_point_1(4.1764471934045394e-01, 3.2788353544125348e-02);
1961  process_point_1(1.7720553241254344e-01, 4.2162588736993016e-02);
1962  process_point_1(2.7347752830883865e-01, 5.1774104507291585e-02);
1963  process_point_3(1.2683309328720416e-03,
1964  1.1897449769695684e-01,
1965  5.0102288385006719e-03);
1966  process_point_3(1.4646950055654417e-02,
1967  2.9837288213625779e-01,
1968  1.4436308113533840e-02);
1969  process_point_3(5.7124757403647919e-02,
1970  1.7226668782135557e-01,
1971  2.4665753212563674e-02);
1972  process_point_3(9.2916249356971847e-02,
1973  3.3686145979634496e-01,
1974  3.8571510787060684e-02);
1975  }
1976  break;
1977  default:
1978  Assert(false, ExcNotImplemented());
1979  }
1980 
1981  Assert(b_point_permutations.size() == b_weights.size(), ExcInternalError());
1982  for (unsigned int permutation_n = 0; permutation_n < b_weights.size();
1983  ++permutation_n)
1984  {
1985  for (const std::array<double, dim + 1> &b_point :
1986  b_point_permutations[permutation_n])
1987  {
1988  const double volume = (dim == 2 ? 1.0 / 2.0 : 1.0 / 6.0);
1989  this->weights.emplace_back(volume * b_weights[permutation_n]);
1990  Point<dim> c_point;
1991  std::copy(b_point.begin(),
1992  b_point.begin() + dim,
1993  c_point.begin_raw());
1994  this->quadrature_points.emplace_back(c_point);
1995  }
1996  }
1997 }
1998 
1999 
2000 
2001 namespace
2002 {
2003  template <int dim>
2005  setup_qiterated_1D(const Quadrature<dim> &, const unsigned int)
2006  {
2007  Assert(false, ExcInternalError());
2008  return Quadrature<dim>();
2009  }
2010 
2011 
2012 
2014  setup_qiterated_1D(const Quadrature<1> &base_quad,
2015  const unsigned int n_copies)
2016  {
2017  return QIterated<1>(base_quad, n_copies);
2018  }
2019 } // namespace
2020 
2021 
2022 
2023 template <int dim>
2025  const unsigned int n_copies)
2026 {
2027  switch (dim)
2028  {
2029  case 1:
2030  static_cast<Quadrature<dim> &>(*this) =
2031  setup_qiterated_1D(base_quad, n_copies);
2032  break;
2033  case 2:
2034  case 3:
2035  {
2036  const auto n_refinements =
2037  static_cast<unsigned int>(std::round(std::log2(n_copies)));
2038  Assert((1u << n_refinements) == n_copies,
2039  ExcMessage("The number of copies must be a power of 2."));
2041  const auto reference_cell = ReferenceCells::get_simplex<dim>();
2043  tria.refine_global(n_refinements);
2044  const Mapping<dim> &mapping =
2045  reference_cell.template get_default_linear_mapping<dim>();
2047 
2048  FEValues<dim> fe_values(mapping,
2049  fe,
2050  base_quad,
2052  std::vector<Point<dim>> points;
2053  std::vector<double> weights;
2054  for (const auto &cell : tria.active_cell_iterators())
2055  {
2056  fe_values.reinit(cell);
2057  for (unsigned int qp = 0; qp < base_quad.size(); ++qp)
2058  {
2059  points.push_back(fe_values.quadrature_point(qp));
2060  weights.push_back(fe_values.JxW(qp));
2061  }
2062  }
2063 
2064  static_cast<Quadrature<dim> &>(*this) =
2065  Quadrature<dim>(points, weights);
2066 
2067  break;
2068  }
2069  default:
2070  Assert(false, ExcNotImplemented());
2071  }
2072 }
2073 
2074 
2075 
2076 template <int dim>
2077 QGaussWedge<dim>::QGaussWedge(const unsigned int n_points)
2078  : Quadrature<dim>()
2079 {
2080  AssertDimension(dim, 3);
2081 
2082  QGaussSimplex<2> quad_tri(n_points);
2083  QGauss<1> quad_line(n_points);
2084 
2085  for (unsigned int i = 0; i < quad_line.size(); ++i)
2086  for (unsigned int j = 0; j < quad_tri.size(); ++j)
2087  {
2088  this->quadrature_points.emplace_back(quad_tri.point(j)[0],
2089  quad_tri.point(j)[1],
2090  quad_line.point(i)[0]);
2091  this->weights.emplace_back(quad_tri.weight(j) * quad_line.weight(i));
2092  }
2093 
2094  AssertDimension(this->quadrature_points.size(), this->weights.size());
2095  Assert(this->quadrature_points.size() > 0,
2096  ExcMessage("No valid quadrature points!"));
2097 }
2098 
2099 
2100 
2101 template <int dim>
2102 QGaussPyramid<dim>::QGaussPyramid(const unsigned int n_points_1D)
2103  : Quadrature<dim>()
2104 {
2105  AssertDimension(dim, 3);
2106 
2107  if (n_points_1D == 1)
2108  {
2109  const double Q14 = 1.0 / 4.0;
2110  const double Q43 = 4.0 / 3.0;
2111 
2112  this->quadrature_points.emplace_back(0, 0, Q14);
2113  this->weights.emplace_back(Q43);
2114  }
2115  else if (n_points_1D == 2)
2116  {
2117  // clang-format off
2118  this->quadrature_points.emplace_back(-0.26318405556971, -0.26318405556971, 0.54415184401122);
2119  this->quadrature_points.emplace_back(-0.50661630334979, -0.50661630334979, 0.12251482265544);
2120  this->quadrature_points.emplace_back(-0.26318405556971, +0.26318405556971, 0.54415184401122);
2121  this->quadrature_points.emplace_back(-0.50661630334979, +0.50661630334979, 0.12251482265544);
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  // clang-format on
2127 
2128  this->weights.emplace_back(0.10078588207983);
2129  this->weights.emplace_back(0.23254745125351);
2130  this->weights.emplace_back(0.10078588207983);
2131  this->weights.emplace_back(0.23254745125351);
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  }
2137 
2138  AssertDimension(this->quadrature_points.size(), this->weights.size());
2139  Assert(this->quadrature_points.size() > 0,
2140  ExcMessage("No valid quadrature points!"));
2141 }
2142 
2143 
2144 
2145 // explicit specialization
2146 // note that 1d formulae are specialized by implementation above
2147 template class QGauss<2>;
2148 template class QGaussLobatto<2>;
2149 template class QMidpoint<2>;
2150 template class QTrapezoid<2>;
2151 template class QSimpson<2>;
2152 template class QMilne<2>;
2153 template class QWeddle<2>;
2154 
2155 template class QGauss<3>;
2156 template class QGaussLobatto<3>;
2157 template class QMidpoint<3>;
2158 template class QTrapezoid<3>;
2159 template class QSimpson<3>;
2160 template class QMilne<3>;
2161 template class QWeddle<3>;
2162 
2163 template class QSorted<1>;
2164 template class QSorted<2>;
2165 template class QSorted<3>;
2166 
2167 template class QTelles<1>;
2168 template class QTelles<2>;
2169 template class QTelles<3>;
2170 
2171 template class QGaussChebyshev<1>;
2172 template class QGaussChebyshev<2>;
2173 template class QGaussChebyshev<3>;
2174 
2175 template class QGaussRadauChebyshev<1>;
2176 template class QGaussRadauChebyshev<2>;
2177 template class QGaussRadauChebyshev<3>;
2178 
2179 template class QGaussLobattoChebyshev<1>;
2180 template class QGaussLobattoChebyshev<2>;
2181 template class QGaussLobattoChebyshev<3>;
2182 
2183 template class QSimplex<1>;
2184 template class QSimplex<2>;
2185 template class QSimplex<3>;
2186 
2187 template class QIteratedSimplex<1>;
2188 template class QIteratedSimplex<2>;
2189 template class QIteratedSimplex<3>;
2190 
2191 template class QSplit<1>;
2192 template class QSplit<2>;
2193 template class QSplit<3>;
2194 
2195 template class QGaussSimplex<0>;
2196 template class QGaussSimplex<1>;
2197 template class QGaussSimplex<2>;
2198 template class QGaussSimplex<3>;
2199 template class QGaussWedge<1>;
2200 template class QGaussWedge<2>;
2201 template class QGaussWedge<3>;
2202 template class QGaussPyramid<1>;
2203 template class QGaussPyramid<2>;
2204 template class QGaussPyramid<3>;
2205 
2206 template class QWitherdenVincentSimplex<1>;
2207 template class QWitherdenVincentSimplex<2>;
2208 template class QWitherdenVincentSimplex<3>;
2209 
Transformed quadrature weights.
QGaussLogR(const unsigned int n, const Point< dim > &x0=Point< dim >(), const double alpha=1, const bool factor_out_singular_weight=false)
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &)
Number * begin_raw()
QGaussSimplex(const unsigned int n_points_1D)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1655
QGaussLog(const unsigned int n, const bool revert=false)
std::vector< double > weights
Definition: quadrature.h:295
QGaussOneOverR(const unsigned int n, const Point< dim > &singularity, const bool factor_out_singular_weight=false)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
QGaussWedge(const unsigned int n_points_1D)
IteratorRange< active_cell_iterator > active_cell_iterators() const
Definition: tria.cc:13264
const ::Triangulation< dim, spacedim > & tria
#define AssertIndexRange(index, range)
Definition: exceptions.h:1720
const double fraction
QGaussChebyshev(const unsigned int n)
Generate a formula with n quadrature points.
const Point< dim > & point(const unsigned int i) const
std::vector< double > get_quadrature_weights(const unsigned int n)
static std::vector< double > get_quadrature_points(const unsigned int n)
QSimplex(const Quadrature< dim > &quad)
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:472
Transformed quadrature points.
static Point< dim > unit_cell_vertex(const unsigned int vertex)
std::vector< double > get_quadrature_points(const unsigned int n)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1571
QGauss(const unsigned int n)
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:482
Quadrature & operator=(const Quadrature< dim > &)
Definition: quadrature.cc:281
void split_point(const Point< dim1+dim2 > &source, Point< dim1 > &p1, Point< dim2 > &p2)
Number jacobi_polynomial_value(const unsigned int degree, const int alpha, const int beta, const Number x)
Definition: polynomial.h:976
QIteratedSimplex(const Quadrature< dim > &base_quadrature, const unsigned int n_copies)
const Point< spacedim > & quadrature_point(const unsigned int q) const
Number * end_raw()
QGaussPyramid(const unsigned int n_points_1D)
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1461
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
static std::vector< double > get_quadrature_weights(const unsigned int n)
QGaussRadauChebyshev(const unsigned int n, EndPoint ep=QGaussRadauChebyshev::left)
Generate a formula with n quadrature points.
QDuffy(const Quadrature< 1 > &radial_quadrature, const Quadrature< 1 > &angular_quadrature, const double beta=1.0)
std::vector< double > get_quadrature_weights(const unsigned int n)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:401
bool compare_weights(const unsigned int a, const unsigned int b) const
Point< 3 > vertices[4]
QGaussLobatto(const unsigned int n)
static constexpr double PI_2
Definition: numbers.h:236
QWitherdenVincentSimplex(const unsigned int n_points_1D, const bool use_odd_order=true)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
bool is_tensor_product_flag
Definition: quadrature.h:304
std::vector< Point< dim > > quadrature_points
Definition: quadrature.h:289
unsigned int size() const
static unsigned int quad_size(const Point< dim > &singularity, const unsigned int n)
QTelles(const Quadrature< 1 > &base_quad, const Point< dim > &singularity)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
void reference_cell(Triangulation< dim, spacedim > &tria, const ReferenceCell &reference_cell)
std::vector< long double > compute_quadrature_weights(const std::vector< long double > &x, const int alpha, const int beta)
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
std::vector< double > get_quadrature_points(const unsigned int n, ::QGaussRadauChebyshev< 1 >::EndPoint ep)
Definition: tensor.h:506
static constexpr double PI
Definition: numbers.h:231
double JxW(const unsigned int quadrature_point) const
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:400
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:138
QGaussLobattoChebyshev(const unsigned int n)
Generate a formula with n quadrature points.
void refine_global(const unsigned int times=1)
Definition: tria.cc:11890
static ::ExceptionBase & ExcNotImplemented()
QSorted(const Quadrature< dim > &quad)
long double gamma(const unsigned int n)
QSplit(const QSimplex< dim > &base, const Point< dim > &split_point)
std::vector< double > get_quadrature_weights(const unsigned int n, ::QGaussRadauChebyshev< 1 >::EndPoint ep)
void copy(const T *begin, const T *end, U *dest)
Quadrature< dim > compute_affine_transformation(const std::array< Point< dim >, dim+1 > &vertices) const
std::vector< double > get_quadrature_points(const unsigned int n)
double weight(const unsigned int i) const
static ::ExceptionBase & ExcInternalError()
QTrianglePolar(const Quadrature< 1 > &radial_quadrature, const Quadrature< 1 > &angular_quadrature)
void reinit(const TriaIterator< DoFCellAccessor< dim, spacedim, level_dof_access >> &cell)
Definition: fe_values.cc:4515