dune/tectonic/spatial-solving/preconditioners/CMakeLists.txt

+add_custom_target(tectonic_dune_spatial-solving_preconditioners SOURCES
+  hierarchicleveliterator.hh
+  levelpatchpreconditioner.hh
+  patchproblem.hh
+  multilevelpatchpreconditioner.hh
+  supportpatchfactory.hh
+)
+
+#install headers
+install(FILES
+  hierarchicleveliterator.hh
+  levelpatchpreconditioner.hh
+  patchproblem.hh
+  multilevelpatchpreconditioner.hh
+  supportpatchfactory.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)

dune/tectonic/spatial-solving/preconditioners/hierarchicleveliterator.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set ts=8 sw=4 et sts=4:
+#ifndef HIERARCHIC_LEVEL_ITERATOR_HH
+#define HIERARCHIC_LEVEL_ITERATOR_HH
+
+#include <dune/common/fvector.hh>
+#include <dune/common/iteratorfacades.hh>
+
+#include <dune/geometry/multilineargeometry.hh>
+#include <dune/geometry/referenceelements.hh>
+
+/** \brief Hierarchic leaf iterator.
+ *
+ *  This iterator loops over all children of a given coarse element and only returns the ones that are leaf.
+ *  If the starting element is leaf itself, then the returned iterator returns the element itself.
+ *  This class also provides a geometry, mapping local coordinates of the children to local coordinates
+ *  in the coarse element. 
+*/
+template <class GridImp>
+class HierarchicLevelIterator :
+    public Dune::ForwardIteratorFacade<HierarchicLevelIterator<GridImp>, const typename GridImp::template Codim<0>::Entity>
+{
+    typedef typename GridImp::template Codim<0>::Entity Element;
+    typedef typename GridImp::HierarchicIterator HierarchicIterator;
+    enum {dim = GridImp::dimension};
+    enum {dimworld = GridImp::dimensionworld};
+public:
+    typedef Dune::CachedMultiLinearGeometry<typename GridImp::ctype, dim, dimworld> LocalGeometry;
+
+    enum PositionFlag {begin, end};
+
+    HierarchicLevelIterator(const GridImp& grid, const Element& element, 
+                            PositionFlag flag, int maxlevel, bool nested = true)
+        : element_(element), maxlevel_(maxlevel), hIt_(element.hend(maxlevel_)),
+            flag_(flag), nested_(nested)
+    {
+
+        // if the element itself is leaf, then we don't have to iterate over the children
+        if (flag==begin && !(element_.level()==maxlevel_)) {
+            hIt_ = element_.hbegin(maxlevel_);
+
+            //NOTE This class by now assumes that possible non-nestedness of the grid levels only arises
+            // due to boundary parametrisation
+            if (!nested_)  {
+                // Check if the element is a boundary element, and set the nested flag correspondingly
+                // If the element is not at the boundary, then we can neglect the nested flag
+                typename GridImp::LevelIntersectionIterator lIt = grid.levelGridView(element.level()).ibegin(element);
+                typename GridImp::LevelIntersectionIterator lEndIt = grid.levelGridView(element.level()).ibegin(element);
+                nested_ = true;
+                for (; lIt != lEndIt; ++lIt)
+                    if (lIt->boundary()) {
+                        nested_ = false;
+                        break;
+                    }
+            }
+
+            if (!(hIt_->level()==maxlevel_))
+                increment();
+        }
+    }
+
+    //! Copy constructor
+    HierarchicLevelIterator(const HierarchicLevelIterator<GridImp>& other) :
+        element_(other.element_),maxlevel_(other.maxlevel_), hIt_(other.hIt_),
+        flag_(other.flag_), nested_(other.nested_)
+    {}
+
+    //! Equality
+    bool equals(const HierarchicLevelIterator<GridImp>& other) const {
+        return  (element_ == other.element_)  && 
+                ((flag_==other.flag_ && flag_==end) || (hIt_ == other.hIt_ && flag_==begin && other.flag_==flag_));
+    }
+
+    //! Prefix increment
+    void increment() {
+
+        HierarchicIterator hEndIt = element_.hend(maxlevel_);
+
+        if (hIt_ == hEndIt) {
+            flag_ = end;
+            return;
+        }
+
+        // Increment until we hit a leaf child element
+        do {
+            ++hIt_;
+        } while(hIt_ != hEndIt && (!(hIt_->level()==maxlevel_)));
+    
+        if (hIt_ == hEndIt)
+            flag_ = end;
+    }
+
+    //! Dereferencing
+    const Element& dereference() const {
+        if (flag_ == end)
+            DUNE_THROW(Dune::Exception,"HierarchicLevelIterator: Cannot dereference end iterator!");
+        return (element_.level()==maxlevel_) ? element_ : *hIt_;
+    }
+
+    //! Compute the local geometry mapping from the leaf child to the starting element
+    LocalGeometry geometryInAncestor() {
+
+        //NOTE: We assume here that the type of the child and the ancestors are the same!
+        const Dune::ReferenceElement<typename GridImp::ctype, dim>& ref = Dune::ReferenceElements<typename GridImp::ctype,dim>::general(element_.type());
+
+        // store local coordinates of the leaf child element within the coarse starting element
+        std::vector<Dune::FieldVector<typename GridImp::ctype,dim> > corners(ref.size(dim));
+        for (int i=0; i<corners.size(); i++)
+            corners[i] = ref.position(i,dim);
+
+        // If the element is leaf, then return an Identity mapping
+        if (element_.level()==maxlevel_)
+            return LocalGeometry(ref,corners);
+
+        if (nested_) {
+            const typename Element::Geometry leafGeom = hIt_->geometry();
+            const typename Element::Geometry coarseGeom = element_.geometry();
+
+            for (int i=0; i<corners.size();i++) 
+                corners[i] = coarseGeom.local(leafGeom.corner(i));
+            return LocalGeometry(ref,corners);
+        }
+ 
+        Element father(*hIt_);
+        while (father != element_) {
+
+            const typename Element::LocalGeometry fatherGeom = hIt_->geometryInFather();
+            father = father->father();
+
+            for (int i=0; i<corners.size(); i++)
+                corners[i] = fatherGeom.global(corners[i]);
+        }
+
+        return LocalGeometry(ref,corners);
+    }
+
+private:
+
+    //! The starting element
+    const Element element_;
+
+    //! The grids maxlevel
+    int maxlevel_;
+
+    //! The actual hierarchic iterator
+    HierarchicIterator hIt_;
+
+    //! Position flag
+    PositionFlag flag_;
+
+    //! Flag that is set true if the grid levels are conforming (i.e. no parametrised boundaries)
+    bool nested_;
+
+};
+#endif

dune/tectonic/spatial-solving/preconditioners/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 <dune/solvers/common/numproc.hh>
+
+#include "../../data-structures/network/levelcontactnetwork.hh"
+
+#include "patchproblem.hh"
+#include "supportpatchfactory.hh"
+
+#include <dune/localfunctions/lagrange/pqkfactory.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/functiontools/boundarydofs.hh>
+
+enum MPPMode {additive, multiplicative};
+inline
+std::istream& operator>>(std::istream& lhs, MPPMode& e)
+{
+  std::string s;
+  lhs >> s;
+
+  if (s == "additive" || s == "1")
+    e = MPPMode::additive;
+  else if (s == "multiplicative" || s == "0")
+    e = MPPMode::multiplicative;
+  else
+    lhs.setstate(std::ios_base::failbit);
+
+  return lhs;
+}
+
+template <class LevelContactNetwork, class PatchSolver, class MatrixType, class VectorType>
+class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorType> {
+
+public:
+    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;
+    using PatchProblem = PatchProblem<MatrixType, VectorType>;
+
+    const MPPMode mode_;
+
+    const LevelContactNetwork& levelContactNetwork_;
+    const LevelContactNetwork& fineContactNetwork_;
+
+    const int level_;
+
+    PatchFactory patchFactory_;
+    std::vector<Patch> patches_;
+    std::vector<std::unique_ptr<PatchProblem>> patchProblems_;
+    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 MPPMode mode = additive) :
+            mode_(mode),
+            levelContactNetwork_(levelContactNetwork),
+            fineContactNetwork_(fineContactNetwork),
+            level_(levelContactNetwork_.level()),
+            patchFactory_(levelContactNetwork_, fineContactNetwork_) {
+
+        setPatchDepth();
+        this->verbosity_ = NumProc::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_ == NumProc::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();
+
+            //printDofLocation(gridView);
+
+            for (const auto& e : elements(gridView)) {
+                const auto& refElement = Dune::ReferenceElements<double, dim>::general(e.type());
+
+                for (int 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_ == NumProc::FULL) {
+            std::cout << std::endl;
+            std::cout << "Total setup time: " << timer.elapsed() << " seconds." << std::endl;
+            std::cout << "---------------------------------------------" << std::endl;
+            timer.stop();
+        }
+
+        if (this->verbosity_ != NumProc::QUIET) {
+            std::cout << "Level: " << level_ << " LevelPatchPreconditioner::build() - SUCCESS" << std::endl;
+        }
+    }
+
+    void setPatchSolver(std::shared_ptr<PatchSolver> patchSolver) {
+        patchSolver_ = patchSolver;
+    }
+
+    void setPatchDepth(const size_t patchDepth = 0) {
+        patchDepth_ = patchDepth;
+    }
+
+    void setVerbosity(NumProc::VerbosityMode verbosity) {
+        this->verbosity_ = verbosity;
+    }
+
+    void setMatrix(const MatrixType& mat) override {
+        Base::setMatrix(mat);
+
+        patchProblems_.resize(patches_.size());
+        for (size_t i=0; i<patches_.size(); i++) {
+            patchProblems_[i] = std::make_unique<PatchProblem>(mat, patches_[i]);
+        }
+
+        //std::cout << "matrix set!" << std::endl;
+    }
+
+    void iterate() override {
+        if (mode_ == additive)
+            iterateAdd();
+        else
+            iterateMult();
+    }
+
+    void iterateAdd() {
+        *(this->x_) = 0;	
+        VectorType x = *(this->x_);
+
+        Dune::Timer timer;
+        timer.start();
+
+        size_t systemSize = 0;
+        size_t count = 0;
+
+        //std::cout << "level::iterate() ... patches: " << patches_.size() << " level size: " << x.size() << std::endl;
+
+        for (size_t i=0; i<patches_.size(); i++) {
+            x = 0;
+
+            /*
+            const auto& p = patches_[i];
+
+            auto ignore = this->ignore();
+            for (size_t i=0; i<ignore.size(); i++) {
+                for (size_t d=0; d<dim; d++) {
+                   if (p[i][d])
+                       ignore[i][d] = true;
+                }
+            }
+
+            auto& step = patchSolver_->getIterationStep();
+            dynamic_cast<LinearIterationStep<MatrixType, VectorType>&>(step).setProblem(*this->mat_, x, *this->rhs_);
+            step.setIgnore(ignore);*/
+
+            const auto& patchMat = patchProblems_[i]->mat();
+
+            patchProblems_[i]->setRhs(*this->rhs_);
+            const auto& patchRhs = patchProblems_[i]->rhs();
+
+            VectorType patchX(patchMat.M());
+            patchX = 0;
+
+            auto& step = patchSolver_->getIterationStep();
+            dynamic_cast<LinearIterationStep<MatrixType, VectorType>&>(step).setProblem(patchMat, patchX, patchRhs);
+
+            // empty ignore
+            Dune::Solvers::DefaultBitVector_t<VectorType> ignore(patchX.size());
+            ignore.unsetAll();
+            step.setIgnore(ignore);
+
+            patchSolver_->check();
+            patchSolver_->preprocess();
+            patchSolver_->solve();
+
+            patchProblems_[i]->prolong(patchX, x);
+
+            *(this->x_) += x;
+
+            /*if (count*1.0/patches_.size() >= 0.1) {
+                std::cout << (int) (i*1.0/patches_.size()*100) << " %. Elapsed time: " << timer.elapsed() << std::endl;
+                count = 0;
+            }
+            count++;
+            systemSize += patchX.size();*/
+        }
+
+       /* timer.stop();
+
+        std::cout << "Total elapsed time: " << timer.elapsed() << std::endl;
+        std::cout << "Average time per patch: " << timer.elapsed()*1.0/patches_.size() << std::endl;
+        std::cout << "Average patch size: " << systemSize*1.0/patches_.size() << std::endl;
+        std::cout << "-------------------------------" << std::endl << std::endl;*/
+    }
+
+    void iterateMult() {
+        *(this->x_) = 0;
+
+        DUNE_THROW(Dune::Exception, "Not implemented!");
+
+/*
+        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
+

dune/tectonic/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh

+#ifndef MULTILEVEL_PATCH_PRECONDITIONER_HH
+#define MULTILEVEL_PATCH_PRECONDITIONER_HH
+
+#include <string>
+
+#include <dune/common/timer.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/parametertree.hh>
+
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+#include <dune/solvers/iterationsteps/blockgssteps.hh>
+#include <dune/solvers/iterationsteps/cgstep.hh>
+#include <dune/solvers/iterationsteps/lineariterationstep.hh>
+#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
+
+#include "../contact/nbodycontacttransfer.hh"
+#include "levelpatchpreconditioner.hh"
+
+template <class ContactNetwork, class MatrixType, class VectorType>
+class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorType> {
+private:
+    using Base = LinearIterationStep<MatrixType, VectorType>;
+
+    using PatchSolver = LoopSolver<VectorType>;
+    using PatchSolverStep = TruncatedBlockGSStep<MatrixType, VectorType>;
+    using Norm = EnergyNorm<MatrixType, VectorType>;
+
+    using LevelContactNetwork = typename ContactNetwork::LevelContactNetwork;
+    using LevelPatchPreconditioner = LevelPatchPreconditioner<LevelContactNetwork, PatchSolver, MatrixType, VectorType>;
+    using BitVector = typename LinearIterationStep<MatrixType, VectorType>::BitVector;
+
+    using MGTransfer = NBodyContactTransfer<ContactNetwork, VectorType>;
+
+    static const int dim = ContactNetwork::dim;
+
+    const ContactNetwork& contactNetwork_;
+    const Dune::BitSetVector<1>& activeLevels_;
+
+    const MPPMode mode_;
+
+    std::vector<std::shared_ptr<LevelPatchPreconditioner>> levelPatchPreconditioners_;
+
+    std::vector<MatrixType> levelMat_;
+    std::vector<VectorType> levelX_;
+    std::vector<VectorType> levelRhs_;
+
+    std::vector<BitVector> ignoreHierarchy_;
+
+    std::vector<std::shared_ptr<EnergyNorm<MatrixType, VectorType>>> levelErrorNorms_;
+    std::vector<std::shared_ptr<LinearIterationStep<MatrixType, VectorType>>> levelItSteps_;
+    std::vector<std::shared_ptr<PatchSolver>> levelSolvers_;
+
+    //std::vector<BitVector> recompute_;
+    std::vector<std::shared_ptr<MGTransfer>> mgTransfer_;
+
+public:
+    MultilevelPatchPreconditioner(const Dune::ParameterTree& parset,
+                                  const ContactNetwork& contactNetwork,
+                                  const Dune::BitSetVector<1>& activeLevels) :
+          contactNetwork_(contactNetwork),
+          activeLevels_(activeLevels),
+          mode_(parset.get<MPPMode>("mode")){
+
+        assert(activeLevels.size() == contactNetwork_.nLevels());
+
+        // init level patch preconditioners and multigrid transfer
+        levelPatchPreconditioners_.resize(1);
+        levelPatchPreconditioners_[0] = nullptr; // blank level representing global coarse problem to make further indexing easier
+
+        int maxLevel = -1;
+        for (size_t i=0; i<activeLevels_.size(); i++) {
+            if (activeLevels_[i][0]) {
+                maxLevel = i;
+                break;
+            }
+        }
+
+        for (size_t i=maxLevel+1; i<activeLevels_.size(); i++) {
+            if (activeLevels_[i][0]) {
+                // init local patch preconditioners on each level
+                const auto& fineNetwork = *contactNetwork_.level(i);
+                levelPatchPreconditioners_.push_back(std::make_shared<LevelPatchPreconditioner>(*contactNetwork_.level(maxLevel), fineNetwork, mode_));
+                levelPatchPreconditioners_.back()->setPatchDepth(parset.get<size_t>("patchDepth"));
+                maxLevel = i;
+            }
+        }
+
+        levelMat_.resize(size()-1);
+        levelX_.resize(size());
+        levelRhs_.resize(size());
+
+        // init patch solvers
+        levelSolvers_.resize(size());
+        levelItSteps_.resize(size());
+        levelErrorNorms_.resize(size());
+
+        // set basesolver
+        levelItSteps_[0] = std::make_shared<PatchSolverStep>();
+        levelErrorNorms_[0] = std::make_shared<Norm>(*levelItSteps_[0].get());
+        levelSolvers_[0] = std::make_shared<PatchSolver>(*levelItSteps_[0].get(),
+                                   parset.get<size_t>("basesolver.maximumIterations"),
+                                   parset.get<double>("basesolver.tolerance"),
+                                   *levelErrorNorms_[0].get(),
+                                   parset.get<Solver::VerbosityMode>("basesolver.verbosity"));
+
+        for (size_t i=1; i<levelSolvers_.size(); i++) {
+            levelItSteps_[i] = std::make_shared<PatchSolverStep>();
+            levelErrorNorms_[i] = std::make_shared<Norm>(*levelItSteps_[i].get());
+            levelSolvers_[i] = std::make_shared<PatchSolver>(*levelItSteps_[i].get(),
+                                       parset.get<size_t>("patchsolver.maximumIterations"),
+                                       parset.get<double>("patchsolver.tolerance"),
+                                       *levelErrorNorms_[i].get(),
+                                       parset.get<Solver::VerbosityMode>("patchsolver.verbosity"));
+            levelPatchPreconditioners_[i]->setPatchSolver(levelSolvers_[i]);
+        }
+
+        // init multigrid transfer
+        mgTransfer_.resize(levelPatchPreconditioners_.size()-1);
+        for (size_t i=0; i<mgTransfer_.size(); i++) {
+            mgTransfer_[i] = std::make_shared<MGTransfer>();
+            mgTransfer_[i]->setup(contactNetwork_, levelPatchPreconditioners_[i+1]->level(), levelPatchPreconditioners_[i+1]->fineLevel());
+        }
+
+        // from dune/contact/solvers/nsnewtonmgstep.cc
+        /*// Remove any recompute filed so that initially the full transferoperator is assembled
+        for (size_t i=0; i<mgTransfer_.size(); i++)
+            std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(mgTransfer_[i])->setRecomputeBitField(nullptr);
+
+        // specify the subset of entries to be reassembled at each iteration
+        recompute_.resize(size());
+        recompute_[recompute_.size()-1] = contactNetwork_.nBodyAssembler().totalHasObstacle_;
+        for (int i=mgTransfer_.size(); i>=1; i--) {
+            std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(mgTransfer_[i-1])->restrict(recompute_[i], recompute_[i-1]);
+          //  std::dynamic_pointer_cast<DenseMultigridTransfer<VectorType, BitVector, MatrixType> >(mgTransfer_[i])->restrict(recompute_[i+1], recompute_[i]);
+
+            //print(recompute_[i], "recompute: ");
+        }
+
+        for (size_t i=0; i<this->mgTransfer_.size(); i++)
+            std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(mgTransfer_[i])->setRecomputeBitField(&recompute_[i]); */
+    }
+
+   void setMatrix(const MatrixType& mat) override {
+       Base::setMatrix(mat);
+
+       ignoreHierarchy_.resize(size());
+       ignoreHierarchy_.back() = Base::ignore();
+       for (int i=mgTransfer_.size()-1; i>=0; i--) {
+           std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(mgTransfer_[i])->restrictToFathers(ignoreHierarchy_[i+1], ignoreHierarchy_[i]);
+       }
+
+        size_t maxLevel = levelPatchPreconditioners_.size()-1;
+        levelPatchPreconditioners_[maxLevel]->setMatrix(mat);
+
+        for (int i=levelPatchPreconditioners_.size()-2; i>0; i--) {
+            const auto& transfer = *mgTransfer_[i];
+
+            transfer.galerkinRestrictSetOccupation(*levelPatchPreconditioners_[i+1]->getMatrix(), levelMat_[i]);
+            transfer.galerkinRestrict(*levelPatchPreconditioners_[i+1]->getMatrix(), levelMat_[i]);
+            levelPatchPreconditioners_[i]->setMatrix(levelMat_[i]);
+
+            // Set solution vector sizes for the lower levels
+            Dune::MatrixVector::resize(levelX_[i], levelMat_[i]);
+            //Dune::MatrixVector::resize(levelRhs_[i], levelMat_[i]);
+        }
+
+        // set coarse global problem
+        const auto& coarseTransfer = *mgTransfer_[0];
+        const auto& fineMat = *levelPatchPreconditioners_[1]->getMatrix();
+        coarseTransfer.galerkinRestrictSetOccupation(fineMat, levelMat_[0]);
+        coarseTransfer.galerkinRestrict(fineMat, levelMat_[0]);
+
+        Dune::MatrixVector::resize(levelX_[0], levelMat_[0]);
+   }
+
+    void iterate() override {
+        //std::cout << "multi::iterate()" << std::endl;
+
+        size_t maxLevel = levelPatchPreconditioners_.size()-1;
+        levelX_[maxLevel] = *this->getIterate();
+        levelRhs_[maxLevel] = *Base::rhs_;
+
+        for (int i=levelPatchPreconditioners_.size()-2; i>0; i--) {
+            const auto& transfer = *mgTransfer_[i];
+
+            transfer.restrict(levelX_[i+1], levelX_[i]);
+            transfer.restrict(levelRhs_[i+1], levelRhs_[i]);
+        }
+
+        // set coarse global problem
+        const auto& coarseTransfer = *mgTransfer_[0];
+        coarseTransfer.restrict(levelX_[1], levelX_[0]);
+        coarseTransfer.restrict(levelRhs_[1], levelRhs_[0]);
+
+
+        if (mode_ == additive)
+            iterateAdd();
+        else
+            iterateMult();
+    }
+
+    void iterateAdd() {
+        *(this->x_) = 0;	
+
+        VectorType x;
+
+        // solve coarse global problem
+        /*LocalProblem<MatrixType, VectorType> localProblem(levelMat_[0], levelRhs_[0], ignoreHierarchy_[0]);
+        VectorType newR;
+        localProblem.getLocalRhs(levelX_[0], newR);
+
+        coarseSolver_.setProblem(localProblem.getMat(), levelX_[0], newR);
+        coarseSolver_.preprocess();
+        coarseSolver_.solve(); */
+
+        auto& step = levelSolvers_[0]->getIterationStep();
+        dynamic_cast<LinearIterationStep<MatrixType, VectorType>&>(step).setProblem(levelMat_[0], levelX_[0], levelRhs_[0]);
+        step.setIgnore(ignoreHierarchy_[0]);
+
+        levelSolvers_[0]->check();
+        levelSolvers_[0]->preprocess();
+        levelSolvers_[0]->solve();
+
+        mgTransfer_[0]->prolong(levelX_[0], x);
+
+        for (size_t i=1; i<levelPatchPreconditioners_.size()-1; i++) {
+            const auto& transfer = *mgTransfer_[i];
+            auto& preconditioner = *levelPatchPreconditioners_[i];
+
+            preconditioner.setIgnore(ignoreHierarchy_[i]);
+            preconditioner.apply(levelX_[i], levelRhs_[i]);
+
+            x += preconditioner.getSol();
+
+            VectorType newX;
+            transfer.prolong(x, newX);
+
+            x = newX;
+        }
+
+        auto& preconditioner = *levelPatchPreconditioners_.back();
+
+        preconditioner.setIgnore(ignoreHierarchy_.back());
+        preconditioner.apply(levelX_.back(), levelRhs_.back());
+
+        x += preconditioner.getSol();
+
+        *(this->x_) = x;
+    }
+
+
+    void iterateMult() {
+        *(this->x_) = 0;
+
+        DUNE_THROW(Dune::Exception, "Not implemented!");
+
+        /*
+        VectorType x;
+
+        for (int i=(levelPatchPreconditioners_.size()-1); i>=0; i--) {
+            VectorType updatedRhs(*(this->rhs_));
+            this->mat_->mmv(*(this->x_), updatedRhs);
+
+            mgTransfer_[i]->restrict(*(this->x_), levelX_[i]);
+            mgTransfer_[i]->restrict(updatedRhs, levelRhs_[i]);
+
+            levelPatchPreconditioners_[i]->setProblem(*(this->mat_), levelX_[i], levelRhs_[i]);
+
+            levelPatchPreconditioners_[i]->iterate();
+            const VectorType& it = levelPatchPreconditioners_[i]->getSol();
+
+            mgTransfer_[i]->prolong(it, x);
+
+            *(this->x_) += x;
+        }
+
+        VectorType updatedCoarseRhs(*(this->rhs_));
+        this->mat_->mmv(*(this->x_), updatedCoarseRhs);
+
+        size_t j = levelPatchPreconditioners_.size();
+        mgTransfer_[j]->restrict(*(this->x_), levelX_[j]);
+        mgTransfer_[j]->restrict(updatedCoarseRhs, levelRhs_[j]);
+
+        coarseGlobalPreconditioner_->setProblem(*(this->mat_), levelX_[j], levelRhs_[j]);
+        coarseGlobalPreconditioner_->iterate();
+        const VectorType& it = coarseGlobalPreconditioner_->getSol();
+
+        mgTransfer_[j]->prolong(it, x);
+
+        *(this->x_) += x; */
+    }
+
+    void build() {
+        for (size_t i=1; i<levelPatchPreconditioners_.size(); i++) {
+            levelPatchPreconditioners_[i]->build();
+        }
+    }
+
+    void setPatchDepth(const size_t patchDepth = 0) {
+        for (size_t i=1; i<levelPatchPreconditioners_.size(); i++) {
+            levelPatchPreconditioners_[i]->setPatchDepth(patchDepth);
+        }
+    }
+
+    auto contactNetwork() const {
+        return contactNetwork_;
+    }
+
+    size_t size() const {
+        return levelPatchPreconditioners_.size();
+    }
+};
+
+#endif
+

dune/tectonic/spatial-solving/preconditioners/patchproblem.hh

+#ifndef SRC_SPATIAL_SOLVING_PRECONDITIONERS_PATCH_PROBLEM_HH
+#define SRC_SPATIAL_SOLVING_PRECONDITIONERS_PATCH_PROBLEM_HH
+
+#include <math.h>   
+#include <dune/common/fmatrix.hh>
+#include <dune/common/function.hh>
+#include <dune/common/timer.hh>
+
+#include <dune/istl/matrixindexset.hh>
+//#include <dune/istl/superlu.hh>
+#include <dune/istl/umfpack.hh>
+
+#include <dune/fufem/assemblers/localoperatorassembler.hh>
+
+#include "../../utils/debugutils.hh"
+
+template <class MatrixType, class DomainType, class RangeType = DomainType>
+class PatchProblem {
+  
+private:    
+    const static size_t dim = DomainType::block_type::dimension;
+
+    using BitVector = Dune::BitSetVector<dim>;
+
+    const MatrixType& mat_;
+
+    std::vector<size_t> localToGlobal_;
+
+    MatrixType localMat_;
+    RangeType localRhs_;
+
+
+public:
+    PatchProblem(const MatrixType& mat, const Dune::BitSetVector<1>& patch) :
+        mat_(mat) {
+	  
+        // construct localToGlobal map
+        localToGlobal_.clear();
+        for (size_t i=0; i<patch.size(); ++i) {
+            if (!patch[i][0]) {
+                localToGlobal_.push_back(i);
+            }
+        }
+
+        // build local matrix
+        auto localDim = localToGlobal_.size();
+        Dune::MatrixIndexSet localIdxSet(localDim, localDim);
+
+        for(size_t rowIdx=0; rowIdx<localDim; rowIdx++) {
+            const auto globalRowIdx = localToGlobal_[rowIdx];
+            const auto& row = mat_[globalRowIdx];
+
+            const auto cEndIt = row.end();
+            for(auto cIt=row.begin(); cIt!=cEndIt; ++cIt) {
+                const auto globalColIdx = cIt.index();
+
+                auto localColIdx = std::find(localToGlobal_.begin(), localToGlobal_.end(), globalColIdx);
+                if (localColIdx!=localToGlobal_.end()) {
+                    localIdxSet.add(rowIdx, localColIdx-localToGlobal_.begin());
+                }
+            }
+        }
+
+        localIdxSet.exportIdx(localMat_);
+
+        for(size_t rowIdx=0; rowIdx<localMat_.N(); rowIdx++) {
+            auto& row = localMat_[rowIdx];
+            const auto& globalRow = mat_[localToGlobal_[rowIdx]];
+
+
+            const auto cEndIt = row.end();
+            for(auto cIt=row.begin(); cIt!=cEndIt; ++cIt) {
+                row[cIt.index()] = globalRow[localToGlobal_[cIt.index()]];
+            }
+        }
+
+        // init local rhs
+        localRhs_.resize(localDim);
+        localRhs_ = 0;
+    }
+
+    const MatrixType& mat() {
+        return localMat_;
+    }
+    
+    const RangeType& rhs() {
+        return localRhs_;
+    }
+
+    void setRhs(const RangeType& rhs){
+        for (size_t i=0; i<localRhs_.size(); i++) {
+            localRhs_[i] = rhs[localToGlobal_[i]];
+        }
+    }
+
+    void prolong(const DomainType& x, DomainType& res){
+        res.resize(mat_.N());
+        res = 0;
+
+        for (size_t i=0; i<x.size(); i++) {
+                res[localToGlobal_[i]] = x[i];
+        }
+    }
+
+    void restrict(const RangeType& x, RangeType& res){
+        res.resize(localToGlobal_.size());
+        res = 0;
+
+        for (size_t i=0; i<res.size(); i++) {
+            res[i] = x[localToGlobal_[i]];
+        }
+    }
+};
+
+#endif

dune/tectonic/spatial-solving/preconditioners/supportpatchfactory.hh

+#ifndef SUPPORT_PATCH_FACTORY_HH
+#define SUPPORT_PATCH_FACTORY_HH
+
+#include <set>
+#include <queue>
+#include <memory>
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/hybridutilities.hh>
+
+#include <dune/fufem/referenceelementhelper.hh>
+#include <dune/grid/common/mcmgmapper.hh>
+
+#include <dune/localfunctions/lagrange/pqkfactory.hh>
+
+//#include "../../data-structures/levelcontactnetwork.hh"
+#include "hierarchicleveliterator.hh"
+
+#include "../../utils/debugutils.hh"
+
+template <class LevelContactNetwork>
+class NetworkIndexMapper {
+private:
+    static const int dim = LevelContactNetwork::dim; //TODO
+    using ctype = typename LevelContactNetwork::ctype;
+
+    using FiniteElementCache = Dune::PQkLocalFiniteElementCache<ctype, ctype, dim, 1>; // P1 finite element cache
+    using Element = typename LevelContactNetwork::GridView::template Codim<0>::Entity;
+
+    const FiniteElementCache cache_;
+
+    const LevelContactNetwork& levelContactNetwork_;
+    const size_t nBodies_;
+
+    std::vector<size_t> vertexOffSet_;
+    std::vector<size_t> faceOffSet_;
+    std::vector<size_t> elementOffSet_;
+
+public:
+    NetworkIndexMapper(const LevelContactNetwork& levelContactNetwork) :
+        levelContactNetwork_(levelContactNetwork),
+        nBodies_(levelContactNetwork_.nBodies()),
+        vertexOffSet_(nBodies_),
+        faceOffSet_(nBodies_),
+        elementOffSet_(nBodies_) {
+
+        size_t vertexOffSet = 0;
+        size_t faceOffSet = 0;
+        size_t elementOffSet = 0;
+
+        for (size_t i=0; i<nBodies_; i++) {
+            const auto& gridView = levelContactNetwork_.body(i)->gridView();
+
+            vertexOffSet_[i] = vertexOffSet;
+            vertexOffSet = vertexOffSet_[i] + gridView.size(dim);
+
+            faceOffSet_[i] = faceOffSet;
+            faceOffSet = faceOffSet_[i] + gridView.size(1);
+
+            elementOffSet_[i] = elementOffSet;
+            elementOffSet = elementOffSet_[i] + gridView.size(0);
+        }
+    }
+
+    size_t vertexIndex(const size_t bodyID, const Element& elem, const size_t localVertex) const {
+        const auto& gridView = levelContactNetwork_.body(bodyID)->gridView();
+
+        auto localIdx = cache().get(elem.type()).localCoefficients().localKey(localVertex).subEntity();
+        return vertexOffSet_[bodyID] + gridView.indexSet().subIndex(elem, localIdx, dim);
+    }
+
+    template <class Intersection>
+    size_t faceIndex(const size_t bodyID, const Intersection& is) const {
+        const auto& gridView = levelContactNetwork_.body(bodyID)->gridView();
+        return faceOffSet_[bodyID] + gridView.indexSet().subIndex(is.inside(), is.indexInInside(), 1);
+    }
+
+    size_t faceIndex(const size_t bodyID, const Element& elem, const size_t indexInInside) const {
+        const auto& gridView = levelContactNetwork_.body(bodyID)->gridView();
+        return faceOffSet_[bodyID] + gridView.indexSet().subIndex(elem, indexInInside, 1);
+    }
+
+    size_t elementIndex(const size_t bodyID, const Element& elem) const {
+        const auto& gridView = levelContactNetwork_.body(bodyID)->gridView();
+        return elementOffSet_[bodyID] + gridView.indexSet().index(elem);
+    }
+
+    const auto& cache() const {
+        return cache_;
+    }
+};
+
+
+
+/**
+ *  constructs BitSetVector whose entries are
+ *      true,   if vertex is outside of patch or part of patch boundary
+ *      false,  if vertex is in interior of patch
+ */
+template <class LevelContactNetwork>
+class SupportPatchFactory
+{
+public:
+    static const int dim = LevelContactNetwork::dim; //TODO
+    using ctype = typename LevelContactNetwork::ctype;
+
+    using Patch = Dune::BitSetVector<1>;
+
+private:
+    using Element = typename LevelContactNetwork::GridView::template Codim<0>::Entity;
+
+    struct RemoteIntersection {
+        size_t bodyID; // used to store bodyID of outside elem
+        size_t contactCouplingID;
+        size_t intersectionID;
+
+        bool flip = false;
+
+        void set(const size_t _bodyID, const size_t _contactCouplingID, const size_t _intersectionID, const bool _flip = false) {
+            bodyID = _bodyID;
+            contactCouplingID = _contactCouplingID;
+            intersectionID = _intersectionID;
+            flip = _flip;
+        }
+
+        auto get(const LevelContactNetwork& levelContactNetwork) const {
+            return std::move(levelContactNetwork.glue(contactCouplingID)->getIntersection(intersectionID));
+        }
+    };
+
+    struct BodyElement {
+        size_t bodyID;
+        Element element;
+
+        void set(const size_t _bodyID, const Element& _elem) {
+            bodyID = _bodyID;
+            element = _elem;
+        }
+
+        void set(const LevelContactNetwork& levelContactNetwork, const RemoteIntersection& rIs) {
+            const auto& is = rIs.get(levelContactNetwork);
+
+            if (rIs.flip) {
+                set(rIs.bodyID, is.inside());
+            } else {
+                set(rIs.bodyID, is.outside());
+            }
+        }
+    };
+
+    using NetworkIndexMapper = NetworkIndexMapper<LevelContactNetwork>;
+
+    const size_t nBodies_;
+
+    const LevelContactNetwork& coarseContactNetwork_;
+    const LevelContactNetwork& fineContactNetwork_;
+
+    size_t nVertices_ = 0;
+
+    std::map<size_t, std::vector<RemoteIntersection>> coarseIntersectionMapper_; // map faceID to neighbor elements
+    std::map<size_t, std::vector<RemoteIntersection>> fineIntersectionMapper_; // map faceID to neighbor elements
+    std::map<size_t, std::vector<size_t>> neighborFaceDofs_; // map faceID to neighbor element dofs
+
+    std::set<size_t> interiorContactDofs_;
+
+    const NetworkIndexMapper coarseIndices_;
+    const NetworkIndexMapper fineIndices_;
+
+public:
+    SupportPatchFactory(const LevelContactNetwork& coarseContactNetwork, const LevelContactNetwork& fineContactNetwork) :
+        nBodies_(coarseContactNetwork.nBodies()),
+        coarseContactNetwork_(coarseContactNetwork),
+        fineContactNetwork_(fineContactNetwork),
+        coarseIndices_(coarseContactNetwork_),
+        fineIndices_(fineContactNetwork_) {
+
+        assert(nBodies_ == fineContactNetwork_.nBodies());
+
+        for (size_t i=0; i<nBodies_; i++) {
+            nVertices_ += fineContactNetwork_.body(i)->nVertices();
+        }
+
+        setup(coarseContactNetwork_, coarseIndices_, coarseIntersectionMapper_);
+        setup(fineContactNetwork_, fineIndices_, fineIntersectionMapper_);
+
+        // setup neighborFaceDofs_
+        for (size_t i=0; i<coarseContactNetwork_.nCouplings(); i++) {
+                const auto& coupling = *coarseContactNetwork_.coupling(i);
+                const auto& glue = *coarseContactNetwork_.glue(i);
+
+                const size_t nmBody = coupling.gridIdx_[0];
+                const size_t mBody = coupling.gridIdx_[1];
+
+                for (const auto& rIs : intersections(glue)) {
+                    size_t nmFaceIdx = coarseIndices_.faceIndex(nmBody, rIs.inside(), rIs.indexInInside());
+                    size_t mFaceIdx = coarseIndices_.faceIndex(mBody, rIs.outside(), rIs.indexInOutside());
+
+                    std::vector<std::set<size_t>> dofs;
+                    remoteIntersectionDofs(coarseIndices_, rIs, nmBody, mBody, dofs);
+
+                    if (dofs[1].size()>0) {
+                        neighborFaceDofs_[nmFaceIdx].insert(neighborFaceDofs_[nmFaceIdx].end(), dofs[1].begin(), dofs[1].end());
+                    }
+
+                    if (dofs[0].size()>0) {
+                        neighborFaceDofs_[mFaceIdx].insert(neighborFaceDofs_[mFaceIdx].end(), dofs[0].begin(), dofs[0].end());
+                    }
+                }
+        }
+
+        for (size_t i=0; i<nBodies_; i++) {
+            //std::cout << "Coarse body" << i << ":" << std::endl;
+
+            for (const auto& e : elements(coarseContactNetwork_.body(i)->gridView())) {
+                //std::cout << "elemID: " << coarseIndices_.elementIndex(i, e) << std::endl;
+                //std::cout << "vertexIDs: ";
+                const int dimElement = Element::dimension;
+
+                /*for (int j=0; j<refElement.size(dim); j++) {
+                    std::cout << coarseIndices_.vertexIndex(i, e, j) << " ";
+                }
+                std::cout << std::endl;*/
+
+                //std::cout << "intersectionIDs: ";
+                for (const auto& is : intersections(coarseContactNetwork_.body(i)->gridView(), e)) {
+                    //std::cout << coarseIndices_.faceIndex(i, e, is.indexInInside()) << " (";
+
+                    std::set<size_t> dofs;
+                    intersectionDofs(coarseIndices_, is, i, dofs);
+                    /*for (auto& d : dofs) {
+                        std::cout << d << " ";
+                    }
+                    std::cout << "); ";*/
+
+                }
+                //std::cout << std::endl  << "--------------" << std::endl;
+            }
+        }
+
+        /*std::cout << std::endl;
+        for (auto& kv : neighborFaceDofs_) {
+            std::cout << "faceID: " << kv.first << " dofs: ";
+
+            const auto& dofs = kv.second;
+            for (auto& d : dofs) {
+                std::cout << d << " ";
+            }
+            std::cout << std::endl;
+        }*/
+
+        /*
+            // neighborElementDofs_ contains dofs of connected intersections that are neighbors to a face given by faceID,
+            // boundary dofs are contained once, interior dofs multiple times
+            // task: remove boundary dofs and redundant multiples of interior dofs
+            // TODO: only works in 2D
+            for (auto& it : neighborElementDofs) {
+                auto& dofs = it.second;
+
+                std::sort(dofs.begin(), dofs.end());
+                auto& dof = dofs.begin();
+                while (dof != dofs.end()) {
+                    auto next = dof;
+                    if (*dof != *(++next))
+                        dof = dofs.erase(dof);
+                }
+                dofs.erase(std::unique(dofs.begin(), dofs.end()), dofs.end());
+            }
+        */
+    }
+
+    void build(const size_t bodyID, const Element& coarseElement, const size_t localVertex, Patch& patchDofs, const size_t patchDepth = 0) {
+        auto isPatchIntersection = [this](auto& is, const size_t bodyID, const std::set<size_t>& patchVertices, std::set<size_t>& newPatchVertices) {
+            std::set<size_t> dofs;
+            intersectionDofs(coarseIndices_, is, bodyID, dofs);
+
+            for (auto& dof : dofs) {
+                if (patchVertices.count(dof)) {
+                    newPatchVertices.insert(dofs.begin(), dofs.end());
+                    return true;
+                }
+            }
+            return false;
+        };
+
+        auto isRPatchIntersection = [this](std::vector<std::set<size_t>>& dofs, const std::set<size_t>& patchVertices) {
+            for (auto& dof : dofs[0]) {
+                if (patchVertices.count(dof)) {
+                    return true;
+                }
+            }
+            return false;
+        };
+
+        patchDofs.resize(nVertices_);
+        patchDofs.setAll();
+
+        std::vector<BodyElement> patchElements;
+        std::set<size_t> visitedElements;
+
+        std::set<size_t> patchVertices;
+        patchVertices.insert(coarseIndices_.vertexIndex(bodyID, coarseElement, localVertex));
+
+        std::set<size_t> visitedBodies;
+
+        BodyElement seedElement;
+        seedElement.set(bodyID, coarseElement);
+
+        std::queue<BodyElement> patchSeeds;
+        patchSeeds.push(seedElement);
+
+        for (size_t depth=0; depth<=patchDepth; depth++) {
+            std::set<size_t> newPatchVertices;
+            std::vector<BodyElement> nextSeeds;
+
+            while (!patchSeeds.empty()) {
+                const auto& patchSeed = patchSeeds.front();
+
+                size_t elemIdx = coarseIndices_.elementIndex(patchSeed.bodyID, patchSeed.element);
+
+                if (visitedElements.count(elemIdx)) {
+                    patchSeeds.pop();
+                    continue;
+                }
+
+                visitedBodies.insert(patchSeed.bodyID);
+                visitedElements.insert(elemIdx);
+                patchElements.emplace_back(patchSeed);
+
+                size_t bodyIdx = patchSeed.bodyID;
+                const auto& elem = patchSeed.element;
+                const auto& gridView = coarseContactNetwork_.body(bodyIdx)->gridView();
+
+                for (const auto& it : intersections(gridView, elem)) {
+                    if (isPatchIntersection(it, bodyIdx, patchVertices, newPatchVertices)) {
+                        if (it.neighbor()) {
+                            BodyElement neighbor;
+                            neighbor.set(bodyIdx, it.outside());
+
+                            size_t neighborIdx = coarseIndices_.elementIndex(neighbor.bodyID, neighbor.element);
+
+                            if (visitedElements.count(neighborIdx))
+                                continue;
+
+                            patchSeeds.push(neighbor);
+                        } else {
+                            size_t faceIdx = coarseIndices_.faceIndex(bodyIdx, it);
+
+                            if (coarseIntersectionMapper_.count(faceIdx)) {
+                                const auto& intersections = coarseIntersectionMapper_[faceIdx];
+                                for (size_t i=0; i<intersections.size(); i++) {
+                                    const auto& intersection = intersections[i];
+                                    const auto& rIs = intersection.get(coarseContactNetwork_);
+
+                                    BodyElement neighbor;
+                                    neighbor.set(coarseContactNetwork_, intersection);
+
+                                    size_t neighborBody = neighbor.bodyID;
+                                    size_t neighborIdx = coarseIndices_.elementIndex(neighborBody, neighbor.element);
+
+                                    if (visitedElements.count(neighborIdx) || visitedBodies.count(neighborBody))
+                                        continue;
+
+                                    std::vector<std::set<size_t>> dofs;
+                                    remoteIntersectionDofs(coarseIndices_, rIs, bodyIdx, neighborBody, dofs, intersection.flip);
+
+                                    if (isRPatchIntersection(dofs, patchVertices)) {
+                                        patchVertices.insert(dofs[1].begin(), dofs[1].end());
+                                        patchSeeds.push(neighbor);
+                                    } else {
+                                        newPatchVertices.insert(dofs[1].begin(), dofs[1].end());
+                                        nextSeeds.emplace_back(neighbor);
+                                    }
+                                }
+                            }
+                        }
+                    } else {
+                        if (it.neighbor()) {
+                            BodyElement neighbor;
+                            neighbor.set(bodyIdx, it.outside());
+
+                            size_t neighborIdx = coarseIndices_.elementIndex(bodyIdx, neighbor.element);
+
+                            if (visitedElements.count(neighborIdx))
+                                continue;
+
+                            nextSeeds.emplace_back(neighbor);
+                        } else {
+                            size_t faceIdx = coarseIndices_.faceIndex(bodyIdx, it);
+
+                            if (coarseIntersectionMapper_.count(faceIdx)) {
+                                const auto& intersections = coarseIntersectionMapper_[faceIdx];
+                                for (size_t i=0; i<intersections.size(); i++) {
+                                    BodyElement neighbor;
+                                    neighbor.set(coarseContactNetwork_, intersections[i]);
+
+                                    size_t neighborBody = neighbor.bodyID;
+                                    size_t neighborIdx = coarseIndices_.elementIndex(neighbor.bodyID, neighbor.element);
+
+                                    if (visitedElements.count(neighborIdx) || visitedBodies.count(neighborBody))
+                                        continue;
+
+                                    newPatchVertices.insert(neighborFaceDofs_[faceIdx].begin(), neighborFaceDofs_[faceIdx].end());
+                                    nextSeeds.emplace_back(neighbor);
+                                }
+                            }
+                        }
+                    }
+                }
+                patchSeeds.pop();
+            }
+
+            for (size_t j=0; j<nextSeeds.size(); j++) {
+                patchSeeds.push(nextSeeds[j]);
+            }
+            patchVertices.insert(newPatchVertices.begin(), newPatchVertices.end());
+        }
+
+        // construct fine patch
+        using HierarchicLevelIteratorType = HierarchicLevelIterator<typename LevelContactNetwork::GridType>;
+        std::set<size_t> patchBoundary;
+        for (size_t i=0; i<patchElements.size(); ++i) {
+            const auto& coarseElem = patchElements[i];
+
+            const auto& grid = coarseContactNetwork_.body(coarseElem.bodyID)->gridView().grid();
+            const auto fineLevel = fineContactNetwork_.body(coarseElem.bodyID)->level();
+
+            // add fine patch elements
+            HierarchicLevelIteratorType endHierIt(grid, coarseElem.element, HierarchicLevelIteratorType::PositionFlag::end, fineLevel);
+            HierarchicLevelIteratorType hierIt(grid, coarseElem.element, HierarchicLevelIteratorType::PositionFlag::begin, fineLevel);
+            for (; hierIt!=endHierIt; ++hierIt) {
+                const Element& fineElem = *hierIt;
+                addLocalDofs(coarseElem, fineElem, visitedElements, patchBoundary, patchDofs);
+            }
+        }
+    }
+
+    const auto& coarseContactNetwork() {
+        return coarseContactNetwork_;
+    }
+
+    const auto& coarseIndices() const {
+        return coarseIndices_;
+    }
+
+private:
+    void setup(const LevelContactNetwork& levelContactNetwork, const NetworkIndexMapper& indices, std::map<size_t, std::vector<RemoteIntersection>>& faceToRIntersections) {
+        for (size_t i=0; i<levelContactNetwork.nCouplings(); i++) {
+            const auto& coupling = *levelContactNetwork.coupling(i);
+            const auto& glue = *levelContactNetwork.glue(i);
+
+            const size_t nmBody = coupling.gridIdx_[0];
+            const size_t mBody = coupling.gridIdx_[1];
+
+            size_t rIsID = 0;
+            for (const auto& rIs : intersections(glue)) {
+                size_t nmFaceIdx = indices.faceIndex(nmBody, rIs.inside(), rIs.indexInInside());
+
+                RemoteIntersection nmIntersection;
+                nmIntersection.set(mBody, i, rIsID, false);
+                faceToRIntersections[nmFaceIdx].emplace_back(nmIntersection);
+
+                if (rIs.neighbor()) {
+                    size_t mFaceIdx = indices.faceIndex(mBody, rIs.outside(), rIs.indexInOutside());
+
+                    RemoteIntersection mIntersection;
+                    mIntersection.set(nmBody, i, rIsID, true);
+                    faceToRIntersections[mFaceIdx].emplace_back(mIntersection);
+                }
+
+                rIsID++;
+            }
+        }
+    }
+
+    template <class Intersection>
+    bool containsInsideSubentity(const Element& elem, const Intersection& intersection, int subEntity, int codim) {
+	    return ReferenceElementHelper<double, dim>::subEntityContainsSubEntity(elem.type(), intersection.indexInInside(), 1, subEntity, codim);
+	}
+      
+    void addLocalDofs(const BodyElement& coarseElem, const Element& fineElem, const std::set<size_t>& coarsePatch, std::set<size_t>& patchBoundary, Patch& patchDofs) {
+	    // insert local dofs
+        const auto& gridView = fineContactNetwork_.body(coarseElem.bodyID)->gridView();
+        const size_t bodyID = coarseElem.bodyID;
+
+        const auto coarseLevel = coarseContactNetwork_.body(coarseElem.bodyID)->level();
+
+	    // search for parts of the global and inner boundary 
+        for (const auto& is : intersections(gridView, fineElem)) {
+            if (is.neighbor()) {
+                std::set<size_t> isDofs;
+                intersectionDofs(fineIndices_, is, bodyID, isDofs);
+
+                auto outsideFather = coarseFather(is.outside(), coarseLevel);
+		    
+                if (!coarsePatch.count(coarseIndices_.elementIndex(coarseElem.bodyID, outsideFather))) {
+                    addToPatchBoundary(isDofs, patchBoundary, patchDofs);
+                } else {
+                    addToPatch(isDofs, patchBoundary, patchDofs);
+                }
+            } else {
+                size_t faceIdx = fineIndices_.faceIndex(bodyID, is);
+                if (fineIntersectionMapper_.count(faceIdx)) {
+                    const auto& intersections = fineIntersectionMapper_[faceIdx];
+
+                    for (size_t i=0; i<intersections.size(); i++) {
+                        const auto& intersection = intersections[i];
+                        const auto& rIs = intersection.get(fineContactNetwork_);
+                        const size_t outsideID = intersection.bodyID;
+
+                        std::vector<std::set<size_t>> rIsDofs;
+                        remoteIntersectionDofs(fineIndices_, rIs, bodyID, outsideID, rIsDofs, intersection.flip);
+
+                        if (rIs.neighbor()) {
+                            Element outsideFather;
+                            const size_t outsideLevel = coarseContactNetwork_.body(outsideID)->level();
+
+                            if (intersection.flip)
+                                outsideFather = coarseFather(rIs.inside(), outsideLevel);
+                            else
+                                outsideFather = coarseFather(rIs.outside(), outsideLevel);
+
+                            std::set<size_t> totalRIsDofs(rIsDofs[0]);
+                            totalRIsDofs.insert(rIsDofs[1].begin(), rIsDofs[1].end());
+
+                            if (!coarsePatch.count(coarseIndices_.elementIndex(outsideID, outsideFather))) {
+                                addToPatchBoundary(totalRIsDofs, patchBoundary, patchDofs);
+                            } else {
+                                addToPatch(totalRIsDofs, patchBoundary, patchDofs);
+                            }
+                        } else {
+                            addToPatchBoundary(rIsDofs[0], patchBoundary, patchDofs);
+                        }
+                    }
+                } else {
+                    std::set<size_t> isDofs;
+                    intersectionDofs(fineIndices_, is, bodyID, isDofs);
+
+                    addToPatch(isDofs, patchBoundary, patchDofs);
+                }
+            }
+        }
+	}
+
+    auto coarseFather(const Element& fineElem, const int coarseLevel) const {
+        Element coarseElem = fineElem;
+        while (coarseElem.level() > coarseLevel) {
+            coarseElem = coarseElem.father();
+        }
+        return coarseElem;
+    }
+
+    void addToPatchBoundary(const std::set<size_t>& dofs, std::set<size_t>& patchBoundary, Patch& patchDofs) {
+        for (auto& dof : dofs) {
+            patchBoundary.insert(dof);
+            patchDofs[dof] = true;
+        }
+    }
+
+    void addToPatch(const std::set<size_t>& dofs, const std::set<size_t>& patchBoundary, Patch& patchDofs) {
+        for (auto& dof : dofs) {
+            if (!patchBoundary.count(dof)) {
+                patchDofs[dof] = false;
+            }
+        }
+    }
+
+    template <class RIntersection>
+    void remoteIntersectionDofs(const NetworkIndexMapper& indices, const RIntersection& is, const size_t insideID, const size_t outsideID, std::vector<std::set<size_t>>& dofs, bool flip = false) {
+        assert(!flip || (flip && is.neighbor()));
+
+        dofs.resize(2);
+        dofs[0].clear();
+        dofs[1].clear();
+
+        const auto& isGeo = is.geometry();
+        const auto& isRefElement = Dune::ReferenceElements<ctype, dim-1>::general(is.type());
+
+        if (!flip) {
+             const auto& inside = is.inside();
+             const auto& insideGeo = inside.geometry();
+             const auto& insideRefElement = Dune::ReferenceElements<ctype, dim>::general(inside.type());
+             const int dimElement = Dune::ReferenceElements<ctype, dim>::dimension;
+
+             for (int i=0; i<insideRefElement.size(is.indexInInside(), 1, dimElement); i++) {
+                 size_t idxInElement = insideRefElement.subEntity(is.indexInInside(), 1, i, dimElement);
+
+                 const auto& localCorner = isGeo.local(insideGeo.corner(idxInElement));
+
+                 if (isRefElement.checkInside(localCorner)) {
+                     dofs[0].insert(indices.vertexIndex(insideID, inside, idxInElement));
+                 }
+             }
+
+             if (is.neighbor()) {
+                 const auto& outside = is.outside();
+                 const auto& outsideGeo = outside.geometry();
+                 const auto& outsideRefElement = Dune::ReferenceElements<ctype, dim>::general(outside.type());
+
+                 for (int i=0; i<outsideRefElement.size(is.indexInOutside(), 1, dimElement); i++) {
+                     size_t idxInElement = outsideRefElement.subEntity(is.indexInOutside(), 1, i, dimElement);
+
+                     const auto& localCorner = isGeo.local(outsideGeo.corner(idxInElement));
+
+                     if (isRefElement.checkInside(localCorner)) {
+                         dofs[1].insert(indices.vertexIndex(outsideID, outside, idxInElement));
+                     }
+                 }
+             }
+        } else {
+            const auto& inside = is.outside();
+            const auto& insideGeo = inside.geometry();
+            const auto& insideRefElement = Dune::ReferenceElements<ctype, dim>::general(inside.type());
+            const int dimElement = Dune::ReferenceElements<ctype, dim>::dimension;
+
+            for (int i=0; i<insideRefElement.size(is.indexInOutside(), 1, dimElement); i++) {
+                size_t idxInElement = insideRefElement.subEntity(is.indexInOutside(), 1, i, dimElement);
+
+                const auto& localCorner = isGeo.local(insideGeo.corner(idxInElement));
+
+                if (isRefElement.checkInside(localCorner)) {
+                    dofs[0].insert(indices.vertexIndex(insideID, inside, idxInElement));
+                }
+            }
+
+            const auto& outside = is.inside();
+            const auto& outsideGeo = outside.geometry();
+            const auto& outsideRefElement = Dune::ReferenceElements<ctype, dim>::general(outside.type());
+
+            for (int i=0; i<outsideRefElement.size(is.indexInInside(), 1, dimElement); i++) {
+                size_t idxInElement = outsideRefElement.subEntity(is.indexInInside(), 1, i, dimElement);
+
+                const auto& localCorner = isGeo.local(outsideGeo.corner(idxInElement));
+
+                if (isRefElement.checkInside(localCorner)) {
+                    dofs[1].insert(indices.vertexIndex(outsideID, outside, idxInElement));
+                }
+            }
+        }
+
+    }
+
+    template <class Intersection>
+    void intersectionDofs(const NetworkIndexMapper& indices, const Intersection& is, const size_t bodyID, std::set<size_t>& dofs) {
+        dofs.clear();
+
+        const Element& elem = is.inside();
+        size_t intersectionIndex = is.indexInInside();
+
+        const int dimElement = Element::dimension;
+        const auto& refElement = Dune::ReferenceElements<double, dimElement>::general(elem.type());
+
+        for (int i=0; i<refElement.size(intersectionIndex, 1, dimElement); i++) {
+            size_t idxInElement = refElement.subEntity(intersectionIndex, 1, i, dimElement);
+            dofs.insert(indices.vertexIndex(bodyID, elem, idxInElement));
+        }
+    }
+};
+#endif

dune/tectonic/spatial-solving/solverfactory.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/solvers/common/wrapownshare.hh>
+#include <dune/solvers/iterationsteps/blockgssteps.hh>
+#include <dune/solvers/solvers/umfpacksolver.hh>
+
+#include "solverfactory.hh"
+
+#include "../utils/debugutils.hh"
+
+template <class Functional, class BitVector>
+SolverFactory<Functional, BitVector>::SolverFactory(
+    const Dune::ParameterTree& parset,
+    Functional& J,
+    const BitVector& ignoreNodes) :
+        parset_(parset),
+        J_(Dune::Solvers::wrap_own_share<const Functional>(std::forward<Functional>(J))),
+        ignoreNodes_(ignoreNodes)
+{}
+
+template <class Functional, class BitVector>
+SolverFactory<Functional, BitVector>::SolverFactory(
+    const Dune::ParameterTree& parset,
+    std::shared_ptr<Functional> J,
+    const BitVector& ignoreNodes) :
+        parset_(parset),
+        J_(Dune::Solvers::wrap_own_share<const Functional>(J)),
+        ignoreNodes_(ignoreNodes)
+{}
+
+template <class Functional, class BitVector>
+template <class LinearSolver>
+void SolverFactory<Functional, BitVector>::build(std::shared_ptr<LinearSolver>& linearSolver) {
+    nonlinearSmoother_ = std::make_shared<NonlinearSmoother>(*J_, dummyIterate_, LocalSolver());
+
+    tnnmgStep_ = std::make_shared<Step>(*J_, dummyIterate_, nonlinearSmoother_, linearSolver, DefectProjection(), LineSearchSolver());
+    tnnmgStep_->setPreSmoothingSteps(parset_.get<int>("main.pre"));
+    tnnmgStep_->setIgnore(ignoreNodes_);
+}
+
+template <class Functional, class BitVector>
+void SolverFactory<Functional, BitVector>::setProblem(Vector& x) {
+    nonlinearSmoother_->setProblem(x);
+    tnnmgStep_->setProblem(x);
+}
+
+
+template <class Functional, class BitVector>
+auto SolverFactory<Functional, BitVector>::step()
+-> std::shared_ptr<Step> {
+    return tnnmgStep_;
+}
+
+#include "solverfactory_ex.cc"

dune/tectonic/spatial-solving/solverfactory.hh

+#ifndef SRC_SPATIAL_SOLVING_SOLVERFACTORY_HH
+#define SRC_SPATIAL_SOLVING_SOLVERFACTORY_HH
+
+//#define NEW_TNNMG_COMPUTE_ITERATES_DIRECTLY
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/parametertree.hh>
+
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+#include <dune/solvers/iterationsteps/cgstep.hh>
+
+#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
+#include <dune/tnnmg/iterationsteps/nonlineargsstep.hh>
+#include <dune/tnnmg/projections/obstacledefectprojection.hh>
+#include <dune/tnnmg/localsolvers/scalarobstaclesolver.hh>
+
+#include "tnnmg/linearization.hh"
+#include "tnnmg/linesearchsolver.hh"
+#include "tnnmg/localbisectionsolver.hh"
+
+template <class FunctionalTEMPLATE, class BitVectorType>
+class SolverFactory {
+public:
+    using Functional = FunctionalTEMPLATE;
+    using Matrix = typename Functional::Matrix;
+    using Vector = typename Functional::Vector;
+    using BitVector = BitVectorType;
+
+    using LocalSolver = LocalBisectionSolver;
+    using NonlinearSmoother = Dune::TNNMG::NonlinearGSStep<Functional, LocalBisectionSolver, BitVector>;
+    using Linearization = Linearization<Functional, BitVector>;
+    using DefectProjection = typename Dune::TNNMG::ObstacleDefectProjection;
+
+    using Step = Dune::TNNMG::TNNMGStep<Functional, BitVector, Linearization, DefectProjection, LineSearchSolver>;
+
+  SolverFactory(const Dune::ParameterTree&,
+                Functional&,
+                const BitVector&);
+
+  SolverFactory(const Dune::ParameterTree&,
+                std::shared_ptr<Functional>,
+                const BitVector&);
+
+  template <class LinearSolver>
+  void build(std::shared_ptr<LinearSolver>& linearSolver);
+
+  void setProblem(Vector& x);
+
+  auto step() -> std::shared_ptr<Step>;
+
+private:
+  const Dune::ParameterTree& parset_;
+
+  Vector dummyIterate_;
+  std::shared_ptr<const Functional> J_;
+  const BitVector& ignoreNodes_;
+
+  // nonlinear smoother
+  std::shared_ptr<NonlinearSmoother> nonlinearSmoother_;
+
+  // TNNMGStep
+  std::shared_ptr<Step> tnnmgStep_;
+};
+
+#endif

dune/tectonic/spatial-solving/solverfactory_ex.cc

+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include <dune/solvers/solvers/loopsolver.hh>
+
+#include "solverfactory_tmpl.cc"
+
+using MyLinearSolver = Dune::Solvers::LoopSolver<Vector, BitVector>;
+
+template class SolverFactory<MyFunctional, BitVector>;
+template void SolverFactory<MyFunctional, BitVector>::build<MyLinearSolver>(std::shared_ptr<MyLinearSolver>&);
+
+template class SolverFactory<MyZeroFunctional, BitVector>;
+template void SolverFactory<MyZeroFunctional, BitVector>::build<MyLinearSolver>(std::shared_ptr<MyLinearSolver>&);
+
+

dune/tectonic/spatial-solving/solverfactory_tmpl.cc

+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include "../explicitgrid.hh"
+#include "../explicitvectors.hh"
+
+#include "../data-structures/friction/globalfriction.hh"
+#include "tnnmg/functional.hh"
+#include "tnnmg/zerononlinearity.hh"
+
+#include "solverfactory.hh"
+
+using MyLinearSolver = Dune::Solvers::LoopSolver<Vector, BitVector>;
+using MyGlobalFriction = GlobalFriction<Matrix, Vector>;
+
+using MyFunctional = Functional<Matrix&, Vector&, MyGlobalFriction&, Vector&, Vector&, double>;
+using MySolverFactory = SolverFactory<MyFunctional, BitVector>;
+
+using MyZeroFunctional = Functional<Matrix&, Vector&, ZeroNonlinearity&, Vector&, Vector&, double>;
+using MyZeroSolverFactory = SolverFactory<MyZeroFunctional, BitVector>;

dune/tectonic/spatial-solving/tnnmg/CMakeLists.txt

+add_custom_target(tectonic_dune_spatial-solving_tnnmg SOURCES
+  functional.hh
+  linearcorrection.hh
+  linearization.hh
+  linesearchsolver.hh
+  localbisectionsolver.hh
+  zerononlinearity.hh
+)
+
+#install headers
+install(FILES
+  functional.hh
+  linearcorrection.hh
+  linearization.hh
+  linesearchsolver.hh
+  localbisectionsolver.hh
+  zerononlinearity.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)

dune/tectonic/spatial-solving/tnnmg/functional.hh

+// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set ts=8 sw=2 et sts=2:
+#ifndef DUNE_TECTONIC_SPATIAL_SOLVING_FUNCTIONAL_HH
+#define DUNE_TECTONIC_SPATIAL_SOLVING_FUNCTIONAL_HH
+
+#include <cstddef>
+#include <type_traits>
+
+#include <dune/solvers/common/wrapownshare.hh>
+#include <dune/solvers/common/copyorreference.hh>
+#include <dune/solvers/common/interval.hh>
+#include <dune/solvers/common/resize.hh>
+
+#include <dune/matrix-vector/algorithm.hh>
+#include <dune/matrix-vector/axpy.hh>
+
+#include <dune/tnnmg/functionals/boxconstrainedquadraticfunctional.hh>
+
+#include "../../utils/maxeigenvalue.hh"
+#include "../../utils/debugutils.hh"
+
+template<class M, class V, class N, class L, class U, class R>
+class DirectionalRestriction;
+
+template<class M, class E, class V, class N, class L, class U, class O, class R>
+class ShiftedFunctional;
+
+template <class M, class V, class N, class L, class U, class R>
+class Functional;
+
+/** \brief Coordinate restriction of box constrained quadratic functional with nonlinearity;
+ *         mainly used for presmoothing in TNNMG algorithm
+ *
+ *  \tparam M global matrix type
+ *  \tparam V global vector type
+ *  \tparam V nonlinearity type
+ *  \tparam R Range type
+ */
+template<class M, class E, class V, class N, class L, class U, class O, class R>
+class ShiftedFunctional {
+
+public:
+  using Matrix = std::decay_t<M>;
+  using Eigenvalues = E;
+  using Vector = std::decay_t<V>;
+  using Nonlinearity = std::decay_t<N>;
+  using LowerObstacle = std::decay_t<L>;
+  using UpperObstacle = std::decay_t<U>;
+  using Origin = std::decay_t<O>;
+  using Range = R;
+
+  template <class MM, class VV, class LL, class UU, class OO>
+  ShiftedFunctional(MM&& matrix, const Eigenvalues& maxEigenvalues, VV&& linearPart, const Nonlinearity& phi,
+                    LL&& lower, UU&& upper, OO&& origin) :
+      quadraticPart_(std::forward<MM>(matrix)),
+      maxEigenvalues_(maxEigenvalues),
+      originalLinearPart_(std::forward<VV>(linearPart)),
+      phi_(phi),
+      originalLowerObstacle_(std::forward<LL>(lower)),
+      originalUpperObstacle_(std::forward<UU>(upper)),
+      origin_(std::forward<OO>(origin))
+  {}
+
+  Range operator()(const Vector& v) const
+  {
+    auto w = origin();
+    w += v;
+    if (Dune::TNNMG::checkInsideBox(w, originalLowerObstacle(), originalUpperObstacle())) {
+        Vector temp;
+        Dune::Solvers::resizeInitializeZero(temp,v);
+        quadraticPart().umv(w, temp);
+        temp *= 0.5;
+        temp -= originalLinearPart();
+        return temp*w + phi()(w);
+      }
+    return std::numeric_limits<Range>::max();
+
+    /*Vector temp;
+    Dune::Solvers::resizeInitializeZero(temp,v);
+    quadraticPart().umv(w, temp);
+    temp *= 0.5;
+    temp -= originalLinearPart();
+    return temp*w + phi()(w);*/
+  }
+
+  const Matrix& quadraticPart() const
+  {
+    return quadraticPart_.get();
+  }
+
+  const Vector& originalLinearPart() const
+  {
+    return originalLinearPart_.get();
+  }
+
+  const Origin& origin() const
+  {
+    return origin_.get();
+  }
+
+  void updateOrigin()
+  {}
+
+  void updateOrigin(std::size_t i)
+  {}
+
+  const LowerObstacle& originalLowerObstacle() const
+  {
+    return originalLowerObstacle_.get();
+  }
+
+  const UpperObstacle& originalUpperObstacle() const
+  {
+    return originalUpperObstacle_.get();
+  }
+
+  const auto& phi() const {
+      return phi_;
+  }
+
+  const auto& maxEigenvalues() const {
+    return maxEigenvalues_;
+  }
+
+protected:
+  Dune::Solvers::ConstCopyOrReference<M> quadraticPart_;
+  const Eigenvalues& maxEigenvalues_;
+
+  Dune::Solvers::ConstCopyOrReference<V> originalLinearPart_;
+  const Nonlinearity& phi_;
+
+  Dune::Solvers::ConstCopyOrReference<L> originalLowerObstacle_;
+  Dune::Solvers::ConstCopyOrReference<U> originalUpperObstacle_;
+
+  Dune::Solvers::ConstCopyOrReference<O> origin_;
+};
+
+/** \brief Coordinate restriction of box constrained quadratic functional with nonlinearity;
+ *         mainly used for line search step in TNNMG algorithm
+ *
+ *  \tparam M Global matrix type
+ *  \tparam V Global vector type
+ *  \tparam N nonlinearity type
+ *  \tparam L Global lower obstacle type
+ *  \tparam U Global upper obstacle type
+ *  \tparam R Range type
+ */
+template<class M, class V, class N, class L, class U, class R=double>
+class DirectionalRestriction :
+  public Dune::TNNMG::BoxConstrainedQuadraticDirectionalRestriction<M,V,L,U,R>
+{
+  using Base = Dune::TNNMG::BoxConstrainedQuadraticDirectionalRestriction<M,V,L,U,R>;
+  using Nonlinearity = N;
+  using Interval = typename Dune::Solvers::Interval<R>;
+
+public:
+  using GlobalMatrix = typename Base::GlobalMatrix;
+  using GlobalVector = typename Base::GlobalVector;
+  using GlobalLowerObstacle = typename Base::GlobalLowerObstacle;
+  using GlobalUpperObstacle = typename Base::GlobalUpperObstacle;
+
+  using Matrix = typename Base::Matrix;
+  using Vector = typename Base::Vector;
+
+  DirectionalRestriction(const GlobalMatrix& matrix, const GlobalVector& linearTerm, const Nonlinearity& phi, const GlobalLowerObstacle& lower, const GlobalLowerObstacle& upper, const GlobalVector& origin, const GlobalVector& direction, const double scaling) :
+    Base(matrix, linearTerm, lower, upper, origin, direction),
+    origin_(origin),
+    direction_(direction),
+    phi_(phi),
+    scaling_(scaling)
+  {
+    phi_.directionalDomain(origin_, direction_, domain_);
+
+    //std::cout << domain_[0] << " " << domain_[1] << "Phi domain:" << std::endl;
+    //std::cout << this->defectLower_ << " " << this->defectUpper_ << "defect obstacles:" << std::endl;
+
+    if (domain_[0] < this->defectLower_) {
+       domain_[0] = this->defectLower_;
+    }
+
+    if (domain_[1] > this->defectUpper_) {
+       domain_[1] = this->defectUpper_;
+    }
+  }
+
+    auto subDifferential(double x) const {
+      Interval D;
+      GlobalVector uxv = origin_;
+
+      Dune::MatrixVector::addProduct(uxv, x, direction_);
+      phi_.directionalSubDiff(uxv, direction_, D);
+      auto const Axmb = this->quadraticPart_ * x - this->linearPart_;
+      D[0] += Axmb;
+      D[1] += Axmb;
+
+      return D;
+    }
+
+    auto domain() const {
+      return domain_;
+    }
+
+    const GlobalVector& origin() const {
+        return origin_;
+    }
+
+    const GlobalVector& direction() const {
+        return direction_;
+    }
+
+    double scaling() const {
+        return scaling_;
+    }
+protected:
+  const GlobalVector& origin_;
+  GlobalVector direction_;
+
+  const Nonlinearity& phi_;
+  Interval domain_;
+
+  const double scaling_;
+};
+
+
+/** \brief Box constrained quadratic functional with nonlinearity
+ *         Note: call setIgnore() to set up functional in affine subspace with ignore information
+ *
+ *  \tparam V Vector type
+ *  \tparam N Nonlinearity type
+ *  \tparam L Lower obstacle type
+ *  \tparam U Upper obstacle type
+ *  \tparam R Range type
+ */
+template<class N, class V, class R>
+class FirstOrderFunctional {
+    using Interval = typename Dune::Solvers::Interval<R>;
+
+public:
+    using Nonlinearity = std::decay_t<N>;
+    using Vector = V;
+    using LowerObstacle = V;
+    using UpperObstacle = V;
+    using Range = R;
+
+    using field_type = typename V::field_type;
+public:
+    template <class LL, class UU, class OO, class DD>
+    FirstOrderFunctional(
+            const Range& maxEigenvalue,
+            const Range& linearPart,
+            const Nonlinearity& phi,
+            LL&& lower,
+            UU&& upper,
+            OO&& origin,
+            DD&& direction) :
+        quadraticPart_(maxEigenvalue),
+        linearPart_(linearPart),
+        lower_(std::forward<LL>(lower)),
+        upper_(std::forward<UU>(upper)),
+        origin_(std::forward<OO>(origin)),
+        direction_(std::forward<DD>(direction)),
+        phi_(phi) {
+
+        // set defect obstacles
+        defectLower_ = -std::numeric_limits<field_type>::max();
+        defectUpper_ = std::numeric_limits<field_type>::max();
+        Dune::TNNMG::directionalObstacles(origin_.get(), direction_.get(), lower_.get(), upper_.get(), defectLower_, defectUpper_);
+
+        // set domain
+        phi_.directionalDomain(origin_.get(), direction_.get(), domain_);
+
+        if (domain_[0] < defectLower_) {
+           domain_[0] = defectLower_;
+        }
+
+        if (domain_[1] > defectUpper_) {
+           domain_[1] = defectUpper_;
+        }
+    }
+
+  Range operator()(const Vector& v) const
+  {
+    DUNE_THROW(Dune::NotImplemented, "Evaluation of FirstOrderFunctional not implemented");
+  }
+
+  auto subDifferential(double x) const {
+    Interval Di;
+    Vector uxv = origin_.get();
+
+    Dune::MatrixVector::addProduct(uxv, x, direction_.get());
+    phi_.directionalSubDiff(uxv, direction_.get(), Di);
+
+    const auto Axmb = quadraticPart_ * x - linearPart_;
+    Di[0] += Axmb;
+    Di[1] += Axmb;
+
+    return Di;
+  }
+
+  const Interval& domain() const {
+    return domain_;
+  }
+
+  const auto& origin() const {
+      return origin_.get();
+  }
+
+  const auto& direction() const {
+      return direction_.get();
+  }
+
+  const auto& lowerObstacle() const {
+      return defectLower_;
+  }
+
+  const auto& upperObstacle() const {
+      return defectUpper_;
+  }
+
+  const auto& quadraticPart() const {
+      return quadraticPart_;
+  }
+
+  const auto& linearPart() const {
+      return linearPart_;
+  }
+private:
+  const Range quadraticPart_;
+  const Range linearPart_;
+
+  Dune::Solvers::ConstCopyOrReference<LowerObstacle> lower_;
+  Dune::Solvers::ConstCopyOrReference<UpperObstacle> upper_;
+
+  Dune::Solvers::ConstCopyOrReference<Vector> origin_;
+  Dune::Solvers::ConstCopyOrReference<Vector> direction_;
+
+  const Nonlinearity& phi_;
+
+  Interval domain_;
+
+  field_type defectLower_;
+  field_type defectUpper_;
+};
+
+// \ToDo This should be an inline friend of ShiftedBoxConstrainedQuadraticFunctional
+// but gcc-4.9.2 shipped with debian jessie does not accept
+// inline friends with auto return type due to bug-59766.
+// Notice, that this is fixed in gcc-4.9.3.
+template<class GlobalShiftedFunctional, class Index>
+auto coordinateRestriction(const GlobalShiftedFunctional& f, const Index& i)
+{
+  using Range = typename GlobalShiftedFunctional::Range;
+  using LocalMatrix = std::decay_t<decltype(f.quadraticPart()[i][i])>;
+  using LocalVector = std::decay_t<decltype(f.originalLinearPart()[i])>;
+  using LocalLowerObstacle = std::decay_t<decltype(f.originalLowerObstacle()[i])>;
+  using LocalUpperObstacle = std::decay_t<decltype(f.originalUpperObstacle()[i])>;
+
+  using namespace Dune::MatrixVector;
+  namespace H = Dune::Hybrid;
+
+  const LocalMatrix* Aii_p = nullptr;
+
+  //print(f.originalLinearPart(), "f.linearPart: ");
+  //print(f.quadraticPart(), "f.quadraticPart: ");
+
+  LocalVector ri = f.originalLinearPart()[i];
+  const auto& Ai = f.quadraticPart()[i];
+  sparseRangeFor(Ai, [&](auto&& Aij, auto&& j) {
+      // TODO Here we must implement a wrapper to guarantee that this will work with proxy matrices!
+      H::ifElse(H::equals(j, i), [&](auto&& id){
+        Aii_p = id(&Aij);
+      });
+      Dune::TNNMG::Imp::mmv(Aij, f.origin()[j], ri, Dune::PriorityTag<1>());
+  });
+
+  //print(*Aii_p, "Aii_p:");
+  //print(ri, "ri:");
+
+  //print(f.originalLowerObstacle()[i], "lower:");
+  //print(f.originalUpperObstacle()[i], "upper:");
+
+  auto& phii = f.phi().restriction(i);
+
+  return ShiftedFunctional<LocalMatrix&, Range, LocalVector, std::decay_t<decltype(phii)>, LocalLowerObstacle&, LocalUpperObstacle&, LocalVector&, Range>(*Aii_p, f.maxEigenvalues()[i], std::move(ri), phii, f.originalLowerObstacle()[i], f.originalUpperObstacle()[i], f.origin()[i]);
+}
+
+
+/** \brief Box constrained quadratic functional with nonlinearity
+ *
+ *  \tparam M Matrix type
+ *  \tparam V Vector type
+ *  \tparam L Lower obstacle type
+ *  \tparam U Upper obstacle type
+ *  \tparam R Range type
+ */
+template<class M, class V, class N, class L, class U, class R>
+class Functional : public Dune::TNNMG::BoxConstrainedQuadraticFunctional<M, V, L, U, R>
+{
+private:
+  using Base = Dune::TNNMG::BoxConstrainedQuadraticFunctional<M, V, L, U, R>;
+
+
+public:
+  using Nonlinearity = std::decay_t<N>;
+
+  using Matrix = typename Base::Matrix;
+  using Vector = typename Base::Vector;
+  using Range = typename Base::Range;
+  using LowerObstacle = typename Base::LowerObstacle;
+  using UpperObstacle = typename Base::UpperObstacle;
+
+  using Eigenvalues = std::vector<Range>;
+
+private:
+    Dune::Solvers::ConstCopyOrReference<N> phi_;
+    Eigenvalues maxEigenvalues_;
+
+public:
+
+  template <class MM, class VV, class NN, class LL, class UU>
+  Functional(
+          MM&& matrix,
+          VV&& linearPart,
+          NN&& phi,
+          LL&& lower,
+          UU&& upper) :
+    Base(std::forward<MM>(matrix), std::forward<VV>(linearPart), std::forward<LL>(lower), std::forward<UU>(upper)),
+    phi_(std::forward<NN>(phi)),
+    maxEigenvalues_(this->linearPart().size())
+  {
+     for (size_t i=0; i<this->quadraticPart().N(); ++i) {
+       const auto& Aii = this->quadraticPart()[i][i];
+       maxEigenvalues_[i] = maxEigenvalue(Aii);
+
+       // due to numerical imprecision, Aii might not be symmetric leading to complex eigenvalues
+       // consider Toeplitz decomposition of Aii and take only symmetric part
+
+       /*auto symBlock = Aii;
+
+       for (size_t j=0; j<symBlock.N(); j++) {
+           for (size_t k=0; k<symBlock.M(); k++) {
+               if (symBlock[j][k]/symBlock[j][j] < 1e-14)
+                symBlock[j][k] = 0;
+           }
+       }
+
+       try {
+       typename Vector::block_type eigenvalues;
+       Dune::FMatrixHelp::eigenValues(symBlock, eigenvalues);
+       maxEigenvalues_[i] =
+           *std::max_element(std::begin(eigenvalues), std::end(eigenvalues));
+       } catch (Dune::Exception &e) {
+            print(symBlock, "symBlock");
+            typename Vector::block_type eigenvalues;
+            Dune::FMatrixHelp::eigenValues(symBlock, eigenvalues);
+            std::cout << "complex eig in dof: " << i << std::endl;
+       } */
+     }
+  }
+
+  const Nonlinearity& phi() const {
+    return phi_.get();
+  }
+
+  const auto& maxEigenvalues() const {
+    return maxEigenvalues_;
+  }
+
+  Range operator()(const Vector& v) const
+  { 
+    //print(v, "v:");
+    //print(this->lowerObstacle(), "lower: ");
+    //print(this->upperObstacle(), "upper: ");
+    //std::cout << Dune::TNNMG::checkInsideBox(v, this->lowerObstacle(), this->upperObstacle()) << " " << Dune::TNNMG::QuadraticFunctional<M,V,R>::operator()(v) << " " << phi_.get().operator()(v) << std::endl;
+    if (Dune::TNNMG::checkInsideBox(v, this->lowerObstacle(), this->upperObstacle()))
+      return Dune::TNNMG::QuadraticFunctional<M,V,R>::operator()(v) + phi_.get()(v);
+    return std::numeric_limits<Range>::max();
+  }
+
+  friend auto directionalRestriction(const Functional& f, const Vector& origin, const Vector& direction)
+    -> DirectionalRestriction<Matrix, Vector, Nonlinearity, LowerObstacle, UpperObstacle, Range>
+  {
+    // rescale direction for numerical stability
+    auto dir = direction;
+    auto dirNorm = dir.two_norm();
+    if (dirNorm > 0.0)
+        dir /= dirNorm;
+    else {
+        dirNorm = 0.0;
+        dir = 0.0;
+    }
+
+
+    return DirectionalRestriction<Matrix, Vector, Nonlinearity, LowerObstacle, UpperObstacle, Range>(f.quadraticPart(), f.linearPart(), f.phi(), f.lowerObstacle(), f.upperObstacle(), origin, dir, dirNorm);
+  }
+
+  friend auto shift(const Functional& f, const Vector& origin)
+    -> ShiftedFunctional<Matrix&, Eigenvalues, Vector&, Nonlinearity&, LowerObstacle&, UpperObstacle&, Vector&, Range>
+  {
+    return ShiftedFunctional<Matrix&, Eigenvalues, Vector&, Nonlinearity&, LowerObstacle&, UpperObstacle&, Vector&, Range>(f.quadraticPart(), f.maxEigenvalues(), f.linearPart(), f.phi(), f.lowerObstacle(), f.upperObstacle(), origin);
+  }
+};
+
+
+#endif

dune/tectonic/spatial-solving/tnnmg/linearcorrection.hh

+#ifndef DUNE_TECTONIC_ITERATIONSTEPS_LINEARCORRECTION_HH
+#define DUNE_TECTONIC_ITERATIONSTEPS_LINEARCORRECTION_HH
+
+#include <functional>
+#include <memory>
+
+#include <dune/solvers/iterationsteps/lineariterationstep.hh>
+#include <dune/solvers/solvers/iterativesolver.hh>
+#include <dune/solvers/solvers/linearsolver.hh>
+#include <dune/solvers/common/canignore.hh>
+
+#include <dune/solvers/common/resize.hh>
+#include <dune/solvers/solvers/umfpacksolver.hh>
+
+/**
+ * \brief linear correction step for use by \ref TNNMGStep
+ *
+ * The function object should solve the linear equation \f$ A x = b \f$.
+ * The ignore bitfield identifies the dofs to be truncated by the
+ * linear solver. The truncation dofs are passed by setting its ignore
+ * member.
+ * Be aware that full rows or columns of `A` might contain only zeroes.
+ */
+template<typename Matrix, typename Vector>
+using LinearCorrection = std::function<void(const Matrix& A, Vector& x, const Vector& b, const Dune::Solvers::DefaultBitVector_t<Vector>& ignore)>;
+
+template<typename Vector>
+Dune::Solvers::DefaultBitVector_t<Vector>
+emptyIgnore(const Vector& v)
+{
+  // TNNMGStep assumes that the linearization and the solver for the
+  // linearized problem will not use the ignoreNodes field
+  Dune::Solvers::DefaultBitVector_t<Vector> ignore;
+  Dune::Solvers::resizeInitialize(ignore, v, false);
+  return ignore;
+}
+
+template<typename Matrix, typename Vector>
+LinearCorrection<Matrix, Vector>
+buildLinearCorrection(std::shared_ptr< Dune::Solvers::LinearSolver<Matrix, Vector> > linearSolver)
+{
+  return [=](const Matrix& A, Vector& x, const Vector& b, const Dune::Solvers::DefaultBitVector_t<Vector>& ignore) {
+
+    auto canIgnoreCast = std::dynamic_pointer_cast<CanIgnore<Dune::Solvers::DefaultBitVector_t<Vector>>>( linearSolver );
+    if (canIgnoreCast)
+      canIgnoreCast->setIgnore(ignore);
+    else
+      DUNE_THROW(Dune::Exception, "LinearCorrection: linearSolver cannot set ignore member!");
+
+    linearSolver->setProblem(A, x, b);
+    linearSolver->preprocess();
+    linearSolver->solve();
+  };
+}
+
+template<typename Matrix, typename Vector>
+LinearCorrection<Matrix, Vector>
+buildLinearCorrection(std::shared_ptr< Dune::Solvers::IterativeSolver<Vector> > iterativeSolver)
+{
+  return [=](const Matrix& A, Vector& x, const Vector& b, const Dune::Solvers::DefaultBitVector_t<Vector>& ignore) {
+
+      // compute reference solution directly
+      LocalProblem<Matrix, Vector> localProblem(A, b, ignore);
+      Vector newR, directX;
+      Dune::Solvers::resizeInitializeZero(directX, b);
+      localProblem.getLocalRhs(directX, newR);
+
+      /*print(*this->ignoreNodes_, "ignoreNodes:");
+      print(A, "A:");
+      print(localProblem.getMat(), "localMat:");*/
+
+      auto directSolver = std::make_shared<Dune::Solvers::UMFPackSolver<Matrix, Vector>>();
+
+      directSolver->setProblem(localProblem.getMat(), directX, newR);
+      directSolver->preprocess();
+      directSolver->solve();
+
+      //x = directX;
+    using LinearIterationStep = Dune::Solvers::LinearIterationStep<Matrix, Vector>;
+
+    auto linearIterationStep = dynamic_cast<LinearIterationStep*>(&iterativeSolver->getIterationStep());
+    if (not linearIterationStep)
+      DUNE_THROW(Dune::Exception, "iterative solver must use a linear iteration step");
+
+    auto empty = emptyIgnore(x);
+    linearIterationStep->setIgnore(empty);
+    //linearIterationStep->setIgnore(ignore);
+    linearIterationStep->setProblem(A, x, b);
+    iterativeSolver->preprocess();
+    iterativeSolver->solve();
+
+
+    const auto& norm = iterativeSolver->getErrorNorm();
+    auto error = norm.diff(linearIterationStep->getSol(), directX);
+
+    std::cout << "Linear solver iterations: " << iterativeSolver->getResult().iterations << " Error: " << error << std::endl;
+  };
+}
+
+template<typename Matrix, typename Vector>
+LinearCorrection<Matrix, Vector>
+buildLinearCorrection(std::shared_ptr< Dune::Solvers::LinearIterationStep<Matrix, Vector> > linearIterationStep, int nIterationSteps = 1)
+{
+  return [=](const Matrix& A, Vector& x, const Vector& b, const Dune::Solvers::DefaultBitVector_t<Vector>& ignore) {
+    linearIterationStep->setIgnore(ignore);
+    linearIterationStep->setProblem(A, x, b);
+    linearIterationStep->preprocess();
+
+    for (int i = 0; i < nIterationSteps; ++i)
+      linearIterationStep->iterate();
+  };
+}
+
+#endif

dune/tectonic/spatial-solving/tnnmg/linearization.hh

+#ifndef DUNE_TECTONIC_LINEARIZATION_HH
+#define DUNE_TECTONIC_LINEARIZATION_HH
+
+#include <cstddef>
+
+#include <dune/common/typetraits.hh>
+#include <dune/common/hybridutilities.hh>
+
+#include <dune/matrix-vector/algorithm.hh>
+
+#include <dune/istl/matrixindexset.hh>
+
+//#include <dune/tnnmg/functionals/bcqfconstrainedlinearization.hh>
+
+#include "../../utils/debugutils.hh"
+
+/**
+ * derived from Dune::TNNMG::BoxConstrainedQuadraticFunctionalConstrainedLinearization<F,BV>
+ */
+template<class F, class BV>
+class Linearization {
+public:
+  using Matrix = typename F::Matrix;
+  using Vector = typename F::Vector;
+  using BitVector = BV;
+  using ConstrainedMatrix = Matrix;
+  using ConstrainedVector = Vector;
+  using ConstrainedBitVector = BitVector;
+
+private:
+  using This = Linearization<F,BV>;
+
+  template <class T>
+  void determineRegularityTruncation(const Vector& x, BitVector& truncationFlags, const T& truncationTolerance)
+  {
+    //std::cout << "x: " << x << std::endl;
+    size_t blocksize = truncationFlags[0].size();
+
+    size_t count = 0;
+
+    for (size_t i = 0; i < x.size(); ++i) {
+      //if (truncationFlags[i].all())
+        //  continue;
+
+      //std::cout << f_.phi().regularity(i, x[i]) << " xnorm: " << x[i].two_norm() << std::endl;
+
+      if (f_.phi().regularity(i, x[i]) > truncationTolerance) {
+        for (size_t j=1; j<blocksize; j++) {
+            truncationFlags[i][j] = true;
+        }
+        count++;
+      }
+    }
+
+     std::cout << "regularityTruncation: " << count << std::endl;
+  }
+
+  template<class NV, class NBV, class T>
+  static void determineTruncation(const NV& x, const NV& lower, const NV& upper, NBV&& truncationFlags, const T& truncationTolerance)
+  {
+    namespace H = Dune::Hybrid;
+    H::ifElse(Dune::IsNumber<NV>(), [&](auto id){
+      if ((id(x) <= id(lower)+truncationTolerance) || (id(x) >= id(upper) - truncationTolerance)) {
+        id(truncationFlags) = true;
+      }
+    }, [&](auto id){
+      H::forEach(H::integralRange(H::size(id(x))), [&](auto&& i) {
+        This::determineTruncation(x[i], lower[i], upper[i], truncationFlags[i], truncationTolerance);
+      });
+    });
+  }
+
+  template<class NV, class NBV>
+  static void truncateVector(NV& x, const NBV& truncationFlags)
+  {
+    namespace H = Dune::Hybrid;
+    H::ifElse(Dune::IsNumber<NV>(), [&](auto id){
+      if (id(truncationFlags))
+        id(x) = 0;
+    }, [&](auto id){
+      H::forEach(H::integralRange(H::size(id(x))), [&](auto&& i) {
+        This::truncateVector(x[i], truncationFlags[i]);
+      });
+    });
+  }
+
+  template<class NM, class RBV, class CBV>
+  static void truncateMatrix(NM& A, const RBV& rowTruncationFlags, const CBV& colTruncationFlags)
+  {
+    namespace H = Dune::Hybrid;
+    using namespace Dune::MatrixVector;
+    H::ifElse(Dune::IsNumber<NM>(), [&](auto id){
+      if(id(rowTruncationFlags) or id(colTruncationFlags))
+        A = 0;
+    }, [&](auto id){
+      H::forEach(H::integralRange(H::size(id(rowTruncationFlags))), [&](auto&& i) {
+        auto&& Ai = A[i];
+        sparseRangeFor(Ai, [&](auto&& Aij, auto&& j) {
+            This::truncateMatrix(Aij, rowTruncationFlags[i], colTruncationFlags[j]);
+        });
+      });
+    });
+  }
+
+public:
+  Linearization(const F& f, const BitVector& ignore) :
+      f_(f),
+      ignore_(ignore),
+      obstacleTruncationTolerance_(1e-10),
+      regularityTruncationTolerance_(1e8)
+  {}
+
+  void bind(const Vector& x) {
+      //std::cout << "Linearization::bind()" << std::endl;
+
+      auto&& A = f_.quadraticPart();
+      auto&& phi = f_.phi();
+
+      // determine which components to truncate
+      // ----------------
+      BitVector obstacleTruncationFlags = ignore_;
+      BitVector regularityTruncationFlags = ignore_;
+
+      //std::cout << "ignore truncation: " << truncationFlags_.count();
+
+      determineTruncation(x, f_.lowerObstacle(), f_.upperObstacle(), obstacleTruncationFlags, obstacleTruncationTolerance_); // obstacle truncation
+
+      //print(obstacleTruncationFlags, "obstacleTruncationFlags:");
+
+      //std::cout << "    obstacle truncation: " << truncationFlags_.count();
+      //std::cout << " " << obstacleTruncationTolerance_;
+
+      determineRegularityTruncation(x, regularityTruncationFlags, regularityTruncationTolerance_); // truncation due to regularity deficit of nonlinearity
+
+      //std::cout << "    regularity truncation: " << truncationFlags_.count() << std::endl;
+      //std::cout << " " << regularityTruncationTolerance_;
+
+      //print(regularityTruncationFlags, "regularityTruncationFlags:");
+
+      truncationFlags_ = obstacleTruncationFlags;
+      size_t blocksize = truncationFlags_[0].size();
+      for (size_t i=0; i<truncationFlags_.size(); i++) {
+          for (size_t j=0; j<blocksize; j++) {
+              truncationFlags_[i][j] = truncationFlags_[i][j] or regularityTruncationFlags[i][j];
+          }
+      }
+
+      // compute hessian
+      // ---------------
+
+      // construct sparsity pattern for linearization
+      Dune::MatrixIndexSet indices(A.N(), A.M());
+      indices.import(A);
+      phi.addHessianIndices(indices);
+
+      // construct matrix from pattern and initialize it
+      indices.exportIdx(hessian_);
+      hessian_ = 0.0;
+
+      // quadratic part
+      for (size_t i = 0; i < A.N(); ++i) {
+        auto const end = std::end(A[i]);
+        for (auto it = std::begin(A[i]); it != end; ++it)
+          hessian_[i][it.index()] += *it;
+      }
+
+      //truncateMatrix(hessian_, obstacleTruncationFlags, obstacleTruncationFlags);
+
+      //print(hessian_, "Hessian: ");
+      //std::cout << "--------" << (double) hessian_[21][21] << "------------" << std::endl;
+
+      //std::cout << hessian_[0][0] << std::endl;
+
+      auto B = hessian_;
+
+      //print(x, "x: ");
+
+      // nonlinearity part
+      phi.addHessian(x, hessian_);
+      //truncateMatrix(hessian_, regularityTruncationFlags, regularityTruncationFlags);
+
+      //std::cout << x[0] << " " << hessian_[0][0] << std::endl;
+
+      B -= hessian_;
+      //print(B, "phi hessian: ");
+
+      // compute gradient
+      // ----------------
+
+      // quadratic part
+      negativeGradient_ = derivative(f_)(x); // A*x - b
+
+      //print(negativeGradient_, "gradient: ");
+      auto C = negativeGradient_;
+
+
+      // nonlinearity part
+      phi.addGradient(x, negativeGradient_);
+
+      C -= negativeGradient_;
+      //print(C, "phi gradient: ");
+
+      // -grad is needed for Newton step
+      negativeGradient_ *= -1;
+
+      // truncate gradient
+      truncateVector(negativeGradient_, truncationFlags_);
+      truncateMatrix(hessian_, truncationFlags_, truncationFlags_);
+
+      //print(hessian_, "hessian: ");
+      //print(negativeGradient_, "negativeGradient: ");
+
+      //print(truncationFlags_, "truncationFlags:");
+  }
+
+
+  void extendCorrection(ConstrainedVector& cv, Vector& v) const {
+    v = cv;
+    truncateVector(v, truncationFlags_);
+  }
+
+  /*const BitVector& truncated() const {
+    return truncationFlags_;
+  }*/
+
+  const auto& negativeGradient() const {
+    return negativeGradient_;
+  }
+
+  const auto& hessian() const {
+    return hessian_;
+  }
+
+private:
+  const F& f_;
+  const BitVector& ignore_;
+
+  double obstacleTruncationTolerance_;
+  double regularityTruncationTolerance_;
+
+  Vector negativeGradient_;
+  Matrix hessian_;
+  BitVector truncationFlags_;
+};
+
+#endif

dune/tectonic/spatial-solving/tnnmg/linesearchsolver.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_TECTONIC_LINESEARCHSOLVER_HH
+#define DUNE_TECTONIC_LINESEARCHSOLVER_HH
+
+#include <dune/tnnmg/problem-classes/bisection.hh>
+
+#include "../../utils/almostequal.hh"
+
+/**
+ * \brief A local solver for scalar quadratic obstacle problems with nonlinearity
+ * using bisection
+ */
+class LineSearchSolver
+{
+public:
+  template<class Vector, class Functional, class BitVector>
+  void operator()(Vector& x, const Functional& f, const BitVector& ignore) const {
+    x = 1.0;
+
+    if (ignore)
+        return;
+
+    Dune::Solvers::Interval<double> D;
+    D = f.subDifferential(0);
+
+    std::cout << "f.A " << f.quadraticPart() << " f.b " << f.linearPart() << std::endl;
+
+    std::cout << D[0] << " " << D[1] << std::endl;
+    std::cout << "domain: " << f.domain()[0] << " " << f.domain()[1] << std::endl;
+
+
+    if (D[1] > 0) // NOTE: Numerical instability can actually get us here
+        return;
+
+
+    if (almost_equal(f.domain()[0], f.domain()[1], 2)) {
+        x = f.domain()[0];
+        std::cout << "no interval: " << x << std::endl;
+        return;
+    }
+
+    int bisectionsteps = 0;
+    const Bisection globalBisection; //(0.0, 1.0, 0.0, 0.0);
+
+    x = globalBisection.minimize(f, f.scaling(), 0.0, bisectionsteps);
+    std::cout << "x: " << x << "scaling: " << f.scaling();
+    x /= f.scaling();
+    std::cout << "final x: "  << x << std::endl;
+    //x = f.domain().projectIn(x);
+  }
+};
+
+#endif
+

dune/tectonic/spatial-solving/tnnmg/localbisectionsolver.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_TECTONIC_LOCALBISECTIONSOLVER_HH
+#define DUNE_TECTONIC_LOCALBISECTIONSOLVER_HH
+
+#include <dune/matrix-vector/axpy.hh>
+
+#include <dune/tnnmg/localsolvers/scalarobstaclesolver.hh>
+#include <dune/tnnmg/problem-classes/bisection.hh>
+
+#include "functional.hh"
+
+#include "../../utils/debugutils.hh"
+
+/**
+ * \brief A local solver for quadratic obstacle problems with nonlinearity
+ * using bisection
+ */
+class LocalBisectionSolver
+{
+public:
+  template<class Vector, class Functional, class BitVector>
+  void operator()(Vector& x, const Functional& f, const BitVector& ignore) const {
+      double safety = 1e-14; //or (f.upperObstacle()-f.lowerObstacle()<safety)
+      if (ignore.all()) {
+          return;
+      }
+
+      auto maxEigenvalue = f.maxEigenvalues();
+
+      auto origin = f.origin();
+      auto linearPart = f.originalLinearPart();
+
+      Dune::MatrixVector::addProduct(linearPart, maxEigenvalue, origin);
+
+      for (size_t j = 0; j < ignore.size(); ++j)
+        if (ignore[j])
+          linearPart[j] = 0; // makes sure result remains in subspace after correction
+        else
+          origin[j] = 0; // shift affine offset
+
+      double const linearPartNorm = linearPart.two_norm();
+      if (linearPartNorm <= 0.0)
+        return;
+      linearPart /= linearPartNorm;
+
+      //std::cout << "maxEigenvalue " << maxEigenvalue << std::endl;
+      //std::cout << "linearPartNorm " << linearPartNorm << std::endl;
+      //print(origin, "origin:");
+      //print(linearPart, "linearPart:");
+
+      using Nonlinearity = std::decay_t<decltype(f.phi())>;
+      using Range = std::decay_t<decltype(f.maxEigenvalues())>;
+      using FirstOrderFunctional = FirstOrderFunctional<Nonlinearity, Vector, Range>;
+
+      auto lower = f.originalLowerObstacle();
+      auto upper = f.originalUpperObstacle();
+
+      //print(lower, "lower:");
+      //print(upper, "upper:");
+
+      FirstOrderFunctional firstOrderFunctional(maxEigenvalue, linearPartNorm, f.phi(),
+                                                lower, upper,
+                                                origin, linearPart);
+
+      //std::cout << "lower: " << firstOrderFunctional.lowerObstacle() << " upper: " << firstOrderFunctional.upperObstacle() << std::endl;
+
+      // scalar obstacle solver without nonlinearity
+      typename Functional::Range alpha;
+      //Dune::TNNMG::ScalarObstacleSolver obstacleSolver;
+      //obstacleSolver(alpha, firstOrderFunctional, false);
+
+      //direction *= alpha;
+
+      /*const auto& A = f.quadraticPart();
+      std::cout << "f.quadratic(): " << std::endl;
+      for (size_t i=0; i<A.N(); i++) {
+            for (size_t j=0; j<A.M(); j++) {
+                std::cout << A[i][j] << " ";
+            }
+            std::cout << std::endl;
+      }*/
+      //std::cout << f.quadraticPart() << std::endl;
+      //std::cout << "A: " << directionalF.quadraticPart() << " " << (directionalF.quadraticPart()==f.quadraticPart()[0][0]) << std::endl;
+      /*std::cout << "b: " << directionalF.linearPart() << std::endl;
+      std::cout << "domain: " << directionalF.domain()[0] << " " << directionalF.domain()[1] << std::endl;
+      auto D = directionalF.subDifferential(0);
+      std::cout << "subDiff at x: " << D[0] << " " << D[1] << std::endl;*/
+
+      //std::cout << "domain: " << firstOrderFunctional.domain()[0] << " " << firstOrderFunctional.domain()[1] << std::endl;
+
+
+      const auto& domain = firstOrderFunctional.domain();
+      if (std::abs(domain[0]-domain[1])>safety) {
+        int bisectionsteps = 0;
+        const Bisection globalBisection(0.0, 1.0, 0.0, 0.0);
+
+        alpha = globalBisection.minimize(firstOrderFunctional, 0.0, 0.0, bisectionsteps);
+      } else {
+          alpha = domain[0];
+      }
+
+      linearPart *= alpha;
+
+      //std::cout << linearPart << std::endl;
+
+#ifdef NEW_TNNMG_COMPUTE_ITERATES_DIRECTLY
+      if (std::abs(alpha)> safety) {
+        x = origin;
+        x += linearPart;
+      } else {
+        x = f.origin();
+      }
+#else
+      if (std::abs(alpha)> safety) {
+        x = origin;
+        x += linearPart;
+        x -= f.origin();
+      }
+#endif
+  }
+};
+
+#endif
+
+
+
+

dune/tectonic/spatial-solving/tnnmg/zerononlinearity.hh

+#ifndef DUNE_TECTONIC_ZERO_NONLINEARITY_HH
+#define DUNE_TECTONIC_ZERO_NONLINEARITY_HH
+
+#include <dune/solvers/common/interval.hh>
+
+/** \file
+ * \brief A dummy nonlinearity class representing the zero function
+ */
+
+class ZeroNonlinearity
+{
+public:
+    ZeroNonlinearity()
+    {}
+
+    const auto& restriction(size_t i) const {
+        return *this;
+    }
+
+    /** \brief Returns zero */
+    template <class VectorType>
+    double operator()(const VectorType& v) const
+    {
+        return 0.0;
+    }
+
+    template <class VectorType>
+    void addGradient(const VectorType& v, VectorType& gradient) const {}
+
+    template <class VectorType, class MatrixType>
+    void addHessian(const VectorType& v, MatrixType& hessian) const {}
+
+    template <class IndexSet>
+    void addHessianIndices(IndexSet& indices) const {}
+
+    /** \brief Returns the interval \f$ [0,0]\f$ */
+    void subDiff(int i, double x, Dune::Solvers::Interval<double>& D) const
+    {
+        D[0] = 0.0;
+        D[1] = 0.0;
+    }
+
+    template <class VectorType>
+    void directionalSubDiff(const VectorType& u, const VectorType& v, Dune::Solvers::Interval<double>& subdifferential) const
+    {
+        subdifferential[0] = 0.0;
+        subdifferential[1] = 0.0;
+    }
+
+    /** \brief Returns 0 */
+    template <class VectorType>
+    double regularity(int i, const VectorType& x) const
+    {
+        return 0.0;
+    }
+
+    void domain(int i, Dune::Solvers::Interval<double>& dom) const
+    {
+        dom[0] = -std::numeric_limits<double>::max();
+        dom[1] = std::numeric_limits<double>::max();
+        return;
+    }
+
+    template <class VectorType>
+    void directionalDomain(const VectorType& u, const VectorType& v, Dune::Solvers::Interval<double>& dom) const
+    {
+        dom[0] = -std::numeric_limits<double>::max();
+        dom[1] = std::numeric_limits<double>::max();
+        return;
+    }
+
+    template <class BitVector>
+    void setIgnore(const BitVector& ignore) {}
+
+    template <class StateVector>
+    void updateAlpha(const StateVector& alpha) {}
+};
+
+#endif
+

dune/tectonic/tectonic.hh

-#ifndef DUNE_TECTONIC_TECTONIC_HH
-#define DUNE_TECTONIC_TECTONIC_HH
-
-// add your classes here
-#endif

dune/tectonic/tests/CMakeLists.txt

+add_subdirectory("contacttest")
+add_subdirectory("hdf5test")
+#add_subdirectory("tnnmgtest")
+
+dune_add_test(SOURCES globalfrictioncontainertest.cc)
+dune_add_test(SOURCES gridgluefrictiontest.cc)
+dune_add_test(SOURCES nodalweightstest.cc)
+dune_add_test(SOURCES supportpatchfactorytest.cc)
+dune_add_test(SOURCES solverfactorytest.cc)

dune/tectonic/tests/common.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_TECTONIC_SRC_TESTS_COMMON_HH
+#define DUNE_TECTONIC_SRC_TESTS_COMMON_HH
+
+#include <dune/grid/common/gridfactory.hh>
+#include <dune/grid/utility/structuredgridfactory.hh>
+
+#include <dune/fufem/referenceelementhelper.hh>
+
+// utility functions
+bool isClose(double a, double b) {
+  return std::abs(a - b) < 1e-14;
+}
+
+template <class Vector>
+bool xyCollinear(Vector const &a, Vector const &b, Vector const &c) {
+  return isClose((b[0] - a[0]) * (c[1] - a[1]), (b[1] - a[1]) * (c[0] - a[0]));
+}
+
+template <class Vector>
+bool xyBoxed(Vector const &v1, Vector const &v2, Vector const &x) {
+  auto const minmax0 = std::minmax(v1[0], v2[0]);
+  auto const minmax1 = std::minmax(v1[1], v2[1]);
+
+  if (minmax0.first - 1e-14 > x[0] or
+      x[0] > minmax0.second + 1e-14)
+    return false;
+  if (minmax1.first - 1e-14 > x[1] or
+      x[1] > minmax1.second + 1e-14)
+    return false;
+
+  return true;
+}
+
+template <class Vector>
+bool xyBetween(Vector const &v1, Vector const &v2, Vector const &x) {
+  return xyCollinear(v1, v2, x) && xyBoxed(v1, v2, x);
+}
+
+// build cube grid given by lowerLeft and upperRight point in global coordinates
+template <class GlobalCoords, class GridType>
+void buildGrid(const GlobalCoords& lowerLeft, const GlobalCoords& upperRight, const int n, std::shared_ptr<GridType>& grid, bool simplexGrid = true) {
+    std::array<unsigned int, GridType::dimension> elements;
+    std::fill(elements.begin(), elements.end(), n);
+
+    Dune::GridFactory<GridType> factory;
+    if (!simplexGrid) {
+        Dune::StructuredGridFactory<GridType>::createCubeGrid(factory, lowerLeft, upperRight, elements);
+        grid = std::move(factory.createGrid());
+    } else {
+        Dune::StructuredGridFactory<GridType>::createSimplexGrid(factory, lowerLeft, upperRight, elements);
+        grid = std::move(factory.createGrid());
+    }
+}
+
+template <class Intersection>
+bool containsInsideSubentity(const Intersection& nIt, int subEntity, int codim) {
+    return ReferenceElementHelper<double, Intersection::GridView::dim>::subEntityContainsSubEntity(nIt.inside().type(), nIt.indexInInside(), 1, subEntity, codim);
+}
+
+#endif