Reference documentation for deal.II version 9.3.3
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
quadrature_lib.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 1998 - 2021 by the deal.II authors
4//
5// This file is part of the deal.II library.
6//
7// The deal.II library is free software; you can use it, redistribute
8// it, and/or modify it under the terms of the GNU Lesser General
9// Public License as published by the Free Software Foundation; either
10// version 2.1 of the License, or (at your option) any later version.
11// The full text of the license can be found in the file LICENSE.md at
12// the top level directory of deal.II.
13//
14// ---------------------------------------------------------------------
15
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
26DEAL_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
35template <>
36QGauss<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
48template <>
49QGaussLobatto<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
61template <>
62QGauss<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
87namespace 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
140template <>
141QGaussLobatto<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
163template <>
164QMidpoint<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
173template <>
174QTrapezoid<1>::QTrapezoid()
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
189template <>
190QSimpson<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
205template <>
206QMilne<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
222template <>
223QWeddle<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
244template <>
245QGaussLog<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
262template <>
263std::vector<double>
264QGaussLog<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
395template <>
396std::vector<double>
397QGaussLog<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
527template <>
528QGaussLogR<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
607template <>
608unsigned int
609QGaussOneOverR<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
629template <>
630QGaussOneOverR<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
678template <>
679QGaussOneOverR<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
768template <int dim>
769QSorted<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
799template <int dim>
800bool
801QSorted<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
810template <int dim>
811QGauss<dim>::QGauss(const unsigned int n)
812 : Quadrature<dim>(QGauss<dim - 1>(n), QGauss<1>(n))
813{}
814
815
816
817template <int dim>
818QGaussLobatto<dim>::QGaussLobatto(const unsigned int n)
819 : Quadrature<dim>(QGaussLobatto<dim - 1>(n), QGaussLobatto<1>(n))
820{}
821
822
823
824template <int dim>
825QMidpoint<dim>::QMidpoint()
826 : Quadrature<dim>(QMidpoint<dim - 1>(), QMidpoint<1>())
827{}
828
829
830
831template <int dim>
832QTrapezoid<dim>::QTrapezoid()
833 : Quadrature<dim>(QTrapezoid<dim - 1>(), QTrapezoid<1>())
834{}
835
836
837
838template <int dim>
839QSimpson<dim>::QSimpson()
840 : Quadrature<dim>(QSimpson<dim - 1>(), QSimpson<1>())
841{}
842
843
844
845template <int dim>
846QMilne<dim>::QMilne()
847 : Quadrature<dim>(QMilne<dim - 1>(), QMilne<1>())
848{}
849
850
851template <int dim>
852QWeddle<dim>::QWeddle()
853 : Quadrature<dim>(QWeddle<dim - 1>(), QWeddle<1>())
854{}
855
856template <int dim>
857QTelles<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
873template <int dim>
874QTelles<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
882template <>
883QTelles<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
956namespace 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
1002template <>
1003QGaussChebyshev<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
1018template <int dim>
1019QGaussChebyshev<dim>::QGaussChebyshev(const unsigned int n)
1020 : Quadrature<dim>(QGaussChebyshev<1>(n))
1021{}
1022
1023
1024namespace 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
1098template <>
1099QGaussRadauChebyshev<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
1117template <int dim>
1118QGaussRadauChebyshev<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
1128namespace 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
1170template <>
1171QGaussLobattoChebyshev<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
1190template <int dim>
1191QGaussLobattoChebyshev<dim>::QGaussLobattoChebyshev(const unsigned int n)
1192 : Quadrature<dim>(QGaussLobattoChebyshev<1>(n))
1193{}
1194
1195
1196
1197template <int dim>
1198QSimplex<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 /* Use "int d" instead of the more natural "unsigned int d" to work
1207 * around a wrong diagnostic in gcc-10.3.0 that warns about that the
1208 * comparison "d < dim" is always false in case of "dim == 0".
1209 * MM 2021 */
1210 for (int d = 0; d < dim; ++d)
1211 r += quad.point(i)[d];
1212 if (r <= 1 + 1e-10)
1213 {
1214 this->quadrature_points.push_back(quad.point(i));
1215 this->weights.push_back(quad.weight(i));
1216 }
1217 }
1218}
1219
1220
1221
1222template <int dim>
1223Quadrature<dim>
1224QSimplex<dim>::compute_affine_transformation(
1225 const std::array<Point<dim>, dim + 1> &vertices) const
1226{
1227 Tensor<2, dim> B;
1228 for (unsigned int d = 0; d < dim; ++d)
1229 B[d] = vertices[d + 1] - vertices[0];
1230
1231 B = transpose(B);
1232 const double J = std::abs(determinant(B));
1233
1234 // if the determinant is zero, we return an empty quadrature
1235 if (J < 1e-12)
1236 return Quadrature<dim>();
1237
1238 std::vector<Point<dim>> qp(this->size());
1239 std::vector<double> w(this->size());
1240
1241 for (unsigned int i = 0; i < this->size(); ++i)
1242 {
1243 qp[i] = Point<dim>(vertices[0] + B * this->point(i));
1244 w[i] = J * this->weight(i);
1245 }
1246
1247 return Quadrature<dim>(qp, w);
1248}
1249
1250
1251
1252QTrianglePolar::QTrianglePolar(const Quadrature<1> &radial_quadrature,
1253 const Quadrature<1> &angular_quadrature)
1254 : QSimplex<2>(Quadrature<2>())
1255{
1256 const QAnisotropic<2> base(radial_quadrature, angular_quadrature);
1257 this->quadrature_points.resize(base.size());
1258 this->weights.resize(base.size());
1259 for (unsigned int i = 0; i < base.size(); ++i)
1260 {
1261 const auto &q = base.point(i);
1262 const auto w = base.weight(i);
1263
1264 const auto xhat = q[0];
1265 const auto yhat = q[1];
1266
1267 const double t = numbers::PI_2 * yhat;
1268 const double pi = numbers::PI;
1269 const double st = std::sin(t);
1270 const double ct = std::cos(t);
1271 const double r = xhat / (st + ct);
1272
1273 const double J = pi * xhat / (2 * (std::sin(pi * yhat) + 1));
1274
1275 this->quadrature_points[i] = Point<2>(r * ct, r * st);
1276 this->weights[i] = w * J;
1277 }
1278}
1279
1280
1281
1282QTrianglePolar::QTrianglePolar(const unsigned int n)
1283 : QTrianglePolar(QGauss<1>(n), QGauss<1>(n))
1284{}
1285
1286
1287
1288QDuffy::QDuffy(const Quadrature<1> &radial_quadrature,
1289 const Quadrature<1> &angular_quadrature,
1290 const double beta)
1291 : QSimplex<2>(Quadrature<2>())
1292{
1293 const QAnisotropic<2> base(radial_quadrature, angular_quadrature);
1294 this->quadrature_points.resize(base.size());
1295 this->weights.resize(base.size());
1296 for (unsigned int i = 0; i < base.size(); ++i)
1297 {
1298 const auto &q = base.point(i);
1299 const auto w = base.weight(i);
1300
1301 const auto xhat = q[0];
1302 const auto yhat = q[1];
1303
1304 const double x = std::pow(xhat, beta) * (1 - yhat);
1305 const double y = std::pow(xhat, beta) * yhat;
1306
1307 const double J = beta * std::pow(xhat, 2. * beta - 1.);
1308
1309 this->quadrature_points[i] = Point<2>(x, y);
1310 this->weights[i] = w * J;
1311 }
1312}
1313
1314
1315
1316QDuffy::QDuffy(const unsigned int n, const double beta)
1317 : QDuffy(QGauss<1>(n), QGauss<1>(n), beta)
1318{}
1319
1320
1321
1322template <int dim>
1323QSplit<dim>::QSplit(const QSimplex<dim> &base, const Point<dim> &split_point)
1324{
1325 AssertThrow(GeometryInfo<dim>::is_inside_unit_cell(split_point, 1e-12),
1326 ExcMessage(
1327 "The split point should be inside the unit reference cell."));
1328
1329 std::array<Point<dim>, dim + 1> vertices;
1330 vertices[0] = split_point;
1331
1332 // Make a simplex from the split_point and the first dim vertices of each
1333 // face. In dimension three, we need to split the face in two triangles, so
1334 // we use once the first dim vertices of each face, and the second time the
1335 // the dim vertices of each face starting from 1.
1336 for (auto f : GeometryInfo<dim>::face_indices())
1337 for (unsigned int start = 0; start < (dim > 2 ? 2 : 1); ++start)
1338 {
1339 for (unsigned int i = 0; i < dim; ++i)
1340 vertices[i + 1] = GeometryInfo<dim>::unit_cell_vertex(
1341 GeometryInfo<dim>::face_to_cell_vertices(f, i + start));
1342 const auto quad = base.compute_affine_transformation(vertices);
1343 if (quad.size())
1344 {
1345 this->quadrature_points.insert(this->quadrature_points.end(),
1346 quad.get_points().begin(),
1347 quad.get_points().end());
1348 this->weights.insert(this->weights.end(),
1349 quad.get_weights().begin(),
1350 quad.get_weights().end());
1351 }
1352 }
1353}
1354
1355
1356
1357template <int dim>
1358QGaussSimplex<dim>::QGaussSimplex(const unsigned int n_points_1D)
1359 : QSimplex<dim>(Quadrature<dim>())
1360{
1361 // fill quadrature points and quadrature weights
1362 if (dim == 0 || dim == 1)
1363 {
1364 const ::QGauss<dim> quad(n_points_1D);
1365
1366 this->quadrature_points = quad.get_points();
1367 this->weights = quad.get_weights();
1368 }
1369 else if (dim == 2)
1370 {
1371 if (n_points_1D == 1)
1372 {
1373 const double p = 1.0 / 3.0;
1374 this->quadrature_points.emplace_back(p, p);
1375 this->weights.emplace_back(0.5);
1376 }
1377 else if (n_points_1D == 2)
1378 {
1379 const double Q23 = 2.0 / 3.0;
1380 const double Q16 = 1.0 / 6.0;
1381
1382 this->quadrature_points.emplace_back(Q23, Q16);
1383 this->quadrature_points.emplace_back(Q16, Q23);
1384 this->quadrature_points.emplace_back(Q16, Q16);
1385 this->weights.emplace_back(Q16);
1386 this->weights.emplace_back(Q16);
1387 this->weights.emplace_back(Q16);
1388 }
1389 else if (n_points_1D == 3)
1390 {
1391 const double q12 = 0.5;
1392
1393 // clang-format off
1394 this->quadrature_points.emplace_back(0.3333333333330, 0.3333333333330);
1395 this->quadrature_points.emplace_back(0.7974269853530, 0.1012865073230);
1396 this->quadrature_points.emplace_back(0.1012865073230, 0.7974269853530);
1397 this->quadrature_points.emplace_back(0.1012865073230, 0.1012865073230);
1398 this->quadrature_points.emplace_back(0.0597158717898, 0.4701420641050);
1399 this->quadrature_points.emplace_back(0.4701420641050, 0.0597158717898);
1400 this->quadrature_points.emplace_back(0.4701420641050, 0.4701420641050);
1401 // clang-format on
1402
1403 this->weights.emplace_back(q12 * 0.225);
1404 this->weights.emplace_back(q12 * 0.125939180545);
1405 this->weights.emplace_back(q12 * 0.125939180545);
1406 this->weights.emplace_back(q12 * 0.125939180545);
1407 this->weights.emplace_back(q12 * 0.132394152789);
1408 this->weights.emplace_back(q12 * 0.132394152789);
1409 this->weights.emplace_back(q12 * 0.132394152789);
1410 }
1411 else if (n_points_1D == 4)
1412 {
1413 Quadrature<dim>::operator=(
1414 QWitherdenVincentSimplex<dim>(n_points_1D));
1415 }
1416 }
1417 else if (dim == 3)
1418 {
1419 if (n_points_1D == 1)
1420 {
1421 const double Q14 = 1.0 / 4.0;
1422 const double Q16 = 1.0 / 6.0;
1423
1424 this->quadrature_points.emplace_back(Q14, Q14, Q14);
1425 this->weights.emplace_back(Q16);
1426 }
1427 else if (n_points_1D == 2)
1428 {
1429 const double Q124 = 1.0 / 6.0 / 4.0;
1430
1431 const double palpha = (5.0 + 3.0 * sqrt(5.0)) / 20.0;
1432 const double pbeta = (5.0 - sqrt(5.0)) / 20.0;
1433 this->quadrature_points.emplace_back(pbeta, pbeta, pbeta);
1434 this->quadrature_points.emplace_back(palpha, pbeta, pbeta);
1435 this->quadrature_points.emplace_back(pbeta, palpha, pbeta);
1436 this->quadrature_points.emplace_back(pbeta, pbeta, palpha);
1437 this->weights.emplace_back(Q124);
1438 this->weights.emplace_back(Q124);
1439 this->weights.emplace_back(Q124);
1440 this->weights.emplace_back(Q124);
1441 }
1442 else if (n_points_1D == 3)
1443 {
1444 const double Q16 = 1.0 / 6.0;
1445
1446 // clang-format off
1447 this->quadrature_points.emplace_back(0.5684305841968444, 0.1438564719343852, 0.1438564719343852);
1448 this->quadrature_points.emplace_back(0.1438564719343852, 0.1438564719343852, 0.1438564719343852);
1449 this->quadrature_points.emplace_back(0.1438564719343852, 0.1438564719343852, 0.5684305841968444);
1450 this->quadrature_points.emplace_back(0.1438564719343852, 0.5684305841968444, 0.1438564719343852);
1451 this->quadrature_points.emplace_back(0.0000000000000000, 0.5000000000000000, 0.5000000000000000);
1452 this->quadrature_points.emplace_back(0.5000000000000000, 0.0000000000000000, 0.5000000000000000);
1453 this->quadrature_points.emplace_back(0.5000000000000000, 0.5000000000000000, 0.0000000000000000);
1454 this->quadrature_points.emplace_back(0.5000000000000000, 0.0000000000000000, 0.0000000000000000);
1455 this->quadrature_points.emplace_back(0.0000000000000000, 0.5000000000000000, 0.0000000000000000);
1456 this->quadrature_points.emplace_back(0.0000000000000000, 0.0000000000000000, 0.5000000000000000);
1457 // clang-format on
1458
1459 this->weights.emplace_back(0.2177650698804054 * Q16);
1460 this->weights.emplace_back(0.2177650698804054 * Q16);
1461 this->weights.emplace_back(0.2177650698804054 * Q16);
1462 this->weights.emplace_back(0.2177650698804054 * Q16);
1463 this->weights.emplace_back(0.0214899534130631 * Q16);
1464 this->weights.emplace_back(0.0214899534130631 * Q16);
1465 this->weights.emplace_back(0.0214899534130631 * Q16);
1466 this->weights.emplace_back(0.0214899534130631 * Q16);
1467 this->weights.emplace_back(0.0214899534130631 * Q16);
1468 this->weights.emplace_back(0.0214899534130631 * Q16);
1469 }
1470 else if (n_points_1D == 4)
1471 {
1472 Quadrature<dim>::operator=(
1473 QWitherdenVincentSimplex<dim>(n_points_1D));
1474 }
1475 }
1476
1477 AssertDimension(this->quadrature_points.size(), this->weights.size());
1478 Assert(this->quadrature_points.size() > 0,
1479 ExcNotImplemented(
1480 "QGaussSimplex is currently only implemented for "
1481 "n_points_1D = 1, 2, 3, and 4 while you are asking for "
1482 "n_points_1D = " +
1483 Utilities::to_string(n_points_1D)));
1484}
1485
1486namespace
1487{
1488 template <std::size_t b_dim>
1489 std::vector<std::array<double, b_dim>>
1490 all_permutations(const std::array<double, b_dim> &b_point)
1491 {
1492 std::vector<std::array<double, b_dim>> output;
1493
1494 // We want all possible permutations of the barycentric coordinates.
1495 // The easiest way to get all of them is to sort the input first and
1496 // then use next_permutation to cycle through them all.
1497 std::array<double, b_dim> temp = b_point;
1498 std::sort(temp.begin(), temp.end());
1499 do
1500 {
1501 output.push_back(temp);
1502 }
1503 while (std::next_permutation(temp.begin(), temp.end()));
1504
1505 return output;
1506 }
1507} // namespace
1508
1509
1510
1511template <int dim>
1512QWitherdenVincentSimplex<dim>::QWitherdenVincentSimplex(
1513 const unsigned int n_points_1D)
1514 : QSimplex<dim>(Quadrature<dim>())
1515{
1516 Assert(1 <= dim && dim <= 3, ExcNotImplemented());
1517 // Just use Gauss in 1D: this is a high-order open rule so this is a
1518 // reasonable equivalent for generic programming.
1519 if (dim == 1)
1520 {
1521 Quadrature<dim>::operator=(::QGauss<dim>(n_points_1D));
1522 return;
1523 }
1524
1525 std::array<double, dim + 1> centroid;
1526 std::fill(centroid.begin(), centroid.end(), 1.0 / (dim + 1.0));
1527 std::vector<std::vector<std::array<double, dim + 1>>> b_point_permutations;
1528 std::vector<double> b_weights;
1529
1530 // We can simplify the implementation of these quadrature rules
1531 // by quite a bit by exploiting symmetry - we do essentially the
1532 // same thing for each barycentric coordinate, so we can express
1533 // our quadrature rule as permutations of barycentric points
1534 // instead of writing things out explicitly.
1535
1536 // Apply a Barycentric permutation where one point is different.
1537 auto process_point_1 = [&](const double a, const double w) {
1538 const double b = 1.0 - dim * a;
1539 std::array<double, dim + 1> b_point;
1540 std::fill(b_point.begin(), b_point.begin() + dim, a);
1541 b_point[dim] = b;
1542
1543 b_weights.push_back(w);
1544 b_point_permutations.push_back(all_permutations(b_point));
1545 };
1546
1547 // Apply a Barycentric permutation where two points (in 3D) are different.
1548 auto process_point_2 = [&](const double a, const double w) {
1549 Assert(dim == 3, ExcInternalError());
1550 const double b = (1.0 - 2.0 * a) / 2.0;
1551 std::array<double, dim + 1> b_point;
1552 std::fill(b_point.begin(), b_point.begin() + dim - 1, a);
1553 b_point[dim - 1] = b;
1554 b_point[dim] = b;
1555
1556 b_weights.push_back(w);
1557 b_point_permutations.push_back(all_permutations(b_point));
1558 };
1559
1560 // Apply a Barycentric permutation where three (or four) points
1561 // are different (since there are two inputs).
1562 auto process_point_3 = [&](const double a, const double b, const double w) {
1563 const double c = 1.0 - (dim - 1.0) * a - b;
1564 std::array<double, dim + 1> b_point;
1565 std::fill(b_point.begin(), b_point.begin() + dim - 1, a);
1566 b_point[dim - 1] = b;
1567 b_point[dim] = c;
1568
1569 b_weights.push_back(w);
1570 b_point_permutations.push_back(all_permutations(b_point));
1571 };
1572
1573 if (n_points_1D == 1)
1574 {
1575 b_point_permutations.push_back({centroid});
1576 b_weights.push_back(1.0);
1577 }
1578 else if (n_points_1D == 2)
1579 {
1580 // This is WV-4 in 2D and WV-3 in 3D
1581 if (dim == 2)
1582 {
1583 process_point_1(9.1576213509770743e-02, 1.0995174365532187e-01);
1584 process_point_1(4.4594849091596489e-01, 2.2338158967801147e-01);
1585 }
1586 else if (dim == 3)
1587 {
1588 process_point_1(3.281633025163817e-01, 1.362178425370874e-01);
1589 process_point_1(1.080472498984286e-01, 1.137821574629126e-01);
1590 }
1591 }
1592 else if (n_points_1D == 3)
1593 {
1594 // This is the WV-5 rule in both 2D and 3D
1595 if (dim == 2)
1596 {
1597 b_weights.push_back(0.225);
1598 b_point_permutations.push_back({centroid});
1599
1600 process_point_1(1.0128650732345634e-01, 1.2593918054482714e-01);
1601 process_point_1(4.7014206410511511e-01, 1.3239415278850619e-01);
1602 }
1603 else if (dim == 3)
1604 {
1605 process_point_1(3.108859192633006e-01, 1.126879257180159e-01);
1606 process_point_1(9.273525031089125e-02, 7.349304311636196e-02);
1607
1608 process_point_2(4.550370412564964e-02, 4.254602077708147e-02);
1609 }
1610 }
1611 else if (n_points_1D == 4)
1612 {
1613 // This is the WV-7 rule in both 2D and 3D
1614 if (dim == 2)
1615 {
1616 process_point_1(3.3730648554587850e-02, 1.6545050110792131e-02);
1617 process_point_1(4.7430969250471822e-01, 7.7086646185986069e-02);
1618 process_point_1(2.4157738259540357e-01, 1.2794417123015558e-01);
1619 process_point_3(4.7036644652595216e-02,
1620 1.9868331479735168e-01,
1621 5.5878732903199779e-02);
1622 }
1623 else if (dim == 3)
1624 {
1625 b_point_permutations.push_back({centroid});
1626 b_weights.push_back(9.548528946413085e-02);
1627
1628 process_point_1(3.157011497782028e-01, 4.232958120996703e-02);
1629 process_point_2(5.048982259839635e-02, 3.189692783285758e-02);
1630
1631 process_point_3(1.888338310260010e-01,
1632 5.751716375870000e-01,
1633 3.720713072833462e-02);
1634 process_point_3(2.126547254148314e-02,
1635 8.108302410985486e-01,
1636 8.110770829903342e-03);
1637 }
1638 }
1639 else if (n_points_1D == 5)
1640 {
1641 // This is the WV-9 rule in both 2D and 3D
1642 if (dim == 2)
1643 {
1644 b_point_permutations.push_back({centroid});
1645 b_weights.push_back(9.7135796282798836e-02);
1646
1647 process_point_1(4.4729513394452691e-02, 2.5577675658698031e-02);
1648 process_point_1(4.8968251919873762e-01, 3.1334700227139071e-02);
1649 process_point_1(4.3708959149293664e-01, 7.7827541004774278e-02);
1650 process_point_1(1.8820353561903275e-01, 7.9647738927210249e-02);
1651
1652 process_point_3(3.6838412054736258e-02,
1653 2.2196298916076568e-01,
1654 4.3283539377289376e-02);
1655 }
1656 else if (dim == 3)
1657 {
1658 b_point_permutations.push_back({centroid});
1659 b_weights.push_back(5.801054891248025e-02);
1660
1661 process_point_1(6.198169755222693e-10, 6.431928175925639e-05);
1662 process_point_1(1.607745353952616e-01, 2.317333846242546e-02);
1663 process_point_1(3.222765218214210e-01, 2.956291233542929e-02);
1664 process_point_1(4.510891834541358e-02, 8.063979979616182e-03);
1665
1666 process_point_2(1.122965460043761e-01, 3.813408010370246e-02);
1667
1668 process_point_3(4.588714487524592e-01,
1669 2.554579233041310e-03,
1670 8.384422198298552e-03);
1671 process_point_3(3.377587068533860e-02,
1672 7.183503264420745e-01,
1673 1.023455935274533e-02);
1674 process_point_3(1.836413698099279e-01,
1675 3.441591057817528e-02,
1676 2.052491596798814e-02);
1677 }
1678 }
1679 else if (n_points_1D == 6)
1680 {
1681 // There is no WV-11 rule in 3D yet
1682 if (dim == 2)
1683 {
1684 b_point_permutations.push_back({centroid});
1685 b_weights.push_back(8.5761179732224219e-02);
1686
1687 process_point_1(2.8485417614371900e-02, 1.0431870512894697e-02);
1688 process_point_1(4.9589190096589092e-01, 1.6606273054585369e-02);
1689 process_point_1(1.0263548271224643e-01, 3.8630759237019321e-02);
1690 process_point_1(4.3846592676435220e-01, 6.7316154079468296e-02);
1691 process_point_1(2.1021995670317828e-01, 7.0515684111716576e-02);
1692
1693 process_point_3(7.3254276860644785e-03,
1694 1.4932478865208237e-01,
1695 1.0290289572953278e-02);
1696 process_point_3(4.6010500165429957e-02,
1697 2.8958112563770588e-01,
1698 4.0332476640500554e-02);
1699 }
1700 else if (dim == 3)
1701 {
1702 Assert(false, ExcNotImplemented());
1703 }
1704 }
1705 else
1706 {
1707 Assert(false, ExcNotImplemented());
1708 }
1709
1710 Assert(b_point_permutations.size() == b_weights.size(), ExcInternalError());
1711 for (unsigned int permutation_n = 0; permutation_n < b_weights.size();
1712 ++permutation_n)
1713 {
1714 for (const std::array<double, dim + 1> &b_point :
1715 b_point_permutations[permutation_n])
1716 {
1717 const double volume = (dim == 2 ? 1.0 / 2.0 : 1.0 / 6.0);
1718 this->weights.emplace_back(volume * b_weights[permutation_n]);
1719 Point<dim> c_point;
1720 std::copy(b_point.begin(),
1721 b_point.begin() + dim,
1722 c_point.begin_raw());
1723 this->quadrature_points.emplace_back(c_point);
1724 }
1725 }
1726}
1727
1728
1729
1730template <int dim>
1731QGaussWedge<dim>::QGaussWedge(const unsigned int n_points)
1732 : Quadrature<dim>()
1733{
1734 AssertDimension(dim, 3);
1735
1736 QGaussSimplex<2> quad_tri(n_points);
1737 QGauss<1> quad_line(n_points);
1738
1739 for (unsigned int i = 0; i < quad_line.size(); ++i)
1740 for (unsigned int j = 0; j < quad_tri.size(); ++j)
1741 {
1742 this->quadrature_points.emplace_back(quad_tri.point(j)[0],
1743 quad_tri.point(j)[1],
1744 quad_line.point(i)[0]);
1745 this->weights.emplace_back(quad_tri.weight(j) * quad_line.weight(i));
1746 }
1747
1748 AssertDimension(this->quadrature_points.size(), this->weights.size());
1749 Assert(this->quadrature_points.size() > 0,
1750 ExcMessage("No valid quadrature points!"));
1751}
1752
1753
1754
1755template <int dim>
1756QGaussPyramid<dim>::QGaussPyramid(const unsigned int n_points_1D)
1757 : Quadrature<dim>()
1758{
1759 AssertDimension(dim, 3);
1760
1761 if (n_points_1D == 1)
1762 {
1763 const double Q14 = 1.0 / 4.0;
1764 const double Q43 = 4.0 / 3.0;
1765
1766 this->quadrature_points.emplace_back(0, 0, Q14);
1767 this->weights.emplace_back(Q43);
1768 }
1769 else if (n_points_1D == 2)
1770 {
1771 // clang-format off
1772 this->quadrature_points.emplace_back(-0.26318405556971, -0.26318405556971, 0.54415184401122);
1773 this->quadrature_points.emplace_back(-0.50661630334979, -0.50661630334979, 0.12251482265544);
1774 this->quadrature_points.emplace_back(-0.26318405556971, +0.26318405556971, 0.54415184401122);
1775 this->quadrature_points.emplace_back(-0.50661630334979, +0.50661630334979, 0.12251482265544);
1776 this->quadrature_points.emplace_back(+0.26318405556971, -0.26318405556971, 0.54415184401122);
1777 this->quadrature_points.emplace_back(+0.50661630334979, -0.50661630334979, 0.12251482265544);
1778 this->quadrature_points.emplace_back(+0.26318405556971, +0.26318405556971, 0.54415184401122);
1779 this->quadrature_points.emplace_back(+0.50661630334979, +0.50661630334979, 0.12251482265544);
1780 // clang-format on
1781
1782 this->weights.emplace_back(0.10078588207983);
1783 this->weights.emplace_back(0.23254745125351);
1784 this->weights.emplace_back(0.10078588207983);
1785 this->weights.emplace_back(0.23254745125351);
1786 this->weights.emplace_back(0.10078588207983);
1787 this->weights.emplace_back(0.23254745125351);
1788 this->weights.emplace_back(0.10078588207983);
1789 this->weights.emplace_back(0.23254745125351);
1790 }
1791
1792 AssertDimension(this->quadrature_points.size(), this->weights.size());
1793 Assert(this->quadrature_points.size() > 0,
1794 ExcMessage("No valid quadrature points!"));
1795}
1796
1797
1798
1799// explicit specialization
1800// note that 1d formulae are specialized by implementation above
1801template class QGauss<2>;
1802template class QGaussLobatto<2>;
1803template class QMidpoint<2>;
1804template class QTrapezoid<2>;
1805template class QSimpson<2>;
1806template class QMilne<2>;
1807template class QWeddle<2>;
1808
1809template class QGauss<3>;
1810template class QGaussLobatto<3>;
1811template class QMidpoint<3>;
1812template class QTrapezoid<3>;
1813template class QSimpson<3>;
1814template class QMilne<3>;
1815template class QWeddle<3>;
1816
1817template class QSorted<1>;
1818template class QSorted<2>;
1819template class QSorted<3>;
1820
1821template class QTelles<1>;
1822template class QTelles<2>;
1823template class QTelles<3>;
1824
1825template class QGaussChebyshev<1>;
1826template class QGaussChebyshev<2>;
1827template class QGaussChebyshev<3>;
1828
1829template class QGaussRadauChebyshev<1>;
1830template class QGaussRadauChebyshev<2>;
1831template class QGaussRadauChebyshev<3>;
1832
1833template class QGaussLobattoChebyshev<1>;
1834template class QGaussLobattoChebyshev<2>;
1835template class QGaussLobattoChebyshev<3>;
1836
1837template class QSimplex<1>;
1838template class QSimplex<2>;
1839template class QSimplex<3>;
1840
1841template class QSplit<1>;
1842template class QSplit<2>;
1843template class QSplit<3>;
1844
1845template class QGaussSimplex<0>;
1846template class QGaussSimplex<1>;
1847template class QGaussSimplex<2>;
1848template class QGaussSimplex<3>;
1849template class QGaussWedge<1>;
1850template class QGaussWedge<2>;
1851template class QGaussWedge<3>;
1852template class QGaussPyramid<1>;
1853template class QGaussPyramid<2>;
1854template class QGaussPyramid<3>;
1855
1856template class QWitherdenVincentSimplex<1>;
1857template class QWitherdenVincentSimplex<2>;
1858template class QWitherdenVincentSimplex<3>;
1859
1860DEAL_II_NAMESPACE_CLOSE
DerivativeForm::transpose
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Definition: derivative_form.h:536
QDuffy::QDuffy
QDuffy(const Quadrature< 1 > &radial_quadrature, const Quadrature< 1 > &angular_quadrature, const double beta=1.0)
Definition: quadrature_lib.cc:1288
QGaussChebyshev::QGaussChebyshev
QGaussChebyshev(const unsigned int n)
Generate a formula with n quadrature points.
Definition: quadrature_lib.cc:1019
QGaussLobattoChebyshev::QGaussLobattoChebyshev
QGaussLobattoChebyshev(const unsigned int n)
Generate a formula with n quadrature points.
Definition: quadrature_lib.cc:1191
QGaussLobatto::QGaussLobatto
QGaussLobatto(const unsigned int n)
Definition: quadrature_lib.cc:818
QGaussLogR::fraction
const double fraction
Definition: quadrature_lib.h:272
QGaussLogR::QGaussLogR
QGaussLogR(const unsigned int n, const Point< dim > &x0=Point< dim >(), const double alpha=1, const bool factor_out_singular_weight=false)
QGaussLog::QGaussLog
QGaussLog(const unsigned int n, const bool revert=false)
QGaussLog::get_quadrature_points
static std::vector< double > get_quadrature_points(const unsigned int n)
QGaussLog::get_quadrature_weights
static std::vector< double > get_quadrature_weights(const unsigned int n)
QGaussOneOverR::QGaussOneOverR
QGaussOneOverR(const unsigned int n, const Point< dim > &singularity, const bool factor_out_singular_weight=false)
QGaussOneOverR::quad_size
static unsigned int quad_size(const Point< dim > &singularity, const unsigned int n)
QGaussPyramid::QGaussPyramid
QGaussPyramid(const unsigned int n_points_1D)
Definition: quadrature_lib.cc:1756
QGaussRadauChebyshev::QGaussRadauChebyshev
QGaussRadauChebyshev(const unsigned int n, EndPoint ep=QGaussRadauChebyshev::left)
Generate a formula with n quadrature points.
Definition: quadrature_lib.cc:1118
QGaussSimplex::QGaussSimplex
QGaussSimplex(const unsigned int n_points_1D)
Definition: quadrature_lib.cc:1358
QGaussWedge::QGaussWedge
QGaussWedge(const unsigned int n_points_1D)
Definition: quadrature_lib.cc:1731
QGauss::QGauss
QGauss(const unsigned int n)
Definition: quadrature_lib.cc:811
QSimplex::QSimplex
QSimplex(const Quadrature< dim > &quad)
Definition: quadrature_lib.cc:1198
QSimplex::compute_affine_transformation
Quadrature< dim > compute_affine_transformation(const std::array< Point< dim >, dim+1 > &vertices) const
Definition: quadrature_lib.cc:1224
QSorted::QSorted
QSorted(const Quadrature< dim > &quad)
Definition: quadrature_lib.cc:769
QSorted::compare_weights
bool compare_weights(const unsigned int a, const unsigned int b) const
Definition: quadrature_lib.cc:801
QSplit::QSplit
QSplit(const QSimplex< dim > &base, const Point< dim > &split_point)
Definition: quadrature_lib.cc:1323
QTelles::QTelles
QTelles(const Quadrature< 1 > &base_quad, const Point< dim > &singularity)
Definition: quadrature_lib.cc:857
QTrianglePolar::QTrianglePolar
QTrianglePolar(const Quadrature< 1 > &radial_quadrature, const Quadrature< 1 > &angular_quadrature)
Definition: quadrature_lib.cc:1252
QWitherdenVincentSimplex::QWitherdenVincentSimplex
QWitherdenVincentSimplex(const unsigned int n_points_1D)
Definition: quadrature_lib.cc:1512
Quadrature::quadrature_points
std::vector< Point< dim > > quadrature_points
Definition: quadrature.h:283
Quadrature::operator=
Quadrature & operator=(const Quadrature< dim > &)
Definition: quadrature.cc:282
Quadrature::point
const Point< dim > & point(const unsigned int i) const
Quadrature::is_tensor_product_flag
bool is_tensor_product_flag
Definition: quadrature.h:298
Quadrature::weight
double weight(const unsigned int i) const
Quadrature::weights
std::vector< double > weights
Definition: quadrature.h:289
Quadrature::size
unsigned int size() const
SymmetricTensor::determinant
constexpr Number determinant(const SymmetricTensor< 2, dim, Number > &t)
Definition: symmetric_tensor.h:2622
Tensor
Definition: tensor.h:472
Tensor::begin_raw
Number * begin_raw()
Tensor::end_raw
Number * end_raw()
DEAL_II_NAMESPACE_OPEN
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:402
DEAL_II_NAMESPACE_CLOSE
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:403
vertices
Point< 3 > vertices[4]
Definition: data_out_base.cc:175
geometry_info.h
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1465
AssertDimension
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:1622
DataOutBase::eps
@ eps
Definition: data_out_base.h:1580
AssertIndexRange
#define AssertIndexRange(index, range)
Definition: exceptions.h:1690
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
AssertThrow
#define AssertThrow(cond, exc)
Definition: exceptions.h:1575
GridTools::volume
double volume(const Triangulation< dim, spacedim > &tria, const Mapping< dim, spacedim > &mapping=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()))
Definition: grid_tools.cc:137
LocalIntegrators::Divergence::norm
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Definition: divergence.h:472
OpenCASCADE::point
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:188
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=(ReferenceCells::get_hypercube< dim >() .template get_default_linear_mapping< dim, spacedim >()), const std::vector< std::vector< double > > &properties={})
Definition: generators.cc:451
Physics::Elasticity::Kinematics::w
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
Physics::Elasticity::Kinematics::b
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Polynomials::jacobi_polynomial_value
Number jacobi_polynomial_value(const unsigned int degree, const int alpha, const int beta, const Number x)
Definition: polynomial.h:976
Utilities::to_string
std::string to_string(const number value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition: utilities.cc:482
internal::QGaussChebyshev::get_quadrature_weights
std::vector< double > get_quadrature_weights(const unsigned int n)
Definition: quadrature_lib.cc:986
internal::QGaussChebyshev::get_quadrature_points
std::vector< double > get_quadrature_points(const unsigned int n)
Definition: quadrature_lib.cc:964
internal::QGaussLobattoChebyshev::get_quadrature_points
std::vector< double > get_quadrature_points(const unsigned int n)
Definition: quadrature_lib.cc:1134
internal::QGaussLobattoChebyshev::get_quadrature_weights
std::vector< double > get_quadrature_weights(const unsigned int n)
Definition: quadrature_lib.cc:1151
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
internal::QGaussLobatto::gamma
long double gamma(const unsigned int n)
Definition: quadrature_lib.cc:96
internal::QGaussRadauChebyshev::get_quadrature_points
std::vector< double > get_quadrature_points(const unsigned int n, ::QGaussRadauChebyshev< 1 >::EndPoint ep)
Definition: quadrature_lib.cc:1030
internal::QGaussRadauChebyshev::get_quadrature_weights
std::vector< double > get_quadrature_weights(const unsigned int n, ::QGaussRadauChebyshev< 1 >::EndPoint ep)
Definition: quadrature_lib.cc:1076
internal::TensorProductManifoldImplementation::split_point
void split_point(const Point< dim1+dim2 > &source, Point< dim1 > &p1, Point< dim2 > &p2)
Definition: tensor_product_manifold.h:160
internal::VectorOperations::copy
void copy(const T *begin, const T *end, U *dest)
Definition: vector_operations_internal.h:68
numbers::PI_2
static constexpr double PI_2
Definition: numbers.h:236
numbers::PI
static constexpr double PI
Definition: numbers.h:231
std::log
::VectorizedArray< Number, width > log(const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5384
std::tan
::VectorizedArray< Number, width > tan(const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5340
std::cos
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5318
std::sin
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5291
std::sqrt
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5406
std::pow
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
Definition: vectorization.h:5422
std::abs
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5468
quadrature_lib.h
GeometryInfo::unit_cell_vertex
static Point< dim > unit_cell_vertex(const unsigned int vertex)
Generated by doxygen 1.9.4