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