Skip to content
Snippets Groups Projects
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