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