Intrepid
http://trilinos.sandia.gov/packages/docs/r10.12/packages/intrepid/test/Discretization/Integration/test_22.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_AdaptiveSparseGrid.hpp"
00053 #include "Intrepid_Utils.hpp"
00054 #include "Teuchos_oblackholestream.hpp"
00055 #include "Teuchos_RCP.hpp"
00056 #include "Teuchos_RefCountPtr.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, int dimension, int order, 
00065                      std::vector<EIntrepidBurkardt> rule,
00066                      std::vector<EIntrepidGrowth> growth) {
00067  
00068   CubatureTensorSorted<long double> lineCub(0,dimension);
00069   AdaptiveSparseGrid<long double,std::vector<long double> >::buildSparseGrid(
00070                                                    lineCub,dimension,
00071                                                    order,rule,
00072                                                    growth,false);
00073 
00074   int size = lineCub.getNumPoints();
00075   FieldContainer<long double> cubPoints(size,dimension);
00076   FieldContainer<long double> cubWeights(size);
00077   lineCub.getCubature(cubPoints,cubWeights);
00078 
00079   int mid  = size/2;
00080   long double Q = 0.0;
00081   if (size%2) {
00082     Q  = cubWeights(mid);
00083     for (int i=0; i<dimension; i++) {
00084       Q *= powl(cubPoints(mid,i),(long double)power[i]);
00085     }
00086   }
00087 
00088   for (int i=0; i<mid; i++) {
00089     long double value1 = cubWeights(i);
00090     long double value2 = cubWeights(size-i-1);
00091     for (int j=0; j<dimension; j++) {
00092       value1 *= powl(cubPoints(i,j),(long double)power[j]);
00093       value2 *= powl(cubPoints(size-i-1,j),(long double)power[j]);
00094     }
00095     Q += value1+value2;
00096   }
00097   return Q;
00098 }
00099 
00100 long double evalInt(int dimension, std::vector<int> power, 
00101                     std::vector<EIntrepidBurkardt> rule) {
00102   long double I = 1.0;
00103 
00104   for (int i=0; i<dimension; i++) {
00105     if (rule[i]==BURK_CLENSHAWCURTIS||rule[i]==BURK_FEJER2||
00106         rule[i]==BURK_LEGENDRE||rule[i]==BURK_PATTERSON || 
00107         rule[i]==BURK_TRAPEZOIDAL) {
00108       if (power[i]%2)
00109         I *= 0.0;
00110       else 
00111         I *= 2.0/((long double)power[i]+1.0);
00112     }
00113     else if (rule[i]==BURK_LAGUERRE) {
00114       I *= tgammal((long double)(power[i]+1));
00115     }
00116     else if (rule[i]==BURK_CHEBYSHEV1) {
00117       long double bot, top;
00118       if (!(power[i]%2)) {
00119         top = 1; bot = 1;
00120         for (int j=2;j<=power[i];j+=2) {
00121           top *= (long double)(j-1);
00122           bot *= (long double)j;
00123         }
00124         I *= M_PI*top/bot;
00125       }
00126       else {
00127         I *= 0.0;
00128       }
00129     }
00130     else if (rule[i]==BURK_CHEBYSHEV2) {
00131       long double bot, top;
00132       if (!(power[i]%2)) {
00133       top = 1; bot = 1;
00134       for (int j=2;j<=power[i];j+=2) {
00135         top *= (long double)(j-1);
00136         bot *= (long double)j;
00137       }
00138       bot *= (long double)(power[i]+2);
00139       I   *= M_PI*top/bot;
00140       }
00141       else {
00142         I *= 0.0;
00143       }
00144     }
00145     else if (rule[i]==BURK_HERMITE||rule[i]==BURK_GENZKEISTER) {
00146       if (power[i]%2) {
00147         I *= 0.0;
00148       }
00149       else {  
00150         long double value = 1.0;
00151         if ((power[i]-1)>=1) {
00152           int n_copy = power[i]-1;
00153           while (1<n_copy) {
00154             value  *= (long double)n_copy;
00155             n_copy -= 2;
00156           }
00157         }
00158         I *= value*sqrt(M_PI)/powl(2.0,(long double)power[i]/2.0);
00159       }
00160     }
00161   }
00162   return I;
00163 }
00164 
00165 int main(int argc, char *argv[]) {
00166   
00167   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00168 
00169   // This little trick lets us print to std::cout only if
00170   // a (dummy) command-line argument is provided.
00171   int iprint     = argc - 1;
00172   Teuchos::RCP<std::ostream> outStream;
00173   Teuchos::oblackholestream bhs; // outputs nothing
00174   if (iprint > 0)
00175     outStream = Teuchos::rcp(&std::cout, false);
00176   else
00177     outStream = Teuchos::rcp(&bhs, false);
00178 
00179   // Save the format state of the original std::cout.
00180   Teuchos::oblackholestream oldFormatState;
00181   oldFormatState.copyfmt(std::cout);
00182  
00183   *outStream \
00184   << "===============================================================================\n" \
00185   << "|                                                                             |\n" \
00186   << "|                         Unit Test (CubatureTensorSorted)                    |\n" \
00187   << "|                                                                             |\n" \
00188   << "|     1) Computing integrals of monomials in 2D                               |\n" \
00189   << "|                                                                             |\n" \
00190   << "|  Questions? Contact  Drew Kouri (dpkouri@sandia.gov) or                     |\n" \
00191   << "|                      Denis Ridzal (dridzal@sandia.gov).                     |\n" \
00192   << "|                                                                             |\n" \
00193   << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00194   << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00195   << "|                                                                             |\n" \
00196   << "===============================================================================\n"\
00197   << "| TEST 22: integrals of monomials in 2D - Anisotropic but no growth rules     |\n"\
00198   << "===============================================================================\n";
00199 
00200 
00201   // internal variables:
00202   int dimension           = 2;
00203   int         errorFlag   = 0;
00204   long double reltol      = 1.0e+02*INTREPID_TOL;
00205   int         maxDeg      = 0;    
00206   long double analyticInt = 0;
00207   long double testInt     = 0;
00208   int         maxOrder    = 6;
00209   std::vector<int> power(2,0);
00210   std::vector<EIntrepidBurkardt> rule1(2,BURK_CLENSHAWCURTIS);
00211   std::vector<EIntrepidGrowth> growth1(2,GROWTH_FULLEXP);
00212 
00213   *outStream << "\nIntegrals of monomials on a reference line (edge):\n";
00214   // compute and compare integrals
00215   try {
00216     for (int i=0; i<=maxOrder; i++) {
00217       maxDeg = i-1;
00218       for (int j=0; j <= maxDeg; j++) {
00219         power[0] = j;
00220         for (int k=0; k <= maxDeg; k++) {
00221           power[1] = k;
00222           if (j+k < maxDeg) {
00223             analyticInt = evalInt(dimension, power, rule1);
00224             testInt     = evalQuad(power,dimension,maxOrder,rule1,growth1);
00225             
00226             long double abstol  = (analyticInt == 0.0 ? reltol : std::fabs(reltol*analyticInt) );
00227             long double absdiff = std::fabs(analyticInt - testInt);
00228             *outStream << "Cubature order " << std::setw(2) << std::left << i 
00229                        << " integrating " << "x^" << std::setw(2) 
00230                        << std::left << j << "y^" << std::setw(2) << std::left 
00231                        << k <<  ":" << "   " << std::scientific 
00232                        << std::setprecision(16) << testInt 
00233                        << "   " << analyticInt << "   " 
00234                        << std::setprecision(4) << absdiff << "   " 
00235                        << "<?" << "   " << abstol << "\n";
00236             if (absdiff > abstol) {
00237               errorFlag++;
00238               *outStream << std::right << std::setw(104) << "^^^^---FAILURE!\n";
00239             }
00240           }
00241         } // end for k
00242         *outStream << "\n";
00243       } // end for j
00244       *outStream << "\n";
00245     } // end for i
00246   }
00247   catch (std::logic_error err) {
00248     *outStream << err.what() << "\n";
00249     errorFlag = -1;
00250   };
00251 
00252 
00253   if (errorFlag != 0)
00254     std::cout << "End Result: TEST FAILED\n";
00255   else
00256     std::cout << "End Result: TEST PASSED\n";
00257 
00258   // reset format state of std::cout
00259   std::cout.copyfmt(oldFormatState);
00260 
00261   return errorFlag;
00262 }