|
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_HCURL_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 HCURL |\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 shards::CellTopology cellType = shards::getCellTopologyData< shards::Hexahedron<> >(); // cell type: hex 00096 00097 /* Related to cubature. */ 00098 DefaultCubatureFactory<double> cubFactory; // create cubature factory 00099 int cubDegree = 20; // cubature degree 00100 Teuchos::RCP<Cubature<double> > myCub = cubFactory.create(cellType, cubDegree); // create default cubature 00101 int spaceDim = myCub->getDimension(); // get spatial dimension 00102 int numCubPoints = myCub->getNumPoints(); // get number of cubature points 00103 00104 /* Related to basis. */ 00105 Basis_HCURL_HEX_I1_FEM<double, FieldContainer<double> > hexBasis; // create H-curl basis on a hex 00106 int numFields = hexBasis.getCardinality(); // get basis cardinality 00107 00108 /* Cell geometries and orientations. */ 00109 int numCells = 4; 00110 int numNodes = 8; 00111 int numCellData = numCells*numNodes*spaceDim; 00112 int numSignData = numCells*numFields; 00113 00114 double hexnodes[] = { 00115 // hex 0 -- affine 00116 -1.0, -1.0, -1.0, 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 // hex 1 -- affine 00125 -3.0, -3.0, 1.0, 00126 6.0, 3.0, 1.0, 00127 7.0, 8.0, 0.0, 00128 -2.0, 2.0, 0.0, 00129 -3.0, -3.0, 4.0, 00130 6.0, 3.0, 4.0, 00131 7.0, 8.0, 3.0, 00132 -2.0, 2.0, 3.0, 00133 // hex 2 -- affine 00134 -3.0, -3.0, 0.0, 00135 9.0, 3.0, 0.0, 00136 15.0, 6.1, 0.0, 00137 3.0, 0.1, 0.0, 00138 9.0, 3.0, 0.1, 00139 21.0, 9.0, 0.1, 00140 27.0, 12.1, 0.1, 00141 15.0, 6.1, 0.1, 00142 // hex 3 -- nonaffine 00143 -2.0, -2.0, 0.0, 00144 2.0, -1.0, 0.0, 00145 1.0, 6.0, 0.0, 00146 -1.0, 1.0, 0.0, 00147 0.0, 0.0, 1.0, 00148 1.0, 0.0, 1.0, 00149 1.0, 1.0, 1.0, 00150 0.0, 1.0, 1.0 00151 }; 00152 00153 short edgesigns[] = { 00154 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 00155 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 1, -1, 00156 -1, -1, -1, 1, 1, 1, -1, -1, 1, 1, -1, 1, 00157 1, 1, 1, 1, -1, -1, -1, -1, 1, 1, 1, 1 00158 }; 00159 00160 /* Computational arrays. */ 00161 FieldContainer<double> cub_points(numCubPoints, spaceDim); 00162 FieldContainer<double> cub_weights(numCubPoints); 00163 FieldContainer<double> cell_nodes(numCells, numNodes, spaceDim); 00164 FieldContainer<short> field_signs(numCells, numFields); 00165 FieldContainer<double> jacobian(numCells, numCubPoints, spaceDim, spaceDim); 00166 FieldContainer<double> jacobian_inv(numCells, numCubPoints, spaceDim, spaceDim); 00167 FieldContainer<double> jacobian_det(numCells, numCubPoints); 00168 FieldContainer<double> weighted_measure(numCells, numCubPoints); 00169 00170 FieldContainer<double> curl_of_basis_at_cub_points(numFields, numCubPoints, spaceDim); 00171 FieldContainer<double> transformed_curl_of_basis_at_cub_points(numCells, numFields, numCubPoints, spaceDim); 00172 FieldContainer<double> weighted_transformed_curl_of_basis_at_cub_points(numCells, numFields, numCubPoints, spaceDim); 00173 FieldContainer<double> stiffness_matrices(numCells, numFields, numFields); 00174 00175 FieldContainer<double> value_of_basis_at_cub_points(numFields, numCubPoints, spaceDim); 00176 FieldContainer<double> transformed_value_of_basis_at_cub_points(numCells, numFields, numCubPoints, spaceDim); 00177 FieldContainer<double> weighted_transformed_value_of_basis_at_cub_points(numCells, numFields, numCubPoints, spaceDim); 00178 FieldContainer<double> mass_matrices(numCells, numFields, numFields); 00179 00180 /******************* START COMPUTATION ***********************/ 00181 00182 // get cubature points and weights 00183 myCub->getCubature(cub_points, cub_weights); 00184 00185 // fill cell vertex array 00186 cell_nodes.setValues(hexnodes, numCellData); 00187 00188 // set basis function signs, for each cell 00189 field_signs.setValues(edgesigns, numSignData); 00190 00191 // compute geometric cell information 00192 CellTools<double>::setJacobian(jacobian, cub_points, cell_nodes, cellType); 00193 CellTools<double>::setJacobianInv(jacobian_inv, jacobian); 00194 CellTools<double>::setJacobianDet(jacobian_det, jacobian); 00195 00196 // compute weighted measure 00197 fst::computeCellMeasure<double>(weighted_measure, jacobian_det, cub_weights); 00198 00199 // Computing stiffness matrices: 00200 // tabulate curls of basis functions at (reference) cubature points 00201 hexBasis.getValues(curl_of_basis_at_cub_points, cub_points, OPERATOR_CURL); 00202 00203 // transform curls of basis functions 00204 fst::HCURLtransformCURL<double>(transformed_curl_of_basis_at_cub_points, 00205 jacobian, 00206 jacobian_det, 00207 curl_of_basis_at_cub_points); 00208 00209 // multiply with weighted measure 00210 fst::multiplyMeasure<double>(weighted_transformed_curl_of_basis_at_cub_points, 00211 weighted_measure, 00212 transformed_curl_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_curl_of_basis_at_cub_points, field_signs); 00216 fst::applyFieldSigns<double>(weighted_transformed_curl_of_basis_at_cub_points, field_signs); 00217 00218 // compute stiffness matrices 00219 fst::integrate<double>(stiffness_matrices, 00220 transformed_curl_of_basis_at_cub_points, 00221 weighted_transformed_curl_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::HCURLtransformVALUE<double>(transformed_value_of_basis_at_cub_points, 00230 jacobian_inv, 00231 value_of_basis_at_cub_points); 00232 00233 // multiply with weighted measure 00234 fst::multiplyMeasure<double>(weighted_transformed_value_of_basis_at_cub_points, 00235 weighted_measure, 00236 transformed_value_of_basis_at_cub_points); 00237 00238 // compute mass matrices 00239 fst::integrate<double>(mass_matrices, 00240 transformed_value_of_basis_at_cub_points, 00241 weighted_transformed_value_of_basis_at_cub_points, 00242 COMP_CPP); 00243 00244 // apply field signs (after the fact, as a post-processing step) 00245 fst::applyLeftFieldSigns<double>(mass_matrices, field_signs); 00246 fst::applyRightFieldSigns<double>(mass_matrices, field_signs); 00247 00248 /******************* STOP COMPUTATION ***********************/ 00249 00250 00251 /******************* START COMPARISON ***********************/ 00252 string basedir = "./testdata"; 00253 for (int cell_id = 0; cell_id < numCells-1; cell_id++) { 00254 00255 stringstream namestream; 00256 string filename; 00257 namestream << basedir << "/mass_HCURL_HEX_I1_FEM" << "_" << "0" << cell_id+1 << ".dat"; 00258 namestream >> filename; 00259 00260 ifstream massfile(&filename[0]); 00261 if (massfile.is_open()) { 00262 if (compareToAnalytic<double>(&mass_matrices(cell_id, 0, 0), massfile, 1e-10, iprint) > 0) 00263 errorFlag++; 00264 massfile.close(); 00265 } 00266 else { 00267 errorFlag = -1; 00268 std::cout << "End Result: TEST FAILED\n"; 00269 return errorFlag; 00270 } 00271 00272 namestream.clear(); 00273 namestream << basedir << "/stiff_HCURL_HEX_I1_FEM" << "_" << "0" << cell_id+1 << ".dat"; 00274 namestream >> filename; 00275 00276 ifstream stifffile(&filename[0]); 00277 if (stifffile.is_open()) 00278 { 00279 if (compareToAnalytic<double>(&stiffness_matrices(cell_id, 0, 0), stifffile, 1e-10, iprint) > 0) 00280 errorFlag++; 00281 stifffile.close(); 00282 } 00283 else { 00284 errorFlag = -1; 00285 std::cout << "End Result: TEST FAILED\n"; 00286 return errorFlag; 00287 } 00288 00289 } 00290 for (int cell_id = 3; cell_id < numCells; cell_id++) { 00291 00292 stringstream namestream; 00293 string filename; 00294 namestream << basedir << "/mass_fp_HCURL_HEX_I1_FEM" << "_" << "0" << cell_id+1 << ".dat"; 00295 namestream >> filename; 00296 00297 ifstream massfile(&filename[0]); 00298 if (massfile.is_open()) { 00299 if (compareToAnalytic<double>(&mass_matrices(cell_id, 0, 0), massfile, 1e-4, iprint, INTREPID_UTILS_SCALAR) > 0) 00300 errorFlag++; 00301 massfile.close(); 00302 } 00303 else { 00304 errorFlag = -1; 00305 std::cout << "End Result: TEST FAILED\n"; 00306 return errorFlag; 00307 } 00308 00309 namestream.clear(); 00310 namestream << basedir << "/stiff_fp_HCURL_HEX_I1_FEM" << "_" << "0" << cell_id+1 << ".dat"; 00311 namestream >> filename; 00312 00313 ifstream stifffile(&filename[0]); 00314 if (stifffile.is_open()) 00315 { 00316 if (compareToAnalytic<double>(&stiffness_matrices(cell_id, 0, 0), stifffile, 1e-4, iprint, INTREPID_UTILS_SCALAR) > 0) 00317 errorFlag++; 00318 stifffile.close(); 00319 } 00320 else { 00321 errorFlag = -1; 00322 std::cout << "End Result: TEST FAILED\n"; 00323 return errorFlag; 00324 } 00325 00326 } 00327 /******************* STOP COMPARISON ***********************/ 00328 00329 *outStream << "\n"; 00330 } 00331 catch (std::logic_error err) { 00332 *outStream << "UNEXPECTED ERROR !!! ----------------------------------------------------------\n"; 00333 *outStream << err.what() << '\n'; 00334 *outStream << "-------------------------------------------------------------------------------" << "\n\n"; 00335 errorFlag = -1000; 00336 }; 00337 00338 00339 if (errorFlag != 0) 00340 std::cout << "End Result: TEST FAILED\n"; 00341 else 00342 std::cout << "End Result: TEST PASSED\n"; 00343 00344 // reset format state of std::cout 00345 std::cout.copyfmt(oldFormatState); 00346 00347 return errorFlag; 00348 }
1.7.4