Intrepid
http://trilinos.sandia.gov/packages/docs/r10.10/packages/intrepid/test/Discretization/FunctionSpaceTools/test_04.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 
00050 #include "Intrepid_FunctionSpaceTools.hpp"
00051 #include "Intrepid_FieldContainer.hpp"
00052 #include "Intrepid_CellTools.hpp"
00053 #include "Intrepid_DefaultCubatureFactory.hpp"
00054 #include "Intrepid_Utils.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 )                                                                                  \
00063 {                                                                                                                   \
00064   try {                                                                                                             \
00065     S ;                                                                                                             \
00066   }                                                                                                                 \
00067   catch (std::logic_error err) {                                                                                    \
00068       *outStream << "Expected Error ----------------------------------------------------------------\n";            \
00069       *outStream << err.what() << '\n';                                                                             \
00070       *outStream << "-------------------------------------------------------------------------------" << "\n\n";    \
00071   };                                                                                                                \
00072 }
00073 
00074 
00075 int main(int argc, char *argv[]) {
00076 
00077   Teuchos::GlobalMPISession mpiSession(&argc, &argv);
00078 
00079   // This little trick lets us print to std::cout only if
00080   // a (dummy) command-line argument is provided.
00081   int iprint     = argc - 1;
00082   Teuchos::RCP<std::ostream> outStream;
00083   Teuchos::oblackholestream bhs; // outputs nothing
00084   if (iprint > 0)
00085     outStream = Teuchos::rcp(&std::cout, false);
00086   else
00087     outStream = Teuchos::rcp(&bhs, false);
00088 
00089   // Save the format state of the original std::cout.
00090   Teuchos::oblackholestream oldFormatState;
00091   oldFormatState.copyfmt(std::cout);
00092 
00093   *outStream \
00094   << "===============================================================================\n" \
00095   << "|                                                                             |\n" \
00096   << "|                      Unit Test (FunctionSpaceTools)                         |\n" \
00097   << "|                                                                             |\n" \
00098   << "|     1) volume integration on tetrahedra, testing dataIntegral               |\n" \
00099   << "|                                                                             |\n" \
00100   << "|  Questions? Contact  Pavel Bochev (pbboche@sandia.gov) or                   |\n" \
00101   << "|                      Denis Ridzal (dridzal@sandia.gov).                     |\n" \
00102   << "|                                                                             |\n" \
00103   << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00104   << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00105   << "|                                                                             |\n" \
00106   << "===============================================================================\n";
00107 
00108 
00109   int errorFlag = 0;
00110 
00111   typedef FunctionSpaceTools fst; 
00112 
00113   *outStream \
00114   << "\n"
00115   << "===============================================================================\n"\
00116   << "| TEST 1: correctness of cell volumes                                         |\n"\
00117   << "===============================================================================\n";
00118 
00119   outStream->precision(20);
00120 
00121   try {
00122       shards::CellTopology cellType = shards::getCellTopologyData< shards::Tetrahedron<> >();   // cell type: tet
00123 
00124       /* Related to cubature. */
00125       DefaultCubatureFactory<double> cubFactory;                                                // create cubature factory
00126       int cubDegree = 0;                                                                        // cubature degree
00127       Teuchos::RCP<Cubature<double> > myCub = cubFactory.create(cellType, cubDegree);           // create default cubature
00128       int spaceDim = myCub->getDimension();                                                     // get spatial dimension 
00129       int numCubPoints = myCub->getNumPoints();                                                 // get number of cubature points
00130 
00131       /* Cell geometries. */
00132       int numCells    = 4;
00133       int numNodes    = 4;
00134       int numCellData = numCells*numNodes*spaceDim;
00135       double tetnodes[] = {
00136         // tet 0
00137         0.0, 0.0, 0.0,
00138         1.0, 0.0, 0.0,
00139         0.0, 1.0, 0.0,
00140         0.0, 0.0, 1.0,
00141         // tet 1
00142         4.0, 5.0, 1.0,
00143        -6.0, 2.0, 0.0,
00144         4.0, -3.0, -1.0,
00145         0.0, 2.0, 5.0,
00146         // tet 2
00147         -6.0, -3.0, 1.0,
00148         9.0, 2.0, 1.0,
00149         8.9, 2.1, 0.9,
00150         8.9, 2.1, 1.1,
00151         // tet 3
00152         -6.0, -3.0, 1.0,
00153         12.0, 3.0, 1.0,
00154         2.9, 0.1, 0.9,
00155         2.9, 0.1, 1.1
00156       };
00157 
00158       /* Analytic volumes. */
00159       double tetvols[] = {1.0/6.0, 194.0/3.0, 1.0/15.0, 2.0/25.0};
00160 
00161       /* Computational arrays. */
00162       FieldContainer<double> cub_points(numCubPoints, spaceDim);
00163       FieldContainer<double> cub_weights(numCubPoints);
00164       FieldContainer<double> cell_nodes(numCells, numNodes, spaceDim);
00165       FieldContainer<double> jacobian(numCells, numCubPoints, spaceDim, spaceDim);
00166       FieldContainer<double> jacobian_det(numCells, numCubPoints);
00167       FieldContainer<double> weighted_measure(numCells, numCubPoints);
00168       FieldContainer<double> data_one(numCells, numCubPoints);
00169       FieldContainer<double> volumes(numCells);
00170 
00171       /******************* START COMPUTATION ***********************/
00172 
00173       // get cubature points and weights
00174       myCub->getCubature(cub_points, cub_weights);
00175 
00176       // fill cell vertex array
00177       cell_nodes.setValues(tetnodes, numCellData);
00178 
00179       // compute geometric cell information
00180       CellTools<double>::setJacobian(jacobian, cub_points, cell_nodes, cellType);
00181       CellTools<double>::setJacobianDet(jacobian_det, jacobian);
00182 
00183       // compute weighted measure
00184       fst::computeCellMeasure<double>(weighted_measure, jacobian_det, cub_weights);
00185 
00186       // set data to 1.0
00187       for (int cell=0; cell<data_one.dimension(0); cell++) { 
00188         for (int qp=0; qp<data_one.dimension(1); qp++) {
00189           data_one(cell,qp) = 1.0;
00190         }
00191       } 
00192 
00193       // compute volumes
00194       fst::integrate<double>(volumes, data_one, weighted_measure, COMP_CPP);
00195 
00196       /******************* STOP COMPUTATION ***********************/
00197 
00198 
00199       /******************* START COMPARISON ***********************/
00200       for (int cell_id = 0; cell_id < numCells; cell_id++) {
00201         *outStream << "Volume of cell " << cell_id << " = " << volumes(cell_id) << "    vs.    Analytic value =  " << tetvols[cell_id] << "\n";
00202         if (std::fabs(volumes(cell_id)-tetvols[cell_id]) > INTREPID_TOL) {
00203           errorFlag++;
00204         }
00205       }
00206       /******************* STOP COMPARISON ***********************/
00207 
00208       *outStream << "\n";
00209   }
00210   catch (std::logic_error err) {
00211     *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n";
00212     *outStream << err.what() << '\n';
00213     *outStream << "-------------------------------------------------------------------------------" << "\n\n";
00214     errorFlag = -1000;
00215   };
00216 
00217 
00218   if (errorFlag != 0)
00219     std::cout << "End Result: TEST FAILED\n";
00220   else
00221     std::cout << "End Result: TEST PASSED\n";
00222 
00223   // reset format state of std::cout
00224   std::cout.copyfmt(oldFormatState);
00225 
00226   return errorFlag;
00227 }