Intrepid
http://trilinos.sandia.gov/packages/docs/r10.12/packages/intrepid/example/Drivers/example_09.cpp
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_CubaturePolylib.hpp"
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 int main(int argc, char *argv[]) {
00115 
00116   //Check number of arguments
00117    if (argc < 4) {
00118       std::cout <<"\n>>> ERROR: Invalid number of arguments.\n\n";
00119       std::cout <<"Usage:\n\n";
00120       std::cout <<"  ./Intrepid_example_Drivers_Example_09.exe deg NX NY verbose\n\n";
00121       std::cout <<" where \n";
00122       std::cout <<"   int deg             - polynomial degree to be used (assumed > 1) \n";
00123       std::cout <<"   int NX              - num intervals in x direction (assumed box domain, 0,1) \n";
00124       std::cout <<"   int NY              - num intervals in y direction (assumed box domain, 0,1) \n";
00125       std::cout <<"   verbose (optional)  - any character, indicates verbose output \n\n";
00126       exit(1);
00127    }
00128   
00129   // This little trick lets us print to std::cout only if
00130   // a (dummy) command-line argument is provided.
00131   int iprint     = argc - 1;
00132   Teuchos::RCP<std::ostream> outStream;
00133   Teuchos::oblackholestream bhs; // outputs nothing
00134   if (iprint > 2)
00135     outStream = Teuchos::rcp(&std::cout, false);
00136   else
00137     outStream = Teuchos::rcp(&bhs, false);
00138   
00139   // Save the format state of the original std::cout.
00140   Teuchos::oblackholestream oldFormatState;
00141   oldFormatState.copyfmt(std::cout);
00142   
00143   *outStream \
00144     << "===============================================================================\n" \
00145     << "|                                                                             |\n" \
00146     << "|  Example: Apply Stiffness Matrix for                                        |\n" \
00147     << "|                   Poisson Equation on Quadrilateral Mesh                    |\n" \
00148     << "|                                                                             |\n" \
00149     << "|  Questions? Contact  Pavel Bochev  (pbboche@sandia.gov),                    |\n" \
00150     << "|                      Denis Ridzal  (dridzal@sandia.gov),                    |\n" \
00151     << "|                      Kara Peterson (kjpeter@sandia.gov).                    |\n" \
00152     << "|                                                                             |\n" \
00153     << "|  Intrepid's website: http://trilinos.sandia.gov/packages/intrepid           |\n" \
00154     << "|  Trilinos website:   http://trilinos.sandia.gov                             |\n" \
00155     << "|                                                                             |\n" \
00156     << "===============================================================================\n";
00157 
00158   
00159   // ************************************ GET INPUTS **************************************
00160   
00161   int deg           = atoi(argv[1]);  // polynomial degree to use
00162   int NX            = atoi(argv[2]);  // num intervals in x direction (assumed box domain, 0,1)
00163   int NY            = atoi(argv[3]);  // num intervals in y direction (assumed box domain, 0,1)
00164   
00165 
00166   // *********************************** CELL TOPOLOGY **********************************
00167   
00168   // Get cell topology for base hexahedron
00169   typedef shards::CellTopology    CellTopology;
00170   CellTopology quad_4(shards::getCellTopologyData<shards::Quadrilateral<4> >() );
00171   
00172   // Get dimensions 
00173   int numNodesPerElem = quad_4.getNodeCount();
00174   int spaceDim = quad_4.getDimension();
00175   
00176   // *********************************** GENERATE MESH ************************************
00177   
00178   *outStream << "Generating mesh ... \n\n";
00179   
00180   *outStream << "   NX" << "   NY\n";
00181   *outStream << std::setw(5) << NX <<
00182     std::setw(5) << NY << "\n\n";
00183   
00184   // Print mesh information
00185   int numElems = NX*NY;
00186   int numNodes = (NX+1)*(NY+1);
00187   *outStream << " Number of Elements: " << numElems << " \n";
00188   *outStream << "    Number of Nodes: " << numNodes << " \n\n";
00189   
00190   // Square
00191   double leftX = 0.0, rightX = 1.0;
00192   double leftY = 0.0, rightY = 1.0;
00193 
00194   // Mesh spacing
00195   double hx = (rightX-leftX)/((double)NX);
00196   double hy = (rightY-leftY)/((double)NY);
00197 
00198   // Get nodal coordinates
00199   FieldContainer<double> nodeCoord(numNodes, spaceDim);
00200   FieldContainer<int> nodeOnBoundary(numNodes);
00201   int inode = 0;
00202   for (int j=0; j<NY+1; j++) {
00203     for (int i=0; i<NX+1; i++) {
00204       nodeCoord(inode,0) = leftX + (double)i*hx;
00205       nodeCoord(inode,1) = leftY + (double)j*hy;
00206       if (j==0 || i==0 || j==NY || i==NX){
00207         nodeOnBoundary(inode)=1;
00208       }
00209       else {
00210         nodeOnBoundary(inode)=0;
00211       }
00212       inode++;
00213     }
00214   }
00215 #define DUMP_DATA
00216 #ifdef DUMP_DATA
00217   // Print nodal coords
00218   ofstream fcoordout("coords.dat");
00219   for (int i=0; i<numNodes; i++) {
00220     fcoordout << nodeCoord(i,0) <<" ";
00221     fcoordout << nodeCoord(i,1) <<"\n";
00222   }
00223   fcoordout.close();
00224 #endif
00225   
00226   
00227   // Element to Node map
00228   // We'll keep it around, but this is only the DOFMap if you are in the lowest order case.
00229   FieldContainer<int> elemToNode(numElems, numNodesPerElem);
00230   int ielem = 0;
00231   for (int j=0; j<NY; j++) {
00232     for (int i=0; i<NX; i++) {
00233       elemToNode(ielem,0) = (NX + 1)*j + i;
00234       elemToNode(ielem,1) = (NX + 1)*j + i + 1;
00235       elemToNode(ielem,2) = (NX + 1)*(j + 1) + i + 1;
00236       elemToNode(ielem,3) = (NX + 1)*(j + 1) + i;
00237       ielem++;
00238     }
00239   }
00240 #ifdef DUMP_DATA
00241   // Output connectivity
00242   ofstream fe2nout("elem2node.dat");
00243   for (int j=0; j<NY; j++) {
00244     for (int i=0; i<NX; i++) {
00245       int ielem = i + j * NX;
00246       for (int m=0; m<numNodesPerElem; m++){
00247         fe2nout << elemToNode(ielem,m) <<"  ";
00248       }
00249       fe2nout <<"\n";
00250     }
00251   }
00252   fe2nout.close();
00253 #endif
00254   
00255 
00256   // ************************************ CUBATURE ************************************** 
00257   *outStream << "Getting cubature ... \n\n";
00258   
00259   // Get numerical integration points and weights
00260   // I only need this on the line since I'm doing tensor products 
00261   DefaultCubatureFactory<double>  cubFactory;                                   
00262 
00263   Teuchos::RCP<Cubature<double,FieldContainer<double>,FieldContainer<double> > > glcub
00264     = Teuchos::rcp(new CubaturePolylib<double,FieldContainer<double>,FieldContainer<double> >(2*deg-1,PL_GAUSS_LOBATTO) );
00265       
00266   const int numCubPoints = glcub->getNumPoints();
00267 
00268   FieldContainer<double> cubPoints1D(numCubPoints, 1);
00269   FieldContainer<double> cubWeights1D(numCubPoints);
00270   
00271   glcub->getCubature(cubPoints1D,cubWeights1D);
00272 
00273   
00274   // ************************************** BASIS ***************************************
00275   *outStream << "Getting basis ... \n\n";
00276   
00277   // Define basis: I only need this on the line also
00278   Basis_HGRAD_LINE_Cn_FEM<double, FieldContainer<double> > lineHGradBasis(deg,POINTTYPE_SPECTRAL);
00279   int numLineFieldsG = lineHGradBasis.getCardinality();
00280   FieldContainer<double> lineGrads(numLineFieldsG, numCubPoints, 1); 
00281   
00282   // Evaluate basis values and gradients at cubature points
00283   lineHGradBasis.getValues(lineGrads, cubPoints1D, OPERATOR_GRAD);
00284 
00285   // ************************************** LTG mapping **********************************
00286   FieldContainer<int> ltgMapping(numElems,numLineFieldsG*numLineFieldsG);
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<numLineFieldsG; m++){
00310         ltgout << ltgMapping(ielem,m) <<"  ";
00311       }
00312       ltgout <<"\n";
00313     }
00314   }
00315   ltgout.close();
00316 #endif
00317   
00318 
00319   // Global arrays in Epetra format 
00320   Epetra_SerialComm Comm;
00321   Epetra_Map globalMapG(numDOF, 0, Comm);
00322 
00323   Epetra_FEVector u(globalMapG);
00324   Epetra_FEVector Ku(globalMapG);
00325 
00326   u.Random();
00327 
00328     
00329   // ************************** Compute element HGrad stiffness matrices *******************************  
00330 //   // Get vertices of all the cells
00331 //   for (int i=0;i<numElems;i++)
00332 //     {
00333 //       for (int j=0;j<4;j++)
00334 //      {
00335 //        const int nodeCur = elemToNode(i,j);
00336 //        for (int k=0;k<spaceDim;k++) 
00337 //          {
00338 //            cellVertices(i,j,k) = nodeCoord(nodeCur,k);
00339 //          }
00340 //      }
00341 //     }
00342 
00343   FieldContainer<double> uScattered(numElems,numLineFieldsG*numLineFieldsG);
00344   FieldContainer<double> KuScattered(numElems,numLineFieldsG*numLineFieldsG);
00345 
00346   // need storage for derivatives of u on each cell
00347   // the number of line dof should be the same as the
00348   // number of cub points.
00349   // This is indexed by Du(q2,q1)
00350   FieldContainer<double> Du(numCubPoints,numCubPoints);
00351 
00352 
00353 
00354   double *uVals = u[0];
00355   double *KuVals = Ku[0];
00356   Epetra_Time scatterTime(Comm);
00357   *outStream << "Scattering\n";
00358   // Scatter
00359   for (int k=0; k<numElems; k++) 
00360     {
00361       for (int i=0;i<numLineFieldsG*numLineFieldsG;i++) 
00362         {
00363           uScattered(k,i) = uVals[ltgMapping(k,i)];
00364         }
00365     }
00366   const double scatTime = scatterTime.ElapsedTime();
00367   *outStream << "Scattered in time " << scatTime << "\n";
00368  
00369   Epetra_Time applyTimer(Comm);
00370   
00371   uScattered.resize(numElems,numLineFieldsG,numLineFieldsG);
00372 
00373   for (int k=0;k<numElems;k++)
00374     {
00375       // local operation: x-derivative term of stiffness matrix
00376       //    evaluate x derivative of u on cell k.
00377       for (int j=0;j<numLineFieldsG;j++)
00378         {
00379           for (int i=0;i<numLineFieldsG;i++)
00380             {
00381               Du(j,i) = 0.0;
00382               for (int q=0;q<numLineFieldsG;q++)
00383                 {
00384                   Du(j,i) += uScattered(k,j,i) * lineGrads(i,q,0);
00385                 }
00386             }
00387         }
00388 
00389       // initialize Ku
00390       for (int i=0;i<numLineFieldsG*numLineFieldsG;i++)
00391         {
00392           KuScattered(k,i) = 0.0;
00393         }
00394 
00395       // loop over basis functions for x term
00396       int cur = 0;
00397       for (int j=0;j<numLineFieldsG;j++)
00398         {
00399           for (int i=0;i<numLineFieldsG;i++)
00400             {
00401               // do the quadrature 
00402               for (int q1=0;q1<numCubPoints;q1++)
00403                 {
00404                   KuScattered(k,cur) += cubWeights1D(j) * cubWeights1D(q1) * Du(j,q1) * lineGrads(i,q1,0);
00405                 }
00406               cur ++;
00407             }
00408         }
00409 
00410       // local operation: y-derivative term of stiffness matrix, store in Du
00411       for (int j=0;j<numLineFieldsG;j++)
00412         {
00413           for (int i=0;i<numLineFieldsG;i++)
00414             {
00415               Du(j,i) = 0.0;
00416               for (int q=0;q<numLineFieldsG;q++)
00417                 {
00418                   Du(j,i) += uScattered(k,j,i) * lineGrads(j,q,0);
00419                 }
00420             }
00421         }
00422 
00423 
00424       // evaluate y-derivatives of u
00425       cur = 0;
00426       for (int j=0;j<numLineFieldsG;j++)
00427         {
00428           for (int i=0;i<numLineFieldsG;i++)
00429             {
00430               for (int q2=0;q2<numCubPoints;q2++)
00431                 {
00432                   KuScattered(k,cur) += cubWeights1D(q2) * Du(q2,i) * lineGrads(j,q2,0) * cubWeights1D(i);
00433                 }
00434             }
00435         }
00436     }
00437   
00438   uScattered.resize(numElems,numLineFieldsG*numLineFieldsG);
00439   
00440   const double applyTime = applyTimer.ElapsedTime();
00441 
00442   *outStream << "Local application: " << applyTime << "\n";
00443       
00444   // gather
00445   Epetra_Time gatherTimer(Comm);
00446   // Gather
00447   for (int k=0;k<numElems;k++)
00448     {
00449       for (int i=0;i<numLineFieldsG*numLineFieldsG;i++)
00450         {
00451           KuVals[ltgMapping(k,i)] += KuScattered(k,i);
00452         }
00453     }
00454 
00455   const double gatherTime = gatherTimer.ElapsedTime();
00456   *outStream << "Gathered in " << gatherTime << "\n";
00457 
00458 
00459 
00460 #ifdef DUMP_DATA
00461    // Dump matrices to disk
00462 //    EpetraExt::RowMatrixToMatlabFile("stiff_matrix.dat",StiffMatrix);
00463 //    EpetraExt::MultiVectorToMatrixMarketFile("rhs_vector.dat",rhs,0,0,false);
00464 #endif
00465 
00466    
00467   std::cout << "End Result: TEST PASSED\n";
00468 
00469    // reset format state of std::cout
00470    std::cout.copyfmt(oldFormatState);
00471    
00472    return 0;
00473 }
00474