|
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_HGRAD_TET_C1_FEM.hpp" 00040 #include "Intrepid_DefaultCubatureFactory.hpp" 00041 #include "Intrepid_Utils.hpp" 00042 #include "Teuchos_oblackholestream.hpp" 00043 #include "Teuchos_RCP.hpp" 00044 #include "Teuchos_GlobalMPISession.hpp" 00045 00046 using namespace std; 00047 using namespace Intrepid; 00048 00049 #define INTREPID_TEST_COMMAND( S ) \ 00050 { \ 00051 try { \ 00052 S ; \ 00053 } \ 00054 catch (std::logic_error err) { \ 00055 *outStream << "Expected Error ----------------------------------------------------------------\n"; \ 00056 *outStream << err.what() << '\n'; \ 00057 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; \ 00058 }; \ 00059 } 00060 00061 00062 int main(int argc, char *argv[]) { 00063 00064 Teuchos::GlobalMPISession mpiSession(&argc, &argv); 00065 00066 // This little trick lets us print to std::cout only if 00067 // a (dummy) command-line argument is provided. 00068 int iprint = argc - 1; 00069 Teuchos::RCP<std::ostream> outStream; 00070 Teuchos::oblackholestream bhs; // outputs nothing 00071 if (iprint > 0) 00072 outStream = Teuchos::rcp(&std::cout, false); 00073 else 00074 outStream = Teuchos::rcp(&bhs, false); 00075 00076 // Save the format state of the original std::cout. 00077 Teuchos::oblackholestream oldFormatState; 00078 oldFormatState.copyfmt(std::cout); 00079 00080 *outStream \ 00081 << "===============================================================================\n" \ 00082 << "| |\n" \ 00083 << "| Unit Test (FunctionSpaceTools) |\n" \ 00084 << "| |\n" \ 00085 << "| 1) basic operator transformations and integration in HGRAD |\n" \ 00086 << "| |\n" \ 00087 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \ 00088 << "| Denis Ridzal (dridzal@sandia.gov). |\n" \ 00089 << "| |\n" \ 00090 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \ 00091 << "| Trilinos website: http://trilinos.sandia.gov |\n" \ 00092 << "| |\n" \ 00093 << "===============================================================================\n"; 00094 00095 00096 int errorFlag = 0; 00097 #ifdef HAVE_INTREPID_DEBUG 00098 int beginThrowNumber = TestForException_getThrowNumber(); 00099 int endThrowNumber = beginThrowNumber + 28; 00100 #endif 00101 00102 typedef FunctionSpaceTools fst; 00103 00104 *outStream \ 00105 << "\n" 00106 << "===============================================================================\n"\ 00107 << "| TEST 1: exceptions |\n"\ 00108 << "===============================================================================\n"; 00109 00110 try{ 00111 #ifdef HAVE_INTREPID_DEBUG 00112 FieldContainer<double, 1> a_2(2); 00113 FieldContainer<double, 2> a_2_2(2, 2); 00114 FieldContainer<double, 3> a_2_3(2, 3); 00115 FieldContainer<double, 4> a_3_2(3, 2); 00116 FieldContainer<double, 5> a_2_2_3(2, 2, 3); 00117 FieldContainer<double, 6> a_2_2_3_3(2, 2, 3, 3); 00118 FieldContainer<double, 7> a_2_2_2(2, 2, 2); 00119 FieldContainer<double, 8> a_2_2_2_3_3(2, 2, 2, 3, 3); 00120 FieldContainer<double, 9> a_2_2_2_2_2(2, 2, 2, 2, 2); 00121 FieldContainer<double, 10> a_2_2_2_2(2, 2, 2, 2); 00122 FieldContainer<double, 11> a_3_2_2_2(3, 2, 2, 2); 00123 FieldContainer<double, 12> a_2_3_2_2(2, 3, 2, 2); 00124 FieldContainer<double, 13> a_2_2_3_2(2, 2, 3, 2); 00125 FieldContainer<double, 14> a_2_2_2_3(2, 2, 2, 3); 00126 00127 *outStream << "\n >>>>> TESTING computeCellMeasure:\n"; 00128 INTREPID_TEST_COMMAND( fst::computeCellMeasure<double>(a_2_2, a_2, a_2) ); 00129 INTREPID_TEST_COMMAND( fst::computeCellMeasure<double>(a_2_2, a_2_2, a_2) ); 00130 00131 *outStream << "\n >>>>> TESTING computeFaceMeasure:\n"; 00132 INTREPID_TEST_COMMAND( fst::computeFaceMeasure<double>(a_2_2, a_2, a_2, 0, shards::getCellTopologyData< shards::Tetrahedron<> >()) ); 00133 INTREPID_TEST_COMMAND( fst::computeFaceMeasure<double>(a_2_2, a_2_2_3_3, a_2, 0, shards::getCellTopologyData< shards::Tetrahedron<> >()) ); 00134 00135 *outStream << "\n >>>>> TESTING computeEdgeMeasure:\n"; 00136 INTREPID_TEST_COMMAND( fst::computeEdgeMeasure<double>(a_2_2, a_2, a_2, 0, shards::getCellTopologyData< shards::Triangle<> >()) ); 00137 INTREPID_TEST_COMMAND( fst::computeEdgeMeasure<double>(a_2_2, a_2_2_2_2, a_2, 0, shards::getCellTopologyData< shards::Triangle<> >()) ); 00138 00139 *outStream << "\n >>>>> TESTING integrate:\n"; 00140 INTREPID_TEST_COMMAND( fst::integrate<double>(a_2_2_2_2, a_2_2_2, a_2_2_2, COMP_CPP) ); 00141 INTREPID_TEST_COMMAND( fst::integrate<double>(a_2, a_2_2, a_2_2, COMP_CPP) ); 00142 INTREPID_TEST_COMMAND( fst::integrate<double>(a_2, a_2_2_3, a_2_2_3, COMP_CPP) ); 00143 INTREPID_TEST_COMMAND( fst::integrate<double>(a_2, a_2_2_3_3, a_2_2_3_3, COMP_CPP) ); 00144 INTREPID_TEST_COMMAND( fst::integrate<double>(a_2_2, a_2_2, a_2_2_2, COMP_CPP) ); 00145 INTREPID_TEST_COMMAND( fst::integrate<double>(a_2_2, a_2_2_3, a_2_2_2_3, COMP_CPP) ); 00146 INTREPID_TEST_COMMAND( fst::integrate<double>(a_2_2, a_2_2_3_3, a_2_2_2_3_3, COMP_CPP) ); 00147 INTREPID_TEST_COMMAND( fst::integrate<double>(a_2_2_2, a_2_2_2, a_2_2_2, COMP_CPP) ); 00148 INTREPID_TEST_COMMAND( fst::integrate<double>(a_2_2_2, a_2_2_2_3, a_2_2_2_3, COMP_CPP) ); 00149 INTREPID_TEST_COMMAND( fst::integrate<double>(a_2_2_2, a_2_2_2_3_3, a_2_2_2_3_3, COMP_CPP) ); 00150 00151 *outStream << "\n >>>>> TESTING operatorIntegral:\n"; 00152 INTREPID_TEST_COMMAND( fst::operatorIntegral<double>(a_2_2_2, a_2_2, a_2_2_2, COMP_CPP) ); 00153 INTREPID_TEST_COMMAND( fst::operatorIntegral<double>(a_2_2_2, a_2_2_2, a_2_2_2, COMP_CPP) ); 00154 INTREPID_TEST_COMMAND( fst::operatorIntegral<double>(a_2_2_2, a_2_2_2_3, a_2_2_2_3, COMP_CPP) ); 00155 INTREPID_TEST_COMMAND( fst::operatorIntegral<double>(a_2_2_2, a_2_2_2_3_3, a_2_2_2_3_3, COMP_CPP) ); 00156 00157 *outStream << "\n >>>>> TESTING functionalIntegral:\n"; 00158 INTREPID_TEST_COMMAND( fst::functionalIntegral<double>(a_2_2, a_2_2_2_3_3, a_2_2_2, COMP_CPP) ); 00159 INTREPID_TEST_COMMAND( fst::functionalIntegral<double>(a_2_2, a_2_2, a_2_2_2, COMP_CPP) ); 00160 INTREPID_TEST_COMMAND( fst::functionalIntegral<double>(a_2_2, a_2_2_3, a_2_2_2_3, COMP_CPP) ); 00161 INTREPID_TEST_COMMAND( fst::functionalIntegral<double>(a_2_2, a_2_2_3_3, a_2_2_2_3_3, COMP_CPP) ); 00162 00163 *outStream << "\n >>>>> TESTING dataIntegral:\n"; 00164 INTREPID_TEST_COMMAND( fst::dataIntegral<double>(a_2, a_2, a_2_2_2, COMP_CPP) ); 00165 INTREPID_TEST_COMMAND( fst::dataIntegral<double>(a_2, a_2_2, a_2_2, COMP_CPP) ); 00166 INTREPID_TEST_COMMAND( fst::dataIntegral<double>(a_2, a_2_2_3, a_2_2_3, COMP_CPP) ); 00167 INTREPID_TEST_COMMAND( fst::dataIntegral<double>(a_2, a_2_2_3_3, a_2_2_3_3, COMP_CPP) ); 00168 00169 *outStream << "\n >>>>> TESTING applyLeftFieldSigns:\n"; 00170 INTREPID_TEST_COMMAND( fst::applyLeftFieldSigns<double>(a_2, a_2) ); 00171 INTREPID_TEST_COMMAND( fst::applyLeftFieldSigns<double>(a_2_2_2, a_2) ); 00172 INTREPID_TEST_COMMAND( fst::applyLeftFieldSigns<double>(a_2_2_2, a_3_2) ); 00173 INTREPID_TEST_COMMAND( fst::applyLeftFieldSigns<double>(a_2_2_2, a_2_3) ); 00174 INTREPID_TEST_COMMAND( fst::applyLeftFieldSigns<double>(a_2_2_2, a_2_2) ); 00175 00176 *outStream << "\n >>>>> TESTING applyRightFieldSigns:\n"; 00177 INTREPID_TEST_COMMAND( fst::applyRightFieldSigns<double>(a_2, a_2) ); 00178 INTREPID_TEST_COMMAND( fst::applyRightFieldSigns<double>(a_2_2_2, a_2) ); 00179 INTREPID_TEST_COMMAND( fst::applyRightFieldSigns<double>(a_2_2_2, a_3_2) ); 00180 INTREPID_TEST_COMMAND( fst::applyRightFieldSigns<double>(a_2_2_2, a_2_3) ); 00181 INTREPID_TEST_COMMAND( fst::applyRightFieldSigns<double>(a_2_2_2, a_2_2) ); 00182 00183 *outStream << "\n >>>>> TESTING applyFieldSigns:\n"; 00184 INTREPID_TEST_COMMAND( fst::applyFieldSigns<double>(a_2, a_2) ); 00185 INTREPID_TEST_COMMAND( fst::applyFieldSigns<double>(a_2_2, a_2) ); 00186 INTREPID_TEST_COMMAND( fst::applyFieldSigns<double>(a_2_2, a_3_2) ); 00187 INTREPID_TEST_COMMAND( fst::applyFieldSigns<double>(a_2_2, a_2_3) ); 00188 INTREPID_TEST_COMMAND( fst::applyFieldSigns<double>(a_2_2_2_3_3, a_2_2) ); 00189 00190 *outStream << "\n >>>>> TESTING evaluate:\n"; 00191 INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2, a_2, a_2_2) ); 00192 INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2, a_2, a_2_2_2_3_3) ); 00193 INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2, a_2_2, a_2_2_2_3_3) ); 00194 INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2_2_3_3, a_3_2, a_2_2_2_3_3) ); 00195 INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2_2_3_3, a_2_3, a_2_2_2_3_3) ); 00196 INTREPID_TEST_COMMAND( fst::evaluate<double>(a_3_2_2_2, a_2_2, a_2_2_2_2_2) ); 00197 INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2_3_2_2, a_2_2, a_2_2_2_2_2) ); 00198 INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2_2_3_2, a_2_2, a_2_2_2_2_2) ); 00199 INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2_2_2_3, a_2_2, a_2_2_2_2_2) ); 00200 INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2_2_2_2, a_2_2, a_2_2_2_2_2) ); 00201 #endif 00202 } 00203 catch (std::logic_error err) { 00204 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n"; 00205 *outStream << err.what() << '\n'; 00206 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00207 errorFlag = -1000; 00208 }; 00209 00210 #ifdef HAVE_INTREPID_DEBUG 00211 if (TestForException_getThrowNumber() != endThrowNumber) 00212 errorFlag++; 00213 #endif 00214 00215 *outStream \ 00216 << "\n" 00217 << "===============================================================================\n"\ 00218 << "| TEST 2: correctness of math operations |\n"\ 00219 << "===============================================================================\n"; 00220 00221 outStream->precision(20); 00222 00223 try { 00224 // cell type: tet 00225 shards::CellTopology cellType = shards::getCellTopologyData< shards::Tetrahedron<> >(); 00226 00227 /* Related to cubature. */ 00228 // create cubature factory 00229 DefaultCubatureFactory<double, FieldContainer<double, 1>, FieldContainer<double, 2> > cubFactory; 00230 // cubature degree 00231 int cubDegree = 2; 00232 // create default cubature 00233 Teuchos::RCP<Cubature<double, FieldContainer<double, 1>, FieldContainer<double, 2> > > myCub = cubFactory.create(cellType, cubDegree); 00234 // get spatial dimension 00235 int spaceDim = myCub->getDimension(); 00236 // get number of cubature points 00237 int numCubPoints = myCub->getNumPoints(); 00238 00239 /* Related to basis. */ 00240 // create tet basis 00241 Basis_HGRAD_TET_C1_FEM<double, FieldContainer<double, 1> > tetBasis; 00242 // get basis cardinality 00243 int numFields = tetBasis.getCardinality(); 00244 00245 /* Cell geometries. */ 00246 int numCells = 4; 00247 int numNodes = 4; 00248 int numCellData = numCells*numNodes*spaceDim; 00249 double tetnodes[] = { 00250 // tet 0 00251 0.0, 0.0, 0.0, 00252 1.0, 0.0, 0.0, 00253 0.0, 1.0, 0.0, 00254 0.0, 0.0, 1.0, 00255 // tet 1 00256 4.0, 5.0, 1.0, 00257 -6.0, 2.0, 0.0, 00258 4.0, -3.0, -1.0, 00259 0.0, 2.0, 5.0, 00260 // tet 2 00261 -6.0, -3.0, 1.0, 00262 9.0, 2.0, 1.0, 00263 8.9, 2.1, 0.9, 00264 8.9, 2.1, 1.1, 00265 // tet 3 00266 -6.0, -3.0, 1.0, 00267 12.0, 3.0, 1.0, 00268 2.9, 0.1, 0.9, 00269 2.9, 0.1, 1.1 00270 }; 00271 00272 /* Computational arrays. */ 00273 FieldContainer<double, 1> cub_points(numCubPoints, spaceDim); 00274 FieldContainer<double, 2> cub_weights(numCubPoints); 00275 FieldContainer<double, 3> cell_nodes(numCells, numNodes, spaceDim); 00276 FieldContainer<double, 4> jacobian(numCells, numCubPoints, spaceDim, spaceDim); 00277 FieldContainer<double, 5> jacobian_inv(numCells, numCubPoints, spaceDim, spaceDim); 00278 FieldContainer<double, 6> jacobian_det(numCells, numCubPoints); 00279 FieldContainer<double, 7> weighted_measure(numCells, numCubPoints); 00280 00281 FieldContainer<double, 1> grad_of_basis_at_cub_points(numFields, numCubPoints, spaceDim); 00282 FieldContainer<double, 9> transformed_grad_of_basis_at_cub_points(numCells, numFields, numCubPoints, spaceDim); 00283 FieldContainer<double, 10> weighted_transformed_grad_of_basis_at_cub_points(numCells, numFields, numCubPoints, spaceDim); 00284 FieldContainer<double, 11> stiffness_matrices(numCells, numFields, numFields); 00285 00286 FieldContainer<double, 1> value_of_basis_at_cub_points(numFields, numCubPoints); 00287 FieldContainer<double, 13> transformed_value_of_basis_at_cub_points(numCells, numFields, numCubPoints); 00288 FieldContainer<double, 14> weighted_transformed_value_of_basis_at_cub_points(numCells, numFields, numCubPoints); 00289 FieldContainer<double, 15> mass_matrices(numCells, numFields, numFields); 00290 00291 /******************* START COMPUTATION ***********************/ 00292 00293 // get cubature points and weights 00294 myCub->getCubature(cub_points, cub_weights); 00295 00296 // fill cell vertex array 00297 cell_nodes.setValues(tetnodes, numCellData); 00298 00299 // compute geometric cell information 00300 CellTools<double>::setJacobian(jacobian, cub_points, cell_nodes, cellType); 00301 CellTools<double>::setJacobianInv(jacobian_inv, jacobian); 00302 CellTools<double>::setJacobianDet(jacobian_det, jacobian); 00303 00304 // compute weighted measure 00305 fst::computeCellMeasure<double>(weighted_measure, jacobian_det, cub_weights); 00306 00307 // Computing stiffness matrices: 00308 // tabulate gradients of basis functions at (reference) cubature points 00309 tetBasis.getValues(grad_of_basis_at_cub_points, cub_points, OPERATOR_GRAD); 00310 00311 // transform gradients of basis functions 00312 fst::HGRADtransformGRAD<double>(transformed_grad_of_basis_at_cub_points, 00313 jacobian_inv, 00314 grad_of_basis_at_cub_points); 00315 00316 // multiply with weighted measure 00317 fst::multiplyMeasure<double>(weighted_transformed_grad_of_basis_at_cub_points, 00318 weighted_measure, 00319 transformed_grad_of_basis_at_cub_points); 00320 00321 // compute stiffness matrices 00322 fst::integrate<double>(stiffness_matrices, 00323 transformed_grad_of_basis_at_cub_points, 00324 weighted_transformed_grad_of_basis_at_cub_points, 00325 COMP_CPP); 00326 00327 00328 // Computing mass matrices: 00329 // tabulate values of basis functions at (reference) cubature points 00330 tetBasis.getValues(value_of_basis_at_cub_points, cub_points, OPERATOR_VALUE); 00331 00332 // transform values of basis functions 00333 fst::HGRADtransformVALUE<double>(transformed_value_of_basis_at_cub_points, 00334 value_of_basis_at_cub_points); 00335 00336 // multiply with weighted measure 00337 fst::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points, 00338 weighted_measure, 00339 transformed_value_of_basis_at_cub_points); 00340 00341 // compute mass matrices 00342 fst::integrate<double>(mass_matrices, 00343 transformed_value_of_basis_at_cub_points, 00344 weighted_transformed_value_of_basis_at_cub_points, 00345 COMP_CPP); 00346 00347 /******************* STOP COMPUTATION ***********************/ 00348 00349 00350 /******************* START COMPARISON ***********************/ 00351 string basedir = "./testdata"; 00352 for (int cell_id = 0; cell_id < numCells; cell_id++) { 00353 00354 stringstream namestream; 00355 string filename; 00356 namestream << basedir << "/mass_TET_FEM_P1" << "_" << "0" << cell_id+1 << ".dat"; 00357 namestream >> filename; 00358 00359 ifstream massfile(&filename[0]); 00360 if (massfile.is_open()) { 00361 if (compareToAnalytic<double>(&mass_matrices(cell_id, 0, 0), massfile, 1e-10, 0) > 0) 00362 errorFlag++; 00363 massfile.close(); 00364 } 00365 else { 00366 errorFlag = -1; 00367 std::cout << "End Result: TEST FAILED\n"; 00368 return errorFlag; 00369 } 00370 00371 namestream.clear(); 00372 namestream << basedir << "/stiff_TET_FEM_P1" << "_" << "0" << cell_id+1 << ".dat"; 00373 namestream >> filename; 00374 00375 ifstream stifffile(&filename[0]); 00376 if (stifffile.is_open()) 00377 { 00378 if (compareToAnalytic<double>(&stiffness_matrices(cell_id, 0, 0), stifffile, 1e-10, 0) > 0) 00379 errorFlag++; 00380 stifffile.close(); 00381 } 00382 else { 00383 errorFlag = -1; 00384 std::cout << "End Result: TEST FAILED\n"; 00385 return errorFlag; 00386 } 00387 00388 } 00389 00390 /******************* STOP COMPARISON ***********************/ 00391 00392 *outStream << "\n"; 00393 } 00394 catch (std::logic_error err) { 00395 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n"; 00396 *outStream << err.what() << '\n'; 00397 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00398 errorFlag = -1000; 00399 }; 00400 00401 00402 if (errorFlag != 0) 00403 std::cout << "End Result: TEST FAILED\n"; 00404 else 00405 std::cout << "End Result: TEST PASSED\n"; 00406 00407 // reset format state of std::cout 00408 std::cout.copyfmt(oldFormatState); 00409 00410 return errorFlag; 00411 }
1.7.4