|
Intrepid
|
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 }
1.7.6.1