52 #include "Teuchos_oblackholestream.hpp" 53 #include "Teuchos_RCP.hpp" 54 #include "Teuchos_ScalarTraits.hpp" 55 #include "Teuchos_GlobalMPISession.hpp" 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);
int main(int argc, char *argv[])
Performs a code-code comparison to FIAT for Nedelec bases on tets (values and curls) ...
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...