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