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