PlayaMatrixLaplacian1D.cpp

00001 /* @HEADER@ */
00002 //   
00003  /* @HEADER@ */
00004 
00005 #include "PlayaMatrixLaplacian1D.hpp"
00006 #include "PlayaEpetraMatrix.hpp"
00007 #include "PlayaIncrementallyConfigurableMatrixFactory.hpp"
00008 
00009 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00010 #include "PlayaVectorImpl.hpp"
00011 #include "PlayaLinearOperatorImpl.hpp"
00012 #endif
00013 
00014 namespace Playa
00015 {
00016 
00017 using namespace Teuchos;
00018 
00019 
00020 MatrixLaplacian1D::MatrixLaplacian1D(int nLocalRows, 
00021   const VectorType<double>& type)
00022   : OperatorBuilder<double>(nLocalRows, type), op_()
00023 {
00024   init(nLocalRows, type);
00025 }
00026 
00027 
00028 
00029 MassMatrix1D::MassMatrix1D(int nLocalRows, 
00030   const VectorType<double>& type)
00031   : OperatorBuilder<double>(nLocalRows, type), op_()
00032 {
00033   init(nLocalRows, type);
00034 }
00035 
00036 
00037 
00038 
00039 void MatrixLaplacian1D::init(int nLocalRows, 
00040   const VectorType<double>& type)
00041 {
00042   int rank = MPIComm::world().getRank();
00043   int nProc = MPIComm::world().getNProc();
00044   RCP<MatrixFactory<double> > mFact;
00045   mFact = vecType().createMatrixFactory(domain(), range());
00046 
00047   if (domain().dim() == domain().numLocalElements())
00048   {
00049     rank = 0;
00050     nProc = 1;
00051   }
00052 
00053   int lowestLocalRow = nLocalRows * rank;
00054 
00055   IncrementallyConfigurableMatrixFactory* icmf 
00056     = dynamic_cast<IncrementallyConfigurableMatrixFactory*>(mFact.get());
00057   if (icmf)
00058   {
00059     for (int i=0; i<nLocalRows; i++)
00060     {
00061       int row = lowestLocalRow + i;
00062       Array<int> colIndices;
00063       if (rank==0 && i==0)
00064       {
00065         colIndices = tuple(row, row+1);
00066       }
00067       else if (rank==nProc-1 && i==nLocalRows-1)
00068       {
00069         colIndices = tuple(row-1, row);
00070       }
00071       else
00072       {
00073         colIndices = tuple(row-1, row, row+1);
00074       }
00075       icmf->initializeNonzerosInRow(row, colIndices.size(),
00076         &(colIndices[0]));
00077     }
00078     icmf->finalize();
00079   }
00080       
00081   op_ = mFact->createMatrix();
00082 
00083   double h = 1.0/((double) domain().dim() + 1);
00084       
00085   RCP<LoadableMatrix<double> > mat = op_.matrix();
00086 
00087   /* fill in with the Laplacian operator */
00088   for (int i=0; i<nLocalRows; i++)
00089   {
00090     int row = lowestLocalRow + i;
00091     Array<int> colIndices;
00092     Array<double> colVals;
00093     if (rank==0 && i==0)
00094     {
00095       colIndices = tuple(row, row+1);
00096       colVals = tuple(2.0/h/h, -1.0/h/h);
00097     }
00098     else if (rank==nProc-1 && i==nLocalRows-1)
00099     {
00100       colIndices = tuple(row-1, row);
00101       colVals = tuple(-1.0/h/h, 2.0/h/h);
00102     }
00103     else
00104     {
00105       colIndices = tuple(row-1, row, row+1);
00106       colVals = tuple(-1.0/h/h, 2.0/h/h, -1.0/h/h);
00107     }
00108     mat->addToRow(row, colIndices.size(), 
00109       &(colIndices[0]), &(colVals[0]));
00110   }
00111 }
00112 
00113 
00114 void MassMatrix1D::init(int nLocalRows, 
00115   const VectorType<double>& type)
00116 {
00117   int rank = MPIComm::world().getRank();
00118   int nProc = MPIComm::world().getNProc();
00119   RCP<MatrixFactory<double> > mFact;
00120   mFact = vecType().createMatrixFactory(domain(), range());
00121 
00122   if (domain().dim() == domain().numLocalElements())
00123   {
00124     rank = 0;
00125     nProc = 1;
00126   }
00127 
00128   int lowestLocalRow = nLocalRows * rank;
00129 
00130   IncrementallyConfigurableMatrixFactory* icmf 
00131     = dynamic_cast<IncrementallyConfigurableMatrixFactory*>(mFact.get());
00132   if (icmf)
00133   {
00134     for (int i=0; i<nLocalRows; i++)
00135     {
00136       int row = lowestLocalRow + i;
00137       Array<int> colIndices;
00138       if (rank==0 && i==0)
00139       {
00140         colIndices = tuple(row, row+1);
00141       }
00142       else if (rank==nProc-1 && i==nLocalRows-1)
00143       {
00144         colIndices = tuple(row-1, row);
00145       }
00146       else
00147       {
00148         colIndices = tuple(row-1, row, row+1);
00149       }
00150       icmf->initializeNonzerosInRow(row, colIndices.size(),
00151         &(colIndices[0]));
00152     }
00153     icmf->finalize();
00154   }
00155       
00156   op_ = mFact->createMatrix();
00157 
00158   double h = 1.0/((double) domain().dim() + 1);
00159       
00160   RCP<LoadableMatrix<double> > mat = op_.matrix();
00161 
00162   /* fill in with the mass operator */
00163   for (int i=0; i<nLocalRows; i++)
00164   {
00165     int row = lowestLocalRow + i;
00166     Array<int> colIndices;
00167     Array<double> colVals;
00168     if (rank==0 && i==0)
00169     {
00170       colIndices = tuple(row, row+1);
00171       colVals = tuple(2.0/3.0, 1.0/3.0);
00172     }
00173     else if (rank==nProc-1 && i==nLocalRows-1)
00174     {
00175       colIndices = tuple(row-1, row);
00176       colVals = tuple(1.0/3.0, 2.0/3.0);
00177     }
00178     else
00179     {
00180       colIndices = tuple(row-1, row, row+1);
00181       colVals = tuple(1.0/3.0, 2.0/3.0, 1.0/3.0);
00182     }
00183     mat->addToRow(row, colIndices.size(), 
00184       &(colIndices[0]), &(colVals[0]));
00185   }
00186 }
00187 
00188 
00189 }

doxygen