Reference documentation for deal.II version 9.2.0
\(\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 - 2020 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 #ifndef DOXYGEN
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 #endif
161 
162 
163 template <>
164 QMidpoint<1>::QMidpoint()
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 <>
174 QTrapez<1>::QTrapez()
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 <>
190 QSimpson<1>::QSimpson()
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 <>
206 QMilne<1>::QMilne()
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 <>
223 QWeddle<1>::QWeddle()
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  const double eps = 1e-8;
612  const bool on_edge =
613  std::any_of(singularity.begin_raw(),
614  singularity.end_raw(),
615  [eps](double coord) {
616  return std::abs(coord) < eps || std::abs(coord - 1.) < eps;
617  });
618  const bool on_vertex =
619  on_edge &&
620  std::abs((singularity - Point<2>(.5, .5)).norm_square() - .5) < eps;
621  if (on_vertex)
622  return 2 * n * n;
623  else if (on_edge)
624  return 4 * n * n;
625  else
626  return 8 * n * n;
627 }
628 
629 template <>
630 QGaussOneOverR<2>::QGaussOneOverR(const unsigned int n,
631  const Point<2> singularity,
632  const bool factor_out_singularity)
633  : Quadrature<2>(quad_size(singularity, n))
634 {
635  // We treat all the cases in the
636  // same way. Split the element in 4
637  // pieces, measure the area, if
638  // it's relevant, add the
639  // quadrature connected to that
640  // singularity.
641  std::vector<QGaussOneOverR<2>> quads;
642  std::vector<Point<2>> origins;
643  // Id of the corner with a
644  // singularity
645  quads.emplace_back(n, 3, factor_out_singularity);
646  quads.emplace_back(n, 2, factor_out_singularity);
647  quads.emplace_back(n, 1, factor_out_singularity);
648  quads.emplace_back(n, 0, factor_out_singularity);
649 
650  origins.emplace_back(0., 0.);
651  origins.emplace_back(singularity[0], 0.);
652  origins.emplace_back(0., singularity[1]);
653  origins.push_back(singularity);
654 
655  // Lexicographical ordering.
656 
657  double eps = 1e-8;
658  unsigned int q_id = 0; // Current quad point index.
659  Tensor<1, 2> dist;
660 
661  for (unsigned int box = 0; box < 4; ++box)
662  {
663  dist = (singularity - GeometryInfo<2>::unit_cell_vertex(box));
664  dist = Point<2>(std::abs(dist[0]), std::abs(dist[1]));
665  double area = dist[0] * dist[1];
666  if (area > eps)
667  for (unsigned int q = 0; q < quads[box].size(); ++q, ++q_id)
668  {
669  const Point<2> &qp = quads[box].point(q);
670  this->quadrature_points[q_id] =
671  origins[box] + Point<2>(dist[0] * qp[0], dist[1] * qp[1]);
672  this->weights[q_id] = quads[box].weight(q) * area;
673  }
674  }
675 }
676 
677 
678 template <>
679 QGaussOneOverR<2>::QGaussOneOverR(const unsigned int n,
680  const unsigned int vertex_index,
681  const bool factor_out_singularity)
682  : Quadrature<2>(2 * n * n)
683 {
684  // This version of the constructor
685  // works only for the 4
686  // vertices. If you need a more
687  // general one, you should use the
688  // one with the Point<2> in the
689  // constructor.
690  AssertIndexRange(vertex_index, 4);
691 
692  // Start with the gauss quadrature formula on the (u,v) reference
693  // element.
694  QGauss<2> gauss(n);
695 
696  Assert(gauss.size() == n * n, ExcInternalError());
697 
698  // For the moment we only implemented this for the vertices of a
699  // quadrilateral. We are planning to do this also for the support
700  // points of arbitrary FE_Q elements, to allow the use of this
701  // class in boundary element programs with higher order mappings.
702  AssertIndexRange(vertex_index, 4);
703 
704  // We create only the first one. All other pieces are rotation of
705  // this one.
706  // In this case the transformation is
707  //
708  // (x,y) = (u, u tan(pi/4 v))
709  //
710  // with Jacobian
711  //
712  // J = pi/4 R / cos(pi/4 v)
713  //
714  // And we get rid of R to take into account the singularity,
715  // unless specified differently in the constructor.
716  std::vector<Point<2>> &ps = this->quadrature_points;
717  std::vector<double> & ws = this->weights;
718  double pi4 = numbers::PI / 4;
719 
720  for (unsigned int q = 0; q < gauss.size(); ++q)
721  {
722  const Point<2> &gp = gauss.point(q);
723  ps[q][0] = gp[0];
724  ps[q][1] = gp[0] * std::tan(pi4 * gp[1]);
725  ws[q] = gauss.weight(q) * pi4 / std::cos(pi4 * gp[1]);
726  if (factor_out_singularity)
727  ws[q] *= (ps[q] - GeometryInfo<2>::unit_cell_vertex(0)).norm();
728  // The other half of the quadrilateral is symmetric with
729  // respect to xy plane.
730  ws[gauss.size() + q] = ws[q];
731  ps[gauss.size() + q][0] = ps[q][1];
732  ps[gauss.size() + q][1] = ps[q][0];
733  }
734 
735  // Now we distribute these vertices in the correct manner
736  double theta = 0;
737  switch (vertex_index)
738  {
739  case 0:
740  theta = 0;
741  break;
742  case 1:
743  //
744  theta = numbers::PI / 2;
745  break;
746  case 2:
747  theta = -numbers::PI / 2;
748  break;
749  case 3:
750  theta = numbers::PI;
751  break;
752  }
753 
754  double R00 = std::cos(theta), R01 = -std::sin(theta);
755  double R10 = std::sin(theta), R11 = std::cos(theta);
756 
757  if (vertex_index != 0)
758  for (unsigned int q = 0; q < size(); ++q)
759  {
760  double x = ps[q][0] - .5, y = ps[q][1] - .5;
761 
762  ps[q][0] = R00 * x + R01 * y + .5;
763  ps[q][1] = R10 * x + R11 * y + .5;
764  }
765 }
766 
767 
768 template <int dim>
769 QSorted<dim>::QSorted(const Quadrature<dim> &quad)
770  : Quadrature<dim>(quad)
771 {
772  std::vector<unsigned int> permutation(quad.size());
773  for (unsigned int i = 0; i < quad.size(); ++i)
774  permutation[i] = i;
775 
776  std::sort(permutation.begin(),
777  permutation.end(),
778  [this](const unsigned int x, const unsigned int y) {
779  return this->compare_weights(x, y);
780  });
781 
782  // At this point, the variable is_tensor_product_flag is set
783  // to the respective value of the given Quadrature in the base
784  // class copy constructor.
785  // We only call a quadrature formula 'tensor product'
786  // if the quadrature points are also sorted lexicographically.
787  // In particular, any reordering destroys that property
788  // and we might need to modify the variable accordingly.
789  for (unsigned int i = 0; i < quad.size(); ++i)
790  {
791  this->weights[i] = quad.weight(permutation[i]);
792  this->quadrature_points[i] = quad.point(permutation[i]);
793  if (permutation[i] != i)
794  this->is_tensor_product_flag = false;
795  }
796 }
797 
798 
799 template <int dim>
800 bool
801 QSorted<dim>::compare_weights(const unsigned int a, const unsigned int b) const
802 {
803  return (this->weights[a] < this->weights[b]);
804 }
805 
806 
807 // construct the quadrature formulae in higher dimensions by
808 // tensor product of lower dimensions
809 
810 template <int dim>
811 QGauss<dim>::QGauss(const unsigned int n)
812  : Quadrature<dim>(QGauss<dim - 1>(n), QGauss<1>(n))
813 {}
814 
815 
816 
817 template <int dim>
818 QGaussLobatto<dim>::QGaussLobatto(const unsigned int n)
819  : Quadrature<dim>(QGaussLobatto<dim - 1>(n), QGaussLobatto<1>(n))
820 {}
821 
822 
823 
824 template <int dim>
825 QMidpoint<dim>::QMidpoint()
826  : Quadrature<dim>(QMidpoint<dim - 1>(), QMidpoint<1>())
827 {}
828 
829 
830 
831 template <int dim>
832 QTrapez<dim>::QTrapez()
833  : Quadrature<dim>(QTrapez<dim - 1>(), QTrapez<1>())
834 {}
835 
836 
837 
838 template <int dim>
839 QSimpson<dim>::QSimpson()
840  : Quadrature<dim>(QSimpson<dim - 1>(), QSimpson<1>())
841 {}
842 
843 
844 
845 template <int dim>
846 QMilne<dim>::QMilne()
847  : Quadrature<dim>(QMilne<dim - 1>(), QMilne<1>())
848 {}
849 
850 
851 template <int dim>
852 QWeddle<dim>::QWeddle()
853  : Quadrature<dim>(QWeddle<dim - 1>(), QWeddle<1>())
854 {}
855 
856 template <int dim>
857 QTelles<dim>::QTelles(const Quadrature<1> &base_quad,
858  const Point<dim> & singularity)
859  : // We need the explicit implementation if dim == 1. If dim > 1 we use the
860  // former implementation and apply a tensorial product to obtain the higher
861  // dimensions.
862  Quadrature<dim>(
863  dim == 2 ?
864  QAnisotropic<dim>(QTelles<1>(base_quad, Point<1>(singularity[0])),
865  QTelles<1>(base_quad, Point<1>(singularity[1]))) :
866  dim == 3 ?
867  QAnisotropic<dim>(QTelles<1>(base_quad, Point<1>(singularity[0])),
868  QTelles<1>(base_quad, Point<1>(singularity[1])),
869  QTelles<1>(base_quad, Point<1>(singularity[2]))) :
870  Quadrature<dim>())
871 {}
872 
873 template <int dim>
874 QTelles<dim>::QTelles(const unsigned int n, const Point<dim> &singularity)
875  : // In this case we map the standard Gauss Legendre formula using the given
876  // singularity point coordinates.
877  Quadrature<dim>(QTelles<dim>(QGauss<1>(n), singularity))
878 {}
879 
880 
881 
882 template <>
883 QTelles<1>::QTelles(const Quadrature<1> &base_quad, const Point<1> &singularity)
884  : // We explicitly implement the Telles' variable change if dim == 1.
885  Quadrature<1>(base_quad)
886 {
887  // We define all the constants to be used in the implementation of
888  // Telles' rule
889  const double eta_bar = singularity[0] * 2. - 1.;
890  const double eta_star = eta_bar * eta_bar - 1.;
891  double gamma_bar;
892 
893  std::vector<Point<1>> quadrature_points_dummy(quadrature_points.size());
894  std::vector<double> weights_dummy(weights.size());
895  unsigned int cont = 0;
896  const double tol = 1e-10;
897  for (unsigned int d = 0; d < quadrature_points.size(); ++d)
898  {
899  if (std::abs(quadrature_points[d][0] - singularity[0]) > tol)
900  {
901  quadrature_points_dummy[d - cont] = quadrature_points[d];
902  weights_dummy[d - cont] = weights[d];
903  }
904  else
905  {
906  // We need to remove the singularity point from the quadrature point
907  // list. To do so we use the variable cont.
908  cont = 1;
909  }
910  }
911  if (cont == 1)
912  {
913  quadrature_points.resize(quadrature_points_dummy.size() - 1);
914  weights.resize(weights_dummy.size() - 1);
915  for (unsigned int d = 0; d < quadrature_points.size(); ++d)
916  {
917  quadrature_points[d] = quadrature_points_dummy[d];
918  weights[d] = weights_dummy[d];
919  }
920  }
921  // We need to check if the singularity is at the boundary of the interval.
922  if (std::abs(eta_star) <= tol)
923  {
924  gamma_bar =
925  std::pow((eta_bar * eta_star + std::abs(eta_star)), 1.0 / 3.0) +
926  std::pow((eta_bar * eta_star - std::abs(eta_star)), 1.0 / 3.0) +
927  eta_bar;
928  }
929  else
930  {
931  gamma_bar = (eta_bar * eta_star + std::abs(eta_star)) /
932  std::abs(eta_bar * eta_star + std::abs(eta_star)) *
933  std::pow(std::abs(eta_bar * eta_star + std::abs(eta_star)),
934  1.0 / 3.0) +
935  (eta_bar * eta_star - std::abs(eta_star)) /
936  std::abs(eta_bar * eta_star - std::abs(eta_star)) *
937  std::pow(std::abs(eta_bar * eta_star - std::abs(eta_star)),
938  1.0 / 3.0) +
939  eta_bar;
940  }
941  for (unsigned int q = 0; q < quadrature_points.size(); ++q)
942  {
943  double gamma = quadrature_points[q][0] * 2 - 1;
944  double eta = (std::pow(gamma - gamma_bar, 3.0) +
945  gamma_bar * (gamma_bar * gamma_bar + 3)) /
946  (1 + 3 * gamma_bar * gamma_bar);
947 
948  double J = 3 * ((gamma - gamma_bar) * (gamma - gamma_bar)) /
949  (1 + 3 * gamma_bar * gamma_bar);
950 
951  quadrature_points[q][0] = (eta + 1) / 2.0;
952  weights[q] = J * weights[q];
953  }
954 }
955 
956 namespace internal
957 {
958  namespace QGaussChebyshev
959  {
963  std::vector<double>
964  get_quadrature_points(const unsigned int n)
965  {
966  std::vector<double> points(n);
967  // n point quadrature: index from 0 to n-1
968  for (unsigned short i = 0; i < n; ++i)
969  // would be cos((2i+1)Pi/(2N+2))
970  // put + Pi so we start from the smallest point
971  // then map from [-1,1] to [0,1]
972  points[i] =
973  1. / 2. *
974  (1. + std::cos(numbers::PI *
975  (1. + double(2 * i + 1) / double(2 * (n - 1) + 2))));
976 
977  return points;
978  }
979 
980 
981 
985  std::vector<double>
986  get_quadrature_weights(const unsigned int n)
987  {
988  std::vector<double> weights(n);
989 
990  for (unsigned short i = 0; i < n; ++i)
991  {
992  // same weights as on [-1,1]
993  weights[i] = numbers::PI / double(n);
994  }
995 
996  return weights;
997  }
998  } // namespace QGaussChebyshev
999 } // namespace internal
1000 
1001 
1002 template <>
1003 QGaussChebyshev<1>::QGaussChebyshev(const unsigned int n)
1004  : Quadrature<1>(n)
1005 {
1006  Assert(n > 0, ExcMessage("Need at least one point for the quadrature rule"));
1007  std::vector<double> p = internal::QGaussChebyshev::get_quadrature_points(n);
1008  std::vector<double> w = internal::QGaussChebyshev::get_quadrature_weights(n);
1009 
1010  for (unsigned int i = 0; i < this->size(); ++i)
1011  {
1012  this->quadrature_points[i] = Point<1>(p[i]);
1013  this->weights[i] = w[i];
1014  }
1015 }
1016 
1017 
1018 template <int dim>
1019 QGaussChebyshev<dim>::QGaussChebyshev(const unsigned int n)
1020  : Quadrature<dim>(QGaussChebyshev<1>(n))
1021 {}
1022 
1023 
1024 namespace internal
1025 {
1026  namespace QGaussRadauChebyshev
1027  {
1028  // Computes the points of the quadrature formula.
1029  std::vector<double>
1030  get_quadrature_points(const unsigned int n,
1031  ::QGaussRadauChebyshev<1>::EndPoint ep)
1032  {
1033  std::vector<double> points(n);
1034  // n point quadrature: index from 0 to n-1
1035  for (unsigned short i = 0; i < n; ++i)
1036  // would be -cos(2i Pi/(2N+1))
1037  // put + Pi so we start from the smallest point
1038  // then map from [-1,1] to [0,1]
1039  switch (ep)
1040  {
1041  case ::QGaussRadauChebyshev<1>::left:
1042  {
1043  points[i] =
1044  1. / 2. *
1045  (1. -
1046  std::cos(numbers::PI *
1047  (1 + 2 * double(i) / (2 * double(n - 1) + 1.))));
1048  break;
1049  }
1050 
1051  case ::QGaussRadauChebyshev<1>::right:
1052  {
1053  points[i] =
1054  1. / 2. *
1055  (1. - std::cos(numbers::PI * (2 * double(n - 1 - i) /
1056  (2 * double(n - 1) + 1.))));
1057  break;
1058  }
1059 
1060  default:
1061  Assert(
1062  false,
1063  ExcMessage(
1064  "This constructor can only be called with either "
1065  "QGaussRadauChebyshev::left or QGaussRadauChebyshev::right as "
1066  "second argument."));
1067  }
1068 
1069  return points;
1070  }
1071 
1072 
1073 
1074  // Computes the weights of the quadrature formula.
1075  std::vector<double>
1076  get_quadrature_weights(const unsigned int n,
1077  ::QGaussRadauChebyshev<1>::EndPoint ep)
1078  {
1079  std::vector<double> weights(n);
1080 
1081  for (unsigned short i = 0; i < n; ++i)
1082  {
1083  // same weights as on [-1,1]
1084  weights[i] = 2. * numbers::PI / double(2 * (n - 1) + 1.);
1085  if (ep == ::QGaussRadauChebyshev<1>::left && i == 0)
1086  weights[i] /= 2.;
1087  else if (ep == ::QGaussRadauChebyshev<1>::right &&
1088  i == (n - 1))
1089  weights[i] /= 2.;
1090  }
1091 
1092  return weights;
1093  }
1094  } // namespace QGaussRadauChebyshev
1095 } // namespace internal
1096 
1097 
1098 template <>
1099 QGaussRadauChebyshev<1>::QGaussRadauChebyshev(const unsigned int n, EndPoint ep)
1100  : Quadrature<1>(n)
1101  , ep(ep)
1102 {
1103  Assert(n > 0, ExcMessage("Need at least one point for quadrature rules"));
1104  std::vector<double> p =
1105  internal::QGaussRadauChebyshev::get_quadrature_points(n, ep);
1106  std::vector<double> w =
1107  internal::QGaussRadauChebyshev::get_quadrature_weights(n, ep);
1108 
1109  for (unsigned int i = 0; i < this->size(); ++i)
1110  {
1111  this->quadrature_points[i] = Point<1>(p[i]);
1112  this->weights[i] = w[i];
1113  }
1114 }
1115 
1116 
1117 template <int dim>
1118 QGaussRadauChebyshev<dim>::QGaussRadauChebyshev(const unsigned int n,
1119  EndPoint ep)
1120  : Quadrature<dim>(QGaussRadauChebyshev<1>(
1121  n,
1122  static_cast<QGaussRadauChebyshev<1>::EndPoint>(ep)))
1123  , ep(ep)
1124 {}
1125 
1126 
1127 
1128 namespace internal
1129 {
1130  namespace QGaussLobattoChebyshev
1131  {
1132  // Computes the points of the quadrature formula.
1133  std::vector<double>
1134  get_quadrature_points(const unsigned int n)
1135  {
1136  std::vector<double> points(n);
1137  // n point quadrature: index from 0 to n-1
1138  for (unsigned short i = 0; i < n; ++i)
1139  // would be cos(i Pi/N)
1140  // put + Pi so we start from the smallest point
1141  // then map from [-1,1] to [0,1]
1142  points[i] =
1143  1. / 2. *
1144  (1. + std::cos(numbers::PI * (1 + double(i) / double(n - 1))));
1145 
1146  return points;
1147  }
1148 
1149  // Computes the weights of the quadrature formula.
1150  std::vector<double>
1151  get_quadrature_weights(const unsigned int n)
1152  {
1153  std::vector<double> weights(n);
1154 
1155  for (unsigned short i = 0; i < n; ++i)
1156  {
1157  // same weights as on [-1,1]
1158  weights[i] = numbers::PI / double((n - 1));
1159  if (i == 0 || i == (n - 1))
1160  weights[i] /= 2.;
1161  }
1162 
1163  return weights;
1164  }
1165  } // namespace QGaussLobattoChebyshev
1166 } // namespace internal
1167 
1168 
1169 
1170 template <>
1171 QGaussLobattoChebyshev<1>::QGaussLobattoChebyshev(const unsigned int n)
1172  : Quadrature<1>(n)
1173 {
1174  Assert(n > 1,
1175  ExcMessage(
1176  "Need at least two points for Gauss-Lobatto quadrature rule"));
1177  std::vector<double> p =
1178  internal::QGaussLobattoChebyshev::get_quadrature_points(n);
1179  std::vector<double> w =
1180  internal::QGaussLobattoChebyshev::get_quadrature_weights(n);
1181 
1182  for (unsigned int i = 0; i < this->size(); ++i)
1183  {
1184  this->quadrature_points[i] = Point<1>(p[i]);
1185  this->weights[i] = w[i];
1186  }
1187 }
1188 
1189 
1190 template <int dim>
1191 QGaussLobattoChebyshev<dim>::QGaussLobattoChebyshev(const unsigned int n)
1192  : Quadrature<dim>(QGaussLobattoChebyshev<1>(n))
1193 {}
1194 
1195 
1196 
1197 template <int dim>
1198 QSimplex<dim>::QSimplex(const Quadrature<dim> &quad)
1199 {
1200  std::vector<Point<dim>> qpoints;
1201  std::vector<double> weights;
1202 
1203  for (unsigned int i = 0; i < quad.size(); ++i)
1204  {
1205  double r = 0;
1206  for (unsigned int d = 0; d < dim; ++d)
1207  r += quad.point(i)[d];
1208  if (r <= 1 + 1e-10)
1209  {
1210  this->quadrature_points.push_back(quad.point(i));
1211  this->weights.push_back(quad.weight(i));
1212  }
1213  }
1214 }
1215 
1216 
1217 
1218 template <int dim>
1219 Quadrature<dim>
1220 QSimplex<dim>::compute_affine_transformation(
1221  const std::array<Point<dim>, dim + 1> &vertices) const
1222 {
1223  Tensor<2, dim> B;
1224  for (unsigned int d = 0; d < dim; ++d)
1225  B[d] = vertices[d + 1] - vertices[0];
1226 
1227  B = transpose(B);
1228  const double J = std::abs(determinant(B));
1229 
1230  // if the determinant is zero, we return an empty quadrature
1231  if (J < 1e-12)
1232  return Quadrature<dim>();
1233 
1234  std::vector<Point<dim>> qp(this->size());
1235  std::vector<double> w(this->size());
1236 
1237  for (unsigned int i = 0; i < this->size(); ++i)
1238  {
1239  qp[i] = Point<dim>(vertices[0] + B * this->point(i));
1240  w[i] = J * this->weight(i);
1241  }
1242 
1243  return Quadrature<dim>(qp, w);
1244 }
1245 
1246 
1247 
1248 QTrianglePolar::QTrianglePolar(const Quadrature<1> &radial_quadrature,
1249  const Quadrature<1> &angular_quadrature)
1250  : QSimplex<2>(Quadrature<2>())
1251 {
1252  const QAnisotropic<2> base(radial_quadrature, angular_quadrature);
1253  this->quadrature_points.resize(base.size());
1254  this->weights.resize(base.size());
1255  for (unsigned int i = 0; i < base.size(); ++i)
1256  {
1257  const auto q = base.point(i);
1258  const auto w = base.weight(i);
1259 
1260  const auto xhat = q[0];
1261  const auto yhat = q[1];
1262 
1263  const double t = numbers::PI_2 * yhat;
1264  const double pi = numbers::PI;
1265  const double st = std::sin(t);
1266  const double ct = std::cos(t);
1267  const double r = xhat / (st + ct);
1268 
1269  const double J = pi * xhat / (2 * (std::sin(pi * yhat) + 1));
1270 
1271  this->quadrature_points[i] = Point<2>(r * ct, r * st);
1272  this->weights[i] = w * J;
1273  }
1274 }
1275 
1276 
1277 
1278 QTrianglePolar::QTrianglePolar(const unsigned int n)
1279  : QTrianglePolar(QGauss<1>(n), QGauss<1>(n))
1280 {}
1281 
1282 
1283 
1284 QDuffy::QDuffy(const Quadrature<1> &radial_quadrature,
1285  const Quadrature<1> &angular_quadrature,
1286  const double beta)
1287  : QSimplex<2>(Quadrature<2>())
1288 {
1289  const QAnisotropic<2> base(radial_quadrature, angular_quadrature);
1290  this->quadrature_points.resize(base.size());
1291  this->weights.resize(base.size());
1292  for (unsigned int i = 0; i < base.size(); ++i)
1293  {
1294  const auto q = base.point(i);
1295  const auto w = base.weight(i);
1296 
1297  const auto xhat = q[0];
1298  const auto yhat = q[1];
1299 
1300  const double x = std::pow(xhat, beta) * (1 - yhat);
1301  const double y = std::pow(xhat, beta) * yhat;
1302 
1303  const double J = beta * std::pow(xhat, 2. * beta - 1.);
1304 
1305  this->quadrature_points[i] = Point<2>(x, y);
1306  this->weights[i] = w * J;
1307  }
1308 }
1309 
1310 
1311 
1312 QDuffy::QDuffy(const unsigned int n, const double beta)
1313  : QDuffy(QGauss<1>(n), QGauss<1>(n), beta)
1314 {}
1315 
1316 
1317 
1318 template <int dim>
1319 QSplit<dim>::QSplit(const QSimplex<dim> &base, const Point<dim> &split_point)
1320 {
1321  AssertThrow(GeometryInfo<dim>::is_inside_unit_cell(split_point, 1e-12),
1322  ExcMessage(
1323  "The split point should be inside the unit reference cell."));
1324 
1325  std::array<Point<dim>, dim + 1> vertices;
1326  vertices[0] = split_point;
1327 
1328  // Make a simplex from the split_point and the first dim vertices of each
1329  // face. In dimension three, we need to split the face in two triangles, so
1330  // we use once the first dim vertices of each face, and the second time the
1331  // the dim vertices of each face starting from 1.
1332  for (auto f : GeometryInfo<dim>::face_indices())
1333  for (unsigned int start = 0; start < (dim > 2 ? 2 : 1); ++start)
1334  {
1335  for (unsigned int i = 0; i < dim; ++i)
1336  vertices[i + 1] = GeometryInfo<dim>::unit_cell_vertex(
1337  GeometryInfo<dim>::face_to_cell_vertices(f, i + start));
1338  const auto quad = base.compute_affine_transformation(vertices);
1339  if (quad.size())
1340  {
1341  this->quadrature_points.insert(this->quadrature_points.end(),
1342  quad.get_points().begin(),
1343  quad.get_points().end());
1344  this->weights.insert(this->weights.end(),
1345  quad.get_weights().begin(),
1346  quad.get_weights().end());
1347  }
1348  }
1349 }
1350 
1351 
1352 
1353 // explicit specialization
1354 // note that 1d formulae are specialized by implementation above
1355 template class QGauss<2>;
1356 template class QGaussLobatto<2>;
1357 template class QMidpoint<2>;
1358 template class QTrapez<2>;
1359 template class QSimpson<2>;
1360 template class QMilne<2>;
1361 template class QWeddle<2>;
1362 
1363 template class QGauss<3>;
1364 template class QGaussLobatto<3>;
1365 template class QMidpoint<3>;
1366 template class QTrapez<3>;
1367 template class QSimpson<3>;
1368 template class QMilne<3>;
1369 template class QWeddle<3>;
1370 
1371 template class QSorted<1>;
1372 template class QSorted<2>;
1373 template class QSorted<3>;
1374 
1375 template class QTelles<1>;
1376 template class QTelles<2>;
1377 template class QTelles<3>;
1378 
1379 template class QGaussChebyshev<1>;
1380 template class QGaussChebyshev<2>;
1381 template class QGaussChebyshev<3>;
1382 
1383 template class QGaussRadauChebyshev<1>;
1384 template class QGaussRadauChebyshev<2>;
1385 template class QGaussRadauChebyshev<3>;
1386 
1387 template class QGaussLobattoChebyshev<1>;
1388 template class QGaussLobattoChebyshev<2>;
1389 template class QGaussLobattoChebyshev<3>;
1390 
1391 template class QSimplex<1>;
1392 template class QSimplex<2>;
1393 template class QSimplex<3>;
1394 
1395 template class QSplit<1>;
1396 template class QSplit<2>;
1397 template class QSplit<3>;
1398 
1399 DEAL_II_NAMESPACE_CLOSE
QAnisotropic
Definition: quadrature.h:320
QSimpson
Definition: quadrature_lib.h:105
internal::QGaussLobatto::gamma
long double gamma(const unsigned int n)
Definition: quadrature_lib.cc:96
QSimplex::QSimplex
QSimplex(const Quadrature< dim > &quad)
Definition: quadrature_lib.cc:1198
QDuffy::QDuffy
QDuffy(const Quadrature< 1 > &radial_quadrature, const Quadrature< 1 > &angular_quadrature, const double beta=1.0)
Definition: quadrature_lib.cc:1284
GeometryInfo::unit_cell_vertex
static Point< dim > unit_cell_vertex(const unsigned int vertex)
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
QGaussRadauChebyshev::QGaussRadauChebyshev
QGaussRadauChebyshev(const unsigned int n, EndPoint ep=QGaussRadauChebyshev::left)
Generate a formula with n quadrature points.
Definition: quadrature_lib.cc:1118
QSorted
Definition: quadrature_lib.h:391
QMilne::QMilne
QMilne()
Definition: quadrature_lib.cc:846
polynomial.h
internal::QGaussRadauChebyshev::get_quadrature_points
std::vector< double > get_quadrature_points(const unsigned int n, ::QGaussRadauChebyshev< 1 >::EndPoint ep)
Definition: quadrature_lib.cc:1030
Quadrature::quadrature_points
std::vector< Point< dim > > quadrature_points
Definition: quadrature.h:284
QSimplex::compute_affine_transformation
Quadrature< dim > compute_affine_transformation(const std::array< Point< dim >, dim+1 > &vertices) const
Definition: quadrature_lib.cc:1220
internal::QGaussLobattoChebyshev::get_quadrature_weights
std::vector< double > get_quadrature_weights(const unsigned int n)
Definition: quadrature_lib.cc:1151
QGaussOneOverR::quad_size
static unsigned int quad_size(const Point< dim > singularity, const unsigned int n)
Tensor::begin_raw
Number * begin_raw()
GeometryInfo
Definition: geometry_info.h:1224
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1649
QGaussLobatto::QGaussLobatto
QGaussLobatto(const unsigned int n)
Definition: quadrature_lib.cc:818
Quadrature::weight
double weight(const unsigned int i) const
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
QGaussLobattoChebyshev::QGaussLobattoChebyshev
QGaussLobattoChebyshev(const unsigned int n)
Generate a formula with n quadrature points.
Definition: quadrature_lib.cc:1191
QGaussLobattoChebyshev
Definition: quadrature_lib.h:571
QGaussLog
Definition: quadrature_lib.h:177
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
QDuffy
Definition: quadrature_lib.h:726
QGaussLog::get_quadrature_points
static std::vector< double > get_quadrature_points(const unsigned int n)
numbers::PI_2
static constexpr double PI_2
Definition: numbers.h:242
QTelles::QTelles
QTelles(const Quadrature< 1 > &base_quad, const Point< dim > &singularity)
Definition: quadrature_lib.cc:857
quadrature_lib.h
OpenCASCADE::point
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:188
QWeddle::QWeddle
QWeddle()
Definition: quadrature_lib.cc:852
QSorted::QSorted
QSorted(const Quadrature< dim > &quad)
Definition: quadrature_lib.cc:769
QSimplex
Definition: quadrature_lib.h:611
Polynomials::jacobi_polynomial_value
Number jacobi_polynomial_value(const unsigned int degree, const int alpha, const int beta, const Number x)
Definition: polynomial.h:845
Quadrature::point
const Point< dim > & point(const unsigned int i) const
QTrapez
Definition: quadrature_lib.h:126
Physics::Elasticity::Kinematics::w
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
QGaussLog::get_quadrature_weights
static std::vector< double > get_quadrature_weights(const unsigned int n)
QSplit::QSplit
QSplit(const QSimplex< dim > &base, const Point< dim > &split_point)
Definition: quadrature_lib.cc:1319
internal::QGaussRadauChebyshev::get_quadrature_weights
std::vector< double > get_quadrature_weights(const unsigned int n, ::QGaussRadauChebyshev< 1 >::EndPoint ep)
Definition: quadrature_lib.cc:1076
QTrapez::QTrapez
QTrapez()
Definition: quadrature_lib.cc:832
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
double
geometry_info.h
QGaussRadauChebyshev::ep
const EndPoint ep
Definition: quadrature_lib.h:552
QMilne
Definition: quadrature_lib.h:141
Tensor
Definition: tensor.h:450
internal::QGaussLobattoChebyshev::get_quadrature_points
std::vector< double > get_quadrature_points(const unsigned int n)
Definition: quadrature_lib.cc:1134
LocalIntegrators::Divergence::norm
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:548
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:358
QSplit
Definition: quadrature_lib.h:764
Physics::Elasticity::Kinematics::b
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
QGaussLog::QGaussLog
QGaussLog(const unsigned int n, const bool revert=false)
Quadrature::weights
std::vector< double > weights
Definition: quadrature.h:290
QTrianglePolar::QTrianglePolar
QTrianglePolar(const Quadrature< 1 > &radial_quadrature, const Quadrature< 1 > &angular_quadrature)
Definition: quadrature_lib.cc:1248
Quadrature::is_tensor_product_flag
bool is_tensor_product_flag
Definition: quadrature.h:299
QSimpson::QSimpson
QSimpson()
Definition: quadrature_lib.cc:839
QGaussLobatto
Definition: quadrature_lib.h:76
QMidpoint::QMidpoint
QMidpoint()
Definition: quadrature_lib.cc:825
QGaussLogR::QGaussLogR
QGaussLogR(const unsigned int n, const Point< dim > x0=Point< dim >(), const double alpha=1, const bool factor_out_singular_weight=false)
internal::QGaussChebyshev::get_quadrature_weights
std::vector< double > get_quadrature_weights(const unsigned int n)
Definition: quadrature_lib.cc:986
QGauss
Definition: quadrature_lib.h:40
vertices
Point< 3 > vertices[4]
Definition: data_out_base.cc:174
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
QSorted::compare_weights
bool compare_weights(const unsigned int a, const unsigned int b) const
Definition: quadrature_lib.cc:801
QMidpoint
Definition: quadrature_lib.h:93
DataOutBase::eps
@ eps
Definition: data_out_base.h:1582
Particles::Generators::quadrature_points
void quadrature_points(const Triangulation< dim, spacedim > &triangulation, const Quadrature< dim > &quadrature, const std::vector< std::vector< BoundingBox< spacedim >>> &global_bounding_boxes, ParticleHandler< dim, spacedim > &particle_handler, const Mapping< dim, spacedim > &mapping=StaticMappingQ1< dim, spacedim >::mapping)
Definition: generators.cc:442
QGaussChebyshev::QGaussChebyshev
QGaussChebyshev(const unsigned int n)
Generate a formula with n quadrature points.
Definition: quadrature_lib.cc:1019
QGaussOneOverR::QGaussOneOverR
QGaussOneOverR(const unsigned int n, const Point< dim > singularity, const bool factor_out_singular_weight=false)
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
QGauss::QGauss
QGauss(const unsigned int n)
Definition: quadrature_lib.cc:811
SymmetricTensor::determinant
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &t)
Definition: symmetric_tensor.h:2645
QWeddle
Definition: quadrature_lib.h:155
Point
Definition: point.h:111
QTelles
Definition: quadrature_lib.h:467
internal::QGaussLobatto::compute_quadrature_weights
std::vector< long double > compute_quadrature_weights(const std::vector< long double > &x, const int alpha, const int beta)
Definition: quadrature_lib.cc:114
Quadrature::size
unsigned int size() const
QGaussChebyshev
Definition: quadrature_lib.h:499
QGaussRadauChebyshev::EndPoint
EndPoint
Definition: quadrature_lib.h:530
internal
Definition: aligned_vector.h:369
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:359
Quadrature
Definition: quadrature.h:85
numbers::PI
static constexpr double PI
Definition: numbers.h:237
QGaussRadauChebyshev
Definition: quadrature_lib.h:524
QTrianglePolar
Definition: quadrature_lib.h:670
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1531
DerivativeForm::transpose
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Definition: derivative_form.h:470
QGaussLogR::fraction
const double fraction
Definition: quadrature_lib.h:265
Tensor::end_raw
Number * end_raw()
internal::TensorProductManifoldImplementation::split_point
void split_point(const Point< dim1+dim2 > &source, Point< dim1 > &p1, Point< dim2 > &p2)
Definition: tensor_product_manifold.h:161
internal::QGaussChebyshev::get_quadrature_points
std::vector< double > get_quadrature_points(const unsigned int n)
Definition: quadrature_lib.cc:964
Generated by   doxygen 1.8.17