PlayaPoissonBoltzmannJacobian.cpp

00001 /* @HEADER@ */
00002 //   
00003 /* @HEADER@ */
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   /* fill in with the Laplacian operator plus exp(-x) */
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 }