Intrepid
http://trilinos.sandia.gov/packages/docs/r10.12/packages/intrepid/test/Discretization/Basis/HGRAD_QUAD_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 
00048 #include "Intrepid_FieldContainer.hpp"
00049 #include "Intrepid_HGRAD_QUAD_C2_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_QUAD_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 (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_QUAD_C2_FEM<double, FieldContainer<double> > quadBasis;
00110   int errorFlag = 0;
00111 
00112  // Initialize throw counter for exception testing
00113   int nException     = 0;
00114   int throwCounter   = 0;
00115 
00116   // Array with the 4 vertices, 4 edge midpoints, center of the reference QUAD and a random point.  
00117   FieldContainer<double> quadNodes(10, 2);
00118   quadNodes(0,0) = -1.0;  quadNodes(0,1) = -1.0;
00119   quadNodes(1,0) =  1.0;  quadNodes(1,1) = -1.0;
00120   quadNodes(2,0) =  1.0;  quadNodes(2,1) =  1.0;
00121   quadNodes(3,0) = -1.0;  quadNodes(3,1) =  1.0;
00122   // edge midpoints
00123   quadNodes(4,0) =  0.0;  quadNodes(4,1) = -1.0;
00124   quadNodes(5,0) =  1.0;  quadNodes(5,1) =  0.0;
00125   quadNodes(6,0) =  0.0;  quadNodes(6,1) =  1.0;
00126   quadNodes(7,0) = -1.0;  quadNodes(7,1) =  0.0;
00127   // center & random point
00128   quadNodes(8,0) =  0.0;  quadNodes(8,1) =  0.0;
00129   quadNodes(9,0) =1./3.;  quadNodes(9,1) =-3./5.;
00130   
00131   
00132   // Generic array for the output values; needs to be properly resized depending on the operator type
00133   FieldContainer<double> vals;
00134 
00135   try{
00136     // exception #1: DIV cannot be applied to scalar functions
00137     // resize vals to rank-2 container with dimensions (num. points, num. basis functions)
00138     vals.resize(quadBasis.getCardinality(), quadNodes.dimension(0));
00139     INTREPID_TEST_COMMAND( quadBasis.getValues(vals, quadNodes, OPERATOR_DIV), throwCounter, nException );
00140         
00141     // Exceptions 2-6: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and 
00142     // getDofTag() to access invalid array elements thereby causing bounds check exception
00143     // exception #2
00144     INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(3,0,0), throwCounter, nException );
00145     // exception #3
00146     INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(1,1,1), throwCounter, nException );
00147     // exception #4
00148     INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(0,4,0), throwCounter, nException );
00149     // exception #5
00150     INTREPID_TEST_COMMAND( quadBasis.getDofTag(10), throwCounter, nException );
00151     // exception #6
00152     INTREPID_TEST_COMMAND( quadBasis.getDofTag(-1), throwCounter, nException );
00153     
00154 #ifdef HAVE_INTREPID_DEBUG
00155     // Exceptions 7- test exception handling with incorrectly dimensioned input/output arrays
00156     // exception #7: input points array must be of rank-2
00157     FieldContainer<double> badPoints1(4, 5, 3);
00158     INTREPID_TEST_COMMAND( quadBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
00159     
00160     // exception #8 dimension 1 in the input point array must equal space dimension of the cell
00161     FieldContainer<double> badPoints2(4, quadBasis.getBaseCellTopology().getDimension() + 1);
00162     INTREPID_TEST_COMMAND( quadBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
00163     
00164     // exception #9 output values must be of rank-2 for OPERATOR_VALUE
00165     FieldContainer<double> badVals1(4, 3, 1);
00166     INTREPID_TEST_COMMAND( quadBasis.getValues(badVals1, quadNodes, OPERATOR_VALUE), throwCounter, nException );
00167 
00168     // exception #10 output values must be of rank-3 for OPERATOR_GRAD
00169     FieldContainer<double> badVals2(4, 3);
00170     INTREPID_TEST_COMMAND( quadBasis.getValues(badVals2, quadNodes, OPERATOR_GRAD), throwCounter, nException );
00171     
00172     // exception #11 output values must be of rank-3 for OPERATOR_CURL
00173     INTREPID_TEST_COMMAND( quadBasis.getValues(badVals2, quadNodes, OPERATOR_CURL), throwCounter, nException );
00174     
00175     // exception #12 output values must be of rank-3 for OPERATOR_D2
00176     INTREPID_TEST_COMMAND( quadBasis.getValues(badVals2, quadNodes, OPERATOR_D2), throwCounter, nException );
00177     
00178     // exception #13 incorrect 0th dimension of output array (must equal number of basis functions)
00179     FieldContainer<double> badVals3(quadBasis.getCardinality() + 1, quadNodes.dimension(0));
00180     INTREPID_TEST_COMMAND( quadBasis.getValues(badVals3, quadNodes, OPERATOR_VALUE), throwCounter, nException );
00181     
00182     // exception #14 incorrect 1st dimension of output array (must equal number of points in quadNodes)
00183     FieldContainer<double> badVals4(quadBasis.getCardinality(), quadNodes.dimension(0) + 1);
00184     INTREPID_TEST_COMMAND( quadBasis.getValues(badVals4, quadNodes, OPERATOR_VALUE), throwCounter, nException );
00185     
00186     // exception #15: incorrect 2nd dimension of output array (must equal the space dimension)
00187     FieldContainer<double> badVals5(quadBasis.getCardinality(), quadNodes.dimension(0), quadBasis.getBaseCellTopology().getDimension() + 1);
00188     INTREPID_TEST_COMMAND( quadBasis.getValues(badVals5, quadNodes, OPERATOR_GRAD), throwCounter, nException );
00189     
00190     // exception #16: incorrect 2nd dimension of output array (must equal D2 cardinality in 2D)
00191     FieldContainer<double> badVals6(quadBasis.getCardinality(), quadNodes.dimension(0), 40);
00192     INTREPID_TEST_COMMAND( quadBasis.getValues(badVals6, quadNodes, OPERATOR_D2), throwCounter, nException );
00193     
00194     // exception #17: incorrect 2nd dimension of output array (must equal D3 cardinality in 2D)
00195     INTREPID_TEST_COMMAND( quadBasis.getValues(badVals6, quadNodes, OPERATOR_D3), throwCounter, nException );
00196 #endif
00197 
00198   }
00199   catch (std::logic_error err) {
00200     *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
00201     *outStream << err.what() << '\n';
00202     *outStream << "-------------------------------------------------------------------------------" << "\n\n";
00203     errorFlag = -1000;
00204   };
00205   
00206   // Check if number of thrown exceptions matches the one we expect 
00207   if (throwCounter != nException) {
00208     errorFlag++;
00209     *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00210   }
00211   
00212   *outStream \
00213     << "\n"
00214     << "===============================================================================\n"\
00215     << "| TEST 2: correctness of tag to enum and enum to tag lookups                  |\n"\
00216     << "===============================================================================\n";
00217   
00218   try{
00219     std::vector<std::vector<int> > allTags = quadBasis.getAllDofTags();
00220     
00221     // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again
00222     for (unsigned i = 0; i < allTags.size(); i++) {
00223       int bfOrd  = quadBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
00224       
00225       std::vector<int> myTag = quadBasis.getDofTag(bfOrd);
00226        if( !( (myTag[0] == allTags[i][0]) &&
00227               (myTag[1] == allTags[i][1]) &&
00228               (myTag[2] == allTags[i][2]) &&
00229               (myTag[3] == allTags[i][3]) ) ) {
00230         errorFlag++;
00231         *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00232         *outStream << " getDofOrdinal( {" 
00233           << allTags[i][0] << ", " 
00234           << allTags[i][1] << ", " 
00235           << allTags[i][2] << ", " 
00236           << allTags[i][3] << "}) = " << bfOrd <<" but \n";   
00237         *outStream << " getDofTag(" << bfOrd << ") = { "
00238           << myTag[0] << ", " 
00239           << myTag[1] << ", " 
00240           << myTag[2] << ", " 
00241           << myTag[3] << "}\n";        
00242       }
00243     }
00244     
00245     // Now do the same but loop over basis functions
00246     for( int bfOrd = 0; bfOrd < quadBasis.getCardinality(); bfOrd++) {
00247       std::vector<int> myTag  = quadBasis.getDofTag(bfOrd);
00248       int myBfOrd = quadBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
00249       if( bfOrd != myBfOrd) {
00250         errorFlag++;
00251         *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00252         *outStream << " getDofTag(" << bfOrd << ") = { "
00253           << myTag[0] << ", " 
00254           << myTag[1] << ", " 
00255           << myTag[2] << ", " 
00256           << myTag[3] << "} but getDofOrdinal({" 
00257           << myTag[0] << ", " 
00258           << myTag[1] << ", " 
00259           << myTag[2] << ", " 
00260           << myTag[3] << "} ) = " << myBfOrd << "\n";
00261       }
00262     }
00263   }
00264   catch (std::logic_error err){
00265     *outStream << err.what() << "\n\n";
00266     errorFlag = -1000;
00267   };
00268   
00269   *outStream \
00270     << "\n"
00271     << "===============================================================================\n"\
00272     << "| TEST 3: correctness of basis function values                                |\n"\
00273     << "===============================================================================\n";
00274   
00275   outStream -> precision(20);
00276   
00277   // VALUE: Correct basis values in (F,P) format:
00278   double basisValues[] = {
00279     1.000000000000000, 0, 0, 0, 0, 0, 0, 0, 0, -0.05333333333333334, 0, \
00280     1.000000000000000, 0, 0, 0, 0, 0, 0, 0, 0.1066666666666667, 0, 0, \
00281     1.000000000000000, 0, 0, 0, 0, 0, 0, -0.02666666666666666, 0, 0, 0, \
00282     1.000000000000000, 0, 0, 0, 0, 0, 0.01333333333333333, 0, 0, 0, 0, \
00283     1.000000000000000, 0, 0, 0, 0, 0.4266666666666667, 0, 0, 0, 0, 0, \
00284     1.000000000000000, 0, 0, 0, 0.1422222222222222, 0, 0, 0, 0, 0, 0, \
00285     1.000000000000000, 0, 0, -0.1066666666666667, 0, 0, 0, 0, 0, 0, 0, \
00286     1.000000000000000, 0, -0.07111111111111112, 0, 0, 0, 0, 0, 0, 0, 0, \
00287     1.000000000000000, 0.5688888888888890};
00288   
00289   // GRAD and D1: Correct gradients and D1 in (F,P,D) format: each row contains 6x2 values of
00290   double basisGrads[] = {
00291     -1.500000000000000, -1.500000000000000, 0.5000000000000000, 0, 0, 0, \
00292     0, 0.5000000000000000, -0.5000000000000000, 0, 0, 0, 0, 0, 0, \
00293     -0.5000000000000000, 0, 0, -0.08000000000000002, 0.1222222222222222, \
00294     -0.5000000000000000, 0, 1.500000000000000, -1.500000000000000, 0, \
00295     0.5000000000000000, 0, 0, 0.5000000000000000, 0, 0, \
00296     -0.5000000000000000, 0, 0, 0, 0, 0, 0, 0.3999999999999999, \
00297     -0.2444444444444444, 0, 0, 0, -0.5000000000000000, 1.500000000000000, \
00298     1.500000000000000, -0.5000000000000000, 0, 0, 0, 0, \
00299     0.5000000000000000, 0.5000000000000000, 0, 0, 0, 0, 0, \
00300     -0.09999999999999998, -0.02222222222222221, 0, -0.5000000000000000, \
00301     0, 0, 0.5000000000000000, 0, -1.500000000000000, 1.500000000000000, \
00302     0, 0, 0, 0, -0.5000000000000000, 0, 0, 0.5000000000000000, 0, 0, \
00303     0.02000000000000000, 0.01111111111111112, 2.000000000000000, 0, \
00304     -2.000000000000000, 0, 0, 0, 0, 0, 0, -1.500000000000000, 0, 0, 0, \
00305     0.5000000000000000, 0, 0, 0, -0.5000000000000000, \
00306     -0.3199999999999999, -0.9777777777777779, 0, 0, 0, 2.000000000000000, \
00307     0, -2.000000000000000, 0, 0, 0, 0, 1.500000000000000, 0, 0, 0, \
00308     -0.5000000000000000, 0, 0.5000000000000000, 0, 0.5333333333333334, \
00309     0.2666666666666666, 0, 0, 0, 0, -2.000000000000000, 0, \
00310     2.000000000000000, 0, 0, -0.5000000000000000, 0, 0, 0, \
00311     1.500000000000000, 0, 0, 0, 0.5000000000000000, 0.07999999999999997, \
00312     -0.08888888888888888, 0, 2.000000000000000, 0, 0, 0, 0, 0, \
00313     -2.000000000000000, 0, 0, 0.5000000000000000, 0, 0, 0, \
00314     -1.500000000000000, 0, -0.5000000000000000, 0, -0.1066666666666667, \
00315     -0.1333333333333333, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2.000000000000000, \
00316     -2.000000000000000, 0, 0, -2.000000000000000, 2.000000000000000, 0, \
00317     0, 0, -0.4266666666666667, 1.066666666666667};
00318   
00319   // D2: Correct multiset of second order partials in (F,P,Dk) format. D2 cardinality = 3 for 2D 
00320   double basisD2[] = {
00321     1.000000000000000, 2.250000000000000, 1.000000000000000, \
00322     1.000000000000000, -0.7500000000000000, 0, 0, 0.2500000000000000, 0, \
00323     0, -0.7500000000000000, 1.000000000000000, 1.000000000000000, \
00324     0.7500000000000000, 0, 0, -0.2500000000000000, 0, 0, \
00325     -0.2500000000000000, 0, 0, 0.7500000000000000, 1.000000000000000, 0, \
00326     0.2500000000000000, 0, 0.4800000000000000, 0.1833333333333334, \
00327     -0.1111111111111111, 1.000000000000000, 0.7500000000000000, 0, \
00328     1.000000000000000, -2.250000000000000, 1.000000000000000, 0, \
00329     0.7500000000000000, 1.000000000000000, 0, -0.2500000000000000, 0, \
00330     1.000000000000000, -0.7500000000000000, 0, 0, -0.7500000000000000, \
00331     1.000000000000000, 0, 0.2500000000000000, 0, 0, 0.2500000000000000, \
00332     0, 0, -0.2500000000000000, 0, 0.4800000000000000, \
00333     -0.9166666666666666, 0.2222222222222222, 0, 0.2500000000000000, 0, 0, \
00334     -0.7500000000000000, 1.000000000000000, 1.000000000000000, \
00335     2.250000000000000, 1.000000000000000, 1.000000000000000, \
00336     -0.7500000000000000, 0, 0, -0.2500000000000000, 0, 0, \
00337     0.7500000000000000, 1.000000000000000, 1.000000000000000, \
00338     0.7500000000000000, 0, 0, -0.2500000000000000, 0, 0, \
00339     0.2500000000000000, 0, -0.1200000000000000, -0.08333333333333331, \
00340     0.2222222222222222, 0, 0.7500000000000000, 1.000000000000000, 0, \
00341     -0.2500000000000000, 0, 1.000000000000000, 0.7500000000000000, 0, \
00342     1.000000000000000, -2.250000000000000, 1.000000000000000, 0, \
00343     0.2500000000000000, 0, 0, 0.2500000000000000, 0, 1.000000000000000, \
00344     -0.7500000000000000, 0, 0, -0.7500000000000000, 1.000000000000000, 0, \
00345     -0.2500000000000000, 0, -0.1200000000000000, 0.01666666666666666, \
00346     -0.1111111111111111, -2.000000000000000, -3.000000000000000, 0, \
00347     -2.000000000000000, 3.000000000000000, 0, 0, -1.000000000000000, 0, \
00348     0, 1.000000000000000, 0, -2.000000000000000, 0, 1.000000000000000, 0, \
00349     1.000000000000000, 0, 0, 0, 1.000000000000000, 0, -1.000000000000000, \
00350     0, 0, 0, 1.000000000000000, -0.9600000000000000, 0.7333333333333332, \
00351     0.8888888888888890, 0, -1.000000000000000, 0, 0, 3.000000000000000, \
00352     -2.000000000000000, 0, -3.000000000000000, -2.000000000000000, 0, \
00353     1.000000000000000, 0, 0, 1.000000000000000, 0, 1.000000000000000, 0, \
00354     -2.000000000000000, 0, -1.000000000000000, 0, 1.000000000000000, 0, \
00355     0, 1.000000000000000, 0, 0, 0.6400000000000001, 1.000000000000000, \
00356     -0.4444444444444444, 0, -1.000000000000000, 0, 0, 1.000000000000000, \
00357     0, -2.000000000000000, -3.000000000000000, 0, -2.000000000000000, \
00358     3.000000000000000, 0, 0, 0, 1.000000000000000, 0, -1.000000000000000, \
00359     0, -2.000000000000000, 0, 1.000000000000000, 0, 1.000000000000000, 0, \
00360     0, 0, 1.000000000000000, 0.2400000000000000, 0.06666666666666665, \
00361     0.8888888888888890, 0, -3.000000000000000, -2.000000000000000, 0, \
00362     1.000000000000000, 0, 0, -1.000000000000000, 0, 0, 3.000000000000000, \
00363     -2.000000000000000, 0, -1.000000000000000, 0, 1.000000000000000, 0, \
00364     0, 0, 1.000000000000000, 0, 1.000000000000000, 0, -2.000000000000000, \
00365     1.000000000000000, 0, 0, 0.6400000000000001, -0.2000000000000001, \
00366     0.2222222222222222, 0, 4.000000000000000, 0, 0, -4.000000000000000, \
00367     0, 0, 4.000000000000000, 0, 0, -4.000000000000000, 0, 0, 0, \
00368     -2.000000000000000, -2.000000000000000, 0, 0, 0, 0, \
00369     -2.000000000000000, -2.000000000000000, 0, 0, -2.000000000000000, 0, \
00370     -2.000000000000000, -1.280000000000000, -0.7999999999999998, \
00371     -1.777777777777778};
00372   
00373   //D3: Correct multiset of second order partials in (F,P,Dk) format. D3 cardinality = 4 for 2D
00374   double basisD3[] = {
00375     0, -1.500000000000000, -1.500000000000000, 0, 0, -1.500000000000000, \
00376     0.5000000000000000, 0, 0, 0.5000000000000000, 0.5000000000000000, 0, \
00377     0, 0.5000000000000000, -1.500000000000000, 0, 0, -1.500000000000000, \
00378     -0.5000000000000000, 0, 0, -0.5000000000000000, 0.5000000000000000, \
00379     0, 0, 0.5000000000000000, -0.5000000000000000, 0, 0, \
00380     -0.5000000000000000, -1.500000000000000, 0, 0, -0.5000000000000000, \
00381     -0.5000000000000000, 0, 0, -1.100000000000000, -0.1666666666666667, \
00382     0, 0, -1.500000000000000, -0.5000000000000000, 0, 0, \
00383     -1.500000000000000, 1.500000000000000, 0, 0, 0.5000000000000000, \
00384     1.500000000000000, 0, 0, 0.5000000000000000, -0.5000000000000000, 0, \
00385     0, -1.500000000000000, 0.5000000000000000, 0, 0, -0.5000000000000000, \
00386     1.500000000000000, 0, 0, 0.5000000000000000, 0.5000000000000000, 0, \
00387     0, -0.5000000000000000, -0.5000000000000000, 0, 0, \
00388     -0.5000000000000000, 0.5000000000000000, 0, 0, -1.100000000000000, \
00389     0.8333333333333333, 0, 0, -0.5000000000000000, -0.5000000000000000, \
00390     0, 0, -0.5000000000000000, 1.500000000000000, 0, 0, \
00391     1.500000000000000, 1.500000000000000, 0, 0, 1.500000000000000, \
00392     -0.5000000000000000, 0, 0, -0.5000000000000000, 0.5000000000000000, \
00393     0, 0, 0.5000000000000000, 1.500000000000000, 0, 0, 1.500000000000000, \
00394     0.5000000000000000, 0, 0, 0.5000000000000000, -0.5000000000000000, 0, \
00395     0, 0.5000000000000000, 0.5000000000000000, 0, 0, \
00396     -0.09999999999999998, 0.8333333333333333, 0, 0, -0.5000000000000000, \
00397     -1.500000000000000, 0, 0, -0.5000000000000000, 0.5000000000000000, 0, \
00398     0, 1.500000000000000, 0.5000000000000000, 0, 0, 1.500000000000000, \
00399     -1.500000000000000, 0, 0, -0.5000000000000000, -0.5000000000000000, \
00400     0, 0, 0.5000000000000000, 0.5000000000000000, 0, 0, \
00401     1.500000000000000, -0.5000000000000000, 0, 0, 0.5000000000000000, \
00402     -1.500000000000000, 0, 0, 0.5000000000000000, -0.5000000000000000, 0, \
00403     0, -0.09999999999999998, -0.1666666666666667, 0, 0, \
00404     3.000000000000000, 2.000000000000000, 0, 0, 3.000000000000000, \
00405     -2.000000000000000, 0, 0, -1.000000000000000, -2.000000000000000, 0, \
00406     0, -1.000000000000000, 2.000000000000000, 0, 0, 3.000000000000000, 0, \
00407     0, 0, 1.000000000000000, -2.000000000000000, 0, 0, \
00408     -1.000000000000000, 0, 0, 0, 1.000000000000000, 2.000000000000000, 0, \
00409     0, 1.000000000000000, 0, 0, 0, 2.200000000000000, \
00410     -0.6666666666666665, 0, 0, 2.000000000000000, 1.000000000000000, 0, \
00411     0, 2.000000000000000, -3.000000000000000, 0, 0, -2.000000000000000, \
00412     -3.000000000000000, 0, 0, -2.000000000000000, 1.000000000000000, 0, \
00413     0, 2.000000000000000, -1.000000000000000, 0, 0, 0, \
00414     -3.000000000000000, 0, 0, -2.000000000000000, -1.000000000000000, 0, \
00415     0, 0, 1.000000000000000, 0, 0, 0, -1.000000000000000, 0, 0, \
00416     1.200000000000000, -1.666666666666667, 0, 0, 1.000000000000000, \
00417     2.000000000000000, 0, 0, 1.000000000000000, -2.000000000000000, 0, 0, \
00418     -3.000000000000000, -2.000000000000000, 0, 0, -3.000000000000000, \
00419     2.000000000000000, 0, 0, 1.000000000000000, 0, 0, 0, \
00420     -1.000000000000000, -2.000000000000000, 0, 0, -3.000000000000000, 0, \
00421     0, 0, -1.000000000000000, 2.000000000000000, 0, 0, \
00422     -1.000000000000000, 0, 0, 0, 0.2000000000000000, -0.6666666666666665, \
00423     0, 0, 2.000000000000000, 3.000000000000000, 0, 0, 2.000000000000000, \
00424     -1.000000000000000, 0, 0, -2.000000000000000, -1.000000000000000, 0, \
00425     0, -2.000000000000000, 3.000000000000000, 0, 0, 2.000000000000000, \
00426     1.000000000000000, 0, 0, 0, -1.000000000000000, 0, 0, \
00427     -2.000000000000000, 1.000000000000000, 0, 0, 0, 3.000000000000000, 0, \
00428     0, 0, 1.000000000000000, 0, 0, 1.200000000000000, 0.3333333333333334, \
00429     0, 0, -4.000000000000000, -4.000000000000000, 0, 0, \
00430     -4.000000000000000, 4.000000000000000, 0, 0, 4.000000000000000, \
00431     4.000000000000000, 0, 0, 4.000000000000000, -4.000000000000000, 0, 0, \
00432     -4.000000000000000, 0, 0, 0, 0, 4.000000000000000, 0, 0, \
00433     4.000000000000000, 0, 0, 0, 0, -4.000000000000000, 0, 0, 0, 0, 0, 0, \
00434     -2.400000000000000, 1.333333333333333, 0};
00435   //D4: Correct multiset of second order partials in (F,P,Dk) format. D4 cardinality = 5 for 2D
00436   double basisD4[] = {
00437     0, 0, 1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
00438     1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
00439     1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
00440     1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
00441     1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
00442     1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
00443     1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
00444     1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
00445     1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
00446     1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
00447     1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
00448     1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
00449     1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
00450     1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
00451     1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
00452     1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
00453     1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
00454     1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
00455     1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
00456     1.000000000000000, 0, 0, 0, 0, 1.000000000000000, 0, 0, 0, 0, \
00457     -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
00458     -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
00459     -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
00460     -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
00461     -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
00462     -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
00463     -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
00464     -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
00465     -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
00466     -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
00467     -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
00468     -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
00469     -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
00470     -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
00471     -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
00472     -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
00473     -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
00474     -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
00475     -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
00476     -2.000000000000000, 0, 0, 0, 0, -2.000000000000000, 0, 0, 0, 0, \
00477     4.000000000000000, 0, 0, 0, 0, 4.000000000000000, 0, 0, 0, 0, \
00478     4.000000000000000, 0, 0, 0, 0, 4.000000000000000, 0, 0, 0, 0, \
00479     4.000000000000000, 0, 0, 0, 0, 4.000000000000000, 0, 0, 0, 0, \
00480     4.000000000000000, 0, 0, 0, 0, 4.000000000000000, 0, 0, 0, 0, \
00481     4.000000000000000, 0, 0, 0, 0, 4.000000000000000, 0, 0};
00482   
00483   try{
00484         
00485     // Dimensions for the output arrays:
00486     int numFields = quadBasis.getCardinality();
00487     int numPoints = quadNodes.dimension(0);
00488     int spaceDim  = quadBasis.getBaseCellTopology().getDimension();
00489     
00490     // Generic array for values, grads, curls, etc. that will be properly sized before each call
00491     FieldContainer<double> vals;
00492     
00493     // Check VALUE of basis functions: resize vals to rank-2 container:
00494     vals.resize(numFields, numPoints);
00495     quadBasis.getValues(vals, quadNodes, OPERATOR_VALUE);
00496     for (int i = 0; i < numFields; i++) {
00497       for (int j = 0; j < numPoints; j++) {
00498         
00499         // Compute offset for (F,P) container
00500         int l =  j + i * numPoints;
00501            if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
00502              errorFlag++;
00503              *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00504 
00505              // Output the multi-index of the value where the error is:
00506              *outStream << " At multi-index { ";
00507              *outStream << i << " ";*outStream << j << " ";
00508              *outStream << "}  computed value: " << vals(i,j)
00509                << " but reference value: " << basisValues[l] << "\n";
00510          }
00511       }
00512     }
00513 
00514     // Check GRAD of basis function: resize vals to rank-3 container
00515     vals.resize(numFields, numPoints, spaceDim);
00516     quadBasis.getValues(vals, quadNodes, OPERATOR_GRAD);
00517     for (int i = 0; i < numFields; i++) {
00518       for (int j = 0; j < numPoints; j++) {
00519         for (int k = 0; k < spaceDim; k++) {
00520           
00521           // basisGrads is (F,P,D), compute offset:
00522           int l = k + j * spaceDim + i * spaceDim * numPoints;
00523            if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
00524              errorFlag++;
00525              *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00526 
00527              // Output the multi-index of the value where the error is:
00528              *outStream << " At multi-index { ";
00529              *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
00530              *outStream << "}  computed grad component: " << vals(i,j,k)
00531                << " but reference grad component: " << basisGrads[l] << "\n";
00532             }
00533          }
00534       }
00535     }
00536 
00537     
00538     // Check D1 of basis function (do not resize vals because it has the correct size: D1 = GRAD)
00539     quadBasis.getValues(vals, quadNodes, OPERATOR_D1);
00540     for (int i = 0; i < numFields; i++) {
00541       for (int j = 0; j < numPoints; j++) {
00542         for (int k = 0; k < spaceDim; k++) {
00543           
00544           // basisGrads is (F,P,D), compute offset:
00545           int l = k + j * spaceDim + i * spaceDim * numPoints;
00546            if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
00547              errorFlag++;
00548              *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00549 
00550              // Output the multi-index of the value where the error is:
00551              *outStream << " At multi-index { ";
00552              *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
00553              *outStream << "}  computed D1 component: " << vals(i,j,k)
00554                << " but reference D1 component: " << basisGrads[l] << "\n";
00555             }
00556          }
00557       }
00558     }
00559 
00560 
00561     // Check CURL of basis function: resize vals just for illustration! 
00562     vals.resize(numFields, numPoints, spaceDim);
00563     quadBasis.getValues(vals, quadNodes, OPERATOR_CURL);
00564     for (int i = 0; i < numFields; i++) {
00565       for (int j = 0; j < numPoints; j++) {
00566         // We will use "rotated" basisGrads to check CURL: get offsets to extract (u_y, -u_x)
00567         int curl_0 = 1 + j * spaceDim + i * spaceDim * numPoints;               // position of y-derivative
00568         int curl_1 = 0 + j * spaceDim + i * spaceDim * numPoints;               // position of x-derivative
00569         
00570         double curl_value_0 = basisGrads[curl_0];
00571         double curl_value_1 =-basisGrads[curl_1];
00572         if (std::abs(vals(i,j,0) - curl_value_0) > INTREPID_TOL) {
00573           errorFlag++;
00574           *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00575           // Output the multi-index of the value where the error is:
00576           *outStream << " At multi-index { ";
00577           *outStream << i << " ";*outStream << j << " ";*outStream << 0 << " ";
00578           *outStream << "}  computed curl component: " << vals(i,j,0)
00579             << " but reference curl component: " << curl_value_0 << "\n";
00580         }
00581         if (std::abs(vals(i,j,1) - curl_value_1) > INTREPID_TOL) {
00582           errorFlag++;
00583           *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00584           // Output the multi-index of the value where the error is:
00585           *outStream << " At multi-index { ";
00586           *outStream << i << " ";*outStream << j << " ";*outStream << 1 << " ";
00587           *outStream << "}  computed curl component: " << vals(i,j,1)
00588             << " but reference curl component: " << curl_value_1 << "\n";
00589         }
00590       }
00591     }
00592     
00593     // Check D2 of basis function
00594     int D2cardinality = Intrepid::getDkCardinality(OPERATOR_D2, spaceDim);
00595     vals.resize(numFields, numPoints, D2cardinality);    
00596     quadBasis.getValues(vals, quadNodes, OPERATOR_D2);
00597     for (int i = 0; i < numFields; i++) {
00598       for (int j = 0; j < numPoints; j++) {
00599         for (int k = 0; k < D2cardinality; k++) {
00600           
00601           // basisD2 is (F,P,Dk), compute offset:
00602           int l = k + j * D2cardinality + i * D2cardinality * numPoints;
00603            if (std::abs(vals(i,j,k) - basisD2[l]) > INTREPID_TOL) {
00604              errorFlag++;
00605              *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00606 
00607              // Output the multi-index of the value where the error is:
00608              *outStream << " At multi-index { ";
00609              *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
00610              *outStream << "}  computed D2 component: " << vals(i,j,k)
00611                << " but reference D2 component: " << basisD2[l] << "\n";
00612             }
00613          }
00614       }
00615     }
00616 
00617     
00618     // Check D3 of basis function
00619     int D3cardinality = Intrepid::getDkCardinality(OPERATOR_D3, spaceDim);
00620     vals.resize(numFields, numPoints, D3cardinality);    
00621     quadBasis.getValues(vals, quadNodes, OPERATOR_D3);
00622     for (int i = 0; i < numFields; i++) {
00623       for (int j = 0; j < numPoints; j++) {
00624         for (int k = 0; k < D3cardinality; k++) {
00625           
00626           // basisD3 is (F,P,Dk), compute offset:
00627           int l = k + j * D3cardinality + i * D3cardinality * numPoints;
00628           if (std::abs(vals(i,j,k) - basisD3[l]) > INTREPID_TOL) {
00629             errorFlag++;
00630             *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00631             
00632             // Output the multi-index of the value where the error is:
00633             *outStream << " At multi-index { ";
00634             *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
00635             *outStream << "}  computed D3 component: " << vals(i,j,k)
00636               << " but reference D3 component: " << basisD2[l] << "\n";
00637           }
00638         }
00639       }
00640     }
00641 
00642     
00643     // Check D4 of basis function
00644     int D4cardinality = Intrepid::getDkCardinality(OPERATOR_D4, spaceDim);
00645     vals.resize(numFields, numPoints, D4cardinality);    
00646     quadBasis.getValues(vals, quadNodes, OPERATOR_D4);
00647     for (int i = 0; i < numFields; i++) {
00648       for (int j = 0; j < numPoints; j++) {
00649         for (int k = 0; k < D4cardinality; k++) {
00650           
00651           // basisD4 is (F,P,Dk), compute offset:
00652           int l = k + j * D4cardinality + i * D4cardinality * numPoints;
00653           if (std::abs(vals(i,j,k) - basisD4[l]) > INTREPID_TOL) {
00654             errorFlag++;
00655             *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00656             
00657             // Output the multi-index of the value where the error is:
00658             *outStream << " At multi-index { ";
00659             *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
00660             *outStream << "}  computed D4 component: " << vals(i,j,k)
00661               << " but reference D4 component: " << basisD2[l] << "\n";
00662           }
00663         }
00664       }
00665     }
00666     
00667 
00668     // Check all higher derivatives - must be zero. 
00669     for(EOperator op = OPERATOR_D5; op < OPERATOR_MAX; op++) {
00670       
00671       // The last dimension is the number of kth derivatives and needs to be resized for every Dk
00672       int DkCardin  = Intrepid::getDkCardinality(op, spaceDim);
00673       vals.resize(numFields, numPoints, DkCardin);    
00674 
00675       quadBasis.getValues(vals, quadNodes, op);
00676       for (int i = 0; i < vals.size(); i++) {
00677         if (std::abs(vals[i]) > INTREPID_TOL) {
00678           errorFlag++;
00679           *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00680           
00681           // Get the multi-index of the value where the error is and the operator order
00682           std::vector<int> myIndex;
00683           vals.getMultiIndex(myIndex,i);
00684           int ord = Intrepid::getOperatorOrder(op);
00685           *outStream << " At multi-index { ";
00686           for(int j = 0; j < vals.rank(); j++) {
00687             *outStream << myIndex[j] << " ";
00688           }
00689           *outStream << "}  computed D"<< ord <<" component: " << vals[i] 
00690             << " but reference D" << ord << " component:  0 \n";
00691         }
00692       }
00693     }    
00694   }
00695   // Catch unexpected errors
00696   catch (std::logic_error err) {
00697     *outStream << err.what() << "\n\n";
00698     errorFlag = -1000;
00699   };
00700 
00701 
00702   *outStream \
00703     << "\n"
00704     << "===============================================================================\n"\
00705     << "| TEST 4: correctness of DoF locations                                        |\n"\
00706     << "===============================================================================\n";
00707 
00708   try{
00709     Teuchos::RCP<Basis<double, FieldContainer<double> > > basis =
00710       Teuchos::rcp(new Basis_HGRAD_QUAD_C2_FEM<double, FieldContainer<double> >);
00711     Teuchos::RCP<DofCoordsInterface<FieldContainer<double> > > coord_iface =
00712       Teuchos::rcp_dynamic_cast<DofCoordsInterface<FieldContainer<double> > >(basis);
00713 
00714     FieldContainer<double> cvals;
00715     FieldContainer<double> bvals(basis->getCardinality(), basis->getCardinality());
00716 
00717     // Check exceptions.
00718 #ifdef HAVE_INTREPID_DEBUG
00719     cvals.resize(1,2,3);
00720     INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
00721     cvals.resize(8,2);
00722     INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
00723     cvals.resize(9,1);
00724     INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
00725 #endif
00726     cvals.resize(9,2);
00727     INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException ); nException--;
00728     // Check if number of thrown exceptions matches the one we expect
00729     if (throwCounter != nException) {
00730       errorFlag++;
00731       *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00732     }
00733 
00734     // Check mathematical correctness.
00735     basis->getValues(bvals, cvals, OPERATOR_VALUE);
00736     char buffer[120];
00737     for (int i=0; i<bvals.dimension(0); i++) {
00738       for (int j=0; j<bvals.dimension(1); j++) {
00739         if ((i != j) && (std::abs(bvals(i,j) - 0.0) > INTREPID_TOL)) {
00740           errorFlag++;
00741           sprintf(buffer, "\nValue of basis function %d at (%6.4e, %6.4e) is %6.4e but should be %6.4e!\n", i, cvals(i,0), cvals(i,1), bvals(i,j), 0.0);
00742           *outStream << buffer;
00743         }
00744         else if ((i == j) && (std::abs(bvals(i,j) - 1.0) > INTREPID_TOL)) {
00745           errorFlag++;
00746           sprintf(buffer, "\nValue of basis function %d at (%6.4e, %6.4e) is %6.4e but should be %6.4e!\n", i, cvals(i,0), cvals(i,1), bvals(i,j), 1.0);
00747           *outStream << buffer;
00748         }
00749       }
00750     }
00751 
00752   }
00753   catch (std::logic_error err){
00754     *outStream << err.what() << "\n\n";
00755     errorFlag = -1000;
00756   };
00757   
00758   if (errorFlag != 0)
00759     std::cout << "End Result: TEST FAILED\n";
00760   else
00761     std::cout << "End Result: TEST PASSED\n";
00762   
00763   // reset format state of std::cout
00764   std::cout.copyfmt(oldFormatState);
00765   
00766   return errorFlag;
00767 }