|
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 00038 #include "Intrepid_FunctionSpaceTools.hpp" 00039 #include "Intrepid_FieldContainer.hpp" 00040 #include "Shards_Array.hpp" 00041 #include "Intrepid_CellTools.hpp" 00042 #include "Intrepid_HDIV_HEX_I1_FEM.hpp" 00043 #include "Intrepid_DefaultCubatureFactory.hpp" 00044 #include "Intrepid_Utils.hpp" 00045 #include "Teuchos_oblackholestream.hpp" 00046 #include "Teuchos_RCP.hpp" 00047 #include "Teuchos_GlobalMPISession.hpp" 00048 00049 using namespace std; 00050 using namespace Intrepid; 00051 00052 SHARDS_ARRAY_DIM_TAG_SIMPLE_DECLARATION( Cell ) 00053 SHARDS_ARRAY_DIM_TAG_SIMPLE_IMPLEMENTATION( Cell ) 00054 00055 SHARDS_ARRAY_DIM_TAG_SIMPLE_DECLARATION( Field ) 00056 SHARDS_ARRAY_DIM_TAG_SIMPLE_IMPLEMENTATION( Field ) 00057 00058 SHARDS_ARRAY_DIM_TAG_SIMPLE_DECLARATION( Point ) 00059 SHARDS_ARRAY_DIM_TAG_SIMPLE_IMPLEMENTATION( Point ) 00060 00061 SHARDS_ARRAY_DIM_TAG_SIMPLE_DECLARATION( Dim ) 00062 SHARDS_ARRAY_DIM_TAG_SIMPLE_IMPLEMENTATION( Dim ) 00063 00064 SHARDS_ARRAY_DIM_TAG_SIMPLE_DECLARATION( Node ) 00065 SHARDS_ARRAY_DIM_TAG_SIMPLE_IMPLEMENTATION( Node ) 00066 00067 00068 typedef FieldContainer<double> FC; 00069 00070 00071 int main(int argc, char *argv[]) { 00072 00073 Teuchos::GlobalMPISession mpiSession(&argc, &argv); 00074 00075 // This little trick lets us print to std::cout only if 00076 // a (dummy) command-line argument is provided. 00077 int iprint = argc - 1; 00078 Teuchos::RCP<std::ostream> outStream; 00079 Teuchos::oblackholestream bhs; // outputs nothing 00080 if (iprint > 0) 00081 outStream = Teuchos::rcp(&std::cout, false); 00082 else 00083 outStream = Teuchos::rcp(&bhs, false); 00084 00085 // Save the format state of the original std::cout. 00086 Teuchos::oblackholestream oldFormatState; 00087 oldFormatState.copyfmt(std::cout); 00088 00089 *outStream \ 00090 << "===============================================================================\n" \ 00091 << "| |\n" \ 00092 << "| Unit Test (FunctionSpaceTools) |\n" \ 00093 << "| |\n" \ 00094 << "| 1) basic operator transformations and integration in HDIV |\n" \ 00095 << "| via shards::Array wrappers |\n" \ 00096 << "| |\n" \ 00097 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov) or |\n" \ 00098 << "| Denis Ridzal (dridzal@sandia.gov). |\n" \ 00099 << "| |\n" \ 00100 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \ 00101 << "| Trilinos website: http://trilinos.sandia.gov |\n" \ 00102 << "| |\n" \ 00103 << "===============================================================================\n"; 00104 00105 00106 int errorFlag = 0; 00107 00108 typedef FunctionSpaceTools fst; 00109 00110 *outStream \ 00111 << "\n" 00112 << "===============================================================================\n"\ 00113 << "| TEST 1: correctness of math operations |\n"\ 00114 << "===============================================================================\n"; 00115 00116 outStream->precision(20); 00117 00118 try { 00119 shards::CellTopology cellType = shards::getCellTopologyData< shards::Hexahedron<> >(); // cell type: hex 00120 00121 /* Related to cubature. */ 00122 DefaultCubatureFactory<double> cubFactory; // create cubature factory 00123 int cubDegree = 20; // cubature degree 00124 Teuchos::RCP<Cubature<double> > myCub = cubFactory.create(cellType, cubDegree); // create default cubature 00125 int spaceDim = myCub->getDimension(); // get spatial dimension 00126 int numCubPoints = myCub->getNumPoints(); // get number of cubature points 00127 00128 /* Related to basis. */ 00129 Basis_HDIV_HEX_I1_FEM<double, FieldContainer<double> > hexBasis; // create H-div basis on a hex 00130 int numFields = hexBasis.getCardinality(); // get basis cardinality 00131 00132 /* Cell geometries and orientations. */ 00133 int numCells = 4; 00134 int numNodes = 8; 00135 int numCellData = numCells*numNodes*spaceDim; 00136 int numSignData = numCells*numFields; 00137 00138 double hexnodes[] = { 00139 // hex 0 -- affine 00140 -1.0, -1.0, -1.0, 00141 1.0, -1.0, -1.0, 00142 1.0, 1.0, -1.0, 00143 -1.0, 1.0, -1.0, 00144 -1.0, -1.0, 1.0, 00145 1.0, -1.0, 1.0, 00146 1.0, 1.0, 1.0, 00147 -1.0, 1.0, 1.0, 00148 // hex 1 -- affine 00149 -3.0, -3.0, 1.0, 00150 6.0, 3.0, 1.0, 00151 7.0, 8.0, 0.0, 00152 -2.0, 2.0, 0.0, 00153 -3.0, -3.0, 4.0, 00154 6.0, 3.0, 4.0, 00155 7.0, 8.0, 3.0, 00156 -2.0, 2.0, 3.0, 00157 // hex 2 -- affine 00158 -3.0, -3.0, 0.0, 00159 9.0, 3.0, 0.0, 00160 15.0, 6.1, 0.0, 00161 3.0, 0.1, 0.0, 00162 9.0, 3.0, 0.1, 00163 21.0, 9.0, 0.1, 00164 27.0, 12.1, 0.1, 00165 15.0, 6.1, 0.1, 00166 // hex 3 -- nonaffine 00167 -2.0, -2.0, 0.0, 00168 2.0, -1.0, 0.0, 00169 1.0, 6.0, 0.0, 00170 -1.0, 1.0, 0.0, 00171 0.0, 0.0, 1.0, 00172 1.0, 0.0, 1.0, 00173 1.0, 1.0, 1.0, 00174 0.0, 1.0, 1.0 00175 }; 00176 00177 short facesigns[] = { 00178 1, 1, 1, 1, 1, 1, 00179 1, -1, 1, -1, 1, -1, 00180 -1, -1, 1, 1, -1, 1, 00181 -1, -1, 1, 1, -1, -1 00182 }; 00183 00184 /* Computational arrays. */ 00185 // First allocate one very large work space. 00186 Teuchos::Array<double> work_space(numCubPoints*spaceDim + 00187 numCubPoints + 00188 numCells*numNodes*spaceDim + 00189 numCells*numCubPoints*spaceDim*spaceDim + 00190 2*numCells*numCubPoints + 00191 numFields*numCubPoints + 00192 2*numCells*numFields*numCubPoints + 00193 numCells*numFields*numFields + 00194 numFields*numCubPoints*spaceDim + 00195 2*numCells*numFields*numCubPoints*spaceDim + 00196 numCells*numFields*numFields 00197 ); 00198 00199 int offset = 0; 00200 shards::Array<double,shards::NaturalOrder,Point,Dim> cub_points(&work_space[offset], numCubPoints, spaceDim); 00201 FC cub_points_FC(cub_points); 00202 offset += numCubPoints*spaceDim; 00203 shards::Array<double,shards::NaturalOrder,Point> cub_weights(&work_space[offset], numCubPoints); 00204 FC cub_weights_FC(cub_weights); 00205 offset += numCubPoints; 00206 shards::Array<double,shards::NaturalOrder,Cell,Node,Dim> cell_nodes(&work_space[offset], numCells, numNodes, spaceDim); 00207 FC cell_nodes_FC(cell_nodes); 00208 offset += numCells*numNodes*spaceDim; 00209 FieldContainer<short> field_signs(numCells, numFields); 00210 shards::Array<double,shards::NaturalOrder,Cell,Point,Dim,Dim> jacobian(&work_space[offset], numCells, numCubPoints, spaceDim, spaceDim); 00211 FC jacobian_FC(jacobian); 00212 offset += numCells*numCubPoints*spaceDim*spaceDim; 00213 //shards::Array<double,shards::NaturalOrder> jacobian_inv(&work_space[offset]numCells, numCubPoints, spaceDim, spaceDim); 00214 shards::Array<double,shards::NaturalOrder,Cell,Point> jacobian_det(&work_space[offset], numCells, numCubPoints); 00215 FC jacobian_det_FC(jacobian_det); 00216 offset += numCells*numCubPoints; 00217 shards::Array<double,shards::NaturalOrder,Cell,Point> weighted_measure(&work_space[offset], numCells, numCubPoints); 00218 FC weighted_measure_FC(weighted_measure); 00219 offset += numCells*numCubPoints; 00220 00221 shards::Array<double,shards::NaturalOrder,Field,Point> div_of_basis_at_cub_points(&work_space[offset], numFields, numCubPoints); 00222 FC div_of_basis_at_cub_points_FC(div_of_basis_at_cub_points); 00223 offset += numFields*numCubPoints; 00224 shards::Array<double,shards::NaturalOrder,Cell,Field,Point> 00225 transformed_div_of_basis_at_cub_points(&work_space[offset], numCells, numFields, numCubPoints); 00226 FC transformed_div_of_basis_at_cub_points_FC(transformed_div_of_basis_at_cub_points); 00227 offset += numCells*numFields*numCubPoints; 00228 shards::Array<double,shards::NaturalOrder,Cell,Field,Point> 00229 weighted_transformed_div_of_basis_at_cub_points(&work_space[offset], numCells, numFields, numCubPoints); 00230 FC weighted_transformed_div_of_basis_at_cub_points_FC(weighted_transformed_div_of_basis_at_cub_points); 00231 offset += numCells*numFields*numCubPoints; 00232 shards::Array<double,shards::NaturalOrder,Cell,Field,Field> stiffness_matrices(&work_space[offset], numCells, numFields, numFields); 00233 FC stiffness_matrices_FC(stiffness_matrices); 00234 offset += numCells*numFields*numFields; 00235 00236 shards::Array<double,shards::NaturalOrder,Field,Point,Dim> value_of_basis_at_cub_points(&work_space[offset], numFields, numCubPoints, spaceDim); 00237 FC value_of_basis_at_cub_points_FC(value_of_basis_at_cub_points); 00238 offset += numFields*numCubPoints*spaceDim; 00239 shards::Array<double,shards::NaturalOrder,Cell,Field,Point,Dim> 00240 transformed_value_of_basis_at_cub_points(&work_space[offset], numCells, numFields, numCubPoints, spaceDim); 00241 FC transformed_value_of_basis_at_cub_points_FC(transformed_value_of_basis_at_cub_points); 00242 offset += numCells*numFields*numCubPoints*spaceDim; 00243 shards::Array<double,shards::NaturalOrder,Cell,Field,Point,Dim> 00244 weighted_transformed_value_of_basis_at_cub_points(&work_space[offset], numCells, numFields, numCubPoints, spaceDim); 00245 FC weighted_transformed_value_of_basis_at_cub_points_FC(weighted_transformed_value_of_basis_at_cub_points); 00246 offset += numCells*numFields*numCubPoints*spaceDim; 00247 shards::Array<double,shards::NaturalOrder,Cell,Field,Field> mass_matrices(&work_space[offset], numCells, numFields, numFields); 00248 FC mass_matrices_FC(mass_matrices); 00249 00250 00251 /******************* START COMPUTATION ***********************/ 00252 00253 // get cubature points and weights 00254 myCub->getCubature(cub_points_FC, cub_weights_FC); 00255 00256 // fill cell vertex array 00257 cell_nodes_FC.setValues(hexnodes, numCellData); 00258 00259 // set basis function signs, for each cell 00260 field_signs.setValues(facesigns, numSignData); 00261 00262 // compute geometric cell information 00263 CellTools<double>::setJacobian(jacobian_FC, cub_points_FC, cell_nodes_FC, cellType); 00264 //CellTools<double>::setJacobianInv(jacobian_inv, jacobian); 00265 CellTools<double>::setJacobianDet(jacobian_det_FC, jacobian_FC); 00266 00267 // compute weighted measure 00268 fst::computeCellMeasure<double>(weighted_measure_FC, jacobian_det_FC, cub_weights_FC); 00269 00270 00271 // Computing stiffness matrices: 00272 // tabulate divergences of basis functions at (reference) cubature points 00273 hexBasis.getValues(div_of_basis_at_cub_points_FC, cub_points_FC, OPERATOR_DIV); 00274 00275 // transform divergences of basis functions 00276 fst::HDIVtransformDIV<double>(transformed_div_of_basis_at_cub_points_FC, 00277 jacobian_det_FC, 00278 div_of_basis_at_cub_points_FC); 00279 00280 // multiply with weighted measure 00281 fst::multiplyMeasure<double>(weighted_transformed_div_of_basis_at_cub_points_FC, 00282 weighted_measure_FC, 00283 transformed_div_of_basis_at_cub_points_FC); 00284 00285 // compute stiffness matrices 00286 fst::integrate<double>(stiffness_matrices_FC, 00287 transformed_div_of_basis_at_cub_points_FC, 00288 weighted_transformed_div_of_basis_at_cub_points_FC, 00289 COMP_CPP); 00290 00291 // apply field signs 00292 fst::applyLeftFieldSigns<double>(stiffness_matrices_FC, field_signs); 00293 fst::applyRightFieldSigns<double>(stiffness_matrices_FC, field_signs); 00294 00295 00296 // Computing mass matrices: 00297 // tabulate values of basis functions at (reference) cubature points 00298 hexBasis.getValues(value_of_basis_at_cub_points_FC, cub_points_FC, OPERATOR_VALUE); 00299 00300 // transform values of basis functions 00301 fst::HDIVtransformVALUE<double>(transformed_value_of_basis_at_cub_points_FC, 00302 jacobian_FC, 00303 jacobian_det_FC, 00304 value_of_basis_at_cub_points_FC); 00305 00306 // multiply with weighted measure 00307 fst::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points_FC, 00308 weighted_measure_FC, 00309 transformed_value_of_basis_at_cub_points_FC); 00310 00311 // compute mass matrices 00312 fst::integrate<double>(mass_matrices_FC, 00313 transformed_value_of_basis_at_cub_points_FC, 00314 weighted_transformed_value_of_basis_at_cub_points_FC, 00315 COMP_CPP); 00316 00317 // apply field signs 00318 fst::applyLeftFieldSigns<double>(mass_matrices_FC, field_signs); 00319 fst::applyRightFieldSigns<double>(mass_matrices_FC, field_signs); 00320 00321 /******************* STOP COMPUTATION ***********************/ 00322 00323 00324 /******************* START COMPARISON ***********************/ 00325 string basedir = "./testdata"; 00326 for (int cell_id = 0; cell_id < numCells-1; cell_id++) { 00327 00328 stringstream namestream; 00329 string filename; 00330 namestream << basedir << "/mass_HDIV_HEX_I1_FEM" << "_" << "0" << cell_id+1 << ".dat"; 00331 namestream >> filename; 00332 00333 ifstream massfile(&filename[0]); 00334 if (massfile.is_open()) { 00335 if (compareToAnalytic<double>(&mass_matrices(cell_id, 0, 0), massfile, 1e-10, iprint) > 0) 00336 errorFlag++; 00337 massfile.close(); 00338 } 00339 else { 00340 errorFlag = -1; 00341 std::cout << "End Result: TEST FAILED\n"; 00342 return errorFlag; 00343 } 00344 00345 namestream.clear(); 00346 namestream << basedir << "/stiff_HDIV_HEX_I1_FEM" << "_" << "0" << cell_id+1 << ".dat"; 00347 namestream >> filename; 00348 00349 ifstream stifffile(&filename[0]); 00350 if (stifffile.is_open()) 00351 { 00352 if (compareToAnalytic<double>(&stiffness_matrices(cell_id, 0, 0), stifffile, 1e-10, iprint) > 0) 00353 errorFlag++; 00354 stifffile.close(); 00355 } 00356 else { 00357 errorFlag = -1; 00358 std::cout << "End Result: TEST FAILED\n"; 00359 return errorFlag; 00360 } 00361 00362 } 00363 for (int cell_id = 3; cell_id < numCells; cell_id++) { 00364 00365 stringstream namestream; 00366 string filename; 00367 namestream << basedir << "/mass_fp_HDIV_HEX_I1_FEM" << "_" << "0" << cell_id+1 << ".dat"; 00368 namestream >> filename; 00369 00370 ifstream massfile(&filename[0]); 00371 if (massfile.is_open()) { 00372 if (compareToAnalytic<double>(&mass_matrices(cell_id, 0, 0), massfile, 1e-4, iprint, INTREPID_UTILS_SCALAR) > 0) 00373 errorFlag++; 00374 massfile.close(); 00375 } 00376 else { 00377 errorFlag = -1; 00378 std::cout << "End Result: TEST FAILED\n"; 00379 return errorFlag; 00380 } 00381 00382 namestream.clear(); 00383 namestream << basedir << "/stiff_fp_HDIV_HEX_I1_FEM" << "_" << "0" << cell_id+1 << ".dat"; 00384 namestream >> filename; 00385 00386 ifstream stifffile(&filename[0]); 00387 if (stifffile.is_open()) 00388 { 00389 if (compareToAnalytic<double>(&stiffness_matrices(cell_id, 0, 0), stifffile, 1e-4, iprint, INTREPID_UTILS_SCALAR) > 0) 00390 errorFlag++; 00391 stifffile.close(); 00392 } 00393 else { 00394 errorFlag = -1; 00395 std::cout << "End Result: TEST FAILED\n"; 00396 return errorFlag; 00397 } 00398 00399 } 00400 /******************* STOP COMPARISON ***********************/ 00401 00402 *outStream << "\n"; 00403 } 00404 catch (std::logic_error err) { 00405 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n"; 00406 *outStream << err.what() << '\n'; 00407 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00408 errorFlag = -1000; 00409 }; 00410 00411 00412 if (errorFlag != 0) 00413 std::cout << "End Result: TEST FAILED\n"; 00414 else 00415 std::cout << "End Result: TEST PASSED\n"; 00416 00417 // reset format state of std::cout 00418 std::cout.copyfmt(oldFormatState); 00419 00420 return errorFlag; 00421 }
1.7.4