Intrepid
http://trilinos.sandia.gov/packages/docs/r10.10/packages/intrepid/test/Shared/IntrepidBurkardtRules/test_01.cpp
00001 // @HEADER
00002 // ************************************************************************
00003 //
00004 //                           Intrepid Package
00005 //                 Copyright (2007) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Pavel Bochev  (pbboche@sandia.gov)
00038 //                    Denis Ridzal  (dridzal@sandia.gov), or
00039 //                    Kara Peterson (kjpeter@sandia.gov)
00040 //
00041 // ************************************************************************
00042 // @HEADER
00043 
00051 #include "Intrepid_BurkardtRules.hpp"
00052 #include "Teuchos_oblackholestream.hpp"
00053 #include "Teuchos_RCP.hpp"
00054 #include "Teuchos_ScalarTraits.hpp"
00055 #include "Teuchos_GlobalMPISession.hpp"
00056 #include "Teuchos_Array.hpp"
00057 
00058 //# include <cstdlib>
00059 //# include <iostream>
00060 //# include <cmath>
00061 //# include <iomanip>
00062 
00063 using namespace std;
00064 using namespace Intrepid;
00065 
00066 #define INTREPID_TEST_COMMAND( S )                                                                                  \
00067 {                                                                                                                   \
00068   try {                                                                                                             \
00069     S ;                                                                                                             \
00070   }                                                                                                                 \
00071   catch (std::logic_error err) {                                                                                    \
00072       *outStream << "Expected Error ----------------------------------------------------------------\n";            \
00073       *outStream << err.what() << '\n';                                                                             \
00074       *outStream << "-------------------------------------------------------------------------------" << "\n\n";    \
00075   };                                                                                                                \
00076 }
00077 
00078 template<class Scalar>
00079 Scalar evalQuad(int order, int power, Scalar x[], Scalar w[]) {
00080  
00081   int mid  = order/2;
00082   Scalar Q = 0.0;
00083   if (order%2) 
00084     Q = w[mid]*powl(x[mid],power);
00085 
00086   for (int i=0; i<mid; i++) {
00087     Q += w[i]*powl(x[i],power)+w[order-i-1]*powl(x[order-i-1],power);
00088   }
00089 
00090   return Q;
00091   /* 
00092   Scalar Q = 0.0;
00093   for (int i=0; i<order; i++) {
00094     Q += w[i]*powl(x[i],power);
00095   }
00096   return Q;
00097   */
00098 }
00099 
00100 template<class Scalar>
00101 Scalar factorial2 (int n) {
00102   Scalar value = 1.0;
00103   if (n<1)
00104     return value;
00105 
00106   int n_copy = n;
00107   while (1<n_copy) {
00108     value  *= (Scalar)n_copy;
00109     n_copy -= 2;
00110   }
00111   return value;
00112 }
00113 
00114 template<class Scalar>
00115 Scalar chebyshev1(int power) {
00116   Scalar bot, exact, top;
00117   if (!(power%2)) {
00118     top = 1; bot = 1;
00119     for (int i=2;i<=power;i+=2) {
00120       top *= (Scalar)(i-1);
00121       bot *= (Scalar)i;
00122     }
00123     exact = M_PI*top/bot;
00124   }
00125   else {
00126     exact = 0.0;
00127   }
00128   return exact;
00129 }
00130 
00131 template<class Scalar>
00132 Scalar chebyshev2(int power) {
00133   Scalar bot, exact, top;
00134   if (!(power%2)) {
00135     top = 1; bot = 1;
00136     for (int i=2;i<=power;i+=2) {
00137       top *= (Scalar)(i-1);
00138       bot *= (Scalar)i;
00139     }
00140     bot *= (Scalar)(power+2);
00141     exact = M_PI*top/bot;
00142   }
00143   else {
00144     exact = 0.0;
00145   }
00146   return exact;
00147 }
00148 
00149 int main(int argc, char *argv[]) {
00150   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00151 
00152   // This little trick lets us print to std::cout only if
00153   // a (dummy) command-line argument is provided.
00154   int iprint     = argc - 1;
00155   Teuchos::RCP<std::ostream> outStream;
00156   Teuchos::oblackholestream bhs; // outputs nothing
00157   if (iprint > 0)
00158     outStream = Teuchos::rcp(&std::cout, false);
00159   else
00160     outStream = Teuchos::rcp(&bhs, false);
00161 
00162   // Save the format state of the original std::cout.
00163   Teuchos::oblackholestream oldFormatState;
00164   oldFormatState.copyfmt(std::cout);
00165 
00166   *outStream \
00167   << "===============================================================================\n" \
00168   << "|                                                                             |\n" \
00169   << "|                       Unit Test (IntrepidBurkardtRules)                     |\n" \
00170   << "|                                                                             |\n" \
00171   << "|     1) the Burkardt rule tests                                              |\n" \
00172   << "|                                                                             |\n" \
00173   << "|  Questions? Contact  Drew Kouri (dpkouri@sandia.gov) or                     |\n" \
00174   << "|                      Denis Ridzal (dridzal@sandia.gov).                     |\n" \
00175   << "|                                                                             |\n" \
00176   << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00177   << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00178   << "|                                                                             |\n" \
00179   << "===============================================================================\n";
00180 
00181 
00182   int errorFlag = 0;
00183 
00184   int maxOrder  = 30;
00185   long double reltol = 1e-8;
00186   long double analyticInt = 0, testInt = 0;
00187   // compute and compare integrals
00188   try {
00189     
00190     *outStream << "Gauss-Legendre Cubature \n";
00191     *outStream << "Integrates functions on [-1,1] weighted by w(x) = 1\n";
00192     for (int i = 1; i<=maxOrder; i++) {
00193       Teuchos::Array<long double> nodes(i), weights(i);
00194       IntrepidBurkardtRules::legendre_compute(i,nodes.getRawPtr(),weights.getRawPtr());
00195       for (int j=0; j<=2*i-1; j++) {
00196         if (j%2)
00197           analyticInt = 0.0;
00198         else 
00199           analyticInt = 2.0/((long double)j+1.0);
00200         testInt = evalQuad(i,j,nodes.getRawPtr(),weights.getRawPtr());
00201         long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) );
00202         long double absdiff = std::fabs(analyticInt - testInt);
00203         *outStream << "Cubature order " << std::setw(2) << std::left << i << " integrating "
00204                    << "x^" << std::setw(2) << std::left << j << ":" << "   "
00205                    << std::scientific << std::setprecision(16) << testInt << "   " 
00206                    << analyticInt << "   " << std::setprecision(4) << absdiff << "   " << "<?" 
00207                    << "   " << abstol << "\n";
00208         if (absdiff > abstol) {
00209           errorFlag++;
00210           *outStream << std::right << std::setw(111) << "^^^^---FAILURE!\n";
00211         }
00212       }
00213     }
00214     *outStream << "\n";    
00215 
00216     *outStream << "Clenshaw-Curtis Cubature \n";
00217     *outStream << "Integrates functions on [-1,1] weighted by w(x) = 1\n";
00218     for (int i = 1; i<=maxOrder; i++) {
00219       Teuchos::Array<long double> nodes(i), weights(i);
00220       IntrepidBurkardtRules::clenshaw_curtis_compute(i,nodes.getRawPtr(),weights.getRawPtr());
00221       for (int j=0; j<i; j++) {
00222         if (j%2)
00223           analyticInt = 0.0;
00224         else 
00225           analyticInt = 2.0/((long double)j+1.0);
00226         testInt = evalQuad(i,j,nodes.getRawPtr(),weights.getRawPtr());
00227         long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) );
00228         long double absdiff = std::fabs(analyticInt - testInt);
00229         *outStream << "Cubature order " << std::setw(2) << std::left << i << " integrating "
00230                    << "x^" << std::setw(2) << std::left << j << ":" << "   "
00231                    << std::scientific << std::setprecision(16) << testInt << "   " 
00232                    << analyticInt << "   " << std::setprecision(4) << absdiff << "   " << "<?" 
00233                    << "   " << abstol << "\n";
00234         if (absdiff > abstol) {
00235           errorFlag++;
00236           *outStream << std::right << std::setw(111) << "^^^^---FAILURE!\n";
00237         }
00238       }
00239     }
00240     *outStream << "\n";    
00241 
00242     *outStream << "Fejer Type 2 Cubature \n";
00243     *outStream << "Integrates functions on [-1,1] weighted by w(x) = 1\n";
00244     for (int i = 1; i<=maxOrder; i++) {
00245       Teuchos::Array<long double> nodes(i), weights(i);
00246       IntrepidBurkardtRules::fejer2_compute(i,nodes.getRawPtr(),weights.getRawPtr());
00247       for (int j=0; j<i; j++) {
00248         if (j%2)
00249           analyticInt = 0.0;
00250         else
00251           analyticInt = 2.0/((long double)j+1.0);
00252         testInt = evalQuad(i,j,nodes.getRawPtr(),weights.getRawPtr());
00253         long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) );
00254         long double absdiff = std::fabs(analyticInt - testInt);
00255         *outStream << "Cubature order " << std::setw(2) << std::left << i << " integrating "
00256                    << "x^" << std::setw(2) << std::left << j << ":" << "   "
00257                    << std::scientific << std::setprecision(16) << testInt << "   " 
00258                    << analyticInt << "   " << std::setprecision(4) << absdiff << "   " << "<?" 
00259                    << "   " << abstol << "\n";
00260         if (absdiff > abstol) {
00261           errorFlag++;
00262           *outStream << std::right << std::setw(111) << "^^^^---FAILURE!\n";
00263         }
00264       }
00265     }
00266     *outStream << "\n";    
00267 
00268     *outStream << "Gauss-Patterson Cubature \n";
00269     *outStream << "Integrates functions on [-1,1] weighted by w(x) = 1\n";
00270     for (int l = 1; l<=7; l++) {
00271       int i = (int)pow(2.0,(double)l+1.0)-1;
00272       Teuchos::Array<long double> nodes(i), weights(i);
00273       IntrepidBurkardtRules::patterson_lookup(i,nodes.getRawPtr(),weights.getRawPtr());
00274       for (int j=0; j<=(1.5*i+0.5); j++) {
00275         if (j%2)
00276           analyticInt = 0.0;
00277         else
00278           analyticInt = 2.0/((long double)j+1.0);
00279         testInt = evalQuad(i,j,nodes.getRawPtr(),weights.getRawPtr());
00280         long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) );
00281         long double absdiff = std::fabs(analyticInt - testInt);
00282         *outStream << "Cubature order " << std::setw(2) << std::left << i << " integrating "
00283                    << "x^" << std::setw(2) << std::left << j << ":" << "   "
00284                    << std::scientific << std::setprecision(16) << testInt << "   " 
00285                    << analyticInt << "   " << std::setprecision(4) << absdiff << "   " << "<?" 
00286                    << "   " << abstol << "\n";
00287         if (absdiff > abstol) {
00288           errorFlag++;
00289           *outStream << std::right << std::setw(111) << "^^^^---FAILURE!\n";
00290         }
00291       }
00292     }
00293     *outStream << "\n";
00294 
00295     *outStream << "Gauss-Chebyshev Type 1 Cubature \n";
00296     *outStream << "Integrates functions on [-1,1] weighted by w(x) = 1/sqrt(1-x^2)\n";
00297     for (int i = 1; i<=maxOrder; i++) {
00298       Teuchos::Array<long double> nodes(i), weights(i);
00299       IntrepidBurkardtRules::chebyshev1_compute(i,nodes.getRawPtr(),weights.getRawPtr());      
00300       for (int j=0; j<=2*i-1; j++) {
00301         analyticInt = chebyshev1<long double>(j);
00302         testInt = evalQuad(i,j,nodes.getRawPtr(),weights.getRawPtr());
00303         long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) );
00304         long double absdiff = std::fabs(analyticInt - testInt);
00305         *outStream << "Cubature order " << std::setw(2) << std::left << i << " integrating "
00306                    << "x^" << std::setw(2) << std::left << j << ":" << "   "
00307                    << std::scientific << std::setprecision(16) << testInt << "   " 
00308                    << analyticInt << "   " << std::setprecision(4) << absdiff << "   " << "<?" 
00309                    << "   " << abstol << "\n";
00310         if (absdiff > abstol) {
00311           errorFlag++;
00312           *outStream << std::right << std::setw(111) << "^^^^---FAILURE!\n";
00313         }
00314       }
00315     }
00316     *outStream << "\n";
00317 
00318     *outStream << "Gauss-Chebyshev Type 2 Cubature \n";
00319     *outStream << "Integrates functions on [-1,1] weighted by w(x) = sqrt(1-x^2)\n";
00320     for (int i = 1; i<=maxOrder; i++) {
00321       Teuchos::Array<long double> nodes(i), weights(i);
00322       IntrepidBurkardtRules::chebyshev2_compute(i,nodes.getRawPtr(),weights.getRawPtr());      
00323       for (int j=0; j<=2*i-1; j++) {
00324         analyticInt = chebyshev2<long double>(j);
00325         testInt = evalQuad(i,j,nodes.getRawPtr(),weights.getRawPtr());
00326         long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) );
00327         long double absdiff = std::fabs(analyticInt - testInt);
00328         *outStream << "Cubature order " << std::setw(2) << std::left << i << " integrating "
00329                    << "x^" << std::setw(2) << std::left << j << ":" << "   "
00330                    << std::scientific << std::setprecision(16) << testInt << "   " 
00331                    << analyticInt << "   " << std::setprecision(4) << absdiff << "   " << "<?" 
00332                    << "   " << abstol << "\n";
00333         if (absdiff > abstol) {
00334           errorFlag++;
00335           *outStream << std::right << std::setw(111) << "^^^^---FAILURE!\n";
00336         }
00337       }
00338     }
00339     *outStream << "\n";
00340     
00341      *outStream << "Gauss-Laguerre Cubature \n";
00342     *outStream << "Integrates functions on [0,oo) weighted by w(x) = exp(-x)\n";
00343     for (int i = 1; i<=maxOrder; i++) {
00344       Teuchos::Array<long double> nodes(i), weights(i);
00345       IntrepidBurkardtRules::laguerre_compute(i,nodes.getRawPtr(),weights.getRawPtr());
00346        for (int j=0; j<=2*i-1; j++) {
00347         analyticInt = tgammal((long double)(j+1));
00348         testInt = evalQuad(i,j,nodes.getRawPtr(),weights.getRawPtr());
00349         long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) );
00350         long double absdiff = std::fabs(analyticInt - testInt);
00351         *outStream << "Cubature order " << std::setw(2) << std::left << i << " integrating "
00352                    << "x^" << std::setw(2) << std::left << j << ":" << "   "
00353                    << std::scientific << std::setprecision(16) << testInt << "   " 
00354                    << analyticInt << "   " << std::setprecision(4) << absdiff << "   " << "<?" 
00355                    << "   " << abstol << "\n";
00356         if (absdiff > abstol) {
00357           errorFlag++;
00358           *outStream << std::right << std::setw(111) << "^^^^---FAILURE!\n";
00359         }
00360       }
00361     }
00362     *outStream << "\n";
00363 
00364     maxOrder = 10;
00365 
00366     *outStream << "Gauss-Hermite Cubature \n";
00367     *outStream << "Integrates functions on (-oo,oo) weighted by w(x) = exp(-x^2)\n";
00368     for (int i = 1; i<=maxOrder; i++) {
00369       Teuchos::Array<long double> nodes(i), weights(i);
00370       IntrepidBurkardtRules::hermite_compute(i,nodes.getRawPtr(),
00371                                              weights.getRawPtr());      
00372       for (int j=0; j<=2*i-1; j++) { 
00373         if (j%2)
00374           analyticInt = 0.0;
00375         else
00376           analyticInt = factorial2<long double>(j-1)*sqrt(M_PI)/powl(2.0,(long double)j/2.0);
00377         testInt = evalQuad(i,j,nodes.getRawPtr(),weights.getRawPtr());
00378         long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) );
00379         long double absdiff = std::fabs(analyticInt - testInt);
00380         *outStream << "Cubature order " << std::setw(2) << std::left << i 
00381                    << " integrating "
00382                    << "x^" << std::setw(2) << std::left << j << ":" 
00383                    << "   "
00384                    << std::scientific << std::setprecision(16) << testInt 
00385                    << "   " 
00386                    << analyticInt << "   " << std::setprecision(4) 
00387                    << absdiff << "   " << "<?" 
00388                    << "   " << abstol << "\n";
00389         if (absdiff > abstol) {
00390           errorFlag++;
00391           *outStream << std::right << std::setw(111) << "^^^^---FAILURE!\n";
00392         }
00393       }
00394     }
00395     *outStream << "\n";    
00396 
00397     reltol = 1e-6;
00398 
00399     *outStream << "Hermite-Genz-Keister Cubature \n";
00400     *outStream << "Integrates functions on (-oo,oo) weighted by w(x) = exp(-x^2)\n";
00401     int order[4] = {1,3, 9,19};
00402     int max[4]   = {1,5,15,29};
00403     for (int l = 0; l<4; l++) {
00404       int i = order[l];
00405       int m = max[l];
00406       Teuchos::Array<long double> nodes(i), weights(i);
00407       IntrepidBurkardtRules::hermite_genz_keister_lookup(i,nodes.getRawPtr(),
00408                                                          weights.getRawPtr());  
00409       for (int j=0; j<=m; j++) { 
00410         if (j%2)
00411           analyticInt = 0.0;
00412         else
00413           analyticInt = factorial2<long double>(j-1)*sqrt(M_PI)/powl(2.0,(long double)j/2.0);
00414         if (i>=36)
00415           analyticInt /= sqrt(M_PI);
00416         testInt = evalQuad(i,j,nodes.getRawPtr(),weights.getRawPtr());
00417         long double abstol = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) );
00418         long double absdiff = std::fabs(analyticInt - testInt);
00419         *outStream << "Cubature order " << std::setw(2) << std::left << i 
00420                    << " integrating "
00421                    << "x^" << std::setw(2) << std::left << j << ":" << "   "
00422                    << std::scientific << std::setprecision(16) << testInt 
00423                    << "   " 
00424                    << analyticInt << "   " << std::setprecision(4) << absdiff 
00425                    << "   " << "<?" 
00426                    << "   " << abstol << "\n";
00427         if (absdiff > abstol) {
00428           errorFlag++;
00429           *outStream << std::right << std::setw(111) << "^^^^---FAILURE!\n";
00430         }
00431       }
00432     }
00433     *outStream << "\n";
00434   }
00435   catch (std::logic_error err) {
00436     *outStream << err.what() << "\n";
00437     errorFlag = -1;
00438   };
00439 
00440 
00441   if (errorFlag != 0)
00442     std::cout << "End Result: TEST FAILED\n";
00443   else
00444     std::cout << "End Result: TEST PASSED\n";
00445 
00446   // reset format state of std::cout
00447   std::cout.copyfmt(oldFormatState);
00448 
00449   return errorFlag;
00450 }
00451