NeoPZ
tpzgaussrule.cpp
Go to the documentation of this file.
1 
5 #include "pz_config.h"
6 #include <math.h>
7 #include <cmath>
8 
9 #include "tpzgaussrule.h"
10 #include "tpzintrulelist.h"
11 #include "pzerror.h"
12 
13 #include "pzvec.h"
14 
15 #ifdef VC
16 #include <algorithm>
17 #endif
18 
19 using namespace std;
20 
21 //***************************************
22 TPZGaussRule::TPZGaussRule(int order,int type,long double alpha,long double beta) {
23  // Cleaning storage
24  fLocation.Resize(0);
25  fWeight.Resize(0);
26 
27  if(alpha <= -1.0L) alpha = 0.0L; // it is required to right computation
28  if(beta <= -1.0L) beta = 0.0L; // it is required to right computation
29  fAlpha = alpha;
30  fBeta = beta;
31 
32  fType = type;
33 
34  if(order < 0) {
35  PZError << "TPZGaussRule criation, order = " << order << " not allowable\n";
36  fNumInt = 0;
37  fLocation.Resize(0);
38  fWeight.Resize(0);
39  return;
40  }
41  switch (fType) {
42  case 1: // Lobatto
43  {
44  // Computing the points and weights for the symmetric Gaussian quadrature (using Legendre polynomials)
45  fOrder = ComputingGaussLobattoQuadrature(order);
46  }
47  break;
48  case 2: // Jacobi
49  if(IsZero(fAlpha) && IsZero(fBeta)) fAlpha = fBeta = 1.0L; // Because fAlpha = fBeta = 0.0 is same as Gauss Legendre case
50  // Computing the points and weights for the symmetric Gaussian quadrature (using Legendre polynomials)
51  fOrder = ComputingGaussJacobiQuadrature(order,fAlpha,fBeta);
52  break;
53  case 3: // Chebyshev
54  fAlpha = fBeta = -0.5L;
55  fOrder = ComputingGaussChebyshevQuadrature(order);
56  break;
57  default: // Legendre (default)
58  // Computing the points and weights for the symmetric Gaussian quadrature (using Legendre polynomials)
59  fOrder = ComputingGaussLegendreQuadrature(order);
60  break;
61  }
62  if(fOrder<0)
63  fOrder = order;
64 }
65 
66 TPZGaussRule::TPZGaussRule(const TPZGaussRule &copy) : fType(copy.fType), fNumInt(copy.fNumInt), fLocation(copy.fLocation),
67  fWeight(copy.fWeight), fAlpha(copy.fAlpha), fBeta(copy.fBeta), fOrder(copy.fOrder)
68 {
69 
70 }
71 
72 TPZGaussRule &TPZGaussRule::operator=(const TPZGaussRule &copy)
73 {
74  fType = copy.fType;
75  fNumInt = copy.fNumInt;
76  fLocation = copy.fLocation;
77  fWeight = copy.fWeight;
78  fAlpha = copy.fAlpha;
79  fBeta = copy.fBeta;
80  fOrder = copy.fOrder;
81  return *this;
82 }
83 
84 
85 
86 //***************************************
87 //***************************************
88 TPZGaussRule::~TPZGaussRule() {
89  fNumInt = 0;
90  fLocation.Resize(0);
91  fWeight.Resize(0);
92 }
93 
94 //***************************************
95 //***************************************
96 long double TPZGaussRule::Loc(int i) const {
97  if (i>=0 && i<fNumInt)
98  return fLocation[i];
99  else {
100  PZError << "ERROR(TPZGaussRule::loc) Out of bounds!!\n";
101  return 0.0L;
102  }
103 }
104 
105 //***************************************
106 //***************************************
107 long double TPZGaussRule::W(int i) const {
108  if (i>=0 && i<fNumInt)
109  return fWeight[i];
110  else {
111  PZError << "ERROR(TPZGaussRule::w) Out of bounds!!\n";
112  return 0.0L;
113  }
114 }
115 
116 void TPZGaussRule::Print(std::ostream &out) {
117  int i;
118  out << "\nGaussian Quadrature ";
119  if(!fType) out << "(Legendre) : " << std::endl;
120  else if(fType == 1) out << "(Lobatto) : " << std::endl;
121  else if(fType == 2) out << "(Jacobi) : " << std::endl;
122  else if(fType == 3) out << "(Chebyshev) : " << std::endl;
123  else out << "(Unkown) : " << std::endl;
124  out << "\nNumber of points : " << fNumInt << std::endl;
125  for(i=0;i<fNumInt;i++) {
126  out << i << " :\t" << fLocation[i] << " \t" << fWeight[i] << std::endl;
127  }
128  out << std::endl;
129 }
130 
131 void TPZGaussRule::SetType(int &type,int order) {
132  if(order < 0) order = 1;
133  if(type < 0 || type > 3) type = 0;
134  fLocation.Resize(0);
135  fWeight.Resize(0);
136  fType = type;
137  int neworder;
138  switch(type) {
139  case 1: // Lobatto
140  {
141  // Computing the points and weights for the symmetric Gaussian quadrature (using Legendre polynomials)
142  neworder=ComputingGaussLobattoQuadrature(order);
143  }
144  break;
145  case 2: // Jacobi
146  {
147  // Computing the points and weights for the symmetric Gaussian quadrature (using Legendre polynomials)
148  neworder=ComputingGaussJacobiQuadrature(order,fAlpha,fBeta);
149  }
150  break;
151  case 3: // Chebyshev
152  {
153  neworder=ComputingGaussChebyshevQuadrature(order);
154  }
155  break;
156  default: // Legendre (default)
157  {
158  // Computing the points and weights for the symmetric Gaussian quadrature (using Legendre polynomials)
159  neworder=ComputingGaussLegendreQuadrature(order);
160  }
161  break;
162  }
163  if(neworder < 0)
164  DebugStop();
165 }
166 
169 // Compute the points and weights for Gauss Legendre Quadrature over the parametric 1D element [-1.0, 1.0] - This quadrature is symmetric
170 int TPZGaussRule::ComputingGaussLegendreQuadrature(int order) {
171  // Computing number of points appropriated to the wished order = 2*npoints - 1. Note: If odd we need increment one point more
172  fNumInt = (int)(0.51*(order+2));
173  if(fNumInt < 1)
174  fNumInt = 1;
175  if(fNumInt > TPZGaussRule::NRULESLEGENDRE_ORDER)
176  fNumInt = TPZGaussRule::NRULESLEGENDRE_ORDER-1;
177 
178  order = (fNumInt-1)*2;
179  ComputingGaussLegendreQuadrature(&fNumInt,fLocation,fWeight);
180  return order;
181 }
182 // Compute the points and weights for Gauss Legendre Quadrature over the parametric 1D element [-1.0, 1.0] - This quadrature is symmetric
183 void TPZGaussRule::ComputingGaussLegendreQuadrature(int *npoints,TPZVec<long double> &Location,TPZVec<long double> &Weight) {
184  long double tol = machinePrecision();
185  long double z1, z, pp, p3, p2, p1, dif, den;
186  int i, j;
187  int64_t iteration;
188 
189  // Cleaning vector to storage
190  Location.Resize(0);
191  Weight.Resize(0);
192  Location.Resize((*npoints),0.0L);
193  Weight.Resize((*npoints),0.0L);
194 
195  int m = ((*npoints)+1)/2;
196  long double weight;
197 
198  for(i=0;i<m;i++) {
199  p1 = ((long double)(i))+(0.75L);
200  p2 = ((long double)((*npoints)))+0.5L;
201  z = cosl((M_PI*p1)/p2); // p2 never is zero
202  iteration = 0L;
203  do {
204  iteration++;
205  p1 = 1.0L;
206  p2 = 0.0L;
207  for(j=0;j<(*npoints);j++) {
208  p3 = p2;
209  p2 = p1;
210  p1 = (((2.0L)*((long double)(j))+(1.0L))*z*p2-(((long double)(j))*p3))/(((long double)(j))+(1.0L)); // denominator never will be zero
211  }
212  den = (z*z)-(1.0L);
213  if(fabs(den)<1.e-16L)
214  z = 0.5L;
215  pp = ((long double)(*npoints))*(z*p1-p2)/den;
216  z1 = z;
217  if(fabs(pp)<1.e-16L)
218  z = 0.5L;
219  else z = z1-p1/pp;
220  dif = fabs(z - z1);
221  }while(fabs(dif) > tol && iteration < PZINTEGRAL_MAXITERATIONS_ALLOWED);
222 
223  // If the maxime number of iterations is reached then the computing is stopped
224  if(iteration == PZINTEGRAL_MAXITERATIONS_ALLOWED) {
225  std::cout << "Reached maxime number of iterations for NPOINTS = " << (*npoints) << ".\n";
226  (*npoints) = 0;
227  Location.Resize(0);
228  Weight.Resize(0);
229  break;
230  }
231  // Philippe Order
232  weight = 2.0L/((1.0L-z*z)*pp*pp);
233  Location[2*i] = -z;
234  Weight[2*i] = weight;
235  if((2*i+1)<(*npoints)) {
236  Location[2*i+1] = z;
237  Weight[2*i+1] = weight;
238  }
239  }
240 }
241 
242 int TPZGaussRule::ComputingGaussChebyshevQuadrature(int order) {
243  // Computing number of points appropriated to the wished order = 2*npoints - 1. Note: If odd we need increment one point more
244  fNumInt = (int)(0.51*(order+2));
245  if(fNumInt < 0)
246  fNumInt = 1;
247 
248  // Cleaning vector to storage
249  fLocation.Resize(0);
250  fWeight.Resize(0);
251  fLocation.Resize(fNumInt, 0.0L);
252  fWeight.Resize(fNumInt, 0.0L);
253 
254  for(int i=0;i<fNumInt;i++) {
255  fLocation[i] = cosl(M_PI*(long double)(i+0.5L)/((long double)(fNumInt)));
256  fWeight[i] = M_PI/((long double)(fNumInt));
257  }
258  return order;
259 }
260 
261 // To compute the polinomial Jacobi on x point, n is the order, alpha and beta are the exponents as knowed
262 long double TPZGaussRule::JacobiPolinomial(long double x,int alpha,int beta,unsigned int n) {
263  // the Jacobi polynomial is evaluated
264  // using a recursion formula.
265  TPZVec<long double> p(n+1);
266  int v, a1, a2, a3, a4;
267 
268  // initial values P_0(x), P_1(x):
269  p[0] = 1.0L;
270  if(n==0) return p[0];
271  p[1] = (((long double)(alpha+beta+2.0L))*x + ((long double)(alpha-beta)))/2.0L;
272  if(n==1) return p[1];
273 
274  for(unsigned int i=1; i<=(n-1); ++i) {
275  v = 2*i + alpha + beta;
276  a1 = 2*(i+1)*(i + alpha + beta + 1)*v;
277  a2 = (v + 1)*(alpha*alpha - beta*beta);
278  a3 = v*(v + 1)*(v + 2);
279  a4 = 2*(i+alpha)*(i+beta)*(v + 2);
280 
281  p[i+1] = static_cast<long double>( (((long double)(a2)) + ((long double)(a3))*x)*p[i] - ((long double)(a4))*p[i-1])/((long double)(a1));
282  } // for
283  return p[n];
284 }
285 
286 using namespace std;
288 int TPZGaussRule::ComputingGaussLobattoQuadrature(int order) {
289  // Computing number of points appropriated to the wished order = 2*npoints - 1. Note: If odd we need increment one point more
290  fNumInt = (int)(0.51*(order+2));
291  fNumInt++;
292 
293  if(fNumInt < 2)
294  fNumInt = 2;
295  if(fNumInt > NRULESLOBATTO_ORDER)
296  fNumInt = NRULESLOBATTO_ORDER-1;
297 
298  order = 2*(fNumInt-1);
299  // Cleaning vector to storage
300  fLocation.Resize(0);
301  fWeight.Resize(0);
302  fLocation.Resize(fNumInt, 0.0L);
303  fWeight.Resize(fNumInt, 0.0L);
304 
305  // Process to compute inner integration points at [-1.,1.]
306  const unsigned int m = fNumInt - 2; // no. of inner points
307 
308  unsigned int counter = 0;
309  // compute quadrature points with a Newton algorithm.
310 
311  // Set tolerance.
312  const long double eps = 1.09e-19L;
313  long double deps = 2.23e-16L;
314 
315  // check whether long double is more accurate than double, and
316  // set tolerances accordingly
317  const long double tol = ((1.L + eps) != 1.L ? max(deps/(100.0L), eps*5.0L) : deps*5.0L);
318 
319  // it take the zeros of the Chebyshev polynomial (alpha=beta=-0.5) as initial values:
320  unsigned int i, k;
321  for(i=0; i<m; ++i)
322  fLocation[i] = - cosl(((2.0L*((long double)(i))+1.0L)/(2.0L*((long double)(m))))*M_PI);
323 
324  long double r, s, J_x, f, delta;
325 
326  for(k=0; k<m; ++k) {
327  r = fLocation[k];
328  if(k>0) r = (r + fLocation[k-1])/2.0L;
329  do {
330  counter++;
331  s = 0.0L;
332  for (i=0; i<k; ++i)
333  s += 1.0L/(r - fLocation[i]);
334 
335  J_x = 0.5L*((long double)(m + 3.0L))*JacobiPolinomial(r,2,2, m-1);
336  f = JacobiPolinomial(r,1,1, m);
337  delta = f/(f*s- J_x);
338  r += delta;
339  }
340  while (fabs(delta) >= tol && counter < PZINTEGRAL_MAXITERATIONS_ALLOWED);
341 
342  // If maxime iterations is reached, clean all points
343  if(counter == PZINTEGRAL_MAXITERATIONS_ALLOWED) {
344  fNumInt = 0;
345  fWeight.Resize(0);
346  fLocation.Resize(0);
347  cout << "Maxime number of iterations allowed in quadrature GaussLobatto. \n" ;
348  return -1; // to compile in linux
349  }
350 
351  fLocation[k] = r;
352  }
353 
354  int ii;
355  for(ii=m-1;ii>-1;ii--)
356  fLocation[ii+1] = fLocation[ii];
357 
358  fLocation[0] = -1.0L;
359  fLocation[fNumInt-1] = 1.0L;
360 
361  // Process to compute the weights to corresponding integration points
362  s = 0.L;
363  i = (unsigned int) fNumInt;
364  f = gamma(i);
365  r = gamma(i+1);
366  const long double factor = 2.L*f*f / (((long double)(fNumInt-1))*f*r);
367 
368  fWeight.Resize(fNumInt,0.0L);
369 
370  for(ii=0; ii < fNumInt; ++ii)
371  {
372  s = JacobiPolinomial(fLocation[ii],0,0, fNumInt-1);
373  fWeight[ii] = factor/(s*s);
374  }
375  return order;
376 }
377 
379 int TPZGaussRule::ComputingGaussJacobiQuadrature(int order,long double alpha, long double beta) {
380  // Computing number of points appropriated to the wished order = 2*npoints - 1. Note: If odd we return -1;need increment one point more
381  fNumInt = (int)(0.51*(order+2));
382  if(fNumInt < 0)
383  fNumInt = 1;
384 
385  ComputingGaussJacobiQuadrature(&fNumInt,alpha,beta,fLocation,fWeight);
386  return order;
387 }
388 
390 void TPZGaussRule::ComputingGaussJacobiQuadrature(int *npoints,long double alpha, long double beta,TPZVec<long double> &Location,TPZVec<long double> &Weight) {
391  // Computing number of points appropriated to the wished order = 2*npoints - 1. Note: If odd we need increment one point more
392  long double an, bn;
393  TPZVec<long double> b;
394  TPZVec<long double> c;
395  long double cc, delta, dp2;
396  int i;
397  long double p1, prod, r1, r2, r3, temp, x0 = 0.0;
398 
399  Location.Resize((*npoints),0.0L);
400  Weight.Resize((*npoints),0.0L);
401 
402  b.Resize((*npoints),0.0L);
403  c.Resize((*npoints),0.0L);
404 
405  // This method permit only alpha > -1.0 and beta > -1.0
406  if(alpha <= -1.0L) {
407  std::cerr << "\nJACOBI_COMPUTE - Fatal error!\n -1.0 < ALPHA is required.\n";
408  std::exit ( 1 );
409  }
410  if(beta <= -1.0L) {
411  std::cerr << "\nJACOBI_COMPUTE - Fatal error!\n -1.0 < BETA is required.\n";
412  std::exit ( 1 );
413  }
414 
415  // Set the recursion coefficients.
416  for(i=1;i<=(*npoints);i++) {
417  if(alpha + beta == 0.0L || beta - alpha == 0.0L) {
418  b[i-1] = 0.0L;
419  }
420  else {
421  b[i-1] = ( alpha + beta ) * ( beta - alpha ) /
422  ( ( alpha + beta + ( 2.0L * ((long double)(i)) ) )
423  * ( alpha + beta + ( 2.0L * ((long double)(i)) - 2.0L ) ) );
424  }
425  if(i==1) {
426  c[i-1] = 0.0L;
427  }
428  else {
429  c[i-1] = 4.0L * ( ((long double)(i)) - 1.0L ) * ( alpha + ( ((long double)(i)) - 1.0L ) ) * ( beta + ( ((long double)(i)) - 1.0L ) )
430  * ( alpha + beta + ( ((long double)(i)) - 1.0L ) ) /
431  ( ( alpha + beta + ( 2.0L * ((long double)(i)) - 1.0L ) ) * (std::pow(alpha + beta + (2.0L*((long double)(i)) - 2.0L),2.0L))
432  * ( alpha + beta + ( 2.0L * ((long double)(i)) - 3.0L ) ) );
433  }
434  }
435 
436  delta = (gamma(alpha + 1.0L) * gamma(beta + 1.0L)) / gamma(alpha + beta + 2.0L);
437 
438  prod = 1.0L;
439  for(i=2;i <= (*npoints);i++)
440  prod = prod * c[i-1];
441  cc = delta * std::pow( 2.0L, alpha + beta + 1.0L ) * prod;
442 
443  for(i=1;i<=(*npoints);i++) {
444  if(i==1) {
445  an = alpha / (long double)((*npoints));
446  bn = beta / (long double)((*npoints));
447 
448  r1 = ( 1.0L + alpha ) * ( 2.78L / ( 4.0L + (long double)((*npoints)*(*npoints))) + 0.768L * an / (long double)((*npoints)));
449  r2 = 1.0L + 1.48L * an + 0.96L * bn + 0.452L * an * an + 0.83L * an * bn;
450 
451  x0 = ( r2 - r1 ) / r2;
452  }
453  else if(i==2) {
454  r1 = ( 4.1L + alpha ) / ( ( 1.0L + alpha ) * ( 1.0L + 0.156L * alpha ) );
455 
456  r2 = 1.0L + 0.06L * ((long double)((*npoints)) - 8.0L) * (1.0L + 0.12L * alpha) / (long double)((*npoints));
457  r3 = 1.0L + 0.012L * beta * (1.0L + 0.25L * fabsl(alpha)) / (long double)((*npoints));
458 
459  x0 = x0 - r1*r2*r3*(1.0L - x0);
460  }
461  else if(i==3) {
462  r1 = ( 1.67L + 0.28L * alpha ) / ( 1.0L + 0.37L * alpha );
463  r2 = 1.0L + 0.22L * ( (long double)((*npoints)) - 8.0L) / (long double)((*npoints));
464  r3 = 1.0L + 8.0L * beta / ( ( 6.28L + beta ) * (long double)((*npoints)*(*npoints)));
465 
466  x0 = x0 - r1 * r2 * r3 * ( Location[0] - x0 );
467  }
468  else if(i<(*npoints)-1) {
469  x0 = 3.0L * Location[i-2] - 3.0L * Location[i-3] + Location[i-4];
470  }
471  else if(i==(*npoints)-1) {
472  r1 = ( 1.0L + 0.235L * beta ) / (0.766L + 0.119L * beta);
473  r2 = 1.0L / ( 1.0L + 0.639L * ((long double)((*npoints)) - 4.0L)
474  / ( 1.0L + 0.71L * ((long double)((*npoints)) - 4.0L)));
475 
476  r3 = 1.0L / ( 1.0L + 20.0L * alpha / ( ( 7.5L + alpha ) *
477  (long double)((*npoints)*(*npoints))));
478 
479  x0 = x0 + r1 * r2 * r3 * ( x0 - Location[i-3] );
480  }
481  else if(i==(*npoints)) {
482  r1 = ( 1.0L + 0.37L * beta ) / ( 1.67L + 0.28L * beta );
483 
484  r2 = 1.0L / ( 1.0L + 0.22L * ( (long double)((*npoints)) - 8.0L )
485  / (long double)((*npoints)));
486 
487  r3 = 1.0L / ( 1.0L + 8.0L * alpha /
488  ( ( 6.28L + alpha ) * (long double)((*npoints)*(*npoints))));
489 
490  x0 = x0 + r1 * r2 * r3 * ( x0 - Location[i-3] );
491  }
492 
493  // Improving the root of the Jacobi polinomial
494  long double d;
495  long double precision;
496  int itera, itera_max = 10;
497 
498  precision = machinePrecision();
499 
500  for(itera=1;itera<=itera_max;itera++) {
501  d = JacobiPolinomial(x0,(*npoints),alpha,beta,b,c,&dp2,&p1);
502  d /= dp2;
503  x0 = x0 - d;
504  if(fabsl(d) <= precision*(fabsl(x0) + 1.0L)) {
505  break;
506  }
507  }
508  Location[i-1] = x0;
509  Weight[i-1] = cc/(dp2*p1);
510  }
511 
512  // inverting the position of the points
513  for(i=1;i<=(*npoints)/2;i++) {
514  temp = Location[i-1];
515  Location[i-1] = Location[(*npoints)-i];
516  Location[(*npoints)-i] = temp;
517  }
518  // inverting the position of the weights
519  for(i=1;i<=(*npoints)/2;i++) {
520  temp = Weight[i-1];
521  Weight[i-1] = Weight[(*npoints)-i];
522  Weight[(*npoints)-i] = temp;
523  }
524 }
525 
527 long double TPZGaussRule::JacobiPolinomial(long double x, int order,long double alpha, long double beta,
528  TPZVec<long double> &b, TPZVec<long double> &c,long double *dp2, long double *p1 )
529 {
530  long double p2, p0, dp0, dp1;
531  int i;
532 
533  *p1 = 1.0L;
534  dp1 = 0.0L;
535  p2 = x + (alpha - beta)/(alpha + beta + 2.0L);
536  *dp2 = 1.0L;
537 
538  for(i=2;i<=order;i++) {
539  p0 = *p1;
540  dp0 = dp1;
541 
542  *p1 = p2;
543  dp1 = *dp2;
544 
545  p2 = (x - b[i-1])*(*p1) - c[i-1]*p0;
546  *dp2 = (x - b[i-1])*dp1 + (*p1) - c[i-1]*dp0;
547  }
548  return p2;
549 }
550 
552 long double machinePrecision() {
553  long double precision = 1.0L;
554  int iterations = 0;
555 
556  while ( 1.0L < (long double) (1.0L + precision) && iterations < PZINTEGRAL_MAXITERATIONS_ALLOWED) {
557  precision = precision / 2.0L;
558  if(iterations == PZINTEGRAL_MAXITERATIONS_ALLOWED) {
559  cout << "\nMaxime iterations reached - epsilon() \n";
560  return 2.0*precision;
561  }
562  }
563  return 2.0L*precision;
564 }
565 
567 long double gamma(unsigned int n) {
568  long double fatorial = (long double)(n - 1);
569  for (int i=n-2; i>1; --i)
570  fatorial *= (long double)i;
571  return fatorial;
572 }
573 
575 long double gamma(long double x) {
576  // Coefficients for minimax approximation over (see reference).
577  long double c[7] = { -1.910444077728E-03L, 8.4171387781295E-04L, -5.952379913043012E-04L,
578  7.93650793500350248E-04L, -2.777777777777681622553E-03L, 8.333333333333333331554247E-02L,
579  5.7083835261E-03L };
580  long double precision = 2.22E-16L;
581  long double fact, half = 0.5L;
582  int i, n;
583  long double p[8] = { -1.71618513886549492533811E+00L, 2.47656508055759199108314E+01L,
584  -3.79804256470945635097577E+02L, 6.29331155312818442661052E+02L, 8.66966202790413211295064E+02L,
585  -3.14512729688483675254357E+04L, -3.61444134186911729807069E+04L, 6.64561438202405440627855E+04L };
586  bool parity;
587  long double q[8] = { -3.08402300119738975254353E+01L, 3.15350626979604161529144E+02L,
588  -1.01515636749021914166146E+03L, -3.10777167157231109440444E+03L, 2.25381184209801510330112E+04L,
589  4.75584627752788110767815E+03L, -1.34659959864969306392456E+05L, -1.15132259675553483497211E+05L };
590  long double res, xnum, y, y1, ysq, z;
591  long double sqrtpi = 0.9189385332046727417803297L;
592  long double sum, value, xbig = 171.624L;
593  long double xden, xinf = 1.79E+308L;
594  long double xminin = 2.23E-308L;
595 
596  parity = false;
597  fact = 1.0L;
598  n = 0;
599  y = x;
600 
601  // When the argument is negative
602  if(y <= 0.0) {
603  y = - x;
604  y1 = (long double)((int)(y)); // take the integer part
605  res = y - y1;
606  if(res != 0.0) {
607  if( y1 != (long double)(((int)(y1*half))*2.0L)) {
608  parity = true;
609  }
610  fact = - M_PI/sinl(M_PI*res);
611  y = y + 1.0L;
612  }
613  else {
614  res = xinf;
615  value = res;
616  return value;
617  }
618  }
619 
620  // When the argument is positive but lower than precision
621  if(y < precision) {
622  if(xminin <= y) res = 1.0L/y;
623  else {
624  res = xinf;
625  value = res;
626  return value;
627  }
628  }
629  // When the argument is positive but lower than 12.0L
630  else if(y < 12.0L) {
631  y1 = y;
632  if(y < 1.0L) {
633  z = y;
634  y = y + 1.0L;
635  }
636  else { // note the value is not higher than 12.0L
637  n = (int)(y) - 1;
638  y = y - (long double)(n);
639  z = y - 1.0L;
640  }
641  xnum = 0.0;
642  xden = 1.0L;
643  for(i=0;i<8;i++) {
644  xnum = (xnum + p[i])*z;
645  xden = xden*z + q[i];
646  }
647  res = xnum/xden + 1.0L;
648  if(y1 < y) {
649  res = res/y1;
650  }
651  else if(y < y1) {
652  for(i=1;i<=n;i++) {
653  res = res*y;
654  y = y + 1.0L;
655  }
656  }
657  }
658  // When the argument is higher than 12.0L
659  else {
660  if(y <= xbig) {
661  ysq = y * y;
662  sum = c[6];
663  for(i=0;i<6;i++) {
664  sum = sum / ysq + c[i];
665  }
666  sum = sum / y - y + sqrtpi;
667  sum = sum + (y-half)*log(y);
668  res = exp(sum);
669  }
670  else {
671  res = xinf;
672  value = res;
673  return value;
674  }
675  }
676  // Adjusting
677  if(parity) res = - res;
678  if(fact != 1.0L) res = fact/res;
679 
680  return res;
681 }
TPZGaussRule::JacobiPolinomial
long double JacobiPolinomial(long double x, int order, long double alpha, long double beta, TPZVec< long double > &b, TPZVec< long double > &c, long double *dp2, long double *p1)
Evaluates the Jacobi polinomial for real x.
Definition: tpzgaussrule.cpp:527
fabs
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ fabs
Definition: tfadfunc.h:140
TPZGaussRule::SetType
void SetType(int &type, int order)
Sets a gaussian quadrature: 0 - Gauss Legendre, 1 - Gauss Lobatto, 2 - Gauss Jacobi with alpha = beta...
Definition: tpzgaussrule.cpp:131
IsZero
bool IsZero(long double a)
Returns if the value a is close Zero as the allowable tolerance.
Definition: pzreal.h:668
TPZGaussRule::W
long double W(int i) const
Returns weight for the ith point.
Definition: tpzgaussrule.cpp:107
TPZGaussRule::fWeight
TPZVec< long double > fWeight
Weight of the integration point.
Definition: tpzgaussrule.h:39
pzvec.h
Templated vector implementation.
pzerror.h
Defines PZError.
TPZGaussRule::ComputingGaussLobattoQuadrature
int ComputingGaussLobattoQuadrature(int order)
Computes the points and weights for Gauss Lobatto Quadrature over the parametric 1D element It is to...
Definition: tpzgaussrule.cpp:288
TPZGaussRule::fOrder
int fOrder
order of the integration rule
Definition: tpzgaussrule.h:47
TPZGaussRule::ComputingGaussJacobiQuadrature
int ComputingGaussJacobiQuadrature(int order, long double alpha, long double beta)
Computes the points and weights for Gauss Lobatto Quadrature over the parametric 1D element It is to...
Definition: tpzgaussrule.cpp:379
TPZGaussRule::ComputingGaussChebyshevQuadrature
int ComputingGaussChebyshevQuadrature(int order)
Computes the points and weights for Gauss Chebyshev Quadrature over the parametric 1D element [-1...
Definition: tpzgaussrule.cpp:242
PZINTEGRAL_MAXITERATIONS_ALLOWED
#define PZINTEGRAL_MAXITERATIONS_ALLOWED
Definition: tpzgaussrule.h:11
TPZGaussRule::operator=
TPZGaussRule & operator=(const TPZGaussRule &copy)
Definition: tpzgaussrule.cpp:72
TPZVec::Resize
virtual void Resize(const int64_t newsize, const T &object)
Resizes the vector object reallocating the necessary storage, copying the existing objects to the new...
Definition: pzvec.h:373
pzgeom::tol
static const double tol
Definition: pzgeoprism.cpp:23
test.f
f
Definition: test.py:287
DebugStop
#define DebugStop()
Returns a message to user put a breakpoint in.
Definition: pzerror.h:20
TPZGaussRule::fNumInt
int fNumInt
Number of integration points for this object.
Definition: tpzgaussrule.h:35
test.res
string res
Definition: test.py:151
TPZGaussRule::fBeta
long double fBeta
BETA is the exponent of (1+X) in the quadrature rule: weight = (1-X)^ALPHA * (1+X)^BETA.
Definition: tpzgaussrule.h:44
TPZGaussRule::TPZGaussRule
TPZGaussRule(int order, int type=0, long double alpha=0.0L, long double beta=0.0L)
Constructor of quadrature rule.
Definition: tpzgaussrule.cpp:22
gamma
long double gamma(unsigned int n)
Evaluate the factorial of a integer.
Definition: tpzgaussrule.cpp:567
pow
TPZFlopCounter pow(const TPZFlopCounter &orig, const TPZFlopCounter &xp)
Returns the power and increments the counter of the power.
Definition: pzreal.h:487
TPZGaussRule::Print
void Print(std::ostream &out=std::cout)
Prints the number of integration points, all points and weights (as one dimension) ...
Definition: tpzgaussrule.cpp:116
tpzintrulelist.h
Contains the TPZIntRuleList class which creates instances of all integration rules for rapid selectio...
log
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ log
Definition: tfadfunc.h:130
TPZGaussRule::~TPZGaussRule
~TPZGaussRule()
Default destructor.
Definition: tpzgaussrule.cpp:88
exp
expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ expr_ exp
Definition: tfadfunc.h:125
TPZGaussRule::fAlpha
long double fAlpha
ALPHA is the exponent of (1-X) in the quadrature rule: weight = (1-X)^ALPHA * (1+X)^BETA.
Definition: tpzgaussrule.h:42
TPZGaussRule::fLocation
TPZVec< long double > fLocation
Location of the integration point.
Definition: tpzgaussrule.h:37
machinePrecision
long double machinePrecision()
Definition: tpzgaussrule.cpp:552
TPZGaussRule::fType
int fType
fType = 1 (Gauss Lobatto quadrature), fType = 2 (Gauss Jacobi quadrature) for any alpha and beta the ...
Definition: tpzgaussrule.h:26
tpzgaussrule.h
Contains the TPZGaussRule class which implements the Gaussian quadrature.
m
clarg::argString m("-m", "input matrix file name (text format)", "matrix.txt")
TPZGaussRule::Loc
long double Loc(int i) const
Returns location of the ith point.
Definition: tpzgaussrule.cpp:96
TPZGaussRule
Implements the Gaussian quadrature. Numerical Integration Abstract class.
Definition: tpzgaussrule.h:19
TPZGaussRule::ComputingGaussLegendreQuadrature
int ComputingGaussLegendreQuadrature(int order)
Computes the points and weights for Gauss Legendre Quadrature over the parametric 1D element [-1...
Definition: tpzgaussrule.cpp:170
PZError
#define PZError
Defines the output device to error messages and the DebugStop() function.
Definition: pzerror.h:15
Generated on Sun Aug 16 2020 11:17:01 for NeoPZ by   doxygen 1.8.13