.

src/data-structures/program_state.hh

@@ -10,11 +10,17 @@
 
 #include <dune/contact/assemblers/nbodyassembler.hh>
 
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+#include <dune/solvers/iterationsteps/cgstep.hh>
+#include <dune/solvers/iterationsteps/blockgssteps.hh>
+
 #include <dune/tectonic/bodydata.hh>
 
 #include "../assemblers.hh"
 #include "contactnetwork.hh"
 #include "matrices.hh"
+#include "../spatial-solving/preconditioners/multilevelpatchpreconditioner.hh"
 #include "../spatial-solving/solverfactory.hh"
 #include "../spatial-solving/tnnmg/functional.hh"
 #include "../spatial-solving/tnnmg/zerononlinearity.hh"
@@ -51,6 +57,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
     using LocalVector = typename Vector::block_type;
     //using LocalMatrix = typename Matrix::block_type;
     const static int dims = LocalVector::dimension;
+    using BitVector = Dune::BitSetVector<dims>;
 
 public:
   ProgramState(const std::vector<size_t>& leafVertexCounts)
@@ -85,14 +92,43 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
   }
 
   // Set up initial conditions
-  template <class GridType>
-  void setupInitialConditions(const Dune::ParameterTree& parset, const ContactNetwork<GridType, Vector>& contactNetwork) {
-    using Matrix = typename ContactNetwork<GridType, Vector>::Matrix;
+  template <class ContactNetwork>
+  void setupInitialConditions(const Dune::ParameterTree& parset, const ContactNetwork& contactNetwork) {
+    using Matrix = typename ContactNetwork::Matrix;
     const auto& nBodyAssembler = contactNetwork.nBodyAssembler();
 
+    // set patch preconditioner for linear correction in TNNMG method
+    using PatchSolver = typename Dune::Solvers::LoopSolver<Vector, BitVector>;
+    using Preconditioner = MultilevelPatchPreconditioner<ContactNetwork, PatchSolver, Matrix, Vector>;
+
+    const auto& preconditionerParset = parset.sub("solver.tnnmg.linear.preconditioner");
+
+    auto gsStep = Dune::Solvers::BlockGSStepFactory<Matrix, Vector>::create(Dune::Solvers::BlockGS::LocalSolvers::direct(0.0));
+    PatchSolver patchSolver(gsStep,
+                               preconditionerParset.get<size_t>("maximumIterations"),
+                               preconditionerParset.get<double>("tolerance"),
+                               nullptr,
+                               preconditionerParset.get<Solver::VerbosityMode>("verbosity"),
+                               false); // absolute error
+
+    Dune::BitSetVector<1> activeLevels(contactNetwork.nLevels(), true);
+    Preconditioner preconditioner(contactNetwork, activeLevels, preconditionerParset.get<MPPMode>("mode"));
+    preconditioner.setPatchSolver(std::forward<PatchSolver>(patchSolver));
+    preconditioner.setPatchDepth(preconditionerParset.get<size_t>("patchDepth"));
+    preconditioner.build();
+
+    using LinearSolver = typename Dune::Solvers::LoopSolver<Vector, BitVector>;
+    using LinearSolverStep = typename Dune::Solvers::CGStep<Matrix, Vector, BitVector>;
+
+    LinearSolverStep linearSolverStep;
+    linearSolverStep.setPreconditioner(preconditioner);
+
+    EnergyNorm<Matrix, Vector> energyNorm(linearSolverStep);
+    LinearSolver linearSolver(linearSolverStep, parset.get<size_t>("linear.maximumIterations"), parset.get<double>("linear.tolerance"), energyNorm, Solver::QUIET);
+
     // Solving a linear problem with a multigrid solver
     auto const solveLinearProblem = [&](
-        const Dune::BitSetVector<dims>& _dirichletNodes, const std::vector<std::shared_ptr<Matrix>>& _matrices,
+        const BitVector& _dirichletNodes, const std::vector<std::shared_ptr<Matrix>>& _matrices,
         const std::vector<Vector>& _rhs, std::vector<Vector>& _x,
         Dune::ParameterTree const &_localParset) {
 
@@ -149,14 +185,15 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
       print(lower, "lower");
       print(upper, "upper");*/
 
-      // set up functional and TNMMG solver
+      // set up functional
       using Functional = Functional<Matrix&, Vector&, ZeroNonlinearity&, Vector&, Vector&, typename Matrix::field_type>;
       Functional J(bilinearForm, totalRhs, ZeroNonlinearity(), lower, upper); //TODO
 
-    //  using Factory = SolverFactory<Functional, Dune::BitSetVector<dims>>;
-    //  Factory factory(parset.sub("solver.tnnmg"), J, _dirichletNodes);
+      // set up TNMMG solver
+      using Factory = SolverFactory<Functional, BitVector>;
+      Factory factory(parset.sub("solver.tnnmg"), J, linearSolver, _dirichletNodes);
 
-   /*   std::vector<Dune::BitSetVector<dims>> bodyDirichletNodes;
+     /* std::vector<BitVector> bodyDirichletNodes;
       nBodyAssembler.postprocess(_dirichletNodes, bodyDirichletNodes);
       for (size_t i=0; i<bodyDirichletNodes.size(); i++) {
         print(bodyDirichletNodes[i], "bodyDirichletNodes_" + std::to_string(i) + ": ");
@@ -166,7 +203,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
       print(totalX, "totalX: ");
       print(totalRhs, "totalRhs: ");*/
 
-     /* auto tnnmgStep = factory.step();
+      auto tnnmgStep = factory.step();
       factory.setProblem(totalX);
 
       const EnergyNorm<Matrix, Vector> norm(bilinearForm);
@@ -184,7 +221,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
 
       Vector res = totalRhs;
       bilinearForm.mmv(tnnmgStep->getSol(), res);
-      std::cout << "TNNMG Res - energy norm: " << norm.operator ()(res) << std::endl; */
+      std::cout << "TNNMG Res - energy norm: " << norm.operator ()(res) << std::endl;
 
       // TODO: remove after debugging
   /*    using DeformedGridType = typename LevelContactNetwork<GridType, dims>::DeformedGridType;
@@ -233,7 +270,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
     // Initial displacement: Start from a situation of minimal stress,
     // which is automatically attained in the case [v = 0 = a].
     // Assuming dPhi(v = 0) = 0, we thus only have to solve Au = ell0
-    Dune::BitSetVector<dims> dirichletNodes;
+    BitVector dirichletNodes;
     contactNetwork.totalNodes("dirichlet", dirichletNodes);
     /*for (size_t i=0; i<dirichletNodes.size(); i++) {
         bool val = false;
@@ -262,7 +299,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
       // Initial normal stress
 
       const auto& body = contactNetwork.body(i);
-      std::vector<std::shared_ptr<typename ContactNetwork<GridType, Vector>::LeafBody::BoundaryCondition>> frictionBoundaryConditions;
+      std::vector<std::shared_ptr<typename ContactNetwork::LeafBody::BoundaryCondition>> frictionBoundaryConditions;
       body->boundaryConditions("friction", frictionBoundaryConditions);
       for (size_t j=0; j<frictionBoundaryConditions.size(); j++) {
           ScalarVector frictionBoundaryStress(weightedNormalStress[i].size());
@@ -277,7 +314,7 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
       Dune::MatrixVector::subtractProduct(accelerationRHS[i], *body->matrices().elasticity, u[i]);
     }
 
-    Dune::BitSetVector<dims> noNodes(dirichletNodes.size(), false);
+    BitVector noNodes(dirichletNodes.size(), false);
     solveLinearProblem(noNodes, contactNetwork.matrices().mass, accelerationRHS, a,
                        parset.sub("a0.solver"));
 

src/spatial-solving/preconditioners/levelpatchpreconditioner.hh

@@ -8,6 +8,7 @@
 #include <dune/common/bitsetvector.hh>
 
 #include <dune/solvers/iterationsteps/lineariterationstep.hh>
+#include <dune/solvers/common/numproc.hh>
 
 #include "../../data-structures/levelcontactnetwork.hh"
 
@@ -60,7 +61,7 @@ class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorTy
             patchFactory_(levelContactNetwork_, fineContactNetwork_) {
 
         setPatchDepth();
-        this->verbosity_ = QUIET;
+        this->verbosity_ = NumProc::QUIET;
     }
 
     // build support patches
@@ -74,7 +75,7 @@ class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorTy
 
         // init local fine level corrections
         Dune::Timer timer;
-        if (this->verbosity_ == FULL) {
+        if (this->verbosity_ == NumProc::FULL) {
             std::cout << std::endl;
             std::cout << "---------------------------------------------" << std::endl;
             std::cout << "Initializing local fine grid corrections! Level: " << level_ << std::endl;
@@ -106,14 +107,14 @@ class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorTy
             }
         }
 
-        if (this->verbosity_ == FULL) {
+        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_ != QUIET) {
+        if (this->verbosity_ != NumProc::QUIET) {
             std::cout << "Level: " << level_ << " LevelPatchPreconditioner::build() - SUCCESS" << std::endl;
         }
     }
@@ -126,7 +127,7 @@ class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorTy
         patchDepth_ = patchDepth;
     }
 
-    void setVerbosity(VerbosityMode verbosity) {
+    void setVerbosity(NumProc::VerbosityMode verbosity) {
         this->verbosity_ = verbosity;
     }
 

src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh

@@ -12,15 +12,16 @@
 #include "nbodycontacttransfer.hh"
 #include "levelpatchpreconditioner.hh"
 
+enum MPPMode {additive, multiplicative};
+
 template <class ContactNetwork, class PatchSolver, class MatrixType, class VectorType>
 class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorType> {
 public:
-    enum Mode {additive, multiplicative};
+    //struct Mode { enum ModeType {additive, multiplicative}; };
 
 private:
     using LevelContactNetwork = typename ContactNetwork::LevelContactNetwork;
     using LevelPatchPreconditioner = LevelPatchPreconditioner<LevelContactNetwork, PatchSolver, MatrixType, VectorType>;
-    using LevelGlobalPreconditioner = LevelPatchPreconditioner<LevelContactNetwork, PatchSolver, MatrixType, VectorType>;
     using BitVector = typename LinearIterationStep<MatrixType, VectorType>::BitVector;
 
     using MGTransfer = NBodyContactTransfer<ContactNetwork, VectorType>;
@@ -30,7 +31,7 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
     const ContactNetwork& contactNetwork_;
     const Dune::BitSetVector<1>& activeLevels_;
 
-    const Mode mode_;
+    const MPPMode mode_;
 
     std::vector<std::shared_ptr<LevelPatchPreconditioner>> levelPatchPreconditioners_;
 
@@ -39,15 +40,16 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
     std::vector<VectorType> levelRhs_;
     std::vector<BitVector> levelIgnore_;
 
-    PatchSolver coarseSolver_;
+    std::unique_ptr<PatchSolver> coarseSolver_;
 
     std::vector<BitVector> recompute_;
+    std::vector<BitVector> truncation_;
     std::vector<std::shared_ptr<MGTransfer>> mgTransfer_; //TODO
 
 public:
     MultilevelPatchPreconditioner(const ContactNetwork& contactNetwork,
                                   const Dune::BitSetVector<1>& activeLevels,
-                                  const Mode mode = MultilevelPatchPreconditioner::Mode::additive) :
+                                  const MPPMode mode = additive) :
           contactNetwork_(contactNetwork),
           activeLevels_(activeLevels),
           mode_(mode) {
@@ -55,15 +57,12 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
         assert(activeLevels.size() == contactNetwork_.nLevels());
 
         // init level patch preconditioners and multigrid transfer
-        levelPatchPreconditioners_.resize(0);
+        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]) {
-                // init global problem on coarsest level
-                const auto& levelNetwork = *contactNetwork_.level(i);
-                coarseGlobalPreconditioner_ = std::make_shared<LevelGlobalPreconditioner>(levelNetwork);
-
                 maxLevel = i;
                 break;
             }
@@ -71,16 +70,16 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
 
         const int coarsestLevel = maxLevel;
 
-        typename LevelPatchPreconditioner::Mode levelMode = LevelPatchPreconditioner::Mode::additive;
+        MPPMode levelMode = additive;
         if (mode_ != additive){
-            levelMode = LevelPatchPreconditioner::Mode::multiplicative;
+            levelMode = multiplicative;
         }
 
         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, levelMode));
+                levelPatchPreconditioners_.push_back(std::make_shared<LevelPatchPreconditioner>(*contactNetwork_.level(maxLevel), fineNetwork, levelMode));
 
                 maxLevel = i;
             }
@@ -90,14 +89,14 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
         levelRhs_.resize(size());
 
         // init multigrid transfer
-        mgTransfer_.resize(levelPatchPreconditioners_.size());
+        mgTransfer_.resize(levelPatchPreconditioners_.size()-1);
         for (size_t i=0; i<mgTransfer_.size(); i++) {
             mgTransfer_[i] = std::make_shared<MGTransfer>();
-            mgTransfer_.setup(contactNetwork_, levelPatchPreconditioners_[i]->level(), levelPatchPreconditioners_[i]->fineLevel());
+            mgTransfer_[i]->setup(contactNetwork_, levelPatchPreconditioners_[i+1]->level(), levelPatchPreconditioners_[i+1]->fineLevel());
         }
 
         //TODO
-        mgTransfer_.push_back(std::make_shared<MGTransfer>(coarseGlobalPreconditioner_->basis(), *allFaultBasis_));
+        //mgTransfer_.push_back(std::make_shared<MGTransfer>(coarseGlobalPreconditioner_->basis(), *allFaultBasis_));
 
         // Remove any recompute filed so that initially the full transferoperator is assembled
         for (size_t i=0; i<mgTransfer_.size(); i++)
@@ -106,49 +105,51 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
         // specify the subset of entries to be reassembled at each iteration
         recompute_.resize(size());
         recompute_[recompute_.size()-1] = contactNetwork_.nBodyAssembler().totalHasObstacle_;
-        for (size_t i=mgTransfer_.size(); i>=1; i--)
-            mgTransfer_[i-1]->restrict(recompute_[i], recompute_[i-1]);
-
-        for (size_t i=0; i<this->mgTransfer_.size(); i++)
+        for (int i=mgTransfer_.size()-1; i>=0; i--) {
+            std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(mgTransfer_[i])->restrict(recompute_[i+1], recompute_[i]);
+          //  std::dynamic_pointer_cast<DenseMultigridTransfer<VectorType, BitVector, MatrixType> >(mgTransfer_[i])->restrict(recompute_[i+1], recompute_[i]);
             std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(mgTransfer_[i])->setRecomputeBitField(&recompute_[i]);
-
+        }
     }
 
 
-   void setIgnore(const BitVector& i) override {
-       this->setIgnore(i);
+   void setIgnore(const BitVector& ignore) {
+       this->setIgnore(ignore);
 
-       auto mgTransfer = std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(this->mgTransfer_[level-1]);
-             mgTransfer->setCriticalBitField(&critical);
+       truncation_.resize(size());
+       truncation_[size()-1] = ignore;
+       for (int i=mgTransfer_.size()-1; i>=0; i--) {
+           std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(mgTransfer_[i])->restrict(truncation_[i+1], truncation_[i]);
+           std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(mgTransfer_[i])->setCriticalBitField(&truncation_[i]);
+       }
    }
 
    void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs) override {
         this->setProblem(mat, x, rhs);
 
-        for (size_t i=0; i<levelPatchPreconditioners_.size(); i++) {
+        size_t maxLevel = levelPatchPreconditioners_.size()-1;
+        levelX_[maxLevel] = x;
+        levelRhs_[maxLevel] = rhs;
+        levelPatchPreconditioners_[maxLevel]->setProblem(mat, levelX_[maxLevel], levelRhs_[maxLevel]);
+
+        for (int i=levelPatchPreconditioners_.size()-2; i>0; i--) {
             const auto& transfer = *mgTransfer_[i];
-            auto& _x = levelX_[i];
-            auto& _rhs = levelRhs_[i];
-            auto& _mat = levelMat_[i];
 
-            transfer.restrict(x, _x);
-            transfer.restrict(rhs, _rhs);
-            transfer.galerkinRestrict(mat, _mat);
+            transfer.restrict(levelX_[i+1], levelX_[i]);
+            transfer.restrict(levelRhs_[i+1], levelRhs_[i]);
+            transfer.galerkinRestrict(*levelPatchPreconditioners_[i+1]->getMatrix(), levelMat_[i]);
 
-            levelPatchPreconditioners_[i]->setProblem(_mat, _x, _rhs);
+            levelPatchPreconditioners_[i]->setProblem(levelMat_[i], levelX_[i], levelRhs_[i]);
         }
 
         // set coarse global problem
-        const auto& coarseTransfer = *mgTransfer_.back();
-        auto& coarseX = levelX_.back();
-        auto& coarseRhs = levelRhs_.back();
-        auto& coarseMat = levelMat_.back();
+        const auto& coarseTransfer = *mgTransfer_[0];
 
-        coarseTransfer.restrict(x, coarseX);
-        coarseTransfer.restrict(rhs, coarseRhs);
-        coarseTransfer.galerkinRestrict(mat, coarseMat);
+        coarseTransfer.restrict(levelX_[1], levelX_[0]);
+        coarseTransfer.restrict(levelRhs_[1], levelRhs_[0]);
+        coarseTransfer.galerkinRestrict(levelMat_[1], levelMat_[0]);
 
-        coarseSolver_.getIterationStep().setProblem(coarseMat, coarseX, coarseRhs);
+        coarseSolver_.getIterationStep().setProblem(levelMat_[0], levelX_[0], levelRhs_[0]);
     }
 
 
@@ -237,19 +238,18 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
     }
 
     void build() {
-        coarseGlobalPreconditioner_->build();
-
         for (size_t i=0; i<levelPatchPreconditioners_.size(); i++) {
             levelPatchPreconditioners_[i]->build();
         }
     }
 
-    void setPatchSolver(PatchSolver&& patchSolver) {
+    void setPatchSolver(PatchSolver&& patchSolver) { //TODO create copies of patch solver for each level
         for (size_t i=0; i<levelPatchPreconditioners_.size(); i++) {
             levelPatchPreconditioners_[i]->setPatchSolver(std::forward<PatchSolver>(patchSolver));
         }
-        coarseGlobalPreconditioner_->setPatchSolver(std::forward<PatchSolver>(patchSolver));
-    };
+
+        coarseSolver_ = Dune::Solvers::wrap_own_share<PatchSolver>(std::forward<PatchSolver>(patchSolver));
+    }
 
     void setPatchDepth(const size_t patchDepth = 0) {
         for (size_t i=0; i<levelPatchPreconditioners_.size(); i++) {
@@ -258,7 +258,7 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
     }
 
     size_t size() const {
-        return levelPatchPreconditioners_.size()+1;
+        return levelPatchPreconditioners_.size();
     }
 };
 

src/spatial-solving/preconditioners/nbodycontacttransfer.cc

 #include <dune/common/fvector.hh>
 #include <dune/common/exceptions.hh>
 #include <dune/istl/bdmatrix.hh>
+#include <dune/istl/matrixindexset.hh>
 
 #include <dune/fufem/assemblers/transferoperatorassembler.hh>
 #include <dune/matrix-vector/blockmatrixview.hh>
@@ -41,7 +42,7 @@ void NBodyContactTransfer<ContactNetwork, VectorType>::setup(const ContactNetwor
 
         } else {
             // gridTransfer is identity if coarse and fine level coincide for body
-            BDMatrix<MatrixBlock>* newMatrix = new BDMatrix<MatrixBlock>(coarseContactNetwork.body(i)->nVertices());
+            Dune::BDMatrix<MatrixBlock>* newMatrix = new Dune::BDMatrix<MatrixBlock>(coarseContactNetwork.body(i)->nVertices());
 
             for (size_t j=0; j<newMatrix->N(); j++)
                 for (int k=0; k<blocksize; k++)
@@ -75,7 +76,7 @@ void NBodyContactTransfer<ContactNetwork, VectorType>::setup(const ContactNetwor
         const auto& contactCouplings = nBodyAssembler.getContactCouplings();
 
         const std::vector<const MatrixType*> mortarTransferOperators(nBodyAssembler.nCouplings());
-        const std::vector<const BitSetVector<1>*> fineHasObstacle(nBodyAssembler.nCouplings());
+        const std::vector<const Dune::BitSetVector<1>*> fineHasObstacle(nBodyAssembler.nCouplings());
         const std::vector<std::array<int,2> > gridIdx(nBodyAssembler.nCouplings());
 
         for (size_t i=0; i<nBodyAssembler.nCouplings(); i++) {
@@ -173,11 +174,11 @@ void NBodyContactTransfer<VectorType>::setupHierarchy(std::vector<std::shared_pt
 }
 */
 
-template<class VectorType>
-void NBodyContactTransfer<VectorType>::combineSubMatrices(const std::vector<const MatrixType*>& submat,
+template<class ContactNetwork, class VectorType>
+void NBodyContactTransfer<ContactNetwork, VectorType>::combineSubMatrices(const std::vector<const MatrixType*>& submat,
                                                                       const std::vector<const MatrixType*>& mortarTransferOperator,
                                                                       const CoordSystemVector& fineLocalCoordSystems,
-                                                                      const std::vector<const BitSetVector<1>*>& fineHasObstacle,
+                                                                      const std::vector<const Dune::BitSetVector<1>*>& fineHasObstacle,
                                                                       const std::vector<std::array<int,2> >& gridIdx)
 {
     // ///////////////////////////////////
@@ -189,7 +190,7 @@ void NBodyContactTransfer<VectorType>::combineSubMatrices(const std::vector<cons
 
     Dune::MatrixVector::BlockMatrixView<MatrixType> view(submat);
 
-    MatrixIndexSet totalIndexSet(view.nRows(), view.nCols());
+    Dune::MatrixIndexSet totalIndexSet(view.nRows(), view.nCols());
 
     // import indices of canonical transfer operator
     for (size_t i=0; i<nGrids; i++)

src/spatial-solving/preconditioners/nbodycontacttransfer.hh

@@ -58,12 +58,8 @@ class NBodyContactTransfer : public TruncatedDenseMGTransfer<VectorType> {
      *  \param fineHasObstacle  Bitfields determining for each coupling which fine grid nodes belong to the nonmortar boundary.
      *  \param gridIdx  For each coupling store the indices of the nonmortar and mortar grid.
     */
-    template <class GridType>
-    void setup(const std::vector<const GridType*>& grids, int colevel,
-               const std::vector<const MatrixType*>& mortarTransferOper,
-               const CoordSystemVector& fineLocalCoordSystems,
-               const std::vector<const BitSetVector<1>*>& fineHasObstacle,
-               const std::vector<std::array<int,2> >& gridIdx);
+    void setup(const ContactNetwork& contactNetwork, const int coarseLevel, const int fineLevel);
+
 
 protected:
     /** \brief Combine all the standard prolongation and mortar operators to one big matrix.
@@ -78,7 +74,7 @@ class NBodyContactTransfer : public TruncatedDenseMGTransfer<VectorType> {
     void combineSubMatrices(const std::vector<const MatrixType*>& submat,
                             const std::vector<const MatrixType*>& mortarTransferOperator,
                             const CoordSystemVector& fineLocalCoordSystems,
-                            const std::vector<const BitSetVector<1>*>& fineHasObstacle,
+                            const std::vector<const Dune::BitSetVector<1>*>& fineHasObstacle,
                             const std::vector<std::array<int,2> >& gridIdx);
 
 };

src/spatial-solving/solverfactory.cc

@@ -11,22 +11,18 @@
 #include "../utils/debugutils.hh"
 
 template <class Functional, class BitVector>
-template <class Preconditioner>
+template <class LinearSolver>
 SolverFactory<Functional, BitVector>::SolverFactory(
     const Dune::ParameterTree& parset,
     Functional& J,
-    Preconditioner&& preconditioner,
+    LinearSolver&& linearSolver,
     const BitVector& ignoreNodes) :
         J_(Dune::Solvers::wrap_own_share<const Functional>(std::forward<Functional>(J))) {
 
     nonlinearSmoother_ = std::make_shared<NonlinearSmoother>(*J_, dummyIterate_, LocalSolver());
 
-    linearSolverStep_ = std::make_shared<LinearSolverStep>();
-    linearSolverStep_.setPreconditioner(std::forward<Preconditioner>(preconditioner));
-    energyNorm_ = std::make_shared<EnergyNorm<Matrix, Vector>>(*linearSolverStep_);
-    linearSolver_ = std::make_shared<LinearSolver>(*linearSolverStep_, parset.get<size_t>("linear.maximumIterations"), parset.get<double>("linear.tolerance"), *energyNorm_, Solver::QUIET);
-
-    tnnmgStep_ = std::make_shared<Step>(*J_, dummyIterate_, nonlinearSmoother_, linearSolver_, DefectProjection(), LineSearchSolver());
+    auto linearSolver = Dune::Solvers::wrap_own_share<LinearSolver>(std::forward<LinearSolver>(linearSolver));
+    tnnmgStep_ = std::make_shared<Step>(*J_, dummyIterate_, nonlinearSmoother_, linearSolver, DefectProjection(), LineSearchSolver());
     tnnmgStep_->setPreSmoothingSteps(parset.get<int>("main.pre"));
     tnnmgStep_->setIgnore(ignoreNodes);
 }

src/spatial-solving/solverfactory.hh

@@ -13,7 +13,6 @@
 //#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
 #include <dune/tnnmg/iterationsteps/nonlineargsstep.hh>
 #include <dune/tnnmg/projections/obstacledefectprojection.hh>
-#include <dune/tnnmg/linearsolvers/fastamgsolver.hh>
 
 #include "tnnmg/tnnmgstep.hh"
 #include "tnnmg/linearization.hh"
@@ -35,13 +34,11 @@ class SolverFactory {
 
     using Step = typename Dune::TNNMG::TNNMGStep<Functional, BitVector, Linearization, DefectProjection, LineSearchSolver>;
 
-    using LinearSolverStep = typename Dune::Solvers::CGStep<Matrix, Vector, BitVector>;
-    using LinearSolver = typename Dune::Solvers::LoopSolver<Vector, BitVector>;
 
-  template <class Preconditioner>
+  template <class LinearSolver>
   SolverFactory(Dune::ParameterTree const &,
                 Functional&,
-                Preconditioner&& preconditioner,
+                LinearSolver&&,
                 const BitVector&);
 
   void setProblem(Vector& x);
@@ -55,11 +52,6 @@ class SolverFactory {
   // nonlinear smoother
   std::shared_ptr<NonlinearSmoother> nonlinearSmoother_;
 
-  // linear solver
-  std::shared_ptr<LinearSolverStep> linearSolverStep_;
-  std::shared_ptr<EnergyNorm<Matrix, Vector>> energyNorm_;
-  std::shared_ptr<LinearSolver> linearSolver_;
-
   // TNNMGStep
   std::shared_ptr<Step> tnnmgStep_;
 };

src/spatial-solving/solverfactory_tmpl.cc

@@ -5,6 +5,8 @@
 #include "../explicitgrid.hh"
 #include "../explicitvectors.hh"
 
+#include <dune/solvers/solvers/loopsolver.hh>
+
 #include "../../dune/tectonic/globalfriction.hh"
 #include "tnnmg/functional.hh"
 #include "tnnmg/zerononlinearity.hh"
@@ -15,5 +17,16 @@ using MyFunctional = Functional<Matrix&, Vector&, MyGlobalFriction&, Vector&, Ve
 using MyZeroFunctional = Functional<Matrix&, Vector&, ZeroNonlinearity&, Vector&, Vector&, double>;
 using MySolverFactory = SolverFactory<MyFunctional, BitVector>;
 
+using MyLinearSolver = Dune::Solvers::LoopSolver<Vector, BitVector>;
+
 template class SolverFactory<MyFunctional, BitVector>;
+template SolverFactory<MyFunctional, BitVector>::SolverFactory(Dune::ParameterTree const &,
+                                                               MyFunctional&,
+                                                               MyLinearSolver&&,
+                                                               const BitVector&);
+
 template class SolverFactory<MyZeroFunctional, BitVector>;
+template SolverFactory<MyZeroFunctional, BitVector>::SolverFactory(Dune::ParameterTree const &,
+                                                               MyZeroFunctional&,
+                                                               MyLinearSolver&&,
+                                                               const BitVector&);