Intrepid
http://trilinos.sandia.gov/packages/docs/r10.10/packages/intrepid/test/Discretization/Basis/HCURL_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_HCURL_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_HCURL_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 CURL 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_HCURL_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     // exception #1: GRAD cannot be applied to HCURL functions 
00134     // resize vals to rank-3 container with dimensions (num. basis functions, num. points, arbitrary)
00135     vals.resize(quadBasis.getCardinality(), quadNodes.dimension(0), 4 );
00136     INTREPID_TEST_COMMAND( quadBasis.getValues(vals, quadNodes, OPERATOR_GRAD), throwCounter, nException );
00137 
00138     // exception #2: DIV cannot be applied to HCURL functions
00139     // resize vals to rank-2 container with dimensions (num. points, num. basis functions)
00140     vals.resize(quadBasis.getCardinality(), quadNodes.dimension(0) );
00141     INTREPID_TEST_COMMAND( quadBasis.getValues(vals, quadNodes, OPERATOR_DIV), throwCounter, nException );
00142         
00143     // Exceptions 3-7: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and 
00144     // getDofTag() to access invalid array elements thereby causing bounds check exception
00145     // exception #3
00146     INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(3,0,0), throwCounter, nException );
00147     // exception #4
00148     INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(1,1,1), throwCounter, nException );
00149     // exception #5
00150     INTREPID_TEST_COMMAND( quadBasis.getDofOrdinal(0,4,1), throwCounter, nException );
00151     // exception #6
00152     INTREPID_TEST_COMMAND( quadBasis.getDofTag(12), throwCounter, nException );
00153     // exception #7
00154     INTREPID_TEST_COMMAND( quadBasis.getDofTag(-1), throwCounter, nException );
00155     
00156 #ifdef HAVE_INTREPID_DEBUG
00157     // Exceptions 8-15 test exception handling with incorrectly dimensioned input/output arrays
00158     // exception #8: input points array must be of rank-2
00159     FieldContainer<double> badPoints1(4, 5, 3);
00160     INTREPID_TEST_COMMAND( quadBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
00161     
00162     // exception #9 dimension 1 in the input point array must equal space dimension of the cell
00163     FieldContainer<double> badPoints2(4, quadBasis.getBaseCellTopology().getDimension() + 1);
00164     INTREPID_TEST_COMMAND( quadBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
00165     
00166     // exception #10 output values must be of rank-3 for OPERATOR_VALUE in 2D
00167     FieldContainer<double> badVals1(4, 3);
00168     INTREPID_TEST_COMMAND( quadBasis.getValues(badVals1, quadNodes, OPERATOR_VALUE), throwCounter, nException );
00169  
00170     FieldContainer<double> badCurls1(4,3,2);
00171     // exception #11 output values must be of rank-2 for OPERATOR_CURL
00172     INTREPID_TEST_COMMAND( quadBasis.getValues(badCurls1, quadNodes, OPERATOR_CURL), throwCounter, nException );
00173     
00174     // exception #12 incorrect 0th dimension of output array (must equal number of basis functions)
00175     FieldContainer<double> badVals2(quadBasis.getCardinality() + 1, quadNodes.dimension(0), quadBasis.getBaseCellTopology().getDimension());
00176     INTREPID_TEST_COMMAND( quadBasis.getValues(badVals2, quadNodes, OPERATOR_VALUE), throwCounter, nException ) ;
00177     
00178     // exception #13 incorrect 1st  dimension of output array (must equal number of points)
00179     FieldContainer<double> badVals3(quadBasis.getCardinality(), quadNodes.dimension(0) + 1, quadBasis.getBaseCellTopology().getDimension() );
00180     INTREPID_TEST_COMMAND( quadBasis.getValues(badVals3, quadNodes, OPERATOR_VALUE), throwCounter, nException ) ;
00181 
00182     // exception #14: incorrect 2nd dimension of output array for VALUE (must equal the space dimension)
00183     FieldContainer<double> badVals4(quadBasis.getCardinality(), quadNodes.dimension(0), quadBasis.getBaseCellTopology().getDimension() - 1);
00184     INTREPID_TEST_COMMAND( quadBasis.getValues(badVals4, quadNodes, OPERATOR_VALUE), throwCounter, nException ) ;
00185         
00186     // exception #15: D2 cannot be applied to HCURL functions 
00187     // resize vals to rank-3 container with dimensions (num. basis functions, num. points, arbitrary)
00188     vals.resize(quadBasis.getCardinality(), 
00189                 quadNodes.dimension(0),  
00190                 Intrepid::getDkCardinality(OPERATOR_D2, quadBasis.getBaseCellTopology().getDimension()));
00191     INTREPID_TEST_COMMAND( quadBasis.getValues(vals, quadNodes, OPERATOR_D2), 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   // Note Teuchos throw number will not pick up exceptions 3-7 and therefore will not match.
00204   if (throwCounter != nException) {
00205     errorFlag++;
00206     *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00207   }
00208 //#endif
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 4x2 correct basis set values at an evaluation point: (P,F,D) layout
00276   double basisValues[] = {
00277     0.500000, 0, 0, 0, 0, 0, 0, -0.500000, 0.500000, 0, 0, 0.500000, 0, \
00278     0, 0, 0, 0, 0, 0, 0.500000, -0.500000, 0, 0, 0, 0, 0, 0, 0, \
00279     -0.500000, 0, 0, -0.500000, 0.250000, 0, 0, 0.250000, -0.250000, 0, \
00280     0, -0.250000, 0.375000, 0, 0, 0.250000, -0.125000, 0, 0, -0.250000, \
00281     0.125000, 0, 0, 0.250000, -0.375000, 0, 0, -0.250000, 0.250000, 0, 0, \
00282     0.125000, -0.250000, 0, 0, -0.375000, 0.250000, 0, 0, 0.375000, \
00283     -0.250000, 0, 0, -0.125000
00284   };
00285   
00286   // CURL: correct values in (F,P) format
00287   double basisCurls[] = {
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     0.25, 0.25, 0.25, 0.25
00297   };
00298 
00299   
00300   try{
00301         
00302     // Dimensions for the output arrays:
00303     int numFields = quadBasis.getCardinality();
00304     int numPoints = quadNodes.dimension(0);
00305     int spaceDim  = quadBasis.getBaseCellTopology().getDimension();
00306     
00307     // Generic array for values and curls that will be properly sized before each call
00308     FieldContainer<double> vals;
00309     
00310     // Check VALUE of basis functions: resize vals to rank-3 container:
00311     vals.resize(numFields, numPoints, spaceDim);
00312     quadBasis.getValues(vals, quadNodes, OPERATOR_VALUE);
00313     for (int i = 0; i < numFields; i++) {
00314       for (int j = 0; j < numPoints; j++) {
00315         for (int k = 0; k < spaceDim; k++) {
00316            
00317           // compute offset for (P,F,D) data layout: indices are P->j, F->i, D->k
00318           int l = k + i * spaceDim + j * spaceDim * numFields;
00319            if (std::abs(vals(i,j,k) - basisValues[l]) > INTREPID_TOL) {
00320              errorFlag++;
00321              *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00322         
00323              // Output the multi-index of the value where the error is: 
00324              *outStream << " At multi-index { ";
00325              *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
00326              *outStream << "}  computed value: " << vals(i,j,k)
00327                << " but reference value: " << basisValues[l] << "\n";
00328             }
00329          }
00330       }
00331     }
00332 
00333     // Check CURL of basis function: resize vals to rank-2 container
00334     vals.resize(numFields, numPoints);
00335     quadBasis.getValues(vals, quadNodes, OPERATOR_CURL);
00336     for (int i = 0; i < numFields; i++) {
00337       for (int j = 0; j < numPoints; j++) {
00338         int l =  i + j * numFields;
00339         if (std::abs(vals(i,j) - basisCurls[l]) > INTREPID_TOL) {
00340           errorFlag++;
00341           *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00342           
00343           // Output the multi-index of the value where the error is:
00344           *outStream << " At multi-index { ";
00345           *outStream << i << " ";*outStream << j << " ";
00346           *outStream << "}  computed curl component: " << vals(i,j)
00347             << " but reference curl component: " << basisCurls[l] << "\n";
00348         }
00349       }
00350     }  
00351    }    
00352   
00353   // Catch unexpected errors
00354   catch (std::logic_error err) {
00355     *outStream << err.what() << "\n\n";
00356     errorFlag = -1000;
00357   };
00358   
00359   if (errorFlag != 0)
00360     std::cout << "End Result: TEST FAILED\n";
00361   else
00362     std::cout << "End Result: TEST PASSED\n";
00363   
00364   // reset format state of std::cout
00365   std::cout.copyfmt(oldFormatState);
00366   
00367   return errorFlag;
00368 }