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