|
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 u_exact( FieldContainer<double> &, const FieldContainer<double> &, int , int, int, int ); 00055 00056 void u_exact( FieldContainer<double> &result, 00057 const FieldContainer<double> &points, 00058 int comp, 00059 int xd, 00060 int yd, 00061 int zd) 00062 { 00063 for (int cell=0;cell<result.dimension(0);cell++){ 00064 for (int pt=0;pt<result.dimension(1);pt++) { 00065 result(cell,pt,comp) = std::pow(points(cell,pt,0),xd)*std::pow(points(cell,pt,1),yd) 00066 *std::pow(points(cell,pt,2),zd); 00067 } 00068 } 00069 return; 00070 } 00071 00072 void rhsFunc( FieldContainer<double> & result , 00073 const FieldContainer<double> &points , 00074 int comp, 00075 int xd, 00076 int yd, 00077 int zd ) 00078 { 00079 u_exact( result , points , comp , xd , yd , zd ); 00080 } 00081 00082 int main(int argc, char *argv[]) { 00083 Teuchos::GlobalMPISession mpiSession(&argc, &argv); 00084 00085 // This little trick lets us print to std::cout only if 00086 // a (dummy) command-line argument is provided. 00087 int iprint = argc - 1; 00088 Teuchos::RCP<std::ostream> outStream; 00089 Teuchos::oblackholestream bhs; // outputs nothing 00090 if (iprint > 0) 00091 outStream = Teuchos::rcp(&std::cout, false); 00092 else 00093 outStream = Teuchos::rcp(&bhs, false); 00094 00095 // Save the format state of the original std::cout. 00096 Teuchos::oblackholestream oldFormatState; 00097 oldFormatState.copyfmt(std::cout); 00098 00099 *outStream \ 00100 << "===============================================================================\n" \ 00101 << "| |\n" \ 00102 << "| Unit Test (Basis_HCURL_TET_In_FEM) |\n" \ 00103 << "| |\n" \ 00104 << "| 1) Patch test involving H(curl) matrices |\n" \ 00105 << "| |\n" \ 00106 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \ 00107 << "| Robert Kirby (robert.c.kirby@ttu.edu), |\n" \ 00108 << "| Denis Ridzal (dridzal@sandia.gov), |\n" \ 00109 << "| Kara Peterson (kjpeter@sandia.gov). |\n" \ 00110 << "| |\n" \ 00111 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \ 00112 << "| Trilinos website: http://trilinos.sandia.gov |\n" \ 00113 << "| |\n" \ 00114 << "===============================================================================\n" \ 00115 << "| TEST 2: Patch test for mass matrices |\n" \ 00116 << "===============================================================================\n"; 00117 00118 00119 int errorFlag = 0; 00120 00121 outStream -> precision(16); 00122 00123 try { 00124 DefaultCubatureFactory<double> cubFactory; // create cubature factory 00125 shards::CellTopology cell(shards::getCellTopologyData< shards::Tetrahedron<> >()); // create parent cell topology 00126 00127 int cellDim = cell.getDimension(); 00128 00129 int min_order = 1; 00130 int max_order = 5; 00131 00132 int numIntervals = max_order; 00133 int numInterpPoints = ((numIntervals + 1)*(numIntervals + 2)*(numIntervals+3))/6; 00134 FieldContainer<double> interp_points_ref(numInterpPoints, cellDim); 00135 int counter = 0; 00136 for (int j=0; j<=numIntervals; j++) { 00137 for (int i=0; i<=numIntervals-j; i++) { 00138 for (int k=0;k<numIntervals-j-i;k++) { 00139 interp_points_ref(counter,0) = i*(1.0/numIntervals); 00140 interp_points_ref(counter,1) = j*(1.0/numIntervals); 00141 interp_points_ref(counter,2) = k*(1.0/numIntervals); 00142 counter++; 00143 } 00144 } 00145 } 00146 00147 for (int basis_order=min_order;basis_order<=max_order;basis_order++) { 00148 // create basis 00149 Teuchos::RCP<Basis<double,FieldContainer<double> > > basis = 00150 Teuchos::rcp(new Basis_HCURL_TET_In_FEM<double,FieldContainer<double> >(basis_order,POINTTYPE_EQUISPACED) ); 00151 00152 int numFields = basis->getCardinality(); 00153 00154 // create cubatures 00155 Teuchos::RCP<Cubature<double> > cellCub = cubFactory.create(cell, 2*(basis_order+1)); 00156 00157 int numCubPointsCell = cellCub->getNumPoints(); 00158 00159 // hold cubature information 00160 FieldContainer<double> cub_points_cell(numCubPointsCell, cellDim); 00161 FieldContainer<double> cub_weights_cell(numCubPointsCell); 00162 00163 // hold basis function information on refcell 00164 FieldContainer<double> value_of_basis_at_cub_points_cell(numFields, numCubPointsCell, cellDim ); 00165 FieldContainer<double> w_value_of_basis_at_cub_points_cell(1, numFields, numCubPointsCell, cellDim); 00166 00167 // holds rhs data 00168 FieldContainer<double> rhs_at_cub_points_cell(1,numCubPointsCell,cellDim); 00169 00170 // FEM mass matrix 00171 FieldContainer<double> fe_matrix_bak(1,numFields,numFields); 00172 FieldContainer<double> fe_matrix(1,numFields,numFields); 00173 FieldContainer<double> rhs_and_soln_vec(1,numFields); 00174 00175 FieldContainer<int> ipiv(numFields); 00176 FieldContainer<double> value_of_basis_at_interp_points( numFields , numInterpPoints , cellDim); 00177 FieldContainer<double> interpolant( 1, numInterpPoints , cellDim ); 00178 00179 int info = 0; 00180 Teuchos::LAPACK<int, double> solver; 00181 00182 // set test tolerance 00183 double zero = (basis_order+1)*(basis_order+1)*1000*INTREPID_TOL; 00184 00185 // build matrices outside the loop, and then just do the rhs 00186 // for each iteration 00187 cellCub->getCubature(cub_points_cell, cub_weights_cell); 00188 00189 // need the vector basis 00190 basis->getValues(value_of_basis_at_cub_points_cell, 00191 cub_points_cell, 00192 OPERATOR_VALUE); 00193 basis->getValues( value_of_basis_at_interp_points , 00194 interp_points_ref , 00195 OPERATOR_VALUE ); 00196 00197 00198 00199 00200 // construct mass matrix 00201 cub_weights_cell.resize(1,numCubPointsCell); 00202 FunctionSpaceTools::multiplyMeasure<double>(w_value_of_basis_at_cub_points_cell , 00203 cub_weights_cell , 00204 value_of_basis_at_cub_points_cell ); 00205 cub_weights_cell.resize(numCubPointsCell); 00206 00207 00208 value_of_basis_at_cub_points_cell.resize( 1 , numFields , numCubPointsCell , cellDim ); 00209 FunctionSpaceTools::integrate<double>(fe_matrix_bak, 00210 w_value_of_basis_at_cub_points_cell , 00211 value_of_basis_at_cub_points_cell , 00212 COMP_BLAS ); 00213 value_of_basis_at_cub_points_cell.resize( numFields , numCubPointsCell , cellDim ); 00214 00215 00216 //std::cout << fe_matrix_bak << std::endl; 00217 00218 for (int x_order=0;x_order<basis_order;x_order++) { 00219 for (int y_order=0;y_order<basis_order-x_order;y_order++) { 00220 for (int z_order=0;z_order<basis_order-x_order-y_order;z_order++) { 00221 for (int comp=0;comp<cellDim;comp++) { 00222 fe_matrix.initialize(); 00223 // copy mass matrix 00224 for (int i=0;i<numFields;i++) { 00225 for (int j=0;j<numFields;j++) { 00226 fe_matrix(0,i,j) = fe_matrix_bak(0,i,j); 00227 } 00228 } 00229 00230 // clear old vector data 00231 rhs_and_soln_vec.initialize(); 00232 00233 // now get rhs vector 00234 00235 cub_points_cell.resize(1,numCubPointsCell,cellDim); 00236 00237 rhs_at_cub_points_cell.initialize(); 00238 rhsFunc(rhs_at_cub_points_cell, 00239 cub_points_cell, 00240 comp, 00241 x_order, 00242 y_order, 00243 z_order); 00244 00245 cub_points_cell.resize(numCubPointsCell,cellDim); 00246 00247 cub_weights_cell.resize(numCubPointsCell); 00248 00249 FunctionSpaceTools::integrate<double>(rhs_and_soln_vec, 00250 rhs_at_cub_points_cell, 00251 w_value_of_basis_at_cub_points_cell, 00252 COMP_BLAS); 00253 00254 // solve linear system 00255 00256 // solver.GESV(numFields, 1, &fe_matrix[0], numFields, &ipiv(0), &rhs_and_soln_vec[0], 00257 // numFields, &info); 00258 solver.POTRF('L',numFields,&fe_matrix[0],numFields,&info); 00259 solver.POTRS('L',numFields,1,&fe_matrix[0],numFields,&rhs_and_soln_vec[0],numFields,&info); 00260 00261 interp_points_ref.resize(1,numInterpPoints,cellDim); 00262 // get exact solution for comparison 00263 FieldContainer<double> exact_solution(1,numInterpPoints,cellDim); 00264 exact_solution.initialize(); 00265 u_exact( exact_solution , interp_points_ref , comp , x_order, y_order, z_order); 00266 interp_points_ref.resize(numInterpPoints,cellDim); 00267 00268 // compute interpolant 00269 // first evaluate basis at interpolation points 00270 value_of_basis_at_interp_points.resize(1,numFields,numInterpPoints,cellDim); 00271 FunctionSpaceTools::evaluate<double>( interpolant , 00272 rhs_and_soln_vec , 00273 value_of_basis_at_interp_points ); 00274 value_of_basis_at_interp_points.resize(numFields,numInterpPoints,cellDim); 00275 00276 RealSpaceTools<double>::subtract(interpolant,exact_solution); 00277 00278 double nrm= RealSpaceTools<double>::vectorNorm(&interpolant[0],interpolant.dimension(1), NORM_TWO); 00279 00280 *outStream << "\nNorm-2 error between scalar components of exact solution of order (" 00281 << x_order << ", " << y_order << ", " << z_order 00282 << ") in component " << comp 00283 << " and finite element interpolant of order " << basis_order << ": " 00284 << nrm << "\n"; 00285 00286 if (nrm > zero) { 00287 *outStream << "\n\nPatch test failed for solution polynomial order (" 00288 << x_order << ", " << y_order << ", " << z_order << ") and basis order (scalar, vector) (" 00289 << basis_order << ", " << basis_order+1 << ")\n\n"; 00290 errorFlag++; 00291 } 00292 } 00293 } 00294 } 00295 } 00296 } 00297 00298 } 00299 00300 catch (std::logic_error err) { 00301 *outStream << err.what() << "\n\n"; 00302 errorFlag = -1000; 00303 }; 00304 00305 if (errorFlag != 0) 00306 std::cout << "End Result: TEST FAILED\n"; 00307 else 00308 std::cout << "End Result: TEST PASSED\n"; 00309 00310 // reset format state of std::cout 00311 std::cout.copyfmt(oldFormatState); 00312 00313 return errorFlag; 00314 }
1.7.4