Intrepid
http://trilinos.sandia.gov/packages/docs/r10.12/packages/intrepid/test/Discretization/Basis/HGRAD_TET_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_TET_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_TET_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_TET_C1_FEM<double, FieldContainer<double> > tetBasis;
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 4 vertices of the reference TET and its center.  
00117   FieldContainer<double> tetNodes(8, 3);
00118   tetNodes(0,0) =  0.0;  tetNodes(0,1) =  0.0;  tetNodes(0,2) =  0.0;  
00119   tetNodes(1,0) =  1.0;  tetNodes(1,1) =  0.0;  tetNodes(1,2) =  0.0;  
00120   tetNodes(2,0) =  0.0;  tetNodes(2,1) =  1.0;  tetNodes(2,2) =  0.0;
00121   tetNodes(3,0) =  0.0;  tetNodes(3,1) =  0.0;  tetNodes(3,2) =  1.0;  
00122   tetNodes(4,0) =  0.25; tetNodes(4,1) =  0.5;  tetNodes(4,2) =  0.1;
00123   tetNodes(5,0) =  0.5;  tetNodes(5,1) =  0.5;  tetNodes(5,2) =  0.0;  
00124   tetNodes(6,0) =  0.5;  tetNodes(6,1) =  0.0;  tetNodes(6,2) =  0.5;  
00125   tetNodes(7,0) =  0.0;  tetNodes(7,1) =  0.5;  tetNodes(7,2) =  0.5;  
00126 
00127 
00128   // Generic array for the output values; needs to be properly resized depending on the operator type
00129   FieldContainer<double> vals;
00130 
00131   try{
00132     // exception #1: CURL cannot be applied to scalar functions
00133     // resize vals to rank-3 container with dimensions (num. points, num. basis functions, arbitrary)
00134     vals.resize(tetBasis.getCardinality(), tetNodes.dimension(0), 4 );
00135     INTREPID_TEST_COMMAND( tetBasis.getValues(vals, tetNodes, OPERATOR_CURL), throwCounter, nException );
00136 
00137     // exception #2: DIV cannot be applied to scalar functions
00138     // resize vals to rank-2 container with dimensions (num. points, num. basis functions)
00139     vals.resize(tetBasis.getCardinality(), tetNodes.dimension(0) );
00140     INTREPID_TEST_COMMAND( tetBasis.getValues(vals, tetNodes, OPERATOR_DIV), 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,0), throwCounter, nException );
00150     // exception #6
00151     INTREPID_TEST_COMMAND( tetBasis.getDofTag(5), throwCounter, nException );
00152     // exception #7
00153     INTREPID_TEST_COMMAND( tetBasis.getDofTag(-1), throwCounter, nException );
00154     
00155 #ifdef HAVE_INTREPID_DEBUG 
00156     // Exceptions 8-18 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, tetBasis.getBaseCellTopology().getDimension() - 1);
00163     INTREPID_TEST_COMMAND( tetBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
00164     
00165     // exception #10 output values must be of rank-2 for OPERATOR_VALUE
00166     FieldContainer<double> badVals1(4, 3, 1);
00167     INTREPID_TEST_COMMAND( tetBasis.getValues(badVals1, tetNodes, OPERATOR_VALUE), throwCounter, nException );
00168     
00169     // exception #11 output values must be of rank-3 for OPERATOR_GRAD
00170     FieldContainer<double> badVals2(4, 3);
00171     INTREPID_TEST_COMMAND( tetBasis.getValues(badVals2, tetNodes, OPERATOR_GRAD), throwCounter, nException );
00172     
00173     // exception #12 output values must be of rank-3 for OPERATOR_D1
00174     INTREPID_TEST_COMMAND( tetBasis.getValues(badVals2, tetNodes, OPERATOR_D1), throwCounter, nException );
00175     
00176     // exception #13 output values must be of rank-3 for OPERATOR_D2
00177     INTREPID_TEST_COMMAND( tetBasis.getValues(badVals2, tetNodes, OPERATOR_D2), throwCounter, nException );
00178     
00179     // exception #14 incorrect 0th dimension of output array (must equal number of basis functions)
00180     FieldContainer<double> badVals3(tetBasis.getCardinality() + 1, tetNodes.dimension(0));
00181     INTREPID_TEST_COMMAND( tetBasis.getValues(badVals3, tetNodes, OPERATOR_VALUE), throwCounter, nException );
00182     
00183     // exception #15 incorrect 1st dimension of output array (must equal number of points)
00184     FieldContainer<double> badVals4(tetBasis.getCardinality(), tetNodes.dimension(0) + 1);
00185     INTREPID_TEST_COMMAND( tetBasis.getValues(badVals4, tetNodes, OPERATOR_VALUE), throwCounter, nException );
00186     
00187     // exception #16: incorrect 2nd dimension of output array (must equal the space dimension)
00188     FieldContainer<double> badVals5(tetBasis.getCardinality(), tetNodes.dimension(0), tetBasis.getBaseCellTopology().getDimension() + 1);
00189     INTREPID_TEST_COMMAND( tetBasis.getValues(badVals5, tetNodes, OPERATOR_GRAD), throwCounter, nException );
00190     
00191     // exception #17: incorrect 2nd dimension of output array (must equal D2 cardinality in 2D)
00192     FieldContainer<double> badVals6(tetBasis.getCardinality(), tetNodes.dimension(0), 40);
00193     INTREPID_TEST_COMMAND( tetBasis.getValues(badVals6, tetNodes, OPERATOR_D1), throwCounter, nException );
00194     
00195     // exception #18: incorrect 2nd dimension of output array (must equal D3 cardinality in 2D)
00196     INTREPID_TEST_COMMAND( tetBasis.getValues(badVals6, tetNodes, OPERATOR_D2), throwCounter, nException );
00197 #endif
00198     
00199   }
00200   catch (std::logic_error err) {
00201     *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
00202     *outStream << err.what() << '\n';
00203     *outStream << "-------------------------------------------------------------------------------" << "\n\n";
00204     errorFlag = -1000;
00205   };
00206   
00207   // Check if number of thrown exceptions matches the one we expect 
00208   if (throwCounter != nException) {
00209     errorFlag++;
00210     *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00211   }
00212   
00213   *outStream \
00214     << "\n"
00215     << "===============================================================================\n"\
00216     << "| TEST 2: correctness of tag to enum and enum to tag lookups                  |\n"\
00217     << "===============================================================================\n";
00218   
00219   try{
00220     std::vector<std::vector<int> > allTags = tetBasis.getAllDofTags();
00221     
00222     // Loop over all tags, lookup the associated dof enumeration and then lookup the tag again
00223     for (unsigned i = 0; i < allTags.size(); i++) {
00224       int bfOrd  = tetBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
00225       
00226       std::vector<int> myTag = tetBasis.getDofTag(bfOrd);
00227        if( !( (myTag[0] == allTags[i][0]) &&
00228               (myTag[1] == allTags[i][1]) &&
00229               (myTag[2] == allTags[i][2]) &&
00230               (myTag[3] == allTags[i][3]) ) ) {
00231         errorFlag++;
00232         *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00233         *outStream << " getDofOrdinal( {" 
00234           << allTags[i][0] << ", " 
00235           << allTags[i][1] << ", " 
00236           << allTags[i][2] << ", " 
00237           << allTags[i][3] << "}) = " << bfOrd <<" but \n";   
00238         *outStream << " getDofTag(" << bfOrd << ") = { "
00239           << myTag[0] << ", " 
00240           << myTag[1] << ", " 
00241           << myTag[2] << ", " 
00242           << myTag[3] << "}\n";        
00243       }
00244     }
00245     
00246     // Now do the same but loop over basis functions
00247     for( int bfOrd = 0; bfOrd < tetBasis.getCardinality(); bfOrd++) {
00248       std::vector<int> myTag  = tetBasis.getDofTag(bfOrd);
00249       int myBfOrd = tetBasis.getDofOrdinal(myTag[0], myTag[1], myTag[2]);
00250       if( bfOrd != myBfOrd) {
00251         errorFlag++;
00252         *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00253         *outStream << " getDofTag(" << bfOrd << ") = { "
00254           << myTag[0] << ", " 
00255           << myTag[1] << ", " 
00256           << myTag[2] << ", " 
00257           << myTag[3] << "} but getDofOrdinal({" 
00258           << myTag[0] << ", " 
00259           << myTag[1] << ", " 
00260           << myTag[2] << ", " 
00261           << myTag[3] << "} ) = " << myBfOrd << "\n";
00262       }
00263     }
00264   }
00265   catch (std::logic_error err){
00266     *outStream << err.what() << "\n\n";
00267     errorFlag = -1000;
00268   };
00269   
00270   *outStream \
00271     << "\n"
00272     << "===============================================================================\n"\
00273     << "| TEST 3: correctness of basis function values                                |\n"\
00274     << "===============================================================================\n";
00275   
00276   outStream -> precision(20);
00277   
00278   // VALUE: Each row gives the 4 correct basis set values at an evaluation point
00279   double basisValues[] = {
00280     1.0, 0.0, 0.0, 0.0,
00281     0.0, 1.0, 0.0, 0.0,
00282     0.0, 0.0, 1.0, 0.0,
00283     0.0, 0.0, 0.0, 1.0,
00284     0.15,0.25,0.5, 0.1,
00285     0.0, 0.5, 0.5, 0.0,
00286     0.0, 0.5, 0.0, 0.5,
00287     0.0, 0.0, 0.5, 0.5
00288   };
00289   
00290   // GRAD and D1: each row gives the 3 x 4 correct values of the gradients of the 4 basis functions
00291   double basisGrads[] = {
00292     -1.0, -1.0, -1.0,    1.0, 0.0, 0.0,    0.0, 1.0, 0.0,    0.0, 0.0, 1.0, 
00293     -1.0, -1.0, -1.0,    1.0, 0.0, 0.0,    0.0, 1.0, 0.0,    0.0, 0.0, 1.0, 
00294     -1.0, -1.0, -1.0,    1.0, 0.0, 0.0,    0.0, 1.0, 0.0,    0.0, 0.0, 1.0, 
00295     -1.0, -1.0, -1.0,    1.0, 0.0, 0.0,    0.0, 1.0, 0.0,    0.0, 0.0, 1.0, 
00296     -1.0, -1.0, -1.0,    1.0, 0.0, 0.0,    0.0, 1.0, 0.0,    0.0, 0.0, 1.0, 
00297     -1.0, -1.0, -1.0,    1.0, 0.0, 0.0,    0.0, 1.0, 0.0,    0.0, 0.0, 1.0, 
00298     -1.0, -1.0, -1.0,    1.0, 0.0, 0.0,    0.0, 1.0, 0.0,    0.0, 0.0, 1.0, 
00299     -1.0, -1.0, -1.0,    1.0, 0.0, 0.0,    0.0, 1.0, 0.0,    0.0, 0.0, 1.0, 
00300   };
00301   
00302   try{
00303         
00304     // Dimensions for the output arrays:
00305     int numFields = tetBasis.getCardinality();
00306     int numPoints = tetNodes.dimension(0);
00307     int spaceDim  = tetBasis.getBaseCellTopology().getDimension();
00308     
00309     // Generic array for values, grads, curls, etc. that will be properly sized before each call
00310     FieldContainer<double> vals;
00311     
00312     // Check VALUE of basis functions: resize vals to rank-2 container:
00313     vals.resize(numFields, numPoints);
00314     tetBasis.getValues(vals, tetNodes, OPERATOR_VALUE);
00315     for (int i = 0; i < numFields; i++) {
00316       for (int j = 0; j < numPoints; j++) {
00317           int l =  i + j * numFields;
00318            if (std::abs(vals(i,j) - basisValues[l]) > INTREPID_TOL) {
00319              errorFlag++;
00320              *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00321 
00322              // Output the multi-index of the value where the error is:
00323              *outStream << " At multi-index { ";
00324              *outStream << i << " ";*outStream << j << " ";
00325              *outStream << "}  computed value: " << vals(i,j)
00326                << " but reference value: " << basisValues[l] << "\n";
00327          }
00328       }
00329     }
00330 
00331     // Check GRAD of basis function: resize vals to rank-3 container
00332     vals.resize(numFields, numPoints, spaceDim);
00333     tetBasis.getValues(vals, tetNodes, OPERATOR_GRAD);
00334     for (int i = 0; i < numFields; i++) {
00335       for (int j = 0; j < numPoints; j++) {
00336         for (int k = 0; k < spaceDim; k++) {
00337            int l = k + i * spaceDim + j * spaceDim * numFields;
00338            if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
00339              errorFlag++;
00340              *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00341 
00342              // Output the multi-index of the value where the error is:
00343              *outStream << " At multi-index { ";
00344              *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
00345              *outStream << "}  computed grad component: " << vals(i,j,k)
00346                << " but reference grad component: " << basisGrads[l] << "\n";
00347             }
00348          }
00349       }
00350     }
00351 
00352     // Check D1 of basis function (do not resize vals because it has the correct size: D1 = GRAD)
00353     tetBasis.getValues(vals, tetNodes, OPERATOR_D1);
00354     for (int i = 0; i < numFields; i++) {
00355       for (int j = 0; j < numPoints; j++) {
00356         for (int k = 0; k < spaceDim; k++) {
00357            int l = k + i * spaceDim + j * spaceDim * numFields;
00358            if (std::abs(vals(i,j,k) - basisGrads[l]) > INTREPID_TOL) {
00359              errorFlag++;
00360              *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00361 
00362              // Output the multi-index of the value where the error is:
00363              *outStream << " At multi-index { ";
00364              *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
00365              *outStream << "}  computed D1 component: " << vals(i,j,k)
00366                << " but reference D1 component: " << basisGrads[l] << "\n";
00367             }
00368          }
00369       }
00370     }
00371 
00372     // Check all higher derivatives - must be zero. 
00373     for(EOperator op = OPERATOR_D2; op < OPERATOR_MAX; op++) {
00374       
00375       // The last dimension is the number of kth derivatives and needs to be resized for every Dk
00376       int DkCardin  = Intrepid::getDkCardinality(op, spaceDim);
00377       vals.resize(numFields, numPoints, DkCardin);    
00378 
00379       tetBasis.getValues(vals, tetNodes, op);
00380       for (int i = 0; i < vals.size(); i++) {
00381         if (std::abs(vals[i]) > INTREPID_TOL) {
00382           errorFlag++;
00383           *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00384           
00385           // Get the multi-index of the value where the error is and the operator order
00386           std::vector<int> myIndex;
00387           vals.getMultiIndex(myIndex,i);
00388           int ord = Intrepid::getOperatorOrder(op);
00389           *outStream << " At multi-index { ";
00390           for(int j = 0; j < vals.rank(); j++) {
00391             *outStream << myIndex[j] << " ";
00392           }
00393           *outStream << "}  computed D"<< ord <<" component: " << vals[i] 
00394             << " but reference D" << ord << " component:  0 \n";
00395         }
00396       }
00397     }    
00398   }
00399   
00400   // Catch unexpected errors
00401   catch (std::logic_error err) {
00402     *outStream << err.what() << "\n\n";
00403     errorFlag = -1000;
00404   };
00405   
00406   if (errorFlag != 0)
00407     std::cout << "End Result: TEST FAILED\n";
00408   else
00409     std::cout << "End Result: TEST PASSED\n";
00410   
00411   // reset format state of std::cout
00412   std::cout.copyfmt(oldFormatState);
00413   
00414   return errorFlag;
00415 }