Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Macros | Functions
Polylib_test.cpp File Reference
#include <stdio.h>
#include <math.h>
#include "Polylib.h"
Include dependency graph for Polylib_test.cpp:

Go to the source code of this file.

Macros

#define NPLOWER   5
 
#define NPUPPER   15
 
#define EPS   1e-12
 
#define GAUSS_INT   1
 
#define GAUSS_RADAUM_INT   1
 
#define GAUSS_RADAUP_INT   1
 
#define GAUSS_LOBATTO_INT   1
 
#define GAUSS_DIFF   1
 
#define GAUSS_RADAUM_DIFF   1
 
#define GAUSS_RADAUP_DIFF   1
 
#define GAUSS_LOBATTO_DIFF   1
 
#define GAUSS_INTERP   1
 
#define GAUSS_RADAUM_INTERP   1
 
#define GAUSS_RADAUP_INTERP   1
 
#define GAUSS_LOBATTO_INTERP   1
 

Functions

double ddot (int, double *, int, double *, int)
 
double * dvector (int, int)
 
 main ()
 

Macro Definition Documentation

#define EPS   1e-12

Definition at line 100 of file Polylib_test.cpp.

Referenced by main().

#define GAUSS_DIFF   1

Definition at line 106 of file Polylib_test.cpp.

#define GAUSS_INT   1

Definition at line 102 of file Polylib_test.cpp.

#define GAUSS_INTERP   1

Definition at line 110 of file Polylib_test.cpp.

#define GAUSS_LOBATTO_DIFF   1

Definition at line 109 of file Polylib_test.cpp.

#define GAUSS_LOBATTO_INT   1

Definition at line 105 of file Polylib_test.cpp.

#define GAUSS_LOBATTO_INTERP   1

Definition at line 113 of file Polylib_test.cpp.

#define GAUSS_RADAUM_DIFF   1

Definition at line 107 of file Polylib_test.cpp.

#define GAUSS_RADAUM_INT   1

Definition at line 103 of file Polylib_test.cpp.

#define GAUSS_RADAUM_INTERP   1

Definition at line 111 of file Polylib_test.cpp.

#define GAUSS_RADAUP_DIFF   1

Definition at line 108 of file Polylib_test.cpp.

#define GAUSS_RADAUP_INT   1

Definition at line 104 of file Polylib_test.cpp.

#define GAUSS_RADAUP_INTERP   1

Definition at line 112 of file Polylib_test.cpp.

#define NPLOWER   5

Definition at line 98 of file Polylib_test.cpp.

Referenced by main().

#define NPUPPER   15

Definition at line 99 of file Polylib_test.cpp.

Referenced by main().

Function Documentation

double ddot ( int  n,
double *  x,
int  incx,
double *  y,
int  incy 
)

Definition at line 486 of file Polylib_test.cpp.

Referenced by main().

487 {
488  register double sum = 0.;
489 
490  while (n--) {
491  sum += (*x) * (*y);
492  x += incx;
493  y += incy;
494  }
495  return sum;
496 }
double * dvector ( int  nl,
int  nh 
)

Definition at line 499 of file Polylib_test.cpp.

Referenced by main().

500 {
501  double *v;
502 
503  v = (double *)malloc((unsigned) (nh-nl+1)*sizeof(double));
504  return v-nl;
505 }
main ( )

Definition at line 119 of file Polylib_test.cpp.

References ddot(), Polylib::Dgj(), Polylib::Dglj(), Polylib::Dgrjm(), Polylib::Dgrjp(), dvector(), EPS, Polylib::Imgj(), Polylib::Imglj(), Polylib::Imgrjm(), Polylib::Imgrjp(), Polylib::jacobfd(), NPLOWER, NPUPPER, CellMLToNektar.cellml_metadata::p, Polylib::zwgj(), Polylib::zwglj(), Polylib::zwgrjm(), and Polylib::zwgrjp().

119  {
120  int np,n,i;
121  double *z,*w,*p,sum=0,alpha,beta,*d,*dt;
122 
123  z = dvector(0,NPUPPER-1);
124  w = dvector(0,NPUPPER-1);
125  p = dvector(0,NPUPPER-1);
126 
127  d = dvector(0,NPUPPER*NPUPPER-1);
128  dt = dvector(0,NPUPPER*NPUPPER-1);
129 
130 #if GAUSS_INT
131  /* Gauss Integration */
132  printf("Begin checking Gauss integration\n");
133  alpha = -0.5;
134  while(alpha <= 5.0){
135  beta = -0.5;
136  while(beta <= 5.0){
137 
138  for(np = NPLOWER; np <= NPUPPER; ++np){
139  zwgj(z,w,np,alpha,beta);
140  for(n = 2; n < 2*np-1; ++n){
141  jacobfd(np,z,p,NULL,n,alpha,beta);
142  sum = ddot(np,w,1,p,1);
143  if(fabs(sum)>EPS)
144  printf("alpha = %lf, beta = %lf, np = %d, n = %d integal was %lg\n"
145  ,alpha,beta,np,n,sum);
146  }
147  }
148 
149  beta += 0.5;
150  }
151  printf("finished checking all beta values for alpha = %lf\n",alpha);
152  alpha += 0.5;
153  }
154  printf("Finished checking Gauss Integration\n");
155 #endif
156 
157 #if GAUSS_RADAUM_INT
158  /* Gauss Radau Integration */
159  printf("Begin checking Gauss Radau Integration\n");
160  alpha = -0.5;
161  while(alpha <= 5.0){
162  beta = -0.5;
163  while(beta <= 5.0){
164  for(np = NPLOWER; np <= NPUPPER; ++np){
165  zwgrjm(z,w,np,alpha,beta);
166  for(n = 2; n < 2*np-2; ++n){
167  jacobfd(np,z,p,NULL,n,alpha,beta);
168  sum = ddot(np,w,1,p,1);
169  if(fabs(sum)>EPS)
170  printf("alpha = %lf, beta = %lf, np = %d, n = %d integal was %lg\n"
171  ,alpha,beta,np,n,sum);
172  }
173  }
174 
175  beta += 0.5;
176  }
177  printf("finished checking all beta values for alpha = %lf\n",alpha);
178  alpha += 0.5;
179  }
180  printf("Finished checking Gauss Radau (z=-1) Integration\n");
181 #endif
182 
183 
184 #if GAUSS_RADAUP_INT
185  /* Gauss Radau Integration */
186  printf("Begin checking Gauss Radau Integration\n");
187  alpha = -0.5;
188  while(alpha <= 5.0){
189  beta = -0.5;
190  while(beta <= 5.0){
191  for(np = NPLOWER; np <= NPUPPER; ++np){
192  zwgrjp(z,w,np,alpha,beta);
193  for(n = 2; n < 2*np-2; ++n){
194  jacobfd(np,z,p,NULL,n,alpha,beta);
195  sum = ddot(np,w,1,p,1);
196  if(fabs(sum)>EPS)
197  printf("alpha = %lf, beta = %lf, np = %d, n = %d integal was %lg\n"
198  ,alpha,beta,np,n,sum);
199  }
200  }
201 
202  beta += 0.5;
203  }
204  printf("finished checking all beta values for alpha = %lf\n",alpha);
205  alpha += 0.5;
206  }
207  printf("Finished checking Gauss Radau (z=1) Integration\n");
208 #endif
209 
210 #if GAUSS_LOBATTO_INT
211  /* Gauss Lobatto Integration */
212  printf("Begin checking Gauss Lobatto integration\n");
213  alpha = -0.5;
214  while(alpha <= 5.0){
215  beta = -0.5;
216  while(beta <= 5.0){
217 
218  for(np = NPLOWER; np <= NPUPPER; ++np){
219  zwglj(z,w,np,alpha,beta);
220  for(n = 2; n < 2*np-3; ++n){
221  jacobfd(np,z,p,NULL,n,alpha,beta);
222  sum = ddot(np,w,1,p,1);
223  if(fabs(sum)>EPS)
224  printf("alpha = %lf, beta = %lf, np = %d, n = %d integal was %lg\n"
225  ,alpha,beta,np,n,sum);
226  }
227  }
228 
229  beta += 0.5;
230  }
231  printf("finished checking all beta values for alpha = %lf\n",alpha);
232  alpha += 0.5;
233  }
234  printf("Finished checking Gauss Lobatto Integration\n");
235 #endif
236 
237 #if GAUSS_DIFF
238  printf("Begin checking differentiation through Gauss points\n");
239  alpha = -0.5;
240  while(alpha <= 5.0){
241  beta = -0.5;
242  while(beta <= 5.0){
243 
244  for(np = NPLOWER; np <= NPUPPER; ++np){
245  zwgj(z,w,np,alpha,beta);
246  for(n = 2; n < np-1; ++n){
247  Dgj(d,dt,z,np,alpha,beta);
248 
249  for(i = 0; i < np; ++i) p[i] = pow(z[i],n);
250  sum = 0;
251  for(i = 0; i < np; ++i)
252  sum += fabs(ddot(np,d+i*np,1,p,1) - n*pow(z[i],n-1));
253  sum /= np;
254  if(fabs(sum)>EPS)
255  printf("alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
256  ,alpha,beta,np,n,sum);
257  }
258  }
259  beta += 0.5;
260  }
261  printf("finished checking all beta values for alpha = %lf\n",alpha);
262  alpha += 0.5;
263  }
264  printf("Finished checking Gauss Jacobi differentiation\n");
265 #endif
266 
267 #if GAUSS_RADAUM_DIFF
268  printf("Begin checking differentiation through Gauss Radau points\n");
269  alpha = -0.5;
270  while(alpha <= 5.0){
271  beta = -0.5;
272  while(beta <= 5.0){
273 
274  for(np = NPLOWER; np <= NPUPPER; ++np){
275  zwgrjm(z,w,np,alpha,beta);
276  for(n = 2; n < np-1; ++n){
277  Dgrjm(d,dt,z,np,alpha,beta);
278  for(i = 0; i < np; ++i) p[i] = pow(z[i],n);
279  sum = 0;
280  for(i = 0; i < np; ++i)
281  sum += fabs(ddot(np,d+i*np,1,p,1) - n*pow(z[i],n-1));
282  sum /= np;
283  if(fabs(sum)>EPS)
284  printf("alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
285  ,alpha,beta,np,n,sum);
286  }
287  }
288  beta += 0.5;
289  }
290  printf("finished checking all beta values for alpha = %lf\n",alpha);
291  alpha += 0.5;
292  }
293  printf("Finished checking Gauss Radau (z=-1) differentiation\n");
294 #endif
295 
296 #if GAUSS_RADAUP_DIFF
297  printf("Begin checking differentiation through Gauss Radau (z=1) points\n");
298  alpha = -0.5;
299  while(alpha <= 5.0){
300  beta = -0.5;
301  while(beta <= 5.0){
302 
303  for(np = NPLOWER; np <= NPUPPER; ++np){
304  zwgrjp(z,w,np,alpha,beta);
305  for(n = 2; n < np-1; ++n){
306  Dgrjp(d,dt,z,np,alpha,beta);
307  for(i = 0; i < np; ++i) p[i] = pow(z[i],n);
308  sum = 0;
309  for(i = 0; i < np; ++i)
310  sum += fabs(ddot(np,d+i*np,1,p,1) - n*pow(z[i],n-1));
311  sum /= np;
312  if(fabs(sum)>EPS)
313  printf("alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
314  ,alpha,beta,np,n,sum);
315  }
316  }
317  beta += 0.5;
318  }
319  printf("finished checking all beta values for alpha = %lf\n",alpha);
320  alpha += 0.5;
321  }
322  printf("Finished checking Gauss Radau (z=1) differentiation\n");
323 #endif
324 
325 #if GAUSS_LOBATTO_DIFF
326  printf("Begin checking differentiation through Gauss Lobatto points\n");
327  alpha = -0.5;
328  while(alpha <= 5.0){
329  beta = -0.5;
330  while(beta <= 5.0){
331 
332  for(np = NPLOWER; np <= NPUPPER; ++np){
333  zwglj(z,w,np,alpha,beta);
334  for(n = 2; n < np-1; ++n){
335  Dglj(d,dt,z,np,alpha,beta);
336  for(i = 0; i < np; ++i) p[i] = pow(z[i],n);
337  sum = 0;
338  for(i = 0; i < np; ++i)
339  sum += fabs(ddot(np,d+i*np,1,p,1) - n*pow(z[i],n-1));
340  sum /= np;
341  if(fabs(sum)>EPS)
342  printf("alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
343  ,alpha,beta,np,n,sum);
344  }
345  }
346  beta += 0.5;
347  }
348  printf("finished checking all beta values for alpha = %lf\n",alpha);
349  alpha += 0.5;
350  }
351  printf("Finished checking Gauss Lobatto differentiation\n");
352 #endif
353 
354  /* check interpolation routines */
355 #if GAUSS_INTERP
356  printf("Begin checking interpolation through Gauss points\n");
357  alpha = -0.5;
358  while(alpha <= 5.0){
359  beta = -0.5;
360  while(beta <= 5.0){
361 
362  for(np = NPLOWER; np <= NPUPPER; ++np){
363  zwgj(z,w,np,alpha,beta);
364  for(n = 2; n < np-1; ++n){
365  for(i = 0; i < np; ++i) {
366  w[i] = 2.0*i/(double)(np-1)-1.0;
367  p[i] = pow(z[i],n);
368  }
369  Imgj(d,z,w,np,np,alpha,beta);
370  sum = 0;
371  for(i = 0; i < np; ++i)
372  sum += fabs(ddot(np,d+i*np,1,p,1) - pow(w[i],n));
373  sum /= np;
374  if(fabs(sum)>EPS)
375  printf("alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
376  ,alpha,beta,np,n,sum);
377  }
378  }
379  beta += 0.5;
380  }
381  printf("finished checking all beta values for alpha = %lf\n",alpha);
382  alpha += 0.5;
383  }
384  printf("Finished checking Gauss Jacobi interpolation\n");
385 #endif
386 
387 #if GAUSS_RADAUM_INTERP
388  printf("Begin checking Interpolation through Gauss Radau (z=-1) points\n");
389  alpha = -0.5;
390  while(alpha <= 5.0){
391  beta = -0.5;
392  while(beta <= 5.0){
393 
394  for(np = NPLOWER; np <= NPUPPER; ++np){
395  zwgrjm(z,w,np,alpha,beta);
396  for(n = 2; n < np-1; ++n){
397  for(i = 0; i < np; ++i) {
398  w[i] = 2.0*i/(double)(np-1)-1.0;
399  p[i] = pow(z[i],n);
400  }
401  Imgrjm(d,z,w,np,np,alpha,beta);
402  sum = 0;
403  for(i = 0; i < np; ++i)
404  sum += fabs(ddot(np,d+i*np,1,p,1) - pow(w[i],n));
405  sum /= np;
406  if(fabs(sum)>EPS)
407  printf("alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
408  ,alpha,beta,np,n,sum);
409  }
410  }
411  beta += 0.5;
412  }
413  printf("finished checking all beta values for alpha = %lf\n",alpha);
414  alpha += 0.5;
415  }
416  printf("Finished checking Gauss Radua Jacobi (z=-1) interpolation\n");
417 #endif
418 #if GAUSS_RADAUP_INTERP
419  printf("Begin checking Interpolation through Gauss Radau (z=1) points\n");
420  alpha = -0.5;
421  while(alpha <= 5.0){
422  beta = -0.5;
423  while(beta <= 5.0){
424 
425  for(np = NPLOWER; np <= NPUPPER; ++np){
426  zwgrjp(z,w,np,alpha,beta);
427  for(n = 2; n < np-1; ++n){
428  for(i = 0; i < np; ++i) {
429  w[i] = 2.0*i/(double)(np-1)-1.0;
430  p[i] = pow(z[i],n);
431  }
432  Imgrjp(d,z,w,np,np,alpha,beta);
433  sum = 0;
434  for(i = 0; i < np; ++i)
435  sum += fabs(ddot(np,d+i*np,1,p,1) - pow(w[i],n));
436  sum /= np;
437  if(fabs(sum)>EPS)
438  printf("alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
439  ,alpha,beta,np,n,sum);
440  }
441  }
442  beta += 0.5;
443  }
444  printf("finished checking all beta values for alpha = %lf\n",alpha);
445  alpha += 0.5;
446  }
447  printf("Finished checking Gauss Radau (z=1) interpolation\n");
448 #endif
449 
450 #if GAUSS_LOBATTO_INTERP
451  printf("Begin checking Interpolation through Gauss Lobatto points\n");
452  alpha = -0.5;
453  while(alpha <= 5.0){
454  beta = -0.5;
455  while(beta <= 5.0){
456 
457  for(np = NPLOWER; np <= NPUPPER; ++np){
458  zwglj(z,w,np,alpha,beta);
459  for(n = 2; n < np-1; ++n){
460  for(i = 0; i < np; ++i) {
461  w[i] = 2.0*i/(double)(np-1)-1.0;
462  p[i] = pow(z[i],n);
463  }
464  Imglj(d,z,w,np,np,alpha,beta);
465  sum = 0;
466  for(i = 0; i < np; ++i)
467  sum += fabs(ddot(np,d+i*np,1,p,1) - pow(w[i],n));
468  sum /= np;
469  if(fabs(sum)>EPS)
470  printf("alpha = %lf, beta = %lf, np = %d, n = %d difference %lg\n"
471  ,alpha,beta,np,n,sum);
472  }
473  }
474  beta += 0.5;
475  }
476  printf("finished checking all beta values for alpha = %lf\n",alpha);
477  alpha += 0.5;
478  }
479  printf("Finished checking Gauss Lobatto interploation\n");
480 #endif
481 
482 
483  free(z); free(w); free(p); free(d); free(dt);
484 }
NPUPPER
#define NPUPPER
Definition: Polylib_test.cpp:99
Polylib::Dgj
void Dgj(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix and its transpose associated.
Definition: Polylib.cpp:937
EPS
#define EPS
Definition: Polylib_test.cpp:100
ddot
double ddot(int, double *, int, double *, int)
Definition: Polylib_test.cpp:486
Polylib::zwgrjp
void zwgrjp(double *z, double *w, const int np, const double alpha, const double beta)
Gauss-Radau-Jacobi zeros and weights with end point at z=1.
Definition: Polylib.cpp:239
Polylib::Imgrjp
void Imgrjp(double *im, const double *zgrj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
Interpolation Operator from Gauss-Radau-Jacobi points.
Definition: Polylib.cpp:1775
Polylib::zwglj
void zwglj(double *z, double *w, const int np, const double alpha, const double beta)
Gauss-Lobatto-Jacobi zeros and weights with end point at z=-1,1.
Definition: Polylib.cpp:311
dvector
double * dvector(int, int)
Definition: Polylib_test.cpp:499
CellMLToNektar.cellml_metadata.p
tuple p
Definition: cellml_metadata.py:350
Polylib::Imglj
void Imglj(double *im, const double *zglj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
Interpolation Operator from Gauss-Lobatto-Jacobi points.
Definition: Polylib.cpp:1837
Polylib::Dgrjm
void Dgrjm(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix and its transpose associated.
Definition: Polylib.cpp:1023
Polylib::Imgrjm
void Imgrjm(double *im, const double *zgrj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
Interpolation Operator from Gauss-Radau-Jacobi points.
Definition: Polylib.cpp:1715
NPLOWER
#define NPLOWER
Definition: Polylib_test.cpp:98
Polylib::Dgrjp
void Dgrjp(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix associated with the.
Definition: Polylib.cpp:1121
Polylib::zwgrjm
void zwgrjm(double *z, double *w, const int np, const double alpha, const double beta)
Gauss-Radau-Jacobi zeros and weights with end point at z=-1.
Definition: Polylib.cpp:163
Polylib::Dglj
void Dglj(double *D, const double *z, const int np, const double alpha, const double beta)
Compute the Derivative Matrix associated with the.
Definition: Polylib.cpp:1221
Polylib::Imgj
void Imgj(double *im, const double *zgj, const double *zm, const int nz, const int mz, const double alpha, const double beta)
Interpolation Operator from Gauss-Jacobi points to an.
Definition: Polylib.cpp:1657
Polylib::jacobfd
void jacobfd(const int np, const double *z, double *poly_in, double *polyd, const int n, const double alpha, const double beta)
Routine to calculate Jacobi polynomials, , and their first derivative, .
Definition: Polylib.cpp:1951
Polylib::zwgj
void zwgj(double *z, double *w, const int np, const double alpha, const double beta)
Gauss-Jacobi zeros and weights.
Definition: Polylib.cpp:107
Generated on Wed Mar 8 2017 12:29:16 for Nektar++ by   doxygen 1.8.8