Reference documentation for deal.II version Git 8236ec9b21 2019-11-20 14:24:58 -0500
\(\newcommand{\dealcoloneq}{\mathrel{\vcenter{:}}=}\)
quadrature_lib.cc
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2019 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 
16 #include <deal.II/base/geometry_info.h>
17 #include <deal.II/base/polynomial.h>
18 #include <deal.II/base/quadrature_lib.h>
19 
20 #include <algorithm>
21 #include <cmath>
22 #include <functional>
23 #include <limits>
24 
25 
26 DEAL_II_NAMESPACE_OPEN
27 
28 
29 // please note: for a given dimension, we need the quadrature formulae
30 // for all lower dimensions as well. That is why in this file the check
31 // is for deal_II_dimension >= any_number and not for ==
32 
33 
34 
35 template <>
36 QGauss<0>::QGauss(const unsigned int)
37  : // there are n_q^dim == 1
38  // points
39  Quadrature<0>(1)
40 {
41  // the single quadrature point gets unit
42  // weight
43  this->weights[0] = 1;
44 }
45 
46 
47 
48 template <>
49 QGaussLobatto<0>::QGaussLobatto(const unsigned int)
50  : // there are n_q^dim == 1
51  // points
52  Quadrature<0>(1)
53 {
54  // the single quadrature point gets unit
55  // weight
56  this->weights[0] = 1;
57 }
58 
59 
60 
61 template <>
62 QGauss<1>::QGauss(const unsigned int n)
63  : Quadrature<1>(n)
64 {
65  if (n == 0)
66  return;
67 
68  std::vector<long double> points =
69  Polynomials::jacobi_polynomial_roots<long double>(n, 0, 0);
70 
71  for (unsigned int i = 0; i < (points.size() + 1) / 2; ++i)
72  {
73  this->quadrature_points[i][0] = points[i];
74  this->quadrature_points[n - i - 1][0] = 1. - points[i];
75 
76  // derivative of Jacobi polynomial
77  const long double pp =
78  0.5 * (n + 1) *
79  Polynomials::jacobi_polynomial_value(n - 1, 1, 1, points[i]);
80  const long double x = -1. + 2. * points[i];
81  const double w = 1. / ((1. - x * x) * pp * pp);
82  this->weights[i] = w;
83  this->weights[n - i - 1] = w;
84  }
85 }
86 
87 namespace internal
88 {
89  namespace QGaussLobatto
90  {
95  long double
96  gamma(const unsigned int n)
97  {
98  long double result = n - 1;
99  for (int i = n - 2; i > 1; --i)
100  result *= i;
101  return result;
102  }
103 
104 
105 
113  std::vector<long double>
114  compute_quadrature_weights(const std::vector<long double> &x,
115  const int alpha,
116  const int beta)
117  {
118  const unsigned int q = x.size();
119  std::vector<long double> w(q);
120 
121  const long double factor =
122  std::pow(2., alpha + beta + 1) * gamma(alpha + q) * gamma(beta + q) /
123  ((q - 1) * gamma(q) * gamma(alpha + beta + q + 1));
124  for (unsigned int i = 0; i < q; ++i)
125  {
126  const long double s =
127  Polynomials::jacobi_polynomial_value(q - 1, alpha, beta, x[i]);
128  w[i] = factor / (s * s);
129  }
130  w[0] *= (beta + 1);
131  w[q - 1] *= (alpha + 1);
132 
133  return w;
134  }
135  } // namespace QGaussLobatto
136 } // namespace internal
137 
138 
139 
140 template <>
141 QGaussLobatto<1>::QGaussLobatto(const unsigned int n)
142  : Quadrature<1>(n)
143 {
144  Assert(n >= 2, ExcNotImplemented());
145 
146  std::vector<long double> points =
147  Polynomials::jacobi_polynomial_roots<long double>(n - 2, 1, 1);
148  points.insert(points.begin(), 0);
149  points.push_back(1.);
150  std::vector<long double> w =
151  internal::QGaussLobatto::compute_quadrature_weights(points, 0, 0);
152 
153  // scale weights to the interval [0.0, 1.0]:
154  for (unsigned int i = 0; i < points.size(); ++i)
155  {
156  this->quadrature_points[i][0] = points[i];
157  this->weights[i] = 0.5 * w[i];
158  }
159 }
160 
161 
162 
163 template <>
165  : Quadrature<1>(1)
166 {
167  this->quadrature_points[0] = Point<1>(0.5);
168  this->weights[0] = 1.0;
169 }
170 
171 
172 
173 template <>
175  : Quadrature<1>(2)
176 {
177  static const double xpts[] = {0.0, 1.0};
178  static const double wts[] = {0.5, 0.5};
179 
180  for (unsigned int i = 0; i < this->size(); ++i)
181  {
182  this->quadrature_points[i] = Point<1>(xpts[i]);
183  this->weights[i] = wts[i];
184  }
185 }
186 
187 
188 
189 template <>
191  : Quadrature<1>(3)
192 {
193  static const double xpts[] = {0.0, 0.5, 1.0};
194  static const double wts[] = {1. / 6., 2. / 3., 1. / 6.};
195 
196  for (unsigned int i = 0; i < this->size(); ++i)
197  {
198  this->quadrature_points[i] = Point<1>(xpts[i]);
199  this->weights[i] = wts[i];
200  }
201 }
202 
203 
204 
205 template <>
207  : Quadrature<1>(5)
208 {
209  static const double xpts[] = {0.0, .25, .5, .75, 1.0};
210  static const double wts[] = {
211  7. / 90., 32. / 90., 12. / 90., 32. / 90., 7. / 90.};
212 
213  for (unsigned int i = 0; i < this->size(); ++i)
214  {
215  this->quadrature_points[i] = Point<1>(xpts[i]);
216  this->weights[i] = wts[i];
217  }
218 }
219 
220 
221 
222 template <>
224  : Quadrature<1>(7)
225 {
226  static const double xpts[] = {
227  0.0, 1. / 6., 1. / 3., .5, 2. / 3., 5. / 6., 1.0};
228  static const double wts[] = {41. / 840.,
229  216. / 840.,
230  27. / 840.,
231  272. / 840.,
232  27. / 840.,
233  216. / 840.,
234  41. / 840.};
235 
236  for (unsigned int i = 0; i < this->size(); ++i)
237  {
238  this->quadrature_points[i] = Point<1>(xpts[i]);
239  this->weights[i] = wts[i];
240  }
241 }
242 
243 
244 template <>
245 QGaussLog<1>::QGaussLog(const unsigned int n, const bool revert)
246  : Quadrature<1>(n)
247 {
248  std::vector<double> p = get_quadrature_points(n);
249  std::vector<double> w = get_quadrature_weights(n);
250 
251  for (unsigned int i = 0; i < this->size(); ++i)
252  {
253  // Using the change of variables x=1-t, it's possible to show
254  // that int f(x)ln|1-x| = int f(1-t) ln|t|, which implies that
255  // we can use this quadrature formula also with weight ln|1-x|.
256  this->quadrature_points[i] =
257  revert ? Point<1>(1 - p[n - 1 - i]) : Point<1>(p[i]);
258  this->weights[i] = revert ? w[n - 1 - i] : w[i];
259  }
260 }
261 
262 template <>
263 std::vector<double>
264 QGaussLog<1>::get_quadrature_points(const unsigned int n)
265 {
266  std::vector<double> q_points(n);
267 
268  switch (n)
269  {
270  case 1:
271  q_points[0] = 0.3333333333333333;
272  break;
273 
274  case 2:
275  q_points[0] = 0.1120088061669761;
276  q_points[1] = 0.6022769081187381;
277  break;
278 
279  case 3:
280  q_points[0] = 0.06389079308732544;
281  q_points[1] = 0.3689970637156184;
282  q_points[2] = 0.766880303938942;
283  break;
284 
285  case 4:
286  q_points[0] = 0.04144848019938324;
287  q_points[1] = 0.2452749143206022;
288  q_points[2] = 0.5561654535602751;
289  q_points[3] = 0.848982394532986;
290  break;
291 
292  case 5:
293  q_points[0] = 0.02913447215197205;
294  q_points[1] = 0.1739772133208974;
295  q_points[2] = 0.4117025202849029;
296  q_points[3] = 0.6773141745828183;
297  q_points[4] = 0.89477136103101;
298  break;
299 
300  case 6:
301  q_points[0] = 0.02163400584411693;
302  q_points[1] = 0.1295833911549506;
303  q_points[2] = 0.3140204499147661;
304  q_points[3] = 0.5386572173517997;
305  q_points[4] = 0.7569153373774084;
306  q_points[5] = 0.922668851372116;
307  break;
308 
309 
310  case 7:
311  q_points[0] = 0.0167193554082585;
312  q_points[1] = 0.100185677915675;
313  q_points[2] = 0.2462942462079286;
314  q_points[3] = 0.4334634932570557;
315  q_points[4] = 0.6323509880476823;
316  q_points[5] = 0.81111862674023;
317  q_points[6] = 0.940848166743287;
318  break;
319 
320  case 8:
321  q_points[0] = 0.01332024416089244;
322  q_points[1] = 0.07975042901389491;
323  q_points[2] = 0.1978710293261864;
324  q_points[3] = 0.354153994351925;
325  q_points[4] = 0.5294585752348643;
326  q_points[5] = 0.7018145299391673;
327  q_points[6] = 0.849379320441094;
328  q_points[7] = 0.953326450056343;
329  break;
330 
331  case 9:
332  q_points[0] = 0.01086933608417545;
333  q_points[1] = 0.06498366633800794;
334  q_points[2] = 0.1622293980238825;
335  q_points[3] = 0.2937499039716641;
336  q_points[4] = 0.4466318819056009;
337  q_points[5] = 0.6054816627755208;
338  q_points[6] = 0.7541101371585467;
339  q_points[7] = 0.877265828834263;
340  q_points[8] = 0.96225055941096;
341  break;
342 
343  case 10:
344  q_points[0] = 0.00904263096219963;
345  q_points[1] = 0.05397126622250072;
346  q_points[2] = 0.1353118246392511;
347  q_points[3] = 0.2470524162871565;
348  q_points[4] = 0.3802125396092744;
349  q_points[5] = 0.5237923179723384;
350  q_points[6] = 0.6657752055148032;
351  q_points[7] = 0.7941904160147613;
352  q_points[8] = 0.898161091216429;
353  q_points[9] = 0.9688479887196;
354  break;
355 
356 
357  case 11:
358  q_points[0] = 0.007643941174637681;
359  q_points[1] = 0.04554182825657903;
360  q_points[2] = 0.1145222974551244;
361  q_points[3] = 0.2103785812270227;
362  q_points[4] = 0.3266955532217897;
363  q_points[5] = 0.4554532469286375;
364  q_points[6] = 0.5876483563573721;
365  q_points[7] = 0.7139638500230458;
366  q_points[8] = 0.825453217777127;
367  q_points[9] = 0.914193921640008;
368  q_points[10] = 0.973860256264123;
369  break;
370 
371  case 12:
372  q_points[0] = 0.006548722279080035;
373  q_points[1] = 0.03894680956045022;
374  q_points[2] = 0.0981502631060046;
375  q_points[3] = 0.1811385815906331;
376  q_points[4] = 0.2832200676673157;
377  q_points[5] = 0.398434435164983;
378  q_points[6] = 0.5199526267791299;
379  q_points[7] = 0.6405109167754819;
380  q_points[8] = 0.7528650118926111;
381  q_points[9] = 0.850240024421055;
382  q_points[10] = 0.926749682988251;
383  q_points[11] = 0.977756129778486;
384  break;
385 
386  default:
387  Assert(false, ExcNotImplemented());
388  break;
389  }
390 
391  return q_points;
392 }
393 
394 
395 template <>
396 std::vector<double>
397 QGaussLog<1>::get_quadrature_weights(const unsigned int n)
398 {
399  std::vector<double> quadrature_weights(n);
400 
401  switch (n)
402  {
403  case 1:
404  quadrature_weights[0] = -1.0;
405  break;
406  case 2:
407  quadrature_weights[0] = -0.7185393190303845;
408  quadrature_weights[1] = -0.2814606809696154;
409  break;
410 
411  case 3:
412  quadrature_weights[0] = -0.5134045522323634;
413  quadrature_weights[1] = -0.3919800412014877;
414  quadrature_weights[2] = -0.0946154065661483;
415  break;
416 
417  case 4:
418  quadrature_weights[0] = -0.3834640681451353;
419  quadrature_weights[1] = -0.3868753177747627;
420  quadrature_weights[2] = -0.1904351269501432;
421  quadrature_weights[3] = -0.03922548712995894;
422  break;
423 
424  case 5:
425  quadrature_weights[0] = -0.2978934717828955;
426  quadrature_weights[1] = -0.3497762265132236;
427  quadrature_weights[2] = -0.234488290044052;
428  quadrature_weights[3] = -0.0989304595166356;
429  quadrature_weights[4] = -0.01891155214319462;
430  break;
431 
432  case 6:
433  quadrature_weights[0] = -0.2387636625785478;
434  quadrature_weights[1] = -0.3082865732739458;
435  quadrature_weights[2] = -0.2453174265632108;
436  quadrature_weights[3] = -0.1420087565664786;
437  quadrature_weights[4] = -0.05545462232488041;
438  quadrature_weights[5] = -0.01016895869293513;
439  break;
440 
441 
442  case 7:
443  quadrature_weights[0] = -0.1961693894252476;
444  quadrature_weights[1] = -0.2703026442472726;
445  quadrature_weights[2] = -0.239681873007687;
446  quadrature_weights[3] = -0.1657757748104267;
447  quadrature_weights[4] = -0.0889432271377365;
448  quadrature_weights[5] = -0.03319430435645653;
449  quadrature_weights[6] = -0.005932787015162054;
450  break;
451 
452  case 8:
453  quadrature_weights[0] = -0.164416604728002;
454  quadrature_weights[1] = -0.2375256100233057;
455  quadrature_weights[2] = -0.2268419844319134;
456  quadrature_weights[3] = -0.1757540790060772;
457  quadrature_weights[4] = -0.1129240302467932;
458  quadrature_weights[5] = -0.05787221071771947;
459  quadrature_weights[6] = -0.02097907374214317;
460  quadrature_weights[7] = -0.003686407104036044;
461  break;
462 
463  case 9:
464  quadrature_weights[0] = -0.1400684387481339;
465  quadrature_weights[1] = -0.2097722052010308;
466  quadrature_weights[2] = -0.211427149896601;
467  quadrature_weights[3] = -0.1771562339380667;
468  quadrature_weights[4] = -0.1277992280331758;
469  quadrature_weights[5] = -0.07847890261203835;
470  quadrature_weights[6] = -0.0390225049841783;
471  quadrature_weights[7] = -0.01386729555074604;
472  quadrature_weights[8] = -0.002408041036090773;
473  break;
474 
475  case 10:
476  quadrature_weights[0] = -0.12095513195457;
477  quadrature_weights[1] = -0.1863635425640733;
478  quadrature_weights[2] = -0.1956608732777627;
479  quadrature_weights[3] = -0.1735771421828997;
480  quadrature_weights[4] = -0.135695672995467;
481  quadrature_weights[5] = -0.0936467585378491;
482  quadrature_weights[6] = -0.05578772735275126;
483  quadrature_weights[7] = -0.02715981089692378;
484  quadrature_weights[8] = -0.00951518260454442;
485  quadrature_weights[9] = -0.001638157633217673;
486  break;
487 
488 
489  case 11:
490  quadrature_weights[0] = -0.1056522560990997;
491  quadrature_weights[1] = -0.1665716806006314;
492  quadrature_weights[2] = -0.1805632182877528;
493  quadrature_weights[3] = -0.1672787367737502;
494  quadrature_weights[4] = -0.1386970574017174;
495  quadrature_weights[5] = -0.1038334333650771;
496  quadrature_weights[6] = -0.06953669788988512;
497  quadrature_weights[7] = -0.04054160079499477;
498  quadrature_weights[8] = -0.01943540249522013;
499  quadrature_weights[9] = -0.006737429326043388;
500  quadrature_weights[10] = -0.001152486965101561;
501  break;
502 
503  case 12:
504  quadrature_weights[0] = -0.09319269144393;
505  quadrature_weights[1] = -0.1497518275763289;
506  quadrature_weights[2] = -0.166557454364573;
507  quadrature_weights[3] = -0.1596335594369941;
508  quadrature_weights[4] = -0.1384248318647479;
509  quadrature_weights[5] = -0.1100165706360573;
510  quadrature_weights[6] = -0.07996182177673273;
511  quadrature_weights[7] = -0.0524069547809709;
512  quadrature_weights[8] = -0.03007108900074863;
513  quadrature_weights[9] = -0.01424924540252916;
514  quadrature_weights[10] = -0.004899924710875609;
515  quadrature_weights[11] = -0.000834029009809656;
516  break;
517 
518  default:
519  Assert(false, ExcNotImplemented());
520  break;
521  }
522 
523  return quadrature_weights;
524 }
525 
526 
527 template <>
528 QGaussLogR<1>::QGaussLogR(const unsigned int n,
529  const Point<1> origin,
530  const double alpha,
531  const bool factor_out_singularity)
532  : Quadrature<1>(
533  ((origin[0] == 0) || (origin[0] == 1)) ? (alpha == 1 ? n : 2 * n) : 4 * n)
534  , fraction(((origin[0] == 0) || (origin[0] == 1.)) ? 1. : origin[0])
535 {
536  // The three quadrature formulas that make this one up. There are
537  // at most two when the origin is one of the extremes, and there is
538  // only one if the origin is one of the extremes and alpha is
539  // equal to one.
540  //
541  // If alpha is different from one, then we need a correction which
542  // is performed with a standard Gauss quadrature rule on each
543  // segment. This is not needed in the standard case where alpha is
544  // equal to one and the origin is on one of the extremes. We
545  // integrate with weight ln|(x-o)/alpha|. In the easy cases, we
546  // only need n quadrature points. In the most difficult one, we
547  // need 2*n points for the first segment, and 2*n points for the
548  // second segment.
549  QGaussLog<1> quad1(n, origin[0] != 0);
550  QGaussLog<1> quad2(n);
551  QGauss<1> quad(n);
552 
553  // Check that the origin is inside 0,1
554  Assert((fraction >= 0) && (fraction <= 1),
555  ExcMessage("Origin is outside [0,1]."));
556 
557  // Non singular offset. This is the start of non singular quad
558  // points.
559  unsigned int ns_offset = (fraction == 1) ? n : 2 * n;
560 
561  for (unsigned int i = 0, j = ns_offset; i < n; ++i, ++j)
562  {
563  // The first i quadrature points are the same as quad1, and
564  // are by default singular.
565  this->quadrature_points[i] = quad1.point(i) * fraction;
566  this->weights[i] = quad1.weight(i) * fraction;
567 
568  // We need to scale with -log|fraction*alpha|
569  if ((alpha != 1) || (fraction != 1))
570  {
571  this->quadrature_points[j] = quad.point(i) * fraction;
572  this->weights[j] =
573  -std::log(alpha / fraction) * quad.weight(i) * fraction;
574  }
575  // In case we need the second quadrature as well, do it now.
576  if (fraction != 1)
577  {
578  this->quadrature_points[i + n] =
579  quad2.point(i) * (1 - fraction) + Point<1>(fraction);
580  this->weights[i + n] = quad2.weight(i) * (1 - fraction);
581 
582  // We need to scale with -log|fraction*alpha|
583  this->quadrature_points[j + n] =
584  quad.point(i) * (1 - fraction) + Point<1>(fraction);
585  this->weights[j + n] =
586  -std::log(alpha / (1 - fraction)) * quad.weight(i) * (1 - fraction);
587  }
588  }
589  if (factor_out_singularity == true)
590  for (unsigned int i = 0; i < size(); ++i)
591  {
592  Assert(
593  this->quadrature_points[i] != origin,
594  ExcMessage(
595  "The singularity cannot be on a Gauss point of the same order!"));
596  double denominator =
597  std::log(std::abs((this->quadrature_points[i] - origin)[0]) / alpha);
598  Assert(denominator != 0.0,
599  ExcMessage(
600  "The quadrature formula you are using does not allow to "
601  "factor out the singularity, which is zero at one point."));
602  this->weights[i] /= denominator;
603  }
604 }
605 
606 
607 template <>
608 unsigned int
609 QGaussOneOverR<2>::quad_size(const Point<2> singularity, const unsigned int n)
610 {
611  double eps = 1e-8;
612  bool on_edge = false;
613  bool on_vertex = false;
614  for (unsigned int i = 0; i < 2; ++i)
615  if ((std::abs(singularity[i]) < eps) ||
616  (std::abs(singularity[i] - 1) < eps))
617  on_edge = true;
618  if (on_edge &&
619  (std::abs((singularity - Point<2>(.5, .5)).norm_square() - .5) < eps))
620  on_vertex = true;
621  if (on_vertex)
622  return (2 * n * n);
623  if (on_edge)
624  return (4 * n * n);
625  return (8 * n * n);
626 }
627 
628 template <>
629 QGaussOneOverR<2>::QGaussOneOverR(const unsigned int n,
630  const Point<2> singularity,
631  const bool factor_out_singularity)
632  : Quadrature<2>(quad_size(singularity, n))
633 {
634  // We treat all the cases in the
635  // same way. Split the element in 4
636  // pieces, measure the area, if
637  // it's relevant, add the
638  // quadrature connected to that
639  // singularity.
640  std::vector<QGaussOneOverR<2>> quads;
641  std::vector<Point<2>> origins;
642  // Id of the corner with a
643  // singularity
644  quads.emplace_back(n, 3, factor_out_singularity);
645  quads.emplace_back(n, 2, factor_out_singularity);
646  quads.emplace_back(n, 1, factor_out_singularity);
647  quads.emplace_back(n, 0, factor_out_singularity);
648 
649  origins.emplace_back(0., 0.);
650  origins.emplace_back(singularity[0], 0.);
651  origins.emplace_back(0., singularity[1]);
652  origins.push_back(singularity);
653 
654  // Lexicographical ordering.
655 
656  double eps = 1e-8;
657  unsigned int q_id = 0; // Current quad point index.
658  Tensor<1, 2> dist;
659 
660  for (unsigned int box = 0; box < 4; ++box)
661  {
662  dist = (singularity - GeometryInfo<2>::unit_cell_vertex(box));
663  dist = Point<2>(std::abs(dist[0]), std::abs(dist[1]));
664  double area = dist[0] * dist[1];
665  if (area > eps)
666  for (unsigned int q = 0; q < quads[box].size(); ++q, ++q_id)
667  {
668  const Point<2> &qp = quads[box].point(q);
669  this->quadrature_points[q_id] =
670  origins[box] + Point<2>(dist[0] * qp[0], dist[1] * qp[1]);
671  this->weights[q_id] = quads[box].weight(q) * area;
672  }
673  }
674 }
675 
676 
677 template <>
678 QGaussOneOverR<2>::QGaussOneOverR(const unsigned int n,
679  const unsigned int vertex_index,
680  const bool factor_out_singularity)
681  : Quadrature<2>(2 * n * n)
682 {
683  // This version of the constructor
684  // works only for the 4
685  // vertices. If you need a more
686  // general one, you should use the
687  // one with the Point<2> in the
688  // constructor.
689  Assert(vertex_index < 4, ExcIndexRange(vertex_index, 0, 4));
690 
691  // Start with the gauss quadrature formula on the (u,v) reference
692  // element.
693  QGauss<2> gauss(n);
694 
695  Assert(gauss.size() == n * n, ExcInternalError());
696 
697  // For the moment we only implemented this for the vertices of a
698  // quadrilateral. We are planning to do this also for the support
699  // points of arbitrary FE_Q elements, to allow the use of this
700  // class in boundary element programs with higher order mappings.
701  Assert(vertex_index < 4, ExcIndexRange(vertex_index, 0, 4));
702 
703  // We create only the first one. All other pieces are rotation of
704  // this one.
705  // In this case the transformation is
706  //
707  // (x,y) = (u, u tan(pi/4 v))
708  //
709  // with Jacobian
710  //
711  // J = pi/4 R / cos(pi/4 v)
712  //
713  // And we get rid of R to take into account the singularity,
714  // unless specified differently in the constructor.
715  std::vector<Point<2>> &ps = this->quadrature_points;
716  std::vector<double> & ws = this->weights;
717  double pi4 = numbers::PI / 4;
718 
719  for (unsigned int q = 0; q < gauss.size(); ++q)
720  {
721  const Point<2> &gp = gauss.point(q);
722  ps[q][0] = gp[0];
723  ps[q][1] = gp[0] * std::tan(pi4 * gp[1]);
724  ws[q] = gauss.weight(q) * pi4 / std::cos(pi4 * gp[1]);
725  if (factor_out_singularity)
726  ws[q] *= (ps[q] - GeometryInfo<2>::unit_cell_vertex(0)).norm();
727  // The other half of the quadrilateral is symmetric with
728  // respect to xy plane.
729  ws[gauss.size() + q] = ws[q];
730  ps[gauss.size() + q][0] = ps[q][1];
731  ps[gauss.size() + q][1] = ps[q][0];
732  }
733 
734  // Now we distribute these vertices in the correct manner
735  double theta = 0;
736  switch (vertex_index)
737  {
738  case 0:
739  theta = 0;
740  break;
741  case 1:
742  //
743  theta = numbers::PI / 2;
744  break;
745  case 2:
746  theta = -numbers::PI / 2;
747  break;
748  case 3:
749  theta = numbers::PI;
750  break;
751  }
752 
753  double R00 = std::cos(theta), R01 = -std::sin(theta);
754  double R10 = std::sin(theta), R11 = std::cos(theta);
755 
756  if (vertex_index != 0)
757  for (unsigned int q = 0; q < size(); ++q)
758  {
759  double x = ps[q][0] - .5, y = ps[q][1] - .5;
760 
761  ps[q][0] = R00 * x + R01 * y + .5;
762  ps[q][1] = R10 * x + R11 * y + .5;
763  }
764 }
765 
766 
767 template <int dim>
769  : Quadrature<dim>(quad)
770 {
771  std::vector<unsigned int> permutation(quad.size());
772  for (unsigned int i = 0; i < quad.size(); ++i)
773  permutation[i] = i;
774 
775  std::sort(permutation.begin(),
776  permutation.end(),
777  [this](const unsigned int x, const unsigned int y) {
778  return this->compare_weights(x, y);
779  });
780 
781  // At this point, the variable is_tensor_product_flag is set
782  // to the respective value of the given Quadrature in the base
783  // class copy constructor.
784  // We only call a quadrature formula 'tensor product'
785  // if the quadrature points are also sorted lexicographically.
786  // In particular, any reordering destroys that property
787  // and we might need to modify the variable accordingly.
788  for (unsigned int i = 0; i < quad.size(); ++i)
789  {
790  this->weights[i] = quad.weight(permutation[i]);
791  this->quadrature_points[i] = quad.point(permutation[i]);
792  if (permutation[i] != i)
793  this->is_tensor_product_flag = false;
794  }
795 }
796 
797 
798 template <int dim>
799 bool
800 QSorted<dim>::compare_weights(const unsigned int a, const unsigned int b) const
801 {
802  return (this->weights[a] < this->weights[b]);
803 }
804 
805 
806 // construct the quadrature formulae in higher dimensions by
807 // tensor product of lower dimensions
808 
809 template <int dim>
810 QGauss<dim>::QGauss(const unsigned int n)
811  : Quadrature<dim>(QGauss<dim - 1>(n), QGauss<1>(n))
812 {}
813 
814 
815 
816 template <int dim>
818  : Quadrature<dim>(QGaussLobatto<dim - 1>(n), QGaussLobatto<1>(n))
819 {}
820 
821 
822 
823 template <int dim>
825  : Quadrature<dim>(QMidpoint<dim - 1>(), QMidpoint<1>())
826 {}
827 
828 
829 
830 template <int dim>
832  : Quadrature<dim>(QTrapez<dim - 1>(), QTrapez<1>())
833 {}
834 
835 
836 
837 template <int dim>
839  : Quadrature<dim>(QSimpson<dim - 1>(), QSimpson<1>())
840 {}
841 
842 
843 
844 template <int dim>
846  : Quadrature<dim>(QMilne<dim - 1>(), QMilne<1>())
847 {}
848 
849 
850 template <int dim>
852  : Quadrature<dim>(QWeddle<dim - 1>(), QWeddle<1>())
853 {}
854 
855 template <int dim>
857  const Point<dim> & singularity)
858  : // We need the explicit implementation if dim == 1. If dim > 1 we use the
859  // former implementation and apply a tensorial product to obtain the higher
860  // dimensions.
861  Quadrature<dim>(
862  dim == 2 ?
863  QAnisotropic<dim>(QTelles<1>(base_quad, Point<1>(singularity[0])),
864  QTelles<1>(base_quad, Point<1>(singularity[1]))) :
865  dim == 3 ?
866  QAnisotropic<dim>(QTelles<1>(base_quad, Point<1>(singularity[0])),
867  QTelles<1>(base_quad, Point<1>(singularity[1])),
868  QTelles<1>(base_quad, Point<1>(singularity[2]))) :
869  Quadrature<dim>())
870 {}
871 
872 template <int dim>
873 QTelles<dim>::QTelles(const unsigned int n, const Point<dim> &singularity)
874  : // In this case we map the standard Gauss Legendre formula using the given
875  // singularity point coordinates.
876  Quadrature<dim>(QTelles<dim>(QGauss<1>(n), singularity))
877 {}
878 
879 
880 
881 template <>
882 QTelles<1>::QTelles(const Quadrature<1> &base_quad, const Point<1> &singularity)
883  : // We explicitly implement the Telles' variable change if dim == 1.
884  Quadrature<1>(base_quad)
885 {
886  // We define all the constants to be used in the implementation of
887  // Telles' rule
888  const double eta_bar = singularity[0] * 2. - 1.;
889  const double eta_star = eta_bar * eta_bar - 1.;
890  double gamma_bar;
891 
892  std::vector<Point<1>> quadrature_points_dummy(quadrature_points.size());
893  std::vector<double> weights_dummy(weights.size());
894  unsigned int cont = 0;
895  const double tol = 1e-10;
896  for (unsigned int d = 0; d < quadrature_points.size(); ++d)
897  {
898  if (std::abs(quadrature_points[d][0] - singularity[0]) > tol)
899  {
900  quadrature_points_dummy[d - cont] = quadrature_points[d];
901  weights_dummy[d - cont] = weights[d];
902  }
903  else
904  {
905  // We need to remove the singularity point from the quadrature point
906  // list. To do so we use the variable cont.
907  cont = 1;
908  }
909  }
910  if (cont == 1)
911  {
912  quadrature_points.resize(quadrature_points_dummy.size() - 1);
913  weights.resize(weights_dummy.size() - 1);
914  for (unsigned int d = 0; d < quadrature_points.size(); ++d)
915  {
916  quadrature_points[d] = quadrature_points_dummy[d];
917  weights[d] = weights_dummy[d];
918  }
919  }
920  // We need to check if the singularity is at the boundary of the interval.
921  if (std::abs(eta_star) <= tol)
922  {
923  gamma_bar =
924  std::pow((eta_bar * eta_star + std::abs(eta_star)), 1.0 / 3.0) +
925  std::pow((eta_bar * eta_star - std::abs(eta_star)), 1.0 / 3.0) +
926  eta_bar;
927  }
928  else
929  {
930  gamma_bar = (eta_bar * eta_star + std::abs(eta_star)) /
931  std::abs(eta_bar * eta_star + std::abs(eta_star)) *
932  std::pow(std::abs(eta_bar * eta_star + std::abs(eta_star)),
933  1.0 / 3.0) +
934  (eta_bar * eta_star - std::abs(eta_star)) /
935  std::abs(eta_bar * eta_star - std::abs(eta_star)) *
936  std::pow(std::abs(eta_bar * eta_star - std::abs(eta_star)),
937  1.0 / 3.0) +
938  eta_bar;
939  }
940  for (unsigned int q = 0; q < quadrature_points.size(); ++q)
941  {
942  double gamma = quadrature_points[q][0] * 2 - 1;
943  double eta = (std::pow(gamma - gamma_bar, 3.0) +
944  gamma_bar * (gamma_bar * gamma_bar + 3)) /
945  (1 + 3 * gamma_bar * gamma_bar);
946 
947  double J = 3 * ((gamma - gamma_bar) * (gamma - gamma_bar)) /
948  (1 + 3 * gamma_bar * gamma_bar);
949 
950  quadrature_points[q][0] = (eta + 1) / 2.0;
951  weights[q] = J * weights[q];
952  }
953 }
954 
955 namespace internal
956 {
957  namespace QGaussChebyshev
958  {
962  std::vector<double>
963  get_quadrature_points(const unsigned int n)
964  {
965  std::vector<double> points(n);
966  // n point quadrature: index from 0 to n-1
967  for (unsigned short i = 0; i < n; ++i)
968  // would be cos((2i+1)Pi/(2N+2))
969  // put + Pi so we start from the smallest point
970  // then map from [-1,1] to [0,1]
971  points[i] =
972  1. / 2. *
973  (1. + std::cos(numbers::PI *
974  (1. + double(2 * i + 1) / double(2 * (n - 1) + 2))));
975 
976  return points;
977  }
978 
979 
980 
984  std::vector<double>
985  get_quadrature_weights(const unsigned int n)
986  {
987  std::vector<double> weights(n);
988 
989  for (unsigned short i = 0; i < n; ++i)
990  {
991  // same weights as on [-1,1]
992  weights[i] = numbers::PI / double(n);
993  }
994 
995  return weights;
996  }
997  } // namespace QGaussChebyshev
998 } // namespace internal
999 
1000 
1001 template <>
1002 QGaussChebyshev<1>::QGaussChebyshev(const unsigned int n)
1003  : Quadrature<1>(n)
1004 {
1005  Assert(n > 0, ExcMessage("Need at least one point for the quadrature rule"));
1006  std::vector<double> p = internal::QGaussChebyshev::get_quadrature_points(n);
1007  std::vector<double> w = internal::QGaussChebyshev::get_quadrature_weights(n);
1008 
1009  for (unsigned int i = 0; i < this->size(); ++i)
1010  {
1011  this->quadrature_points[i] = Point<1>(p[i]);
1012  this->weights[i] = w[i];
1013  }
1014 }
1015 
1016 
1017 template <int dim>
1019  : Quadrature<dim>(QGaussChebyshev<1>(n))
1020 {}
1021 
1022 
1023 namespace internal
1024 {
1025  namespace QGaussRadauChebyshev
1026  {
1027  // Computes the points of the quadrature formula.
1028  std::vector<double>
1029  get_quadrature_points(const unsigned int n,
1031  {
1032  std::vector<double> points(n);
1033  // n point quadrature: index from 0 to n-1
1034  for (unsigned short i = 0; i < n; ++i)
1035  // would be -cos(2i Pi/(2N+1))
1036  // put + Pi so we start from the smallest point
1037  // then map from [-1,1] to [0,1]
1038  switch (ep)
1039  {
1040  case ::QGaussRadauChebyshev<1>::left:
1041  {
1042  points[i] =
1043  1. / 2. *
1044  (1. -
1045  std::cos(numbers::PI *
1046  (1 + 2 * double(i) / (2 * double(n - 1) + 1.))));
1047  break;
1048  }
1049 
1050  case ::QGaussRadauChebyshev<1>::right:
1051  {
1052  points[i] =
1053  1. / 2. *
1054  (1. - std::cos(numbers::PI * (2 * double(n - 1 - i) /
1055  (2 * double(n - 1) + 1.))));
1056  break;
1057  }
1058 
1059  default:
1060  Assert(
1061  false,
1062  ExcMessage(
1063  "This constructor can only be called with either "
1064  "QGaussRadauChebyshev::left or QGaussRadauChebyshev::right as "
1065  "second argument."));
1066  }
1067 
1068  return points;
1069  }
1070 
1071 
1072 
1073  // Computes the weights of the quadrature formula.
1074  std::vector<double>
1075  get_quadrature_weights(const unsigned int n,
1077  {
1078  std::vector<double> weights(n);
1079 
1080  for (unsigned short i = 0; i < n; ++i)
1081  {
1082  // same weights as on [-1,1]
1083  weights[i] = 2. * numbers::PI / double(2 * (n - 1) + 1.);
1084  if (ep == ::QGaussRadauChebyshev<1>::left && i == 0)
1085  weights[i] /= 2.;
1086  else if (ep == ::QGaussRadauChebyshev<1>::right &&
1087  i == (n - 1))
1088  weights[i] /= 2.;
1089  }
1090 
1091  return weights;
1092  }
1093  } // namespace QGaussRadauChebyshev
1094 } // namespace internal
1095 
1096 
1097 template <>
1098 QGaussRadauChebyshev<1>::QGaussRadauChebyshev(const unsigned int n, EndPoint ep)
1099  : Quadrature<1>(n)
1100  , ep(ep)
1101 {
1102  Assert(n > 0, ExcMessage("Need at least one point for quadrature rules"));
1103  std::vector<double> p =
1104  internal::QGaussRadauChebyshev::get_quadrature_points(n, ep);
1105  std::vector<double> w =
1106  internal::QGaussRadauChebyshev::get_quadrature_weights(n, ep);
1107 
1108  for (unsigned int i = 0; i < this->size(); ++i)
1109  {
1110  this->quadrature_points[i] = Point<1>(p[i]);
1111  this->weights[i] = w[i];
1112  }
1113 }
1114 
1115 
1116 template <int dim>
1118  EndPoint ep)
1119  : Quadrature<dim>(QGaussRadauChebyshev<1>(
1120  n,
1121  static_cast<QGaussRadauChebyshev<1>::EndPoint>(ep)))
1122  , ep(ep)
1123 {}
1124 
1125 
1126 
1127 namespace internal
1128 {
1129  namespace QGaussLobattoChebyshev
1130  {
1131  // Computes the points of the quadrature formula.
1132  std::vector<double>
1133  get_quadrature_points(const unsigned int n)
1134  {
1135  std::vector<double> points(n);
1136  // n point quadrature: index from 0 to n-1
1137  for (unsigned short i = 0; i < n; ++i)
1138  // would be cos(i Pi/N)
1139  // put + Pi so we start from the smallest point
1140  // then map from [-1,1] to [0,1]
1141  points[i] =
1142  1. / 2. *
1143  (1. + std::cos(numbers::PI * (1 + double(i) / double(n - 1))));
1144 
1145  return points;
1146  }
1147 
1148  // Computes the weights of the quadrature formula.
1149  std::vector<double>
1150  get_quadrature_weights(const unsigned int n)
1151  {
1152  std::vector<double> weights(n);
1153 
1154  for (unsigned short i = 0; i < n; ++i)
1155  {
1156  // same weights as on [-1,1]
1157  weights[i] = numbers::PI / double((n - 1));
1158  if (i == 0 || i == (n - 1))
1159  weights[i] /= 2.;
1160  }
1161 
1162  return weights;
1163  }
1164  } // namespace QGaussLobattoChebyshev
1165 } // namespace internal
1166 
1167 
1168 
1169 template <>
1171  : Quadrature<1>(n)
1172 {
1173  Assert(n > 1,
1174  ExcMessage(
1175  "Need at least two points for Gauss-Lobatto quadrature rule"));
1176  std::vector<double> p =
1177  internal::QGaussLobattoChebyshev::get_quadrature_points(n);
1178  std::vector<double> w =
1179  internal::QGaussLobattoChebyshev::get_quadrature_weights(n);
1180 
1181  for (unsigned int i = 0; i < this->size(); ++i)
1182  {
1183  this->quadrature_points[i] = Point<1>(p[i]);
1184  this->weights[i] = w[i];
1185  }
1186 }
1187 
1188 
1189 template <int dim>
1191  : Quadrature<dim>(QGaussLobattoChebyshev<1>(n))
1192 {}
1193 
1194 
1195 
1196 template <int dim>
1198 {
1199  std::vector<Point<dim>> qpoints;
1200  std::vector<double> weights;
1201 
1202  for (unsigned int i = 0; i < quad.size(); ++i)
1203  {
1204  double r = 0;
1205  for (unsigned int d = 0; d < dim; ++d)
1206  r += quad.point(i)[d];
1207  if (r <= 1 + 1e-10)
1208  {
1209  this->quadrature_points.push_back(quad.point(i));
1210  this->weights.push_back(quad.weight(i));
1211  }
1212  }
1213 }
1214 
1215 
1216 
1217 template <int dim>
1220  const std::array<Point<dim>, dim + 1> &vertices) const
1221 {
1222  Tensor<2, dim> B;
1223  for (unsigned int d = 0; d < dim; ++d)
1224  B[d] = vertices[d + 1] - vertices[0];
1225 
1226  B = transpose(B);
1227  const double J = std::abs(determinant(B));
1228 
1229  // if the determinant is zero, we return an empty quadrature
1230  if (J < 1e-12)
1231  return Quadrature<dim>();
1232 
1233  std::vector<Point<dim>> qp(this->size());
1234  std::vector<double> w(this->size());
1235 
1236  for (unsigned int i = 0; i < this->size(); ++i)
1237  {
1238  qp[i] = Point<dim>(vertices[0] + B * this->point(i));
1239  w[i] = J * this->weight(i);
1240  }
1241 
1242  return Quadrature<dim>(qp, w);
1243 }
1244 
1245 
1246 
1248  const Quadrature<1> &angular_quadrature)
1249  : QSimplex<2>(Quadrature<2>())
1250 {
1251  const QAnisotropic<2> base(radial_quadrature, angular_quadrature);
1252  this->quadrature_points.resize(base.size());
1253  this->weights.resize(base.size());
1254  for (unsigned int i = 0; i < base.size(); ++i)
1255  {
1256  const auto q = base.point(i);
1257  const auto w = base.weight(i);
1258 
1259  const auto xhat = q[0];
1260  const auto yhat = q[1];
1261 
1262  const double t = numbers::PI_2 * yhat;
1263  const double pi = numbers::PI;
1264  const double st = std::sin(t);
1265  const double ct = std::cos(t);
1266  const double r = xhat / (st + ct);
1267 
1268  const double J = pi * xhat / (2 * (std::sin(pi * yhat) + 1));
1269 
1270  this->quadrature_points[i] = Point<2>(r * ct, r * st);
1271  this->weights[i] = w * J;
1272  }
1273 }
1274 
1275 
1276 
1277 QTrianglePolar::QTrianglePolar(const unsigned int n)
1278  : QTrianglePolar(QGauss<1>(n), QGauss<1>(n))
1279 {}
1280 
1281 
1282 
1283 QDuffy::QDuffy(const Quadrature<1> &radial_quadrature,
1284  const Quadrature<1> &angular_quadrature,
1285  const double beta)
1286  : QSimplex<2>(Quadrature<2>())
1287 {
1288  const QAnisotropic<2> base(radial_quadrature, angular_quadrature);
1289  this->quadrature_points.resize(base.size());
1290  this->weights.resize(base.size());
1291  for (unsigned int i = 0; i < base.size(); ++i)
1292  {
1293  const auto q = base.point(i);
1294  const auto w = base.weight(i);
1295 
1296  const auto xhat = q[0];
1297  const auto yhat = q[1];
1298 
1299  const double x = std::pow(xhat, beta) * (1 - yhat);
1300  const double y = std::pow(xhat, beta) * yhat;
1301 
1302  const double J = beta * std::pow(xhat, 2. * beta - 1.);
1303 
1304  this->quadrature_points[i] = Point<2>(x, y);
1305  this->weights[i] = w * J;
1306  }
1307 }
1308 
1309 
1310 
1311 QDuffy::QDuffy(const unsigned int n, const double beta)
1312  : QDuffy(QGauss<1>(n), QGauss<1>(n), beta)
1313 {}
1314 
1315 
1316 
1317 template <int dim>
1318 QSplit<dim>::QSplit(const QSimplex<dim> &base, const Point<dim> &split_point)
1319 {
1321  ExcMessage(
1322  "The split point should be inside the unit reference cell."));
1323 
1324  std::array<Point<dim>, dim + 1> vertices;
1325  vertices[0] = split_point;
1326 
1327  // Make a simplex from the split_point and the first dim vertices of each
1328  // face. In dimension three, we need to split the face in two triangles, so
1329  // we use once the first dim vertices of each face, and the second time the
1330  // the dim vertices of each face starting from 1.
1331  for (unsigned int f = 0; f < GeometryInfo<dim>::faces_per_cell; ++f)
1332  for (unsigned int start = 0; start < (dim > 2 ? 2 : 1); ++start)
1333  {
1334  for (unsigned int i = 0; i < dim; ++i)
1335  vertices[i + 1] = GeometryInfo<dim>::unit_cell_vertex(
1337  const auto quad = base.compute_affine_transformation(vertices);
1338  if (quad.size())
1339  {
1340  this->quadrature_points.insert(this->quadrature_points.end(),
1341  quad.get_points().begin(),
1342  quad.get_points().end());
1343  this->weights.insert(this->weights.end(),
1344  quad.get_weights().begin(),
1345  quad.get_weights().end());
1346  }
1347  }
1348 }
1349 
1350 
1351 
1352 // explicit specialization
1353 // note that 1d formulae are specialized by implementation above
1354 template class QGauss<2>;
1355 template class QGaussLobatto<2>;
1356 template class QMidpoint<2>;
1357 template class QTrapez<2>;
1358 template class QSimpson<2>;
1359 template class QMilne<2>;
1360 template class QWeddle<2>;
1361 
1362 template class QGauss<3>;
1363 template class QGaussLobatto<3>;
1364 template class QMidpoint<3>;
1365 template class QTrapez<3>;
1366 template class QSimpson<3>;
1367 template class QMilne<3>;
1368 template class QWeddle<3>;
1369 
1370 template class QSorted<1>;
1371 template class QSorted<2>;
1372 template class QSorted<3>;
1373 
1374 template class QTelles<1>;
1375 template class QTelles<2>;
1376 template class QTelles<3>;
1377 
1378 template class QGaussChebyshev<1>;
1379 template class QGaussChebyshev<2>;
1380 template class QGaussChebyshev<3>;
1381 
1382 template class QGaussRadauChebyshev<1>;
1383 template class QGaussRadauChebyshev<2>;
1384 template class QGaussRadauChebyshev<3>;
1385 
1386 template class QGaussLobattoChebyshev<1>;
1387 template class QGaussLobattoChebyshev<2>;
1388 template class QGaussLobattoChebyshev<3>;
1389 
1390 template class QSimplex<1>;
1391 template class QSimplex<2>;
1392 template class QSimplex<3>;
1393 
1394 template class QSplit<1>;
1395 template class QSplit<2>;
1396 template class QSplit<3>;
1397 
1398 DEAL_II_NAMESPACE_CLOSE
QGaussLog(const unsigned int n, const bool revert=false)
std::vector< double > weights
Definition: quadrature.h:290
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)
QGaussChebyshev(const unsigned int n)
Generate a formula with n quadrature points.
const Point< dim > & point(const unsigned int i) const
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:548
static Point< dim > unit_cell_vertex(const unsigned int vertex)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1523
constexpr SymmetricTensor< rank_, dim, Number > transpose(const SymmetricTensor< rank_, dim, Number > &t)
QGauss(const unsigned int n)
static ::ExceptionBase & ExcIndexRange(int arg1, int arg2, int arg3)
Definition: point.h:110
static ::ExceptionBase & ExcMessage(std::string arg1)
#define Assert(cond, exc)
Definition: exceptions.h:1411
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)
bool compare_weights(const unsigned int a, const unsigned int b) const
QGaussLobatto(const unsigned int n)
static constexpr double PI_2
Definition: numbers.h:243
bool is_tensor_product_flag
Definition: quadrature.h:299
std::vector< Point< dim > > quadrature_points
Definition: quadrature.h:284
unsigned int size() const
QTelles(const Quadrature< 1 > &base_quad, const Point< dim > &singularity)
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Definition: tensor.h:422
static constexpr double PI
Definition: numbers.h:238
DEAL_II_CONSTEXPR Number determinant(const SymmetricTensor< 2, dim, Number > &)
QGaussLobattoChebyshev(const unsigned int n)
Generate a formula with n quadrature points.
static ::ExceptionBase & ExcNotImplemented()
QSorted(const Quadrature< dim > &quad)
QGaussLogR(const unsigned int n, const Point< dim > x0=Point< dim >(), const double alpha=1, const bool factor_out_singular_weight=false)
QSplit(const QSimplex< dim > &base, const Point< dim > &split_point)
Quadrature< dim > compute_affine_transformation(const std::array< Point< dim >, dim+1 > &vertices) const
double weight(const unsigned int i) const
static unsigned int quad_size(const Point< dim > singularity, const unsigned int n)
static ::ExceptionBase & ExcInternalError()
QTrianglePolar(const Quadrature< 1 > &radial_quadrature, const Quadrature< 1 > &angular_quadrature)