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