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