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