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