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