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