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