Intrepid
http://trilinos.sandia.gov/packages/docs/r10.10/packages/intrepid/test/Discretization/Integration/test_15.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_CubatureLineSorted.hpp"
00052 #include "Intrepid_Utils.hpp"
00053 #include "Teuchos_oblackholestream.hpp"
00054 #include "Teuchos_RCP.hpp"
00055 #include "Teuchos_GlobalMPISession.hpp"
00056 
00057 using namespace Intrepid;
00058 
00059 /*
00060   Computes integrals of monomials over a given reference cell.
00061 */
00062 long double evalQuad(int order, int power, EIntrepidBurkardt rule) {
00063  
00064   CubatureLineSorted<long double> lineCub(rule,order,false);
00065   CubatureLineSorted<long double> lineCub2(rule,order-1,false);
00066   lineCub.update(-1.0,lineCub2,1.0);
00067   int size = lineCub.getNumPoints();
00068   FieldContainer<long double> cubPoints(size);
00069   FieldContainer<long double> cubWeights(size);
00070   lineCub.getCubature(cubPoints,cubWeights);
00071 
00072   long double Q = 0.0;
00073   for (int i=0; i<size; i++) {
00074     Q += cubWeights(i)*powl(cubPoints(i),(long double)power);
00075   }
00076   return Q;
00077 
00078   /*
00079   int mid  = size/2;
00080   long double Q = 0.0;
00081   if (size%2) 
00082     Q = cubWeights(mid)*powl(cubPoints(mid),power);
00083 
00084   for (int i=0; i<mid; i++) {
00085     Q += cubWeights(i)*powl(cubPoints(i),power)+
00086       cubWeights(order-i-1)*powl(cubPoints(size-i-1),power);
00087   }
00088   return Q;
00089   */
00090 }
00091 
00092 int main(int argc, char *argv[]) {
00093   
00094   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00095 
00096   // This little trick lets us print to std::cout only if
00097   // a (dummy) command-line argument is provided.
00098   int iprint     = argc - 1;
00099   Teuchos::RCP<std::ostream> outStream;
00100   Teuchos::oblackholestream bhs; // outputs nothing
00101   if (iprint > 0)
00102     outStream = Teuchos::rcp(&std::cout, false);
00103   else
00104     outStream = Teuchos::rcp(&bhs, false);
00105 
00106   // Save the format state of the original std::cout.
00107   Teuchos::oblackholestream oldFormatState;
00108   oldFormatState.copyfmt(std::cout);
00109  
00110   *outStream \
00111   << "===============================================================================\n" \
00112   << "|                                                                             |\n" \
00113   << "|                         Unit Test (CubatureLineSorted)                      |\n" \
00114   << "|                                                                             |\n" \
00115   << "|     1) Computing differential integrals of monomials in 1D                  |\n" \
00116   << "|                                                                             |\n" \
00117   << "|  Questions? Contact  Drew Kouri (dpkouri@sandia.gov) or                     |\n" \
00118   << "|                      Denis Ridzal (dridzal@sandia.gov).                     |\n" \
00119   << "|                                                                             |\n" \
00120   << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00121   << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00122   << "|                                                                             |\n" \
00123   << "===============================================================================\n"\
00124   << "| TEST 15: differential integrals of monomials in 1D                          |\n"\
00125   << "===============================================================================\n";
00126 
00127   // internal variables:
00128   int         errorFlag   = 0;
00129   long double abstol      = 1.0e+05*INTREPID_TOL;
00130   int         maxDeg      = 0;    
00131   long double analyticInt = 0;
00132   long double testInt     = 0;
00133   int         maxOrder    = 60;
00134   EIntrepidBurkardt rule  = BURK_LEGENDRE;
00135 
00136   *outStream << "\nDifferential integrals of monomials on a reference line:\n";
00137   // compute and compare integrals
00138   try {
00139     *outStream << "Testing " << EIntrepidBurkardtToString(rule) << "\n";
00140     // compute integrals
00141     for (int i=2; i <= maxOrder; i++) {
00142       maxDeg = 2*(i-1)-1;  
00143       for (int j=0; j <= maxDeg; j++) {
00144         testInt             = evalQuad(i,j,rule);
00145         long double absdiff = std::fabs(testInt);
00146         *outStream << "Cubature order " << std::setw(2) << std::left << i 
00147                    << " integrating "
00148                    << "x^" << std::setw(2) << std::left << j <<  ":" << "   "
00149                    << std::scientific << std::setprecision(16) << testInt 
00150                    << "   " << analyticInt 
00151                    << "   " << std::setprecision(4) << absdiff << "   " 
00152                    << "<?" << "   " << abstol 
00153                    << "\n";
00154         if (absdiff > abstol) {
00155           errorFlag++;
00156           *outStream << std::right << std::setw(104) << "^^^^---FAILURE!\n";
00157         }
00158       } // end for j
00159       *outStream << "\n";
00160     } // end for i
00161   }
00162   catch (std::logic_error err) {
00163     *outStream << err.what() << "\n";
00164     errorFlag = -1;
00165   };
00166 
00167 
00168   if (errorFlag != 0)
00169     std::cout << "End Result: TEST FAILED\n";
00170   else
00171     std::cout << "End Result: TEST PASSED\n";
00172 
00173   // reset format state of std::cout
00174   std::cout.copyfmt(oldFormatState);
00175 
00176   return errorFlag;
00177 }