Intrepid
http://trilinos.sandia.gov/packages/docs/r10.10/packages/intrepid/test/Discretization/Integration/test_11.cpp
Go to the documentation of this file.
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 
00044 
00051 //#include "Intrepid_CubatureLineSorted.hpp"
00052 #include "Intrepid_CubatureLineSorted.hpp"
00053 #include "Intrepid_Utils.hpp"
00054 #include "Teuchos_oblackholestream.hpp"
00055 #include "Teuchos_RCP.hpp"
00056 #include "Teuchos_GlobalMPISession.hpp"
00057 
00058 using namespace Intrepid;
00059 
00060 /*
00061   Computes integrals of monomials over a given reference cell.
00062 */
00063 long double evalQuad(int order, int power, EIntrepidBurkardt rule) {
00064  
00065   CubatureLineSorted<long double> lineCub(rule,order,false);  
00066   int size = lineCub.getNumPoints();
00067   FieldContainer<long double> cubPoints(size);
00068   FieldContainer<long double> cubWeights(size);
00069   lineCub.getCubature(cubPoints,cubWeights);
00070 
00071   int mid  = size/2;
00072   long double Q = 0.0;
00073   if (size%2) 
00074     Q = cubWeights(mid)*powl(cubPoints(mid),power);
00075 
00076   for (int i=0; i<mid; i++) {
00077     Q += cubWeights(i)*powl(cubPoints(i),power)+
00078       cubWeights(size-i-1)*powl(cubPoints(size-i-1),power);
00079   }
00080   return Q;
00081 }
00082 
00083 long double evalInt(int power, EIntrepidBurkardt rule) {
00084   double I = 0;
00085 
00086   if (rule==BURK_CLENSHAWCURTIS||rule==BURK_FEJER2||
00087       rule==BURK_LEGENDRE||rule==BURK_PATTERSON || 
00088       rule==BURK_TRAPEZOIDAL) {
00089     if (power%2)
00090       I = 0.0;
00091     else 
00092       I = 2.0/((long double)power+1.0);
00093   }
00094   else if (rule==BURK_LAGUERRE) {
00095     I = tgammal((long double)(power+1));
00096   }
00097   else if (rule==BURK_CHEBYSHEV1) {
00098     long double bot, top;
00099     if (!(power%2)) {
00100       top = 1; bot = 1;
00101       for (int i=2;i<=power;i+=2) {
00102         top *= (long double)(i-1);
00103         bot *= (long double)i;
00104       }
00105       I = M_PI*top/bot;
00106     }
00107     else {
00108       I = 0.0;
00109     }
00110   }
00111   else if (rule==BURK_CHEBYSHEV2) {
00112     long double bot, top;
00113     if (!(power%2)) {
00114       top = 1; bot = 1;
00115       for (int i=2;i<=power;i+=2) {
00116         top *= (long double)(i-1);
00117         bot *= (long double)i;
00118       }
00119       bot *= (long double)(power+2);
00120       I    = M_PI*top/bot;
00121     }
00122     else {
00123       I = 0.0;
00124     }
00125   }
00126   else if (rule==BURK_HERMITE||rule==BURK_GENZKEISTER) {
00127     if (power%2) {
00128       I = 0.0;
00129     }
00130     else {  
00131       long double value = 1.0;
00132       if ((power-1)>=1) {
00133         int n_copy = power-1;
00134         while (1<n_copy) {
00135           value  *= (long double)n_copy;
00136           n_copy -= 2;
00137         }
00138       }
00139       I = value*sqrt(M_PI)/powl(2.0,(long double)power/2.0);
00140     }
00141   }
00142   return I;
00143 }
00144 
00145 int main(int argc, char *argv[]) {
00146   
00147   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00148 
00149   // This little trick lets us print to std::cout only if
00150   // a (dummy) command-line argument is provided.
00151   int iprint     = argc - 1;
00152   Teuchos::RCP<std::ostream> outStream;
00153   Teuchos::oblackholestream bhs; // outputs nothing
00154   if (iprint > 0)
00155     outStream = Teuchos::rcp(&std::cout, false);
00156   else
00157     outStream = Teuchos::rcp(&bhs, false);
00158 
00159   // Save the format state of the original std::cout.
00160   Teuchos::oblackholestream oldFormatState;
00161   oldFormatState.copyfmt(std::cout);
00162  
00163   *outStream \
00164   << "===============================================================================\n" \
00165   << "|                                                                             |\n" \
00166   << "|                         Unit Test (CubatureLineSorted)                      |\n" \
00167   << "|                                                                             |\n" \
00168   << "|     1) Computing integrals of monomials in 1D                               |\n" \
00169   << "|                                                                             |\n" \
00170   << "|  Questions? Contact  Drew Kouri (dpkouri@sandia.gov) or                     |\n" \
00171   << "|                      Denis Ridzal (dridzal@sandia.gov).                     |\n" \
00172   << "|                                                                             |\n" \
00173   << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00174   << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00175   << "|                                                                             |\n" \
00176   << "===============================================================================\n"\
00177   << "| TEST 11: integrals of monomials in 1D                                       |\n"\
00178   << "===============================================================================\n";
00179 
00180   // internal variables:
00181   int         errorFlag   = 0;
00182   long double reltol      = 1.0e+05*INTREPID_TOL;
00183   int         maxDeg      = 0;    
00184   long double analyticInt = 0;
00185   long double testInt     = 0;
00186   int         maxOrder    = 30;
00187 
00188   *outStream << "\nIntegrals of monomials on a reference line (edge):\n";
00189   // compute and compare integrals
00190   try {
00191     for (EIntrepidBurkardt rule=BURK_CHEBYSHEV1; rule <= BURK_LAGUERRE; rule++) {
00192       *outStream << "Testing " << EIntrepidBurkardtToString(rule) << "\n";
00193       // compute integrals
00194       if (rule==BURK_HERMITE)
00195         maxOrder = 10;
00196       else if (rule==BURK_TRAPEZOIDAL) 
00197         maxOrder = 2;
00198       else 
00199         maxOrder = 30;
00200 
00201       if (rule!=BURK_PATTERSON&&rule!=BURK_GENZKEISTER) {
00202         for (int i=1; i <= maxOrder; i++) {
00203           if ( rule==BURK_CHEBYSHEV1 ||
00204                rule==BURK_CHEBYSHEV2 ||
00205                rule==BURK_LEGENDRE   ||
00206                rule==BURK_LAGUERRE   ||
00207                rule==BURK_HERMITE      )
00208             maxDeg = 2*i-1;
00209           else if ( rule==BURK_CLENSHAWCURTIS ||
00210                     rule==BURK_FEJER2         ||
00211                     rule==BURK_TRAPEZOIDAL      ) 
00212             maxDeg = i-1;
00213           
00214           for (int j=0; j <= maxDeg; j++) {
00215             analyticInt = evalInt(j,rule);
00216             testInt     = evalQuad(i,j,rule);
00217             
00218             long double abstol  = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) );
00219             long double absdiff = std::fabs(analyticInt - testInt);
00220             *outStream << "Cubature order " << std::setw(2) << std::left 
00221                        << i << " integrating "
00222                        << "x^" << std::setw(2) << std::left << j <<  ":" 
00223                        << "   "
00224                        << std::scientific << std::setprecision(16) << testInt 
00225                        << "   " << analyticInt 
00226                        << "   " << std::setprecision(4) << absdiff << "   " 
00227                        << "<?" << "   " << abstol 
00228                        << "\n";
00229             if (absdiff > abstol) {
00230               errorFlag++;
00231               *outStream << std::right << std::setw(104) 
00232                          << "^^^^---FAILURE!\n";
00233             }
00234           } // end for j
00235           *outStream << "\n";
00236         } // end for i
00237       }
00238       else if (rule==BURK_PATTERSON) {
00239         for (int i=0; i < 8; i++) {
00240           int l = (int)std::pow(2.0,(double)i+1.0)-1;
00241           if (i==0) 
00242             maxDeg = 1;
00243           else
00244             maxDeg = (int)(1.5*(double)l+0.5);
00245           for (int j=0; j <= maxDeg; j++) {
00246             analyticInt = evalInt(j,rule);
00247             testInt     = evalQuad(l,j,rule);
00248             
00249             long double abstol  = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) );
00250             long double absdiff = std::fabs(analyticInt - testInt);
00251             *outStream << "Cubature order " << std::setw(2) << std::left 
00252                        << l << " integrating "
00253                        << "x^" << std::setw(2) << std::left << j <<  ":" 
00254                        << "   "
00255                        << std::scientific << std::setprecision(16) << testInt 
00256                        << "   " << analyticInt 
00257                        << "   " << std::setprecision(4) << absdiff << "   " 
00258                        << "<?" << "   " << abstol 
00259                        << "\n";
00260             if (absdiff > abstol) {
00261               errorFlag++;
00262               *outStream << std::right << std::setw(104) 
00263                          << "^^^^---FAILURE!\n";
00264             }
00265           } // end for j
00266           *outStream << "\n";
00267         } // end for i
00268       }
00269       else if (rule==BURK_GENZKEISTER) {
00270         reltol *= 1.0e+02;
00271         int o_ghk[4] = {1,3, 9,19};
00272         int p_ghk[4] = {1,5,15,29};
00273         for (int i=0; i < 4; i++) {
00274           int l  = o_ghk[i];
00275           maxDeg = p_ghk[i];
00276           for (int j=0; j <= maxDeg; j++) {
00277             analyticInt = evalInt(j,rule);
00278             testInt     = evalQuad(l,j,rule);
00279             
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 
00283                        << l << " integrating "
00284                        << "x^" << std::setw(2) << std::left << j <<  ":" 
00285                        << "   "
00286                        << std::scientific << std::setprecision(16) << testInt 
00287                        << "   " << analyticInt 
00288                        << "   " << std::setprecision(4) << absdiff << "   " 
00289                        << "<?" << "   " << abstol 
00290                        << "\n";
00291             if (absdiff > abstol) {
00292               errorFlag++;
00293               *outStream << std::right << std::setw(104) << "^^^^---FAILURE!\n";
00294             }
00295           } // end for j
00296           *outStream << "\n";
00297         } // end for i
00298       }
00299     } // end for rule
00300   }
00301   catch (std::logic_error err) {
00302     *outStream << err.what() << "\n";
00303     errorFlag = -1;
00304   };
00305 
00306 
00307   if (errorFlag != 0)
00308     std::cout << "End Result: TEST FAILED\n";
00309   else
00310     std::cout << "End Result: TEST PASSED\n";
00311 
00312   // reset format state of std::cout
00313   std::cout.copyfmt(oldFormatState);
00314 
00315   return errorFlag;
00316 }