Intrepid
http://trilinos.sandia.gov/packages/docs/r10.12/packages/intrepid/test/Discretization/Basis/HGRAD_QUAD_Cn_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_QUAD_Cn_FEM.hpp"
00050 #include "Teuchos_oblackholestream.hpp"
00051 #include "Teuchos_RCP.hpp"
00052 #include "Teuchos_GlobalMPISession.hpp"
00053 #include "Intrepid_PointTools.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 
00072 int main(int argc, char *argv[]) {
00073 
00074   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00075   
00076   // This little trick lets us print to std::cout only if
00077   // a (dummy) command-line argument is provided.
00078   int iprint     = argc - 1;
00079   Teuchos::RCP<std::ostream> outStream;
00080   Teuchos::oblackholestream bhs; // outputs nothing
00081   if (iprint > 0)
00082     outStream = Teuchos::rcp(&std::cout, false);
00083   else
00084     outStream = Teuchos::rcp(&bhs, false);
00085   
00086   // Save the format state of the original std::cout.
00087   Teuchos::oblackholestream oldFormatState;
00088   oldFormatState.copyfmt(std::cout);
00089   
00090   *outStream \
00091     << "===============================================================================\n" \
00092     << "|                                                                             |\n" \
00093     << "|                 Unit Test (Basis_HGRAD_QUAD_Cn_FEM)                         |\n" \
00094     << "|                                                                             |\n" \
00095     << "|     1) Conversion of Dof tags into Dof ordinals and back                    |\n" \
00096     << "|     2) Basis values for VALUE, GRAD, CURL, and Dk operators                 |\n" \
00097     << "|                                                                             |\n" \
00098     << "|  Questions? Contact  Pavel Bochev  (pbboche@sandia.gov),                    |\n" \
00099     << "|                      Robert Kirby  (robert.c.kirby@ttu.edu),                |\n" \
00100     << "|                      Denis Ridzal  (dridzal@sandia.gov),                    |\n" \
00101     << "|                      Kara Peterson (kjpeter@sandia.gov).                    |\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     << "| TEST 1: Basis creation, exception testing                                   |\n"\
00108     << "===============================================================================\n";
00109   
00110   // Define basis and error flag
00111   // get points for basis
00112   const int deg=2;
00113   shards::CellTopology line(shards::getCellTopologyData< shards::Line<> >()); 
00114   FieldContainer<double> pts(PointTools::getLatticeSize(line,deg),1);
00115   PointTools::getLattice<double,FieldContainer<double> >(pts,line,deg);
00116 
00117   Basis_HGRAD_QUAD_Cn_FEM<double, FieldContainer<double> > quadBasis(deg,deg,pts,pts);
00118   int errorFlag = 0;
00119 
00120   // Initialize throw counter for exception testing
00121   int nException     = 0;
00122   int throwCounter   = 0;
00123 
00124   // Array with the 4 vertices, 4 edge midpoints, center of the reference QUAD and a random point.  
00125   FieldContainer<double> quadNodes(10, 2);
00126   quadNodes(0,0) = -1.0;  quadNodes(0,1) = -1.0;
00127   quadNodes(1,0) =  1.0;  quadNodes(1,1) = -1.0;
00128   quadNodes(2,0) =  1.0;  quadNodes(2,1) =  1.0;
00129   quadNodes(3,0) = -1.0;  quadNodes(3,1) =  1.0;
00130   // edge midpoints
00131   quadNodes(4,0) =  0.0;  quadNodes(4,1) = -1.0;
00132   quadNodes(5,0) =  1.0;  quadNodes(5,1) =  0.0;
00133   quadNodes(6,0) =  0.0;  quadNodes(6,1) =  1.0;
00134   quadNodes(7,0) = -1.0;  quadNodes(7,1) =  0.0;
00135   // center & random point
00136   quadNodes(8,0) =  0.0;  quadNodes(8,1) =  0.0;
00137   quadNodes(9,0) =1./3.;  quadNodes(9,1) =-3./5.;
00138   
00139   // Generic array for the output values; needs to be properly resized depending on the operator type
00140   FieldContainer<double> vals;
00141 
00142   try{
00143     // exception #1: DIV cannot be applied to scalar functions
00144     // resize vals to rank-2 container with dimensions (num. points, num. basis functions)
00145     vals.resize(quadBasis.getCardinality(), quadNodes.dimension(0));
00146     INTREPID_TEST_COMMAND( quadBasis.getValues(vals, quadNodes, OPERATOR_DIV), throwCounter, nException );
00147         
00148     // Exceptions 2-6: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and 
00149     // getDofTag() to access invalid array elements thereby causing bounds check exception
00150     // exception #2
00151      INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(3,0,0), throwCounter, nException );
00152      // exception #3
00153      INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(1,1,1), throwCounter, nException );
00154     // exception #4
00155     INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(0,4,0), throwCounter, nException );
00156     // exception #5
00157     INTREPID_TEST_COMMAND( quadBasis.getDofTag(10), throwCounter, nException );
00158     // exception #6
00159     INTREPID_TEST_COMMAND( quadBasis.getDofTag(-1), throwCounter, nException );
00160     
00161 #ifdef HAVE_INTREPID_DEBUG
00162     // Exceptions 7- test exception handling with incorrectly dimensioned input/output arrays
00163     // exception #7: input points array must be of rank-2
00164     FieldContainer<double> badPoints1(4, 5, 3);
00165     INTREPID_TEST_COMMAND( quadBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
00166     
00167     // exception #8 dimension 1 in the input point array must equal space dimension of the cell
00168     FieldContainer<double> badPoints2(4, quadBasis.getBaseCellTopology().getDimension() + 1);
00169     INTREPID_TEST_COMMAND( quadBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
00170     
00171     // exception #9 output values must be of rank-2 for OPERATOR_VALUE
00172     FieldContainer<double> badVals1(4, 3, 1);
00173     INTREPID_TEST_COMMAND( quadBasis.getValues(badVals1, quadNodes, OPERATOR_VALUE), throwCounter, nException );
00174 
00175     // exception #10 output values must be of rank-3 for OPERATOR_GRAD
00176     FieldContainer<double> badVals2(4, 3);
00177     INTREPID_TEST_COMMAND( quadBasis.getValues(badVals2, quadNodes, OPERATOR_GRAD), throwCounter, nException );
00178     
00179     // exception #11 output values must be of rank-3 for OPERATOR_CURL
00180     INTREPID_TEST_COMMAND( quadBasis.getValues(badVals2, quadNodes, OPERATOR_CURL), throwCounter, nException );
00181     
00182     // exception #12 output values must be of rank-3 for OPERATOR_D2
00183     INTREPID_TEST_COMMAND( quadBasis.getValues(badVals2, quadNodes, OPERATOR_D2), throwCounter, nException );
00184     
00185     // exception #13 incorrect 0th dimension of output array (must equal number of basis functions)
00186     FieldContainer<double> badVals3(quadBasis.getCardinality() + 1, quadNodes.dimension(0));
00187     INTREPID_TEST_COMMAND( quadBasis.getValues(badVals3, quadNodes, OPERATOR_VALUE), throwCounter, nException );
00188     
00189     // exception #14 incorrect 1st dimension of output array (must equal number of points in quadNodes)
00190     FieldContainer<double> badVals4(quadBasis.getCardinality(), quadNodes.dimension(0) + 1);
00191     INTREPID_TEST_COMMAND( quadBasis.getValues(badVals4, quadNodes, OPERATOR_VALUE), throwCounter, nException );
00192     
00193     // exception #15: incorrect 2nd dimension of output array (must equal the space dimension)
00194     FieldContainer<double> badVals5(quadBasis.getCardinality(), quadNodes.dimension(0), quadBasis.getBaseCellTopology().getDimension() + 1);
00195     INTREPID_TEST_COMMAND( quadBasis.getValues(badVals5, quadNodes, OPERATOR_GRAD), throwCounter, nException );
00196     
00197     // exception #16: incorrect 2nd dimension of output array (must equal D2 cardinality in 2D)
00198     FieldContainer<double> badVals6(quadBasis.getCardinality(), quadNodes.dimension(0), 40);
00199     INTREPID_TEST_COMMAND( quadBasis.getValues(badVals6, quadNodes, OPERATOR_D2), throwCounter, nException );
00200     
00201     // exception #17: incorrect 2nd dimension of output array (must equal D3 cardinality in 2D)
00202     INTREPID_TEST_COMMAND( quadBasis.getValues(badVals6, quadNodes, OPERATOR_D3), throwCounter, nException );
00203 #endif
00204 
00205   }
00206   catch (std::logic_error err) {
00207     *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
00208     *outStream << err.what() << '\n';
00209     *outStream << "-------------------------------------------------------------------------------" << "\n\n";
00210     errorFlag = -1000;
00211   };
00212   
00213   // Check if number of thrown exceptions matches the one we expect 
00214   if (throwCounter != nException) {
00215     errorFlag++;
00216     *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00217   }
00218   
00219   *outStream \
00220     << "\n"
00221     << "===============================================================================\n"\
00222     << "| TEST 2: correctness of tag to enum and enum to tag lookups                  |\n"\
00223     << "===============================================================================\n";
00224   
00225   try{
00226     std::vector<std::vector<int> > allTags = quadBasis.getAllDofTags();
00227     
00228     // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again
00229     for (unsigned i = 0; i < allTags.size(); i++) {
00230       int bfOrd  = quadBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
00231       
00232       std::vector<int> myTag = quadBasis.getDofTag(bfOrd);
00233        if( !( (myTag[0] == allTags[i][0]) &&
00234               (myTag[1] == allTags[i][1]) &&
00235               (myTag[2] == allTags[i][2]) &&
00236               (myTag[3] == allTags[i][3]) ) ) {
00237         errorFlag++;
00238         *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00239         *outStream << " getDofOrdinal( {" 
00240           << allTags[i][0] << ", " 
00241           << allTags[i][1] << ", " 
00242           << allTags[i][2] << ", " 
00243           << allTags[i][3] << "}) = " << bfOrd <<" but \n";   
00244         *outStream << " getDofTag(" << bfOrd << ") = { "
00245           << myTag[0] << ", " 
00246           << myTag[1] << ", " 
00247           << myTag[2] << ", " 
00248           << myTag[3] << "}\n";        
00249       }
00250     }
00251     
00252     // Now do the same but loop over basis functions
00253     for( int bfOrd = 0; bfOrd < quadBasis.getCardinality(); bfOrd++) {
00254       std::vector<int> myTag  = quadBasis.getDofTag(bfOrd);
00255       int myBfOrd = quadBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
00256       if( bfOrd != myBfOrd) {
00257         errorFlag++;
00258         *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00259         *outStream << " getDofTag(" << bfOrd << ") = { "
00260           << myTag[0] << ", " 
00261           << myTag[1] << ", " 
00262           << myTag[2] << ", " 
00263           << myTag[3] << "} but getDofOrdinal({" 
00264           << myTag[0] << ", " 
00265           << myTag[1] << ", " 
00266           << myTag[2] << ", " 
00267           << myTag[3] << "} ) = " << myBfOrd << "\n";
00268       }
00269     }
00270   }
00271   catch (std::logic_error err){
00272     *outStream << err.what() << "\n\n";
00273     errorFlag = -1000;
00274   };
00275   
00276   *outStream \
00277     << "\n"
00278     << "===============================================================================\n"\
00279     << "| TEST 3: correctness of basis function values                                |\n"\
00280     << "===============================================================================\n";
00281   
00282   outStream -> precision(20);
00283   
00284   // VALUE: Correct basis values in (F,P) format:
00285   double basisValues[] = {
00286     1, 0, 0, 0, 0, 0, 0, 0, 0, -0.05333333333333334, \
00287     0, 0, 0, 0, 1, 0, 0, 0, 0, 0.4266666666666667, \
00288     0, 1, 0, 0, 0, 0, 0, 0, 0, 0.1066666666666667, \
00289     0, 0, 0, 0, 0, 0, 0, 1, 0, -0.07111111111111112 , \
00290     0, 0, 0, 0, 0, 0, 0, 0, 1, 0.5688888888888890, \
00291     0, 0, 0, 0, 0, 1, 0, 0, 0, 0.1422222222222222 ,\
00292     0, 0, 0, 1, 0, 0, 0, 0, 0, 0.01333333333333333, \
00293     0, 0, 0, 0, 0, 0, 1, 0, 0, -0.1066666666666667, \
00294     0, 0, 1, 0, 0, 0, 0, 0, 0, -0.02666666666666666 };
00295 
00296   // FIXME HERE: needs to be reordered.
00297   
00298   // GRAD and D1: Correct gradients and D1 in (F,P,D) format
00299   // 9 basis functions, each evaluated at 10 points, with two
00300   // components at each point.
00301   // that looks like 10 per to me.
00302   double basisGrads[] = {
00303     //
00304     -1.500000000000000, -1.500000000000000, 0.5000000000000000, 0, 0, 0, 0, 0.5000000000000000, -0.5000000000000000, 0, \
00305     0, 0, 0, 0, 0, -0.5000000000000000, 0, 0, -0.08000000000000002, 0.1222222222222222, \
00306     //
00307     2.000000000000000, 0, -2.000000000000000, 0, 0, 0, 0, 0, 0, -1.500000000000000, \
00308     0, 0, 0, 0.5000000000000000, 0, 0, 0, -0.5000000000000000,  -0.3199999999999999, -0.9777777777777779, \
00309     //
00310     -0.5000000000000000, 0, 1.500000000000000, -1.500000000000000, 0, 0.5000000000000000, 0, 0, 0.5000000000000000, 0, \
00311     0, -0.5000000000000000, 0, 0, 0, 0, 0, 0, 0.3999999999999999, -0.2444444444444444, \
00312     //
00313     0, 2.0, 0, 0, 0, 0, 0, -2.000000000000000, 0, 0, \
00314     0.5000000000000000, 0, 0, 0, -1.50, 0, -0.50, 0, -0.1066666666666667, -0.1333333333333333, \
00315     //
00316     0, 0, 0, 0, 0, 0, 0, 0, 0, 2.0,\
00317     -2.00, 0, 0, -2.0, 2.0, 0, 0, 0, -0.4266666666666667, 1.066666666666667 ,  \
00318     //
00319     0, 0, 0, 2.000000000000000, 0, -2.000000000000000, 0, 0, 0, 0, \
00320     1.5, 0, 0, 0, -0.5, 0, 0.5000000000000000, 0, 0.5333333333333334, 0.2666666666666666 ,  \
00321     //
00322     0, -0.5000000000000000, 0, 0, 0.5000000000000000, 0, -1.500000000000000, 1.500000000000000, 0, 0, \
00323     0, 0, -0.5000000000000000, 0, 0, 0.5000000000000000, 0, 0, 0.02000000000000000, 0.01111111111111112 , \
00324     //
00325     0, 0, 0, 0, -2.0, 0, 2.0, 0, 0, -0.50,                              \
00326     0, 0, 0, 1.5, 0, 0, 0, 0.5000000000000000, 0.07999999999999997, -0.08888888888888888, \
00327     //
00328     0, 0, 0, -0.5000000000000000, 1.500000000000000, 1.500000000000000, -0.5000000000000000, 0, 0, 0, \
00329     0, 0.5000000000000000, 0.5000000000000000, 0, 0, 0, 0, 0, -0.09999999999999998, -0.02222222222222221 \
00330     //
00331   };
00332   
00333   // D2: Correct multiset of second order partials in (F,P,Dk) format. D2 cardinality = 3 for 2D 
00334   // 10 quad points and 3 values per point, so
00335   // each bf consists of 30 values.
00336   double basisD2[] = {
00337     1.0, 2.25, 1.0, 1.0, -0.75, 0, 0, 0.25, 0,  0, -0.75, 1.0, 1.0,  0.75, 0, 0, -0.25, 0, 0,  -0.25, 0, 0, 0.75, 1.0, 0,  0.25, 0, 0.48, 0.1833333333333334, -0.1111111111111111,
00338     //
00339     -2.0, -3.0, 0, -2.0, 3.0, 0, 0, -1.0, 0, \
00340     0, 1.0, 0, -2.0, 0, 1.0, 0, \
00341     1.0, 0, 0, 0, 1.0, 0, -1.0, \
00342     0, 0, 0, 1.0, -0.96, 0.7333333333333332, \
00343     0.8888888888888890, \
00344     //
00345     1.0, 0.75, 0, 1.0, -2.25, 1.0, 0, 0.75, 1.0, 0, -0.25, 0, \
00346     1.0, -0.75, 0, 0, -0.75, 1.0, 0, 0.25, 0, 0, 0.25, \
00347     0, 0, -0.25, 0, 0.48, -0.9166666666666666, 0.2222222222222222, 
00348     //
00349     0, -3.0, -2.0, 0, 1.0, 0, 0, -1.0, 0, 0, 3.0, \
00350     -2.0, 0, -1.0, 0, 1.0, 0, 0, 0, 1.0, 0, 1.0, 0, -2.0, \
00351     1.0, 0, 0, 0.6400000000000001, -0.2000000000000001,  0.2222222222222222, \
00352     //
00353     0, 4.0, 0, 0, -4.0, 0, 0, 4.0, 0, 0, -4.0, 0, 0, 0, \
00354     -2.0, -2.0, 0, 0, 0, 0, -2.0, -2.0, 0, 0, -2.0, 0, \
00355     -2.0, -1.280000000000000, -0.7999999999999998, -1.777777777777778 ,
00356     //
00357     0, -1.0, 0, 0, 3.0, -2.0, 0, -3.0, -2.0, 0, \
00358     1.0, 0, 0, 1.0, 0, 1.0, 0, -2.0, 0, -1.0, 0, 1.0, 0, \
00359     0, 1.0, 0, 0, 0.6400000000000001, 1.0, -0.4444444444444444, \
00360     //
00361     0, 0.75, 1.0, 0, -0.25, 0, 1.0, 0.75, 0, 1.0, -2.25, 1.0, 0, \
00362     0.25, 0, 0, 0.25, 0, 1.0, -0.75, 0, 0, -0.75, 1.0, 0, \
00363     -0.25, 0, -0.12, 0.01666666666666666,  -0.1111111111111111, \
00364     //
00365     0, -1.0, 0, 0, 1.0, 0, -2.0, -3.0, 0, -2.0, 3.0, 0, 0, 0, 1.0, 0, -1.0, \
00366     0, -2.0, 0, 1.0, 0, 1.0, 0, \
00367     0, 0, 1.0, 0.24, 0.06666666666666665,  0.8888888888888890,  \
00368     //
00369     0, 0.25, 0, 0, -0.75, 1.0, 1.0, 2.25, 1.0, 1.0, \
00370     -0.75, 0, 0, -0.25, 0, 0, 0.75, 1.0, 1.0, \
00371     0.75, 0, 0, -0.25, 0, 0, 0.25, 0, -0.12, -0.08333333333333331, 0.2222222222222222 \
00372   };
00373   
00374   //D3: Correct multiset of second order partials in (F,P,Dk) format. D3 cardinality = 4 for 2D
00375   double basisD3[] = {
00376     0, -1.5, -1.5, 0, 0, -1.5, 0.5, 0, 0, 0.5, 
00377     0.5, 0, 0, 0.5, -1.5, 0, 0, -1.5, -0.5, 0, 
00378     0, -0.5, 0.5, 0, 0, 0.5, -0.5, 0, 0, -0.5, 
00379     -1.5, 0, 0, -0.5, -0.5, 0, 0, -1.1, -0.1666666666666667, 0, 
00380     //
00381     0, 3.0, 2.0, 0, 0, 3.0, -2.0, 0, 0, -1.0, 
00382     -2.0, 0, 0, -1.0, 2.0, 0, 0, 3.0, 0, 0, 
00383     0, 1.0, -2.0, 0, 0, -1.0, 0, 0, 0, 1.0,
00384     2.0, 0, 0, 1.0, 0, 0, 0, 2.2, -0.6666666666666665, 0, 
00385     //
00386     0, -1.5, -0.5, 0, 0, -1.5, 1.5, 0, 0, 0.5, 
00387     1.5, 0, 0, 0.5, -0.5, 0, 0, -1.5, 0.5, 0, 
00388     0, -0.5, 1.5, 0, 0, 0.5, 0.5, 0, 0, -0.5, 
00389     -0.5, 0, 0, -0.5, 0.5, 0, 0, -1.1, 0.8333333333333333, 0,
00390     //
00391     0, 2.0, 3.0, 0, 0, 2.0, -1.0, 0, 0, -2.0,
00392     -1.0, 0, 0, -2.0, 3.0, 0, 0, 2.0,  1.0, 0, 
00393     0, 0, -1.0, 0, 0, -2.0, 1.0, 0, 0, 0, 
00394     3.0, 0, 0, 0, 1.0, 0, 0, 1.2, 0.3333333333333334, 0, 
00395     //
00396     0, -4.0, -4.0, 0, 0, -4.0, 4.0, 0, 0, 4.0, 
00397     4.0, 0, 0, 4.0, -4.0, 0, 0, -4.0, 0, 0, 
00398     0, 0, 4.0, 0, 0, 4.0, 0, 0, 0, 0, 
00399     -4.0, 0, 0, 0, 0, 0, 0, -2.40, 1.333333333333333, 0,
00400     //
00401     0, 2.0, 1.0, 0, 0, 2.0, -3.0, 0, 0, -2.0, 
00402     -3.0, 0, 0, -2.0, 1.0, 0, 0, 2.0, -1.0, 0, 
00403     0, 0, -3.0, 0, 0, -2.0, -1.0, 0, 0, 0, 
00404     1.0, 0, 0, 0, -1.0, 0, 0, 1.2, -1.666666666666667, 0 ,
00405     //
00406     0, -0.5, -1.5, 0, 0, -0.5, 0.5, 0, 0, 1.5, 
00407     0.5, 0, 0, 1.5, -1.5, 0, 0, -0.5, -0.5, 0, 
00408     0, 0.5, 0.5, 0, 0, 1.5, -0.5, 0, 0, 0.5,
00409     -1.5, 0, 0, 0.5, -0.5, 0,  0, -0.09999999999999998, -0.1666666666666667, 0, 
00410     //
00411     0, 1.0, 2.0, 0, 0, 1.0, -2.0, 0, 0, -3.0, 
00412     -2.0, 0, 0, -3.0, 2.0, 0, 0, 1.0, 0, 0, 
00413     0, -1.0, -2.0, 0, 0, -3.0, 0,  0, 0, -1.0, 
00414     2.0, 0, 0, -1.0, 0, 0, 0, 0.2, -0.6666666666666665, 0,
00415     //
00416     0, -0.5, -0.5, 0, 0, -0.5, 1.5, 0, 0, 1.5, 
00417     1.5, 0, 0, 1.5, -0.5, 0, 0, -0.5, 0.5, 0, 
00418     0, 0.5, 1.5, 0, 0, 1.5, 0.5, 0, 0, 0.5, 
00419     -0.5, 0, 0, 0.5, 0.5, 0, 0, -0.09999999999999998, 0.8333333333333333, 0
00420   };
00421   //D4: Correct multiset of second order partials in (F,P,Dk) format. D4 cardinality = 5 for 2D
00422   double basisD4[] = {
00423     0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 
00424     0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 
00425     0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0,
00426     0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 
00427     0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 
00428     //
00429     0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0, 
00430     0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0, 
00431     0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0, 
00432     0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0, 
00433     0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0,
00434     //
00435     0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 
00436     0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 
00437     0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 
00438     0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 
00439     0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 
00440     //
00441     0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0, 
00442     0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0, 
00443     0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0, 
00444     0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0, 
00445     0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0,
00446     //
00447     0, 0, 4.0, 0, 0, 0, 0, 4.0, 0, 0, 
00448     0, 0, 4.0, 0, 0, 0, 0, 4.0, 0, 0, 
00449     0, 0, 4.0, 0, 0, 0, 0, 4.0, 0, 0, 
00450     0, 0, 4.0, 0, 0, 0, 0, 4.0, 0, 0, 
00451     0, 0, 4.0, 0, 0, 0, 0, 4.0, 0, 0,
00452     //
00453     0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0, 
00454     0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0, 
00455     0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0, 
00456     0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0,
00457     0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0,
00458     //
00459     0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 
00460     0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 
00461     0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 
00462     0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 
00463     0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0,
00464     //
00465     0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0, 
00466     0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0, 
00467     0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0, 
00468     0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0, 
00469     0, 0, -2.0, 0, 0, 0, 0, -2.0, 0, 0,
00470     //
00471     0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 
00472     0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 
00473     0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 
00474     0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0, 
00475     0, 0, 1.0, 0, 0, 0, 0, 1.0, 0, 0
00476 };
00477   
00478   try{
00479         
00480     // Dimensions for the output arrays:
00481     int numFields = quadBasis.getCardinality();
00482     int numPoints = quadNodes.dimension(0);
00483     int spaceDim  = quadBasis.getBaseCellTopology().getDimension();
00484     
00485     // Generic array for values, grads, curls, etc. that will be properly sized before each call
00486     FieldContainer<double> vals;
00487     
00488     // Check VALUE of basis functions: resize vals to rank-2 container:
00489     vals.resize(numFields, numPoints);
00490     quadBasis.getValues(vals, quadNodes, OPERATOR_VALUE);
00491     for (int i = 0; i < numFields; i++) {
00492       for (int j = 0; j < numPoints; j++) {
00493         
00494         // Compute offset for (F,P) container
00495         int l =  j + i * numPoints;
00496         if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
00497           errorFlag++;
00498           *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00499           
00500           // Output the multi-index of the value where the error is:
00501           *outStream << " At multi-index { ";
00502           *outStream << i << " ";*outStream << j << " ";
00503           *outStream << "}  computed value: " << vals(i,j)
00504                      << " but reference value: " << basisValues[l] << "\n";
00505         }
00506       }
00507     }
00508 
00509     // Check GRAD of basis function: resize vals to rank-3 container
00510     vals.resize(numFields, numPoints, spaceDim);
00511     quadBasis.getValues(vals, quadNodes, OPERATOR_GRAD);
00512     for (int i = 0; i < numFields; i++) {
00513       for (int j = 0; j < numPoints; j++) {
00514         for (int k = 0; k < spaceDim; k++) {
00515           
00516           // basisGrads is (F,P,D), compute offset:
00517           int l = k + j * spaceDim + i * spaceDim * numPoints;
00518            if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
00519              errorFlag++;
00520              *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00521 
00522              // Output the multi-index of the value where the error is:
00523              *outStream << " At multi-index { ";
00524              *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
00525              *outStream << "}  computed grad component: " << vals(i,j,k)
00526                << " but reference grad component: " << basisGrads[l] << "\n";
00527             }
00528          }
00529       }
00530     }
00531 
00532     
00533     // Check D1 of basis function (do not resize vals because it has the correct size: D1 = GRAD)
00534     quadBasis.getValues(vals, quadNodes, OPERATOR_D1);
00535     for (int i = 0; i < numFields; i++) {
00536       for (int j = 0; j < numPoints; j++) {
00537         for (int k = 0; k < spaceDim; k++) {
00538           
00539           // basisGrads is (F,P,D), compute offset:
00540           int l = k + j * spaceDim + i * spaceDim * numPoints;
00541            if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
00542              errorFlag++;
00543              *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00544 
00545              // Output the multi-index of the value where the error is:
00546              *outStream << " At multi-index { ";
00547              *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
00548              *outStream << "}  computed D1 component: " << vals(i,j,k)
00549                << " but reference D1 component: " << basisGrads[l] << "\n";
00550             }
00551          }
00552       }
00553     }
00554 
00555 
00556     // Check CURL of basis function: resize vals just for illustration! 
00557     vals.resize(numFields, numPoints, spaceDim);
00558     quadBasis.getValues(vals, quadNodes, OPERATOR_CURL);
00559     for (int i = 0; i < numFields; i++) {
00560       for (int j = 0; j < numPoints; j++) {
00561         // We will use "rotated" basisGrads to check CURL: get offsets to extract (u_y, -u_x)
00562         int curl_0 = 1 + j * spaceDim + i * spaceDim * numPoints;               // position of y-derivative
00563         int curl_1 = 0 + j * spaceDim + i * spaceDim * numPoints;               // position of x-derivative
00564         
00565         double curl_value_0 = basisGrads[curl_0];
00566         double curl_value_1 =-basisGrads[curl_1];
00567         if (std::abs(vals(i,j,0) - curl_value_0) > INTREPID_TOL) {
00568           errorFlag++;
00569           *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00570           // Output the multi-index of the value where the error is:
00571           *outStream << " At multi-index { ";
00572           *outStream << i << " ";*outStream << j << " ";*outStream << 0 << " ";
00573           *outStream << "}  computed curl component: " << vals(i,j,0)
00574             << " but reference curl component: " << curl_value_0 << "\n";
00575         }
00576         if (std::abs(vals(i,j,1) - curl_value_1) > INTREPID_TOL) {
00577           errorFlag++;
00578           *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00579           // Output the multi-index of the value where the error is:
00580           *outStream << " At multi-index { ";
00581           *outStream << i << " ";*outStream << j << " ";*outStream << 1 << " ";
00582           *outStream << "}  computed curl component: " << vals(i,j,1)
00583             << " but reference curl component: " << curl_value_1 << "\n";
00584         }
00585       }
00586     }
00587     
00588     // Check D2 of basis function
00589     int D2cardinality = Intrepid::getDkCardinality(OPERATOR_D2, spaceDim);
00590     vals.resize(numFields, numPoints, D2cardinality);    
00591     quadBasis.getValues(vals, quadNodes, OPERATOR_D2);
00592     for (int i = 0; i < numFields; i++) {
00593       for (int j = 0; j < numPoints; j++) {
00594         for (int k = 0; k < D2cardinality; k++) {
00595           
00596           // basisD2 is (F,P,Dk), compute offset:
00597           int l = k + j * D2cardinality + i * D2cardinality * numPoints;
00598            if (std::abs(vals(i,j,k) - basisD2[l]) > INTREPID_TOL) {
00599              errorFlag++;
00600              *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00601 
00602              // Output the multi-index of the value where the error is:
00603              *outStream << " At multi-index { ";
00604              *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
00605              *outStream << "}  computed D2 component: " << vals(i,j,k)
00606                << " but reference D2 component: " << basisD2[l] << "\n";
00607             }
00608          }
00609       }
00610     }
00611 
00612     
00613     // Check D3 of basis function
00614     int D3cardinality = Intrepid::getDkCardinality(OPERATOR_D3, spaceDim);
00615     vals.resize(numFields, numPoints, D3cardinality);    
00616     quadBasis.getValues(vals, quadNodes, OPERATOR_D3);
00617     for (int i = 0; i < numFields; i++) {
00618       for (int j = 0; j < numPoints; j++) {
00619         for (int k = 0; k < D3cardinality; k++) {
00620           
00621           // basisD3 is (F,P,Dk), compute offset:
00622           int l = k + j * D3cardinality + i * D3cardinality * numPoints;
00623           if (std::abs(vals(i,j,k) - basisD3[l]) > INTREPID_TOL) {
00624             errorFlag++;
00625             *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00626             
00627             // Output the multi-index of the value where the error is:
00628             *outStream << " At multi-index { ";
00629             *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
00630             *outStream << "}  computed D3 component: " << vals(i,j,k)
00631               << " but reference D3 component: " << basisD2[l] << "\n";
00632           }
00633         }
00634       }
00635     }
00636     
00637     // Check D4 of basis function
00638     int D4cardinality = Intrepid::getDkCardinality(OPERATOR_D4, spaceDim);
00639     vals.resize(numFields, numPoints, D4cardinality);    
00640     quadBasis.getValues(vals, quadNodes, OPERATOR_D4);
00641     for (int i = 0; i < numFields; i++) {
00642       for (int j = 0; j < numPoints; j++) {
00643         for (int k = 0; k < D4cardinality; k++) {
00644           
00645           // basisD4 is (F,P,Dk), compute offset:
00646           int l = k + j * D4cardinality + i * D4cardinality * numPoints;
00647           if (std::abs(vals(i,j,k) - basisD4[l]) > INTREPID_TOL) {
00648             errorFlag++;
00649             *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00650             
00651             // Output the multi-index of the value where the error is:
00652             *outStream << " At multi-index { ";
00653             *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
00654             *outStream << "}  computed D4 component: " << vals(i,j,k)
00655               << " but reference D4 component: " << basisD4[l] << "\n";
00656           }
00657         }
00658       }
00659     }
00660     
00661 
00662     // Check all higher derivatives - must be zero. 
00663     for(EOperator op = OPERATOR_D5; op <= OPERATOR_D6; op++) {
00664       
00665       // The last dimension is the number of kth derivatives and needs to be resized for every Dk
00666       int DkCardin  = Intrepid::getDkCardinality(op, spaceDim);
00667       vals.resize(numFields, numPoints, DkCardin);    
00668 
00669       quadBasis.getValues(vals, quadNodes, op);
00670       for (int i = 0; i < vals.size(); i++) {
00671         if (std::abs(vals[i]) > INTREPID_TOL) {
00672           errorFlag++;
00673           *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00674           
00675           // Get the multi-index of the value where the error is and the operator order
00676           std::vector<int> myIndex;
00677           vals.getMultiIndex(myIndex,i);
00678           int ord = Intrepid::getOperatorOrder(op);
00679           *outStream << " At multi-index { ";
00680           for(int j = 0; j < vals.rank(); j++) {
00681             *outStream << myIndex[j] << " ";
00682           }
00683           *outStream << "}  computed D"<< ord <<" component: " << vals[i] 
00684             << " but reference D" << ord << " component:  0 \n";
00685         }
00686       }
00687 }
00688   }
00689   
00690   // Catch unexpected errors
00691   catch (std::logic_error err) {
00692     *outStream << err.what() << "\n\n";
00693     errorFlag = -1000;
00694   };
00695   
00696   if (errorFlag != 0)
00697     std::cout << "End Result: TEST FAILED\n";
00698   else
00699     std::cout << "End Result: TEST PASSED\n";
00700   
00701   // reset format state of std::cout
00702   std::cout.copyfmt(oldFormatState);
00703   
00704   return errorFlag;
00705 }