|
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_CubatureDirectLineGauss.hpp" 00052 #include "Intrepid_CubatureDirectLineGaussJacobi20.hpp" 00053 #include "Intrepid_CubatureDirectTriDefault.hpp" 00054 #include "Intrepid_CubatureDirectTetDefault.hpp" 00055 #include "Intrepid_CubatureTensor.hpp" 00056 #include "Intrepid_CubatureTensorPyr.hpp" 00057 #include "Shards_CellTopology.hpp" 00058 #include "Teuchos_oblackholestream.hpp" 00059 #include "Teuchos_RCP.hpp" 00060 #include "Teuchos_GlobalMPISession.hpp" 00061 00062 using namespace Intrepid; 00063 00064 #define INTREPID_TEST_COMMAND( S ) \ 00065 { \ 00066 try { \ 00067 S ; \ 00068 } \ 00069 catch (std::logic_error err) { \ 00070 *outStream << "Expected Error ----------------------------------------------------------------\n"; \ 00071 *outStream << err.what() << '\n'; \ 00072 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \ 00073 }; \ 00074 } 00075 00076 /* 00077 Computes volume of a given reference cell. 00078 */ 00079 double computeRefVolume(shards::CellTopology & cellTopology, int cubDegree) { 00080 Teuchos::RCP< Cubature<double> > myCub; 00081 double vol = 0.0; 00082 00083 switch (cellTopology.getBaseCellTopologyData()->key) { 00084 00085 case shards::Line<>::key: 00086 myCub = Teuchos::rcp(new CubatureDirectLineGauss<double>(cubDegree)); 00087 break; 00088 case shards::Triangle<>::key: 00089 myCub = Teuchos::rcp(new CubatureDirectTriDefault<double>(cubDegree)); 00090 break; 00091 case shards::Tetrahedron<>::key: 00092 myCub = Teuchos::rcp(new CubatureDirectTetDefault<double>(cubDegree)); 00093 break; 00094 case shards::Quadrilateral<>::key: { 00095 std::vector< Teuchos::RCP< Cubature<double> > > lineCubs(2); 00096 lineCubs[0] = Teuchos::rcp(new CubatureDirectLineGauss<double>(cubDegree)); 00097 lineCubs[1] = Teuchos::rcp(new CubatureDirectLineGauss<double>((cubDegree+7) % INTREPID_CUBATURE_LINE_GAUSS_MAX)); 00098 myCub = Teuchos::rcp(new CubatureTensor<double>(lineCubs)); 00099 } 00100 break; 00101 case shards::Hexahedron<>::key: { 00102 std::vector< Teuchos::RCP< Cubature<double> > > lineCubs(2); 00103 std::vector< Teuchos::RCP< Cubature<double> > > miscCubs(2); 00104 lineCubs[0] = Teuchos::rcp(new CubatureDirectLineGauss<double>(cubDegree)); 00105 lineCubs[1] = Teuchos::rcp(new CubatureDirectLineGauss<double>((cubDegree+7) % INTREPID_CUBATURE_LINE_GAUSS_MAX)); 00106 miscCubs[0] = Teuchos::rcp(new CubatureTensor<double>(lineCubs)); 00107 miscCubs[1] = Teuchos::rcp(new CubatureDirectLineGauss<double>((cubDegree+25) % INTREPID_CUBATURE_LINE_GAUSS_MAX)); 00108 myCub = Teuchos::rcp(new CubatureTensor<double>(miscCubs)); 00109 } 00110 break; 00111 case shards::Wedge<>::key: { 00112 Teuchos::RCP<CubatureDirect<double> > triCub = Teuchos::rcp(new CubatureDirectTriDefault<double>(cubDegree)); 00113 Teuchos::RCP<CubatureDirect<double> > lineCub = Teuchos::rcp(new CubatureDirectLineGauss<double>(cubDegree)); 00114 myCub = Teuchos::rcp(new CubatureTensor<double>(triCub,lineCub)); 00115 } 00116 break; 00117 case shards::Pyramid<>::key: { 00118 std::vector< Teuchos::RCP< Cubature<double> > > lineCubs(3); 00119 lineCubs[0] = Teuchos::rcp(new CubatureDirectLineGauss<double>(cubDegree)); 00120 lineCubs[1] = Teuchos::rcp(new CubatureDirectLineGauss<double>(cubDegree)); 00121 lineCubs[2] = Teuchos::rcp(new CubatureDirectLineGaussJacobi20<double>(cubDegree)); 00122 myCub = Teuchos::rcp(new CubatureTensorPyr<double>(lineCubs)); 00123 } 00124 break; 00125 00126 default: 00127 TEUCHOS_TEST_FOR_EXCEPTION( ( (cellTopology.getBaseCellTopologyData()->key != shards::Line<>::key), 00128 (cellTopology.getBaseCellTopologyData()->key != shards::Triangle<>::key), 00129 (cellTopology.getBaseCellTopologyData()->key != shards::Tetrahedron<>::key), 00130 (cellTopology.getBaseCellTopologyData()->key != shards::Quadrilateral<>::key), 00131 (cellTopology.getBaseCellTopologyData()->key != shards::Hexahedron<>::key), 00132 (cellTopology.getBaseCellTopologyData()->key != shards::Wedge<>::key), 00133 (cellTopology.getBaseCellTopologyData()->key != shards::Pyramid<>::key) ), 00134 std::invalid_argument, 00135 ">>> ERROR (Unit Test -- Cubature -- Volume): Invalid cell type."); 00136 } // end switch 00137 00138 int numCubPoints = myCub->getNumPoints(); 00139 00140 FieldContainer<double> cubPoints( numCubPoints, myCub->getDimension() ); 00141 FieldContainer<double> cubWeights( numCubPoints ); 00142 myCub->getCubature(cubPoints, cubWeights); 00143 00144 for (int i=0; i<numCubPoints; i++) 00145 vol += cubWeights[i]; 00146 00147 return vol; 00148 } 00149 00150 00151 int main(int argc, char *argv[]) { 00152 00153 Teuchos::GlobalMPISession mpiSession(&argc, &argv); 00154 00155 // This little trick lets us print to std::cout only if 00156 // a (dummy) command-line argument is provided. 00157 int iprint = argc - 1; 00158 Teuchos::RCP<std::ostream> outStream; 00159 Teuchos::oblackholestream bhs; // outputs nothing 00160 if (iprint > 0) 00161 outStream = Teuchos::rcp(&std::cout, false); 00162 else 00163 outStream = Teuchos::rcp(&bhs, false); 00164 00165 // Save the format state of the original std::cout. 00166 Teuchos::oblackholestream oldFormatState; 00167 oldFormatState.copyfmt(std::cout); 00168 00169 *outStream \ 00170 << "===============================================================================\n" \ 00171 << "| |\n" \ 00172 << "| Unit Test (CubatureDirect,CubatureTensor) |\n" \ 00173 << "| |\n" \ 00174 << "| 1) Computing volumes of reference cells |\n" \ 00175 << "| |\n" \ 00176 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \ 00177 << "| Denis Ridzal (dridzal@sandia.gov). |\n" \ 00178 << "| |\n" \ 00179 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \ 00180 << "| Trilinos website: http://trilinos.sandia.gov |\n" \ 00181 << "| |\n" \ 00182 << "===============================================================================\n"\ 00183 << "| TEST 1: construction and basic functionality |\n"\ 00184 << "===============================================================================\n"; 00185 00186 int errorFlag = 0; 00187 00188 int beginThrowNumber = Teuchos::TestForException_getThrowNumber(); 00189 int endThrowNumber = beginThrowNumber + 7; 00190 00191 try { 00192 /* Line cubature. */ 00193 INTREPID_TEST_COMMAND( CubatureDirectLineGauss<double> lineCub(-1) ); 00194 INTREPID_TEST_COMMAND( CubatureDirectLineGauss<double> lineCub(INTREPID_CUBATURE_LINE_GAUSS_MAX+1) ); 00195 INTREPID_TEST_COMMAND( CubatureDirectLineGauss<double> lineCub; 00196 std::string testName = "INTREPID_CUBATURE_LINE_GAUSS"; 00197 std::string lineCubName = lineCub.getName(); 00198 *outStream << "\nComparing strings: " << testName << " and " << lineCubName << "\n\n"; 00199 TEUCHOS_TEST_FOR_EXCEPTION( (testName != lineCubName), std::logic_error, "Name mismatch!" ) ); 00200 INTREPID_TEST_COMMAND( CubatureDirectLineGauss<double> lineCub; 00201 std::vector<int> accuracy; 00202 lineCub.getAccuracy(accuracy); 00203 TEUCHOS_TEST_FOR_EXCEPTION( (accuracy[0] != 0), std::logic_error, "Check member getAccuracy!" ) ); 00204 INTREPID_TEST_COMMAND( CubatureDirectLineGauss<double> lineCub(55); 00205 TEUCHOS_TEST_FOR_EXCEPTION( (lineCub.getNumPoints() != 28), std::logic_error, "Check member getNumPoints!" ) ); 00206 INTREPID_TEST_COMMAND( CubatureDirectLineGauss<double> lineCub; 00207 TEUCHOS_TEST_FOR_EXCEPTION( (lineCub.getDimension() != 1), 00208 std::logic_error, 00209 "Check member dimension!" ) ); 00210 /* Triangle cubature. */ 00211 INTREPID_TEST_COMMAND( CubatureDirectTriDefault<double> triCub(-1) ); 00212 INTREPID_TEST_COMMAND( CubatureDirectTriDefault<double> triCub(INTREPID_CUBATURE_TRI_DEFAULT_MAX+1) ); 00213 INTREPID_TEST_COMMAND( CubatureDirectTriDefault<double> triCub; 00214 std::string testName = "INTREPID_CUBATURE_TRI_DEFAULT"; 00215 std::string triCubName = triCub.getName(); 00216 *outStream << "\nComparing strings: " << testName << " and " << triCubName << "\n\n"; 00217 TEUCHOS_TEST_FOR_EXCEPTION( (testName != triCubName), std::logic_error, "Name mismatch!" ) ); 00218 INTREPID_TEST_COMMAND( CubatureDirectTriDefault<double> triCub; 00219 std::vector<int> accuracy; 00220 triCub.getAccuracy(accuracy); 00221 TEUCHOS_TEST_FOR_EXCEPTION( (accuracy[0] != 0), std::logic_error, "Check member getAccuracy!" ) ); 00222 INTREPID_TEST_COMMAND( CubatureDirectTriDefault<double> triCub(17); 00223 TEUCHOS_TEST_FOR_EXCEPTION( (triCub.getNumPoints() != 61), std::logic_error, "Check member getNumPoints!" ) ); 00224 INTREPID_TEST_COMMAND( CubatureDirectTriDefault<double> triCub; 00225 TEUCHOS_TEST_FOR_EXCEPTION( (triCub.getDimension() != 2), 00226 std::logic_error, 00227 "Check member dimension!" ) ); 00228 /* Tetrahedron cubature. */ 00229 INTREPID_TEST_COMMAND( CubatureDirectTetDefault<double> tetCub(-1) ); 00230 INTREPID_TEST_COMMAND( CubatureDirectTetDefault<double> tetCub(INTREPID_CUBATURE_TET_DEFAULT_MAX+1) ); 00231 INTREPID_TEST_COMMAND( CubatureDirectTetDefault<double> tetCub; 00232 std::string testName = "INTREPID_CUBATURE_TET_DEFAULT"; 00233 std::string tetCubName = tetCub.getName(); 00234 *outStream << "\nComparing strings: " << testName << " and " << tetCubName << "\n\n"; 00235 std::vector< Teuchos::RCP< Cubature<double> > > lineCubs(2); 00236 TEUCHOS_TEST_FOR_EXCEPTION( (testName != tetCubName), std::logic_error, "Name mismatch!" ) ); 00237 INTREPID_TEST_COMMAND( CubatureDirectTetDefault<double> tetCub; 00238 std::vector<int> accuracy; 00239 tetCub.getAccuracy(accuracy); 00240 TEUCHOS_TEST_FOR_EXCEPTION( (accuracy[0] != 0), std::logic_error, "Check member getAccuracy!" ) ); 00241 INTREPID_TEST_COMMAND( CubatureDirectTetDefault<double> tetCub(17); 00242 TEUCHOS_TEST_FOR_EXCEPTION( (tetCub.getNumPoints() != 495), std::logic_error, "Check member getNumPoints!" ) ); 00243 INTREPID_TEST_COMMAND( CubatureDirectTetDefault<double> tetCub; 00244 TEUCHOS_TEST_FOR_EXCEPTION( (tetCub.getDimension() != 3), 00245 std::logic_error, 00246 "Check member getCellTopology!" ) ); 00247 /* Tensor cubature. */ 00248 INTREPID_TEST_COMMAND( std::vector< Teuchos::RCP< Cubature<double> > > lineCubs(0); 00249 CubatureTensor<double> quadCub(lineCubs) ); 00250 INTREPID_TEST_COMMAND( std::vector< Teuchos::RCP< Cubature<double> > > miscCubs(3); 00251 std::vector< Teuchos::RCP< Cubature<double> > > lineCubs(2); 00252 lineCubs[0] = Teuchos::rcp(new CubatureDirectLineGauss<double>(3)); 00253 lineCubs[1] = Teuchos::rcp(new CubatureDirectLineGauss<double>(16)); 00254 miscCubs[0] = Teuchos::rcp(new CubatureTensor<double>(lineCubs)); 00255 miscCubs[1] = Teuchos::rcp(new CubatureDirectTriDefault<double>); 00256 miscCubs[2] = Teuchos::rcp(new CubatureDirectTetDefault<double>(19)); 00257 CubatureTensor<double> tensorCub(miscCubs); 00258 std::vector<int> a(4); a[0]=3; a[1]=16; a[2]=0; a[3]=19; 00259 std::vector<int> atest(4); 00260 tensorCub.getAccuracy(atest); 00261 TEUCHOS_TEST_FOR_EXCEPTION( (a != atest), std::logic_error, "Check member getAccuracy!" ) ); 00262 INTREPID_TEST_COMMAND( std::vector< Teuchos::RCP< Cubature<double> > > lineCubs(2); 00263 lineCubs[0] = Teuchos::rcp(new CubatureDirectLineGauss<double>(15)); 00264 lineCubs[1] = Teuchos::rcp(new CubatureDirectLineGauss<double>(11)); 00265 CubatureTensor<double> tensorCub(lineCubs); 00266 TEUCHOS_TEST_FOR_EXCEPTION( (tensorCub.getNumPoints() != 48), std::logic_error, "Check member getNumPoints!" ) ); 00267 INTREPID_TEST_COMMAND( std::vector< Teuchos::RCP< Cubature<double> > > miscCubs(3); 00268 miscCubs[0] = Teuchos::rcp(new CubatureDirectLineGauss<double>); 00269 miscCubs[1] = Teuchos::rcp(new CubatureDirectTriDefault<double>); 00270 miscCubs[2] = Teuchos::rcp(new CubatureDirectTetDefault<double>); 00271 CubatureTensor<double> tensorCub(miscCubs); 00272 TEUCHOS_TEST_FOR_EXCEPTION( (tensorCub.getDimension() != 6), std::logic_error, "Check member dimension!" ) ); 00273 INTREPID_TEST_COMMAND( std::vector< Teuchos::RCP< Cubature<double> > > miscCubs(3); 00274 miscCubs[0] = Teuchos::rcp(new CubatureDirectLineGauss<double>(3)); 00275 miscCubs[1] = Teuchos::rcp(new CubatureDirectLineGauss<double>(7)); 00276 miscCubs[2] = Teuchos::rcp(new CubatureDirectLineGauss<double>(5)); 00277 CubatureTensor<double> tensorCub(miscCubs); 00278 FieldContainer<double> points(tensorCub.getNumPoints(), tensorCub.getDimension()); 00279 FieldContainer<double> weights(tensorCub.getNumPoints()); 00280 tensorCub.getCubature(points, weights) 00281 ) 00282 INTREPID_TEST_COMMAND( Teuchos::RCP< CubatureDirect<double> > lineCub = Teuchos::rcp(new CubatureDirectLineGauss<double>(15)); 00283 Teuchos::RCP< CubatureDirect<double> > triCub = Teuchos::rcp(new CubatureDirectTriDefault<double>(12)); 00284 CubatureTensor<double> tensorCub(lineCub, triCub); 00285 std::vector<int> a(2); a[0] = 15; a[1] = 12; 00286 std::vector<int> atest(2); 00287 tensorCub.getAccuracy(atest); 00288 TEUCHOS_TEST_FOR_EXCEPTION( (tensorCub.getDimension() != 3) || (a != atest), 00289 std::logic_error, 00290 "Check constructormembers dimension and getAccuracy!" ) ); 00291 INTREPID_TEST_COMMAND( Teuchos::RCP< CubatureDirect<double> > lineCub = Teuchos::rcp(new CubatureDirectLineGauss<double>(15)); 00292 Teuchos::RCP< CubatureDirect<double> > triCub = Teuchos::rcp(new CubatureDirectTriDefault<double>(12)); 00293 CubatureTensor<double> tensorCub(triCub, lineCub, triCub); 00294 std::vector<int> a(3); a[0] = 12; a[1] = 15; a[2] = 12; 00295 std::vector<int> atest(3); 00296 tensorCub.getAccuracy(atest); 00297 TEUCHOS_TEST_FOR_EXCEPTION( (tensorCub.getDimension() != 5) || (a != atest), 00298 std::logic_error, 00299 "Check constructor and members dimension and getAccuracy!" ) ); 00300 INTREPID_TEST_COMMAND( Teuchos::RCP< CubatureDirect<double> > triCub = Teuchos::rcp(new CubatureDirectTriDefault<double>(12)); 00301 CubatureTensor<double> tensorCub(triCub, 5); 00302 std::vector<int> a(5); a[0] = 12; a[1] = 12; a[2] = 12; a[3] = 12; a[4] = 12; 00303 std::vector<int> atest(5); 00304 tensorCub.getAccuracy(atest); 00305 TEUCHOS_TEST_FOR_EXCEPTION( (tensorCub.getDimension() != 10) || (a != atest), 00306 std::logic_error, 00307 "Check constructor and members dimension and getAccuracy!" ) ); 00308 if (Teuchos::TestForException_getThrowNumber() != endThrowNumber) { 00309 errorFlag = -1000; 00310 } 00311 } 00312 catch (std::logic_error err) { 00313 *outStream << err.what() << "\n"; 00314 errorFlag = -1000; 00315 }; 00316 00317 *outStream \ 00318 << "===============================================================================\n"\ 00319 << "| TEST 2: volume computations |\n"\ 00320 << "===============================================================================\n"; 00321 00322 double testVol = 0.0; 00323 double tol = 100.0 * INTREPID_TOL; 00324 00325 // list of analytic volume values, listed in the enumerated reference cell order up to CELL_HEXAPRISM 00326 double volumeList[] = {0.0, 2.0, 1.0/2.0, 4.0, 1.0/6.0, 8.0, 1.0, 4.0/3.0, 32.0}; 00327 00328 *outStream << "\nReference cell volumes:\n\n"; 00329 00330 try { 00331 shards::CellTopology line(shards::getCellTopologyData< shards::Line<> >()); 00332 for (int deg=0; deg<=INTREPID_CUBATURE_LINE_GAUSS_MAX; deg++) { 00333 testVol = computeRefVolume(line, deg); 00334 *outStream << std::setw(30) << "Line volume --> " << std::setw(10) << std::scientific << testVol << 00335 std::setw(10) << "diff = " << std::setw(10) << std::scientific << std::abs(testVol - volumeList[1]) << "\n"; 00336 if (std::abs(testVol - volumeList[1]) > tol) { 00337 errorFlag = 1; 00338 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00339 } 00340 } 00341 *outStream << "\n\n"; 00342 shards::CellTopology tri(shards::getCellTopologyData< shards::Triangle<> >()); 00343 for (int deg=0; deg<=INTREPID_CUBATURE_TRI_DEFAULT_MAX; deg++) { 00344 testVol = computeRefVolume(tri, deg); 00345 *outStream << std::setw(30) << "Triangle volume --> " << std::setw(10) << std::scientific << testVol << 00346 std::setw(10) << "diff = " << std::setw(10) << std::scientific << std::abs(testVol - volumeList[2]) << "\n"; 00347 if (std::abs(testVol - volumeList[2]) > tol) { 00348 errorFlag = 1; 00349 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00350 } 00351 } 00352 *outStream << "\n\n"; 00353 shards::CellTopology quad(shards::getCellTopologyData< shards::Quadrilateral<> >()); 00354 for (int deg=0; deg<=INTREPID_CUBATURE_LINE_GAUSS_MAX; deg++) { 00355 testVol = computeRefVolume(quad, deg); 00356 *outStream << std::setw(30) << "Quadrilateral volume --> " << std::setw(10) << std::scientific << testVol << 00357 std::setw(10) << "diff = " << std::setw(10) << std::scientific << std::abs(testVol - volumeList[3]) << "\n"; 00358 if (std::abs(testVol - volumeList[3]) > tol) { 00359 errorFlag = 1; 00360 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00361 } 00362 } 00363 *outStream << "\n\n"; 00364 shards::CellTopology tet(shards::getCellTopologyData< shards::Tetrahedron<> >()); 00365 for (int deg=0; deg<=INTREPID_CUBATURE_TET_DEFAULT_MAX; deg++) { 00366 testVol = computeRefVolume(tet, deg); 00367 *outStream << std::setw(30) << "Tetrahedron volume --> " << std::setw(10) << std::scientific << testVol << 00368 std::setw(10) << "diff = " << std::setw(10) << std::scientific << std::abs(testVol - volumeList[4]) << "\n"; 00369 if (std::abs(testVol - volumeList[4]) > tol) { 00370 errorFlag = 1; 00371 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00372 } 00373 } 00374 *outStream << "\n\n"; 00375 shards::CellTopology hex(shards::getCellTopologyData< shards::Hexahedron<> >()); 00376 for (int deg=0; deg<=INTREPID_CUBATURE_LINE_GAUSS_MAX; deg++) { 00377 testVol = computeRefVolume(hex, deg); 00378 *outStream << std::setw(30) << "Hexahedron volume --> " << std::setw(10) << std::scientific << testVol << 00379 std::setw(10) << "diff = " << std::setw(10) << std::scientific << std::abs(testVol - volumeList[5]) << "\n"; 00380 if (std::abs(testVol - volumeList[5]) > tol) { 00381 errorFlag = 1; 00382 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00383 } 00384 } 00385 *outStream << "\n\n"; 00386 shards::CellTopology wedge(shards::getCellTopologyData< shards::Wedge<> >()); 00387 for (int deg=0; deg<=std::min(INTREPID_CUBATURE_LINE_GAUSS_MAX,INTREPID_CUBATURE_TRI_DEFAULT_MAX); deg++) { 00388 testVol = computeRefVolume(wedge, deg); 00389 *outStream << std::setw(30) << "Wedge volume --> " << std::setw(10) << std::scientific << testVol << 00390 std::setw(10) << "diff = " << std::setw(10) << std::scientific << std::abs(testVol - volumeList[6]) << "\n"; 00391 if (std::abs(testVol - volumeList[6]) > tol) { 00392 errorFlag = 1; 00393 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00394 } 00395 } 00396 *outStream << "\n\n"; 00397 shards::CellTopology pyr(shards::getCellTopologyData< shards::Pyramid<> >()); 00398 for (int deg=0; deg<=std::min(INTREPID_CUBATURE_LINE_GAUSS_MAX,INTREPID_CUBATURE_LINE_GAUSSJACOBI20_MAX); deg++) { 00399 testVol = computeRefVolume(pyr, deg); 00400 *outStream << std::setw(30) << "Pyramid volume --> " << std::setw(10) << std::scientific << testVol << 00401 std::setw(10) << "diff = " << std::setw(10) << std::scientific << std::abs(testVol - volumeList[7]) << "\n"; 00402 if (std::abs(testVol - volumeList[7]) > tol) { 00403 errorFlag = 1; 00404 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00405 } 00406 } 00407 *outStream << "\n\n"; 00408 for (int deg=0; deg<=20; deg++) { 00409 Teuchos::RCP<CubatureDirectLineGauss<double> > lineCub = Teuchos::rcp(new CubatureDirectLineGauss<double>(deg)); 00410 CubatureTensor<double> hypercubeCub(lineCub, 5); 00411 int numCubPoints = hypercubeCub.getNumPoints(); 00412 FieldContainer<double> cubPoints( numCubPoints, hypercubeCub.getDimension() ); 00413 FieldContainer<double> cubWeights( numCubPoints ); 00414 hypercubeCub.getCubature(cubPoints, cubWeights); 00415 testVol = 0; 00416 for (int i=0; i<numCubPoints; i++) 00417 testVol += cubWeights[i]; 00418 *outStream << std::setw(30) << "5-D Hypercube volume --> " << std::setw(10) << std::scientific << testVol << 00419 std::setw(10) << "diff = " << std::setw(10) << std::scientific << std::abs(testVol - volumeList[8]) << "\n"; 00420 if (std::abs(testVol - volumeList[8])/std::abs(testVol) > tol) { 00421 errorFlag = 1; 00422 *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n"; 00423 } 00424 } 00425 } 00426 catch (std::logic_error err) { 00427 *outStream << err.what() << "\n"; 00428 errorFlag = -1; 00429 }; 00430 00431 00432 if (errorFlag != 0) 00433 std::cout << "End Result: TEST FAILED\n"; 00434 else 00435 std::cout << "End Result: TEST PASSED\n"; 00436 00437 // reset format state of std::cout 00438 std::cout.copyfmt(oldFormatState); 00439 00440 return errorFlag; 00441 }
1.7.6.1