|
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 00049 #include "Intrepid_FunctionSpaceTools.hpp" 00050 #include "Intrepid_FieldContainer.hpp" 00051 #include "Intrepid_CellTools.hpp" 00052 #include "Intrepid_HGRAD_TET_C1_FEM.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) basic operator transformations and integration in HGRAD |\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 #ifdef HAVE_INTREPID_DEBUG 00111 int beginThrowNumber = Teuchos::TestForException_getThrowNumber(); 00112 int endThrowNumber = beginThrowNumber + 28; 00113 #endif 00114 00115 typedef FunctionSpaceTools fst; 00116 00117 *outStream \ 00118 << "\n" 00119 << "===============================================================================\n"\ 00120 << "| TEST 1: exceptions |\n"\ 00121 << "===============================================================================\n"; 00122 00123 try{ 00124 #ifdef HAVE_INTREPID_DEBUG 00125 FieldContainer<double> a_2(2); 00126 FieldContainer<double> a_2_2(2, 2); 00127 FieldContainer<double> a_2_3(2, 3); 00128 FieldContainer<double> a_3_2(3, 2); 00129 FieldContainer<double> a_2_2_3(2, 2, 3); 00130 FieldContainer<double> a_2_2_3_3(2, 2, 3, 3); 00131 FieldContainer<double> a_2_2_2(2, 2, 2); 00132 FieldContainer<double> a_2_2_2_3_3(2, 2, 2, 3, 3); 00133 FieldContainer<double> a_2_2_2_2_2(2, 2, 2, 2, 2); 00134 FieldContainer<double> a_2_2_2_2(2, 2, 2, 2); 00135 FieldContainer<double> a_3_2_2_2(3, 2, 2, 2); 00136 FieldContainer<double> a_2_3_2_2(2, 3, 2, 2); 00137 FieldContainer<double> a_2_2_3_2(2, 2, 3, 2); 00138 FieldContainer<double> a_2_2_2_3(2, 2, 2, 3); 00139 00140 *outStream << "\n >>>>> TESTING computeCellMeasure:\n"; 00141 INTREPID_TEST_COMMAND( fst::computeCellMeasure<double>(a_2_2, a_2, a_2) ); 00142 INTREPID_TEST_COMMAND( fst::computeCellMeasure<double>(a_2_2, a_2_2, a_2) ); 00143 00144 *outStream << "\n >>>>> TESTING computeFaceMeasure:\n"; 00145 INTREPID_TEST_COMMAND( fst::computeFaceMeasure<double>(a_2_2, a_2, a_2, 0, shards::getCellTopologyData< shards::Tetrahedron<> >()) ); 00146 INTREPID_TEST_COMMAND( fst::computeFaceMeasure<double>(a_2_2, a_2_2_3_3, a_2, 0, shards::getCellTopologyData< shards::Tetrahedron<> >()) ); 00147 00148 *outStream << "\n >>>>> TESTING computeEdgeMeasure:\n"; 00149 INTREPID_TEST_COMMAND( fst::computeEdgeMeasure<double>(a_2_2, a_2, a_2, 0, shards::getCellTopologyData< shards::Triangle<> >()) ); 00150 INTREPID_TEST_COMMAND( fst::computeEdgeMeasure<double>(a_2_2, a_2_2_2_2, a_2, 0, shards::getCellTopologyData< shards::Triangle<> >()) ); 00151 00152 *outStream << "\n >>>>> TESTING integrate:\n"; 00153 INTREPID_TEST_COMMAND( fst::integrate<double>(a_2_2_2_2, a_2_2_2, a_2_2_2, COMP_CPP) ); 00154 INTREPID_TEST_COMMAND( fst::integrate<double>(a_2, a_2_2, a_2_2, COMP_CPP) ); 00155 INTREPID_TEST_COMMAND( fst::integrate<double>(a_2, a_2_2_3, a_2_2_3, COMP_CPP) ); 00156 INTREPID_TEST_COMMAND( fst::integrate<double>(a_2, a_2_2_3_3, a_2_2_3_3, COMP_CPP) ); 00157 INTREPID_TEST_COMMAND( fst::integrate<double>(a_2_2, a_2_2, a_2_2_2, COMP_CPP) ); 00158 INTREPID_TEST_COMMAND( fst::integrate<double>(a_2_2, a_2_2_3, a_2_2_2_3, COMP_CPP) ); 00159 INTREPID_TEST_COMMAND( fst::integrate<double>(a_2_2, a_2_2_3_3, a_2_2_2_3_3, COMP_CPP) ); 00160 INTREPID_TEST_COMMAND( fst::integrate<double>(a_2_2_2, a_2_2_2, a_2_2_2, COMP_CPP) ); 00161 INTREPID_TEST_COMMAND( fst::integrate<double>(a_2_2_2, a_2_2_2_3, a_2_2_2_3, COMP_CPP) ); 00162 INTREPID_TEST_COMMAND( fst::integrate<double>(a_2_2_2, a_2_2_2_3_3, a_2_2_2_3_3, COMP_CPP) ); 00163 00164 *outStream << "\n >>>>> TESTING operatorIntegral:\n"; 00165 INTREPID_TEST_COMMAND( fst::operatorIntegral<double>(a_2_2_2, a_2_2, a_2_2_2, COMP_CPP) ); 00166 INTREPID_TEST_COMMAND( fst::operatorIntegral<double>(a_2_2_2, a_2_2_2, a_2_2_2, COMP_CPP) ); 00167 INTREPID_TEST_COMMAND( fst::operatorIntegral<double>(a_2_2_2, a_2_2_2_3, a_2_2_2_3, COMP_CPP) ); 00168 INTREPID_TEST_COMMAND( fst::operatorIntegral<double>(a_2_2_2, a_2_2_2_3_3, a_2_2_2_3_3, COMP_CPP) ); 00169 00170 *outStream << "\n >>>>> TESTING functionalIntegral:\n"; 00171 INTREPID_TEST_COMMAND( fst::functionalIntegral<double>(a_2_2, a_2_2_2_3_3, a_2_2_2, COMP_CPP) ); 00172 INTREPID_TEST_COMMAND( fst::functionalIntegral<double>(a_2_2, a_2_2, a_2_2_2, COMP_CPP) ); 00173 INTREPID_TEST_COMMAND( fst::functionalIntegral<double>(a_2_2, a_2_2_3, a_2_2_2_3, COMP_CPP) ); 00174 INTREPID_TEST_COMMAND( fst::functionalIntegral<double>(a_2_2, a_2_2_3_3, a_2_2_2_3_3, COMP_CPP) ); 00175 00176 *outStream << "\n >>>>> TESTING dataIntegral:\n"; 00177 INTREPID_TEST_COMMAND( fst::dataIntegral<double>(a_2, a_2, a_2_2_2, COMP_CPP) ); 00178 INTREPID_TEST_COMMAND( fst::dataIntegral<double>(a_2, a_2_2, a_2_2, COMP_CPP) ); 00179 INTREPID_TEST_COMMAND( fst::dataIntegral<double>(a_2, a_2_2_3, a_2_2_3, COMP_CPP) ); 00180 INTREPID_TEST_COMMAND( fst::dataIntegral<double>(a_2, a_2_2_3_3, a_2_2_3_3, COMP_CPP) ); 00181 00182 *outStream << "\n >>>>> TESTING applyLeftFieldSigns:\n"; 00183 INTREPID_TEST_COMMAND( fst::applyLeftFieldSigns<double>(a_2, a_2) ); 00184 INTREPID_TEST_COMMAND( fst::applyLeftFieldSigns<double>(a_2_2_2, a_2) ); 00185 INTREPID_TEST_COMMAND( fst::applyLeftFieldSigns<double>(a_2_2_2, a_3_2) ); 00186 INTREPID_TEST_COMMAND( fst::applyLeftFieldSigns<double>(a_2_2_2, a_2_3) ); 00187 INTREPID_TEST_COMMAND( fst::applyLeftFieldSigns<double>(a_2_2_2, a_2_2) ); 00188 00189 *outStream << "\n >>>>> TESTING applyRightFieldSigns:\n"; 00190 INTREPID_TEST_COMMAND( fst::applyRightFieldSigns<double>(a_2, a_2) ); 00191 INTREPID_TEST_COMMAND( fst::applyRightFieldSigns<double>(a_2_2_2, a_2) ); 00192 INTREPID_TEST_COMMAND( fst::applyRightFieldSigns<double>(a_2_2_2, a_3_2) ); 00193 INTREPID_TEST_COMMAND( fst::applyRightFieldSigns<double>(a_2_2_2, a_2_3) ); 00194 INTREPID_TEST_COMMAND( fst::applyRightFieldSigns<double>(a_2_2_2, a_2_2) ); 00195 00196 *outStream << "\n >>>>> TESTING applyFieldSigns:\n"; 00197 INTREPID_TEST_COMMAND( fst::applyFieldSigns<double>(a_2, a_2) ); 00198 INTREPID_TEST_COMMAND( fst::applyFieldSigns<double>(a_2_2, a_2) ); 00199 INTREPID_TEST_COMMAND( fst::applyFieldSigns<double>(a_2_2, a_3_2) ); 00200 INTREPID_TEST_COMMAND( fst::applyFieldSigns<double>(a_2_2, a_2_3) ); 00201 INTREPID_TEST_COMMAND( fst::applyFieldSigns<double>(a_2_2_2_3_3, a_2_2) ); 00202 00203 *outStream << "\n >>>>> TESTING evaluate:\n"; 00204 INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2, a_2, a_2_2) ); 00205 INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2, a_2, a_2_2_2_3_3) ); 00206 INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2, a_2_2, a_2_2_2_3_3) ); 00207 INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2_2_3_3, a_3_2, a_2_2_2_3_3) ); 00208 INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2_2_3_3, a_2_3, a_2_2_2_3_3) ); 00209 INTREPID_TEST_COMMAND( fst::evaluate<double>(a_3_2_2_2, a_2_2, a_2_2_2_2_2) ); 00210 INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2_3_2_2, a_2_2, a_2_2_2_2_2) ); 00211 INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2_2_3_2, a_2_2, a_2_2_2_2_2) ); 00212 INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2_2_2_3, a_2_2, a_2_2_2_2_2) ); 00213 INTREPID_TEST_COMMAND( fst::evaluate<double>(a_2_2_2_2, a_2_2, a_2_2_2_2_2) ); 00214 #endif 00215 } 00216 catch (std::logic_error err) { 00217 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n"; 00218 *outStream << err.what() << '\n'; 00219 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00220 errorFlag = -1000; 00221 }; 00222 00223 #ifdef HAVE_INTREPID_DEBUG 00224 if (Teuchos::TestForException_getThrowNumber() != endThrowNumber) 00225 errorFlag++; 00226 #endif 00227 00228 *outStream \ 00229 << "\n" 00230 << "===============================================================================\n"\ 00231 << "| TEST 2: correctness of math operations |\n"\ 00232 << "===============================================================================\n"; 00233 00234 outStream->precision(20); 00235 00236 try { 00237 shards::CellTopology cellType = shards::getCellTopologyData< shards::Tetrahedron<> >(); // cell type: tet 00238 00239 /* Related to cubature. */ 00240 DefaultCubatureFactory<double> cubFactory; // create cubature factory 00241 int cubDegree = 2; // cubature degree 00242 Teuchos::RCP<Cubature<double> > myCub = cubFactory.create(cellType, cubDegree); // create default cubature 00243 int spaceDim = myCub->getDimension(); // get spatial dimension 00244 int numCubPoints = myCub->getNumPoints(); // get number of cubature points 00245 00246 /* Related to basis. */ 00247 Basis_HGRAD_TET_C1_FEM<double, FieldContainer<double> > tetBasis; // create tet basis 00248 int numFields = tetBasis.getCardinality(); // get basis cardinality 00249 00250 /* Cell geometries. */ 00251 int numCells = 4; 00252 int numNodes = 4; 00253 int numCellData = numCells*numNodes*spaceDim; 00254 double tetnodes[] = { 00255 // tet 0 00256 0.0, 0.0, 0.0, 00257 1.0, 0.0, 0.0, 00258 0.0, 1.0, 0.0, 00259 0.0, 0.0, 1.0, 00260 // tet 1 00261 4.0, 5.0, 1.0, 00262 -6.0, 2.0, 0.0, 00263 4.0, -3.0, -1.0, 00264 0.0, 2.0, 5.0, 00265 // tet 2 00266 -6.0, -3.0, 1.0, 00267 9.0, 2.0, 1.0, 00268 8.9, 2.1, 0.9, 00269 8.9, 2.1, 1.1, 00270 // tet 3 00271 -6.0, -3.0, 1.0, 00272 12.0, 3.0, 1.0, 00273 2.9, 0.1, 0.9, 00274 2.9, 0.1, 1.1 00275 }; 00276 00277 /* Computational arrays. */ 00278 FieldContainer<double> cub_points(numCubPoints, spaceDim); 00279 FieldContainer<double> cub_weights(numCubPoints); 00280 FieldContainer<double> cell_nodes(numCells, numNodes, spaceDim); 00281 FieldContainer<double> jacobian(numCells, numCubPoints, spaceDim, spaceDim); 00282 FieldContainer<double> jacobian_inv(numCells, numCubPoints, spaceDim, spaceDim); 00283 FieldContainer<double> jacobian_det(numCells, numCubPoints); 00284 FieldContainer<double> weighted_measure(numCells, numCubPoints); 00285 00286 FieldContainer<double> grad_of_basis_at_cub_points(numFields, numCubPoints, spaceDim); 00287 FieldContainer<double> transformed_grad_of_basis_at_cub_points(numCells, numFields, numCubPoints, spaceDim); 00288 FieldContainer<double> weighted_transformed_grad_of_basis_at_cub_points(numCells, numFields, numCubPoints, spaceDim); 00289 FieldContainer<double> stiffness_matrices(numCells, numFields, numFields); 00290 00291 FieldContainer<double> value_of_basis_at_cub_points(numFields, numCubPoints); 00292 FieldContainer<double> transformed_value_of_basis_at_cub_points(numCells, numFields, numCubPoints); 00293 FieldContainer<double> weighted_transformed_value_of_basis_at_cub_points(numCells, numFields, numCubPoints); 00294 FieldContainer<double> mass_matrices(numCells, numFields, numFields); 00295 00296 /******************* START COMPUTATION ***********************/ 00297 00298 // get cubature points and weights 00299 myCub->getCubature(cub_points, cub_weights); 00300 00301 // fill cell vertex array 00302 cell_nodes.setValues(tetnodes, numCellData); 00303 00304 // compute geometric cell information 00305 CellTools<double>::setJacobian(jacobian, cub_points, cell_nodes, cellType); 00306 CellTools<double>::setJacobianInv(jacobian_inv, jacobian); 00307 CellTools<double>::setJacobianDet(jacobian_det, jacobian); 00308 00309 // compute weighted measure 00310 fst::computeCellMeasure<double>(weighted_measure, jacobian_det, cub_weights); 00311 00312 // Computing stiffness matrices: 00313 // tabulate gradients of basis functions at (reference) cubature points 00314 tetBasis.getValues(grad_of_basis_at_cub_points, cub_points, OPERATOR_GRAD); 00315 00316 // transform gradients of basis functions 00317 fst::HGRADtransformGRAD<double>(transformed_grad_of_basis_at_cub_points, 00318 jacobian_inv, 00319 grad_of_basis_at_cub_points); 00320 00321 // multiply with weighted measure 00322 fst::multiplyMeasure<double>(weighted_transformed_grad_of_basis_at_cub_points, 00323 weighted_measure, 00324 transformed_grad_of_basis_at_cub_points); 00325 00326 // compute stiffness matrices 00327 fst::integrate<double>(stiffness_matrices, 00328 transformed_grad_of_basis_at_cub_points, 00329 weighted_transformed_grad_of_basis_at_cub_points, 00330 COMP_CPP); 00331 00332 00333 // Computing mass matrices: 00334 // tabulate values of basis functions at (reference) cubature points 00335 tetBasis.getValues(value_of_basis_at_cub_points, cub_points, OPERATOR_VALUE); 00336 00337 // transform values of basis functions 00338 fst::HGRADtransformVALUE<double>(transformed_value_of_basis_at_cub_points, 00339 value_of_basis_at_cub_points); 00340 00341 // multiply with weighted measure 00342 fst::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points, 00343 weighted_measure, 00344 transformed_value_of_basis_at_cub_points); 00345 00346 // compute mass matrices 00347 fst::integrate<double>(mass_matrices, 00348 transformed_value_of_basis_at_cub_points, 00349 weighted_transformed_value_of_basis_at_cub_points, 00350 COMP_CPP); 00351 00352 /******************* STOP COMPUTATION ***********************/ 00353 00354 00355 /******************* START COMPARISON ***********************/ 00356 string basedir = "./testdata"; 00357 for (int cell_id = 0; cell_id < numCells; cell_id++) { 00358 00359 stringstream namestream; 00360 string filename; 00361 namestream << basedir << "/mass_TET_FEM_P1" << "_" << "0" << cell_id+1 << ".dat"; 00362 namestream >> filename; 00363 00364 ifstream massfile(&filename[0]); 00365 if (massfile.is_open()) { 00366 if (compareToAnalytic<double>(&mass_matrices(cell_id, 0, 0), massfile, 1e-10, 0) > 0) 00367 errorFlag++; 00368 massfile.close(); 00369 } 00370 else { 00371 errorFlag = -1; 00372 std::cout << "End Result: TEST FAILED\n"; 00373 return errorFlag; 00374 } 00375 00376 namestream.clear(); 00377 namestream << basedir << "/stiff_TET_FEM_P1" << "_" << "0" << cell_id+1 << ".dat"; 00378 namestream >> filename; 00379 00380 ifstream stifffile(&filename[0]); 00381 if (stifffile.is_open()) 00382 { 00383 if (compareToAnalytic<double>(&stiffness_matrices(cell_id, 0, 0), stifffile, 1e-10, 0) > 0) 00384 errorFlag++; 00385 stifffile.close(); 00386 } 00387 else { 00388 errorFlag = -1; 00389 std::cout << "End Result: TEST FAILED\n"; 00390 return errorFlag; 00391 } 00392 00393 } 00394 00395 /******************* STOP COMPARISON ***********************/ 00396 00397 *outStream << "\n"; 00398 } 00399 catch (std::logic_error err) { 00400 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n"; 00401 *outStream << err.what() << '\n'; 00402 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00403 errorFlag = -1000; 00404 }; 00405 00406 00407 if (errorFlag != 0) 00408 std::cout << "End Result: TEST FAILED\n"; 00409 else 00410 std::cout << "End Result: TEST PASSED\n"; 00411 00412 // reset format state of std::cout 00413 std::cout.copyfmt(oldFormatState); 00414 00415 return errorFlag; 00416 }
1.7.6.1