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