Intrepid
http://trilinos.sandia.gov/packages/docs/r10.12/packages/intrepid/test/Discretization/Integration/test_09.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_CubatureGenSparse.hpp"
00052 #include "Intrepid_Utils.hpp"
00053 #include "Teuchos_oblackholestream.hpp"
00054 #include "Teuchos_RCP.hpp"
00055 #include "Teuchos_BLAS.hpp"
00056 #include "Teuchos_GlobalMPISession.hpp"
00057 
00058 using namespace Intrepid;
00059 
00060 
00061 /*
00062   Monomial evaluation.
00063     in 1D, for point p(x)    : x^xDeg
00064     in 2D, for point p(x,y)  : x^xDeg * y^yDeg
00065     in 3D, for point p(x,y,z): x^xDeg * y^yDeg * z^zDeg
00066 */
00067 double computeMonomial(FieldContainer<double> & p, int xDeg, int yDeg=0, int zDeg=0) {
00068   double val = 1.0;
00069   int polydeg[3];
00070   polydeg[0] = xDeg; polydeg[1] = yDeg; polydeg[2] = zDeg;
00071   for (int i=0; i<p.dimension(0); i++) {
00072     val *= std::pow(p(i),polydeg[i]);
00073   }
00074   return val;
00075 }
00076 
00077 
00078 /*
00079   Computes integrals of monomials over a given reference cell.
00080 */
00081 void computeIntegral(Teuchos::Array<double>& testIntFixDeg, int cubDegree) {
00082 
00083   CubatureGenSparse<double,3> myCub(cubDegree);
00084 
00085   int cubDim       = myCub.getDimension();
00086   int numCubPoints = myCub.getNumPoints();
00087   int numPolys     = (cubDegree+1)*(cubDegree+2)*(cubDegree+3)/6;
00088 
00089   FieldContainer<double> point(cubDim);
00090   FieldContainer<double> cubPoints(numCubPoints, cubDim);
00091   FieldContainer<double> cubWeights(numCubPoints);
00092   FieldContainer<double> functValues(numCubPoints, numPolys);
00093 
00094   myCub.getCubature(cubPoints, cubWeights);
00095 
00096   int polyCt = 0;
00097   for (int xDeg=0; xDeg <= cubDegree; xDeg++) {
00098     for (int yDeg=0; yDeg <= cubDegree-xDeg; yDeg++) {
00099       for (int zDeg=0; zDeg <= cubDegree-xDeg-yDeg; zDeg++) {
00100         for (int i=0; i<numCubPoints; i++) {
00101           for (int j=0; j<cubDim; j++) {
00102             point(j) = cubPoints(i,j);
00103           }
00104           functValues(i,polyCt) = computeMonomial(point, xDeg, yDeg, zDeg);
00105         }
00106         polyCt++;
00107       }
00108     }
00109   }
00110 
00111   Teuchos::BLAS<int, double> myblas;
00112   int inc = 1;
00113   double alpha = 1.0;
00114   double beta  = 0.0;
00115   myblas.GEMV(Teuchos::NO_TRANS, numPolys, numCubPoints, alpha, &functValues[0], numPolys,
00116               &cubWeights[0], inc, beta, &testIntFixDeg[0], inc);
00117 }
00118 
00119 
00120 int main(int argc, char *argv[]) {
00121 
00122   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00123 
00124   // This little trick lets us print to std::cout only if
00125   // a (dummy) command-line argument is provided.
00126   int iprint     = argc - 1;
00127   Teuchos::RCP<std::ostream> outStream;
00128   Teuchos::oblackholestream bhs; // outputs nothing
00129   if (iprint > 0)
00130     outStream = Teuchos::rcp(&std::cout, false);
00131   else
00132     outStream = Teuchos::rcp(&bhs, false);
00133 
00134   // Save the format state of the original std::cout.
00135   Teuchos::oblackholestream oldFormatState;
00136   oldFormatState.copyfmt(std::cout);
00137  
00138   *outStream \
00139   << "===============================================================================\n" \
00140   << "|                                                                             |\n" \
00141   << "|                      Unit Test (CubatureGenSparse)                          |\n" \
00142   << "|                                                                             |\n" \
00143   << "|     1) Computing integrals of monomials on reference cells in 3D            |\n" \
00144   << "|                         - using Level 2 BLAS -                              |\n" \
00145   << "|                                                                             |\n" \
00146   << "|  Questions? Contact  Pavel Bochev (pbboche@sandia.gov),                     |\n" \
00147   << "|                      Denis Ridzal (dridzal@sandia.gov) or                   |\n" \
00148   << "|                      Matthew Keegan (mskeega@sandia.gov).                   |\n" \
00149   << "|                                                                             |\n" \
00150   << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00151   << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00152   << "|                                                                             |\n" \
00153   << "===============================================================================\n"\
00154   << "| TEST 1: integrals of monomials in 3D (Level 2 BLAS version) using           |\n"\
00155   << "|           Generalized Sparse Grid Construction                              |\n"\
00156   << "===============================================================================\n";
00157 
00158   // internal variables:
00159   int                                      errorFlag = 0;
00160   int                                      polyCt = 0;
00161   int                                      offset = 0;
00162   Teuchos::Array< Teuchos::Array<double> > testInt;
00163   Teuchos::Array< Teuchos::Array<double> > analyticInt;
00164   Teuchos::Array<double>                   tmparray(1);
00165   double                                   reltol = 1.0e+04 * INTREPID_TOL;
00166   int maxDeg                             = 20; // can be as large as INTREPID_CUBATURE_SPARSE3D_GAUSS_MAX, but runtime is excessive
00167   int maxOffset                          = INTREPID_CUBATURE_LINE_GAUSS_MAX;
00168   int numPoly                            = (maxDeg+1)*(maxDeg+2)*(maxDeg+3)/6;
00169   int numAnalytic                        = (maxOffset+1)*(maxOffset+2)*(maxOffset+3)/6;
00170   testInt.assign(numPoly, tmparray);
00171   analyticInt.assign(numAnalytic, tmparray);
00172 
00173   // get names of files with analytic values
00174   std::string basedir = "./data";
00175   std::stringstream namestream;
00176   std::string filename;
00177   namestream << basedir << "/HEX_integrals" << ".dat";
00178   namestream >> filename;
00179 
00180   // format of data files with analytic values
00181   TypeOfExactData dataFormat = INTREPID_UTILS_FRACTION;
00182 
00183   // compute and compare integrals
00184   try {
00185       *outStream << "\nIntegrals of monomials:\n";
00186       std::ifstream filecompare(&filename[0]);
00187       // compute integrals
00188       for (int cubDeg=0; cubDeg <= maxDeg; cubDeg++) {
00189         int numMonomials = (cubDeg+1)*(cubDeg+2)*(cubDeg+3)/6; 
00190         testInt[cubDeg].resize(numMonomials);
00191         computeIntegral(testInt[cubDeg], cubDeg);
00192       }
00193       // get analytic values
00194       if (filecompare.is_open()) {
00195         getAnalytic(analyticInt, filecompare, dataFormat);
00196         // close file
00197         filecompare.close();
00198       }
00199       // perform comparison
00200       for (int cubDeg=0; cubDeg <= maxDeg; cubDeg++) {
00201         polyCt = 0;
00202         offset = 0;
00203         int oldErrorFlag = errorFlag;
00204         for (int xDeg=0; xDeg <= cubDeg; xDeg++) {
00205           for (int yDeg=0; yDeg <= cubDeg-xDeg; yDeg++) {
00206             for (int zDeg=0; zDeg <= cubDeg-xDeg-yDeg; zDeg++) {
00207               double abstol = ( analyticInt[polyCt+offset][0] == 0.0 ? reltol : std::fabs(reltol*analyticInt[polyCt+offset][0]) );
00208               double absdiff = std::fabs(analyticInt[polyCt+offset][0] - testInt[cubDeg][polyCt]);
00209               if (absdiff > abstol) {
00210                 *outStream << "Cubature order " << std::setw(2) << std::left << cubDeg << " integrating "
00211                            << "x^" << std::setw(2) << std::left << xDeg << " * y^" << std::setw(2) << yDeg
00212                            << " * z^" << std::setw(2) << zDeg << ":" << "   "
00213                            << std::scientific << std::setprecision(16)
00214                            << testInt[cubDeg][polyCt] << "   " << analyticInt[polyCt+offset][0] << "   "
00215                            << std::setprecision(4) << absdiff << "   " << "<?" << "   " << abstol << "\n";
00216                 errorFlag++;
00217                 *outStream << std::right << std::setw(118) << "^^^^---FAILURE!\n";
00218               }
00219               polyCt++;
00220             }
00221             offset = offset + maxOffset - cubDeg;
00222           }
00223           offset = offset + (maxOffset - cubDeg)*(maxOffset - cubDeg + 1)/2;
00224         }
00225         *outStream << "Cubature order " << std::setw(2) << std::left << cubDeg;
00226         if (errorFlag == oldErrorFlag)
00227          *outStream << " passed.\n";
00228         else
00229          *outStream << " failed.\n";
00230       }
00231       *outStream << "\n";
00232   }
00233   catch (std::logic_error err) {
00234     *outStream << err.what() << "\n";
00235     errorFlag = -1;
00236   };
00237 
00238 
00239   if (errorFlag != 0)
00240     std::cout << "End Result: TEST FAILED\n";
00241   else
00242     std::cout << "End Result: TEST PASSED\n";
00243 
00244   // reset format state of std::cout
00245   std::cout.copyfmt(oldFormatState);
00246   return errorFlag;
00247 }