projectedblockgsstep.hh

// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=8 sw=4 sts=4:
#ifndef DUNE_PROJECTED_BLOCK_GAUSS_SEIDEL_STEP_HH
#define DUNE_PROJECTED_BLOCK_GAUSS_SEIDEL_STEP_HH

#include <vector>

#include <dune/matrix-vector/axpy.hh>

#include <dune/solvers/common/boxconstraint.hh>
#include <dune/solvers/iterationsteps/lineariterationstep.hh>

template<class MatrixType, class VectorType>
class ProjectedBlockGSStep : public LinearIterationStep<MatrixType, VectorType>
{
    using Base = LinearIterationStep<MatrixType, VectorType>;
    using VectorBlock = typename VectorType::block_type;
    using Field = typename VectorType::field_type;

    enum {BlockSize = VectorBlock::dimension};

public:
    //! Default constructor.  Doesn't init anything
    ProjectedBlockGSStep() {}

    //! Constructor with a linear problem
    ProjectedBlockGSStep(const MatrixType& mat, VectorType& x, const VectorType& rhs)
        : Base(mat, x, rhs)
    {}

    //! Perform one iteration
    virtual void iterate();

    /** \brief Compute the residual of the current iterate of a (block) degree of freedom
     *
     *  \param index Global index of the (block) degree of freedom
     *  \param r Write residual in this vector
     */
    virtual void residual(int index, VectorBlock& r) const {
        r = (*this->rhs_)[index];

        const auto& row = (*this->mat_)[index];
        for (auto cIt = row.begin(); cIt!=row.end(); ++cIt) {
            // r_i -= A_ij x_j
            Dune::MatrixVector::subtractProduct(r, *cIt, (*this->x_)[cIt.index()]);
        }
    }

    using HasObstacle = Dune::BitSetVector<BlockSize>;
    using Obstacle = BoxConstraint<Field, BlockSize>;
    const HasObstacle* hasObstacle_ = nullptr;
    const std::vector<Obstacle>* obstacles_ = nullptr;
};

#include "projectedblockgsstep.cc"

#endif