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