Intrepid
http://trilinos.sandia.gov/packages/docs/r10.12/packages/intrepid/test/Discretization/Basis/HGRAD_WEDGE_C1_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 
00048 #include "Intrepid_FieldContainer.hpp"
00049 #include "Intrepid_HGRAD_WEDGE_C1_FEM.hpp"
00050 #include "Teuchos_oblackholestream.hpp"
00051 #include "Teuchos_RCP.hpp"
00052 #include "Teuchos_GlobalMPISession.hpp"
00053 
00054 using namespace std;
00055 using namespace Intrepid;
00056 
00057 #define INTREPID_TEST_COMMAND( S , throwCounter, nException )                                                              \
00058 {                                                                                                                          \
00059   ++nException;                                                                                                            \
00060   try {                                                                                                                    \
00061     S ;                                                                                                                    \
00062   }                                                                                                                        \
00063   catch (std::logic_error err) {                                                                                           \
00064       ++throwCounter;                                                                                                      \
00065       *outStream << "Expected Error " << nException << " -------------------------------------------------------------\n"; \
00066       *outStream << err.what() << '\n';                                                                                    \
00067       *outStream << "-------------------------------------------------------------------------------" << "\n\n";           \
00068   };                                                                                                                       \
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_WEDGE_C1_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 (kjpeter@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   // Define basis and error flag
00109   Basis_HGRAD_WEDGE_C1_FEM<double, FieldContainer<double> > wedgeBasis;
00110   int errorFlag = 0;
00111 
00112   // Initialize throw counter for exception testing
00113   int nException     = 0;
00114   int throwCounter   = 0;
00115 
00116   // Define array containing the 6 vertices of the reference TRIPRISM and some other points.  
00117   FieldContainer<double> wedgeNodes(12, 3);
00118   wedgeNodes(0,0) =  0.0;  wedgeNodes(0,1) =  0.0;  wedgeNodes(0,2) = -1.0;  
00119   wedgeNodes(1,0) =  1.0;  wedgeNodes(1,1) =  0.0;  wedgeNodes(1,2) = -1.0;  
00120   wedgeNodes(2,0) =  0.0;  wedgeNodes(2,1) =  1.0;  wedgeNodes(2,2) = -1.0;
00121   wedgeNodes(3,0) =  0.0;  wedgeNodes(3,1) =  0.0;  wedgeNodes(3,2) =  1.0;  
00122   wedgeNodes(4,0) =  1.0;  wedgeNodes(4,1) =  0.0;  wedgeNodes(4,2) =  1.0;  
00123   wedgeNodes(5,0) =  0.0;  wedgeNodes(5,1) =  1.0;  wedgeNodes(5,2) =  1.0;
00124     
00125   wedgeNodes(6,0) =  0.25; wedgeNodes(6,1) =  0.5;  wedgeNodes(6,2) = -1.0;  
00126   wedgeNodes(7,0) =  0.5;  wedgeNodes(7,1) =  0.25; wedgeNodes(7,2) =  0.0;  
00127   wedgeNodes(8,0) =  0.25; wedgeNodes(8,1) =  0.30; wedgeNodes(8,2) =  1.0;
00128   wedgeNodes(9,0) =  0.3;  wedgeNodes(9,1) =  0.0;  wedgeNodes(9,2) =  0.75;
00129   wedgeNodes(10,0)=  0.0;  wedgeNodes(10,1)=  0.44; wedgeNodes(10,2)= -0.23;  
00130   wedgeNodes(11,0)=  0.4;  wedgeNodes(11,1)=  0.6;  wedgeNodes(11,2)=  0.0;  
00131 
00132 
00133   // Generic array for the output values; needs to be properly resized depending on the operator type
00134   FieldContainer<double> vals;
00135 
00136   try{
00137     // exception #1: CURL cannot be applied to scalar functions
00138     // resize vals to rank-3 container with dimensions (num. points, num. basis functions)
00139     vals.resize(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), 3 );
00140     INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_DIV), throwCounter, nException );
00141 
00142     // exception #2: DIV cannot be applied to scalar functions
00143     // resize vals to rank-2 container with dimensions (num. points, num. basis functions)
00144     vals.resize(wedgeBasis.getCardinality(), wedgeNodes.dimension(0) );
00145     INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_DIV), throwCounter, nException );
00146         
00147     // Exceptions 3-7: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and 
00148     // getDofTag() to access invalid array elements thereby causing bounds check exception
00149     // exception #3
00150     INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(3,0,0), throwCounter, nException );
00151     // exception #4
00152     INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(1,1,1), throwCounter, nException );
00153     // exception #5
00154     INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(0,6,0), throwCounter, nException );
00155     // exception #6
00156     INTREPID_TEST_COMMAND( wedgeBasis.getDofTag(7), throwCounter, nException );
00157     // exception #7
00158     INTREPID_TEST_COMMAND( wedgeBasis.getDofTag(-1), throwCounter, nException );
00159     
00160 #ifdef HAVE_INTREPID_DEBUG
00161     // Exceptions 8-18 test exception handling with incorrectly dimensioned input/output arrays
00162     // exception #8: input points array must be of rank-2
00163     FieldContainer<double> badPoints1(4, 5, 3);
00164     INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
00165     
00166     // exception #9 dimension 1 in the input point array must equal space dimension of the cell
00167     FieldContainer<double> badPoints2(4, 2);
00168     INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
00169     
00170     // exception #10 output values must be of rank-2 for OPERATOR_VALUE
00171     FieldContainer<double> badVals1(4, 3, 1);
00172     INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals1, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
00173     
00174     // exception #11 output values must be of rank-3 for OPERATOR_GRAD
00175     FieldContainer<double> badVals2(4, 3);
00176     INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals2, wedgeNodes, OPERATOR_GRAD), throwCounter, nException );
00177     
00178     // exception #12 output values must be of rank-3 for OPERATOR_D1
00179     INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals2, wedgeNodes, OPERATOR_D1), throwCounter, nException );
00180     
00181     // exception #13 output values must be of rank-3 for OPERATOR_D2
00182     INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals2, wedgeNodes, OPERATOR_D2), throwCounter, nException );
00183     
00184     // exception #14 incorrect 0th dimension of output array (must equal number of basis functions)
00185     FieldContainer<double> badVals3(wedgeBasis.getCardinality() + 1, wedgeNodes.dimension(0));
00186     INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals3, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
00187     
00188     // exception #15 incorrect 1st dimension of output array (must equal number of points)
00189     FieldContainer<double> badVals4(wedgeBasis.getCardinality(), wedgeNodes.dimension(0) + 1);
00190     INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals4, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
00191     
00192     // exception #16: incorrect 2nd dimension of output array (must equal the space dimension)
00193     FieldContainer<double> badVals5(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), 4);
00194     INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals5, wedgeNodes, OPERATOR_GRAD), throwCounter, nException );
00195     
00196     // exception #17: incorrect 2nd dimension of output array (must equal D2 cardinality in 3D)
00197     FieldContainer<double> badVals6(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), 40);
00198     INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals6, wedgeNodes, OPERATOR_D2), throwCounter, nException );
00199     
00200     // exception #18: incorrect 2nd dimension of output array (must equal D3 cardinality in 3D)
00201     INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals6, wedgeNodes, OPERATOR_D3), throwCounter, nException );
00202 #endif
00203 
00204   }
00205   catch (std::logic_error err) {
00206     *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
00207     *outStream << err.what() << '\n';
00208     *outStream << "-------------------------------------------------------------------------------" << "\n\n";
00209     errorFlag = -1000;
00210   };
00211   
00212   // Check if number of thrown exceptions matches the one we expect - 18
00213   if (throwCounter != nException) {
00214     errorFlag++;
00215     *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00216   }
00217   
00218   *outStream \
00219     << "\n"
00220     << "===============================================================================\n"\
00221     << "| TEST 2: correctness of tag to enum and enum to tag lookups                  |\n"\
00222     << "===============================================================================\n";
00223   
00224   try{
00225     std::vector<std::vector<int> > allTags = wedgeBasis.getAllDofTags();
00226     
00227     // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again
00228     for (unsigned i = 0; i < allTags.size(); i++) {
00229       int bfOrd  = wedgeBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
00230       
00231       std::vector<int> myTag = wedgeBasis.getDofTag(bfOrd);
00232        if( !( (myTag[0] == allTags[i][0]) &&
00233               (myTag[1] == allTags[i][1]) &&
00234               (myTag[2] == allTags[i][2]) &&
00235               (myTag[3] == allTags[i][3]) ) ) {
00236         errorFlag++;
00237         *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00238         *outStream << " getDofOrdinal( {" 
00239           << allTags[i][0] << ", " 
00240           << allTags[i][1] << ", " 
00241           << allTags[i][2] << ", " 
00242           << allTags[i][3] << "}) = " << bfOrd <<" but \n";   
00243         *outStream << " getDofTag(" << bfOrd << ") = { "
00244           << myTag[0] << ", " 
00245           << myTag[1] << ", " 
00246           << myTag[2] << ", " 
00247           << myTag[3] << "}\n";        
00248       }
00249     }
00250     
00251     // Now do the same but loop over basis functions
00252     for( int bfOrd = 0; bfOrd < wedgeBasis.getCardinality(); bfOrd++) {
00253       std::vector<int> myTag  = wedgeBasis.getDofTag(bfOrd);
00254       int myBfOrd = wedgeBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
00255       if( bfOrd != myBfOrd) {
00256         errorFlag++;
00257         *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00258         *outStream << " getDofTag(" << bfOrd << ") = { "
00259           << myTag[0] << ", " 
00260           << myTag[1] << ", " 
00261           << myTag[2] << ", " 
00262           << myTag[3] << "} but getDofOrdinal({" 
00263           << myTag[0] << ", " 
00264           << myTag[1] << ", " 
00265           << myTag[2] << ", " 
00266           << myTag[3] << "} ) = " << myBfOrd << "\n";
00267       }
00268     }
00269   }
00270   catch (std::logic_error err){
00271     *outStream << err.what() << "\n\n";
00272     errorFlag = -1000;
00273   };
00274   
00275   *outStream \
00276     << "\n"
00277     << "===============================================================================\n"\
00278     << "| TEST 3: correctness of basis function values                                |\n"\
00279     << "===============================================================================\n";
00280   
00281   outStream -> precision(20);
00282   
00283   // VALUE: Each row gives the 4 correct basis set values at an evaluation point
00284   double basisValues[] = {
00285     1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 
00286     0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 
00287     0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 
00288     0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 
00289     0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 
00290     0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 
00291     //
00292     0.25,     0.25,     0.5,    0.,     0.,     0.,
00293     0.125,    0.25,     0.125,  0.125,  0.25,   0.125,
00294     0.,       0.,       0.,     0.45,   0.25,   0.3,
00295     0.0875,   0.0375,   0,      0.6125, 0.2625, 0,
00296     0.3444,   0,        0.2706, 0.2156, 0,      0.1694,
00297     0.,       0.2,      0.3,    0.,     0.2,    0.3
00298   };
00299   
00300   // GRAD and D1: each row gives the 3 x 4 correct values of the gradients of the 4 basis functions
00301   double basisGrads[] = {
00302     -1., -1., -0.5, 1., 0, 0, 0, 1., 0, 0., 0., 0.5, 0., 0, 0, 0, 0., 0, \
00303     -1., -1., 0., 1., 0, -0.5, 0, 1., 0, 0., 0., 0., 0., 0, 0.5, 0, 0., \
00304     0, -1., -1., 0., 1., 0, 0, 0, 1., -0.5, 0., 0., 0., 0., 0, 0, 0, 0., \
00305     0.5, 0., 0., -0.5, 0., 0, 0, 0, 0., 0, -1., -1., 0.5, 1., 0, 0, 0, \
00306     1., 0, 0., 0., 0., 0., 0, -0.5, 0, 0., 0, -1., -1., 0., 1., 0, 0.5, \
00307     0, 1., 0, 0., 0., 0., 0., 0, 0, 0, 0., -0.5, -1., -1., 0., 1., 0, 0, \
00308     0, 1., 0.5, -1., -1., -0.125, 1., 0, -0.125, 0, 1., -0.25, 0., 0., \
00309     0.125, 0., 0, 0.125, 0, 0., 0.25, -0.5, -0.5, -0.125, 0.5, 0, -0.25, \
00310     0, 0.5, -0.125, -0.5, -0.5, 0.125, 0.5, 0, 0.25, 0, 0.5, 0.125, 0., \
00311     0., -0.225, 0., 0, -0.125, 0, 0., -0.15, -1., -1., 0.225, 1., 0, \
00312     0.125, 0, 1., 0.15, -0.125, -0.125, -0.35, 0.125, 0, -0.15, 0, 0.125, \
00313     0, -0.875, -0.875, 0.35, 0.875, 0, 0.15, 0, 0.875, 0, -0.615, -0.615, \
00314     -0.28, 0.615, 0, 0, 0, 0.615, -0.22, -0.385, -0.385, 0.28, 0.385, 0, \
00315     0, 0, 0.385, 0.22, -0.5, -0.5, 0., 0.5, 0, -0.2, 0, 0.5, -0.3, -0.5, \
00316     -0.5, 0., 0.5, 0, 0.2, 0, 0.5, 0.3
00317   };
00318   //D2: flat array with the values of D2 applied to basis functions. Multi-index is (P,F,K)
00319   double basisD2[] = {
00320     0, 0, 0.5, 0, 0.5, 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, \
00321     -0.5, 0, -0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, \
00322     0, 0.5, 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, \
00323     -0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0.5, \
00324     0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, -0.5, 0, \
00325     0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0.5, 0, 0, 0, \
00326     -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, -0.5, 0, 0, 0, \
00327     0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0.5, 0, 0, 0, -0.5, \
00328     0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, -0.5, 0, 0, 0, 0.5, 0, \
00329     0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0.5, 0, 0, 0, -0.5, 0, 0, 0, \
00330     0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, -0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, \
00331     0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0.5, 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, \
00332     0, -0.5, 0, 0, 0, -0.5, 0, -0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, \
00333     0.5, 0, 0, 0, 0.5, 0, 0.5, 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, \
00334     0, 0, 0, -0.5, 0, -0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, \
00335     0, 0.5, 0, 0.5, 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, \
00336     -0.5, 0, -0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, \
00337     0, 0.5, 0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, \
00338     -0.5, 0, 0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0.5, \
00339     0, 0, 0, -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, -0.5, 0, \
00340     0, 0, 0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0, 0, 0, 0.5, 0, 0.5, 0, 0, 0, \
00341     -0.5, 0, 0, 0, 0, 0, 0, 0, -0.5, 0, 0, 0, -0.5, 0, -0.5, 0, 0, 0, \
00342     0.5, 0, 0, 0, 0, 0, 0, 0, 0.5, 0
00343   };
00344   try{
00345         
00346     // Dimensions for the output arrays:
00347     int numFields = wedgeBasis.getCardinality();
00348     int numPoints = wedgeNodes.dimension(0);
00349     int spaceDim  = wedgeBasis.getBaseCellTopology().getDimension();
00350     int D2Cardin  = Intrepid::getDkCardinality(OPERATOR_D2, spaceDim);
00351     
00352     // Generic array for values, grads, curls, etc. that will be properly sized before each call
00353     FieldContainer<double> vals;
00354     
00355     // Check VALUE of basis functions: resize vals to rank-2 container:
00356     vals.resize(numFields, numPoints);
00357     wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_VALUE);
00358     for (int i = 0; i < numFields; i++) {
00359       for (int j = 0; j < numPoints; j++) {
00360           int l =  i + j * numFields;
00361            if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
00362              errorFlag++;
00363              *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00364 
00365              // Output the multi-index of the value where the error is:
00366              *outStream << " At multi-index { ";
00367              *outStream << i << " ";*outStream << j << " ";
00368              *outStream << "}  computed value: " << vals(i,j)
00369                << " but reference value: " << basisValues[l] << "\n";
00370          }
00371       }
00372     }
00373 
00374     // Check GRAD of basis function: resize vals to rank-3 container
00375     vals.resize(numFields, numPoints, spaceDim);
00376     wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_GRAD);
00377     for (int i = 0; i < numFields; i++) {
00378       for (int j = 0; j < numPoints; j++) {
00379         for (int k = 0; k < spaceDim; k++) {
00380            int l = k + i * spaceDim + j * spaceDim * numFields;
00381            if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
00382              errorFlag++;
00383              *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00384 
00385              // Output the multi-index of the value where the error is:
00386              *outStream << " At multi-index { ";
00387              *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
00388              *outStream << "}  computed grad component: " << vals(i,j,k)
00389                << " but reference grad component: " << basisGrads[l] << "\n";
00390             }
00391          }
00392       }
00393     }
00394 
00395     // Check D1 of basis function (do not resize vals because it has the correct size: D1 = GRAD)
00396     wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_D1);
00397     for (int i = 0; i < numFields; i++) {
00398       for (int j = 0; j < numPoints; j++) {
00399         for (int k = 0; k < spaceDim; k++) {
00400            int l = k + i * spaceDim + j * spaceDim * numFields;
00401            if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
00402              errorFlag++;
00403              *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00404 
00405              // Output the multi-index of the value where the error is:
00406              *outStream << " At multi-index { ";
00407              *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
00408              *outStream << "}  computed D1 component: " << vals(i,j,k)
00409                << " but reference D1 component: " << basisGrads[l] << "\n";
00410             }
00411          }
00412       }
00413     }
00414 
00415     // Check D2 of basis function
00416     vals.resize(numFields, numPoints, D2Cardin);    
00417     wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_D2);
00418     for (int i = 0; i < numFields; i++) {
00419       for (int j = 0; j < numPoints; j++) {
00420         for (int k = 0; k < D2Cardin; k++) {
00421            int l = k + i * D2Cardin + j * D2Cardin * numFields;
00422            if (std::abs(vals(i,j,k) - basisD2[l]) > INTREPID_TOL) {
00423              errorFlag++;
00424              *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00425 
00426              // Output the multi-index of the value where the error is:
00427              *outStream << " At multi-index { ";
00428              *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
00429              *outStream << "}  computed D2 component: " << vals(i,j,k)
00430                << " but reference D2 component: " << basisD2[l] << "\n";
00431             }
00432          }
00433       }
00434     }
00435 
00436     // Check all higher derivatives - must be zero. 
00437     for(EOperator op = OPERATOR_D3; op < OPERATOR_MAX; op++) {
00438       
00439       // The last dimension is the number of kth derivatives and needs to be resized for every Dk
00440       int DkCardin  = Intrepid::getDkCardinality(op, spaceDim);
00441       vals.resize(numFields, numPoints, DkCardin);    
00442 
00443       wedgeBasis.getValues(vals, wedgeNodes, op);
00444       for (int i = 0; i < vals.size(); i++) {
00445         if (std::abs(vals[i]) > INTREPID_TOL) {
00446           errorFlag++;
00447           *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00448           
00449           // Get the multi-index of the value where the error is and the operator order
00450           std::vector<int> myIndex;
00451           vals.getMultiIndex(myIndex,i);
00452           int ord = Intrepid::getOperatorOrder(op);
00453           *outStream << " At multi-index { ";
00454           for(int j = 0; j < vals.rank(); j++) {
00455             *outStream << myIndex[j] << " ";
00456           }
00457           *outStream << "}  computed D"<< ord <<" component: " << vals[i] 
00458             << " but reference D" << ord << " component:  0 \n";
00459         }
00460       }
00461     }    
00462   }
00463   
00464   // Catch unexpected errors
00465   catch (std::logic_error err) {
00466     *outStream << err.what() << "\n\n";
00467     errorFlag = -1000;
00468   };
00469   
00470   if (errorFlag != 0)
00471     std::cout << "End Result: TEST FAILED\n";
00472   else
00473     std::cout << "End Result: TEST PASSED\n";
00474   
00475   // reset format state of std::cout
00476   std::cout.copyfmt(oldFormatState);
00477   
00478   return errorFlag;
00479 }