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