levelpatchpreconditioner.hh

#ifndef LEVEL_PATCH_PRECONDITIONER_HH
#define LEVEL_PATCH_PRECONDITIONER_HH

#include <string>

#include <dune/common/timer.hh>
#include <dune/common/fvector.hh>
#include <dune/common/bitsetvector.hh>

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

#include "../../data-structures/levelcontactnetwork.hh"

#include "supportpatchfactory.hh"

#include <dune/localfunctions/lagrange/pqkfactory.hh>

#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/functiontools/boundarydofs.hh>


template <class LevelContactNetwork, class PatchSolver, class MatrixType, class VectorType>
class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorType> {

public:
    enum Mode {additive, multiplicative};
    static const int dim = LevelContactNetwork::dim;

private:
    using Base = LinearIterationStep<MatrixType, VectorType>;

    using ctype = typename LevelContactNetwork::ctype;

    using PatchFactory = SupportPatchFactory<LevelContactNetwork>;
    using Patch = typename PatchFactory::Patch;

    const Mode mode_;

    const LevelContactNetwork& levelContactNetwork_;
    const LevelContactNetwork& fineContactNetwork_;

    const int level_;

    PatchFactory patchFactory_;
    std::vector<Patch> patches_;

    std::shared_ptr<PatchSolver> patchSolver_;
    size_t patchDepth_;

public:

    // for each coarse patch given by levelContactNetwork: set up local problem, compute local correction
    LevelPatchPreconditioner(const LevelContactNetwork& levelContactNetwork,
                             const LevelContactNetwork& fineContactNetwork,
                             const Mode mode = LevelPatchPreconditioner::Mode::additive) :
            mode_(mode),
            levelContactNetwork_(levelContactNetwork),
            fineContactNetwork_(fineContactNetwork),
            level_(levelContactNetwork_.level()),
            patchFactory_(levelContactNetwork_, fineContactNetwork_) {

        setPatchDepth();
        this->verbosity_ = QUIET;
    }

    // build support patches
    void build() {
        size_t totalNVertices = 0;
        for (size_t i=0; i<levelContactNetwork_.nBodies(); i++) {
            totalNVertices += levelContactNetwork_.body(i)->nVertices();
        }

        patches_.resize(totalNVertices);

        // init local fine level corrections
        Dune::Timer timer;
        if (this->verbosity_ == FULL) {
            std::cout << std::endl;
            std::cout << "---------------------------------------------" << std::endl;
            std::cout << "Initializing local fine grid corrections! Level: " << level_ << std::endl;


            timer.reset();
            timer.start();
        }

        Dune::BitSetVector<1> vertexVisited(totalNVertices);
        vertexVisited.unsetAll();

        const auto& levelIndices = patchFactory_.coarseIndices();

        for (size_t bodyIdx=0; bodyIdx<levelContactNetwork_.nBodies(); bodyIdx++) {
            const auto& gridView = levelContactNetwork_.body(bodyIdx)->gridView();

            for (const auto& e : elements(gridView)) {
                const auto& refElement = Dune::ReferenceElements<double, dim>::general(e.type());

                for (size_t i=0; i<refElement.size(dim); i++) {
                    auto globalIdx = levelIndices.vertexIndex(bodyIdx, e, i);

                    if (!vertexVisited[globalIdx][0]) {
                        vertexVisited[globalIdx][0] = true;
                        patchFactory_.build(bodyIdx, e, i, patches_[globalIdx], patchDepth_);
                    }
                }
            }
        }

        if (this->verbosity_ == FULL) {
            std::cout << std::endl;
            std::cout << "Total setup time: " << timer.elapsed() << " seconds." << std::endl;
            std::cout << "---------------------------------------------" << std::endl;
            timer.stop();
        }

        if (this->verbosity_ != QUIET) {
            std::cout << "Level: " << level_ << " LevelPatchPreconditioner::build() - SUCCESS" << std::endl;
        }
    }

    void setPatchSolver(PatchSolver&& patchSolver) {
        patchSolver_ = Dune::Solvers::wrap_own_share<PatchSolver>(std::forward<PatchSolver>(patchSolver));
    }

    void setPatchDepth(const size_t patchDepth = 0) {
        patchDepth_ = patchDepth;
    }

    void setVerbosity(VerbosityMode verbosity) {
        this->verbosity_ = verbosity;
    }

    virtual void iterate() {
        if (mode_ == additive)
            iterateAdd();
        else
            iterateMult();
    }

    void iterateAdd() {
        *(this->x_) = 0;	
        VectorType x = *(this->x_);

        for (const auto& p : patches_) {
            x = 0;

            auto& step = patchSolver_->getIterationStep();
            step.setProblem(this->mat_, x, this->rhs_);
            step.setIgnore(p);

            patchSolver_->check();
            patchSolver_->preprocess();
            patchSolver_->solve();

            *(this->x_) += patchSolver_->getSol();
        }
    }

    void iterateMult() {
        *(this->x_) = 0;
/*
        VectorType x = *(this->x_);
        for (size_t i=0; i<patches_.size(); i++) {
            VectorType updatedRhs(*(this->rhs_));
            this->mat_->mmv(*(this->x_), updatedRhs);

            patchSolver_.setProblem(*this->mat_, *(this->x_), updatedRhs);
            patchSolver_.setIgnore(patches_[i]);
            patchSolver_.solve();

            *(this->x_) += x;
        } */
    }

    size_t size() const {
        return patches_.size();
    }

    size_t level() const {
        return level_;
    }

    size_t fineLevel() const {
        return fineContactNetwork_.level();
    }
};

#endif