|
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 00081 // Intrepid includes 00082 #include "Intrepid_FunctionSpaceTools.hpp" 00083 #include "Intrepid_FunctionSpaceToolsInPlace.hpp" 00084 #include "Intrepid_FieldContainer.hpp" 00085 #include "Intrepid_CellTools.hpp" 00086 #include "Intrepid_CubaturePolylib.hpp" 00087 #include "Intrepid_CubatureTensor.hpp" 00088 #include "Intrepid_HGRAD_HEX_Cn_FEM.hpp" 00089 #include "Intrepid_TensorProductSpaceTools.hpp" 00090 #include "Intrepid_Utils.hpp" 00091 00092 // Epetra includes 00093 #include "Epetra_Time.h" 00094 #include "Epetra_Map.h" 00095 #include "Epetra_FEVector.h" 00096 #include "Epetra_SerialComm.h" 00097 00098 // Teuchos includes 00099 #include "Teuchos_oblackholestream.hpp" 00100 #include "Teuchos_RCP.hpp" 00101 #include "Teuchos_Time.hpp" 00102 00103 // Shards includes 00104 #include "Shards_CellTopology.hpp" 00105 00106 // EpetraExt includes 00107 #include "EpetraExt_MultiVectorOut.h" 00108 00109 using namespace std; 00110 using namespace Intrepid; 00111 using Teuchos::rcp; 00112 00113 int main(int argc, char *argv[]) { 00114 00115 //Check number of arguments 00116 if (argc < 4) { 00117 std::cout <<"\n>>> ERROR: Invalid number of arguments.\n\n"; 00118 std::cout <<"Usage:\n\n"; 00119 std::cout <<" ./Intrepid_example_Drivers_Example_14.exe deg NX NY NZ verbose\n\n"; 00120 std::cout <<" where \n"; 00121 std::cout <<" int deg - polynomial degree to be used (assumed >= 1) \n"; 00122 std::cout <<" int NX - num intervals in x direction (assumed box domain, 0,1) \n"; 00123 std::cout <<" int NY - num intervals in y direction (assumed box domain, 0,1) \n"; 00124 std::cout <<" int NZ - num intervals in y direction (assumed box domain, 0,1) \n"; 00125 std::cout <<" verbose (optional) - any character, indicates verbose output \n\n"; 00126 exit(1); 00127 } 00128 00129 // This little trick lets us print to std::cout only if 00130 // a (dummy) command-line argument is provided. 00131 int iprint = argc - 1; 00132 Teuchos::RCP<std::ostream> outStream; 00133 Teuchos::oblackholestream bhs; // outputs nothing 00134 if (iprint > 2) 00135 outStream = Teuchos::rcp(&std::cout, false); 00136 else 00137 outStream = Teuchos::rcp(&bhs, false); 00138 00139 // Save the format state of the original std::cout. 00140 Teuchos::oblackholestream oldFormatState; 00141 oldFormatState.copyfmt(std::cout); 00142 00143 *outStream \ 00144 << "===============================================================================\n" \ 00145 << "| |\n" \ 00146 << "| Example: Apply Stiffness Matrix for |\n" \ 00147 << "| Poisson Equation on Hexahedral Mesh |\n" \ 00148 << "| |\n" \ 00149 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \ 00150 << "| Denis Ridzal (dridzal@sandia.gov), |\n" \ 00151 << "| Kara Peterson (kjpeter@sandia.gov). |\n" \ 00152 << "| |\n" \ 00153 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \ 00154 << "| Trilinos website: http://trilinos.sandia.gov |\n" \ 00155 << "| |\n" \ 00156 << "===============================================================================\n"; 00157 00158 00159 // ************************************ GET INPUTS ************************************** 00160 00161 int deg = atoi(argv[1]); // polynomial degree to use 00162 int NX = atoi(argv[2]); // num intervals in x direction (assumed box domain, 0,1) 00163 int NY = atoi(argv[3]); // num intervals in y direction (assumed box domain, 0,1) 00164 int NZ = atoi(argv[4]); // num intervals in y direction (assumed box domain, 0,1) 00165 00166 00167 // *********************************** CELL TOPOLOGY ********************************** 00168 00169 // Get cell topology for base hexahedron 00170 typedef shards::CellTopology CellTopology; 00171 CellTopology hex_8(shards::getCellTopologyData<shards::Hexahedron<8> >() ); 00172 00173 // Get dimensions 00174 int numNodesPerElem = hex_8.getNodeCount(); 00175 int spaceDim = hex_8.getDimension(); 00176 00177 // *********************************** GENERATE MESH ************************************ 00178 00179 *outStream << "Generating mesh ... \n\n"; 00180 00181 *outStream << " NX" << " NY" << " NZ\n"; 00182 *outStream << std::setw(5) << NX << 00183 std::setw(5) << NY << std::setw(5) << NZ << "\n\n"; 00184 00185 // Print mesh information 00186 int numElems = NX*NY*NZ; 00187 int numNodes = (NX+1)*(NY+1)*(NZ+1); 00188 *outStream << " Number of Elements: " << numElems << " \n"; 00189 *outStream << " Number of Nodes: " << numNodes << " \n\n"; 00190 00191 // Cube 00192 double leftX = 0.0, rightX = 1.0; 00193 double leftY = 0.0, rightY = 1.0; 00194 double leftZ = 0.0, rightZ = 1.0; 00195 00196 // Mesh spacing 00197 double hx = (rightX-leftX)/((double)NX); 00198 double hy = (rightY-leftY)/((double)NY); 00199 double hz = (rightZ-leftZ)/((double)NZ); 00200 00201 // Get nodal coordinates 00202 FieldContainer<double> nodeCoord(numNodes, spaceDim); 00203 FieldContainer<int> nodeOnBoundary(numNodes); 00204 int inode = 0; 00205 for (int k=0; k<NZ+1; k++) 00206 { 00207 for (int j=0; j<NY+1; j++) 00208 { 00209 for (int i=0; i<NX+1; i++) 00210 { 00211 nodeCoord(inode,0) = leftX + (double)i*hx; 00212 nodeCoord(inode,1) = leftY + (double)j*hy; 00213 nodeCoord(inode,2) = leftZ + (double)k*hz; 00214 if (k==0 || k==NZ || j==0 || i==0 || j==NY || i==NX) 00215 { 00216 nodeOnBoundary(inode)=1; 00217 } 00218 else 00219 { 00220 nodeOnBoundary(inode)=0; 00221 } 00222 inode++; 00223 } 00224 } 00225 } 00226 //#define DUMP_DATA 00227 #ifdef DUMP_DATA 00228 // Print nodal coords 00229 ofstream fcoordout("coords.dat"); 00230 for (int i=0; i<numNodes; i++) { 00231 fcoordout << nodeCoord(i,0) <<" "; 00232 fcoordout << nodeCoord(i,1) <<" "; 00233 fcoordout << nodeCoord(i,2) <<"\n"; 00234 } 00235 fcoordout.close(); 00236 #endif 00237 00238 00239 00240 // ********************************* CUBATURE AND BASIS *********************************** 00241 *outStream << "Getting cubature and basis ... \n\n"; 00242 00243 // Get numerical integration points and weights 00244 // I only need this on the line since I'm doing tensor products 00245 Teuchos::RCP<Cubature<double,FieldContainer<double>,FieldContainer<double> > > glcub 00246 = Teuchos::rcp(new CubaturePolylib<double,FieldContainer<double>,FieldContainer<double> >(2*deg-1,PL_GAUSS_LOBATTO) ); 00247 00248 const int numCubPoints1D = glcub->getNumPoints(); 00249 00250 FieldContainer<double> cubPoints1D(numCubPoints1D, 1); 00251 FieldContainer<double> cubWeights1D(numCubPoints1D); 00252 00253 glcub->getCubature(cubPoints1D,cubWeights1D); 00254 00255 std::vector<Teuchos::RCP<Cubature<double,FieldContainer<double>,FieldContainer<double> > > > 00256 cub_to_tensor; 00257 cub_to_tensor.push_back( glcub ); 00258 cub_to_tensor.push_back( glcub ); 00259 cub_to_tensor.push_back( glcub ); 00260 00261 Array<RCP<FieldContainer<double> > > wts_by_dim(3); 00262 wts_by_dim[0] = rcp( &cubWeights1D , false ); wts_by_dim[1] = wts_by_dim[0]; wts_by_dim[2] = wts_by_dim[1]; 00263 00264 CubatureTensor<double,FieldContainer<double>,FieldContainer<double> > cubhex( cub_to_tensor ); 00265 00266 Basis_HGRAD_HEX_Cn_FEM<double, FieldContainer<double> > hexBasis( deg , POINTTYPE_SPECTRAL ); 00267 00268 Array< Array< RCP< Basis< double , FieldContainer<double> > > > > &bases = hexBasis.getBases(); 00269 00270 // get the bases tabulated at the quadrature points, dimension-by-dimension 00271 00272 Array< RCP< FieldContainer<double> > > basisVals( 3 ); 00273 FieldContainer<double> bvals1D( bases[0][0]->getCardinality() , numCubPoints1D ); 00274 bases[0][0]->getValues( bvals1D , cubPoints1D , OPERATOR_VALUE ); 00275 basisVals[0] = rcp( &bvals1D , false ); basisVals[1] = basisVals[0]; basisVals[2] = basisVals[0]; 00276 00277 Array< RCP< FieldContainer<double> > > basisDVals( 3 ); 00278 FieldContainer<double> bdvals1D( bases[0][0]->getCardinality() , numCubPoints1D , 1); 00279 bases[0][0]->getValues( bdvals1D , cubPoints1D , OPERATOR_D1 ); 00280 basisDVals[0] = rcp( &bdvals1D , false ); basisDVals[1] = basisDVals[0]; basisDVals[2] = basisDVals[0]; 00281 00282 00283 const int numCubPoints = cubhex.getNumPoints(); 00284 FieldContainer<double> cubPoints3D( numCubPoints , 3 ); 00285 FieldContainer<double> cubWeights3D( numCubPoints ); 00286 cubhex.getCubature( cubPoints3D , cubWeights3D ); 00287 00288 00289 FieldContainer<int> elemToNode(numElems, numNodesPerElem); 00290 int ielem = 0; 00291 for (int k=0; k<NZ; k++) 00292 { 00293 for (int j=0; j<NY; j++) 00294 { 00295 for (int i=0; i<NX; i++) 00296 { 00297 elemToNode(ielem,0) = k * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i; 00298 elemToNode(ielem,1) = k * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i + 1; 00299 elemToNode(ielem,2) = k * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i + 1; 00300 elemToNode(ielem,3) = k * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i; 00301 elemToNode(ielem,4) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i; 00302 elemToNode(ielem,5) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i + 1; 00303 elemToNode(ielem,6) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i + 1; 00304 elemToNode(ielem,7) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i; 00305 ielem++; 00306 } 00307 } 00308 } 00309 #ifdef DUMP_DATA 00310 // Output connectivity 00311 ofstream fe2nout("elem2node.dat"); 00312 for (int k=0;k<NZ;k++) 00313 { 00314 for (int j=0; j<NY; j++) 00315 { 00316 for (int i=0; i<NX; i++) 00317 { 00318 int ielem = i + j * NX + k * NY * NY; 00319 for (int m=0; m<numNodesPerElem; m++) 00320 { 00321 fe2nout << elemToNode(ielem,m) <<" "; 00322 } 00323 fe2nout <<"\n"; 00324 } 00325 } 00326 } 00327 fe2nout.close(); 00328 #endif 00329 00330 00331 // ********************************* 3-D LOCAL-TO-GLOBAL MAPPING ******************************* 00332 FieldContainer<int> ltgMapping(numElems,hexBasis.getCardinality()); 00333 const int numDOF = (NX*deg+1)*(NY*deg+1)*(NZ*deg+1); 00334 ielem=0; 00335 for (int k=0;k<NZ;k++) 00336 { 00337 for (int j=0;j<NY;j++) 00338 { 00339 for (int i=0;i<NX;i++) 00340 { 00341 const int start = k * ( NY * deg + 1 ) * ( NX * deg + 1 ) + j * ( NX * deg + 1 ) + i * deg; 00342 // loop over local dof on this cell 00343 int local_dof_cur=0; 00344 for (int kloc=0;kloc<=deg;kloc++) 00345 { 00346 for (int jloc=0;jloc<=deg;jloc++) 00347 { 00348 for (int iloc=0;iloc<=deg;iloc++) 00349 { 00350 ltgMapping(ielem,local_dof_cur) = start 00351 + kloc * ( NX * deg + 1 ) * ( NY * deg + 1 ) 00352 + jloc * ( NX * deg + 1 ) 00353 + iloc; 00354 local_dof_cur++; 00355 } 00356 } 00357 } 00358 ielem++; 00359 } 00360 } 00361 } 00362 #ifdef DUMP_DATA 00363 // Output ltg mapping 00364 ielem = 0; 00365 ofstream ltgout("ltg.dat"); 00366 for (int k=0;k<NZ;k++) 00367 { 00368 for (int j=0; j<NY; j++) 00369 { 00370 for (int i=0; i<NX; i++) 00371 { 00372 int ielem = i + j * NX + k * NX * NY; 00373 for (int m=0; m<hexBasis.getCardinality(); m++) 00374 { 00375 ltgout << ltgMapping(ielem,m) <<" "; 00376 } 00377 ltgout <<"\n"; 00378 } 00379 } 00380 } 00381 ltgout.close(); 00382 #endif 00383 00384 // ********** DECLARE GLOBAL OBJECTS ************* 00385 Epetra_SerialComm Comm; 00386 Epetra_Map globalMapG(numDOF, 0, Comm); 00387 Epetra_FEVector u(globalMapG); u.Random(); 00388 Epetra_FEVector Ku(globalMapG); 00389 00390 // ************* For Jacobians ********************** 00391 FieldContainer<double> cellVertices(numElems,numNodesPerElem,spaceDim); 00392 FieldContainer<double> cellJacobian(numElems,numCubPoints,spaceDim,spaceDim); 00393 FieldContainer<double> cellJacobInv(numElems,numCubPoints,spaceDim,spaceDim); 00394 FieldContainer<double> cellJacobDet(numElems,numCubPoints); 00395 00396 00397 // get vertices of cells (for computing Jacobians) 00398 for (int i=0;i<numElems;i++) 00399 { 00400 for (int j=0;j<numNodesPerElem;j++) 00401 { 00402 const int nodeCur = elemToNode(i,j); 00403 for (int k=0;k<spaceDim;k++) 00404 { 00405 cellVertices(i,j,k) = nodeCoord(nodeCur,k); 00406 } 00407 } 00408 } 00409 00410 // jacobian evaluation 00411 CellTools<double>::setJacobian(cellJacobian,cubPoints3D,cellVertices,hex_8); 00412 CellTools<double>::setJacobianInv(cellJacobInv, cellJacobian ); 00413 CellTools<double>::setJacobianDet(cellJacobDet, cellJacobian ); 00414 00415 00416 // ************* MATRIX-FREE APPLICATION 00417 FieldContainer<double> uScattered(numElems,1,hexBasis.getCardinality()); 00418 FieldContainer<double> KuScattered(numElems,1,hexBasis.getCardinality()); 00419 FieldContainer<double> gradU(numElems,1,hexBasis.getCardinality(),3); 00420 00421 u.GlobalAssemble(); 00422 00423 00424 00425 Ku.PutScalar(0.0); 00426 Ku.GlobalAssemble(); 00427 00428 double *uVals = u[0]; 00429 double *KuVals = Ku[0]; 00430 00431 Teuchos::Time full_timer( "Time to apply operator matrix-free:" ); 00432 Teuchos::Time scatter_timer( "Time to scatter dof:" ); 00433 Teuchos::Time elementwise_timer( "Time to do elementwise computation:" ); 00434 Teuchos::Time grad_timer( "Time to compute gradients:" ); 00435 Teuchos::Time pointwise_timer( "Time to do pointwise transformations:" ); 00436 Teuchos::Time moment_timer( "Time to compute moments:" ); 00437 Teuchos::Time gather_timer( "Time to gather dof:" ); 00438 00439 full_timer.start(); 00440 00441 scatter_timer.start(); 00442 for (int k=0; k<numElems; k++) 00443 { 00444 for (int i=0;i<hexBasis.getCardinality();i++) 00445 { 00446 uScattered(k,0,i) = uVals[ltgMapping(k,i)]; 00447 } 00448 } 00449 scatter_timer.stop(); 00450 00451 elementwise_timer.start(); 00452 00453 grad_timer.start(); 00454 Intrepid::TensorProductSpaceTools::evaluateGradient<double>( gradU , uScattered ,basisVals , basisDVals ); 00455 grad_timer.stop(); 00456 pointwise_timer.start(); 00457 Intrepid::FunctionSpaceToolsInPlace::HGRADtransformGRAD<double>( gradU , cellJacobian ); 00458 Intrepid::FunctionSpaceToolsInPlace::HGRADtransformGRADDual<double>( gradU , cellJacobian ); 00459 Intrepid::FunctionSpaceToolsInPlace::multiplyMeasure<double>( gradU , cellJacobDet ); 00460 pointwise_timer.stop(); 00461 moment_timer.start(); 00462 Intrepid::TensorProductSpaceTools::momentsGrad<double>( KuScattered , gradU , basisVals , basisDVals , wts_by_dim ); 00463 moment_timer.stop(); 00464 elementwise_timer.stop(); 00465 gather_timer.start(); 00466 for (int k=0;k<numElems;k++) 00467 { 00468 for (int i=0;i<hexBasis.getCardinality();i++) 00469 { 00470 KuVals[ltgMapping(k,i)] += KuScattered(k,0,i); 00471 } 00472 } 00473 gather_timer.stop(); 00474 full_timer.stop(); 00475 00476 *outStream << full_timer.name() << " " << full_timer.totalElapsedTime() << " sec\n"; 00477 *outStream << "\t" << scatter_timer.name() << " " << scatter_timer.totalElapsedTime() << " sec\n"; 00478 *outStream << "\t" << elementwise_timer.name() << " " << elementwise_timer.totalElapsedTime() << " sec\n"; 00479 *outStream << "\t\t" << grad_timer.name() << " " << grad_timer.totalElapsedTime() << " sec\n"; 00480 *outStream << "\t\t" << pointwise_timer.name() << " " << pointwise_timer.totalElapsedTime() << " sec\n"; 00481 *outStream << "\t\t" << moment_timer.name() << " " << moment_timer.totalElapsedTime() << " sec\n"; 00482 *outStream << "\t" << gather_timer.name() << " " << gather_timer.totalElapsedTime() << " sec\n"; 00483 00484 00485 *outStream << "End Result: TEST PASSED\n"; 00486 00487 // reset format state of std::cout 00488 std::cout.copyfmt(oldFormatState); 00489 00490 return 0; 00491 } 00492
1.7.6.1