Intrepid: http://trilinos.sandia.gov/packages/docs/r10.0/packages/intrepid/test/Discretization/Basis/HDIV_TET_I1_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_HDIV_TET_I1_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_HDIV_TET_I1_FEM)                           |\n" \
00080     << "|                                                                             |\n" \
00081     << "|     1) Conversion of Dof tags into Dof ordinals and back                    |\n" \
00082     << "|     3) Exception tests on input arguments and invalid operators             |\n" \
00083     << "|     2) Basis values for VALUE and DIV operators                             |\n" \
00084     << "|                                                                             |\n" \
00085     << "|  Questions? Contact  Pavel Bochev  (pbboche@sandia.gov),                    |\n" \
00086     << "|                      Denis Ridzal  (dridzal@sandia.gov),                    |\n" \
00087     << "|                      Kara Peterson (kjpeter@sandia.gov).                    |\n" \
00088     << "|                                                                             |\n" \
00089     << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00090     << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00091     << "|                                                                             |\n" \
00092     << "===============================================================================\n"\
00093     << "| TEST 1: Basis creation, exception testing                                   |\n"\
00094     << "===============================================================================\n";
00095   
00096   // Define basis and error flag
00097   Basis_HDIV_TET_I1_FEM<double, FieldContainer<double> > tetBasis;
00098   int errorFlag = 0;
00099 
00100   // Initialize throw counter for exception testing
00101   int nException     = 0;
00102   int throwCounter   = 0;
00103 
00104   // Define array containing the 4 vertices of the reference TET and its 6 edge midpoints.
00105   FieldContainer<double> tetNodes(10, 3);
00106   tetNodes(0,0) =  0.0;  tetNodes(0,1) =  0.0;  tetNodes(0,2) =  0.0;
00107   tetNodes(1,0) =  1.0;  tetNodes(1,1) =  0.0;  tetNodes(1,2) =  0.0;
00108   tetNodes(2,0) =  0.0;  tetNodes(2,1) =  1.0;  tetNodes(2,2) =  0.0;
00109   tetNodes(3,0) =  0.0;  tetNodes(3,1) =  0.0;  tetNodes(3,2) =  1.0;
00110   tetNodes(4,0) =  0.5;  tetNodes(4,1) =  0.0;  tetNodes(4,2) =  0.0;
00111   tetNodes(5,0) =  0.5;  tetNodes(5,1) =  0.5;  tetNodes(5,2) =  0.0;
00112   tetNodes(6,0) =  0.0;  tetNodes(6,1) =  0.5;  tetNodes(6,2) =  0.0;
00113   tetNodes(7,0) =  0.0;  tetNodes(7,1) =  0.0;  tetNodes(7,2) =  0.5;
00114   tetNodes(8,0) =  0.5;  tetNodes(8,1) =  0.0;  tetNodes(8,2) =  0.5;
00115   tetNodes(9,0) =  0.0;  tetNodes(9,1) =  0.5;  tetNodes(9,2) =  0.5;
00116 
00117   // Generic array for the output values; needs to be properly resized depending on the operator type
00118   FieldContainer<double> vals;
00119 
00120   try{
00121     // exception #1: GRAD cannot be applied to HDIV functions 
00122     // resize vals to rank-3 container with dimensions (num. basis functions, num. points, arbitrary)
00123     vals.resize(tetBasis.getCardinality(), tetNodes.dimension(0), 3 );
00124     INTREPID_TEST_COMMAND( tetBasis.getValues(vals, tetNodes, OPERATOR_GRAD), throwCounter, nException );
00125 
00126     // exception #2: CURL cannot be applied to HDIV functions
00127     INTREPID_TEST_COMMAND( tetBasis.getValues(vals, tetNodes, OPERATOR_CURL), throwCounter, nException );
00128         
00129     // Exceptions 3-7: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and 
00130     // getDofTag() to access invalid array elements thereby causing bounds check exception
00131     // exception #3
00132     INTREPID_TEST_COMMAND( tetBasis.getDofOrdinal(3,0,0), throwCounter, nException );
00133     // exception #4
00134     INTREPID_TEST_COMMAND( tetBasis.getDofOrdinal(1,1,1), throwCounter, nException );
00135     // exception #5
00136     INTREPID_TEST_COMMAND( tetBasis.getDofOrdinal(0,4,1), throwCounter, nException );
00137     // exception #6
00138     INTREPID_TEST_COMMAND( tetBasis.getDofTag(7), throwCounter, nException );
00139     // exception #7
00140     INTREPID_TEST_COMMAND( tetBasis.getDofTag(-1), throwCounter, nException );
00141 
00142 #ifdef HAVE_INTREPID_DEBUG
00143     // Exceptions 8-15 test exception handling with incorrectly dimensioned input/output arrays
00144     // exception #8: input points array must be of rank-2
00145     FieldContainer<double> badPoints1(4, 5, 3);
00146     INTREPID_TEST_COMMAND( tetBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
00147     
00148     // exception #9 dimension 1 in the input point array must equal space dimension of the cell
00149     FieldContainer<double> badPoints2(4, 2);
00150     INTREPID_TEST_COMMAND( tetBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
00151     
00152     // exception #10 output values must be of rank-3 for OPERATOR_VALUE
00153     FieldContainer<double> badVals1(4, 3);
00154     INTREPID_TEST_COMMAND( tetBasis.getValues(badVals1, tetNodes, OPERATOR_VALUE), throwCounter, nException );
00155  
00156     // exception #11 output values must be of rank-2 for OPERATOR_DIV
00157     FieldContainer<double> badVals2(4, 3, 1);
00158     INTREPID_TEST_COMMAND( tetBasis.getValues(badVals2, tetNodes, OPERATOR_VALUE), throwCounter, nException );
00159 
00160     // exception #12 incorrect 0th dimension of output array (must equal number of basis functions)
00161     FieldContainer<double> badVals3(tetBasis.getCardinality() + 1, tetNodes.dimension(0), 3);
00162     INTREPID_TEST_COMMAND( tetBasis.getValues(badVals3, tetNodes, OPERATOR_VALUE), throwCounter, nException );
00163 
00164     // exception #13 incorrect 0th dimension of output array (must equal number of basis functions)
00165     FieldContainer<double> badVals4(tetBasis.getCardinality() + 1, tetNodes.dimension(0));
00166     INTREPID_TEST_COMMAND( tetBasis.getValues(badVals4, tetNodes, OPERATOR_DIV), throwCounter, nException );
00167 
00168     // exception #14 incorrect 1st dimension of output array (must equal number of points)
00169     FieldContainer<double> badVals5(tetBasis.getCardinality(), tetNodes.dimension(0) + 1, 3);
00170     INTREPID_TEST_COMMAND( tetBasis.getValues(badVals5, tetNodes, OPERATOR_VALUE), throwCounter, nException );
00171 
00172     // exception #15 incorrect 1st dimension of output array (must equal number of points)
00173     FieldContainer<double> badVals6(tetBasis.getCardinality(), tetNodes.dimension(0) + 1);
00174     INTREPID_TEST_COMMAND( tetBasis.getValues(badVals6, tetNodes, OPERATOR_DIV), throwCounter, nException );
00175 
00176     // exception #16: incorrect 2nd dimension of output array (must equal the space dimension)
00177     FieldContainer<double> badVals7(tetBasis.getCardinality(), tetNodes.dimension(0), 4);
00178     INTREPID_TEST_COMMAND( tetBasis.getValues(badVals7, tetNodes, OPERATOR_VALUE), throwCounter, nException );
00179 #endif
00180     
00181   }
00182   catch (std::logic_error err) {
00183     *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
00184     *outStream << err.what() << '\n';
00185     *outStream << "-------------------------------------------------------------------------------" << "\n\n";
00186     errorFlag = -1000;
00187   };
00188   
00189   // Check if number of thrown exceptions matches the one we expect
00190   if (throwCounter != nException) {
00191     errorFlag++;
00192     *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00193   }
00194   
00195   *outStream \
00196     << "\n"
00197     << "===============================================================================\n"\
00198     << "| TEST 2: correctness of tag to enum and enum to tag lookups                  |\n"\
00199     << "===============================================================================\n";
00200   
00201   try{
00202     std::vector<std::vector<int> > allTags = tetBasis.getAllDofTags();
00203     
00204     // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again
00205     for (unsigned i = 0; i < allTags.size(); i++) {
00206       int bfOrd  = tetBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
00207       
00208       std::vector<int> myTag = tetBasis.getDofTag(bfOrd);
00209        if( !( (myTag[0] == allTags[i][0]) &&
00210               (myTag[1] == allTags[i][1]) &&
00211               (myTag[2] == allTags[i][2]) &&
00212               (myTag[3] == allTags[i][3]) ) ) {
00213         errorFlag++;
00214         *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00215         *outStream << " getDofOrdinal( {" 
00216           << allTags[i][0] << ", " 
00217           << allTags[i][1] << ", " 
00218           << allTags[i][2] << ", " 
00219           << allTags[i][3] << "}) = " << bfOrd <<" but \n";   
00220         *outStream << " getDofTag(" << bfOrd << ") = { "
00221           << myTag[0] << ", " 
00222           << myTag[1] << ", " 
00223           << myTag[2] << ", " 
00224           << myTag[3] << "}\n";        
00225       }
00226     }
00227     
00228     // Now do the same but loop over basis functions
00229     for( int bfOrd = 0; bfOrd < tetBasis.getCardinality(); bfOrd++) {
00230       std::vector<int> myTag  = tetBasis.getDofTag(bfOrd);
00231       int myBfOrd = tetBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
00232       if( bfOrd != myBfOrd) {
00233         errorFlag++;
00234         *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00235         *outStream << " getDofTag(" << bfOrd << ") = { "
00236           << myTag[0] << ", " 
00237           << myTag[1] << ", " 
00238           << myTag[2] << ", " 
00239           << myTag[3] << "} but getDofOrdinal({" 
00240           << myTag[0] << ", " 
00241           << myTag[1] << ", " 
00242           << myTag[2] << ", " 
00243           << myTag[3] << "} ) = " << myBfOrd << "\n";
00244       }
00245     }
00246   }
00247   catch (std::logic_error err){
00248     *outStream << err.what() << "\n\n";
00249     errorFlag = -1000;
00250   };
00251   
00252   *outStream \
00253     << "\n"
00254     << "===============================================================================\n"\
00255     << "| TEST 3: correctness of basis function values                                |\n"\
00256     << "===============================================================================\n";
00257   
00258   outStream -> precision(20);
00259   
00260   // VALUE: Correct basis values in (P,F,D) format: each row gives the 4x3 correct basis set values 
00261   // at an evaluation point. Note that getValues returns results as an (F,P,D) array.
00262   double basisValues[] = {
00263     // 4 vertices
00264     0.,-2.0,0.,    0.,0.,0.,    -2.0,0.,0.,     0.,0.,-2.0,
00265     2.0,-2.0,0.,   2.0,0.,0.,    0.,0.,0.,      2.0,0.,-2.0,
00266     0.,0.,0.,      0.,2.0,0.,   -2.0,2.0,0.,    0,2.0,-2.0,
00267     0.,-2.0,2.0,   0.,0.,2.0,   -2.0,0.,2.0,    0.,0.,0.,
00268     // 6 edge midpoints
00269     1.0,-2.0,0.,   1.0,0.,0.,    -1.0,0.,0.,     1.0,0.,-2.0,
00270     1.0,-1.0,0.,   1.0,1.0,0.,   -1.0,1.0,0.,    1.0,1.0,-2.0,
00271     0.,-1.0,0.,    0.,1.0,0.,    -2.0,1.0,0.,    0.,1.0,-2.0,
00272     0.,-2.0,1.0,   0.,0.,1.0,    -2.0,0.,1.0,    0.,0.,-1.0,
00273     1.0,-2.0,1.0,  1.0,0.,1.0,   -1.0,0.,1.0,    1.0,0.,-1.0,
00274     0.,-1.0,1.0,   0.,1.0,1.0,   -2.0,1.0,1.0,   0.,1.0,-1.0
00275     // bf0         bf1                bf2            bf3
00276   };
00277   
00278   // DIV: each row gives the 4 correct values of the divergence of the 4 basis functions
00279   double basisDivs[] = {   
00280     // 4 vertices
00281      6.0, 6.0, 6.0, 6.0,
00282      6.0, 6.0, 6.0, 6.0,
00283      6.0, 6.0, 6.0, 6.0,
00284      6.0, 6.0, 6.0, 6.0,
00285     // 6 edge midpoints
00286      6.0, 6.0, 6.0, 6.0,
00287      6.0, 6.0, 6.0, 6.0,
00288      6.0, 6.0, 6.0, 6.0,
00289      6.0, 6.0, 6.0, 6.0,
00290      6.0, 6.0, 6.0, 6.0,
00291      6.0, 6.0, 6.0, 6.0
00292   };
00293   
00294   try{
00295         
00296     // Dimensions for the output arrays:
00297     int numFields = tetBasis.getCardinality();
00298     int numPoints = tetNodes.dimension(0);
00299     int spaceDim  = tetBasis.getBaseCellTopology().getDimension();
00300     
00301     // Generic array for values and curls that will be properly sized before each call
00302     FieldContainer<double> vals;
00303     
00304     // Check VALUE of basis functions: resize vals to rank-3 container:
00305     vals.resize(numFields, numPoints, spaceDim);
00306     tetBasis.getValues(vals, tetNodes, OPERATOR_VALUE);
00307     for (int i = 0; i < numFields; i++) {
00308       for (int j = 0; j < numPoints; j++) {
00309         for (int k = 0; k < spaceDim; k++) {
00310           // basisValues is (P,F,D) array so its multiindex is (j,i,k) and not (i,j,k)!
00311            int l = k + i * spaceDim + j * spaceDim * numFields;
00312            if (std::abs(vals(i,j,k) - basisValues[l]) > INTREPID_TOL) {
00313              errorFlag++;
00314              *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00315 
00316              // Output the multi-index of the value where the error is:
00317              *outStream << " At (Field,Point,Dim) multi-index { ";
00318              *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
00319              *outStream << "}  computed value: " << vals(i,j,k)
00320                << " but reference value: " << basisValues[l] << "\n";
00321             }
00322          }
00323       }
00324     }
00325 
00326     
00327     // Check DIV of basis function: resize vals to rank-2 container
00328     vals.resize(numFields, numPoints);
00329     tetBasis.getValues(vals, tetNodes, OPERATOR_DIV);
00330     for (int i = 0; i < numFields; i++) {
00331       for (int j = 0; j < numPoints; j++) {
00332           int l =  i + j * numFields;
00333            if (std::abs(vals(i,j) - basisDivs[l]) > INTREPID_TOL) {
00334              errorFlag++;
00335              *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00336 
00337              // Output the multi-index of the value where the error is:
00338              *outStream << " At multi-index { ";
00339              *outStream << i << " ";*outStream << j << " ";
00340              *outStream << "}  computed divergence component: " << vals(i,j)
00341                << " but reference divergence component: " << basisDivs[l] << "\n";
00342          }
00343       }
00344     }
00345     
00346    }    
00347   
00348   // Catch unexpected errors
00349   catch (std::logic_error err) {
00350     *outStream << err.what() << "\n\n";
00351     errorFlag = -1000;
00352   };
00353   
00354   if (errorFlag != 0)
00355     std::cout << "End Result: TEST FAILED\n";
00356   else
00357     std::cout << "End Result: TEST PASSED\n";
00358   
00359   // reset format state of std::cout
00360   std::cout.copyfmt(oldFormatState);
00361   
00362   return errorFlag;
00363 }
http://trilinos.sandia.gov/packages/docs/r10.0/packages/intrepid/test/Discretization/Basis/HDIV_TET_I1_FEM/test_01.cpp
