Intrepid
http://trilinos.sandia.gov/packages/docs/r10.10/packages/intrepid/test/Discretization/Basis/HGRAD_HEX_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_HEX_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_HEX_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_HEX_C1_FEM<double, FieldContainer<double> > hexBasis;
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 8 vertices of the reference HEX, its center and 6 face centers
00117   FieldContainer<double> hexNodes(15, 3);
00118   hexNodes(0,0) = -1.0;  hexNodes(0,1) = -1.0;  hexNodes(0,2) = -1.0;
00119   hexNodes(1,0) =  1.0;  hexNodes(1,1) = -1.0;  hexNodes(1,2) = -1.0;
00120   hexNodes(2,0) =  1.0;  hexNodes(2,1) =  1.0;  hexNodes(2,2) = -1.0;
00121   hexNodes(3,0) = -1.0;  hexNodes(3,1) =  1.0;  hexNodes(3,2) = -1.0;
00122   
00123   hexNodes(4,0) = -1.0;  hexNodes(4,1) = -1.0;  hexNodes(4,2) =  1.0;
00124   hexNodes(5,0) =  1.0;  hexNodes(5,1) = -1.0;  hexNodes(5,2) =  1.0;
00125   hexNodes(6,0) =  1.0;  hexNodes(6,1) =  1.0;  hexNodes(6,2) =  1.0;
00126   hexNodes(7,0) = -1.0;  hexNodes(7,1) =  1.0;  hexNodes(7,2) =  1.0;  
00127   
00128   hexNodes(8,0) =  0.0;  hexNodes(8,1) =  0.0;  hexNodes(8,2) =  0.0;
00129   
00130   hexNodes(9,0) =  1.0;  hexNodes(9,1) =  0.0;  hexNodes(9,2) =  0.0;
00131   hexNodes(10,0)= -1.0;  hexNodes(10,1)=  0.0;  hexNodes(10,2)=  0.0;
00132   
00133   hexNodes(11,0)=  0.0;  hexNodes(11,1)=  1.0;  hexNodes(11,2)=  0.0;
00134   hexNodes(12,0)=  0.0;  hexNodes(12,1)= -1.0;  hexNodes(12,2)=  0.0;
00135   
00136   hexNodes(13,0)=  0.0;  hexNodes(13,1)=  0.0;  hexNodes(13,2)=  1.0;
00137   hexNodes(14,0)=  0.0;  hexNodes(14,1)=  0.0;  hexNodes(14,2)= -1.0;
00138 
00139   
00140   // Generic array for the output values; needs to be properly resized depending on the operator type
00141   FieldContainer<double> vals;
00142 
00143   try{
00144     // exception #1: CURL cannot be applied to scalar functions in 3D
00145     // resize vals to rank-3 container with dimensions (num. basis functions, num. points, arbitrary)
00146     vals.resize(hexBasis.getCardinality(), hexNodes.dimension(0), 4 );
00147     INTREPID_TEST_COMMAND( hexBasis.getValues(vals, hexNodes, OPERATOR_CURL), throwCounter, nException );
00148 
00149     // exception #2: DIV cannot be applied to scalar functions in 3D
00150     // resize vals to rank-2 container with dimensions (num. basis functions, num. points)
00151     vals.resize(hexBasis.getCardinality(), hexNodes.dimension(0) );
00152     INTREPID_TEST_COMMAND( hexBasis.getValues(vals, hexNodes, OPERATOR_DIV), throwCounter, nException );
00153         
00154     // Exceptions 3-7: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and 
00155     // getDofTag() to access invalid array elements thereby causing bounds check exception
00156     // exception #3
00157     INTREPID_TEST_COMMAND( hexBasis.getDofOrdinal(3,0,0), throwCounter, nException );
00158     // exception #4
00159     INTREPID_TEST_COMMAND( hexBasis.getDofOrdinal(1,1,1), throwCounter, nException );
00160     // exception #5
00161     INTREPID_TEST_COMMAND( hexBasis.getDofOrdinal(0,4,1), throwCounter, nException );
00162     // exception #6
00163     INTREPID_TEST_COMMAND( hexBasis.getDofTag(8), throwCounter, nException );
00164     // exception #7
00165     INTREPID_TEST_COMMAND( hexBasis.getDofTag(-1), throwCounter, nException );
00166 
00167 #ifdef HAVE_INTREPID_DEBUG 
00168     // Exceptions 8-18 test exception handling with incorrectly dimensioned input/output arrays
00169     // exception #8: input points array must be of rank-2
00170     FieldContainer<double> badPoints1(4, 5, 3);
00171     INTREPID_TEST_COMMAND( hexBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
00172 
00173     // exception #9 dimension 1 in the input point array must equal space dimension of the cell
00174     FieldContainer<double> badPoints2(4, 2);
00175     INTREPID_TEST_COMMAND( hexBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
00176 
00177     // exception #10 output values must be of rank-2 for OPERATOR_VALUE
00178     FieldContainer<double> badVals1(4, 3, 1);
00179     INTREPID_TEST_COMMAND( hexBasis.getValues(badVals1, hexNodes, OPERATOR_VALUE), throwCounter, nException );
00180  
00181     // exception #11 output values must be of rank-3 for OPERATOR_GRAD
00182     FieldContainer<double> badVals2(4, 3);
00183     INTREPID_TEST_COMMAND( hexBasis.getValues(badVals2, hexNodes, OPERATOR_GRAD), throwCounter, nException );
00184     
00185     // exception #12 output values must be of rank-3 for OPERATOR_D1
00186     INTREPID_TEST_COMMAND( hexBasis.getValues(badVals2, hexNodes, OPERATOR_D1), throwCounter, nException );
00187     
00188     // exception #13 output values must be of rank-3 for OPERATOR_D2
00189     INTREPID_TEST_COMMAND( hexBasis.getValues(badVals2, hexNodes, OPERATOR_D2), throwCounter, nException );
00190     
00191     // exception #14 incorrect 0th dimension of output array (must equal number of basis functions)
00192     FieldContainer<double> badVals3(hexBasis.getCardinality() + 1, hexNodes.dimension(0));
00193     INTREPID_TEST_COMMAND( hexBasis.getValues(badVals3, hexNodes, OPERATOR_VALUE), throwCounter, nException );
00194     
00195     // exception #15 incorrect 1st dimension of output array (must equal number of points)
00196     FieldContainer<double> badVals4(hexBasis.getCardinality(), hexNodes.dimension(0) + 1);
00197     INTREPID_TEST_COMMAND( hexBasis.getValues(badVals4, hexNodes, OPERATOR_VALUE), throwCounter, nException );
00198     
00199     // exception #16: incorrect 2nd dimension of output array (must equal the space dimension)
00200     FieldContainer<double> badVals5(hexBasis.getCardinality(), hexNodes.dimension(0), 4);
00201     INTREPID_TEST_COMMAND( hexBasis.getValues(badVals5, hexNodes, OPERATOR_GRAD), throwCounter, nException );
00202     
00203     // exception #17: incorrect 2nd dimension of output array (must equal D2 cardinality in 3D)
00204     FieldContainer<double> badVals6(hexBasis.getCardinality(), hexNodes.dimension(0), 40);
00205     INTREPID_TEST_COMMAND( hexBasis.getValues(badVals6, hexNodes, OPERATOR_D2), throwCounter, nException );
00206     
00207     // exception #18: incorrect 2nd dimension of output array (must equal D3 cardinality in 3D)
00208     FieldContainer<double> badVals7(hexBasis.getCardinality(), hexNodes.dimension(0), 50);
00209     INTREPID_TEST_COMMAND( hexBasis.getValues(badVals7, hexNodes, OPERATOR_D3), throwCounter, nException );
00210 #endif
00211 
00212   }
00213   catch (std::logic_error err) {
00214     *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
00215     *outStream << err.what() << '\n';
00216     *outStream << "-------------------------------------------------------------------------------" << "\n\n";
00217     errorFlag = -1000;
00218   };
00219   
00220   // Check if number of thrown exceptions matches the one we expect 
00221   // Note Teuchos throw number will not pick up exceptions 3-7 and therefore will not match.
00222   if (throwCounter != nException) {
00223     errorFlag++;
00224     *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00225   }
00226   
00227   *outStream \
00228     << "\n"
00229     << "===============================================================================\n"\
00230     << "| TEST 2: correctness of tag to enum and enum to tag lookups                  |\n"\
00231     << "===============================================================================\n";
00232   
00233   try{
00234     std::vector<std::vector<int> > allTags = hexBasis.getAllDofTags();
00235     
00236     // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again
00237     for (unsigned i = 0; i < allTags.size(); i++) {
00238       int bfOrd  = hexBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
00239       
00240       std::vector<int> myTag = hexBasis.getDofTag(bfOrd);
00241        if( !( (myTag[0] == allTags[i][0]) &&
00242               (myTag[1] == allTags[i][1]) &&
00243               (myTag[2] == allTags[i][2]) &&
00244               (myTag[3] == allTags[i][3]) ) ) {
00245         errorFlag++;
00246         *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00247         *outStream << " getDofOrdinal( {" 
00248           << allTags[i][0] << ", " 
00249           << allTags[i][1] << ", " 
00250           << allTags[i][2] << ", " 
00251           << allTags[i][3] << "}) = " << bfOrd <<" but \n";   
00252         *outStream << " getDofTag(" << bfOrd << ") = { "
00253           << myTag[0] << ", " 
00254           << myTag[1] << ", " 
00255           << myTag[2] << ", " 
00256           << myTag[3] << "}\n";        
00257       }
00258     }
00259     
00260     // Now do the same but loop over basis functions
00261     for( int bfOrd = 0; bfOrd < hexBasis.getCardinality(); bfOrd++) {
00262       std::vector<int> myTag  = hexBasis.getDofTag(bfOrd);
00263       int myBfOrd = hexBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
00264       if( bfOrd != myBfOrd) {
00265         errorFlag++;
00266         *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00267         *outStream << " getDofTag(" << bfOrd << ") = { "
00268           << myTag[0] << ", " 
00269           << myTag[1] << ", " 
00270           << myTag[2] << ", " 
00271           << myTag[3] << "} but getDofOrdinal({" 
00272           << myTag[0] << ", " 
00273           << myTag[1] << ", " 
00274           << myTag[2] << ", " 
00275           << myTag[3] << "} ) = " << myBfOrd << "\n";
00276       }
00277     }
00278   }
00279   catch (std::logic_error err){
00280     *outStream << err.what() << "\n\n";
00281     errorFlag = -1000;
00282   };
00283   
00284   *outStream \
00285     << "\n"
00286     << "===============================================================================\n"\
00287     << "| TEST 3: correctness of basis function values                                |\n"\
00288     << "===============================================================================\n";
00289   
00290   outStream -> precision(20);
00291   
00292   // VALUE: Each row gives the 8 correct basis set values at an evaluation point
00293   double basisValues[] = {
00294     // bottom 4 vertices
00295     1.0, 0.0, 0.0, 0.0,  0.0, 0.0, 0.0, 0.0,
00296     0.0, 1.0, 0.0, 0.0,  0.0, 0.0, 0.0, 0.0,
00297     0.0, 0.0, 1.0, 0.0,  0.0, 0.0, 0.0, 0.0,
00298     0.0, 0.0, 0.0, 1.0,  0.0, 0.0, 0.0, 0.0,
00299     // top 4 vertices
00300     0.0, 0.0, 0.0, 0.0,  1.0, 0.0, 0.0, 0.0,
00301     0.0, 0.0, 0.0, 0.0,  0.0, 1.0, 0.0, 0.0,
00302     0.0, 0.0, 0.0, 0.0,  0.0, 0.0, 1.0, 0.0,
00303     0.0, 0.0, 0.0, 0.0,  0.0, 0.0, 0.0, 1.0,
00304     // center {0, 0, 0}
00305     0.125, 0.125, 0.125, 0.125,  0.125, 0.125, 0.125, 0.125,
00306     // faces { 1, 0, 0} and {-1, 0, 0}
00307     0.0,   0.25,  0.25,  0.0,    0.0,   0.25,  0.25,  0.0,
00308     0.25,  0.0,   0.0,   0.25,   0.25,  0.0,   0.0,   0.25,
00309     // faces { 0, 1, 0} and { 0,-1, 0}
00310     0.0,   0.0,   0.25,  0.25,   0.0,   0.0,   0.25,  0.25,
00311     0.25,  0.25,  0.0,   0.0,    0.25,  0.25,  0.0,   0.0,
00312     // faces {0, 0, 1} and {0, 0, -1}
00313     0.0,   0.0,   0.0,   0.0,    0.25,  0.25,  0.25,  0.25,
00314     0.25,  0.25,  0.25,  0.25,   0.0,   0.0,   0.0,   0.0,
00315   };
00316   
00317   // GRAD and D1: each row gives the 3x8 correct values of the gradients of the 8 basis functions
00318   double basisGrads[] = {   
00319     // points 0-3
00320     -0.5,-0.5,-0.5,  0.5, 0.0, 0.0,  0.0, 0.0, 0.0,  0.0, 0.5, 0.0,  0.0, 0.0, 0.5,  0.0, 0.0, 0.0,  0.0, 0.0, 0.0,  0.0, 0.0, 0.0,
00321     -0.5, 0.0, 0.0,  0.5,-0.5,-0.5,  0.0, 0.5, 0.0,  0.0, 0.0, 0.0,  0.0, 0.0, 0.0,  0.0, 0.0, 0.5,  0.0, 0.0, 0.0,  0.0, 0.0, 0.0,
00322      0.0, 0.0, 0.0,  0.0,-0.5, 0.0,  0.5, 0.5,-0.5, -0.5, 0.0, 0.0,  0.0, 0.0, 0.0,  0.0, 0.0, 0.0,  0.0, 0.0, 0.5,  0.0, 0.0, 0.0,
00323      0.0,-0.5, 0.0,  0.0, 0.0, 0.0,  0.5, 0.0, 0.0, -0.5, 0.5,-0.5,  0.0, 0.0, 0.0,  0.0, 0.0, 0.0,  0.0, 0.0, 0.0,  0.0, 0.0, 0.5,
00324      // points 4-7
00325      0.0, 0.0,-0.5,  0.0, 0.0, 0.0,  0.0, 0.0, 0.0,  0.0, 0.0, 0.0, -0.5,-0.5, 0.5,  0.5, 0.0, 0.0,  0.0, 0.0, 0.0,  0.0, 0.5, 0.0,
00326      0.0, 0.0, 0.0,  0.0, 0.0,-0.5,  0.0, 0.0, 0.0,  0.0, 0.0, 0.0, -0.5, 0.0, 0.0,  0.5,-0.5, 0.5,  0.0, 0.5, 0.0,  0.0, 0.0, 0.0,
00327      0.0, 0.0, 0.0,  0.0, 0.0, 0.0,  0.0, 0.0,-0.5,  0.0, 0.0, 0.0,  0.0, 0.0, 0.0,  0.0,-0.5, 0.0,  0.5, 0.5, 0.5, -0.5, 0.0, 0.0,
00328      0.0, 0.0, 0.0,  0.0, 0.0, 0.0,  0.0, 0.0, 0.0,  0.0, 0.0,-0.5,  0.0,-0.5, 0.0,  0.0, 0.0, 0.0,  0.5, 0.0, 0.0, -0.5, 0.5, 0.5,
00329     // point 8
00330     -0.125,-0.125,-0.125,  0.125,-0.125,-0.125,  0.125, 0.125,-0.125, \
00331     -0.125, 0.125,-0.125, -0.125,-0.125, 0.125,  0.125,-0.125, 0.125, \
00332      0.125, 0.125, 0.125, -0.125, 0.125, 0.125,
00333     // point 9
00334     -0.125, 0.0,   0.0,    0.125,-0.25, -0.25,   0.125, 0.25, -0.25,  -0.125, 0.0, 0.0, \
00335     -0.125, 0.0,   0.0,    0.125,-0.25,  0.25,   0.125, 0.25,  0.25,  -0.125, 0.0, 0.0,
00336     // point 10
00337     -0.125,-0.25, -0.25,   0.125, 0.0,   0.0,    0.125, 0.0,   0.0,   -0.125, 0.25, -0.25,\
00338     -0.125,-0.25,  0.25,   0.125, 0.0,   0.0,    0.125, 0.0,   0.0,   -0.125, 0.25,  0.25,
00339     // point 11
00340      0.0,  -0.125, 0.0,    0.0,  -0.125, 0.0,    0.25,  0.125,-0.25,  -0.25,  0.125,-0.25,\
00341      0.0,  -0.125, 0.0,    0.0,  -0.125, 0.0,    0.25,  0.125, 0.25,  -0.25,  0.125, 0.25,
00342     // point 12
00343     -0.25, -0.125,-0.25,   0.25, -0.125,-0.25,   0.0,   0.125, 0.0,    0.0,   0.125, 0.0, \
00344     -0.25, -0.125, 0.25,   0.25, -0.125, 0.25,   0.0,   0.125, 0.0,    0.0,   0.125, 0.0,
00345     // point 13
00346      0.0,   0.0,  -0.125,  0.0,   0.0,  -0.125,  0.0,   0.0,  -0.125,  0.0,   0.0,  -0.125, \
00347     -0.25, -0.25,  0.125,  0.25, -0.25,  0.125,  0.25,  0.25,  0.125, -0.25,  0.25,  0.125,
00348     // point 14
00349     -0.25, -0.25, -0.125,  0.25, -0.25, -0.125,  0.25,  0.25, -0.125, -0.25,  0.25, -0.125, \
00350      0.0,   0.0,   0.125,  0.0,   0.0,   0.125,  0.0,   0.0,   0.125,  0.0,   0.0,   0.125
00351   };
00352   
00353   //D2: flat array with the values of D2 applied to basis functions. Multi-index is (P,F,K)
00354   double basisD2[] = {
00355     // point 0
00356     0, 0.25, 0.25, 0, 0.25, 0, 0, -0.25, -0.25, 0, 0., 0, 0, 0.25, 0., 0, \
00357     0., 0, 0, -0.25, 0., 0, -0.25, 0, 0, 0., -0.25, 0, -0.25, 0, 0, 0., \
00358     0.25, 0, 0., 0, 0, 0., 0., 0, 0., 0, 0, 0., 0., 0, 0.25, 0., \
00359     // point 1
00360     0, 0.25, 0.25, 0, 0., 0, 0, -0.25, -0.25, 0, 0.25, 0, 0, 0.25, 0., 0, \
00361     -0.25, 0, 0, -0.25, 0., 0, 0., 0, 0, 0., -0.25, 0, 0., 0, 0, 0., \
00362     0.25, 0, -0.25, 0, 0, 0., 0., 0, 0.25, 0, 0, 0., 0., 0, 0., 0., \
00363     // Point 2
00364     0, 0.25, 0., 0, 0., 0, 0, -0.25, 0., 0, 0.25, 0, 0, 0.25, -0.25, 0, \
00365     -0.25, 0, 0, -0.25, 0.25, 0, 0., 0, 0, 0., 0., 0, 0., 0, 0, 0., 0., \
00366     0, -0.25, 0, 0, 0., 0.25, 0, 0.25, 0, 0, 0., -0.25, 0, 0., 0., \
00367     // Point 3
00368     0, 0.25, 0., 0, 0.25, 0, 0, -0.25, 0., 0, 0., 0, 0, 0.25, -0.25, 0, \
00369     0., 0, 0, -0.25, 0.25, 0, -0.25, 0, 0, 0., 0., 0, -0.25, 0, 0, 0., \
00370     0., 0, 0., 0, 0, 0., 0.25, 0, 0., 0, 0, 0., -0.25, 0, 0.25, 0.,\
00371     // Point 4
00372     0, 0., 0.25, 0, 0.25, 0, 0, 0., -0.25, 0, 0., 0, 0, 0., 0., 0, 0., 0, \
00373     0, 0., 0., 0, -0.25, 0, 0, 0.25, -0.25, 0, -0.25, 0, 0, -0.25, 0.25, \
00374     0, 0., 0, 0, 0.25, 0., 0, 0., 0, 0, -0.25, 0., 0, 0.25, 0., \
00375     // Point 5
00376     0, 0., 0.25, 0, 0., 0, 0, 0., -0.25, 0, 0.25, 0, 0, 0., 0., 0, -0.25, \
00377     0, 0, 0., 0., 0, 0., 0, 0, 0.25, -0.25, 0, 0., 0, 0, -0.25, 0.25, 0, \
00378     -0.25, 0, 0, 0.25, 0., 0, 0.25, 0, 0, -0.25, 0., 0, 0., 0., \
00379     // Point 6
00380     0, 0., 0., 0, 0., 0, 0, 0., 0., 0, 0.25, 0, 0, 0., -0.25, 0, -0.25, \
00381     0, 0, 0., 0.25, 0, 0., 0, 0, 0.25, 0., 0, 0., 0, 0, -0.25, 0., 0, \
00382     -0.25, 0, 0, 0.25, 0.25, 0, 0.25, 0, 0, -0.25, -0.25, 0, 0., 0., \
00383     // Point 7
00384     0, 0., 0., 0, 0.25, 0, 0, 0., 0., 0, 0., 0, 0, 0., -0.25, 0, 0., 0, \
00385     0, 0., 0.25, 0, -0.25, 0, 0, 0.25, 0., 0, -0.25, 0, 0, -0.25, 0., 0, \
00386     0., 0, 0, 0.25, 0.25, 0, 0., 0, 0, -0.25, -0.25, 0, 0.25, 0., \
00387     // Point 8
00388     0, 0.125, 0.125, 0, 0.125, 0, 0, -0.125, -0.125, 0, 0.125, 0, 0, \
00389     0.125, -0.125, 0, -0.125, 0, 0, -0.125, 0.125, 0, -0.125, 0, 0, \
00390     0.125, -0.125, 0, -0.125, 0, 0, -0.125, 0.125, 0, -0.125, 0, 0, \
00391     0.125, 0.125, 0, 0.125, 0, 0, -0.125, -0.125, 0, 0.125, 0., \
00392     // Point 9
00393     0, 0.125, 0.125, 0, 0., 0, 0, -0.125, -0.125, 0, 0.25, 0, 0, 0.125, \
00394     -0.125, 0, -0.25, 0, 0, -0.125, 0.125, 0, 0., 0, 0, 0.125, -0.125, 0, \
00395     0., 0, 0, -0.125, 0.125, 0, -0.25, 0, 0, 0.125, 0.125, 0, 0.25, 0, 0, \
00396     -0.125, -0.125, 0, 0., 0., \
00397     // Point 10
00398     0, 0.125, 0.125, 0, 0.25, 0, 0, -0.125, -0.125, 0, 0., 0, 0, 0.125, \
00399     -0.125, 0, 0., 0, 0, -0.125, 0.125, 0, -0.25, 0, 0, 0.125, -0.125, 0, \
00400     -0.25, 0, 0, -0.125, 0.125, 0, 0., 0, 0, 0.125, 0.125, 0, 0., 0, 0, \
00401     -0.125, -0.125, 0, 0.25, 0., \
00402     // Point 11
00403     0, 0.125, 0., 0, 0.125, 0, 0, -0.125, 0., 0, 0.125, 0, 0, 0.125, \
00404     -0.25, 0, -0.125, 0, 0, -0.125, 0.25, 0, -0.125, 0, 0, 0.125, 0., 0, \
00405     -0.125, 0, 0, -0.125, 0., 0, -0.125, 0, 0, 0.125, 0.25, 0, 0.125, 0, \
00406     0, -0.125, -0.25, 0, 0.125, 0., \
00407     // Point 12
00408     0, 0.125, 0.25, 0, 0.125, 0, 0, -0.125, -0.25, 0, 0.125, 0, 0, 0.125, \
00409     0., 0, -0.125, 0, 0, -0.125, 0., 0, -0.125, 0, 0, 0.125, -0.25, 0, \
00410     -0.125, 0, 0, -0.125, 0.25, 0, -0.125, 0, 0, 0.125, 0., 0, 0.125, 0, \
00411     0, -0.125, 0., 0, 0.125, 0., \
00412     // Point 13
00413     0, 0., 0.125, 0, 0.125, 0, 0, 0., -0.125, 0, 0.125, 0, 0, 0., -0.125, \
00414     0, -0.125, 0, 0, 0., 0.125, 0, -0.125, 0, 0, 0.25, -0.125, 0, -0.125, \
00415     0, 0, -0.25, 0.125, 0, -0.125, 0, 0, 0.25, 0.125, 0, 0.125, 0, 0, \
00416     -0.25, -0.125, 0, 0.125, 0., \
00417     // Point 14
00418     0, 0.25, 0.125, 0, 0.125, 0, 0, -0.25, -0.125, 0, 0.125, 0, 0, 0.25, \
00419     -0.125, 0, -0.125, 0, 0, -0.25, 0.125, 0, -0.125, 0, 0, 0., -0.125, \
00420     0, -0.125, 0, 0, 0., 0.125, 0, -0.125, 0, 0, 0., 0.125, 0, 0.125, 0, \
00421     0, 0., -0.125, 0, 0.125, 0.
00422   };
00423   
00424   try{
00425         
00426     // Dimensions for the output arrays:
00427     int numFields = hexBasis.getCardinality();
00428     int numPoints = hexNodes.dimension(0);
00429     int spaceDim  = hexBasis.getBaseCellTopology().getDimension();
00430     int D2Cardin  = Intrepid::getDkCardinality(OPERATOR_D2, spaceDim);
00431     
00432     // Generic array for values, grads, curls, etc. that will be properly sized before each call
00433     FieldContainer<double> vals;
00434     
00435     // Check VALUE of basis functions: resize vals to rank-2 container:
00436     vals.resize(numFields, numPoints);
00437     hexBasis.getValues(vals, hexNodes, OPERATOR_VALUE);
00438     for (int i = 0; i < numFields; i++) {
00439       for (int j = 0; j < numPoints; j++) {
00440           int l =  i + j * numFields;
00441            if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
00442              errorFlag++;
00443              *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00444 
00445              // Output the multi-index of the value where the error is:
00446              *outStream << " At multi-index { ";
00447              *outStream << i << " ";*outStream << j << " ";
00448              *outStream << "}  computed value: " << vals(i,j)
00449                << " but reference value: " << basisValues[l] << "\n";
00450          }
00451       }
00452     }
00453     
00454     // Check GRAD of basis function: resize vals to rank-3 container
00455     vals.resize(numFields, numPoints, spaceDim);
00456     hexBasis.getValues(vals, hexNodes, OPERATOR_GRAD);
00457     for (int i = 0; i < numFields; i++) {
00458       for (int j = 0; j < numPoints; j++) {
00459         for (int k = 0; k < spaceDim; k++) {
00460            int l = k + i * spaceDim + j * spaceDim * numFields;
00461            if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
00462              errorFlag++;
00463              *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00464 
00465              // Output the multi-index of the value where the error is:
00466              *outStream << " At multi-index { ";
00467              *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
00468              *outStream << "}  computed grad component: " << vals(i,j,k)
00469                << " but reference grad component: " << basisGrads[l] << "\n";
00470             }
00471          }
00472       }
00473     }
00474     
00475     // Check D1 of basis function (do not resize vals because it has the correct size: D1 = GRAD)
00476     hexBasis.getValues(vals, hexNodes, OPERATOR_D1);
00477     for (int i = 0; i < numFields; i++) {
00478       for (int j = 0; j < numPoints; j++) {
00479         for (int k = 0; k < spaceDim; k++) {
00480            int l = k + i * spaceDim + j * spaceDim * numFields;
00481            if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
00482              errorFlag++;
00483              *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00484 
00485              // Output the multi-index of the value where the error is:
00486              *outStream << " At multi-index { ";
00487              *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
00488              *outStream << "}  computed D1 component: " << vals(i,j,k)
00489                << " but reference D1 component: " << basisGrads[l] << "\n";
00490             }
00491          }
00492       }
00493     }
00494 
00495     
00496     // Check D2 of basis function
00497     vals.resize(numFields, numPoints, D2Cardin);    
00498     hexBasis.getValues(vals, hexNodes, OPERATOR_D2);
00499     for (int i = 0; i < numFields; i++) {
00500       for (int j = 0; j < numPoints; j++) {
00501         for (int k = 0; k < D2Cardin; k++) {
00502            int l = k + i * D2Cardin + j * D2Cardin * numFields;
00503            if (std::abs(vals(i,j,k) - basisD2[l]) > INTREPID_TOL) {
00504              errorFlag++;
00505              *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00506 
00507              // Output the multi-index of the value where the error is:
00508              *outStream << " At multi-index { ";
00509              *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
00510              *outStream << "}  computed D2 component: " << vals(i,j,k)
00511                << " but reference D2 component: " << basisD2[l] << "\n";
00512             }
00513          }
00514       }
00515     }
00516 
00517     // Check all higher derivatives - must be zero. 
00518     for(EOperator op = OPERATOR_D3; op < OPERATOR_MAX; op++) {
00519       
00520       // The last dimension is the number of kth derivatives and needs to be resized for every Dk
00521       int DkCardin  = Intrepid::getDkCardinality(op, spaceDim);
00522       vals.resize(numFields, numPoints, DkCardin);    
00523 
00524       hexBasis.getValues(vals, hexNodes, op);
00525       for (int i = 0; i < vals.size(); i++) {
00526         if (std::abs(vals[i]) > INTREPID_TOL) {
00527           errorFlag++;
00528           *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00529           
00530           // Get the multi-index of the value where the error is and the operator order
00531           std::vector<int> myIndex;
00532           vals.getMultiIndex(myIndex,i);
00533           int ord = Intrepid::getOperatorOrder(op);
00534           *outStream << " At multi-index { ";
00535           for(int j = 0; j < vals.rank(); j++) {
00536             *outStream << myIndex[j] << " ";
00537           }
00538           *outStream << "}  computed D"<< ord <<" component: " << vals[i] 
00539             << " but reference D" << ord << " component:  0 \n";
00540         }
00541       }
00542     }    
00543   }
00544   
00545   // Catch unexpected errors
00546   catch (std::logic_error err) {
00547     *outStream << err.what() << "\n\n";
00548     errorFlag = -1000;
00549   };
00550 
00551   *outStream \
00552     << "\n"
00553     << "===============================================================================\n"\
00554     << "| TEST 4: correctness of DoF locations                                        |\n"\
00555     << "===============================================================================\n";
00556 
00557   try{
00558     Teuchos::RCP<Basis<double, FieldContainer<double> > > basis =
00559       Teuchos::rcp(new Basis_HGRAD_HEX_C1_FEM<double, FieldContainer<double> >);
00560     Teuchos::RCP<DofCoordsInterface<FieldContainer<double> > > coord_iface =
00561       Teuchos::rcp_dynamic_cast<DofCoordsInterface<FieldContainer<double> > >(basis);
00562 
00563     FieldContainer<double> cvals;
00564     FieldContainer<double> bvals(basis->getCardinality(), basis->getCardinality());
00565 
00566     // Check exceptions.
00567 #ifdef HAVE_INTREPID_DEBUG
00568     cvals.resize(1,2,3);
00569     INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
00570     cvals.resize(3,2);
00571     INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
00572     cvals.resize(8,2);
00573     INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException );
00574 #endif
00575     cvals.resize(8,3);
00576     INTREPID_TEST_COMMAND( coord_iface->getDofCoords(cvals), throwCounter, nException ); nException--;
00577     // Check if number of thrown exceptions matches the one we expect
00578     if (throwCounter != nException) {
00579       errorFlag++;
00580       *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00581     }
00582 
00583     // Check mathematical correctness.
00584     basis->getValues(bvals, cvals, OPERATOR_VALUE);
00585     char buffer[120];
00586     for (int i=0; i<bvals.dimension(0); i++) {
00587       for (int j=0; j<bvals.dimension(1); j++) {
00588         if ((i != j) && (std::abs(bvals(i,j) - 0.0) > INTREPID_TOL)) {
00589           errorFlag++;
00590           sprintf(buffer, "\nValue of basis function %d at (%6.4e, %6.4e, %6.4e) is %6.4e but should be %6.4e!\n", i, cvals(i,0), cvals(i,1), cvals(i,2), bvals(i,j), 0.0);
00591           *outStream << buffer;
00592         }
00593         else if ((i == j) && (std::abs(bvals(i,j) - 1.0) > INTREPID_TOL)) {
00594           errorFlag++;
00595           sprintf(buffer, "\nValue of basis function %d at (%6.4e, %6.4e, %6.4e) is %6.4e but should be %6.4e!\n", i, cvals(i,0), cvals(i,1), cvals(i,2), bvals(i,j), 1.0);
00596           *outStream << buffer;
00597         }
00598       }
00599     }
00600 
00601   }
00602   catch (std::logic_error err){
00603     *outStream << err.what() << "\n\n";
00604     errorFlag = -1000;
00605   };
00606   
00607   if (errorFlag != 0)
00608     std::cout << "End Result: TEST FAILED\n";
00609   else
00610     std::cout << "End Result: TEST PASSED\n";
00611   
00612   // reset format state of std::cout
00613   std::cout.copyfmt(oldFormatState);
00614   
00615   return errorFlag;
00616 }