|
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 00069 // Intrepid includes 00070 #include "Intrepid_FunctionSpaceTools.hpp" 00071 #include "Intrepid_FieldContainer.hpp" 00072 #include "Intrepid_CellTools.hpp" 00073 //#include "Intrepid_ArrayTools.hpp" 00074 #include "Intrepid_HGRAD_HEX_Cn_FEM.hpp" 00075 //#include "Intrepid_RealSpaceTools.hpp" 00076 #include "Intrepid_DefaultCubatureFactory.hpp" 00077 #include "Intrepid_Utils.hpp" 00078 00079 // Epetra includes 00080 #include "Epetra_Time.h" 00081 #include "Epetra_Map.h" 00082 #include "Epetra_FEVector.h" 00083 #include "Epetra_SerialComm.h" 00084 00085 // Teuchos includes 00086 #include "Teuchos_oblackholestream.hpp" 00087 #include "Teuchos_RCP.hpp" 00088 //#include "Teuchos_BLAS.hpp" 00089 //#include "Teuchos_BLAS_types.hpp" 00090 00091 // Shards includes 00092 #include "Shards_CellTopology.hpp" 00093 00094 // EpetraExt includes 00095 #include "EpetraExt_MultiVectorOut.h" 00096 00097 using namespace std; 00098 using namespace Intrepid; 00099 00100 int main(int argc, char *argv[]) { 00101 00102 //Check number of arguments 00103 if (argc < 4) { 00104 std::cout <<"\n>>> ERROR: Invalid number of arguments.\n\n"; 00105 std::cout <<"Usage:\n\n"; 00106 std::cout <<" ./Intrepid_example_Drivers_Example_10.exe deg NX NY NZ verbose\n\n"; 00107 std::cout <<" where \n"; 00108 std::cout <<" int deg - polynomial degree to be used (assumed >= 1) \n"; 00109 std::cout <<" int NX - num intervals in x direction (assumed box domain, 0,1) \n"; 00110 std::cout <<" int NY - num intervals in y direction (assumed box domain, 0,1) \n"; 00111 std::cout <<" int NZ - num intervals in y direction (assumed box domain, 0,1) \n"; 00112 std::cout <<" verbose (optional) - any character, indicates verbose output \n\n"; 00113 exit(1); 00114 } 00115 00116 // This little trick lets us print to std::cout only if 00117 // a (dummy) command-line argument is provided. 00118 int iprint = argc - 1; 00119 Teuchos::RCP<std::ostream> outStream; 00120 Teuchos::oblackholestream bhs; // outputs nothing 00121 if (iprint > 2) 00122 outStream = Teuchos::rcp(&std::cout, false); 00123 else 00124 outStream = Teuchos::rcp(&bhs, false); 00125 00126 // Save the format state of the original std::cout. 00127 Teuchos::oblackholestream oldFormatState; 00128 oldFormatState.copyfmt(std::cout); 00129 00130 *outStream \ 00131 << "===============================================================================\n" \ 00132 << "| |\n" \ 00133 << "| Example: Build Stiffness Matrix for |\n" \ 00134 << "| Poisson Equation on Hexahedral Mesh |\n" \ 00135 << "| |\n" \ 00136 << "| Questions? Contact Pavel Bochev (pbboche@sandia.gov), |\n" \ 00137 << "| Denis Ridzal (dridzal@sandia.gov), |\n" \ 00138 << "| Kara Peterson (kjpeter@sandia.gov). |\n" \ 00139 << "| |\n" \ 00140 << "| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \ 00141 << "| Trilinos website: http://trilinos.sandia.gov |\n" \ 00142 << "| |\n" \ 00143 << "===============================================================================\n"; 00144 00145 00146 // ************************************ GET INPUTS ************************************** 00147 00148 int deg = atoi(argv[1]); // polynomial degree to use 00149 int NX = atoi(argv[2]); // num intervals in x direction (assumed box domain, 0,1) 00150 int NY = atoi(argv[3]); // num intervals in y direction (assumed box domain, 0,1) 00151 int NZ = atoi(argv[4]); // num intervals in y direction (assumed box domain, 0,1) 00152 00153 00154 // *********************************** CELL TOPOLOGY ********************************** 00155 00156 // Get cell topology for base hexahedron 00157 typedef shards::CellTopology CellTopology; 00158 CellTopology hex_8(shards::getCellTopologyData<shards::Hexahedron<8> >() ); 00159 00160 // Get dimensions 00161 int numNodesPerElem = hex_8.getNodeCount(); 00162 int spaceDim = hex_8.getDimension(); 00163 00164 // *********************************** GENERATE MESH ************************************ 00165 00166 *outStream << "Generating mesh ... \n\n"; 00167 00168 *outStream << " NX" << " NY" << " NZ\n"; 00169 *outStream << std::setw(5) << NX << 00170 std::setw(5) << NY << std::setw(5) << NZ << "\n\n"; 00171 00172 // Print mesh information 00173 int numElems = NX*NY*NZ; 00174 int numNodes = (NX+1)*(NY+1)*(NZ+1); 00175 *outStream << " Number of Elements: " << numElems << " \n"; 00176 *outStream << " Number of Nodes: " << numNodes << " \n\n"; 00177 00178 // Cube 00179 double leftX = 0.0, rightX = 1.0; 00180 double leftY = 0.0, rightY = 1.0; 00181 double leftZ = 0.0, rightZ = 1.0; 00182 00183 // Mesh spacing 00184 double hx = (rightX-leftX)/((double)NX); 00185 double hy = (rightY-leftY)/((double)NY); 00186 double hz = (rightZ-leftZ)/((double)NZ); 00187 00188 // Get nodal coordinates 00189 FieldContainer<double> nodeCoord(numNodes, spaceDim); 00190 FieldContainer<int> nodeOnBoundary(numNodes); 00191 int inode = 0; 00192 for (int k=0; k<NZ+1; k++) 00193 { 00194 for (int j=0; j<NY+1; j++) 00195 { 00196 for (int i=0; i<NX+1; i++) 00197 { 00198 nodeCoord(inode,0) = leftX + (double)i*hx; 00199 nodeCoord(inode,1) = leftY + (double)j*hy; 00200 nodeCoord(inode,2) = leftZ + (double)k*hz; 00201 if (k==0 || k==NZ || j==0 || i==0 || j==NY || i==NX) 00202 { 00203 nodeOnBoundary(inode)=1; 00204 } 00205 else 00206 { 00207 nodeOnBoundary(inode)=0; 00208 } 00209 inode++; 00210 } 00211 } 00212 } 00213 #define DUMP_DATA 00214 #ifdef DUMP_DATA 00215 // Print nodal coords 00216 ofstream fcoordout("coords.dat"); 00217 for (int i=0; i<numNodes; i++) { 00218 fcoordout << nodeCoord(i,0) <<" "; 00219 fcoordout << nodeCoord(i,1) <<" "; 00220 fcoordout << nodeCoord(i,2) <<"\n"; 00221 } 00222 fcoordout.close(); 00223 #endif 00224 00225 00226 // Element to Node map 00227 // We'll keep it around, but this is only the DOFMap if you are in the lowest order case. 00228 FieldContainer<int> elemToNode(numElems, numNodesPerElem); 00229 int ielem = 0; 00230 for (int k=0; k<NZ; k++) 00231 { 00232 for (int j=0; j<NY; j++) 00233 { 00234 for (int i=0; i<NX; i++) 00235 { 00236 elemToNode(ielem,0) = k * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i; 00237 elemToNode(ielem,1) = k * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i + 1; 00238 elemToNode(ielem,2) = k * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i + 1; 00239 elemToNode(ielem,3) = k * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i; 00240 elemToNode(ielem,4) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i; 00241 elemToNode(ielem,5) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i + 1; 00242 elemToNode(ielem,6) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i + 1; 00243 elemToNode(ielem,7) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i; 00244 ielem++; 00245 } 00246 } 00247 } 00248 #ifdef DUMP_DATA 00249 // Output connectivity 00250 ofstream fe2nout("elem2node.dat"); 00251 for (int k=0;k<NZ;k++) 00252 { 00253 for (int j=0; j<NY; j++) 00254 { 00255 for (int i=0; i<NX; i++) 00256 { 00257 int ielem = i + j * NX + k * NY * NY; 00258 for (int m=0; m<numNodesPerElem; m++) 00259 { 00260 fe2nout << elemToNode(ielem,m) <<" "; 00261 } 00262 fe2nout <<"\n"; 00263 } 00264 } 00265 } 00266 fe2nout.close(); 00267 #endif 00268 00269 // ************************************ CUBATURE ************************************** 00270 *outStream << "Getting cubature ... \n\n"; 00271 00272 // Get numerical integration points and weights 00273 DefaultCubatureFactory<double> cubFactory; 00274 int cubDegree = 2*deg; 00275 Teuchos::RCP<Cubature<double> > quadCub = cubFactory.create(hex_8, cubDegree); 00276 00277 int cubDim = quadCub->getDimension(); 00278 int numCubPoints = quadCub->getNumPoints(); 00279 00280 FieldContainer<double> cubPoints(numCubPoints, cubDim); 00281 FieldContainer<double> cubWeights(numCubPoints); 00282 00283 quadCub->getCubature(cubPoints, cubWeights); 00284 00285 // ************************************** BASIS *************************************** 00286 00287 *outStream << "Getting basis ... \n\n"; 00288 00289 // Define basis 00290 Basis_HGRAD_HEX_Cn_FEM<double, FieldContainer<double> > quadHGradBasis(deg,POINTTYPE_SPECTRAL); 00291 int numFieldsG = quadHGradBasis.getCardinality(); 00292 FieldContainer<double> quadGVals(numFieldsG, numCubPoints); 00293 FieldContainer<double> quadGrads(numFieldsG, numCubPoints, spaceDim); 00294 00295 // Evaluate basis values and gradients at cubature points 00296 quadHGradBasis.getValues(quadGVals, cubPoints, OPERATOR_VALUE); 00297 quadHGradBasis.getValues(quadGrads, cubPoints, OPERATOR_GRAD); 00298 00299 // create the local-global mapping 00300 FieldContainer<int> ltgMapping(numElems,numFieldsG); 00301 const int numDOF = (NX*deg+1)*(NY*deg+1)*(NZ*deg+1); 00302 ielem=0; 00303 for (int k=0;k<NZ;k++) 00304 { 00305 for (int j=0;j<NY;j++) 00306 { 00307 for (int i=0;i<NX;i++) 00308 { 00309 const int start = k * ( NY * deg + 1 ) * ( NX * deg + 1 ) + j * ( NX * deg + 1 ) + i * deg; 00310 // loop over local dof on this cell 00311 int local_dof_cur=0; 00312 for (int kloc=0;kloc<=deg;kloc++) 00313 { 00314 for (int jloc=0;jloc<=deg;jloc++) 00315 { 00316 for (int iloc=0;iloc<=deg;iloc++) 00317 { 00318 ltgMapping(ielem,local_dof_cur) = start 00319 + kloc * ( NX * deg + 1 ) * ( NY * deg + 1 ) 00320 + jloc * ( NX * deg + 1 ) 00321 + iloc; 00322 local_dof_cur++; 00323 } 00324 } 00325 } 00326 ielem++; 00327 } 00328 } 00329 } 00330 #ifdef DUMP_DATA 00331 // Output ltg mapping 00332 ielem = 0; 00333 ofstream ltgout("ltg.dat"); 00334 for (int k=0;k<NZ;k++) 00335 { 00336 for (int j=0; j<NY; j++) 00337 { 00338 for (int i=0; i<NX; i++) 00339 { 00340 int ielem = i + j * NX + k * NX * NY; 00341 for (int m=0; m<numFieldsG; m++) 00342 { 00343 ltgout << ltgMapping(ielem,m) <<" "; 00344 } 00345 ltgout <<"\n"; 00346 } 00347 } 00348 } 00349 ltgout.close(); 00350 #endif 00351 00352 // ********** DECLARE GLOBAL OBJECTS ************* 00353 Epetra_SerialComm Comm; 00354 Epetra_Map globalMapG(numDOF, 0, Comm); 00355 Epetra_FEVector u(globalMapG); u.Random(); 00356 Epetra_FEVector Ku(globalMapG); 00357 00358 00359 00360 // ********** CONSTRUCT AND INSERT LOCAL STIFFNESS MATRICES *********** 00361 *outStream << "Building reference stiffness matrix...\n\n"; 00362 typedef CellTools<double> CellTools; 00363 typedef FunctionSpaceTools fst; 00364 00365 // jacobian information 00366 FieldContainer<double> refCellNodes(1,numNodesPerElem,spaceDim); 00367 FieldContainer<double> cellJacobian(1,numCubPoints,spaceDim,spaceDim); 00368 FieldContainer<double> cellJacobInv(1,numCubPoints,spaceDim,spaceDim); 00369 FieldContainer<double> cellJacobDet(1,numCubPoints); 00370 00371 // element stiffness matrices and supporting storage space 00372 FieldContainer<double> localStiffMatrix(1, numFieldsG, numFieldsG); 00373 FieldContainer<double> transformedBasisGradients(1,numFieldsG,numCubPoints,spaceDim); 00374 FieldContainer<double> weightedTransformedBasisGradients(1,numFieldsG,numCubPoints,spaceDim); 00375 FieldContainer<double> weightedMeasure(1, numCubPoints); 00376 00377 Epetra_Time localConstructTimer( Comm ); 00378 refCellNodes(0,0,0) = 0.0; refCellNodes(0,0,1) = 0.0; refCellNodes(0,0,2) = 0.0; 00379 refCellNodes(0,1,0) = hx; refCellNodes(0,1,1) = 0.0; refCellNodes(0,1,2) = 0.0; 00380 refCellNodes(0,2,0) = hx; refCellNodes(0,2,1) = hy; refCellNodes(0,2,2) = 0.0; 00381 refCellNodes(0,3,0) = 0.0; refCellNodes(0,3,1) = hy; refCellNodes(0,3,2) = 0.0; 00382 refCellNodes(0,4,0) = 0.0; refCellNodes(0,4,1) = 0.0; refCellNodes(0,4,2) = hz; 00383 refCellNodes(0,5,0) = hx; refCellNodes(0,5,1) = 0.0; refCellNodes(0,5,2) = hz; 00384 refCellNodes(0,6,0) = hx; refCellNodes(0,6,1) = hy; refCellNodes(0,6,2) = hz; 00385 refCellNodes(0,7,0) = 0.0; refCellNodes(0,7,1) = hy; refCellNodes(0,7,2) = hz; 00386 00387 // jacobian evaluation 00388 CellTools::setJacobian(cellJacobian,cubPoints,refCellNodes,hex_8); 00389 CellTools::setJacobianInv(cellJacobInv, cellJacobian ); 00390 CellTools::setJacobianDet(cellJacobDet, cellJacobian ); 00391 00392 // transform reference element gradients to each cell 00393 fst::HGRADtransformGRAD<double>(transformedBasisGradients, cellJacobInv, quadGrads); 00394 00395 // compute weighted measure 00396 fst::computeCellMeasure<double>(weightedMeasure, cellJacobDet, cubWeights); 00397 00398 // multiply values with weighted measure 00399 fst::multiplyMeasure<double>(weightedTransformedBasisGradients, 00400 weightedMeasure, transformedBasisGradients); 00401 00402 // integrate to compute element stiffness matrix 00403 fst::integrate<double>(localStiffMatrix, 00404 transformedBasisGradients, weightedTransformedBasisGradients , COMP_BLAS); 00405 00406 const double localConstructTime = localConstructTimer.ElapsedTime(); 00407 00408 // ************* MATRIX-FREE APPLICATION 00409 FieldContainer<double> uScattered(numElems,numFieldsG); 00410 FieldContainer<double> KuScattered(numElems,numFieldsG); 00411 00412 u.GlobalAssemble(); 00413 00414 Epetra_Time multTimer(Comm); 00415 Teuchos::BLAS<int,double> blas; 00416 Ku.PutScalar(0.0); 00417 Ku.GlobalAssemble(); 00418 00419 double *uVals = u[0]; 00420 double *KuVals = Ku[0]; 00421 00422 Epetra_Time scatterTimer(Comm); 00423 std::cout << "Scattering\n"; 00424 // Scatter 00425 for (int k=0; k<numElems; k++) 00426 { 00427 for (int i=0;i<numFieldsG;i++) 00428 { 00429 uScattered(k,i) = uVals[ltgMapping(k,i)]; 00430 } 00431 } 00432 const double scatterTime = scatterTimer.ElapsedTime(); 00433 00434 Epetra_Time blasTimer(Comm); 00435 blas.GEMM(Teuchos::NO_TRANS , Teuchos::NO_TRANS , 00436 numFieldsG , numElems, numFieldsG , 00437 1.0 , 00438 &localStiffMatrix(0,0,0) , 00439 numFieldsG , 00440 &uScattered(0,0) , 00441 numFieldsG , 00442 0.0 , 00443 &KuScattered(0,0) , 00444 numFieldsG ); 00445 const double blasTime = blasTimer.ElapsedTime(); 00446 00447 Epetra_Time gatherTimer(Comm); 00448 // Gather 00449 for (int k=0;k<numElems;k++) 00450 { 00451 for (int i=0;i<numFieldsG;i++) 00452 { 00453 KuVals[ltgMapping(k,i)] += KuScattered(k,i); 00454 } 00455 } 00456 00457 const double gatherTime = gatherTimer.ElapsedTime(); 00458 00459 00460 *outStream << "Time to build local matrix (including Jacobian computation): "<< localConstructTime << "\n"; 00461 *outStream << "Time to scatter " << scatterTime << "\n"; 00462 *outStream << "Time for local application " << blasTime << "\n"; 00463 *outStream << "Time to gather " << gatherTime << "\n"; 00464 *outStream << "Total matrix-free time " << scatterTime + blasTime + gatherTime << "\n"; 00465 00466 00467 *outStream << "End Result: TEST PASSED\n"; 00468 00469 // reset format state of std::cout 00470 std::cout.copyfmt(oldFormatState); 00471 00472 return 0; 00473 } 00474
1.7.4