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