00001
00002
00003
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
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
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 }