|
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 // This library is free software; you can redistribute it and/or modify 00011 // it under the terms of the GNU Lesser General Public License as 00012 // published by the Free Software Foundation; either version 2.1 of the 00013 // License, or (at your option) any later version. 00014 // 00015 // This library is distributed in the hope that it will be useful, but 00016 // WITHOUT ANY WARRANTY; without even the implied warranty of 00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 00018 // Lesser General Public License for more details. 00019 // 00020 // You should have received a copy of the GNU Lesser General Public 00021 // License along with this library; if not, write to the Free Software 00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 00023 // USA 00024 // Questions? Contact Pavel Bochev (pbboche@sandia.gov) or 00025 // Denis Ridzal (dridzal@sandia.gov). 00026 // 00027 // ************************************************************************ 00028 // @HEADER 00029 00036 #include "Intrepid_FunctionSpaceTools.hpp" 00037 #include "Intrepid_FieldContainer.hpp" 00038 #include "Intrepid_CellTools.hpp" 00039 #include "Intrepid_DefaultCubatureFactory.hpp" 00040 #include "Intrepid_Utils.hpp" 00041 #include "Teuchos_oblackholestream.hpp" 00042 #include "Teuchos_RCP.hpp" 00043 #include "Teuchos_GlobalMPISession.hpp" 00044 00045 using namespace std; 00046 using namespace Intrepid; 00047 00048 #define INTREPID_TEST_COMMAND( S ) \ 00049 { \ 00050 try { \ 00051 S ; \ 00052 } \ 00053 catch (std::logic_error err) { \ 00054 *outStream << "Expected Error ----------------------------------------------------------------\n"; \ 00055 *outStream << err.what() << '\n'; \ 00056 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \ 00057 }; \ 00058 } 00059 00060 00061 int main(int argc, char *argv[]) { 00062 00063 Teuchos::GlobalMPISession mpiSession(&argc, &argv); 00064 00065 // This little trick lets us print to std::cout only if 00066 // a (dummy) command-line argument is provided. 00067 int iprint = argc - 1; 00068 Teuchos::RCP<std::ostream> outStream; 00069 Teuchos::oblackholestream bhs; // outputs nothing 00070 if (iprint > 0) 00071 outStream = Teuchos::rcp(&std::cout, false); 00072 else 00073 outStream = Teuchos::rcp(&bhs, false); 00074 00075 // Save the format state of the original std::cout. 00076 Teuchos::oblackholestream oldFormatState; 00077 oldFormatState.copyfmt(std::cout); 00078 00079 *outStream \ 00080 << "===============================================================================\n" \ 00081 << "| |\n" \ 00082 << "| Unit Test (FunctionSpaceTools) |\n" \ 00083 << "| |\n" \ 00084 << "| 1) volume integration on tetrahedra, testing dataIntegral |\n" \ 00085 << "| |\n" \ 00086 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \ 00087 << "| Denis Ridzal (dridzal@sandia.gov). |\n" \ 00088 << "| |\n" \ 00089 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \ 00090 << "| Trilinos website: http://trilinos.sandia.gov |\n" \ 00091 << "| |\n" \ 00092 << "===============================================================================\n"; 00093 00094 00095 int errorFlag = 0; 00096 00097 typedef FunctionSpaceTools fst; 00098 00099 *outStream \ 00100 << "\n" 00101 << "===============================================================================\n"\ 00102 << "| TEST 1: correctness of cell volumes |\n"\ 00103 << "===============================================================================\n"; 00104 00105 outStream->precision(20); 00106 00107 try { 00108 shards::CellTopology cellType = shards::getCellTopologyData< shards::Tetrahedron<> >(); // cell type: tet 00109 00110 /* Related to cubature. */ 00111 DefaultCubatureFactory<double> cubFactory; // create cubature factory 00112 int cubDegree = 0; // cubature degree 00113 Teuchos::RCP<Cubature<double> > myCub = cubFactory.create(cellType, cubDegree); // create default cubature 00114 int spaceDim = myCub->getDimension(); // get spatial dimension 00115 int numCubPoints = myCub->getNumPoints(); // get number of cubature points 00116 00117 /* Cell geometries. */ 00118 int numCells = 4; 00119 int numNodes = 4; 00120 int numCellData = numCells*numNodes*spaceDim; 00121 double tetnodes[] = { 00122 // tet 0 00123 0.0, 0.0, 0.0, 00124 1.0, 0.0, 0.0, 00125 0.0, 1.0, 0.0, 00126 0.0, 0.0, 1.0, 00127 // tet 1 00128 4.0, 5.0, 1.0, 00129 -6.0, 2.0, 0.0, 00130 4.0, -3.0, -1.0, 00131 0.0, 2.0, 5.0, 00132 // tet 2 00133 -6.0, -3.0, 1.0, 00134 9.0, 2.0, 1.0, 00135 8.9, 2.1, 0.9, 00136 8.9, 2.1, 1.1, 00137 // tet 3 00138 -6.0, -3.0, 1.0, 00139 12.0, 3.0, 1.0, 00140 2.9, 0.1, 0.9, 00141 2.9, 0.1, 1.1 00142 }; 00143 00144 /* Analytic volumes. */ 00145 double tetvols[] = {1.0/6.0, 194.0/3.0, 1.0/15.0, 2.0/25.0}; 00146 00147 /* Computational arrays. */ 00148 FieldContainer<double> cub_points(numCubPoints, spaceDim); 00149 FieldContainer<double> cub_weights(numCubPoints); 00150 FieldContainer<double> cell_nodes(numCells, numNodes, spaceDim); 00151 FieldContainer<double> jacobian(numCells, numCubPoints, spaceDim, spaceDim); 00152 FieldContainer<double> jacobian_det(numCells, numCubPoints); 00153 FieldContainer<double> weighted_measure(numCells, numCubPoints); 00154 FieldContainer<double> data_one(numCells, numCubPoints); 00155 FieldContainer<double> volumes(numCells); 00156 00157 /******************* START COMPUTATION ***********************/ 00158 00159 // get cubature points and weights 00160 myCub->getCubature(cub_points, cub_weights); 00161 00162 // fill cell vertex array 00163 cell_nodes.setValues(tetnodes, numCellData); 00164 00165 // compute geometric cell information 00166 CellTools<double>::setJacobian(jacobian, cub_points, cell_nodes, cellType); 00167 CellTools<double>::setJacobianDet(jacobian_det, jacobian); 00168 00169 // compute weighted measure 00170 fst::computeCellMeasure<double>(weighted_measure, jacobian_det, cub_weights); 00171 00172 // set data to 1.0 00173 for (int cell=0; cell<data_one.dimension(0); cell++) { 00174 for (int qp=0; qp<data_one.dimension(1); qp++) { 00175 data_one(cell,qp) = 1.0; 00176 } 00177 } 00178 00179 // compute volumes 00180 fst::integrate<double>(volumes, data_one, weighted_measure, COMP_CPP); 00181 00182 /******************* STOP COMPUTATION ***********************/ 00183 00184 00185 /******************* START COMPARISON ***********************/ 00186 for (int cell_id = 0; cell_id < numCells; cell_id++) { 00187 *outStream << "Volume of cell " << cell_id << " = " << volumes(cell_id) << " vs. Analytic value = " << tetvols[cell_id] << "\n"; 00188 if (std::fabs(volumes(cell_id)-tetvols[cell_id]) > INTREPID_TOL) { 00189 errorFlag++; 00190 } 00191 } 00192 /******************* STOP COMPARISON ***********************/ 00193 00194 *outStream << "\n"; 00195 } 00196 catch (std::logic_error err) { 00197 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n"; 00198 *outStream << err.what() << '\n'; 00199 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00200 errorFlag = -1000; 00201 }; 00202 00203 00204 if (errorFlag != 0) 00205 std::cout << "End Result: TEST FAILED\n"; 00206 else 00207 std::cout << "End Result: TEST PASSED\n"; 00208 00209 // reset format state of std::cout 00210 std::cout.copyfmt(oldFormatState); 00211 00212 return errorFlag; 00213 }
1.7.4