Intrepid
http://trilinos.sandia.gov/packages/docs/r10.10/packages/intrepid/test/Shared/PointTools/test_01.cpp
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 
00049 #include "Intrepid_PointTools.hpp"
00050 #include "Intrepid_FieldContainer.hpp"
00051 #include "Teuchos_oblackholestream.hpp"
00052 #include "Teuchos_GlobalMPISession.hpp"
00053 #include "Shards_CellTopology.hpp"
00054 
00055 using namespace std;
00056 using namespace Intrepid;
00057 
00058 
00059 #define INTREPID_TEST_COMMAND( S , throwCounter, nException )                                                              \
00060 {                                                                                                                          \
00061   ++nException;                                                                                                            \
00062   try {                                                                                                                    \
00063     S ;                                                                                                                    \
00064   }                                                                                                                        \
00065   catch (std::logic_error err) {                                                                                           \
00066       ++throwCounter;                                                                                                      \
00067       *outStream << "Expected Error " << nException << " -------------------------------------------------------------\n"; \
00068       *outStream << err.what() << '\n';                                                                                    \
00069       *outStream << "-------------------------------------------------------------------------------" << "\n\n";           \
00070   };                                                                                                                       \
00071 }
00072 
00073 
00074 int main(int argc, char *argv[]) {
00075 
00076   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00077 
00078   // This little trick lets us print to std::cout only if
00079   // a (dummy) command-line argument is provided.
00080   int iprint     = argc - 1;
00081   Teuchos::RCP<std::ostream> outStream;
00082   Teuchos::oblackholestream bhs; // outputs nothing
00083   if (iprint > 0)
00084     outStream = Teuchos::rcp(&std::cout, false);
00085   else
00086     outStream = Teuchos::rcp(&bhs, false);
00087 
00088   // Save the format state of the original std::cout.
00089   Teuchos::oblackholestream oldFormatState;
00090   oldFormatState.copyfmt(std::cout);
00091 
00092   *outStream \
00093   << "===============================================================================\n" \
00094   << "|                                                                             |\n" \
00095   << "|                       Unit Test (PointTools)                                |\n" \
00096   << "|                                                                             |\n" \
00097   << "|     1) Construction of equispaced and warped lattices on simplices          |\n" \
00098   << "|                                                                             |\n" \
00099   << "|  Questions? Contact  Pavel Bochev (pbboche@sandia.gov) or                   |\n" \
00100   << "|                      Denis Ridzal (dridzal@sandia.gov) or                   |\n" \
00101   << "|                      Robert Kirby (robert.c.kirby@ttu.edu)                  |\n" \
00102   << "|                                                                             |\n" \
00103   << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00104   << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00105   << "|                                                                             |\n" \
00106   << "===============================================================================\n";
00107 
00108 
00109 
00110   int errorFlag = 0;
00111 
00112   shards::CellTopology myLine_2( shards::getCellTopologyData< shards::Line<2> >() );
00113   shards::CellTopology myTri_3( shards::getCellTopologyData< shards::Triangle<3> >() );
00114   shards::CellTopology myTet_4( shards::getCellTopologyData< shards::Tetrahedron<4> >() );
00115 
00116 
00117   *outStream \
00118   << "\n"
00119   << "===============================================================================\n"\
00120   << "| TEST 1: size of lattices                                                   |\n"\
00121   << "===============================================================================\n";
00122 
00123   try{
00124     // first try the lattices with offset = 0.  This is a spot-check
00125 
00126     if (PointTools::getLatticeSize( myLine_2 , 4 , 0 ) != 5) {
00127       errorFlag++;
00128       *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00129       *outStream << " size of 4th order lattice on a line with no offset: " << PointTools::getLatticeSize( myLine_2 , 4 , 0 )  << "\n";
00130       *outStream << " should be 5\n";
00131     }
00132 
00133 
00134     if (PointTools::getLatticeSize( myTri_3 , 3 , 0 ) != 10) {
00135       errorFlag++;
00136       *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00137       *outStream << " size of 3rd order lattice on a line with no offset: " << PointTools::getLatticeSize( myTri_3 , 3 , 0 )  << "\n";
00138       *outStream << " should be 10\n";    
00139     }
00140 
00141 
00142     if (PointTools::getLatticeSize( myTet_4 , 3 , 0 ) != 20) {
00143       errorFlag++;
00144       *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00145       *outStream << " size of 3rd order lattice on a tet with no offset: " << PointTools::getLatticeSize( myTet_4 , 3 , 0 )  << "\n";
00146       *outStream << " should be 20\n"; 
00147     }
00148 
00149                         
00150     // check with the offset = 1
00151     if (PointTools::getLatticeSize( myLine_2 , 3 , 1 ) != 2) {
00152       errorFlag++;
00153       *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00154       *outStream << " size of 3rd order lattice on a line with 1 offset: " << PointTools::getLatticeSize( myLine_2 , 3 , 1 )  << "\n";
00155       *outStream << " should be 2\n";       
00156     }
00157 
00158     if (PointTools::getLatticeSize( myTri_3 , 4 , 1 ) != 3) {
00159       errorFlag++;
00160       *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00161       *outStream << " size of 4th order lattice on a triangle with 1 offset: " << PointTools::getLatticeSize( myTri_3 , 4 , 1 )  << "\n";
00162       *outStream << " should be 3\n";           
00163     }
00164 
00165     if (PointTools::getLatticeSize( myTet_4 , 5 , 1 ) != 4) {
00166       errorFlag++;
00167       *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00168       *outStream << " size of 5th order lattice on a tetrahedron with 1 offset: " << PointTools::getLatticeSize( myTet_4 , 5 , 1 )  << "\n";
00169       *outStream << " should be 4\n";   
00170     }
00171 
00172   }
00173   catch (std::exception &err) {
00174     *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
00175     *outStream << err.what() << '\n';
00176     *outStream << "-------------------------------------------------------------------------------" << "\n\n";
00177     errorFlag = -1000;
00178   };
00179 
00180   // Now verify that we throw an exception on some of the non-supported cell types.
00181 
00182   *outStream \
00183   << "\n"
00184   << "===============================================================================\n"\
00185   << "| TEST 2: check for unsupported cell types                                     \n"\
00186   << "===============================================================================\n";
00187   try{
00188     try {
00189       PointTools::getLatticeSize( shards::getCellTopologyData< shards::Quadrilateral<4> >() , 3 , 0 );
00190     }
00191     catch (std::invalid_argument err) {
00192       *outStream << err.what() << "\n";
00193     }
00194 
00195     try {
00196       PointTools::getLatticeSize( shards::getCellTopologyData< shards::Hexahedron<8> >() , 3 , 0 );
00197     }
00198     catch (std::invalid_argument err) {
00199       *outStream << err.what() << "\n";
00200     }
00201   }
00202   catch (std::exception &err) {
00203     *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
00204     *outStream << err.what() << '\n';
00205     *outStream << "-------------------------------------------------------------------------------" << "\n\n";
00206     errorFlag = -1000;
00207   };
00208 
00209   // Check for malformed point arrays
00210 
00211 
00212   *outStream \
00213   << "\n"
00214   << "===============================================================================\n"\
00215   << "| TEST 2: malformed point arrays                                               \n"\
00216   << "===============================================================================\n";
00217   try{
00218     // line: not enough points allocated
00219     try {
00220       FieldContainer<double> pts(3,1);
00221       PointTools::getLatticeSize( shards::getCellTopologyData< shards::Line<2> >() , 5 , 0 );
00222     }
00223     catch (std::invalid_argument err) {
00224       *outStream << err.what() << "\n";
00225     }
00226     // line: wrong dimension for points
00227     try {
00228       FieldContainer<double> pts(6,2);
00229       PointTools::getLatticeSize( shards::getCellTopologyData< shards::Line<2> >() , 5 , 0 );
00230     }
00231     catch (std::invalid_argument err) {
00232       *outStream << err.what() << "\n";
00233     }
00234     // triangle: too many points allocated
00235     try {
00236       FieldContainer<double> pts(4,2);
00237       PointTools::getLatticeSize( shards::getCellTopologyData< shards::Triangle<3> >() , 3 , 1 );
00238     }
00239     catch (std::invalid_argument err) {
00240       *outStream << err.what() << "\n";
00241     }
00242     // triangle: wrong dimension for points
00243     try {
00244       FieldContainer<double> pts(6,1);
00245       PointTools::getLatticeSize( shards::getCellTopologyData< shards::Triangle<3> >() , 3 , 0 );
00246     }
00247     catch (std::invalid_argument err) {
00248       *outStream << err.what() << "\n";
00249     }
00250     // tetrahedron: not enough points allocated
00251     try {
00252       FieldContainer<double> pts(4,2);
00253       PointTools::getLatticeSize( shards::getCellTopologyData< shards::Tetrahedron<4> >() , 2 , 0 );
00254     }
00255     catch (std::invalid_argument err) {
00256       *outStream << err.what() << "\n";
00257     }
00258     // tetrahedron: wrong dimension for points
00259     try {
00260       FieldContainer<double> pts(4,2);
00261       PointTools::getLatticeSize( shards::getCellTopologyData< shards::Tetrahedron<4> >() , 1 , 0 );
00262     }
00263     catch (std::invalid_argument err) {
00264       *outStream << err.what() << "\n";
00265     }
00266 
00267   }
00268   catch (std::exception &err) {
00269     *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
00270     *outStream << err.what() << '\n';
00271     *outStream << "-------------------------------------------------------------------------------" << "\n\n";
00272     errorFlag = -1000;
00273   };
00274 
00275 
00276   *outStream \
00277   << "\n"
00278   << "===============================================================================\n"\
00279   << "| TEST 3: check values of triangular lattice compared to Warburton's code      \n"\
00280   << "===============================================================================\n";
00281   try {
00282     // triangle case
00283     const int order = 4;
00284     const int offset = 0;
00285     int numPts = PointTools::getLatticeSize( myTri_3 , order , offset );
00286     int ptDim = 2;
00287     FieldContainer<double> warpBlendPts( numPts , ptDim );
00288     PointTools::getLattice<double,FieldContainer<double> >( warpBlendPts ,
00289                                                             myTri_3 ,
00290                                                             order ,
00291                                                             offset ,
00292                                                             POINTTYPE_WARPBLEND );
00293     FieldContainer<double> verts( 1, 3 , 2 );
00294     verts(0,0,0) = -1.0;
00295     verts(0,0,1) = -1.0/sqrt(3.0);
00296     verts(0,1,0) = 1.0;
00297     verts(0,1,1) = -1.0/sqrt(3.0);
00298     verts(0,2,0) = 0.0;
00299     verts(0,2,1) = 2.0/sqrt(3.0);
00300 
00301     // holds points on the equilateral triangle
00302     FieldContainer<double> warpBlendMappedPts( numPts , ptDim );
00303 
00304     CellTools<double>::mapToPhysicalFrame(warpBlendMappedPts ,
00305                                           warpBlendPts ,
00306                                           verts ,
00307                                           myTri_3 ,
00308                                           0 );
00309 
00310     // Values from MATLAB code
00311     double points[] = { -1.000000000000000 , -0.577350269189626 ,
00312                         -0.654653670707977 , -0.577350269189626 ,
00313                         -0.000000000000000 , -0.577350269189626 ,
00314                         0.654653670707977  , -0.577350269189626 ,
00315                         1.000000000000000  , -0.577350269189626 ,
00316                         -0.827326835353989 , -0.278271574919028 ,
00317                         -0.327375261332958 , -0.189010195256608 ,
00318                         0.327375261332958 , -0.189010195256608 ,
00319                         0.827326835353989,  -0.278271574919028,
00320                         -0.500000000000000,   0.288675134594813,
00321                         0.000000000000000,   0.378020390513215,
00322                         0.500000000000000,   0.288675134594813,
00323                         -0.172673164646011,   0.855621844108654,
00324                         0.172673164646011,   0.855621844108654,
00325                         0,   1.154700538379252 };
00326 
00327     // compare
00328     for (int i=0;i<numPts;i++) {
00329       for (int j=0;j<2;j++) {
00330         int l = 2*i+j;
00331         if (std::abs(warpBlendMappedPts(i,j) - points[l]) > INTREPID_TOL) {
00332           errorFlag++;
00333           *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00334           
00335           // Output the multi-index of the value where the error is:
00336           *outStream << " At multi-index { ";
00337           *outStream << i << " ";*outStream << j << " ";
00338           *outStream << "}  computed value: " << warpBlendMappedPts(i,j)
00339                      << " but correct value: " << points[l] << "\n";
00340         }
00341       }
00342     }
00343 
00344 
00345   }
00346   catch (std::exception &err) {
00347     *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
00348     *outStream << err.what() << '\n';
00349     *outStream << "-------------------------------------------------------------------------------" << "\n\n";
00350     errorFlag = -1000;
00351   };
00352 
00353 
00354   *outStream \
00355   << "\n"
00356   << "===============================================================================\n"\
00357   << "| TEST 4: check values of tetrahedral lattice compared to Warburton's code     \n"\
00358   << "===============================================================================\n";
00359   try {
00360     // triangle case
00361     const int order = 6;
00362     const int offset = 0;
00363     int numPts = PointTools::getLatticeSize( myTet_4 , order , offset );
00364     int ptDim = 3;
00365     FieldContainer<double> warpBlendPts( numPts , ptDim );
00366     PointTools::getLattice<double,FieldContainer<double> >( warpBlendPts ,
00367                                                             myTet_4 ,
00368                                                             order ,
00369                                                             offset ,
00370                                                             POINTTYPE_WARPBLEND );
00371 
00372     FieldContainer<double> verts(1,4,3);
00373     verts(0,0,0) = -1.0;
00374     verts(0,0,1) = -1.0/sqrt(3.0);
00375     verts(0,0,2) = -1.0/sqrt(6.0);
00376     verts(0,1,0) = 1.0;
00377     verts(0,1,1) = -1.0/sqrt(3.0);
00378     verts(0,1,2) = -1.0/sqrt(6.0);
00379     verts(0,2,0) = 0.0;
00380     verts(0,2,1) = 2.0 / sqrt(3.0);
00381     verts(0,2,2) = -1.0/sqrt(6.0);
00382     verts(0,3,0) = 0.0;
00383     verts(0,3,1) = 0.0;
00384     verts(0,3,2) = 3.0 / sqrt(6.0);
00385 
00386 
00387     // points on the equilateral tet
00388     FieldContainer<double> warpBlendMappedPts( numPts , ptDim );
00389 
00390     CellTools<double>::mapToPhysicalFrame(warpBlendMappedPts ,
00391                                           warpBlendPts ,
00392                                           verts ,
00393                                           myTet_4 ,
00394                                           0 );
00395 
00396     // Values from MATLAB code
00397     double points[] = {  -1.000000000000000,  -0.577350269189626,  -0.408248290463863,
00398                         -0.830223896278567,  -0.577350269189626,  -0.408248290463863,
00399                         -0.468848793470714,  -0.577350269189626,  -0.408248290463863,
00400                         -0.000000000000000,  -0.577350269189626,  -0.408248290463863,
00401                         0.468848793470714,  -0.577350269189626, -0.408248290463863,
00402                         0.830223896278567,  -0.577350269189626,  -0.408248290463863,
00403                         1.000000000000000,  -0.577350269189626,  -0.408248290463863,
00404                         -0.915111948139283,  -0.430319850411323,  -0.408248290463863,
00405                         -0.660434383303427,  -0.381301968982318,  -0.408248290463863,
00406                         -0.239932664820086,  -0.368405260495326,  -0.408248290463863,
00407                         0.239932664820086,  -0.368405260495326,  -0.408248290463863,
00408                         0.660434383303426, -0.381301968982318,  -0.408248290463863,
00409                         0.915111948139283,  -0.430319850411323,  -0.408248290463863,
00410                         -0.734424396735357,  -0.117359831084509,  -0.408248290463863,
00411                         -0.439014646886819,  -0.023585152684228,  -0.408248290463863,
00412                         -0.000000000000000,  -0.000000000000000,  -0.408248290463863,
00413                         0.439014646886819,  -0.023585152684228,  -0.408248290463863,
00414                         0.734424396735357,  -0.117359831084509,  -0.408248290463863,
00415                         -0.500000000000000,   0.288675134594813,  -0.408248290463863,
00416                         -0.199081982066733,   0.391990413179555,  -0.408248290463863,
00417                         0.199081982066733,   0.391990413179555,  -0.408248290463863,
00418                         0.500000000000000,   0.288675134594813,  -0.408248290463863,
00419                         -0.265575603264643,   0.694710100274135,  -0.408248290463863,
00420                         -0.000000000000000,  0.762603937964635,  -0.408248290463863,
00421                         0.265575603264643,   0.694710100274135,  -0.408248290463863,
00422                         -0.084888051860716,   1.007670119600949,  -0.408248290463863,
00423                         0.084888051860716,   1.007670119600949,  -0.408248290463863,
00424                         0,   1.154700538379252,  -0.408248290463863,
00425                         -0.915111948139284,  -0.528340129596858,  -0.269626682252082,
00426                         -0.660434383303427,  -0.512000835787190,  -0.223412180441618,
00427                         -0.239932664820086,  -0.507701932958193,  -0.211253047073435,
00428                         0.239932664820086,  -0.507701932958193,  -0.211253047073435,
00429                         0.660434383303426,  -0.512000835787190,  -0.223412180441618,
00430                         0.915111948139284,  -0.528340129596858,  -0.269626682252082,
00431                         -0.773622922202284,  -0.315952535579882,  -0.223412180441618,
00432                         -0.421605613935553,  -0.243414114697549,  -0.172119771139157,
00433                         -0.000000000000000,  -0.224211101329670,  -0.158541190167514,
00434                         0.421605613935553,  -0.243414114697549,  -0.172119771139157,
00435                         0.773622922202284,  -0.315952535579882,  -0.223412180441618,
00436                         -0.559649103902302,   0.046063183547205,  -0.211253047073435,
00437                         -0.194172509561981,   0.112105550664835,  -0.158541190167514,
00438                         0.194172509561981,   0.112105550664835,  -0.158541190167514,
00439                         0.559649103902302,   0.046063183547205,  -0.211253047073435,
00440                         -0.319716439082216,   0.461638749410988,  -0.211253047073435,
00441                         -0.000000000000000,   0.486828229395098,  -0.172119771139157,
00442                         0.319716439082216,   0.461638749410988,  -0.211253047073435,
00443                         -0.113188538898858,   0.827953371367071,  -0.223412180441618,
00444                         0.113188538898858,   0.827953371367071,  -0.223412180441618,
00445                         -0.000000000000000,   1.056680259193716,  -0.269626682252082,
00446                         -0.734424396735357,  -0.424020123154587,   0.025434853622935,
00447                         -0.439014646886819,  -0.392761897021160,   0.113846468290170,
00448                         -0.000000000000000,  -0.384900179459751,   0.136082763487954,
00449                         0.439014646886819,  -0.392761897021160,   0.113846468290170,
00450                         0.734424396735357,  -0.424020123154587,   0.025434853622935,
00451                         -0.559649103902302,  -0.183816888326860,   0.113846468290170,
00452                         -0.194172509561981,  -0.112105550664835,   0.158541190167514,
00453                         0.194172509561981,  -0.112105550664835,   0.158541190167514,
00454                         0.559649103902302,  -0.183816888326860,   0.113846468290170,
00455                         -0.333333333333333,   0.192450089729875,   0.136082763487954,
00456                         -0.000000000000000,   0.224211101329670,   0.158541190167514,
00457                         0.333333333333333,   0.192450089729875,   0.136082763487954,
00458                         -0.120634457015483,   0.576578785348020,   0.113846468290170,
00459                         0.120634457015482,   0.576578785348020,   0.113846468290170,
00460                         -0.000000000000000,   0.848040246309174,   0.025434853622935,
00461                         -0.500000000000000,  -0.288675134594813,   0.408248290463863,
00462                         -0.199081982066733,  -0.254236708399899,   0.505654869247127,
00463                         0.199081982066733,  -0.254236708399899,   0.505654869247127,
00464                         0.500000000000000,  -0.288675134594813,   0.408248290463863,
00465                         -0.319716439082216,  -0.045291699705599,   0.505654869247127,
00466                         -0.000000000000000,  -0.000000000000000,   0.516359313417471,
00467                         0.319716439082216,  -0.045291699705599,   0.505654869247127,
00468                         -0.120634457015483,   0.299528408105498,   0.505654869247127,
00469                         0.120634457015483,  0.299528408105498,   0.505654869247127,
00470                         -0.000000000000000,   0.577350269189626,   0.408248290463863,
00471                         -0.265575603264643,  -0.153330146035039,   0.791061727304791,
00472                         -0.000000000000000,  -0.130698866804872,   0.855072651347100,
00473                         0.265575603264643,  -0.153330146035039,   0.791061727304791,
00474                         -0.113188538898858,   0.065349433402436,   0.855072651347100,
00475                         0.113188538898858,   0.065349433402436,   0.855072651347099,
00476                         0,   0.306660292070078,   0.791061727304791,
00477                         -0.084888051860717,  -0.049010139592768,   1.086123263179808,
00478                         0.084888051860717,  -0.049010139592768,   1.086123263179808,
00479                         0.000000000000000,   0.098020279185535,   1.086123263179808,
00480                         0,                   0,   1.224744871391589
00481     };
00482 
00483     // compare
00484     for (int i=0;i<numPts;i++) {
00485       for (int j=0;j<ptDim;j++) {
00486         int l = ptDim*i+j;
00487         if (std::abs(warpBlendMappedPts(i,j) - points[l]) > INTREPID_TOL) {
00488           errorFlag++;
00489           *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00490           
00491           // Output the multi-index of the value where the error is:
00492           *outStream << " At multi-index { ";
00493           *outStream << i << " ";*outStream << j << " ";
00494           *outStream << "}  computed value: " << warpBlendMappedPts(i,j)
00495                      << " but correct value: " << points[l] << "\n";
00496         }
00497       }
00498     }
00499 
00500 
00501   }
00502   catch (std::exception &err) {
00503     *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
00504     *outStream << err.what() << '\n';
00505     *outStream << "-------------------------------------------------------------------------------" << "\n\n";
00506     errorFlag = -1000;
00507   };
00508 
00509 
00510 
00511   if (errorFlag != 0)
00512     std::cout << "End Result: TEST FAILED\n";
00513   else
00514     std::cout << "End Result: TEST PASSED\n";
00515 
00516   // reset format state of std::cout
00517   std::cout.copyfmt(oldFormatState);
00518 
00519   return errorFlag;
00520 }