|
Intrepid
|
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_CubatureSparse.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 CubatureSparse<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 (CubatureSparse) |\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 Sparse |\n"\ 00155 << "| 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 }
1.7.6.1