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