Intrepid
/usr/src/RPM/BUILD/trilinos10-10.6.4/packages/intrepid/test/Discretization/Basis/HGRAD_TRI_C2_FEM/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 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Pavel Bochev  (pbboche@sandia.gov), 
00025 //                    Denis Ridzal  (dridzal@sandia.gov),
00026 //                    Kara Peterson (kjpeter@sandia.gov).
00027 //
00028 // ************************************************************************
00029 // @HEADER
00030 
00036 #include "Intrepid_FieldContainer.hpp"
00037 #include "Intrepid_HGRAD_TRI_C2_FEM.hpp"
00038 #include "Teuchos_oblackholestream.hpp"
00039 #include "Teuchos_RCP.hpp"
00040 #include "Teuchos_GlobalMPISession.hpp"
00041 
00042 using namespace std;
00043 using namespace Intrepid;
00044 
00045 #define INTREPID_TEST_COMMAND( S , throwCounter, nException )                                                              \
00046 {                                                                                                                          \
00047   ++nException;                                                                                                            \
00048   try {                                                                                                                    \
00049     S ;                                                                                                                    \
00050   }                                                                                                                        \
00051   catch (std::logic_error err) {                                                                                           \
00052       ++throwCounter;                                                                                                      \
00053       *outStream << "Expected Error " << nException << " -------------------------------------------------------------\n"; \
00054       *outStream << err.what() << '\n';                                                                                    \
00055       *outStream << "-------------------------------------------------------------------------------" << "\n\n";           \
00056   };                                                                                                                       \
00057 }
00058 int main(int argc, char *argv[]) {
00059 
00060   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00061 
00062   // This little trick lets us print to std::cout only if
00063   // a (dummy) command-line argument is provided.
00064   int iprint     = argc - 1;
00065   Teuchos::RCP<std::ostream> outStream;
00066   Teuchos::oblackholestream bhs; // outputs nothing
00067   if (iprint > 0)
00068     outStream = Teuchos::rcp(&std::cout, false);
00069   else
00070     outStream = Teuchos::rcp(&bhs, false);
00071 
00072   // Save the format state of the original std::cout.
00073   Teuchos::oblackholestream oldFormatState;
00074   oldFormatState.copyfmt(std::cout);
00075 
00076   *outStream \
00077     << "===============================================================================\n" \
00078     << "|                                                                             |\n" \
00079     << "|                 Unit Test (Basis_HGRAD_TRI_C2_FEM)                          |\n" \
00080     << "|                                                                             |\n" \
00081     << "|     1) Conversion of Dof tags into Dof ordinals and back                    |\n" \
00082     << "|     2) Basis values for VALUE, GRAD, CURL, and Dk operators                 |\n" \
00083     << "|                                                                             |\n" \
00084     << "|  Questions? Contact  Pavel Bochev  (pbboche@sandia.gov),                    |\n" \
00085     << "|                      Denis Ridzal  (dridzal@sandia.gov),                    |\n" \
00086     << "|                      Kara Peterson (dridzal@sandia.gov).                    |\n" \
00087     << "|                                                                             |\n" \
00088     << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00089     << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00090     << "|                                                                             |\n" \
00091     << "===============================================================================\n"\
00092     << "| TEST 1: Basis creation, exception testing                                   |\n"\
00093     << "===============================================================================\n";
00094 
00095   
00096   // Define basis and error flag
00097   Basis_HGRAD_TRI_C2_FEM<double, FieldContainer<double> > triBasis;
00098   int errorFlag = 0;
00099 
00100   // Initialize throw counter for exception testing
00101   int nException     = 0;
00102   int throwCounter   = 0;
00103   
00104   // The unisolvent set for the default Lagrangian basis of degree 2 uses the standard equispaced
00105   // lattice on Triangle consisting of the 3 vertices and the 3 edge midpoints. To check basis
00106   // correctness it suffices to evaluate each basis function at the lattice points, defined below,
00107   // but we throw in an additional random point.
00108   FieldContainer<double> triNodes(7, 2);
00109   triNodes(0,0) =  0.0;     triNodes(0,1) =  0.0;
00110   triNodes(1,0) =  1.0;     triNodes(1,1) =  0.0;
00111   triNodes(2,0) =  0.0;     triNodes(2,1) =  1.0;
00112   triNodes(3,0) =  0.5;     triNodes(3,1) =  0.0;
00113   triNodes(4,0) =  0.5;     triNodes(4,1) =  0.5;
00114   triNodes(5,0) =  0.0;     triNodes(5,1) =  0.5;
00115   triNodes(6,0) =  0.1732;  triNodes(6,1) = 0.6714;
00116   
00117   // Generic array for the output values; needs to be properly resized depending on the operator type
00118   FieldContainer<double> vals;
00119 
00120   try{
00121     // exception #1: DIV cannot be applied to scalar functions
00122     // resize vals to rank-2 container with dimensions (num. points, num. basis functions)
00123     vals.resize(triBasis.getCardinality(), triNodes.dimension(0) );
00124     INTREPID_TEST_COMMAND( triBasis.getValues(vals, triNodes, OPERATOR_DIV), throwCounter, nException );
00125     
00126     // Exceptions 2-6: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and 
00127     // getDofTag() to access invalid array elements thereby causing bounds check exception
00128     // exception #2
00129     INTREPID_TEST_COMMAND( triBasis.getDofOrdinal(2,0,0), throwCounter, nException );
00130     // exception #3
00131     INTREPID_TEST_COMMAND( triBasis.getDofOrdinal(1,1,1), throwCounter, nException );
00132     // exception #4
00133     INTREPID_TEST_COMMAND( triBasis.getDofOrdinal(0,4,0), throwCounter, nException );
00134     // exception #5
00135     INTREPID_TEST_COMMAND( triBasis.getDofTag(6), throwCounter, nException );
00136     // exception #6
00137     INTREPID_TEST_COMMAND( triBasis.getDofTag(-1), throwCounter, nException );
00138 
00139 #ifdef HAVE_INTREPID_DEBUG
00140     // Exceptions 7-17 test exception handling with incorrectly dimensioned input/output arrays
00141     // exception #7: input points array must be of rank-2
00142     FieldContainer<double> badPoints1(4, 5, 3);
00143     INTREPID_TEST_COMMAND( triBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
00144     
00145     // exception #8 dimension 1 in the input point array must equal space dimension of the cell
00146     FieldContainer<double> badPoints2(4, triBasis.getBaseCellTopology().getDimension() + 1);
00147     INTREPID_TEST_COMMAND( triBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
00148     
00149     // exception #9 output values must be of rank-2 for OPERATOR_VALUE
00150     FieldContainer<double> badVals1(4, 3, 1);
00151     INTREPID_TEST_COMMAND( triBasis.getValues(badVals1, triNodes, OPERATOR_VALUE), throwCounter, nException );
00152     
00153     // exception #10 output values must be of rank-3 for OPERATOR_GRAD
00154     FieldContainer<double> badVals2(4, 3);
00155     INTREPID_TEST_COMMAND( triBasis.getValues(badVals2, triNodes, OPERATOR_GRAD), throwCounter, nException );
00156     
00157     // exception #11 output values must be of rank-3 for OPERATOR_CURL
00158     INTREPID_TEST_COMMAND( triBasis.getValues(badVals2, triNodes, OPERATOR_CURL), throwCounter, nException );
00159     
00160     // exception #12 output values must be of rank-3 for OPERATOR_D2
00161     INTREPID_TEST_COMMAND( triBasis.getValues(badVals2, triNodes, OPERATOR_D2), throwCounter, nException );
00162     
00163     // exception #13 incorrect 1st dimension of output array (must equal number of basis functions)
00164     FieldContainer<double> badVals3(triBasis.getCardinality() + 1, triNodes.dimension(0));
00165     INTREPID_TEST_COMMAND( triBasis.getValues(badVals3, triNodes, OPERATOR_VALUE), throwCounter, nException );
00166     
00167     // exception #14 incorrect 0th dimension of output array (must equal number of points)
00168     FieldContainer<double> badVals4(triBasis.getCardinality(), triNodes.dimension(0) + 1);
00169     INTREPID_TEST_COMMAND( triBasis.getValues(badVals4, triNodes, OPERATOR_VALUE), throwCounter, nException );
00170     
00171     // exception #15: incorrect 2nd dimension of output array (must equal the space dimension)
00172     FieldContainer<double> badVals5(triBasis.getCardinality(), triNodes.dimension(0), triBasis.getBaseCellTopology().getDimension() + 1);
00173     INTREPID_TEST_COMMAND( triBasis.getValues(badVals5, triNodes, OPERATOR_GRAD), throwCounter, nException );
00174     
00175     // exception #16: incorrect 2nd dimension of output array (must equal D2 cardinality in 2D)
00176     FieldContainer<double> badVals6(triBasis.getCardinality(), triNodes.dimension(0), 40);
00177     INTREPID_TEST_COMMAND( triBasis.getValues(badVals6, triNodes, OPERATOR_D2), throwCounter, nException );
00178     
00179     // exception #17: incorrect 2nd dimension of output array (must equal D3 cardinality in 2D)
00180     INTREPID_TEST_COMMAND( triBasis.getValues(badVals6, triNodes, OPERATOR_D3), throwCounter, nException );
00181 #endif
00182     
00183   }
00184   catch (std::logic_error err) {
00185     *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
00186     *outStream << err.what() << '\n';
00187     *outStream << "-------------------------------------------------------------------------------" << "\n\n";
00188     errorFlag = -1000;
00189   };
00190   
00191   // Check if number of thrown exceptions matches the one we expect
00192   if (throwCounter != nException) {
00193     errorFlag++;
00194     *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00195   }
00196   
00197   *outStream \
00198     << "\n"
00199     << "===============================================================================\n"\
00200     << "| TEST 2: correctness of tag to enum and enum to tag lookups                  |\n"\
00201     << "===============================================================================\n";
00202 
00203   try{
00204     std::vector<std::vector<int> > allTags = triBasis.getAllDofTags();
00205     
00206     // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again
00207     for (unsigned i = 0; i < allTags.size(); i++) {
00208       int bfOrd  = triBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
00209       
00210       std::vector<int> myTag = triBasis.getDofTag(bfOrd);
00211       if( !( (myTag[0] == allTags[i][0]) &&
00212              (myTag[1] == allTags[i][1]) &&
00213              (myTag[2] == allTags[i][2]) &&
00214              (myTag[3] == allTags[i][3]) ) ) {
00215         errorFlag++;
00216         *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00217         *outStream << " getDofOrdinal( {" 
00218           << allTags[i][0] << ", " 
00219           << allTags[i][1] << ", " 
00220           << allTags[i][2] << ", " 
00221           << allTags[i][3] << "}) = " << bfOrd <<" but \n";   
00222         *outStream << " getDofTag(" << bfOrd << ") = { "
00223           << myTag[0] << ", " 
00224           << myTag[1] << ", " 
00225           << myTag[2] << ", " 
00226           << myTag[3] << "}\n";        
00227       }
00228     }
00229     
00230     // Now do the same but loop over basis functions
00231     for( int bfOrd = 0; bfOrd < triBasis.getCardinality(); bfOrd++) {
00232       std::vector<int> myTag  = triBasis.getDofTag(bfOrd);
00233       int myBfOrd = triBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
00234       if( bfOrd != myBfOrd) {
00235         errorFlag++;
00236         *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00237         *outStream << " getDofTag(" << bfOrd << ") = { "
00238           << myTag[0] << ", " 
00239           << myTag[1] << ", " 
00240           << myTag[2] << ", " 
00241           << myTag[3] << "} but getDofOrdinal({" 
00242           << myTag[0] << ", " 
00243           << myTag[1] << ", " 
00244           << myTag[2] << ", " 
00245           << myTag[3] << "} ) = " << myBfOrd << "\n";
00246       }
00247     }
00248   }
00249   catch (std::logic_error err){
00250     *outStream << err.what() << "\n\n";
00251     errorFlag = -1000;
00252   };
00253   
00254   *outStream \
00255     << "\n"
00256     << "===============================================================================\n"\
00257     << "| TEST 3: correctness of basis function values                                |\n"\
00258     << "===============================================================================\n";
00259 
00260   outStream -> precision(20);
00261 
00262   // VALUE: Correct basis values in (F,P) format: each row gives the 6 correct basis values at
00263   // the lattice points (3 vertices and 3 edge midpoints)
00264   double basisValues[] = {
00265     1.0, 0.0, 0.0,    0.0, 0.0, 0.0,    -0.10710168,
00266     0.0, 1.0, 0.0,    0.0, 0.0, 0.0,    -0.11320352,
00267     0.0, 0.0, 1.0,    0.0, 0.0, 0.0,     0.23015592,
00268     0.0, 0.0, 0.0,    1.0, 0.0, 0.0,     0.10766112,
00269     0.0, 0.0, 0.0,    0.0, 1.0, 0.0,     0.46514592,
00270     0.0, 0.0, 0.0,    0.0, 0.0, 1.0,     0.41734224
00271   };
00272 
00273   // GRAD and D1: Correct gradients and D1 in (F,P,D) format: each row contains 6x2 values of
00274   // gradients of basis functions. Same array is used to check correctness of CURL.
00275   double basisGrads[] = {
00276     //   V0            V1           V2                   M0             M1            M2          RP
00277     -3.0, -3.0,     1.0, 1.0,   1.0, 1.0,           -1.0, -1.0,     1.0, 1.0,     -1.0,-1.0,     0.3784,  0.3784,
00278     -1.0,  0.0,     3.0, 0.0,  -1.0, 0.0,            1.0,  0.0,     1.0, 0.0,     -1.0, 0.0,    -0.3072,  0.0,
00279      0.0, -1.0,     0.0,-1.0,   0.0, 3.0,            0.0, -1.0,     0.0, 1.0,      0.0, 1.0,     0.0,     1.6856,
00280      4.0,  0.0,    -4.0,-4.0,   0.0, 0.0,            0.0, -2.0,    -2.0,-2.0,      2.0, 0.0,    -0.0712, -0.6928, 
00281      0.0,  0.0,     0.0, 4.0,   4.0, 0.0,            0.0,  2.0,     2.0, 2.0,      2.0, 0.0,     2.6856,  0.6928,
00282      0.0,  4.0,     0.0, 0.0,  -4.0,-4.0,            0.0,  2.0,    -2.0,-2.0,     -2.0, 0.0,    -2.6856, -2.064
00283   };
00284   
00285   // D2: Correct multiset of second order partials in (F,P,Dk) format. D2 cardinality = 3 for 2D 
00286   double basisD2[] = {
00287     //   V0                  V1              V2                   M0              M1                M2            RP
00288     4.0, 4.0, 4.0,    4.0, 4.0, 4.0,    4.0, 4.0, 4.0,    4.0, 4.0, 4.0,    4.0, 4.0, 4.0,    4.0, 4.0, 4.0,   4.0, 4.0, 4.0,
00289     4.0, 0.0, 0.0,    4.0, 0.0, 0.0,    4.0, 0.0, 0.0,    4.0, 0.0, 0.0,    4.0, 0.0, 0.0,    4.0, 0.0, 0.0,   4.0, 0.0, 0.0,
00290     0.0, 0.0, 4.0,    0.0, 0.0, 4.0,    0.0, 0.0, 4.0,    0.0, 0.0, 4.0,    0.0, 0.0, 4.0,    0.0, 0.0, 4.0,   0.0, 0.0, 4.0,
00291    -8.0,-4.0, 0.0,   -8.0,-4.0, 0.0,   -8.0,-4.0, 0.0,   -8.0,-4.0, 0.0,   -8.0,-4.0, 0.0,   -8.0,-4.0, 0.0,  -8.0,-4.0, 0.0,
00292     0.0, 4.0, 0.0,    0.0, 4.0, 0.0,    0.0, 4.0, 0.0,    0.0, 4.0, 0.0,    0.0, 4.0, 0.0,    0.0, 4.0, 0.0,   0.0, 4.0, 0.0,
00293     0.0,-4.0,-8.0,    0.0,-4.0,-8.0,    0.0,-4.0,-8.0,    0.0,-4.0,-8.0,    0.0,-4.0,-8.0,    0.0,-4.0,-8.0,   0.0,-4.0,-8.0
00294   };
00295   try{
00296     
00297     // Dimensions for the output arrays:
00298     int numFields = triBasis.getCardinality();
00299     int numPoints = triNodes.dimension(0);
00300     int spaceDim  = triBasis.getBaseCellTopology().getDimension();
00301     
00302     // Generic array for values, grads, curls, etc. that will be properly sized before each call
00303     FieldContainer<double> vals;
00304     
00305     // Check VALUE of basis functions: resize vals to rank-2 container:
00306     vals.resize(numFields, numPoints);
00307     triBasis.getValues(vals, triNodes, OPERATOR_VALUE);
00308     for (int i = 0; i < numFields; i++) {
00309       for (int j = 0; j < numPoints; j++) {
00310         
00311         // Compute offset for (F,P) container
00312           int l =  j + i * numPoints;
00313            if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
00314              errorFlag++;
00315              *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00316 
00317              // Output the multi-index of the value where the error is:
00318              *outStream << " At multi-index { ";
00319              *outStream << i << " ";*outStream << j << " ";
00320              *outStream << "}  computed value: " << vals(i,j)
00321                << " but reference value: " << basisValues[l] << "\n";
00322          }
00323       }
00324     }
00325     
00326     // Check GRAD of basis function: resize vals to rank-3 container
00327     vals.resize(numFields, numPoints, spaceDim);
00328     triBasis.getValues(vals, triNodes, OPERATOR_GRAD);
00329     for (int i = 0; i < numFields; i++) {
00330       for (int j = 0; j < numPoints; j++) {
00331         for (int k = 0; k < spaceDim; k++) {
00332           // basisGrads is (F,P,D), compute offset:
00333            int l = k + j * spaceDim + i * spaceDim * numPoints;
00334            if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
00335              errorFlag++;
00336              *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00337 
00338              // Output the multi-index of the value where the error is:
00339              *outStream << " At multi-index { ";
00340              *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
00341              *outStream << "}  computed grad component: " << vals(i,j,k)
00342                << " but reference grad component: " << basisGrads[l] << "\n";
00343             }
00344          }
00345       }
00346     }
00347 
00348     // Check D1 of basis function (do not resize vals because it has the correct size: D1 = GRAD)
00349     triBasis.getValues(vals, triNodes, OPERATOR_D1);
00350     for (int i = 0; i < numFields; i++) {
00351       for (int j = 0; j < numPoints; j++) {
00352         for (int k = 0; k < spaceDim; k++) {
00353           // basisGrads is (F,P,D), compute offset:
00354           int l = k + j * spaceDim + i * spaceDim * numPoints;
00355            if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
00356              errorFlag++;
00357              *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00358 
00359              // Output the multi-index of the value where the error is:
00360              *outStream << " At multi-index { ";
00361              *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
00362              *outStream << "}  computed D1 component: " << vals(i,j,k)
00363                << " but reference D1 component: " << basisGrads[l] << "\n";
00364             }
00365          }
00366       }
00367     }
00368 
00369     // Check CURL of basis function: resize vals just for illustration! 
00370     vals.resize(numFields, numPoints, spaceDim);
00371     triBasis.getValues(vals, triNodes, OPERATOR_CURL);
00372     for (int i = 0; i < numFields; i++) {
00373       for (int j = 0; j < numPoints; j++) {
00374         // We will use "rotated" basisGrads to check CURL: get offsets to extract (u_y, -u_x)
00375         int curl_0 = 1 + j * spaceDim + i * spaceDim * numPoints;               // position of y-derivative
00376         int curl_1 = 0 + j * spaceDim + i * spaceDim * numPoints;               // position of x-derivative
00377         
00378         double curl_value_0 = basisGrads[curl_0];
00379         double curl_value_1 =-basisGrads[curl_1];
00380         if (std::abs(vals(i,j,0) - curl_value_0) > INTREPID_TOL) {
00381           errorFlag++;
00382           *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00383           // Output the multi-index of the value where the error is:
00384           *outStream << " At multi-index { ";
00385           *outStream << i << " ";*outStream << j << " ";*outStream << 0 << " ";
00386           *outStream << "}  computed curl component: " << vals(i,j,0)
00387             << " but reference curl component: " << curl_value_0 << "\n";
00388         }
00389         if (std::abs(vals(i,j,1) - curl_value_1) > INTREPID_TOL) {
00390           errorFlag++;
00391           *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00392           // Output the multi-index of the value where the error is:
00393           *outStream << " At multi-index { ";
00394           *outStream << i << " ";*outStream << j << " ";*outStream << 1 << " ";
00395           *outStream << "}  computed curl component: " << vals(i,j,1)
00396             << " but reference curl component: " << curl_value_1 << "\n";
00397         }
00398       }
00399     }
00400     
00401     // Check D2 of basis function
00402     int D2cardinality = Intrepid::getDkCardinality(OPERATOR_D2, spaceDim);
00403     vals.resize(numFields, numPoints, D2cardinality);
00404     triBasis.getValues(vals, triNodes, OPERATOR_D2);
00405     for (int i = 0; i < numFields; i++) {
00406       for (int j = 0; j < numPoints; j++) {
00407         for (int k = 0; k < D2cardinality; k++) {
00408           
00409           // basisD2 is (F,P,Dk), compute offset:
00410           int l = k + j * D2cardinality + i * D2cardinality * numPoints;
00411           if (std::abs(vals(i,j,k) - basisD2[l]) > INTREPID_TOL) {
00412             errorFlag++;
00413             *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00414             
00415             // Output the multi-index of the value where the error is:
00416             *outStream << " At multi-index { ";
00417             *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
00418             *outStream << "}  computed D2 component: " << vals(i,j,k)
00419               << " but reference D2 component: " << basisD2[l] << "\n";
00420           }
00421         }
00422       }
00423     }
00424           
00425     // Check all higher derivatives - must be zero. 
00426     for(EOperator op = OPERATOR_D3; op < OPERATOR_MAX; op++) {
00427       
00428       // The last dimension is the number of kth derivatives and needs to be resized for every Dk
00429       int DkCardin  = Intrepid::getDkCardinality(op, spaceDim);
00430       vals.resize(numFields, numPoints, DkCardin);    
00431       
00432       triBasis.getValues(vals, triNodes, op);
00433       for (int i = 0; i < vals.size(); i++) {
00434         if (std::abs(vals[i]) > INTREPID_TOL) {
00435           errorFlag++;
00436           *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00437           
00438           // Get the multi-index of the value where the error is and the operator order
00439           std::vector<int> myIndex;
00440           vals.getMultiIndex(myIndex,i);
00441           int ord = Intrepid::getOperatorOrder(op);
00442           *outStream << " At multi-index { ";
00443           for(int j = 0; j < vals.rank(); j++) {
00444             *outStream << myIndex[j] << " ";
00445           }
00446           *outStream << "}  computed D"<< ord <<" component: " << vals[i] 
00447             << " but reference D" << ord << " component:  0 \n";
00448         }
00449       }
00450     }    
00451   }
00452 
00453   // Catch unexpected errors
00454   catch (std::logic_error err) {
00455     *outStream << err.what() << "\n\n";
00456     errorFlag = -1000;
00457   };
00458 
00459   if (errorFlag != 0)
00460     std::cout << "End Result: TEST FAILED\n";
00461   else
00462     std::cout << "End Result: TEST PASSED\n";
00463 
00464   // reset format state of std::cout
00465   std::cout.copyfmt(oldFormatState);
00466 
00467   return errorFlag;
00468 }