Intrepid
http://trilinos.sandia.gov/packages/docs/r10.10/packages/intrepid/example/Drivers/example_07.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_FECrsMatrix.h"
00096 #include "Epetra_FEVector.h"
00097 #include "Epetra_SerialComm.h"
00098 
00099 // Teuchos includes
00100 #include "Teuchos_oblackholestream.hpp"
00101 #include "Teuchos_RCP.hpp"
00102 #include "Teuchos_BLAS.hpp"
00103 
00104 // Shards includes
00105 #include "Shards_CellTopology.hpp"
00106 
00107 // EpetraExt includes
00108 #include "EpetraExt_RowMatrixOut.h"
00109 #include "EpetraExt_MultiVectorOut.h"
00110 
00111 using namespace std;
00112 using namespace Intrepid;
00113 
00114 // Functions to evaluate exact solution and derivatives
00115 double evalu(double & x, double & y, double & z);
00116 int evalGradu(double & x, double & y, double & z, double & gradu1, double & gradu2, double & gradu3);
00117 double evalDivGradu(double & x, double & y, double & z);
00118 
00119 int main(int argc, char *argv[]) {
00120 
00121   //Check number of arguments
00122    if (argc < 4) {
00123       std::cout <<"\n>>> ERROR: Invalid number of arguments.\n\n";
00124       std::cout <<"Usage:\n\n";
00125       std::cout <<"  ./Intrepid_example_Drivers_Example_05.exe deg NX NY verbose\n\n";
00126       std::cout <<" where \n";
00127       std::cout <<"   int deg             - polynomial degree to be used (assumed > 1) \n";
00128       std::cout <<"   int NX              - num intervals in x direction (assumed box domain, 0,1) \n";
00129       std::cout <<"   int NY              - num intervals in y direction (assumed box domain, 0,1) \n";
00130       std::cout <<"   verbose (optional)  - any character, indicates verbose output \n\n";
00131       exit(1);
00132    }
00133   
00134   // This little trick lets us print to std::cout only if
00135   // a (dummy) command-line argument is provided.
00136   int iprint     = argc - 1;
00137   Teuchos::RCP<std::ostream> outStream;
00138   Teuchos::oblackholestream bhs; // outputs nothing
00139   if (iprint > 2)
00140     outStream = Teuchos::rcp(&std::cout, false);
00141   else
00142     outStream = Teuchos::rcp(&bhs, false);
00143   
00144   // Save the format state of the original std::cout.
00145   Teuchos::oblackholestream oldFormatState;
00146   oldFormatState.copyfmt(std::cout);
00147   
00148   *outStream \
00149     << "===============================================================================\n" \
00150     << "|                                                                             |\n" \
00151     << "|  Example: Generate Stiffness Matrix and Right Hand Side Vector for          |\n" \
00152     << "|                   Poisson Equation on Quadrilateral Mesh                    |\n" \
00153     << "|                                                                             |\n" \
00154     << "|  Questions? Contact  Pavel Bochev  (pbboche@sandia.gov),                    |\n" \
00155     << "|                      Denis Ridzal  (dridzal@sandia.gov),                    |\n" \
00156     << "|                      Kara Peterson (kjpeter@sandia.gov).                    |\n" \
00157     << "|                                                                             |\n" \
00158     << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00159     << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00160     << "|                                                                             |\n" \
00161     << "===============================================================================\n";
00162 
00163   
00164   // ************************************ GET INPUTS **************************************
00165   
00166   int deg          = atoi(argv[1]);  // polynomial degree to use
00167   int NX            = atoi(argv[2]);  // num intervals in x direction (assumed box domain, 0,1)
00168   int NY            = atoi(argv[3]);  // num intervals in y direction (assumed box domain, 0,1)
00169   
00170 
00171   // *********************************** CELL TOPOLOGY **********************************
00172   
00173   // Get cell topology for base hexahedron
00174   typedef shards::CellTopology    CellTopology;
00175   CellTopology quad_4(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
00176   
00177   // Get dimensions 
00178   int numNodesPerElem = quad_4.getNodeCount();
00179   int spaceDim = quad_4.getDimension();
00180   
00181   // *********************************** GENERATE MESH ************************************
00182   
00183   *outStream << "Generating mesh ... \n\n";
00184   
00185   *outStream << "   NX" << "   NY\n";
00186   *outStream << std::setw(5) << NX <<
00187     std::setw(5) << NY << "\n\n";
00188   
00189   // Print mesh information
00190   int numElems = NX*NY;
00191   int numNodes = (NX+1)*(NY+1);
00192   *outStream << " Number of Elements: " << numElems << " \n";
00193   *outStream << "    Number of Nodes: " << numNodes << " \n\n";
00194   
00195   // Square
00196   double leftX = 0.0, rightX = 1.0;
00197   double leftY = 0.0, rightY = 1.0;
00198 
00199   // Mesh spacing
00200   double hx = (rightX-leftX)/((double)NX);
00201   double hy = (rightY-leftY)/((double)NY);
00202 
00203   // Get nodal coordinates
00204   FieldContainer<double> nodeCoord(numNodes, spaceDim);
00205   FieldContainer<int> nodeOnBoundary(numNodes);
00206   int inode = 0;
00207   for (int j=0; j<NY+1; j++) {
00208     for (int i=0; i<NX+1; i++) {
00209       nodeCoord(inode,0) = leftX + (double)i*hx;
00210       nodeCoord(inode,1) = leftY + (double)j*hy;
00211       if (j==0 || i==0 || j==NY || i==NX){
00212         nodeOnBoundary(inode)=1;
00213       }
00214       else {
00215         nodeOnBoundary(inode)=0;
00216       }
00217       inode++;
00218     }
00219   }
00220 #define DUMP_DATA
00221 #ifdef DUMP_DATA
00222   // Print nodal coords
00223   ofstream fcoordout("coords.dat");
00224   for (int i=0; i<numNodes; i++) {
00225     fcoordout << nodeCoord(i,0) <<" ";
00226     fcoordout << nodeCoord(i,1) <<"\n";
00227   }
00228   fcoordout.close();
00229 #endif
00230   
00231   
00232   // Element to Node map
00233   // We'll keep it around, but this is only the DOFMap if you are in the lowest order case.
00234   FieldContainer<int> elemToNode(numElems, numNodesPerElem);
00235   int ielem = 0;
00236   for (int j=0; j<NY; j++) {
00237     for (int i=0; i<NX; i++) {
00238       elemToNode(ielem,0) = (NX + 1)*j + i;
00239       elemToNode(ielem,1) = (NX + 1)*j + i + 1;
00240       elemToNode(ielem,2) = (NX + 1)*(j + 1) + i + 1;
00241       elemToNode(ielem,3) = (NX + 1)*(j + 1) + i;
00242       ielem++;
00243     }
00244   }
00245 #ifdef DUMP_DATA
00246   // Output connectivity
00247   ofstream fe2nout("elem2node.dat");
00248   for (int j=0; j<NY; j++) {
00249     for (int i=0; i<NX; i++) {
00250       int ielem = i + j * NX;
00251       for (int m=0; m<numNodesPerElem; m++){
00252         fe2nout << elemToNode(ielem,m) <<"  ";
00253       }
00254       fe2nout <<"\n";
00255     }
00256   }
00257   fe2nout.close();
00258 #endif
00259   
00260 
00261   // ************************************ CUBATURE ************************************** 
00262   *outStream << "Getting cubature ... \n\n";
00263   
00264   // Get numerical integration points and weights
00265   DefaultCubatureFactory<double>  cubFactory;                                   
00266   int cubDegree = 2*deg;
00267   Teuchos::RCP<Cubature<double> > quadCub = cubFactory.create(quad_4, cubDegree); 
00268   
00269   int cubDim       = quadCub->getDimension();
00270   int numCubPoints = quadCub->getNumPoints();
00271   
00272   FieldContainer<double> cubPoints(numCubPoints, cubDim);
00273   FieldContainer<double> cubWeights(numCubPoints);
00274   
00275   quadCub->getCubature(cubPoints, cubWeights);
00276   
00277 
00278   // ************************************** BASIS ***************************************
00279   
00280   *outStream << "Getting basis ... \n\n";
00281   
00282   // Define basis 
00283   Basis_HGRAD_QUAD_Cn_FEM<double, FieldContainer<double> > quadHGradBasis(deg,POINTTYPE_SPECTRAL);
00284   int numFieldsG = quadHGradBasis.getCardinality();
00285   FieldContainer<double> quadGVals(numFieldsG, numCubPoints); 
00286   FieldContainer<double> quadGrads(numFieldsG, numCubPoints, spaceDim); 
00287   
00288   // Evaluate basis values and gradients at cubature points
00289   quadHGradBasis.getValues(quadGVals, cubPoints, OPERATOR_VALUE);
00290   quadHGradBasis.getValues(quadGrads, cubPoints, OPERATOR_GRAD);
00291 
00292   // create the local-global mapping for higher order elements
00293   FieldContainer<int> ltgMapping(numElems,numFieldsG);
00294   const int numDOF = (NX*deg+1)*(NY*deg+1);
00295   ielem=0;
00296   for (int j=0;j<NY;j++) {
00297     for (int i=0;i<NX;i++) {
00298       const int start = deg * j * ( NX * deg + 1 ) + i * deg;
00299       // loop over local dof on this cell
00300       int local_dof_cur=0;
00301       for (int vertical=0;vertical<=deg;vertical++) {
00302         for (int horizontal=0;horizontal<=deg;horizontal++) {
00303           ltgMapping(ielem,local_dof_cur) = start + vertical*(NX*deg+1)+horizontal;
00304           local_dof_cur++;
00305         }
00306       }
00307       ielem++;
00308     }
00309   }
00310 #ifdef DUMP_DATA
00311   // Output ltg mapping
00312   ofstream ltgout("ltg.dat");
00313   for (int j=0; j<NY; j++) {
00314     for (int i=0; i<NX; i++) {
00315       int ielem = i + j * NX;
00316       for (int m=0; m<numFieldsG; m++){
00317         ltgout << ltgMapping(ielem,m) <<"  ";
00318       }
00319       ltgout <<"\n";
00320     }
00321   }
00322   ltgout.close();
00323 #endif
00324   
00325   // ******** CREATE A SINGLE STIFFNESS MATRIX, WHICH IS REPLICATED ON ALL ELEMENTS *********
00326   *outStream << "Building stiffness matrix and right hand side ... \n\n";
00327 
00328   // Settings and data structures for mass and stiffness matrices
00329   typedef CellTools<double>  CellTools;
00330   typedef FunctionSpaceTools fst;
00331   int numCells = 1; 
00332 
00333   // Container for nodes
00334   FieldContainer<double> refQuadNodes(numCells, numNodesPerElem, spaceDim);
00335   // Containers for Jacobian
00336   FieldContainer<double> refQuadJacobian(numCells, numCubPoints, spaceDim, spaceDim);
00337   FieldContainer<double> refQuadJacobInv(numCells, numCubPoints, spaceDim, spaceDim);
00338   FieldContainer<double> refQuadJacobDet(numCells, numCubPoints);
00339   // Containers for element HGRAD stiffness matrix
00340   FieldContainer<double> localStiffMatrix(numCells, numFieldsG, numFieldsG);
00341   FieldContainer<double> weightedMeasure(numCells, numCubPoints);
00342   FieldContainer<double> quadGradsTransformed(numCells, numFieldsG, numCubPoints, spaceDim);
00343   FieldContainer<double> quadGradsTransformedWeighted(numCells, numFieldsG, numCubPoints, spaceDim);
00344   // Containers for right hand side vectors
00345   FieldContainer<double> rhsData(numCells, numCubPoints);
00346   FieldContainer<double> localRHS(numCells, numFieldsG);
00347   FieldContainer<double> quadGValsTransformed(numCells, numFieldsG, numCubPoints);
00348   FieldContainer<double> quadGValsTransformedWeighted(numCells, numFieldsG, numCubPoints);
00349   // Container for cubature points in physical space
00350   FieldContainer<double> physCubPoints(numCells, numCubPoints, cubDim);
00351   
00352   // Global arrays in Epetra format 
00353   // we will explicitly build the sparsity pattern before instantiating the matrix later.
00354   Epetra_SerialComm Comm;
00355   Epetra_Map globalMapG(numDOF, 0, Comm);
00356   Epetra_FEVector u(globalMapG);
00357   Epetra_FEVector Ku(globalMapG);
00358   u.Random();
00359     
00360   // ************************** Compute element HGrad stiffness matrices *******************************  
00361   refQuadNodes(0,0,0) = 0.0;
00362   refQuadNodes(0,0,1) = 0.0;
00363   refQuadNodes(0,1,0) = hx;
00364   refQuadNodes(0,1,1) = 0.0;
00365   refQuadNodes(0,2,0) = hx;
00366   refQuadNodes(0,2,1) = hy;
00367   refQuadNodes(0,3,0) = 0.0;
00368   refQuadNodes(0,3,1) = hy;
00369 
00370   // Compute cell Jacobians, their inverses and their determinants
00371   CellTools::setJacobian(refQuadJacobian, cubPoints, refQuadNodes, quad_4);
00372   CellTools::setJacobianInv(refQuadJacobInv, refQuadJacobian );
00373   CellTools::setJacobianDet(refQuadJacobDet, refQuadJacobian );
00374   
00375   // transform from [-1,1]^2 to [0,hx]x[0,hy]
00376   fst::HGRADtransformGRAD<double>(quadGradsTransformed, refQuadJacobInv, quadGrads);
00377       
00378   // compute weighted measure
00379   fst::computeCellMeasure<double>(weightedMeasure, refQuadJacobDet, cubWeights);
00380 
00381   // multiply values with weighted measure
00382   fst::multiplyMeasure<double>(quadGradsTransformedWeighted,
00383                                weightedMeasure, quadGradsTransformed);
00384 
00385   // integrate to compute element stiffness matrix
00386   fst::integrate<double>(localStiffMatrix,
00387                          quadGradsTransformed, quadGradsTransformedWeighted, COMP_BLAS);
00388 
00389   Epetra_Time graphTimer(Comm);
00390   Epetra_CrsGraph grph( Copy , globalMapG , 4 * numFieldsG );
00391   for (int k=0;k<numElems;k++) 
00392     {
00393       for (int i=0;i<numFieldsG;i++)
00394         {
00395           grph.InsertGlobalIndices(ltgMapping(k,i),numFieldsG,&ltgMapping(k,0));
00396         }
00397     }
00398   grph.FillComplete();
00399   const double graphTime = graphTimer.ElapsedTime();
00400   std::cout << "Graph computed in " << graphTime << "\n";
00401 
00402   Epetra_Time instantiateTimer( Comm );
00403   Epetra_FECrsMatrix StiffMatrix( Copy , grph );
00404   const double instantiateTime = instantiateTimer.ElapsedTime(  );
00405   std::cout << "Matrix instantiated in " << instantiateTime << "\n";
00406 
00407   Epetra_Time assemblyTimer(Comm);
00408 
00409   // *** Element loop ***
00410    for (int k=0; k<numElems; k++) 
00411      {
00412        // assemble into global matrix
00413        StiffMatrix.InsertGlobalValues(numFieldsG,&ltgMapping(k,0),numFieldsG,&ltgMapping(k,0),&localStiffMatrix(0,0,0));
00414 
00415      }
00416 
00417 
00418   // Assemble global matrices
00419    StiffMatrix.GlobalAssemble(); StiffMatrix.FillComplete();
00420 
00421    double assembleTime = assemblyTimer.ElapsedTime();
00422    std::cout << "Time to insert reference element matrix into global matrix: " << assembleTime << std::endl;
00423    std::cout << "Total matrix construction time: " << assembleTime + instantiateTime + graphTime << "\n";
00424    std::cout << "There are " << StiffMatrix.NumGlobalNonzeros() << " nonzeros in the matrix.\n";
00425    std::cout << "There are " << numDOF << " global degrees of freedom.\n";
00426  
00427    Epetra_Time multTimer(Comm);
00428    StiffMatrix.Apply(u,Ku);
00429    double multTime = multTimer.ElapsedTime();
00430    std::cout << "Time to apply: " << multTime << std::endl;
00431 
00432 
00433 #ifdef DUMP_DATA
00434    // Dump matrices to disk
00435 //    EpetraExt::RowMatrixToMatlabFile("stiff_matrix.dat",StiffMatrix);
00436 //    EpetraExt::MultiVectorToMatrixMarketFile("rhs_vector.dat",rhs,0,0,false);
00437 #endif
00438 
00439    std::cout << "End Result: TEST PASSED\n";   
00440 
00441    // reset format state of std::cout
00442    std::cout.copyfmt(oldFormatState);
00443    
00444    return 0;
00445 }
00446 
00447 
00448 // Calculates value of exact solution u
00449  double evalu(double & x, double & y, double & z)
00450  {
00451  /*
00452    // function1
00453     double exactu = sin(M_PI*x)*sin(M_PI*y)*sin(M_PI*z);
00454  */
00455 
00456    // function2
00457    double exactu = sin(M_PI*x)*sin(M_PI*y)*sin(M_PI*z)*exp(x+y+z);
00458 
00459    return exactu;
00460  }
00461 
00462 // Calculates gradient of exact solution u
00463  int evalGradu(double & x, double & y, double & z, double & gradu1, double & gradu2, double & gradu3)
00464  {
00465  /*
00466    // function 1
00467        gradu1 = M_PI*cos(M_PI*x)*sin(M_PI*y)*sin(M_PI*z);
00468        gradu2 = M_PI*sin(M_PI*x)*cos(M_PI*y)*sin(M_PI*z);
00469        gradu3 = M_PI*sin(M_PI*x)*sin(M_PI*y)*cos(M_PI*z);
00470  */
00471 
00472    // function2
00473        gradu1 = (M_PI*cos(M_PI*x)+sin(M_PI*x))
00474                   *sin(M_PI*y)*sin(M_PI*z)*exp(x+y+z);
00475        gradu2 = (M_PI*cos(M_PI*y)+sin(M_PI*y))
00476                   *sin(M_PI*x)*sin(M_PI*z)*exp(x+y+z);
00477        gradu3 = (M_PI*cos(M_PI*z)+sin(M_PI*z))
00478                   *sin(M_PI*x)*sin(M_PI*y)*exp(x+y+z);
00479   
00480    return 0;
00481  }
00482 
00483 // Calculates Laplacian of exact solution u
00484  double evalDivGradu(double & x, double & y, double & z)
00485  {
00486  /*
00487    // function 1
00488     double divGradu = -3.0*M_PI*M_PI*sin(M_PI*x)*sin(M_PI*y)*sin(M_PI*z);
00489  */
00490 
00491    // function 2
00492    double divGradu = -3.0*M_PI*M_PI*sin(M_PI*x)*sin(M_PI*y)*sin(M_PI*z)*exp(x+y+z)
00493                     + 2.0*M_PI*cos(M_PI*x)*sin(M_PI*y)*sin(M_PI*z)*exp(x+y+z)
00494                     + 2.0*M_PI*cos(M_PI*y)*sin(M_PI*x)*sin(M_PI*z)*exp(x+y+z)
00495                     + 2.0*M_PI*cos(M_PI*z)*sin(M_PI*x)*sin(M_PI*y)*exp(x+y+z)
00496                     + 3.0*sin(M_PI*x)*sin(M_PI*y)*sin(M_PI*z)*exp(x+y+z);
00497    
00498    return divGradu;
00499  }
00500