|
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_HGRAD_TRI_Cn_FEM_ORTH.hpp" 00038 #include "Intrepid_HDIV_TRI_In_FEM.hpp" 00039 #include "Intrepid_DefaultCubatureFactory.hpp" 00040 #include "Intrepid_RealSpaceTools.hpp" 00041 #include "Intrepid_ArrayTools.hpp" 00042 #include "Intrepid_FunctionSpaceTools.hpp" 00043 #include "Intrepid_CellTools.hpp" 00044 #include "Teuchos_oblackholestream.hpp" 00045 #include "Teuchos_RCP.hpp" 00046 #include "Teuchos_GlobalMPISession.hpp" 00047 #include "Teuchos_SerialDenseMatrix.hpp" 00048 #include "Teuchos_SerialDenseVector.hpp" 00049 #include "Teuchos_LAPACK.hpp" 00050 00051 using namespace std; 00052 using namespace Intrepid; 00053 00054 void rhsFunc( FieldContainer<double> &, const FieldContainer<double> &, int, int ); 00055 void u_exact( FieldContainer<double> &, const FieldContainer<double> &, int, int ); 00056 00057 // This is the rhs for (div tau,w) = (f,w), 00058 // which makes f the negative Laplacian of scalar solution 00059 void rhsFunc( FieldContainer<double> &result, 00060 const FieldContainer<double> &points, 00061 int xd, 00062 int yd ) 00063 { 00064 for (int cell=0;cell<result.dimension(0);cell++) { 00065 for (int pt=0;pt<result.dimension(1);pt++) { 00066 result(cell,pt) = 0.0; 00067 if (xd >=2) { 00068 result(cell,pt) += xd*(xd-1)*pow(points(cell,pt,0),xd-2)*pow(points(cell,pt,1),yd); 00069 } 00070 if (yd >=2) { 00071 result(cell,pt) += yd*(yd-1)*pow(points(cell,pt,0),xd)*pow(points(cell,pt,1),yd-2); 00072 } 00073 } 00074 } 00075 } 00076 00077 void u_exact( FieldContainer<double> &result, 00078 const FieldContainer<double> &points, 00079 int xd, 00080 int yd) 00081 { 00082 for (int cell=0;cell<result.dimension(0);cell++){ 00083 for (int pt=0;pt<result.dimension(1);pt++) { 00084 result(cell,pt) = std::pow(points(cell,pt,0),xd)*std::pow(points(cell,pt,1),yd); 00085 } 00086 } 00087 return; 00088 } 00089 00090 int main(int argc, char *argv[]) { 00091 Teuchos::GlobalMPISession mpiSession(&argc, &argv); 00092 00093 // This little trick lets us print to std::cout only if 00094 // a (dummy) command-line argument is provided. 00095 int iprint = argc - 1; 00096 Teuchos::RCP<std::ostream> outStream; 00097 Teuchos::oblackholestream bhs; // outputs nothing 00098 if (iprint > 0) 00099 outStream = Teuchos::rcp(&std::cout, false); 00100 else 00101 outStream = Teuchos::rcp(&bhs, false); 00102 00103 // Save the format state of the original std::cout. 00104 Teuchos::oblackholestream oldFormatState; 00105 oldFormatState.copyfmt(std::cout); 00106 00107 *outStream \ 00108 << "===============================================================================\n" \ 00109 << "| |\n" \ 00110 << "| Unit Test (Basis_HGRAD_TRI_Cn_FEM_ORTH) |\n" \ 00111 << "| |\n" \ 00112 << "| 1) Patch test involving H(div) matrices |\n" \ 00113 << "| for the Dirichlet problem on a triangular patch |\n" \ 00114 << "| Omega with boundary Gamma. |\n" \ 00115 << "| |\n" \ 00116 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \ 00117 << "| Robert Kirby (robert.c.kirby@ttu.edu), |\n" \ 00118 << "| Denis Ridzal (dridzal@sandia.gov), |\n" \ 00119 << "| Kara Peterson (kjpeter@sandia.gov). |\n" \ 00120 << "| |\n" \ 00121 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \ 00122 << "| Trilinos website: http://trilinos.sandia.gov |\n" \ 00123 << "| |\n" \ 00124 << "===============================================================================\n"\ 00125 << "| TEST 1: Patch test |\n"\ 00126 << "===============================================================================\n"; 00127 00128 00129 int errorFlag = 0; 00130 00131 outStream -> precision(16); 00132 00133 try { 00134 DefaultCubatureFactory<double> cubFactory; // create cubature factory 00135 shards::CellTopology cell(shards::getCellTopologyData< shards::Triangle<> >()); // create parent cell topology 00136 shards::CellTopology side(shards::getCellTopologyData< shards::Line<> >()); // create relevant subcell (side) topology 00137 00138 int cellDim = cell.getDimension(); 00139 int sideDim = side.getDimension(); 00140 00141 int min_order = 0; 00142 int max_order = 9; 00143 00144 int numIntervals = max_order; 00145 int numInterpPoints = ((numIntervals + 1)*(numIntervals + 2))/2; 00146 FieldContainer<double> interp_points_ref(numInterpPoints, 2); 00147 int counter = 0; 00148 for (int j=0; j<=numIntervals; j++) { 00149 for (int i=0; i<=numIntervals; i++) { 00150 if (i <= numIntervals-j) { 00151 interp_points_ref(counter,0) = i*(1.0/numIntervals); 00152 interp_points_ref(counter,1) = j*(1.0/numIntervals); 00153 counter++; 00154 } 00155 } 00156 } 00157 00158 interp_points_ref.resize(numInterpPoints,2); 00159 00160 for (int basis_order=min_order;basis_order<=max_order;basis_order++) { 00161 // create bases 00162 Teuchos::RCP<Basis<double,FieldContainer<double> > > vectorBasis = 00163 Teuchos::rcp(new Basis_HDIV_TRI_In_FEM<double,FieldContainer<double> >(basis_order+1,POINTTYPE_EQUISPACED) ); 00164 Teuchos::RCP<Basis<double,FieldContainer<double> > > scalarBasis = 00165 Teuchos::rcp(new Basis_HGRAD_TRI_Cn_FEM_ORTH<double,FieldContainer<double> >(basis_order) ); 00166 00167 int numVectorFields = vectorBasis->getCardinality(); 00168 int numScalarFields = scalarBasis->getCardinality(); 00169 int numTotalFields = numVectorFields + numScalarFields; 00170 00171 // create cubatures 00172 Teuchos::RCP<Cubature<double> > cellCub = cubFactory.create(cell, 2*(basis_order+1)); 00173 Teuchos::RCP<Cubature<double> > sideCub = cubFactory.create(side, 2*(basis_order+1)); 00174 00175 int numCubPointsCell = cellCub->getNumPoints(); 00176 int numCubPointsSide = sideCub->getNumPoints(); 00177 00178 // hold cubature information 00179 FieldContainer<double> cub_points_cell(numCubPointsCell, cellDim); 00180 FieldContainer<double> cub_weights_cell(numCubPointsCell); 00181 FieldContainer<double> cub_points_side( numCubPointsSide, sideDim ); 00182 FieldContainer<double> cub_weights_side( numCubPointsSide ); 00183 FieldContainer<double> cub_points_side_refcell( numCubPointsSide , cellDim ); 00184 00185 // hold basis function information on refcell 00186 FieldContainer<double> value_of_v_basis_at_cub_points_cell(numVectorFields, numCubPointsCell, cellDim ); 00187 FieldContainer<double> w_value_of_v_basis_at_cub_points_cell(1, numVectorFields, numCubPointsCell, cellDim); 00188 FieldContainer<double> div_of_v_basis_at_cub_points_cell( numVectorFields, numCubPointsCell ); 00189 FieldContainer<double> w_div_of_v_basis_at_cub_points_cell( 1, numVectorFields , numCubPointsCell ); 00190 FieldContainer<double> value_of_s_basis_at_cub_points_cell(numScalarFields,numCubPointsCell); 00191 FieldContainer<double> w_value_of_s_basis_at_cub_points_cell(1,numScalarFields,numCubPointsCell); 00192 00193 // containers for side integration: 00194 // I just need the normal component of the vector basis 00195 // and the exact solution at the cub points 00196 FieldContainer<double> value_of_v_basis_at_cub_points_side(numVectorFields,numCubPointsSide,cellDim); 00197 FieldContainer<double> n_of_v_basis_at_cub_points_side(numVectorFields,numCubPointsSide); 00198 FieldContainer<double> w_n_of_v_basis_at_cub_points_side(1,numVectorFields,numCubPointsSide); 00199 FieldContainer<double> diri_data_at_cub_points_side(1,numCubPointsSide); 00200 FieldContainer<double> side_normal(cellDim); 00201 00202 // holds rhs data 00203 FieldContainer<double> rhs_at_cub_points_cell(1,numCubPointsCell); 00204 00205 // FEM matrices and vectors 00206 FieldContainer<double> fe_matrix_M(1,numVectorFields,numVectorFields); 00207 FieldContainer<double> fe_matrix_B(1,numVectorFields,numScalarFields); 00208 FieldContainer<double> fe_matrix(1,numTotalFields,numTotalFields); 00209 00210 FieldContainer<double> rhs_vector_vec(1,numVectorFields); 00211 FieldContainer<double> rhs_vector_scal(1,numScalarFields); 00212 FieldContainer<double> rhs_and_soln_vec(1,numTotalFields); 00213 00214 FieldContainer<int> ipiv(numTotalFields); 00215 FieldContainer<double> value_of_s_basis_at_interp_points( numScalarFields , numInterpPoints); 00216 FieldContainer<double> interpolant( 1 , numInterpPoints ); 00217 00218 // set test tolerance 00219 double zero = (basis_order+1)*(basis_order+1)*100*INTREPID_TOL; 00220 00221 // build matrices outside the loop, and then just do the rhs 00222 // for each iteration 00223 00224 cellCub->getCubature(cub_points_cell, cub_weights_cell); 00225 sideCub->getCubature(cub_points_side, cub_weights_side); 00226 00227 // need the vector basis & its divergences 00228 vectorBasis->getValues(value_of_v_basis_at_cub_points_cell, 00229 cub_points_cell, 00230 OPERATOR_VALUE); 00231 vectorBasis->getValues(div_of_v_basis_at_cub_points_cell, 00232 cub_points_cell, 00233 OPERATOR_DIV); 00234 00235 // need the scalar basis as well 00236 scalarBasis->getValues(value_of_s_basis_at_cub_points_cell, 00237 cub_points_cell, 00238 OPERATOR_VALUE); 00239 00240 // construct mass matrix 00241 cub_weights_cell.resize(1,numCubPointsCell); 00242 FunctionSpaceTools::multiplyMeasure<double>(w_value_of_v_basis_at_cub_points_cell , 00243 cub_weights_cell , 00244 value_of_v_basis_at_cub_points_cell ); 00245 cub_weights_cell.resize(numCubPointsCell); 00246 00247 00248 value_of_v_basis_at_cub_points_cell.resize( 1 , numVectorFields , numCubPointsCell , cellDim ); 00249 FunctionSpaceTools::integrate<double>(fe_matrix_M, 00250 w_value_of_v_basis_at_cub_points_cell , 00251 value_of_v_basis_at_cub_points_cell , 00252 COMP_BLAS ); 00253 value_of_v_basis_at_cub_points_cell.resize( numVectorFields , numCubPointsCell , cellDim ); 00254 00255 // div matrix 00256 cub_weights_cell.resize(1,numCubPointsCell); 00257 FunctionSpaceTools::multiplyMeasure<double>(w_div_of_v_basis_at_cub_points_cell, 00258 cub_weights_cell, 00259 div_of_v_basis_at_cub_points_cell); 00260 cub_weights_cell.resize(numCubPointsCell); 00261 00262 value_of_s_basis_at_cub_points_cell.resize(1,numScalarFields,numCubPointsCell); 00263 FunctionSpaceTools::integrate<double>(fe_matrix_B, 00264 w_div_of_v_basis_at_cub_points_cell , 00265 value_of_s_basis_at_cub_points_cell , 00266 COMP_BLAS ); 00267 value_of_s_basis_at_cub_points_cell.resize(numScalarFields,numCubPointsCell); 00268 00269 00270 // construct div matrix 00271 00272 for (int x_order=0;x_order<=basis_order;x_order++) { 00273 for (int y_order=0;y_order<=basis_order-x_order;y_order++) { 00274 00275 // reset global matrix since I destroyed it in LU factorization. 00276 fe_matrix.initialize(); 00277 // insert mass matrix into global matrix 00278 for (int i=0;i<numVectorFields;i++) { 00279 for (int j=0;j<numVectorFields;j++) { 00280 fe_matrix(0,i,j) = fe_matrix_M(0,i,j); 00281 } 00282 } 00283 00284 // insert div matrix into global matrix 00285 for (int i=0;i<numVectorFields;i++) { 00286 for (int j=0;j<numScalarFields;j++) { 00287 fe_matrix(0,i,numVectorFields+j)=-fe_matrix_B(0,i,j); 00288 fe_matrix(0,j+numVectorFields,i)=fe_matrix_B(0,i,j); 00289 } 00290 } 00291 00292 // clear old vector data 00293 rhs_vector_vec.initialize(); 00294 rhs_vector_scal.initialize(); 00295 rhs_and_soln_vec.initialize(); 00296 00297 // now get rhs vector 00298 // rhs_vector_scal is just (rhs,w) for w in the scalar basis 00299 // I already have the scalar basis tabulated. 00300 cub_points_cell.resize(1,numCubPointsCell,cellDim); 00301 rhsFunc(rhs_at_cub_points_cell, 00302 cub_points_cell, 00303 x_order, 00304 y_order); 00305 00306 cub_points_cell.resize(numCubPointsCell,cellDim); 00307 00308 cub_weights_cell.resize(1,numCubPointsCell); 00309 FunctionSpaceTools::multiplyMeasure<double>(w_value_of_s_basis_at_cub_points_cell, 00310 cub_weights_cell, 00311 value_of_s_basis_at_cub_points_cell); 00312 cub_weights_cell.resize(numCubPointsCell); 00313 FunctionSpaceTools::integrate<double>(rhs_vector_scal, 00314 rhs_at_cub_points_cell, 00315 w_value_of_s_basis_at_cub_points_cell, 00316 COMP_BLAS); 00317 00318 for (int i=0;i<numScalarFields;i++) { 00319 rhs_and_soln_vec(0,numVectorFields+i) = rhs_vector_scal(0,i); 00320 } 00321 00322 00323 // now get <u,v.n> on boundary 00324 for (unsigned side_cur=0;side_cur<3;side_cur++) { 00325 // map side cubature to current side 00326 CellTools<double>::mapToReferenceSubcell( cub_points_side_refcell , 00327 cub_points_side , 00328 sideDim , 00329 (int)side_cur , 00330 cell ); 00331 00332 // Evaluate dirichlet data 00333 cub_points_side_refcell.resize(1,numCubPointsSide,cellDim); 00334 u_exact(diri_data_at_cub_points_side, 00335 cub_points_side_refcell,x_order,y_order); 00336 cub_points_side_refcell.resize(numCubPointsSide,cellDim); 00337 00338 // get normal direction, this has the edge weight factored into it already 00339 CellTools<double>::getReferenceSideNormal(side_normal , 00340 (int)side_cur,cell ); 00341 00342 // v.n at cub points on side 00343 vectorBasis->getValues(value_of_v_basis_at_cub_points_side , 00344 cub_points_side_refcell , 00345 OPERATOR_VALUE ); 00346 00347 00348 for (int i=0;i<numVectorFields;i++) { 00349 for (int j=0;j<numCubPointsSide;j++) { 00350 n_of_v_basis_at_cub_points_side(i,j) = 0.0; 00351 for (int k=0;k<cellDim;k++) { 00352 n_of_v_basis_at_cub_points_side(i,j) += side_normal(k) * 00353 value_of_v_basis_at_cub_points_side(i,j,k); 00354 } 00355 } 00356 } 00357 00358 cub_weights_side.resize(1,numCubPointsSide); 00359 FunctionSpaceTools::multiplyMeasure<double>(w_n_of_v_basis_at_cub_points_side, 00360 cub_weights_side, 00361 n_of_v_basis_at_cub_points_side); 00362 cub_weights_side.resize(numCubPointsSide); 00363 00364 FunctionSpaceTools::integrate<double>(rhs_vector_vec, 00365 diri_data_at_cub_points_side, 00366 w_n_of_v_basis_at_cub_points_side, 00367 COMP_BLAS, 00368 false); 00369 for (int i=0;i<numVectorFields;i++) { 00370 rhs_and_soln_vec(0,i) -= rhs_vector_vec(0,i); 00371 } 00372 00373 } 00374 00375 // solve linear system 00376 int info = 0; 00377 Teuchos::LAPACK<int, double> solver; 00378 solver.GESV(numTotalFields, 1, &fe_matrix[0], numTotalFields, &ipiv(0), &rhs_and_soln_vec[0], 00379 numTotalFields, &info); 00380 00381 // compute interpolant; the scalar entries are last 00382 scalarBasis->getValues(value_of_s_basis_at_interp_points, 00383 interp_points_ref, 00384 OPERATOR_VALUE); 00385 for (int pt=0;pt<numInterpPoints;pt++) { 00386 interpolant(0,pt)=0.0; 00387 for (int i=0;i<numScalarFields;i++) { 00388 interpolant(0,pt) += rhs_and_soln_vec(0,numVectorFields+i) 00389 * value_of_s_basis_at_interp_points(i,pt); 00390 } 00391 } 00392 00393 interp_points_ref.resize(1,numInterpPoints,cellDim); 00394 // get exact solution for comparison 00395 FieldContainer<double> exact_solution(1,numInterpPoints); 00396 u_exact( exact_solution , interp_points_ref , x_order, y_order); 00397 interp_points_ref.resize(numInterpPoints,cellDim); 00398 00399 RealSpaceTools<double>::add(interpolant,exact_solution); 00400 00401 double nrm= RealSpaceTools<double>::vectorNorm(&interpolant[0],interpolant.dimension(1), NORM_TWO); 00402 00403 *outStream << "\nNorm-2 error between scalar components of exact solution polynomial of order (" 00404 << x_order << ", " << y_order << ") and finite element interpolant of order " << basis_order << ": " 00405 << nrm << "\n"; 00406 00407 if (nrm > zero) { 00408 *outStream << "\n\nPatch test failed for solution polynomial order (" 00409 << x_order << ", " << y_order << ") and basis order (scalar / vector) (" 00410 << basis_order << ", " << basis_order + 1 << ")\n\n"; 00411 errorFlag++; 00412 } 00413 00414 } 00415 } 00416 } 00417 00418 } 00419 catch (std::logic_error err) { 00420 *outStream << err.what() << "\n\n"; 00421 errorFlag = -1000; 00422 }; 00423 00424 if (errorFlag != 0) 00425 std::cout << "End Result: TEST FAILED\n"; 00426 else 00427 std::cout << "End Result: TEST PASSED\n"; 00428 00429 // reset format state of std::cout 00430 std::cout.copyfmt(oldFormatState); 00431 00432 return errorFlag; 00433 }
1.7.4