|
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), 00025 // Denis Ridzal (dridzal@sandia.gov), 00026 // Kara Peterson (kjpeter@sandia.gov). 00027 // 00028 // ************************************************************************ 00029 // @HEADER 00030 00036 #include "Intrepid_FieldContainer.hpp" 00037 #include "Intrepid_HCURL_TET_In_FEM.hpp" 00038 #include "Intrepid_DefaultCubatureFactory.hpp" 00039 #include "Intrepid_RealSpaceTools.hpp" 00040 #include "Intrepid_ArrayTools.hpp" 00041 #include "Intrepid_FunctionSpaceTools.hpp" 00042 #include "Intrepid_CellTools.hpp" 00043 #include "Teuchos_oblackholestream.hpp" 00044 #include "Teuchos_RCP.hpp" 00045 #include "Teuchos_GlobalMPISession.hpp" 00046 #include "Teuchos_SerialDenseMatrix.hpp" 00047 #include "Teuchos_SerialDenseVector.hpp" 00048 #include "Teuchos_LAPACK.hpp" 00049 00050 using namespace std; 00051 using namespace Intrepid; 00052 00053 void rhsFunc( FieldContainer<double> &, const FieldContainer<double> &, int, int, int, int ); 00054 void bcFunc( FieldContainer<double> &, const FieldContainer<double> & , 00055 const shards::CellTopology & , 00056 int , int , int , int , int ); 00057 void u_exact( FieldContainer<double> &, const FieldContainer<double> &, int , int, int, int ); 00058 00059 void u_exact( FieldContainer<double> &result, 00060 const FieldContainer<double> &points, 00061 int comp, 00062 int xd, 00063 int yd, 00064 int zd) 00065 { 00066 for (int cell=0;cell<result.dimension(0);cell++){ 00067 for (int pt=0;pt<result.dimension(1);pt++) { 00068 result(cell,pt,comp) = std::pow(points(cell,pt,0),xd)*std::pow(points(cell,pt,1),yd) 00069 *std::pow(points(cell,pt,2),zd); 00070 } 00071 } 00072 return; 00073 } 00074 00075 void bcFunc( FieldContainer<double> & result, 00076 const FieldContainer<double> & points , 00077 const shards::CellTopology & parentCell , 00078 int sideOrdinal , int comp , int a, int b, int c ) 00079 { 00080 00081 int numCells = result.dimension(0); 00082 int numPoints = result.dimension(1); 00083 00084 // reference face normal 00085 FieldContainer<double> normal(3); 00086 CellTools<double>::getReferenceFaceNormal(normal,sideOrdinal,parentCell); 00087 00088 result.initialize(); 00089 00090 if (comp == 0) { 00091 for (int cell=0;cell<numCells;cell++) { 00092 for (int pt=0;pt<numPoints;pt++) { 00093 // first component 00094 if (c > 0) { 00095 result(cell,pt,0) -= c * normal(2) * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c-1); 00096 } 00097 if (b > 0) { 00098 result(cell,pt,0) -= b * normal(1) * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b-1) * std::pow(points(cell,pt,2),c); 00099 } 00100 // second component 00101 if (b > 0) { 00102 result(cell,pt,1) = b * normal(0) * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b-1) * std::pow(points(cell,pt,2),c); 00103 } 00104 // third component 00105 if (c > 0) { 00106 result(cell,pt,2) = c * normal(0) * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c-1); 00107 } 00108 } 00109 } 00110 } 00111 else if (comp == 1) { 00112 for (int cell=0;cell<numCells;cell++) { 00113 for (int pt=0;pt<numPoints;pt++) { 00114 // first component 00115 if (a > 0) { 00116 result(cell,pt,0) = a * normal(1) * std::pow(points(cell,pt,0),a-1) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c); 00117 } 00118 // second component 00119 if (c > 0) { 00120 result(cell,pt,1) -= c * normal(2) * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c-1); 00121 } 00122 if (a > 0) { 00123 result(cell,pt,1) -= a * normal(0) * std::pow(points(cell,pt,0),a-1) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c); 00124 } 00125 // third component 00126 if (c > 0) { 00127 result(cell,pt,2) = c * normal(1) * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c-1); 00128 } 00129 } 00130 } 00131 } 00132 else if (comp == 2) { 00133 for (int cell=0;cell<numCells;cell++) { 00134 for (int pt=0;pt<numPoints;pt++) { 00135 // first component 00136 if (a > 0) { 00137 result(cell,pt,0) = a * normal(2) * std::pow(points(cell,pt,0),a-1) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c); 00138 } 00139 // second component 00140 if (b > 0) { 00141 result(cell,pt,1) = b * normal(2) * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b-1) * std::pow(points(cell,pt,2),c); 00142 } 00143 // third component 00144 if (b > 0) { 00145 result(cell,pt,2) -= b * normal(1) * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b-1)* std::pow(points(cell,pt,2),c); 00146 } 00147 if (a > 0) { 00148 result(cell,pt,2) -= a * normal(0) * std::pow(points(cell,pt,0),a-1) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c); 00149 } 00150 } 00151 } 00152 } 00153 00154 return; 00155 } 00156 00157 void rhsFunc( FieldContainer<double> & result , 00158 const FieldContainer<double> &points , 00159 int comp, 00160 int a, 00161 int b, 00162 int c ) 00163 { 00164 // rhs is curl(curl(E))+E, so load up the exact solution 00165 00166 u_exact(result,points,comp,a,b,c); 00167 00168 if (comp == 0) { 00169 // component 0 00170 if (b>=2) { 00171 for (int cell=0;cell<result.dimension(0);cell++) { 00172 for (int pt=0;pt<result.dimension(1);pt++) { 00173 result(cell,pt,0) -= b*(b-1)*std::pow(points(cell,pt,0),a)*std::pow(points(cell,pt,1),b-2)*std::pow(points(cell,pt,2),c); 00174 } 00175 } 00176 } 00177 if (c>=2) { 00178 for (int cell=0;cell<result.dimension(0);cell++) { 00179 for (int pt=0;pt<result.dimension(1);pt++) { 00180 result(cell,pt,0)-=c*(c-1)*std::pow(points(cell,pt,0),a)*std::pow(points(cell,pt,1),b)*std::pow(points(cell,pt,2),c-2); 00181 } 00182 } 00183 } 00184 // component one 00185 if ( (a>=1) && (b>=1) ) { 00186 for (int cell=0;cell<result.dimension(0);cell++) { 00187 for (int pt=0;pt<result.dimension(1);pt++) { 00188 result(cell,pt,1)+=a*b*std::pow(points(cell,pt,0),a-1)*std::pow(points(cell,pt,1),b-1)*std::pow(points(cell,pt,2),c); 00189 } 00190 } 00191 } 00192 // component two 00193 if ( (a>=1) && (c>=1) ) { 00194 for (int cell=0;cell<result.dimension(0);cell++) { 00195 for (int pt=0;pt<result.dimension(1);pt++) { 00196 result(cell,pt,2)+=a*c*std::pow(points(cell,pt,0),a-1)*std::pow(points(cell,pt,1),b)*std::pow(points(cell,pt,2),c-1); 00197 } 00198 } 00199 } 00200 } 00201 if (comp == 1) { 00202 for (int cell=0;cell<result.dimension(0);cell++) { 00203 for (int pt=0;pt<result.dimension(1);pt++) { 00204 // first component 00205 if (a > 0 && b > 0) { 00206 result(cell,pt,0) += a * b * std::pow(points(cell,pt,0),a-1) * std::pow(points(cell,pt,1),b-1) * std::pow(points(cell,pt,2),c); 00207 } 00208 // second component 00209 if (c > 1) { 00210 result(cell,pt,1) -= c*(c-1.0)*std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b)*std::pow(points(cell,pt,2),c-2); 00211 } 00212 if (a > 1) { 00213 result(cell,pt,1) -= a*(a-1.0)*std::pow(points(cell,pt,0),a-2) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c); 00214 } 00215 // third component 00216 if (b > 0 && c > 0) { 00217 result(cell,pt,2) += b * c * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b-1) * std::pow(points(cell,pt,2),c-1); 00218 } 00219 } 00220 } 00221 } 00222 else if (comp == 2) { 00223 for (int cell=0;cell<result.dimension(0);cell++) { 00224 for (int pt=0;pt<result.dimension(1);pt++) { 00225 // first component 00226 if (a>0 && c>0) { 00227 result(cell,pt,0) += a * c * std::pow(points(cell,pt,0),a-1) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c-1); 00228 } 00229 // second component 00230 if (b>0 && c>0) { 00231 result(cell,pt,1) += b * c * std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b-1) * std::pow(points(cell,pt,2),c-1); 00232 } 00233 // third component 00234 if (b>1) { 00235 result(cell,pt,2) -= b*(b-1.0)*std::pow(points(cell,pt,0),a) * std::pow(points(cell,pt,1),b-2) * std::pow(points(cell,pt,2),c); 00236 } 00237 if (a>1) { 00238 result(cell,pt,2) -= a*(a-1.0)*std::pow(points(cell,pt,0),a-2) * std::pow(points(cell,pt,1),b) * std::pow(points(cell,pt,2),c); 00239 } 00240 } 00241 } 00242 } 00243 return; 00244 } 00245 00246 int main(int argc, char *argv[]) { 00247 Teuchos::GlobalMPISession mpiSession(&argc, &argv); 00248 00249 // This little trick lets us print to std::cout only if 00250 // a (dummy) command-line argument is provided. 00251 int iprint = argc - 1; 00252 Teuchos::RCP<std::ostream> outStream; 00253 Teuchos::oblackholestream bhs; // outputs nothing 00254 if (iprint > 0) 00255 outStream = Teuchos::rcp(&std::cout, false); 00256 else 00257 outStream = Teuchos::rcp(&bhs, false); 00258 00259 // Save the format state of the original std::cout. 00260 Teuchos::oblackholestream oldFormatState; 00261 oldFormatState.copyfmt(std::cout); 00262 00263 *outStream \ 00264 << "===============================================================================\n" \ 00265 << "| |\n" \ 00266 << "| Unit Test (Basis_HGRAD_TRI_In_FEM) |\n" \ 00267 << "| |\n" \ 00268 << "| 1) Patch test involving H(curl) matrices |\n" \ 00269 << "| |\n" \ 00270 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \ 00271 << "| Robert Kirby (robert.c.kirby@ttu.edu), |\n" \ 00272 << "| Denis Ridzal (dridzal@sandia.gov), |\n" \ 00273 << "| Kara Peterson (kjpeter@sandia.gov). |\n" \ 00274 << "| |\n" \ 00275 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \ 00276 << "| Trilinos website: http://trilinos.sandia.gov |\n" \ 00277 << "| |\n" \ 00278 << "===============================================================================\n" \ 00279 << "| TEST 1: Patch test for mass + curl-curl matrices |\n" \ 00280 << "===============================================================================\n"; 00281 00282 00283 int errorFlag = 0; 00284 00285 outStream -> precision(16); 00286 00287 try { 00288 DefaultCubatureFactory<double> cubFactory; // create cubature factory 00289 shards::CellTopology cell(shards::getCellTopologyData< shards::Tetrahedron<> >()); // create parent cell topology 00290 00291 shards::CellTopology side(shards::getCellTopologyData< shards::Triangle<> >()); 00292 00293 int cellDim = cell.getDimension(); 00294 int sideDim = side.getDimension(); 00295 00296 int min_order = 1; 00297 int max_order = 4; 00298 00299 int numIntervals = max_order; 00300 int numInterpPoints = ((numIntervals + 1)*(numIntervals + 2)*(numIntervals+3))/6; 00301 FieldContainer<double> interp_points_ref(numInterpPoints, cellDim); 00302 int counter = 0; 00303 for (int j=0; j<=numIntervals; j++) { 00304 for (int i=0; i<=numIntervals-j; i++) { 00305 for (int k=0;k<numIntervals-j-i;k++) { 00306 interp_points_ref(counter,0) = i*(1.0/numIntervals); 00307 interp_points_ref(counter,1) = j*(1.0/numIntervals); 00308 interp_points_ref(counter,2) = k*(1.0/numIntervals); 00309 counter++; 00310 } 00311 } 00312 } 00313 00314 for (int basis_order=min_order;basis_order<=max_order;basis_order++) { 00315 // create basis 00316 Teuchos::RCP<Basis<double,FieldContainer<double> > > basis = 00317 Teuchos::rcp(new Basis_HCURL_TET_In_FEM<double,FieldContainer<double> >(basis_order,POINTTYPE_EQUISPACED) ); 00318 00319 int numFields = basis->getCardinality(); 00320 00321 // create cubatures 00322 Teuchos::RCP<Cubature<double> > cellCub = cubFactory.create(cell, 2*(basis_order+1)); 00323 Teuchos::RCP<Cubature<double> > sideCub = cubFactory.create(side, 2*(basis_order+1)); 00324 00325 int numCubPointsCell = cellCub->getNumPoints(); 00326 int numCubPointsSide = sideCub->getNumPoints(); 00327 00328 // hold cubature information 00329 FieldContainer<double> cub_points_cell(numCubPointsCell, cellDim); 00330 FieldContainer<double> cub_weights_cell(numCubPointsCell); 00331 FieldContainer<double> cub_points_side(numCubPointsSide,sideDim); 00332 FieldContainer<double> cub_weights_side(numCubPointsSide); 00333 FieldContainer<double> cub_points_side_refcell(numCubPointsSide,cellDim); 00334 00335 // hold basis function information on refcell 00336 FieldContainer<double> value_of_basis_at_cub_points_cell(numFields, numCubPointsCell, cellDim ); 00337 FieldContainer<double> w_value_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell, cellDim); 00338 FieldContainer<double> curl_of_basis_at_cub_points_cell(numFields, numCubPointsCell, cellDim ); 00339 FieldContainer<double> w_curl_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell, cellDim); 00340 FieldContainer<double> value_of_basis_at_cub_points_side_refcell(numFields,numCubPointsSide,cellDim ); 00341 FieldContainer<double> w_value_of_basis_at_cub_points_side_refcell(1,numFields,numCubPointsSide,cellDim ); 00342 00343 00344 // holds rhs data 00345 FieldContainer<double> rhs_at_cub_points_cell(1,numCubPointsCell,cellDim); 00346 FieldContainer<double> bc_func_at_cub_points_side_refcell(1,numCubPointsSide,cellDim); 00347 FieldContainer<double> bc_fields_per_side(1,numFields); 00348 00349 // FEM mass matrix 00350 FieldContainer<double> fe_matrix_bak(1,numFields,numFields); 00351 FieldContainer<double> fe_matrix(1,numFields,numFields); 00352 FieldContainer<double> rhs_and_soln_vec(1,numFields); 00353 00354 FieldContainer<double> value_of_basis_at_interp_points( numFields , numInterpPoints , cellDim); 00355 FieldContainer<double> interpolant( 1, numInterpPoints , cellDim ); 00356 00357 int info = 0; 00358 Teuchos::LAPACK<int, double> solver; 00359 00360 // set test tolerance 00361 double zero = (basis_order+1)*(basis_order+1)*1000*INTREPID_TOL; 00362 00363 // build matrices outside the loop, and then just do the rhs 00364 // for each iteration 00365 cellCub->getCubature(cub_points_cell, cub_weights_cell); 00366 sideCub->getCubature(cub_points_side, cub_weights_side); 00367 00368 // need the vector basis 00369 basis->getValues(value_of_basis_at_cub_points_cell, 00370 cub_points_cell, 00371 OPERATOR_VALUE); 00372 basis->getValues(curl_of_basis_at_cub_points_cell, 00373 cub_points_cell, 00374 OPERATOR_CURL); 00375 00376 basis->getValues( value_of_basis_at_interp_points , 00377 interp_points_ref , 00378 OPERATOR_VALUE ); 00379 00380 00381 // construct mass and curl-curl matrices 00382 cub_weights_cell.resize(1,numCubPointsCell); 00383 FunctionSpaceTools::multiplyMeasure<double>(w_value_of_basis_at_cub_points_cell , 00384 cub_weights_cell , 00385 value_of_basis_at_cub_points_cell ); 00386 FunctionSpaceTools::multiplyMeasure<double>(w_curl_of_basis_at_cub_points_cell , 00387 cub_weights_cell , 00388 curl_of_basis_at_cub_points_cell ); 00389 cub_weights_cell.resize(numCubPointsCell); 00390 00391 00392 value_of_basis_at_cub_points_cell.resize( 1 , numFields , numCubPointsCell , cellDim ); 00393 curl_of_basis_at_cub_points_cell.resize( 1 , numFields , numCubPointsCell , cellDim ); 00394 FunctionSpaceTools::integrate<double>(fe_matrix_bak, 00395 w_value_of_basis_at_cub_points_cell , 00396 value_of_basis_at_cub_points_cell , 00397 COMP_BLAS ); 00398 FunctionSpaceTools::integrate<double>(fe_matrix_bak, 00399 w_curl_of_basis_at_cub_points_cell , 00400 curl_of_basis_at_cub_points_cell , 00401 COMP_BLAS , 00402 true ); 00403 value_of_basis_at_cub_points_cell.resize( numFields , numCubPointsCell , cellDim ); 00404 curl_of_basis_at_cub_points_cell.resize( numFields , numCubPointsCell , cellDim ); 00405 00406 for (int comp=0;comp<cellDim;comp++) { 00407 for (int x_order=0;x_order<basis_order;x_order++) { 00408 for (int y_order=0;y_order<basis_order-x_order;y_order++) { 00409 for (int z_order=0;z_order<basis_order-x_order-y_order;z_order++) { 00410 fe_matrix.initialize(); 00411 // copy mass matrix 00412 for (int i=0;i<numFields;i++) { 00413 for (int j=0;j<numFields;j++) { 00414 fe_matrix(0,i,j) = fe_matrix_bak(0,i,j); 00415 } 00416 } 00417 00418 // clear old vector data 00419 rhs_and_soln_vec.initialize(); 00420 00421 // now get rhs vector 00422 cub_points_cell.resize(1,numCubPointsCell,cellDim); 00423 00424 rhs_at_cub_points_cell.initialize(); 00425 rhsFunc(rhs_at_cub_points_cell, 00426 cub_points_cell, 00427 comp, 00428 x_order, 00429 y_order, 00430 z_order); 00431 00432 cub_points_cell.resize(numCubPointsCell,cellDim); 00433 00434 cub_weights_cell.resize(numCubPointsCell); 00435 00436 FunctionSpaceTools::integrate<double>(rhs_and_soln_vec, 00437 rhs_at_cub_points_cell, 00438 w_value_of_basis_at_cub_points_cell, 00439 COMP_BLAS); 00440 00441 // now I need to get the boundary condition terms in place 00442 00443 for (unsigned i=0;i<4;i++) { 00444 // map side quadrature to reference cell side 00445 CellTools<double>::mapToReferenceSubcell(cub_points_side_refcell,cub_points_side,sideDim, 00446 (int) i, cell); 00447 00448 //evaluate basis at these points 00449 basis->getValues( value_of_basis_at_cub_points_side_refcell , 00450 cub_points_side_refcell , 00451 OPERATOR_VALUE ); 00452 00453 // evaluate imposed current on surface, which is n x curl(u_exact), on the quad points 00454 cub_points_side_refcell.resize( 1 , numCubPointsSide , cellDim ); 00455 bcFunc(bc_func_at_cub_points_side_refcell, 00456 cub_points_side_refcell, 00457 cell , 00458 i, 00459 comp , 00460 x_order, 00461 y_order, 00462 z_order); 00463 cub_points_side_refcell.resize( numCubPointsSide , cellDim ); 00464 00465 // now I need to integrate the bc function against the test basis 00466 // need to weight the basis functions with quadrature weights 00467 // the size of the face is embedded in the normal 00468 cub_weights_side.resize(1,numCubPointsSide); 00469 FunctionSpaceTools::multiplyMeasure<double>(w_value_of_basis_at_cub_points_side_refcell , 00470 cub_weights_side , 00471 value_of_basis_at_cub_points_side_refcell ); 00472 cub_weights_side.resize(numCubPointsSide); 00473 00474 FunctionSpaceTools::integrate<double>( bc_fields_per_side , 00475 bc_func_at_cub_points_side_refcell , 00476 w_value_of_basis_at_cub_points_side_refcell , 00477 COMP_BLAS ); 00478 00479 RealSpaceTools<double>::subtract(rhs_and_soln_vec, bc_fields_per_side ); 00480 00481 00482 } 00483 00484 // solve linear system 00485 solver.POTRF('L',numFields,&fe_matrix[0],numFields,&info); 00486 solver.POTRS('L',numFields,1,&fe_matrix[0],numFields,&rhs_and_soln_vec[0],numFields,&info); 00487 00488 interp_points_ref.resize(1,numInterpPoints,cellDim); 00489 // get exact solution for comparison 00490 FieldContainer<double> exact_solution(1,numInterpPoints,cellDim); 00491 exact_solution.initialize(); 00492 u_exact( exact_solution , interp_points_ref , comp , x_order, y_order, z_order); 00493 interp_points_ref.resize(numInterpPoints,cellDim); 00494 00495 // compute interpolant 00496 // first evaluate basis at interpolation points 00497 value_of_basis_at_interp_points.resize(1,numFields,numInterpPoints,cellDim); 00498 FunctionSpaceTools::evaluate<double>( interpolant , 00499 rhs_and_soln_vec , 00500 value_of_basis_at_interp_points ); 00501 value_of_basis_at_interp_points.resize(numFields,numInterpPoints,cellDim); 00502 00503 RealSpaceTools<double>::subtract(interpolant,exact_solution); 00504 00505 double nrm= RealSpaceTools<double>::vectorNorm(&interpolant[0],interpolant.dimension(1), NORM_TWO); 00506 00507 *outStream << "\nNorm-2 error between scalar components of exact solution of order (" 00508 << x_order << ", " << y_order << ", " << z_order 00509 << ") in component " << comp 00510 << " and finite element interpolant of order " << basis_order << ": " 00511 << nrm << "\n"; 00512 00513 if (nrm > zero) { 00514 *outStream << "\n\nPatch test failed for solution polynomial order (" 00515 << x_order << ", " << y_order << ", " << z_order << ") and basis order " 00516 << basis_order << "\n\n"; 00517 errorFlag++; 00518 } 00519 } 00520 } 00521 } 00522 } 00523 } 00524 00525 } 00526 00527 catch (std::logic_error err) { 00528 *outStream << err.what() << "\n\n"; 00529 errorFlag = -1000; 00530 }; 00531 00532 if (errorFlag != 0) 00533 std::cout << "End Result: TEST FAILED\n"; 00534 else 00535 std::cout << "End Result: TEST PASSED\n"; 00536 00537 // reset format state of std::cout 00538 std::cout.copyfmt(oldFormatState); 00539 00540 return errorFlag; 00541 }
1.7.4