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