Intrepid: http://trilinos.sandia.gov/packages/docs/r10.0/packages/intrepid/test/Discretization/Basis/HGRAD_WEDGE_C2_FEM/test

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 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Pavel Bochev  (pbboche@sandia.gov),
00025 //                    Denis Ridzal  (dridzal@sandia.gov),
00026 //                    Kara Peterson (kjpeter@sandia.gov).
00027 //
00028 // ************************************************************************
00029 // @HEADER
00030 
00035 #include "Intrepid_FieldContainer.hpp"
00036 #include "Intrepid_HGRAD_WEDGE_C2_FEM.hpp"
00037 #include "Teuchos_oblackholestream.hpp"
00038 #include "Teuchos_RCP.hpp"
00039 #include "Teuchos_GlobalMPISession.hpp"
00040 
00041 using namespace std;
00042 using namespace Intrepid;
00043 
00044 #define INTREPID_TEST_COMMAND( S , throwCounter, nException )                                                              \
00045 {                                                                                                                          \
00046   ++nException;                                                                                                            \
00047   try {                                                                                                                    \
00048     S ;                                                                                                                    \
00049   }                                                                                                                        \
00050   catch (std::logic_error err) {                                                                                           \
00051       ++throwCounter;                                                                                                      \
00052       *outStream << "Expected Error " << nException << " -------------------------------------------------------------\n"; \
00053       *outStream << err.what() << '\n';                                                                                    \
00054       *outStream << "-------------------------------------------------------------------------------" << "\n\n";           \
00055   };                                                                                                                       \
00056 }
00057 
00058 int main(int argc, char *argv[]) {
00059   
00060   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00061 
00062   // This little trick lets us print to std::cout only if
00063   // a (dummy) command-line argument is provided.
00064   int iprint     = argc - 1;
00065   Teuchos::RCP<std::ostream> outStream;
00066   Teuchos::oblackholestream bhs; // outputs nothing
00067   if (iprint > 0)
00068     outStream = Teuchos::rcp(&std::cout, false);
00069   else
00070     outStream = Teuchos::rcp(&bhs, false);
00071   
00072   // Save the format state of the original std::cout.
00073   Teuchos::oblackholestream oldFormatState;
00074   oldFormatState.copyfmt(std::cout);
00075   
00076   *outStream \
00077     << "===============================================================================\n" \
00078     << "|                                                                             |\n" \
00079     << "|                 Unit Test (Basis_HGRAD_WEDGE_C2_FEM)                        |\n" \
00080     << "|                                                                             |\n" \
00081     << "|     1) Conversion of Dof tags into Dof ordinals and back                    |\n" \
00082     << "|     2) Basis values for VALUE, GRAD, and Dk operators                       |\n" \
00083     << "|                                                                             |\n" \
00084     << "|  Questions? Contact  Pavel Bochev  (pbboche@sandia.gov),                    |\n" \
00085     << "|                      Denis Ridzal  (dridzal@sandia.gov),                    |\n" \
00086     << "|                      Kara Peterson (kjpeter@sandia.gov).                    |\n" \
00087     << "|                                                                             |\n" \
00088     << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00089     << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00090     << "|                                                                             |\n" \
00091     << "===============================================================================\n"\
00092     << "| TEST 1: Basis creation, exception testing                                   |\n"\
00093     << "===============================================================================\n";
00094   
00095   // Define basis and error flag
00096   Basis_HGRAD_WEDGE_C2_FEM<double, FieldContainer<double> > wedgeBasis;
00097   int errorFlag = 0;
00098 
00099   // Initialize throw counter for exception testing
00100   int nException     = 0;
00101   int throwCounter   = 0;
00102 
00103   // Nodes of Wedde<18>: vertices, edge midpoints, Quadrilateral face centers 
00104   FieldContainer<double> wedgeNodes(18, 3);
00105   wedgeNodes(0,0) =  0.0;  wedgeNodes(0,1) =  0.0;  wedgeNodes(0,2) = -1.0;  
00106   wedgeNodes(1,0) =  1.0;  wedgeNodes(1,1) =  0.0;  wedgeNodes(1,2) = -1.0;  
00107   wedgeNodes(2,0) =  0.0;  wedgeNodes(2,1) =  1.0;  wedgeNodes(2,2) = -1.0;
00108   wedgeNodes(3,0) =  0.0;  wedgeNodes(3,1) =  0.0;  wedgeNodes(3,2) =  1.0;  
00109   wedgeNodes(4,0) =  1.0;  wedgeNodes(4,1) =  0.0;  wedgeNodes(4,2) =  1.0;  
00110   wedgeNodes(5,0) =  0.0;  wedgeNodes(5,1) =  1.0;  wedgeNodes(5,2) =  1.0;
00111     
00112   wedgeNodes(6,0) =  0.5;  wedgeNodes(6,1) =  0.0;  wedgeNodes(6,2) = -1.0;  
00113   wedgeNodes(7,0) =  0.5;  wedgeNodes(7,1) =  0.5;  wedgeNodes(7,2) = -1.0;  
00114   wedgeNodes(8,0) =  0.0;  wedgeNodes(8,1) =  0.5;  wedgeNodes(8,2) = -1.0;
00115   wedgeNodes(9,0) =  0.0;  wedgeNodes(9,1) =  0.0;  wedgeNodes(9,2) =  0.0;
00116   wedgeNodes(10,0)=  1.0;  wedgeNodes(10,1)=  0.0;  wedgeNodes(10,2)=  0.0;  
00117   wedgeNodes(11,0)=  0.0;  wedgeNodes(11,1)=  1.0;  wedgeNodes(11,2)=  0.0;  
00118  
00119   wedgeNodes(12,0)=  0.5;  wedgeNodes(12,1)=  0.0;  wedgeNodes(12,2)=  1.0;  
00120   wedgeNodes(13,0)=  0.5;  wedgeNodes(13,1)=  0.5;  wedgeNodes(13,2)=  1.0;  
00121   wedgeNodes(14,0)=  0.0;  wedgeNodes(14,1)=  0.5;  wedgeNodes(14,2)=  1.0;  
00122   wedgeNodes(15,0)=  0.5;  wedgeNodes(15,1)=  0.0;  wedgeNodes(15,2)=  0.0;  
00123   wedgeNodes(16,0)=  0.5;  wedgeNodes(16,1)=  0.5;  wedgeNodes(16,2)=  0.0;  
00124   wedgeNodes(17,0)=  0.0;  wedgeNodes(17,1)=  0.5;  wedgeNodes(17,2)=  0.0;  
00125 
00126 
00127   // Generic array for the output values; needs to be properly resized depending on the operator type
00128   FieldContainer<double> vals;
00129 
00130   try{
00131     // exception #1: CURL cannot be applied to scalar functions
00132     // resize vals to rank-3 container with dimensions (num. points, num. basis functions)
00133     vals.resize(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), 3 );
00134     INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_DIV), throwCounter, nException );
00135 
00136     // exception #2: DIV cannot be applied to scalar functions
00137     // resize vals to rank-2 container with dimensions (num. points, num. basis functions)
00138     vals.resize(wedgeBasis.getCardinality(), wedgeNodes.dimension(0) );
00139     INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_DIV), throwCounter, nException );
00140         
00141     // Exceptions 3-7: 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 #3
00144     INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(3,0,0), throwCounter, nException );
00145     // exception #4
00146     INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(1,1,1), throwCounter, nException );
00147     // exception #5
00148     INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(0,9,0), throwCounter, nException );
00149     // exception #6
00150     INTREPID_TEST_COMMAND( wedgeBasis.getDofTag(18), throwCounter, nException );
00151     // exception #7
00152     INTREPID_TEST_COMMAND( wedgeBasis.getDofTag(-1), throwCounter, nException );
00153     
00154 #ifdef HAVE_INTREPID_DEBUG
00155     // Exceptions 8-18 test exception handling with incorrectly dimensioned input/output arrays
00156     // exception #8: input points array must be of rank-2
00157     FieldContainer<double> badPoints1(4, 5, 3);
00158     INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
00159     
00160     // exception #9 dimension 1 in the input point array must equal space dimension of the cell
00161     FieldContainer<double> badPoints2(4, wedgeBasis.getBaseCellTopology().getDimension() + 1);
00162     INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
00163     
00164     // exception #10 output values must be of rank-2 for OPERATOR_VALUE
00165     FieldContainer<double> badVals1(4, 3, 1);
00166     INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals1, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
00167     
00168     // exception #11 output values must be of rank-3 for OPERATOR_GRAD
00169     FieldContainer<double> badVals2(4, 3);
00170     INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals2, wedgeNodes, OPERATOR_GRAD), throwCounter, nException );
00171     
00172     // exception #12 output values must be of rank-3 for OPERATOR_D1
00173     INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals2, wedgeNodes, OPERATOR_D1), throwCounter, nException );
00174     
00175     // exception #13 output values must be of rank-3 for OPERATOR_D2
00176     INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals2, wedgeNodes, OPERATOR_D2), throwCounter, nException );
00177     
00178     // exception #14 incorrect 0th dimension of output array (must equal number of basis functions)
00179     FieldContainer<double> badVals3(wedgeBasis.getCardinality() + 1, wedgeNodes.dimension(0));
00180     INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals3, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
00181     
00182     // exception #15 incorrect 1st dimension of output array (must equal number of points)
00183     FieldContainer<double> badVals4(wedgeBasis.getCardinality(), wedgeNodes.dimension(0) + 1);
00184     INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals4, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
00185     
00186     // exception #16: incorrect 2nd dimension of output array (must equal the space dimension)
00187     FieldContainer<double> badVals5(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), wedgeBasis.getBaseCellTopology().getDimension() - 1);
00188     INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals5, wedgeNodes, OPERATOR_GRAD), throwCounter, nException );
00189     
00190     // exception #17: incorrect 2nd dimension of output array (must equal D2 cardinality in 3D)
00191     FieldContainer<double> badVals6(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), 40);
00192     INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals6, wedgeNodes, OPERATOR_D2), throwCounter, nException );
00193     
00194     // exception #18: incorrect 2nd dimension of output array (must equal D3 cardinality in 3D)
00195     INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals6, wedgeNodes, 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 - 18
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 = wedgeBasis.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  = wedgeBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
00224       
00225       std::vector<int> myTag = wedgeBasis.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 < wedgeBasis.getCardinality(); bfOrd++) {
00247       std::vector<int> myTag  = wedgeBasis.getDofTag(bfOrd);
00248       int myBfOrd = wedgeBasis.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 function values in (F,P) format
00278   double basisValues[] = {
00279     1.00, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.00, 0, \
00280     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.00, 0, 0, 0, 0, \
00281     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.00, 0, 0, 0, 0, 0, 0, 0, \
00282     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.00, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
00283     0, 0, 0, 0, 0, 0, 0, 0, 1.00, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
00284     0, 0, 0, 0, 0, 1.00, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
00285     0, 0, 1.00, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
00286     1.00, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.00, 0, \
00287     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.00, 0, 0, 0, 0, \
00288     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.00, 0, 0, 0, 0, 0, 0, 0, \
00289     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.00, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
00290     0, 0, 0, 0, 0, 0, 0, 0, 1.00, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
00291     0, 0, 0, 0, 0, 1.00, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
00292     0, 0, 1.00, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, \
00293     1.00, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.00  };
00294   
00295   // GRAD, D1, D2, D3 and D4 test values are stored in files due to their large size
00296   std::string     fileName;
00297   std::ifstream   dataFile;
00298   
00299   // GRAD and D1 values are stored in (F,P,D) format in a data file. Read file and do the test
00300   std::vector<double> basisGrads;           // Flat array for the gradient values.
00301   
00302   fileName = "./testdata/WEDGE_C2_GradVals.dat";
00303   dataFile.open(fileName.c_str());
00304   TEST_FOR_EXCEPTION( !dataFile.good(), std::logic_error,
00305                       ">>> ERROR (HGRAD_WEDGE_C2/test01): could not open GRAD values data file, test aborted.");
00306   while (!dataFile.eof() ){
00307     double temp;
00308     string line;                            // string for one line of input file
00309     std::getline(dataFile, line);           // get next line from file
00310     stringstream data_line(line);           // convert to stringstream
00311     while(data_line >> temp){               // extract value from line
00312       basisGrads.push_back(temp);           // push into vector
00313     }
00314   }
00315   // It turns out that just closing and then opening the ifstream variable does not reset it
00316   // and subsequent open() command fails. One fix is to explicitely clear the ifstream, or
00317   // scope the variables.
00318   dataFile.close();
00319   dataFile.clear();
00320   
00321   
00322   //D2: flat array with the values of D2 applied to basis functions. Multi-index is (F,P,D2cardinality)
00323   std::vector<double> basisD2;
00324   
00325   fileName = "./testdata/WEDGE_C2_D2Vals.dat";
00326   dataFile.open(fileName.c_str());
00327   TEST_FOR_EXCEPTION( !dataFile.good(), std::logic_error,
00328                       ">>> ERROR (HGRAD_WEDGE_C2/test01): could not open D2 values data file, test aborted.");
00329   
00330   while (!dataFile.eof() ){
00331     double temp;
00332     string line;                            // string for one line of input file
00333     std::getline(dataFile, line);           // get next line from file
00334     stringstream data_line(line);           // convert to stringstream
00335     while(data_line >> temp){               // extract value from line
00336       basisD2.push_back(temp);              // push into vector
00337     }
00338   }
00339   dataFile.close();
00340   dataFile.clear();
00341   
00342   
00343   //D3: flat array with the values of D3 applied to basis functions. Multi-index is (F,P,D3cardinality)
00344   std::vector<double> basisD3;
00345   
00346   fileName = "./testdata/WEDGE_C2_D3Vals.dat";
00347   dataFile.open(fileName.c_str());
00348   TEST_FOR_EXCEPTION( !dataFile.good(), std::logic_error,
00349                       ">>> ERROR (HGRAD_WEDGE_C2/test01): could not open D3 values data file, test aborted.");
00350   
00351   while (!dataFile.eof() ){
00352     double temp;
00353     string line;                            // string for one line of input file
00354     std::getline(dataFile, line);           // get next line from file
00355     stringstream data_line(line);           // convert to stringstream
00356     while(data_line >> temp){               // extract value from line
00357       basisD3.push_back(temp);              // push into vector
00358     }
00359   }
00360   dataFile.close();
00361   dataFile.clear();
00362   
00363   
00364   //D4: flat array with the values of D3 applied to basis functions. Multi-index is (F,P,D4cardinality)
00365   std::vector<double> basisD4;
00366   
00367   fileName = "./testdata/WEDGE_C2_D4Vals.dat";
00368   dataFile.open(fileName.c_str());
00369   TEST_FOR_EXCEPTION( !dataFile.good(), std::logic_error,
00370                       ">>> ERROR (HGRAD_WEDGE_C2/test01): could not open D4 values data file, test aborted.");
00371   
00372   while (!dataFile.eof() ){
00373     double temp;
00374     string line;                            // string for one line of input file
00375     std::getline(dataFile, line);           // get next line from file
00376     stringstream data_line(line);           // convert to stringstream
00377     while(data_line >> temp){               // extract value from line
00378       basisD4.push_back(temp);              // push into vector
00379     }
00380   }
00381   dataFile.close();
00382   dataFile.clear();
00383 
00384   
00385   try{
00386         
00387     // Dimensions for the output arrays:
00388     int numFields = wedgeBasis.getCardinality();
00389     int numPoints = wedgeNodes.dimension(0);
00390     int spaceDim  = wedgeBasis.getBaseCellTopology().getDimension();
00391     
00392     // Generic array for values, grads, curls, etc. that will be properly sized before each call
00393     FieldContainer<double> vals;
00394     
00395     // Check VALUE of basis functions: resize vals to rank-2 container:
00396     vals.resize(numFields, numPoints);
00397     wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_VALUE);
00398     for (int i = 0; i < numFields; i++) {
00399       for (int j = 0; j < numPoints; j++) {
00400           int l =  i + j * numFields;
00401            if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
00402              errorFlag++;
00403              *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00404 
00405              // Output the multi-index of the value where the error is:
00406              *outStream << " At multi-index { ";
00407              *outStream << i << " ";*outStream << j << " ";
00408              *outStream << "}  computed value: " << vals(i,j)
00409                << " but reference value: " << basisValues[l] << "\n";
00410          }
00411       }
00412     }
00413 
00414     
00415     
00416     // Check GRAD of basis function: resize vals to rank-3 container
00417     vals.resize(numFields, numPoints, spaceDim);
00418     wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_GRAD);
00419     for (int i = 0; i < numFields; i++) {
00420       for (int j = 0; j < numPoints; j++) {
00421         for (int k = 0; k < spaceDim; k++) {
00422           
00423           // basisGrads is (F,P,D), compute offset:
00424           int l = k + j * spaceDim + i * spaceDim * numPoints;
00425           if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
00426             errorFlag++;
00427             *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00428             
00429             // Output the multi-index of the value where the error is:
00430             *outStream << " At multi-index { ";
00431             *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
00432             *outStream << "}  computed grad component: " << vals(i,j,k)
00433               << " but reference grad component: " << basisGrads[l] << "\n";
00434           }
00435         }
00436       }
00437     }
00438     
00439     // Check D1 of basis function (do not resize vals because it has the correct size: D1 = GRAD)
00440     wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_D1);
00441     for (int i = 0; i < numFields; i++) {
00442       for (int j = 0; j < numPoints; j++) {
00443         for (int k = 0; k < spaceDim; k++) {
00444           
00445           // basisGrads is (F,P,D), compute offset:
00446           int l = k + j * spaceDim + i * spaceDim * numPoints;
00447           if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
00448             errorFlag++;
00449             *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00450             
00451             // Output the multi-index of the value where the error is:
00452             *outStream << " At multi-index { ";
00453             *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
00454             *outStream << "}  computed D1 component: " << vals(i,j,k)
00455               << " but reference D1 component: " << basisGrads[l] << "\n";
00456           }
00457         }
00458       }
00459     }
00460     
00461     
00462     // Check D2 of basis function
00463     int D2cardinality = Intrepid::getDkCardinality(OPERATOR_D2, spaceDim);
00464     vals.resize(numFields, numPoints, D2cardinality);    
00465     wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_D2);
00466     for (int i = 0; i < numFields; i++) {
00467       for (int j = 0; j < numPoints; j++) {
00468         for (int k = 0; k < D2cardinality; k++) {
00469           
00470           // basisD2 is (F,P,Dk), compute offset:
00471           int l = k + j * D2cardinality + i * D2cardinality * numPoints;
00472           if (std::abs(vals(i,j,k) - basisD2[l]) > INTREPID_TOL) {
00473             errorFlag++;
00474             *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00475             
00476             // Output the multi-index of the value where the error is:
00477             *outStream << " At multi-index { ";
00478             *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
00479             *outStream << "}  computed D2 component: " << vals(i,j,k)
00480               << " but reference D2 component: " << basisD2[l] << "\n";
00481           }
00482         }
00483       }
00484     }
00485     
00486     
00487     // Check D3 of basis function
00488     int D3cardinality = Intrepid::getDkCardinality(OPERATOR_D3, spaceDim);
00489     vals.resize(numFields, numPoints, D3cardinality);    
00490     wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_D3);
00491     
00492     for (int i = 0; i < numFields; i++) {
00493       for (int j = 0; j < numPoints; j++) {
00494         for (int k = 0; k < D3cardinality; k++) {
00495           
00496           // basisD3 is (F,P,Dk), compute offset:
00497           int l = k + j * D3cardinality + i * D3cardinality * numPoints;
00498           if (std::abs(vals(i,j,k) - basisD3[l]) > INTREPID_TOL) {
00499             errorFlag++;
00500             *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00501             
00502             // Output the multi-index of the value where the error is:
00503             *outStream << " At multi-index { ";
00504             *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
00505             *outStream << "}  computed D3 component: " << vals(i,j,k)
00506               << " but reference D3 component: " << basisD3[l] << "\n";
00507           }
00508         }
00509       }
00510     }
00511     
00512     
00513     // Check D4 of basis function
00514     int D4cardinality = Intrepid::getDkCardinality(OPERATOR_D4, spaceDim);
00515     vals.resize(numFields, numPoints, D4cardinality);    
00516     wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_D4);
00517     for (int i = 0; i < numFields; i++) {
00518       for (int j = 0; j < numPoints; j++) {
00519         for (int k = 0; k < D4cardinality; k++) {
00520           
00521           // basisD4 is (F,P,Dk), compute offset:
00522           int l = k + j * D4cardinality + i * D4cardinality * numPoints;
00523           if (std::abs(vals(i,j,k) - basisD4[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 D4 component: " << vals(i,j,k)
00531               << " but reference D4 component: " << basisD2[l] << "\n";
00532           }
00533         }
00534       }
00535     }
00536     
00537         
00538     // Check all higher derivatives - must be zero. 
00539     for(EOperator op = OPERATOR_D5; op < OPERATOR_MAX; op++) {
00540       
00541       // The last dimension is the number of kth derivatives and needs to be resized for every Dk
00542       int DkCardin  = Intrepid::getDkCardinality(op, spaceDim);
00543       vals.resize(numFields, numPoints, DkCardin);    
00544 
00545       wedgeBasis.getValues(vals, wedgeNodes, op);
00546       for (int i = 0; i < vals.size(); i++) {
00547         if (std::abs(vals[i]) > INTREPID_TOL) {
00548           errorFlag++;
00549           *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00550           
00551           // Get the multi-index of the value where the error is and the operator order
00552           std::vector<int> myIndex;
00553           vals.getMultiIndex(myIndex,i);
00554           int ord = Intrepid::getOperatorOrder(op);
00555           *outStream << " At multi-index { ";
00556           for(int j = 0; j < vals.rank(); j++) {
00557             *outStream << myIndex[j] << " ";
00558           }
00559           *outStream << "}  computed D"<< ord <<" component: " << vals[i] 
00560             << " but reference D" << ord << " component:  0 \n";
00561         }
00562       }
00563     }    
00564   }
00565   
00566   // Catch unexpected errors
00567   catch (std::logic_error err) {
00568     *outStream << err.what() << "\n\n";
00569     errorFlag = -1000;
00570   };
00571   
00572   if (errorFlag != 0)
00573     std::cout << "End Result: TEST FAILED\n";
00574   else
00575     std::cout << "End Result: TEST PASSED\n";
00576   
00577   // reset format state of std::cout
00578   std::cout.copyfmt(oldFormatState);
00579   
00580   return errorFlag;
00581 }
http://trilinos.sandia.gov/packages/docs/r10.0/packages/intrepid/test/Discretization/Basis/HGRAD_WEDGE_C2_FEM/test_01.cpp
