projectedblockgsstep.cc

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

template<class MatrixType, class VectorType>
inline
void ProjectedBlockGSStep<MatrixType, VectorType>::iterate()
{
    if (hasObstacle_->size()!= (unsigned int)this->x_->size())
        DUNE_THROW(SolverError, "Size of hasObstacle (" << hasObstacle_->size()
                   << ") doesn't match solution vector (" << this->x_->size() << ")");

    const MatrixType& mat = *this->mat_;

    for (size_t i=0; i<this->x_->size(); i++) {
        auto ignoreCount = (*this->ignoreNodes_)[i].count();
        if (ignoreCount == BlockSize)
            continue;
        if (ignoreCount > 0)
          DUNE_THROW(Dune::NotImplemented,
                     "Individual blocks must be either ignored completely, or not at all");

        bool zeroDiagonal = false;
        for (size_t j=0; j<BlockSize; j++) {
            // When using this solver as part of a truncated multigrid solver,
            // the diagonal entries of the matrix may get completely truncated
            // away.  In this case we just do nothing here.
            zeroDiagonal |= (mat[i][i][j][j] < 1e-10);
        }

        if (zeroDiagonal)
            continue;

        VectorBlock r;
        this->residual(i, r);

        // Compute x_i += A_{i,i}^{-1} r[i]
        VectorBlock v;
        VectorBlock& x = (*this->x_)[i];


        if ((*hasObstacle_)[i].none()) {

            // No obstacle --> solve linear problem
            mat[i][i].solve(v, r);

        } else {

            // Solve the local constraint minimization problem
            // We use a projected Gauss-Seidel, for lack of anything better
            Obstacle defectObstacle = (*obstacles_)[i];
            defectObstacle -= x;

            // Initial correction
            v = 0;

            for (size_t j=0; j< ((BlockSize==1) ? 1 : 20); j++) {

                for (size_t k=0; k<BlockSize; k++) {

                    // Compute residual
                    Field sr = 0;
                    for (size_t l=0; l<BlockSize; l++)
                        sr += mat[i][i][k][l] * v[l];

                    sr = r[k] - sr;
                    v[k] += sr / mat[i][i][k][k];

                    // Project
                    if ((*hasObstacle_)[i][k]==true) {
                        if (v[k] < defectObstacle.lower(k))
                            v[k] = defectObstacle.lower(k);
                        else if (v[k] > defectObstacle.upper(k))
                            v[k] = defectObstacle.upper(k);
                    }

                }

            }

        }

        // Add correction
        x += v;

    }

}