Intrepid
http://trilinos.sandia.gov/packages/docs/r10.12/packages/intrepid/test/Discretization/Basis/HGRAD_TET_COMP12_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_COMP12_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_COMP12_FEM)                          |\n" \
00093     << "|                                                                             |\n" \
00094     << "|     1) Evaluation of Basis Function Values                                  |\n" \
00095     << "|                                                                             |\n" \
00096     << "|  Questions? Contact  Pavel Bochev  (pbboche@sandia.gov),                    |\n" \
00097     << "|                      Denis Ridzal  (dridzal@sandia.gov),                    |\n" \
00098     << "|                      Kara Peterson (kjpeter@sandia.gov),                    |\n" \
00099     << "|                      Jake Ostien   (jtostie@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, values                                              |\n"\
00106     << "===============================================================================\n";
00107   
00108   // Define basis and error flag
00109   Basis_HGRAD_TET_COMP12_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 10 vertices of the reference TET
00117   FieldContainer<double> tetNodes(10, 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.5;  tetNodes(4,1) = 0.0;  tetNodes(4,2) = 0.0;
00123   tetNodes(5,0) = 0.5;  tetNodes(5,1) = 0.5;  tetNodes(5,2) = 0.0;
00124   tetNodes(6,0) = 0.0;  tetNodes(6,1) = 0.5;  tetNodes(6,2) = 0.0;
00125   tetNodes(7,0) = 0.0;  tetNodes(7,1) = 0.0;  tetNodes(7,2) = 0.5;  
00126   tetNodes(8,0) = 0.5;  tetNodes(8,1) = 0.0;  tetNodes(8,2) = 0.5;
00127   tetNodes(9,0) = 0.0;  tetNodes(9,1) = 0.5;  tetNodes(9,2) = 0.5;
00128   // Define array containing 5 integration points
00129   FieldContainer<double> tetPoints(5, 3);
00130   tetPoints(0,0) = 0.25;     tetPoints(0,1) = 0.25;     tetPoints(0,2) = 0.25;
00131   tetPoints(1,0) = 0.5;      tetPoints(1,1) = (1./6.);  tetPoints(1,2) = (1./6.);
00132   tetPoints(2,0) = (1./6.);  tetPoints(2,1) = 0.5;      tetPoints(2,2) = (1./6.);
00133   tetPoints(3,0) = (1./6.);  tetPoints(3,1) = (1./6.);  tetPoints(3,2) = 0.5;
00134   tetPoints(4,0) = (1./6.);  tetPoints(4,1) = (1./6.);  tetPoints(4,2) = (1./6.);
00135 
00136   // output precision
00137   outStream -> precision(20);
00138   
00139   // VALUE: Each row gives the 10 correct basis set values at an evaluation point
00140   double nodalBasisValues[] = {
00141     // first 4 vertices
00142     1.0, 0.0, 0.0, 0.0,  0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
00143     0.0, 1.0, 0.0, 0.0,  0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
00144     0.0, 0.0, 1.0, 0.0,  0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
00145     0.0, 0.0, 0.0, 1.0,  0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
00146     // second 6 vertices
00147     0.0, 0.0, 0.0, 0.0,  1.0, 0.0, 0.0, 0.0, 0.0, 0.0,
00148     0.0, 0.0, 0.0, 0.0,  0.0, 1.0, 0.0, 0.0, 0.0, 0.0,
00149     0.0, 0.0, 0.0, 0.0,  0.0, 0.0, 1.0, 0.0, 0.0, 0.0,
00150     0.0, 0.0, 0.0, 0.0,  0.0, 0.0, 0.0, 1.0, 0.0, 0.0,
00151     0.0, 0.0, 0.0, 0.0,  0.0, 0.0, 0.0, 0.0, 1.0, 0.0,
00152     0.0, 0.0, 0.0, 0.0,  0.0, 0.0, 0.0, 0.0, 0.0, 1.0
00153   };
00154   double pointBasisValues[] = {
00155     // pt 0 {0.25, 0.25, 0.25}
00156     0.0, 0.0, 0.0, 0.0,  1./6., 1./6., 1./6., 1./6., 1./6., 1./6.,
00157     // pt 1 {0.5, 1/6, 1/6}
00158     0.0, 0.0, 0.0, 0.0,  1./3., 1./3., 0.0, 0.0, 1./3., 0.0,
00159     // pt 2 {1/6, 0.5, 0.1/6}
00160     0.0, 0.0, 0.0, 0.0,  0.0, 1./3., 1./3., 0.0, 0.0, 1./3.,
00161     // pt 3 {1/6, 1/6, 0.5}
00162     0.0, 0.0, 0.0, 0.0,  0.0, 0.0, 0.0, 1./3., 1./3., 1./3.,
00163     // pt 4 {1/6, 1/6, 1/6}
00164     0.0, 0.0, 0.0, 0.0,  1./3., 0.0, 1./3., 1./3., 0.0, 0.0
00165   };
00166 
00167   // GRAD and D1: each row gives the 3x10 correct values of the gradients of the 10 basis functions
00168   double pointBasisGrads[] = {   
00169     // point 0
00170     -0.25, -0.25, -0.25,  0.25, 0.00,  0.00,   0.00,  0.25, 0.0,  0.00, 0.0, 0.25,  0.0, -0.75, -0.75,  \
00171      0.75,  0.75,  0.00, -0.75, 0.00, -0.75,  -0.75, -0.75, 0.0,  0.75, 0.0, 0.75,  0.0,  0.75,  0.75,
00172 
00173     // point 1
00174      0.15, 0.15, 0.15,     1.45, 0.00, 0.00,        0.00, -0.15, 0.0,    0.00, 0.0, -0.15,    -104./60., -2.15, -2.15,  \
00175      5./12., 2.15, 0.00,  -17./60., 0.00, -0.15,  -17./60., -0.15, 0.0,  5./12., 0.0,  2.15,  -2./15.,    0.15,  0.15,
00176 
00177     // point 2
00178     0.15, 0.15, 0.15,   -0.15, 0.0, 0.0,           0.0, 1.45, 0.0,       0.0, 0.0, -0.15,      0.0, -17./60., -0.15,  \
00179     2.15, 5./12., 0.0,  -2.15, -104./60., -2.15,  -0.15, -17./60., 0.0,  0.15, -2./15., 0.15,  0.0, 5./12., 2.15,
00180 
00181     // point 3
00182     0.15, 0.15, 0.15,     -0.15, 0.0, 0.0,        0.0, -0.15, 0.0,         0.0, 0.0, 1.45,     0.0, -0.15, -17./60.,  \
00183     0.15, 0.15, -2./15.,  -0.15, 0.0, -17./60.,  -2.15, -2.15, -104./60.,  2.15, 0.0, 5./12.,  0.0, 2.15, 5./12.,
00184 
00185     // point 4
00186     -1.45, -1.45, -1.45,       -0.15, 0.0, 0.0,              0.0, -0.15, 0.0,            0.0, 0.0, -0.15,           104./60., -5./12., -5./12.,  \
00187     17./60., 17./60., 2./15.,  -5./12., 104./60., -5./12.,  -5./12., -5./12., 104./60.,  17./60., 2./15., 17./60.,  2./15., 17./60., 17./60.
00188   };
00189 
00190   try{
00191         
00192     // Dimensions for the output arrays:
00193     int numFields = tetBasis.getCardinality();
00194     int numNodes  = tetNodes.dimension(0);
00195     int spaceDim  = tetBasis.getBaseCellTopology().getDimension();
00196     
00197     // Generic array for values, grads, curls, etc. that will be properly sized before each call
00198     FieldContainer<double> vals;
00199     
00200     // Check VALUE of basis functions at nodes: resize vals to rank-2 container:\n";
00201     vals.resize(numFields, numNodes);
00202     tetBasis.getValues(vals, tetNodes, OPERATOR_VALUE);
00203 
00204     for (int i = 0; i < numFields; i++) {
00205       for (int j = 0; j < numNodes; j++) {
00206           int l =  i + j * numFields;
00207           if (std::abs(vals(i,j) - nodalBasisValues[l]) > INTREPID_TOL) {
00208              errorFlag++;
00209              *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00210 
00211              // Output the multi-index of the value where the error is:
00212              *outStream << " At multi-index { ";
00213              *outStream << i << " ";*outStream << j << " ";
00214              *outStream << "}  computed value: " << vals(i,j)
00215                << " but reference value: " << nodalBasisValues[l] << "\n";
00216          }
00217       }
00218     }
00219 
00220     // Check VALUE of basis functions at points: resize vals to rank-2 container:\n";
00221     int numPoints = tetPoints.dimension(0);
00222     vals.resize(numFields, numPoints);
00223     vals.initialize(0.0);
00224     tetBasis.getValues(vals, tetPoints, OPERATOR_VALUE);
00225 
00226     for (int i = 0; i < numFields; i++) {
00227       for (int j = 0; j < numPoints; j++) {
00228           int l =  i + j * numFields;
00229           if (std::abs(vals(i,j) - pointBasisValues[l]) > INTREPID_TOL) {
00230              errorFlag++;
00231              *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00232 
00233              // Output the multi-index of the value where the error is:
00234              *outStream << " At multi-index { ";
00235              *outStream << i << " ";*outStream << j << " ";
00236              *outStream << "}  computed value: " << vals(i,j)
00237                << " but reference value: " << pointBasisValues[l] << "\n";
00238          }
00239       }
00240     }
00241 
00242     // Check GRAD of basis functions at points: resize vals to rank-3 container:\n";
00243     numPoints = tetPoints.dimension(0);
00244     vals.resize(numFields, numPoints, spaceDim);
00245     vals.initialize(0.0);
00246     tetBasis.getValues(vals, tetPoints, OPERATOR_GRAD);
00247 
00248      for (int i = 0; i < numFields; i++) {
00249        for (int j = 0; j < numPoints; j++) {
00250          for (int k = 0; k < spaceDim; k++) {
00251            int l = k + i * spaceDim + j * spaceDim * numFields;
00252            if (std::abs(vals(i,j,k) - pointBasisGrads[l]) > INTREPID_TOL) {
00253              errorFlag++;
00254              *outStream << std::setw(70) << "^^^^----FAILURE!" << "\n";
00255              
00256              // Output the multi-index of the value where the error is:
00257              *outStream << " At multi-index { ";
00258              *outStream << i << " ";*outStream << j << " ";*outStream << k << " ";
00259              *outStream << "}  computed grad component: " << vals(i,j,k)
00260                         << " but reference grad component: " << pointBasisGrads[l] << "\n";
00261            }
00262          }
00263        }
00264      }
00265   }
00266   // Catch unexpected errors
00267   catch (std::logic_error err) {
00268     *outStream << err.what() << "\n\n";
00269     errorFlag = -1000;
00270   };
00271 
00272   
00273   if (errorFlag != 0)
00274     std::cout << "End Result: TEST FAILED\n";
00275   else
00276     std::cout << "End Result: TEST PASSED\n";
00277   
00278   // reset format state of std::cout
00279   std::cout.copyfmt(oldFormatState);
00280   
00281   return errorFlag;
00282 }