00001
00002
00003
00004
00005 #include "PlayaHeatOperator1D.hpp"
00006 #include "PlayaIncrementallyConfigurableMatrixFactory.hpp"
00007
00008 #ifndef HAVE_TEUCHOS_EXPLICIT_INSTANTIATION
00009 #include "PlayaVectorImpl.hpp"
00010 #include "PlayaLinearOperatorImpl.hpp"
00011 #endif
00012
00013 using namespace Playa;
00014 using namespace Teuchos;
00015
00016
00017 HeatOperator1D::HeatOperator1D(int nLocalRows, const VectorType<double>& type)
00018 : OperatorBuilder<double>(nLocalRows, type),
00019 matrixFactory_(),
00020 cj_(0.0),
00021 h_(0.0),
00022 nLocalRows_(nLocalRows),
00023 lowestLocalRow_(0)
00024 {
00025 int rank = MPIComm::world().getRank();
00026 int nProc = MPIComm::world().getNProc();
00027 matrixFactory_ = vecType().createMatrixFactory(domain(), range());
00028
00029 int n = domain().dim();
00030 h_ = 1.0/(n - 1.0);
00031
00032 lowestLocalRow_ = nLocalRows * rank;
00033
00034 IncrementallyConfigurableMatrixFactory* icmf
00035 = dynamic_cast<IncrementallyConfigurableMatrixFactory*>(matrixFactory_.get());
00036 for (int i=0; i<nLocalRows; i++)
00037 {
00038 int row = lowestLocalRow_ + i;
00039 Array<int> colIndices;
00040 if ((rank==0 && i==0) || (rank==nProc-1 && i==nLocalRows-1))
00041 {
00042 colIndices = tuple(row);
00043 }
00044 else
00045 {
00046 colIndices = tuple(row-1, row, row+1);
00047 }
00048 icmf->initializeNonzerosInRow(row, colIndices.size(),
00049 &(colIndices[0]));
00050 }
00051 icmf->finalize();
00052 }
00053
00054
00055
00056 LinearOperator<double> HeatOperator1D::getOp() const
00057 {
00058
00059 int rank = MPIComm::world().getRank();
00060 int nProc = MPIComm::world().getNProc();
00061 LinearOperator<double> rtn = matrixFactory_->createMatrix();
00062
00063 RCP<LoadableMatrix<double> > mat = rtn.matrix();
00064
00065
00066 for (int i=0; i<nLocalRows_; i++)
00067 {
00068 int row = lowestLocalRow_ + i;
00069 Array<int> colIndices;
00070 Array<double> colVals;
00071 if ((rank==0 && i==0) || (rank==nProc-1 && i==nLocalRows_-1))
00072 {
00073 colIndices = tuple(row);
00074 colVals = tuple(1.0);
00075 }
00076 else
00077 {
00078 colIndices = tuple(row-1, row, row+1);
00079 colVals = tuple(-1.0/h_/h_, 2.0/h_/h_ + cj_, -1.0/h_/h_);
00080 }
00081 mat->addToRow(row, colIndices.size(),
00082 &(colIndices[0]), &(colVals[0]));
00083 }
00084
00085 return rtn;
00086 }