Intrepid
http://trilinos.sandia.gov/packages/docs/r10.12/packages/intrepid/test/Discretization/Basis/HDIV_WEDGE_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_WEDGE_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_WEDGE_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_WEDGE_I1_FEM<double, FieldContainer<double> > wedgeBasis;
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 6 vertices of the reference WEDGE and 6 other points.
00117   FieldContainer<double> wedgeNodes(12, 3);
00118   wedgeNodes(0,0) =  0.0;  wedgeNodes(0,1) =  0.0;  wedgeNodes(0,2) = -1.0;
00119   wedgeNodes(1,0) =  1.0;  wedgeNodes(1,1) =  0.0;  wedgeNodes(1,2) = -1.0;
00120   wedgeNodes(2,0) =  0.0;  wedgeNodes(2,1) =  1.0;  wedgeNodes(2,2) = -1.0;
00121   wedgeNodes(3,0) =  0.0;  wedgeNodes(3,1) =  0.0;  wedgeNodes(3,2) =  1.0;
00122   wedgeNodes(4,0) =  1.0;  wedgeNodes(4,1) =  0.0;  wedgeNodes(4,2) =  1.0;
00123   wedgeNodes(5,0) =  0.0;  wedgeNodes(5,1) =  1.0;  wedgeNodes(5,2) =  1.0;
00124 
00125   wedgeNodes(6,0) =  0.25; wedgeNodes(6,1) =  0.5;  wedgeNodes(6,2) = -1.0;
00126   wedgeNodes(7,0) =  0.5;  wedgeNodes(7,1) =  0.25; wedgeNodes(7,2) =  0.0;
00127   wedgeNodes(8,0) =  0.25; wedgeNodes(8,1) =  0.25; wedgeNodes(8,2) =  1.0;
00128   wedgeNodes(9,0) =  0.25; wedgeNodes(9,1) =  0.0;  wedgeNodes(9,2) =  0.75;
00129   wedgeNodes(10,0)=  0.0;  wedgeNodes(10,1)=  0.5;  wedgeNodes(10,2)= -0.25;
00130   wedgeNodes(11,0)=  0.5;  wedgeNodes(11,1)=  0.5;  wedgeNodes(11,2)=  0.0;
00131 
00132   // Generic array for the output values; needs to be properly resized depending on the operator type
00133   FieldContainer<double> vals;
00134 
00135   try{
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(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), 3 );
00139     INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_GRAD), throwCounter, nException );
00140 
00141     // exception #2: CURL cannot be applied to HDIV functions
00142     INTREPID_TEST_COMMAND( wedgeBasis.getValues(vals, wedgeNodes, 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( wedgeBasis.getDofOrdinal(3,0,0), throwCounter, nException );
00148     // exception #4
00149     INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(1,1,1), throwCounter, nException );
00150     // exception #5
00151     INTREPID_TEST_COMMAND( wedgeBasis.getDofOrdinal(0,4,1), throwCounter, nException );
00152     // exception #6
00153     INTREPID_TEST_COMMAND( wedgeBasis.getDofTag(11), throwCounter, nException );
00154     // exception #7
00155     INTREPID_TEST_COMMAND( wedgeBasis.getDofTag(-1), throwCounter, nException );
00156 
00157 #ifdef HAVE_INTREPID_DEBUG
00158     // Exceptions 8-15 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( wedgeBasis.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, 2);
00165     INTREPID_TEST_COMMAND( wedgeBasis.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, 3);
00169     INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals1, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
00170  
00171     // exception #11 output values must be of rank-2 for OPERATOR_DIV
00172     FieldContainer<double> badVals2(4, 3, 1);
00173     INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals2, wedgeNodes, OPERATOR_DIV), throwCounter, nException );
00174     
00175     // exception #12 incorrect 0th dimension of output array (must equal number of basis functions)
00176     FieldContainer<double> badVals3(wedgeBasis.getCardinality() + 1, wedgeNodes.dimension(0), 3);
00177     INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals3, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
00178 
00179     // exception #13 incorrect 0th dimension of output array (must equal number of basis functions)
00180     FieldContainer<double> badVals4(wedgeBasis.getCardinality() + 1, wedgeNodes.dimension(0));
00181     INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals4, wedgeNodes, OPERATOR_DIV), throwCounter, nException );
00182 
00183     // exception #14 incorrect 1st dimension of output array (must equal number of points)
00184     FieldContainer<double> badVals5(wedgeBasis.getCardinality(), wedgeNodes.dimension(0) + 1, 3);
00185     INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals5, wedgeNodes, OPERATOR_VALUE), throwCounter, nException );
00186 
00187     // exception #15 incorrect 1st dimension of output array (must equal number of points)
00188     FieldContainer<double> badVals6(wedgeBasis.getCardinality(), wedgeNodes.dimension(0) + 1);
00189     INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals6, wedgeNodes, OPERATOR_DIV), throwCounter, nException );
00190 
00191     // exception #16: incorrect 2nd dimension of output array (must equal the space dimension)
00192     FieldContainer<double> badVals7(wedgeBasis.getCardinality(), wedgeNodes.dimension(0), 4);
00193     INTREPID_TEST_COMMAND( wedgeBasis.getValues(badVals7, wedgeNodes, 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 = wedgeBasis.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  = wedgeBasis.getDofOrdinal(allTags[i][0], allTags[i][1], allTags[i][2]);
00222       
00223       std::vector<int> myTag = wedgeBasis.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 < wedgeBasis.getCardinality(); bfOrd++) {
00245       std::vector<int> myTag  = wedgeBasis.getDofTag(bfOrd);
00246       int myBfOrd = wedgeBasis.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 5x3 correct basis set values at an evaluation point
00276   double basisValues[] = {
00277     0, -0.500000, 0, 0, 0, 0, -0.500000, 0, 0, 0, 0, -2.00000, 0, 0, 0, \
00278     0.500000, -0.500000, 0, 0.500000, 0, 0, 0, 0, 0, 0, 0, -2.00000, 0, \
00279     0, 0, 0, 0, 0, 0, 0.500000, 0, -0.500000, 0.500000, 0, 0, 0, \
00280     -2.00000, 0, 0, 0, 0, -0.500000, 0, 0, 0, 0, -0.500000, 0, 0, 0, 0, \
00281     0, 0, 0, 2.00000, 0.500000, -0.500000, 0, 0.500000, 0, 0, 0, 0, 0, 0, \
00282     0, 0, 0, 0, 2.00000, 0, 0, 0, 0, 0.500000, 0, -0.500000, 0.500000, 0, \
00283     0, 0, 0, 0, 0, 2.00000, 0.125000, -0.250000, 0, 0.125000, 0.250000, \
00284     0, -0.375000, 0.250000, 0, 0, 0, -2.00000, 0, 0, 0, 0.250000, \
00285     -0.375000, 0, 0.250000, 0.125000, 0, -0.250000, 0.125000, 0, 0, 0, \
00286     -1.00000, 0, 0, 1.00000, 0.125000, -0.375000, 0, 0.125000, 0.125000, \
00287     0, -0.375000, 0.125000, 0, 0, 0, 0, 0, 0, 2.00000, 0.125000, \
00288     -0.500000, 0, 0.125000, 0, 0, -0.375000, 0, 0, 0, 0, -0.250000, 0, 0, \
00289     1.75000, 0, -0.250000, 0, 0, 0.250000, 0, -0.500000, 0.250000, 0, 0, \
00290     0, -1.25000, 0, 0, 0.750000, 0.250000, -0.250000, 0, 0.250000, \
00291     0.250000, 0, -0.250000, 0.250000, 0, 0, 0, -1.00000, 0, 0, 1.00000};
00292   
00293   // DIV: each row pair gives the 5 correct values of the divergence of the 5 basis functions
00294   double basisDivs[] = {   
00295     // 6 vertices
00296     1.0, 1.0, 1.0, 1.0, 1.0,
00297     1.0, 1.0, 1.0, 1.0, 1.0,
00298     1.0, 1.0, 1.0, 1.0, 1.0,
00299     1.0, 1.0, 1.0, 1.0, 1.0,
00300     1.0, 1.0, 1.0, 1.0, 1.0,
00301     1.0, 1.0, 1.0, 1.0, 1.0,
00302     // 6 other points
00303     1.0, 1.0, 1.0, 1.0, 1.0,
00304     1.0, 1.0, 1.0, 1.0, 1.0,
00305     1.0, 1.0, 1.0, 1.0, 1.0,
00306     1.0, 1.0, 1.0, 1.0, 1.0,
00307     1.0, 1.0, 1.0, 1.0, 1.0,
00308     1.0, 1.0, 1.0, 1.0, 1.0
00309   };
00310   
00311   try{
00312         
00313     // Dimensions for the output arrays:
00314     int numFields = wedgeBasis.getCardinality();
00315     int numPoints = wedgeNodes.dimension(0);
00316     int spaceDim  = wedgeBasis.getBaseCellTopology().getDimension();
00317     
00318     // Generic array for values and curls that will be properly sized before each call
00319     FieldContainer<double> vals;
00320     
00321     // Check VALUE of basis functions: resize vals to rank-3 container:
00322     vals.resize(numFields, numPoints, spaceDim);
00323     wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_VALUE);
00324     for (int i = 0; i < numFields; i++) {
00325       for (int j = 0; j < numPoints; j++) {
00326         for (int k = 0; k < spaceDim; k++) {
00327            int l = k + i * spaceDim + j * spaceDim * numFields;
00328            if (std::abs(vals(i,j,k) - basisValues[l]) > INTREPID_TOL) {
00329              errorFlag++;
00330              *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00331 
00332              // Output the multi-index of the value where the error is:
00333              *outStream << " At multi-index { ";
00334              *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
00335              *outStream << "}  computed value: " << vals(i,j,k)
00336                << " but reference value: " << basisValues[l] << "\n";
00337             }
00338          }
00339       }
00340     }
00341 
00342     
00343     // Check DIV of basis function: resize vals to rank-2 container
00344     vals.resize(numFields, numPoints);
00345     wedgeBasis.getValues(vals, wedgeNodes, OPERATOR_DIV);
00346     for (int i = 0; i < numFields; i++) {
00347       for (int j = 0; j < numPoints; j++) {
00348           int l =  i + j * numFields;
00349            if (std::abs(vals(i,j) - basisDivs[l]) > INTREPID_TOL) {
00350              errorFlag++;
00351              *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00352 
00353              // Output the multi-index of the value where the error is:
00354              *outStream << " At multi-index { ";
00355              *outStream << i << " ";*outStream << j << " ";
00356              *outStream << "}  computed divergence component: " << vals(i,j)
00357                << " but reference divergence component: " << basisDivs[l] << "\n";
00358          }
00359       }
00360     }
00361     
00362    }    
00363   
00364   // Catch unexpected errors
00365   catch (std::logic_error err) {
00366     *outStream << err.what() << "\n\n";
00367     errorFlag = -1000;
00368   };
00369   
00370   if (errorFlag != 0)
00371     std::cout << "End Result: TEST FAILED\n";
00372   else
00373     std::cout << "End Result: TEST PASSED\n";
00374   
00375   // reset format state of std::cout
00376   std::cout.copyfmt(oldFormatState);
00377   
00378   return errorFlag;
00379 }