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