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