-
Oliver Sander authored
[[Imported from SVN: r3137]]
Oliver Sander authored[[Imported from SVN: r3137]]
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
truncatedsaddlepointgsstep.hh 2.99 KiB
#ifndef TRUNCATED_SADDLE_POINT_GAUSS_SEIDEL_STEP_HH
#define TRUNCATED_SADDLE_POINT_GAUSS_SEIDEL_STEP_HH
#include <dune/common/bitsetvector.hh>
#include <dune/common/fmatrix.hh>
#include <dune/istl/matrix.hh>
#include <dune/solvers/iterationsteps/lineariterationstep.hh>
//! Point block Gauß-Seidel step for saddle point problems.
// Currently only works for scalar problems.
template<class MatrixType, class VectorType, class BitVectorType = Dune::BitSetVector<VectorType::block_type::dimension> >
class TruncatedSaddlePointGSStep
: public LinearIterationStep<MatrixType, VectorType, BitVectorType>
{
public:
//! Default constructor. Doesn't init anything
TruncatedSaddlePointGSStep() {}
//! Constructor with a linear problem
TruncatedSaddlePointGSStep(MatrixType& mat, VectorType& x, VectorType& rhs)
: LinearIterationStep<MatrixType,VectorType>(mat, x, rhs)
{}
virtual VectorType getSol()
{
return *(this->x_);
}
//! Perform one iteration
virtual void iterate()
{
const MatrixType& mat = *this->mat_;
const VectorType& rhs = *this->rhs_;
const BitVectorType& ignore = *this->ignoreNodes_;
VectorType& x = *this->x_;
typedef typename VectorType::block_type VectorBlock;
typedef typename MatrixType::block_type MatrixBlock;
typedef typename MatrixBlock::row_type::const_iterator ColumnIterator;
typedef typename VectorBlock::block_type VectorBlockBlock;
typedef typename MatrixBlock::block_type MatrixBlockBlock;
int P = 2;
// TODO: How must a proper smoother for general problems look like ?
int N = mat[0][0].N();
for(int row=0; row<N; ++row)
{
Dune::BlockVector<VectorBlockBlock> r(P);
Dune::Matrix<MatrixBlockBlock> A(P,P);
Dune::Matrix<Dune::FieldMatrix<char,1,1> > ANZ(P,P);
ANZ = false;
for(int i=0; i<P; ++i)
{
r[i] = rhs[i][row];
for(int j=0; j<P; ++j)
{
ColumnIterator it = mat[i][j][row].begin();
ColumnIterator end = mat[i][j][row].end();
for(; it!=end; ++it)
{
int col = it.index();
if (col == row)
{
A[i][j] = *it;
ANZ[i][j] = true;
}
else
it->mmv(x[j][col],r[i]);
}
}
}
// TODO: How do we generalize this for blocks ?
if (ANZ[0][0])
{
double Z = A[1][0]/A[0][0];
x[1][row] = (r[1] - Z*r[0])/(A[1][1] - Z*A[0][1]);
x[0][row] = (r[0] - A[0][1]*x[1][row])/A[0][0];
}
else
{
x[1][row] = r[1]/A[1][1];
x[0][row] = 0;
}
}
}
};
#endif