Intrepid
http://trilinos.sandia.gov/packages/docs/r10.12/packages/intrepid/test/Discretization/Basis/HGRAD_LINE_Cn_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 
00049 #include "Intrepid_FieldContainer.hpp"
00050 #include "Intrepid_HGRAD_LINE_Cn_FEM.hpp"
00051 #include "Intrepid_DefaultCubatureFactory.hpp"
00052 #include "Intrepid_ArrayTools.hpp"
00053 #include "Intrepid_PointTools.hpp"
00054 #include "Intrepid_FunctionSpaceTools.hpp"
00055 #include "Teuchos_oblackholestream.hpp"
00056 #include "Teuchos_RCP.hpp"
00057 #include "Teuchos_GlobalMPISession.hpp"
00058 
00059 using namespace std;
00060 using namespace Intrepid;
00061 
00062 #define INTREPID_TEST_COMMAND( S , throwCounter, nException )                                                              \
00063 {                                                                                                                          \
00064   ++nException;                                                                                                            \
00065   try {                                                                                                                    \
00066     S ;                                                                                                                    \
00067   }                                                                                                                        \
00068   catch (std::logic_error err) {                                                                                           \
00069       ++throwCounter;                                                                                                      \
00070       *outStream << "Expected Error " << nException << " -------------------------------------------------------------\n"; \
00071       *outStream << err.what() << '\n';                                                                                    \
00072       *outStream << "-------------------------------------------------------------------------------" << "\n\n";           \
00073   };                                                                                                                       \
00074 }
00075 
00076 
00077 int main(int argc, char *argv[]) {
00078 
00079   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00080 
00081   // This little trick lets us print to std::cout only if
00082   // a (dummy) command-line argument is provided.
00083   int iprint     = argc - 1;
00084   Teuchos::RCP<std::ostream> outStream;
00085   Teuchos::oblackholestream bhs; // outputs nothing
00086   if (iprint > 0)
00087     outStream = Teuchos::rcp(&std::cout, false);
00088   else
00089     outStream = Teuchos::rcp(&bhs, false);
00090 
00091   // Save the format state of the original std::cout.
00092   Teuchos::oblackholestream oldFormatState;
00093   oldFormatState.copyfmt(std::cout);
00094 
00095   *outStream \
00096     << "===============================================================================\n" \
00097     << "|                                                                             |\n" \
00098     << "|               Unit Test (Basis_HGRAD_LINE_Cn_FEM)                           |\n" \
00099     << "|                                                                             |\n" \
00100     << "|     1) Conversion of Dof tags into Dof ordinals and back                    |\n" \
00101     << "|     2) Basis values for VALUE, GRAD, CURL, and Dk operators                 |\n" \
00102     << "|                                                                             |\n" \
00103     << "|  Questions? Contact  Pavel Bochev  (pbboche@sandia.gov),                    |\n" \
00104     << "|                      Denis Ridzal  (dridzal@sandia.gov),                    |\n" \
00105     << "|                      Kara Peterson (kjpeter@sandia.gov).                    |\n" \
00106     << "|                                                                             |\n" \
00107     << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00108     << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00109     << "|                                                                             |\n" \
00110     << "===============================================================================\n"\
00111     << "| TEST 1: Basis creation, exception testing                                   |\n"\
00112     << "===============================================================================\n";
00113 
00114   
00115   // Define basis and error flag
00116   shards::CellTopology line(shards::getCellTopologyData< shards::Line<> >());   // create cell topology
00117   const int deg = 5;
00118 
00119   // get the points for the basis
00120   FieldContainer<double> pts(PointTools::getLatticeSize(line,deg),1);
00121   PointTools::getLattice<double,FieldContainer<double> >(pts,line,deg);
00122 
00123   Basis_HGRAD_LINE_Cn_FEM<double, FieldContainer<double> > lineBasis(deg, pts);
00124   int errorFlag = 0;
00125 
00126   // Initialize throw counter for exception testing
00127   int nException     = 0;
00128   int throwCounter   = 0;
00129   
00130   // Define array containing vertices of the reference Line and a few other points   
00131   int numIntervals = 100;
00132   FieldContainer<double> lineNodes(numIntervals+1, 1);
00133   for (int i=0; i<numIntervals+1; i++) {
00134     lineNodes(i,0) = -1.0+(2.0*(double)i)/(double)numIntervals;
00135   }
00136 
00137   // Generic array for the output values; needs to be properly resized depending on the operator type
00138   FieldContainer<double> vals;
00139 
00140   try{
00141     // exception #1: DIV cannot be applied to scalar functions
00142     // resize vals to rank-2 container with dimensions (num. points, num. basis functions)
00143     vals.resize(lineBasis.getCardinality(), lineNodes.dimension(0) );
00144     //INTREPID_TEST_COMMAND( lineBasis.getValues(vals, lineNodes, OPERATOR_DIV), throwCounter, nException );
00145 
00146     // Exceptions 1-5: all bf tags/bf Ids below are wrong and should cause getDofOrdinal() and
00147     // getDofTag() to access invalid array elements thereby causing bounds check exception
00148     // exception #1
00149     INTREPID_TEST_COMMAND( lineBasis.getDofOrdinal(2,0,0), throwCounter, nException );
00150     // exception #2
00151     INTREPID_TEST_COMMAND( lineBasis.getDofOrdinal(1,1,1), throwCounter, nException );
00152      // exception #3
00153     INTREPID_TEST_COMMAND( lineBasis.getDofOrdinal(1,0,7), throwCounter, nException );
00154     // exception #4
00155     INTREPID_TEST_COMMAND( lineBasis.getDofTag(6), throwCounter, nException );
00156     // exception #5
00157     INTREPID_TEST_COMMAND( lineBasis.getDofTag(-1), throwCounter, nException );
00158     // not an exception
00159     INTREPID_TEST_COMMAND( lineBasis.getDofTag(5), throwCounter, nException ); --nException;
00160 
00161 #ifdef HAVE_INTREPID_DEBUG
00162     // Exceptions 6-14 test exception handling with incorrectly dimensioned input/output arrays
00163     // exception #6: input points array must be of rank-2
00164     FieldContainer<double> badPoints1(4, 5, 3);
00165     INTREPID_TEST_COMMAND( lineBasis.getValues(vals, badPoints1, OPERATOR_VALUE), throwCounter, nException );
00166 
00167     // exception #7: dimension 1 in the input point array must equal space dimension of the cell
00168     FieldContainer<double> badPoints2(4, 3);
00169     INTREPID_TEST_COMMAND( lineBasis.getValues(vals, badPoints2, OPERATOR_VALUE), throwCounter, nException );
00170 
00171     // exception #8: output values must be of rank-2 for OPERATOR_VALUE
00172     FieldContainer<double> badVals1(4, 3, 1);
00173     INTREPID_TEST_COMMAND( lineBasis.getValues(badVals1, lineNodes, OPERATOR_VALUE), throwCounter, nException );
00174 
00175     // exception #9: output values must be of rank-3 for OPERATOR_GRAD
00176     FieldContainer<double> badVals2(4, 3);
00177     INTREPID_TEST_COMMAND( lineBasis.getValues(badVals2, lineNodes, OPERATOR_GRAD), throwCounter, nException );
00178 
00179     // exception #10: output values must be of rank-2 for OPERATOR_D1
00180     INTREPID_TEST_COMMAND( lineBasis.getValues(badVals2, lineNodes, OPERATOR_D1), throwCounter, nException );
00181 
00182     // exception #11: incorrect 0th dimension of output array (must equal number of basis functions)
00183     FieldContainer<double> badVals3(lineBasis.getCardinality() + 1, lineNodes.dimension(0));
00184     INTREPID_TEST_COMMAND( lineBasis.getValues(badVals3, lineNodes, OPERATOR_VALUE), throwCounter, nException );
00185     
00186     // exception #12: incorrect 1st dimension of output array (must equal number of points)
00187     FieldContainer<double> badVals4(lineBasis.getCardinality(), lineNodes.dimension(0) + 1);
00188     INTREPID_TEST_COMMAND( lineBasis.getValues(badVals4, lineNodes, OPERATOR_VALUE), throwCounter, nException );
00189 
00190     // exception #13: incorrect 2nd dimension of output array (must equal spatial dimension)
00191     FieldContainer<double> badVals5(lineBasis.getCardinality(), lineNodes.dimension(0), 2);
00192     INTREPID_TEST_COMMAND( lineBasis.getValues(badVals5, lineNodes, OPERATOR_GRAD), throwCounter, nException );
00193 
00194     
00195     // not an exception
00196     FieldContainer<double> goodVals2(lineBasis.getCardinality(), lineNodes.dimension(0));
00197     INTREPID_TEST_COMMAND( lineBasis.getValues(goodVals2, lineNodes, OPERATOR_VALUE), throwCounter, nException ); --nException;
00198 #endif
00199 
00200   }
00201   catch (std::logic_error err) {
00202     *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
00203     *outStream << err.what() << '\n';
00204     *outStream << "-------------------------------------------------------------------------------" << "\n\n";
00205     errorFlag = -1000;
00206   };
00207 
00208   // Check if number of thrown exceptions matches the one we expect
00209   if (throwCounter != nException) {
00210     errorFlag++;
00211     *outStream << std::setw(70) << "FAILURE! Incorrect number of exceptions." << "\n";
00212   }
00213 
00214    *outStream \
00215      << "\n"
00216      << "===============================================================================\n" \
00217      << "| TEST 2: correctness of basis function values                                |\n" \
00218      << "===============================================================================\n";
00219    outStream -> precision(20);
00220 
00221 
00222    try {
00223      // Check Kronecker property for Lagrange polynomials.
00224      int maxorder = 4;
00225 
00226      for (int ordi=1; ordi <= maxorder; ordi++) {
00227        FieldContainer<double> pts(PointTools::getLatticeSize(line,ordi),1);
00228        PointTools::getLattice<double,FieldContainer<double> >(pts,line,ordi);
00229        Basis_HGRAD_LINE_Cn_FEM<double, FieldContainer<double> > lineBasis(ordi, pts);
00230        FieldContainer<double> vals(lineBasis.getCardinality(),pts.dimension(0));
00231        lineBasis.getValues(vals,pts,OPERATOR_VALUE);
00232        for (int i=0;i<lineBasis.getCardinality();i++) {
00233         for (int j=0;j<pts.dimension(0);j++) {
00234           if ( i==j && std::abs( vals(i,j) - 1.0 ) > INTREPID_TOL ) {
00235             errorFlag++;
00236             *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00237             *outStream << " Basis function " << i << " does not have unit value at its node\n";
00238           }
00239           if ( i!=j && std::abs( vals(i,j) ) > INTREPID_TOL ) {
00240             errorFlag++;
00241             *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00242             *outStream << " Basis function " << i << " does not vanish at node " << j << "\n";
00243           }
00244         }
00245        }
00246      }
00247    }
00248   // Catch unexpected errors
00249   catch (std::logic_error err) {
00250     *outStream << err.what() << "\n\n";
00251     errorFlag = -1000;
00252   };
00253 
00254   if (errorFlag != 0)
00255     std::cout << "End Result: TEST FAILED\n";
00256   else
00257     std::cout << "End Result: TEST PASSED\n";
00258 
00259   // reset format state of std::cout
00260   std::cout.copyfmt(oldFormatState);
00261 
00262   return errorFlag;
00263 }