|
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 00082 // Intrepid includes 00083 #include "Intrepid_FunctionSpaceTools.hpp" 00084 #include "Intrepid_FieldContainer.hpp" 00085 #include "Intrepid_CellTools.hpp" 00086 //#include "Intrepid_ArrayTools.hpp" 00087 #include "Intrepid_HGRAD_HEX_Cn_FEM.hpp" 00088 //#include "Intrepid_RealSpaceTools.hpp" 00089 #include "Intrepid_DefaultCubatureFactory.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_BLAS.hpp" 00102 //#include "Teuchos_BLAS_types.hpp" 00103 00104 // Shards includes 00105 #include "Shards_CellTopology.hpp" 00106 00107 // EpetraExt includes 00108 #include "EpetraExt_MultiVectorOut.h" 00109 00110 using namespace std; 00111 using namespace Intrepid; 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_10.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: Build 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 // Element to Node map 00240 // We'll keep it around, but this is only the DOFMap if you are in the lowest order case. 00241 FieldContainer<int> elemToNode(numElems, numNodesPerElem); 00242 int ielem = 0; 00243 for (int k=0; k<NZ; k++) 00244 { 00245 for (int j=0; j<NY; j++) 00246 { 00247 for (int i=0; i<NX; i++) 00248 { 00249 elemToNode(ielem,0) = k * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i; 00250 elemToNode(ielem,1) = k * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i + 1; 00251 elemToNode(ielem,2) = k * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i + 1; 00252 elemToNode(ielem,3) = k * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i; 00253 elemToNode(ielem,4) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i; 00254 elemToNode(ielem,5) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + j * ( NX + 1 ) + i + 1; 00255 elemToNode(ielem,6) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i + 1; 00256 elemToNode(ielem,7) = ( k + 1 ) * ( NX + 1 ) * ( NY + 1 ) + ( j + 1 ) * ( NX + 1 ) + i; 00257 ielem++; 00258 } 00259 } 00260 } 00261 #ifdef DUMP_DATA 00262 // Output connectivity 00263 ofstream fe2nout("elem2node.dat"); 00264 for (int k=0;k<NZ;k++) 00265 { 00266 for (int j=0; j<NY; j++) 00267 { 00268 for (int i=0; i<NX; i++) 00269 { 00270 int ielem = i + j * NX + k * NY * NY; 00271 for (int m=0; m<numNodesPerElem; m++) 00272 { 00273 fe2nout << elemToNode(ielem,m) <<" "; 00274 } 00275 fe2nout <<"\n"; 00276 } 00277 } 00278 } 00279 fe2nout.close(); 00280 #endif 00281 00282 // ************************************ CUBATURE ************************************** 00283 *outStream << "Getting cubature ... \n\n"; 00284 00285 // Get numerical integration points and weights 00286 DefaultCubatureFactory<double> cubFactory; 00287 int cubDegree = 2*deg; 00288 Teuchos::RCP<Cubature<double> > quadCub = cubFactory.create(hex_8, cubDegree); 00289 00290 int cubDim = quadCub->getDimension(); 00291 int numCubPoints = quadCub->getNumPoints(); 00292 00293 FieldContainer<double> cubPoints(numCubPoints, cubDim); 00294 FieldContainer<double> cubWeights(numCubPoints); 00295 00296 quadCub->getCubature(cubPoints, cubWeights); 00297 00298 // ************************************** BASIS *************************************** 00299 00300 *outStream << "Getting basis ... \n\n"; 00301 00302 // Define basis 00303 Basis_HGRAD_HEX_Cn_FEM<double, FieldContainer<double> > quadHGradBasis(deg,POINTTYPE_SPECTRAL); 00304 int numFieldsG = quadHGradBasis.getCardinality(); 00305 FieldContainer<double> quadGVals(numFieldsG, numCubPoints); 00306 FieldContainer<double> quadGrads(numFieldsG, numCubPoints, spaceDim); 00307 00308 // Evaluate basis values and gradients at cubature points 00309 quadHGradBasis.getValues(quadGVals, cubPoints, OPERATOR_VALUE); 00310 quadHGradBasis.getValues(quadGrads, cubPoints, OPERATOR_GRAD); 00311 00312 // create the local-global mapping 00313 FieldContainer<int> ltgMapping(numElems,numFieldsG); 00314 const int numDOF = (NX*deg+1)*(NY*deg+1)*(NZ*deg+1); 00315 ielem=0; 00316 for (int k=0;k<NZ;k++) 00317 { 00318 for (int j=0;j<NY;j++) 00319 { 00320 for (int i=0;i<NX;i++) 00321 { 00322 const int start = k * ( NY * deg + 1 ) * ( NX * deg + 1 ) + j * ( NX * deg + 1 ) + i * deg; 00323 // loop over local dof on this cell 00324 int local_dof_cur=0; 00325 for (int kloc=0;kloc<=deg;kloc++) 00326 { 00327 for (int jloc=0;jloc<=deg;jloc++) 00328 { 00329 for (int iloc=0;iloc<=deg;iloc++) 00330 { 00331 ltgMapping(ielem,local_dof_cur) = start 00332 + kloc * ( NX * deg + 1 ) * ( NY * deg + 1 ) 00333 + jloc * ( NX * deg + 1 ) 00334 + iloc; 00335 local_dof_cur++; 00336 } 00337 } 00338 } 00339 ielem++; 00340 } 00341 } 00342 } 00343 #ifdef DUMP_DATA 00344 // Output ltg mapping 00345 ielem = 0; 00346 ofstream ltgout("ltg.dat"); 00347 for (int k=0;k<NZ;k++) 00348 { 00349 for (int j=0; j<NY; j++) 00350 { 00351 for (int i=0; i<NX; i++) 00352 { 00353 int ielem = i + j * NX + k * NX * NY; 00354 for (int m=0; m<numFieldsG; m++) 00355 { 00356 ltgout << ltgMapping(ielem,m) <<" "; 00357 } 00358 ltgout <<"\n"; 00359 } 00360 } 00361 } 00362 ltgout.close(); 00363 #endif 00364 00365 // ********** DECLARE GLOBAL OBJECTS ************* 00366 Epetra_SerialComm Comm; 00367 Epetra_Map globalMapG(numDOF, 0, Comm); 00368 Epetra_FEVector u(globalMapG); u.Random(); 00369 Epetra_FEVector Ku(globalMapG); 00370 00371 00372 00373 // ********** CONSTRUCT AND INSERT LOCAL STIFFNESS MATRICES *********** 00374 *outStream << "Building reference stiffness matrix...\n\n"; 00375 typedef CellTools<double> CellTools; 00376 typedef FunctionSpaceTools fst; 00377 00378 // jacobian information 00379 FieldContainer<double> refCellNodes(1,numNodesPerElem,spaceDim); 00380 FieldContainer<double> cellJacobian(1,numCubPoints,spaceDim,spaceDim); 00381 FieldContainer<double> cellJacobInv(1,numCubPoints,spaceDim,spaceDim); 00382 FieldContainer<double> cellJacobDet(1,numCubPoints); 00383 00384 // element stiffness matrices and supporting storage space 00385 FieldContainer<double> localStiffMatrix(1, numFieldsG, numFieldsG); 00386 FieldContainer<double> transformedBasisGradients(1,numFieldsG,numCubPoints,spaceDim); 00387 FieldContainer<double> weightedTransformedBasisGradients(1,numFieldsG,numCubPoints,spaceDim); 00388 FieldContainer<double> weightedMeasure(1, numCubPoints); 00389 00390 Epetra_Time localConstructTimer( Comm ); 00391 refCellNodes(0,0,0) = 0.0; refCellNodes(0,0,1) = 0.0; refCellNodes(0,0,2) = 0.0; 00392 refCellNodes(0,1,0) = hx; refCellNodes(0,1,1) = 0.0; refCellNodes(0,1,2) = 0.0; 00393 refCellNodes(0,2,0) = hx; refCellNodes(0,2,1) = hy; refCellNodes(0,2,2) = 0.0; 00394 refCellNodes(0,3,0) = 0.0; refCellNodes(0,3,1) = hy; refCellNodes(0,3,2) = 0.0; 00395 refCellNodes(0,4,0) = 0.0; refCellNodes(0,4,1) = 0.0; refCellNodes(0,4,2) = hz; 00396 refCellNodes(0,5,0) = hx; refCellNodes(0,5,1) = 0.0; refCellNodes(0,5,2) = hz; 00397 refCellNodes(0,6,0) = hx; refCellNodes(0,6,1) = hy; refCellNodes(0,6,2) = hz; 00398 refCellNodes(0,7,0) = 0.0; refCellNodes(0,7,1) = hy; refCellNodes(0,7,2) = hz; 00399 00400 // jacobian evaluation 00401 CellTools::setJacobian(cellJacobian,cubPoints,refCellNodes,hex_8); 00402 CellTools::setJacobianInv(cellJacobInv, cellJacobian ); 00403 CellTools::setJacobianDet(cellJacobDet, cellJacobian ); 00404 00405 // transform reference element gradients to each cell 00406 fst::HGRADtransformGRAD<double>(transformedBasisGradients, cellJacobInv, quadGrads); 00407 00408 // compute weighted measure 00409 fst::computeCellMeasure<double>(weightedMeasure, cellJacobDet, cubWeights); 00410 00411 // multiply values with weighted measure 00412 fst::multiplyMeasure<double>(weightedTransformedBasisGradients, 00413 weightedMeasure, transformedBasisGradients); 00414 00415 // integrate to compute element stiffness matrix 00416 fst::integrate<double>(localStiffMatrix, 00417 transformedBasisGradients, weightedTransformedBasisGradients , COMP_BLAS); 00418 00419 const double localConstructTime = localConstructTimer.ElapsedTime(); 00420 00421 // ************* MATRIX-FREE APPLICATION 00422 FieldContainer<double> uScattered(numElems,numFieldsG); 00423 FieldContainer<double> KuScattered(numElems,numFieldsG); 00424 00425 u.GlobalAssemble(); 00426 00427 Epetra_Time multTimer(Comm); 00428 Teuchos::BLAS<int,double> blas; 00429 Ku.PutScalar(0.0); 00430 Ku.GlobalAssemble(); 00431 00432 double *uVals = u[0]; 00433 double *KuVals = Ku[0]; 00434 00435 Epetra_Time scatterTimer(Comm); 00436 std::cout << "Scattering\n"; 00437 // Scatter 00438 for (int k=0; k<numElems; k++) 00439 { 00440 for (int i=0;i<numFieldsG;i++) 00441 { 00442 uScattered(k,i) = uVals[ltgMapping(k,i)]; 00443 } 00444 } 00445 const double scatterTime = scatterTimer.ElapsedTime(); 00446 00447 Epetra_Time blasTimer(Comm); 00448 blas.GEMM(Teuchos::NO_TRANS , Teuchos::NO_TRANS , 00449 numFieldsG , numElems, numFieldsG , 00450 1.0 , 00451 &localStiffMatrix(0,0,0) , 00452 numFieldsG , 00453 &uScattered(0,0) , 00454 numFieldsG , 00455 0.0 , 00456 &KuScattered(0,0) , 00457 numFieldsG ); 00458 const double blasTime = blasTimer.ElapsedTime(); 00459 00460 Epetra_Time gatherTimer(Comm); 00461 // Gather 00462 for (int k=0;k<numElems;k++) 00463 { 00464 for (int i=0;i<numFieldsG;i++) 00465 { 00466 KuVals[ltgMapping(k,i)] += KuScattered(k,i); 00467 } 00468 } 00469 00470 const double gatherTime = gatherTimer.ElapsedTime(); 00471 00472 00473 *outStream << "Time to build local matrix (including Jacobian computation): "<< localConstructTime << "\n"; 00474 *outStream << "Time to scatter " << scatterTime << "\n"; 00475 *outStream << "Time for local application " << blasTime << "\n"; 00476 *outStream << "Time to gather " << gatherTime << "\n"; 00477 *outStream << "Total matrix-free time " << scatterTime + blasTime + gatherTime << "\n"; 00478 00479 00480 *outStream << "End Result: TEST PASSED\n"; 00481 00482 // reset format state of std::cout 00483 std::cout.copyfmt(oldFormatState); 00484 00485 return 0; 00486 } 00487
1.7.6.1