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