52 #include "Teuchos_oblackholestream.hpp"
53 #include "Teuchos_RCP.hpp"
54 #include "Teuchos_ScalarTraits.hpp"
55 #include "Teuchos_GlobalMPISession.hpp"
59 using namespace Intrepid;
63 #define TEST_EPS 1000*INTREPID_TOL
66 #define GAUSS_RADAUM_INT 1
67 #define GAUSS_RADAUP_INT 1
68 #define GAUSS_LOBATTO_INT 1
70 #define GAUSS_RADAUM_DIFF 1
71 #define GAUSS_RADAUP_DIFF 1
72 #define GAUSS_LOBATTO_DIFF 1
73 #define GAUSS_INTERP 1
74 #define GAUSS_RADAUM_INTERP 1
75 #define GAUSS_RADAUP_INTERP 1
76 #define GAUSS_LOBATTO_INTERP 1
79 #define INTREPID_TEST_COMMAND( S ) \
84 catch (std::logic_error err) { \
85 *outStream << "Expected Error ----------------------------------------------------------------\n"; \
86 *outStream << err.what() << '\n'; \
87 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \
92 template<
class Scalar>
93 Scalar ddot(
int n, Scalar *x,
int incx, Scalar *y,
int incy)
105 template<
class Scalar>
106 Scalar * dvector(
int nl,
int nh)
110 v = (Scalar *)malloc((
unsigned) (nh-nl+1)*
sizeof(Scalar));
193 int main(
int argc,
char *argv[]) {
195 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
199 int iprint = argc - 1;
200 Teuchos::RCP<std::ostream> outStream;
201 Teuchos::oblackholestream bhs;
203 outStream = Teuchos::rcp(&std::cout,
false);
205 outStream = Teuchos::rcp(&bhs,
false);
208 Teuchos::oblackholestream oldFormatState;
209 oldFormatState.copyfmt(std::cout);
212 <<
"===============================================================================\n" \
214 <<
"| Unit Test (IntrepidPolylib) |\n" \
216 <<
"| 1) the original Polylib tests |\n" \
218 <<
"| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \
219 <<
"| Denis Ridzal (dridzal@sandia.gov). |\n" \
221 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
222 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
224 <<
"===============================================================================\n";
228 int beginThrowNumber = Teuchos::TestForException_getThrowNumber();
229 int endThrowNumber = beginThrowNumber + 1;
236 <<
"===============================================================================\n"\
237 <<
"| TEST 1: exceptions |\n"\
238 <<
"===============================================================================\n";
241 INTREPID_TEST_COMMAND( iplib.
gammaF((
double)0.33) );
243 catch (std::logic_error err) {
244 *outStream <<
"UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
245 *outStream << err.what() <<
'\n';
246 *outStream <<
"-------------------------------------------------------------------------------" <<
"\n\n";
250 if (Teuchos::TestForException_getThrowNumber() != endThrowNumber)
255 <<
"===============================================================================\n"\
256 <<
"| TEST 2: correctness of math operations |\n"\
257 <<
"===============================================================================\n";
259 outStream->precision(20);
264 double *z,*w,*p,sum=0,alpha,beta,*d;
266 z = dvector<double>(0,NPUPPER-1);
267 w = dvector<double>(0,NPUPPER-1);
268 p = dvector<double>(0,NPUPPER-1);
270 d = dvector<double>(0,NPUPPER*NPUPPER-1);
274 *outStream <<
"Begin checking Gauss integration\n";
280 for(np = NPLOWER; np <= NPUPPER; ++np){
281 ipl::zwgj(z,w,np,alpha,beta);
282 for(n = 2; n < 2*np-1; ++n){
283 ipl::jacobfd(np,z,p,(
double*)0,n,alpha,beta);
284 sum = ddot(np,w,1,p,1);
285 if(std::abs(sum)>TEST_EPS) {
287 *outStream <<
"ERROR: alpha = " << alpha <<
", beta = " << beta <<
288 ", np = " << np <<
", n = " << n <<
" integral was " << sum <<
"\n";
295 *outStream <<
"finished checking all beta values for alpha = " << alpha <<
"\n";
298 *outStream <<
"Finished checking Gauss Integration\n";
303 *outStream <<
"Begin checking Gauss-Radau (z = -1) integration\n";
309 for(np = NPLOWER; np <= NPUPPER; ++np){
310 ipl::zwgrjm(z,w,np,alpha,beta);
311 for(n = 2; n < 2*np-2; ++n){
312 ipl::jacobfd(np,z,p,(
double*)0,n,alpha,beta);
313 sum = ddot(np,w,1,p,1);
314 if(std::abs(sum)>TEST_EPS) {
316 *outStream <<
"ERROR: alpha = " << alpha <<
", beta = " << beta <<
317 ", np = " << np <<
", n = " << n <<
" integral was " << sum <<
"\n";
324 *outStream <<
"finished checking all beta values for alpha = " << alpha <<
"\n";
327 *outStream <<
"Finished checking Gauss-Radau (z = -1) Integration\n";
332 *outStream <<
"Begin checking Gauss-Radau (z = +1) integration\n";
338 for(np = NPLOWER; np <= NPUPPER; ++np){
339 ipl::zwgrjp(z,w,np,alpha,beta);
340 for(n = 2; n < 2*np-2; ++n){
341 ipl::jacobfd(np,z,p,(
double*)0,n,alpha,beta);
342 sum = ddot(np,w,1,p,1);
343 if(std::abs(sum)>TEST_EPS) {
345 *outStream <<
"ERROR: alpha = " << alpha <<
", beta = " << beta <<
346 ", np = " << np <<
", n = " << n <<
" integral was " << sum <<
"\n";
353 *outStream <<
"finished checking all beta values for alpha = " << alpha <<
"\n";
356 *outStream <<
"Finished checking Gauss-Radau (z = +1) Integration\n";
359 #if GAUSS_LOBATTO_INT
361 *outStream <<
"Begin checking Gauss-Lobatto integration\n";
367 for(np = NPLOWER; np <= NPUPPER; ++np){
368 ipl::zwglj(z,w,np,alpha,beta);
369 for(n = 2; n < 2*np-3; ++n){
370 ipl::jacobfd(np,z,p,(
double*)0,n,alpha,beta);
371 sum = ddot(np,w,1,p,1);
372 if(std::abs(sum)>TEST_EPS) {
374 *outStream <<
"ERROR: alpha = " << alpha <<
", beta = " << beta <<
375 ", np = " << np <<
", n = " << n <<
" integral was " << sum <<
"\n";
382 *outStream <<
"finished checking all beta values for alpha = " << alpha <<
"\n";
385 *outStream <<
"Finished checking Gauss-Lobatto Integration\n";
389 *outStream <<
"Begin checking differentiation through Gauss points\n";
395 for(np = NPLOWER; np <= NPUPPER; ++np){
396 ipl::zwgj(z,w,np,alpha,beta);
397 for(n = 2; n < np-1; ++n){
398 ipl::Dgj(d,z,np,alpha,beta);
400 for(i = 0; i < np; ++i) p[i] = std::pow(z[i],n);
402 for(i = 0; i < np; ++i)
403 sum += std::abs(ddot(np,d+i*np,1,p,1) - n*std::pow(z[i],n-1));
405 if(std::abs(sum)>TEST_EPS) {
407 *outStream <<
"ERROR: alpha = " << alpha <<
", beta = " << beta <<
408 ", np = " << np <<
", n = " << n <<
" difference " << sum <<
"\n";
414 *outStream <<
"finished checking all beta values for alpha = " << alpha <<
"\n";
417 *outStream <<
"Finished checking Gauss Jacobi differentiation\n";
420 #if GAUSS_RADAUM_DIFF
421 *outStream <<
"Begin checking differentiation through Gauss-Radau (z=-1) points\n";
427 for(np = NPLOWER; np <= NPUPPER; ++np){
428 ipl::zwgrjm(z,w,np,alpha,beta);
429 for(n = 2; n < np-1; ++n){
430 ipl::Dgrjm(d,z,np,alpha,beta);
432 for(i = 0; i < np; ++i) p[i] = std::pow(z[i],n);
434 for(i = 0; i < np; ++i)
435 sum += std::abs(ddot(np,d+i*np,1,p,1) - n*std::pow(z[i],n-1));
437 if(std::abs(sum)>TEST_EPS) {
439 *outStream <<
"ERROR: alpha = " << alpha <<
", beta = " << beta <<
440 ", np = " << np <<
", n = " << n <<
" difference " << sum <<
"\n";
446 *outStream <<
"finished checking all beta values for alpha = " << alpha <<
"\n";
449 *outStream <<
"Finished checking Gauss-Radau (z=-1) differentiation\n";
452 #if GAUSS_RADAUP_DIFF
453 *outStream <<
"Begin checking differentiation through Gauss-Radau (z=+1) points\n";
459 for(np = NPLOWER; np <= NPUPPER; ++np){
460 ipl::zwgrjp(z,w,np,alpha,beta);
461 for(n = 2; n < np-1; ++n){
462 ipl::Dgrjp(d,z,np,alpha,beta);
464 for(i = 0; i < np; ++i) p[i] = std::pow(z[i],n);
466 for(i = 0; i < np; ++i)
467 sum += std::abs(ddot(np,d+i*np,1,p,1) - n*std::pow(z[i],n-1));
469 if(std::abs(sum)>TEST_EPS) {
471 *outStream <<
"ERROR: alpha = " << alpha <<
", beta = " << beta <<
472 ", np = " << np <<
", n = " << n <<
" difference " << sum <<
"\n";
478 *outStream <<
"finished checking all beta values for alpha = " << alpha <<
"\n";
481 *outStream <<
"Finished checking Gauss-Radau (z=+1) differentiation\n";
484 #if GAUSS_RADAUP_DIFF
485 *outStream <<
"Begin checking differentiation through Gauss-Lobatto (z=+1) points\n";
491 for(np = NPLOWER; np <= NPUPPER; ++np){
492 ipl::zwglj(z,w,np,alpha,beta);
493 for(n = 2; n < np-1; ++n){
494 ipl::Dglj(d,z,np,alpha,beta);
496 for(i = 0; i < np; ++i) p[i] = std::pow(z[i],n);
498 for(i = 0; i < np; ++i)
499 sum += std::abs(ddot(np,d+i*np,1,p,1) - n*std::pow(z[i],n-1));
501 if(std::abs(sum)>TEST_EPS) {
503 *outStream <<
"ERROR: alpha = " << alpha <<
", beta = " << beta <<
504 ", np = " << np <<
", n = " << n <<
" difference " << sum <<
"\n";
510 *outStream <<
"finished checking all beta values for alpha = " << alpha <<
"\n";
513 *outStream <<
"Finished checking Gauss-Lobatto differentiation\n";
517 *outStream <<
"Begin checking interpolation through Gauss points\n";
523 for(np = NPLOWER; np <= NPUPPER; ++np){
524 ipl::zwgj(z,w,np,alpha,beta);
525 for(n = 2; n < np-1; ++n){
526 for(i = 0; i < np; ++i) {
527 w[i] = 2.0*i/(double)(np-1)-1.0;
528 p[i] = std::pow(z[i],n);
530 ipl::Imgj(d,z,w,np,np,alpha,beta);
532 for(i = 0; i < np; ++i)
533 sum += std::abs(ddot(np,d+i*np,1,p,1) - std::pow(w[i],n));
535 if(std::abs(sum)>TEST_EPS) {
537 *outStream <<
"ERROR: alpha = " << alpha <<
", beta = " << beta <<
538 ", np = " << np <<
", n = " << n <<
" difference " << sum <<
"\n";
544 *outStream <<
"finished checking all beta values for alpha = " << alpha <<
"\n";
547 *outStream <<
"Finished checking Gauss Jacobi interpolation\n";
550 #if GAUSS_RADAUM_INTERP
551 *outStream <<
"Begin checking interpolation through Gauss-Radau (z=-1) points\n";
557 for(np = NPLOWER; np <= NPUPPER; ++np){
558 ipl::zwgrjm(z,w,np,alpha,beta);
559 for(n = 2; n < np-1; ++n){
560 for(i = 0; i < np; ++i) {
561 w[i] = 2.0*i/(double)(np-1)-1.0;
562 p[i] = std::pow(z[i],n);
564 ipl::Imgrjm(d,z,w,np,np,alpha,beta);
566 for(i = 0; i < np; ++i)
567 sum += std::abs(ddot(np,d+i*np,1,p,1) - std::pow(w[i],n));
569 if(std::abs(sum)>TEST_EPS) {
571 *outStream <<
"ERROR: alpha = " << alpha <<
", beta = " << beta <<
572 ", np = " << np <<
", n = " << n <<
" difference " << sum <<
"\n";
578 *outStream <<
"finished checking all beta values for alpha = " << alpha <<
"\n";
581 *outStream <<
"Finished checking Gauss-Radau (z=-1) interpolation\n";
584 #if GAUSS_RADAUP_INTERP
585 *outStream <<
"Begin checking interpolation through Gauss-Radau (z=+1) points\n";
591 for(np = NPLOWER; np <= NPUPPER; ++np){
592 ipl::zwgrjp(z,w,np,alpha,beta);
593 for(n = 2; n < np-1; ++n){
594 for(i = 0; i < np; ++i) {
595 w[i] = 2.0*i/(double)(np-1)-1.0;
596 p[i] = std::pow(z[i],n);
598 ipl::Imgrjp(d,z,w,np,np,alpha,beta);
600 for(i = 0; i < np; ++i)
601 sum += std::abs(ddot(np,d+i*np,1,p,1) - std::pow(w[i],n));
603 if(std::abs(sum)>TEST_EPS) {
605 *outStream <<
"ERROR: alpha = " << alpha <<
", beta = " << beta <<
606 ", np = " << np <<
", n = " << n <<
" difference " << sum <<
"\n";
612 *outStream <<
"finished checking all beta values for alpha = " << alpha <<
"\n";
615 *outStream <<
"Finished checking Gauss-Radau (z=+1) interpolation\n";
618 #if GAUSS_LOBATTO_INTERP
619 *outStream <<
"Begin checking interpolation through Gauss-Lobatto points\n";
625 for(np = NPLOWER; np <= NPUPPER; ++np){
626 ipl::zwglj(z,w,np,alpha,beta);
627 for(n = 2; n < np-1; ++n){
628 for(i = 0; i < np; ++i) {
629 w[i] = 2.0*i/(double)(np-1)-1.0;
630 p[i] = std::pow(z[i],n);
632 ipl::Imglj(d,z,w,np,np,alpha,beta);
634 for(i = 0; i < np; ++i)
635 sum += std::abs(ddot(np,d+i*np,1,p,1) - std::pow(w[i],n));
637 if(std::abs(sum)>TEST_EPS) {
639 *outStream <<
"ERROR: alpha = " << alpha <<
", beta = " << beta <<
640 ", np = " << np <<
", n = " << n <<
" difference " << sum <<
"\n";
646 *outStream <<
"finished checking all beta values for alpha = " << alpha <<
"\n";
649 *outStream <<
"Finished checking Gauss-Lobatto interpolation\n";
652 free(z); free(w); free(p); free(d);
659 catch (std::logic_error err) {
660 *outStream <<
"UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
661 *outStream << err.what() <<
'\n';
662 *outStream <<
"-------------------------------------------------------------------------------" <<
"\n\n";
668 std::cout <<
"End Result: TEST FAILED\n";
670 std::cout <<
"End Result: TEST PASSED\n";
673 std::cout.copyfmt(oldFormatState);
Providing orthogonal polynomial calculus and interpolation, created by Spencer Sherwin, Aeronautics, Imperial College London, modified and redistributed by D. Ridzal.
static Scalar gammaF(const Scalar x)
Calculate the Gamma function , , for integer values x and halves.
Header file for a set of functions providing orthogonal polynomial polynomial calculus and interpolat...