Intrepid
http://trilinos.sandia.gov/packages/docs/r10.12/packages/intrepid/example/Drivers/example_06.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 
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_QUAD_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_06.exe deg NX NY 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 <<"   verbose (optional)  - any character, indicates verbose output \n\n";
00125       exit(1);
00126    }
00127   
00128   // This little trick lets us print to std::cout only if
00129   // a (dummy) command-line argument is provided.
00130   int iprint     = argc - 1;
00131   Teuchos::RCP<std::ostream> outStream;
00132   Teuchos::oblackholestream bhs; // outputs nothing
00133   if (iprint > 2)
00134     outStream = Teuchos::rcp(&std::cout, false);
00135   else
00136     outStream = Teuchos::rcp(&bhs, false);
00137   
00138   // Save the format state of the original std::cout.
00139   Teuchos::oblackholestream oldFormatState;
00140   oldFormatState.copyfmt(std::cout);
00141   
00142   *outStream \
00143     << "===============================================================================\n" \
00144     << "|                                                                             |\n" \
00145     << "|  Example: Apply Stiffness Matrix for                                        |\n" \
00146     << "|                   Poisson Equation on Quadrilateral Mesh                    |\n" \
00147     << "|                                                                             |\n" \
00148     << "|  Questions? Contact  Pavel Bochev  (pbboche@sandia.gov),                    |\n" \
00149     << "|                      Denis Ridzal  (dridzal@sandia.gov),                    |\n" \
00150     << "|                      Kara Peterson (kjpeter@sandia.gov).                    |\n" \
00151     << "|                                                                             |\n" \
00152     << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00153     << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00154     << "|                                                                             |\n" \
00155     << "===============================================================================\n";
00156 
00157   
00158   // ************************************ GET INPUTS **************************************
00159   
00160   int deg          = atoi(argv[1]);  // polynomial degree to use
00161   int NX            = atoi(argv[2]);  // num intervals in x direction (assumed box domain, 0,1)
00162   int NY            = atoi(argv[3]);  // num intervals in y direction (assumed box domain, 0,1)
00163   
00164 
00165   // *********************************** CELL TOPOLOGY **********************************
00166   
00167   // Get cell topology for base hexahedron
00168   typedef shards::CellTopology    CellTopology;
00169   CellTopology quad_4(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
00170   
00171   // Get dimensions 
00172   int numNodesPerElem = quad_4.getNodeCount();
00173   int spaceDim = quad_4.getDimension();
00174   
00175   // *********************************** GENERATE MESH ************************************
00176   
00177   *outStream << "Generating mesh ... \n\n";
00178   
00179   *outStream << "   NX" << "   NY\n";
00180   *outStream << std::setw(5) << NX <<
00181     std::setw(5) << NY << "\n\n";
00182   
00183   // Print mesh information
00184   int numElems = NX*NY;
00185   int numNodes = (NX+1)*(NY+1);
00186   *outStream << " Number of Elements: " << numElems << " \n";
00187   *outStream << "    Number of Nodes: " << numNodes << " \n\n";
00188   
00189   // Square
00190   double leftX = 0.0, rightX = 1.0;
00191   double leftY = 0.0, rightY = 1.0;
00192 
00193   // Mesh spacing
00194   double hx = (rightX-leftX)/((double)NX);
00195   double hy = (rightY-leftY)/((double)NY);
00196 
00197   // Get nodal coordinates
00198   FieldContainer<double> nodeCoord(numNodes, spaceDim);
00199   FieldContainer<int> nodeOnBoundary(numNodes);
00200   int inode = 0;
00201   for (int j=0; j<NY+1; j++) {
00202     for (int i=0; i<NX+1; i++) {
00203       nodeCoord(inode,0) = leftX + (double)i*hx;
00204       nodeCoord(inode,1) = leftY + (double)j*hy;
00205       if (j==0 || i==0 || j==NY || i==NX){
00206         nodeOnBoundary(inode)=1;
00207       }
00208       else {
00209         nodeOnBoundary(inode)=0;
00210       }
00211       inode++;
00212     }
00213   }
00214 #define DUMP_DATA
00215 #ifdef DUMP_DATA
00216   // Print nodal coords
00217   ofstream fcoordout("coords.dat");
00218   for (int i=0; i<numNodes; i++) {
00219     fcoordout << nodeCoord(i,0) <<" ";
00220     fcoordout << nodeCoord(i,1) <<"\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 j=0; j<NY; j++) {
00231     for (int i=0; i<NX; i++) {
00232       elemToNode(ielem,0) = (NX + 1)*j + i;
00233       elemToNode(ielem,1) = (NX + 1)*j + i + 1;
00234       elemToNode(ielem,2) = (NX + 1)*(j + 1) + i + 1;
00235       elemToNode(ielem,3) = (NX + 1)*(j + 1) + i;
00236       ielem++;
00237     }
00238   }
00239 #ifdef DUMP_DATA
00240   // Output connectivity
00241   ofstream fe2nout("elem2node.dat");
00242   for (int j=0; j<NY; j++) {
00243     for (int i=0; i<NX; i++) {
00244       int ielem = i + j * NX;
00245       for (int m=0; m<numNodesPerElem; m++){
00246         fe2nout << elemToNode(ielem,m) <<"  ";
00247       }
00248       fe2nout <<"\n";
00249     }
00250   }
00251   fe2nout.close();
00252 #endif
00253   
00254   // ************************************ CUBATURE ************************************** 
00255   *outStream << "Getting cubature ... \n\n";
00256   
00257   // Get numerical integration points and weights
00258   DefaultCubatureFactory<double>  cubFactory;                                   
00259   int cubDegree = 2*deg;
00260   Teuchos::RCP<Cubature<double> > quadCub = cubFactory.create(quad_4, cubDegree); 
00261   
00262   int cubDim       = quadCub->getDimension();
00263   int numCubPoints = quadCub->getNumPoints();
00264   
00265   FieldContainer<double> cubPoints(numCubPoints, cubDim);
00266   FieldContainer<double> cubWeights(numCubPoints);
00267   
00268   quadCub->getCubature(cubPoints, cubWeights);
00269   
00270 
00271   // ************************************** BASIS ***************************************
00272   
00273   *outStream << "Getting basis ... \n\n";
00274   
00275   // Define basis 
00276   Basis_HGRAD_QUAD_Cn_FEM<double, FieldContainer<double> > quadHGradBasis(deg,POINTTYPE_SPECTRAL);
00277   int numFieldsG = quadHGradBasis.getCardinality();
00278   FieldContainer<double> quadGVals(numFieldsG, numCubPoints); 
00279   FieldContainer<double> quadGrads(numFieldsG, numCubPoints, spaceDim); 
00280   
00281   // Evaluate basis values and gradients at cubature points
00282   quadHGradBasis.getValues(quadGVals, cubPoints, OPERATOR_VALUE);
00283   quadHGradBasis.getValues(quadGrads, cubPoints, OPERATOR_GRAD);
00284 
00285   // create the local-global mapping for higher order elements
00286   FieldContainer<int> ltgMapping(numElems,numFieldsG);
00287   const int numDOF = (NX*deg+1)*(NY*deg+1);
00288   ielem=0;
00289   for (int j=0;j<NY;j++) {
00290     for (int i=0;i<NX;i++) {
00291       const int start = deg * j * ( NX * deg + 1 ) + i * deg;
00292       // loop over local dof on this cell
00293       int local_dof_cur=0;
00294       for (int vertical=0;vertical<=deg;vertical++) {
00295         for (int horizontal=0;horizontal<=deg;horizontal++) {
00296           ltgMapping(ielem,local_dof_cur) = start + vertical*(NX*deg+1)+horizontal;
00297           local_dof_cur++;
00298         }
00299       }
00300       ielem++;
00301     }
00302   }
00303 #ifdef DUMP_DATA
00304   // Output ltg mapping
00305 //   ofstream ltgout("ltg.dat");
00306 //   for (int j=0; j<NY; j++) {
00307 //     for (int i=0; i<NX; i++) {
00308 //       int ielem = i + j * NX;
00309 //       for (int m=0; m<numFieldsG; m++){
00310 //      ltgout << ltgMapping(ielem,m) <<"  ";
00311 //       }
00312 //       ltgout <<"\n";
00313 //     }
00314 //   }
00315 //   ltgout.close();
00316 #endif
00317   
00318   // ******** CREATE A SINGLE STIFFNESS MATRIX, WHICH IS REPLICATED ON ALL ELEMENTS *********
00319   *outStream << "Applying stiffness matrix and right hand side ... \n\n";
00320 
00321   // Settings and data structures for mass and stiffness matrices
00322   typedef CellTools<double>  CellTools;
00323   typedef FunctionSpaceTools fst;
00324   int numCells = 1; 
00325 
00326   // Container for nodes
00327   FieldContainer<double> refQuadNodes(numCells, numNodesPerElem, spaceDim);
00328   // Containers for Jacobian
00329   FieldContainer<double> refQuadJacobian(numCells, numCubPoints, spaceDim, spaceDim);
00330   FieldContainer<double> refQuadJacobInv(numCells, numCubPoints, spaceDim, spaceDim);
00331   FieldContainer<double> refQuadJacobDet(numCells, numCubPoints);
00332   // Containers for element HGRAD stiffness matrix
00333   FieldContainer<double> localStiffMatrix(numCells, numFieldsG, numFieldsG);
00334   FieldContainer<double> weightedMeasure(numCells, numCubPoints);
00335   FieldContainer<double> quadGradsTransformed(numCells, numFieldsG, numCubPoints, spaceDim);
00336   FieldContainer<double> quadGradsTransformedWeighted(numCells, numFieldsG, numCubPoints, spaceDim);
00337   // Containers for right hand side vectors
00338   FieldContainer<double> rhsData(numCells, numCubPoints);
00339   FieldContainer<double> localRHS(numCells, numFieldsG);
00340   FieldContainer<double> quadGValsTransformed(numCells, numFieldsG, numCubPoints);
00341   FieldContainer<double> quadGValsTransformedWeighted(numCells, numFieldsG, numCubPoints);
00342   // Container for cubature points in physical space
00343   FieldContainer<double> physCubPoints(numCells, numCubPoints, cubDim);
00344   
00345   // Global arrays in Epetra format 
00346   Epetra_SerialComm Comm;
00347   Epetra_Map globalMapG(numDOF, 0, Comm);
00348   Epetra_FEVector u(globalMapG);
00349   Epetra_FEVector Ku(globalMapG);
00350   u.Random();
00351 
00352   std::cout << "About to start ref element matrix\n";
00353 
00354   // ************************** Compute element HGrad stiffness matrices *******************************  
00355   refQuadNodes(0,0,0) = 0.0;
00356   refQuadNodes(0,0,1) = 0.0;
00357   refQuadNodes(0,1,0) = hx;
00358   refQuadNodes(0,1,1) = 0.0;
00359   refQuadNodes(0,2,0) = hx;
00360   refQuadNodes(0,2,1) = hy;
00361   refQuadNodes(0,3,0) = 0.0;
00362   refQuadNodes(0,3,1) = hy;
00363 
00364   // Compute cell Jacobians, their inverses and their determinants
00365   CellTools::setJacobian(refQuadJacobian, cubPoints, refQuadNodes, quad_4);
00366   CellTools::setJacobianInv(refQuadJacobInv, refQuadJacobian );
00367   CellTools::setJacobianDet(refQuadJacobDet, refQuadJacobian );
00368   
00369   // transform from [-1,1]^2 to [0,hx]x[0,hy]
00370   fst::HGRADtransformGRAD<double>(quadGradsTransformed, refQuadJacobInv, quadGrads);
00371       
00372   // compute weighted measure
00373   fst::computeCellMeasure<double>(weightedMeasure, refQuadJacobDet, cubWeights);
00374 
00375   // multiply values with weighted measure
00376   fst::multiplyMeasure<double>(quadGradsTransformedWeighted,
00377                                weightedMeasure, quadGradsTransformed);
00378 
00379   // integrate to compute element stiffness matrix
00380   fst::integrate<double>(localStiffMatrix,
00381                          quadGradsTransformed, quadGradsTransformedWeighted, COMP_BLAS);
00382 
00383   std::cout << "Finished with reference element matrix\n";
00384 
00385   
00386   // now we will scatter global degrees of freedom, apply the local stiffness matrix 
00387   // with BLAS, and then gather the results
00388   FieldContainer<double> uScattered(numElems,numFieldsG);
00389   FieldContainer<double> KuScattered(numElems,numFieldsG);
00390 
00391   // to extract info from u
00392 
00393   u.GlobalAssemble();
00394 
00395   Epetra_Time multTimer(Comm);
00396 
00397   Ku.PutScalar(0.0);
00398   Ku.GlobalAssemble();
00399 
00400   double *uVals = u[0];
00401   double *KuVals = Ku[0];
00402 
00403   Teuchos::BLAS<int,double> blas;
00404   Epetra_Time scatterTime(Comm);
00405   std::cout << "Scattering\n";
00406   // Scatter
00407   for (int k=0; k<numElems; k++) 
00408     {
00409       for (int i=0;i<numFieldsG;i++) 
00410         {
00411           uScattered(k,i) = uVals[ltgMapping(k,i)];
00412         }
00413     }
00414   const double scatTime = scatterTime.ElapsedTime();
00415   std::cout << "Scattered in time " << scatTime << "\n";
00416 
00417   Epetra_Time blasTimer(Comm);
00418   blas.GEMM(Teuchos::NO_TRANS , Teuchos::NO_TRANS , 
00419             numFieldsG , numElems, numFieldsG  , 
00420             1.0 , 
00421             &localStiffMatrix(0,0,0) , 
00422             numFieldsG ,
00423             &uScattered(0,0) , 
00424             numFieldsG , 
00425             0.0 , 
00426              &KuScattered(0,0) , 
00427             numFieldsG );
00428   const double blasTime = blasTimer.ElapsedTime();
00429   std::cout << "Element matrices applied in " << blasTime << "\n";
00430 
00431   Epetra_Time gatherTimer(Comm);
00432   // Gather
00433   for (int k=0;k<numElems;k++)
00434     {
00435       for (int i=0;i<numFieldsG;i++)
00436         {
00437           KuVals[ltgMapping(k,i)] += KuScattered(k,i);
00438         }
00439     }
00440 
00441   const double gatherTime = gatherTimer.ElapsedTime();
00442   std::cout << "Gathered in " << gatherTime << "\n";
00443   
00444 
00445   const double applyTime = gatherTime + blasTime + scatTime;
00446   std::cout << "Time to do matrix-free product: " << applyTime << std::endl;
00447   
00448   
00449   std::cout << "End Result: TEST PASSED\n";
00450   
00451   // reset format state of std::cout
00452   std::cout.copyfmt(oldFormatState);
00453   
00454   return 0;
00455 }
00456