.

src/data-structures/program_state.hh

@@ -25,6 +25,12 @@
 #include "../spatial-solving/tnnmg/zerononlinearity.hh"
 #include "../utils/debugutils.hh"
 
+
+#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
+#include <dune/solvers/iterationsteps/multigridstep.hh>
+
+#include "../spatial-solving/preconditioners/nbodycontacttransfer.hh"
+
 template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class BodyState {
 public:
   using Vector = VectorTEMPLATE;
@@ -98,8 +104,17 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
 
     const auto& preconditionerParset = parset.sub("solver.tnnmg.linear.preconditioner");
 
+    std::cout << "-- setupInitialConditions --" << std::endl;
+    std::cout << "----------------------------" << std::endl;
+
+    /*
+    std::cout << "Building preconditioner..." << std::endl;
     using Preconditioner = MultilevelPatchPreconditioner<ContactNetwork, Matrix, Vector>;
     Dune::BitSetVector<1> activeLevels(contactNetwork.nLevels(), true);
+    /*Dune::BitSetVector<1> activeLevels(contactNetwork.nLevels(), false);
+    activeLevels.back() = true;
+    activeLevels[activeLevels.size()-2] = true;*/
+    /*print(activeLevels, "activeLevels:");
     Preconditioner preconditioner(preconditionerParset, contactNetwork, activeLevels);
     preconditioner.build();
 
@@ -110,7 +125,31 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
     linearSolverStep.setPreconditioner(preconditioner);
 
     EnergyNorm<Matrix, Vector> energyNorm(linearSolverStep);
-    LinearSolver linearSolver(linearSolverStep, parset.get<size_t>("solver.tnnmg.linear.maximumIterations"), parset.get<double>("solver.tnnmg.linear.tolerance"), energyNorm, Solver::QUIET);
+    LinearSolver linearSolver(linearSolverStep, parset.get<size_t>("solver.tnnmg.linear.maximumIterations"), parset.get<double>("solver.tnnmg.linear.tolerance"), energyNorm, Solver::FULL);
+    */
+    using LinearSolver = typename Dune::Solvers::LoopSolver<Vector, BitVector>;
+    using TransferOperator = NBodyContactTransfer<ContactNetwork, Vector>;
+    using TransferOperators = std::vector<std::shared_ptr<TransferOperator>>;
+
+    TransferOperators transfer(contactNetwork.nLevels()-1);
+    for (size_t i=0; i<transfer.size(); i++) {
+        transfer[i] = std::make_shared<TransferOperator>();
+        transfer[i]->setup(contactNetwork, i, i+1);
+    }
+
+    // Remove any recompute filed so that initially the full transferoperator is assembled
+    for (size_t i=0; i<transfer.size(); i++)
+        std::dynamic_pointer_cast<TruncatedMGTransfer<Vector> >(transfer[i])->setRecomputeBitField(nullptr);
+
+    auto smoother = TruncatedBlockGSStep<Matrix, Vector>{};
+    auto linearMultigridStep = std::make_shared<Dune::Solvers::MultigridStep<Matrix, Vector> >();
+    linearMultigridStep->setMGType(1, 3, 3);
+    linearMultigridStep->setSmoother(&smoother);
+    linearMultigridStep->setTransferOperators(transfer);
+
+    EnergyNorm<Matrix, Vector> mgNorm(*linearMultigridStep);
+    LinearSolver mgSolver(linearMultigridStep, parset.get<size_t>("solver.tnnmg.linear.maximumIterations"), parset.get<double>("solver.tnnmg.linear.tolerance"), mgNorm, Solver::FULL);
+
 
     // Solving a linear problem with a multigrid solver
     auto const solveLinearProblem = [&](
@@ -177,7 +216,8 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
 
       // set up TNMMG solver
       using Factory = SolverFactory<Functional, BitVector>;
-      Factory factory(parset.sub("solver.tnnmg"), J, linearSolver, _dirichletNodes);
+      //Factory factory(parset.sub("solver.tnnmg"), J, linearSolver, _dirichletNodes);
+      Factory factory(parset.sub("solver.tnnmg"), J, mgSolver, _dirichletNodes);
 
      /* std::vector<BitVector> bodyDirichletNodes;
       nBodyAssembler.postprocess(_dirichletNodes, bodyDirichletNodes);
@@ -270,10 +310,12 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
         }
     }*/
 
+    std::cout << "solving linear problem for u..." << std::endl;
+
     solveLinearProblem(dirichletNodes, contactNetwork.matrices().elasticity, ell0, u,
                        parset.sub("u0.solver"));
 
-    print(u, "initial u:");
+    //print(u, "initial u:");
 
     // Initial acceleration: Computed in agreement with Ma = ell0 - Au
     // (without Dirichlet constraints), again assuming dPhi(v = 0) = 0
@@ -300,11 +342,13 @@ template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
       Dune::MatrixVector::subtractProduct(accelerationRHS[i], *body->matrices().elasticity, u[i]);
     }
 
+    std::cout << "solving linear problem for a..." << std::endl;
+
     BitVector noNodes(dirichletNodes.size(), false);
     solveLinearProblem(noNodes, contactNetwork.matrices().mass, accelerationRHS, a,
                        parset.sub("a0.solver"));
 
-    print(a, "initial a:");
+    //print(a, "initial a:");
   }
 
 private:

src/multi-body-problem-2D.cfg

@@ -21,4 +21,4 @@ tolerance         = 1e-5
 tolerance          = 1e-10
 
 [solver.tnnmg.linear.preconditioner]
-tolerance          = 1e-10
+tolerance          = 1e-8

src/multi-body-problem.cc

@@ -190,7 +190,7 @@ int main(int argc, char *argv[]) {
     }
     contactNetwork.setDeformation(def);*/
 
-    const auto & coarseContactNetwork = *contactNetwork.level(3);
+   /* const auto & coarseContactNetwork = *contactNetwork.level(3);
     const auto & fineContactNetwork = *contactNetwork.level(4);
    // SupportPatchFactory<typename ContactNetwork::LevelContactNetwork> supportPatchFactory(coarseContactNetwork, fineContactNetwork);
     size_t bodyID = 2;
@@ -222,7 +222,7 @@ int main(int argc, char *argv[]) {
             dofs[j] = j;
         }
         writeToVTK(gv, dofs, "", "body_" + std::to_string(i) + "_fine");
-    }
+    } */
 
     /*
     using Patch = typename SupportPatchFactory<typename ContactNetwork::LevelContactNetwork>::Patch;
@@ -267,9 +267,6 @@ int main(int argc, char *argv[]) {
         }
     } */
 
-    std::cout << "Done! :)" << std::endl;
-    //return 1;
-
     //printMortarBasis<Vector>(contactNetwork.nBodyAssembler());
 
     // -----------------
@@ -417,7 +414,7 @@ int main(int argc, char *argv[]) {
         }
     }*/
 
-    print(totalDirichletNodes, "totalDirichletNodes:");
+    //print(totalDirichletNodes, "totalDirichletNodes:");
 
     using Functional = Functional<Matrix&, Vector&, GlobalFriction<Matrix, Vector>&, Vector&, Vector&, field_type>;
     using NonlinearFactory = SolverFactory<Functional, BitVector>;
@@ -433,18 +430,18 @@ int main(int argc, char *argv[]) {
     BoundaryNodes dirichletNodes;
     contactNetwork.boundaryNodes("dirichlet", dirichletNodes);
 
-    for (size_t i=0; i<dirichletNodes.size(); i++) {
+    /*for (size_t i=0; i<dirichletNodes.size(); i++) {
         for (size_t j=0; j<dirichletNodes[i].size(); j++) {
         print(*dirichletNodes[i][j], "dirichletNodes_body_" + std::to_string(i) + "_boundary_" + std::to_string(j));
         }
-    }
+    }*/
 
     std::vector<const Dune::BitSetVector<1>*> frictionNodes;
     contactNetwork.frictionNodes(frictionNodes);
 
-    for (size_t i=0; i<frictionNodes.size(); i++) {
+    /*for (size_t i=0; i<frictionNodes.size(); i++) {
         print(*frictionNodes[i], "frictionNodes_body_" + std::to_string(i));
-    }
+    }*/
 
     Updaters current(
         initRateUpdater(
@@ -527,21 +524,21 @@ int main(int argc, char *argv[]) {
     Preconditioner preconditioner(contactNetwork, activeLevels, preconditionerParset.get<Preconditioner::Mode>("mode"));
     preconditioner.setPatchSolver(patchSolver);
     preconditioner.setPatchDepth(preconditionerParset.get<size_t>("patchDepth"));
-
+*/
     // set adaptive time stepper
     typename ContactNetwork::ExternalForces externalForces;
     contactNetwork.externalForces(externalForces);
 
-    AdaptiveTimeStepper<NonlinearFactory, std::decay_t<decltype(nBodyAssembler)>, Updaters, std::decay_t<decltype(stateEnergyNorms)>>
-        adaptiveTimeStepper(parset, nBodyAssembler, totalDirichletNodes, globalFriction, frictionNodes, current,
+    AdaptiveTimeStepper<NonlinearFactory, std::decay_t<decltype(contactNetwork)>, Updaters, std::decay_t<decltype(stateEnergyNorms)>>
+        adaptiveTimeStepper(parset, contactNetwork, totalDirichletNodes, globalFriction, frictionNodes, current,
                             programState.relativeTime, programState.relativeTau,
                             externalForces, stateEnergyNorms, mustRefine);
 
     while (!adaptiveTimeStepper.reachedEnd()) {
       programState.timeStep++;
 
-      preconditioner.build();
-      iterationCount = adaptiveTimeStepper.advance(preconditioner);
+      //preconditioner.build();
+      iterationCount = adaptiveTimeStepper.advance();
 
       programState.relativeTime = adaptiveTimeStepper.relativeTime_;
       programState.relativeTau = adaptiveTimeStepper.relativeTau_;
@@ -559,6 +556,8 @@ int main(int argc, char *argv[]) {
 
       report();
 
+      break;
+
       if (programState.timeStep==50) {
         std::cout << "limit of timeSteps reached!" << std::endl;
         break; // TODO remove after debugging
@@ -569,7 +568,7 @@ int main(int argc, char *argv[]) {
         break;
       }
     }
-    */
+
 
     std::cout.rdbuf(coutbuf); //reset to standard output again
 

src/multi-body-problem.cfg

@@ -11,7 +11,7 @@ vtk.write       = true
 
 [problem]
 finalTime       = 100  # [s] #1000
-bodyCount       = 4
+bodyCount       = 2
 
 [body]
 bulkModulus     = 0.5e5 # [Pa]
@@ -49,11 +49,11 @@ relativeTau = 1e-4 # 1e-6
 scheme = newmark
 
 [u0.solver]
-maximumIterations = 10000
+maximumIterations = 100
 verbosity         = full
 
 [a0.solver]
-maximumIterations = 10000
+maximumIterations = 100
 verbosity         = full
 
 [v.solver]

src/spatial-solving/fixedpointiterator.cc

@@ -30,6 +30,11 @@
 #include "../utils/tobool.hh"
 #include "../utils/debugutils.hh"
 
+#include <dune/solvers/solvers/loopsolver.hh>
+#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
+#include <dune/solvers/iterationsteps/multigridstep.hh>
+
+#include "../spatial-solving/preconditioners/nbodycontacttransfer.hh"
 
 void FixedPointIterationCounter::operator+=(
     FixedPointIterationCounter const &other) {
@@ -37,16 +42,16 @@ void FixedPointIterationCounter::operator+=(
   multigridIterations += other.multigridIterations;
 }
 
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
-FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::FixedPointIterator(
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
+FixedPointIterator<Factory, ContactNetwork, Updaters, ErrorNorms>::FixedPointIterator(
     Dune::ParameterTree const &parset,
-    NBodyAssembler& nBodyAssembler,
+    const ContactNetwork& contactNetwork,
     const IgnoreVector& ignoreNodes,
     GlobalFriction& globalFriction,
     const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
     const ErrorNorms& errorNorms)
     : parset_(parset),
-      nBodyAssembler_(nBodyAssembler),
+      contactNetwork_(contactNetwork),
       ignoreNodes_(ignoreNodes),
       globalFriction_(globalFriction),
       bodywiseNonmortarBoundaries_(bodywiseNonmortarBoundaries),
@@ -58,17 +63,18 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::FixedPointIte
       verbosity_(parset.get<Solver::VerbosityMode>("v.solver.verbosity")),
       errorNorms_(errorNorms) {}
 
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
-template <class Preconditioner>
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
 FixedPointIterationCounter
-FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::run(
-    Updaters updaters, Preconditioner&& preconditioner,
+FixedPointIterator<Factory, ContactNetwork, Updaters, ErrorNorms>::run(
+    Updaters updaters,
     const std::vector<Matrix>& velocityMatrices, const std::vector<Vector>& velocityRHSs,
     std::vector<Vector>& velocityIterates) {
 
   std::cout << "FixedPointIterator::run()" << std::endl;
 
-  const auto nBodies = nBodyAssembler_.nGrids();
+  const auto& nBodyAssembler = contactNetwork_.nBodyAssembler();
+
+  const auto nBodies = nBodyAssembler.nGrids();
 
   std::vector<const Matrix*> matrices_ptr(nBodies);
   for (int i=0; i<nBodies; i++) {
@@ -77,17 +83,17 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::run(
 
   // assemble full global contact problem
   Matrix bilinearForm;
-  nBodyAssembler_.assembleJacobian(matrices_ptr, bilinearForm);
+  nBodyAssembler.assembleJacobian(matrices_ptr, bilinearForm);
 
   //print(bilinearForm, "bilinearForm:");
 
   Vector totalRhs;
-  nBodyAssembler_.assembleRightHandSide(velocityRHSs, totalRhs);
+  nBodyAssembler.assembleRightHandSide(velocityRHSs, totalRhs);
 
   //print(totalRhs, "totalRhs:");
 
   // get lower and upper obstacles
-  const auto& totalObstacles = nBodyAssembler_.totalObstacles_;
+  const auto& totalObstacles = nBodyAssembler.totalObstacles_;
   Vector lower(totalObstacles.size());
   Vector upper(totalObstacles.size());
 
@@ -101,10 +107,10 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::run(
     }
   }
 
-  print(totalObstacles, "totalObstacles:");
+  //print(totalObstacles, "totalObstacles:");
 
-  print(lower, "lower obstacles:");
-  print(upper, "upper obstacles:");
+  //print(lower, "lower obstacles:");
+  //print(upper, "upper obstacles:");
 
   // comput velocity obstacles
   Vector vLower, vUpper;
@@ -113,20 +119,44 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::run(
   updaters.rate_->extractOldDisplacement(u0);
 
   Vector totalu0, totalv0;
-  nBodyAssembler_.concatenateVectors(u0, totalu0);
-  nBodyAssembler_.concatenateVectors(v0, totalv0);
+  nBodyAssembler.concatenateVectors(u0, totalu0);
+  nBodyAssembler.concatenateVectors(v0, totalv0);
 
   updaters.rate_->velocityObstacles(totalu0, lower, totalv0, vLower);
   updaters.rate_->velocityObstacles(totalu0, upper, totalv0, vUpper);
 
-  print(vLower, "vLower obstacles:");
-  print(vUpper, "vUpper obstacles:");
+  //print(vLower, "vLower obstacles:");
+  //print(vUpper, "vUpper obstacles:");
 
   std::cout << "- Problem assembled: success" << std::endl;
 
+  using LinearSolver = typename Dune::Solvers::LoopSolver<Vector, IgnoreVector>;
+  using TransferOperator = NBodyContactTransfer<ContactNetwork, Vector>;
+  using TransferOperators = std::vector<std::shared_ptr<TransferOperator>>;
+
+  TransferOperators transfer(contactNetwork_.nLevels()-1);
+  for (size_t i=0; i<transfer.size(); i++) {
+      transfer[i] = std::make_shared<TransferOperator>();
+      transfer[i]->setup(contactNetwork_, i, i+1);
+  }
+
+  // Remove any recompute filed so that initially the full transferoperator is assembled
+  for (size_t i=0; i<transfer.size(); i++)
+      std::dynamic_pointer_cast<TruncatedMGTransfer<Vector> >(transfer[i])->setRecomputeBitField(nullptr);
+
+  auto smoother = TruncatedBlockGSStep<Matrix, Vector>{};
+  auto linearMultigridStep = std::make_shared<Dune::Solvers::MultigridStep<Matrix, Vector> >();
+  linearMultigridStep->setMGType(1, 3, 3);
+  linearMultigridStep->setSmoother(&smoother);
+  linearMultigridStep->setTransferOperators(transfer);
+
+  EnergyNorm<Matrix, Vector> mgNorm(*linearMultigridStep);
+  LinearSolver mgSolver(linearMultigridStep, parset_.get<size_t>("solver.tnnmg.linear.maximumIterations"), parset_.get<double>("solver.tnnmg.linear.tolerance"), mgNorm, Solver::QUIET);
+
+
   // set up functional and TNMMG solver
   Functional J(bilinearForm, totalRhs, globalFriction_, vLower, vUpper);
-  Factory solverFactory(parset_.sub("solver.tnnmg"), J, std::forward<Preconditioner>(preconditioner), ignoreNodes_);
+  Factory solverFactory(parset_.sub("solver.tnnmg"), J, mgSolver, ignoreNodes_);
   auto step = solverFactory.step();
 
   std::cout << "- Functional and TNNMG step setup: success" << std::endl;
@@ -147,7 +177,7 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::run(
     globalFriction_.updateAlpha(alpha);
 
     Vector totalVelocityIterate;
-    nBodyAssembler_.nodalToTransformed(velocityIterates, totalVelocityIterate);
+    nBodyAssembler.nodalToTransformed(velocityIterates, totalVelocityIterate);
 
     //print(velocityIterates, "velocityIterates:");
     //print(totalVelocityIterate, "totalVelocityIterate:");
@@ -164,8 +194,8 @@ FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::run(
     std::cout << "-- Solved" << std::endl;
 
     const auto& tnnmgSol = step->getSol();
-    nBodyAssembler_.postprocess(tnnmgSol, velocityIterates);
-    //nBodyAssembler_.postprocess(totalVelocityIterate, velocityIterates);
+    nBodyAssembler.postprocess(tnnmgSol, velocityIterates);
+    //nBodyAssembler.postprocess(totalVelocityIterate, velocityIterates);
 
     //print(totalVelocityIterate, "totalVelocityIterate:");
     //print(velocityIterates, "velocityIterates:");
@@ -239,13 +269,14 @@ std::ostream &operator<<(std::ostream &stream,
                 << ")";
 }
 
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
-void FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::relativeVelocities(
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
+void FixedPointIterator<Factory, ContactNetwork, Updaters, ErrorNorms>::relativeVelocities(
     const Vector& v,
     std::vector<Vector>& v_rel) const {
 
-    const size_t nBodies = nBodyAssembler_.nGrids();
-   // const auto& contactCouplings = nBodyAssembler_.getContactCouplings();
+    const auto& nBodyAssembler = contactNetwork_.nBodyAssembler();
+    const size_t nBodies = nBodyAssembler.nGrids();
+   // const auto& contactCouplings = nBodyAssembler.getContactCouplings();
 
     size_t globalIdx = 0;
     v_rel.resize(nBodies);

src/spatial-solving/fixedpointiterator.hh

@@ -23,7 +23,7 @@ struct FixedPointIterationCounter {
 std::ostream &operator<<(std::ostream &stream,
                          FixedPointIterationCounter const &fpic);
 
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
 class FixedPointIterator {
   using Functional = typename Factory::Functional;
   using ScalarVector = typename Updaters::StateUpdater::ScalarVector;
@@ -46,22 +46,20 @@ class FixedPointIterator {
 
 public:
   FixedPointIterator(const Dune::ParameterTree& parset,
-                     NBodyAssembler& nBodyAssembler,
+                     const ContactNetwork& contactNetwork,
                      const IgnoreVector& ignoreNodes,
                      GlobalFriction& globalFriction,
                      const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
                      const ErrorNorms& errorNorms);
 
-  template <class Preconditioner>
   FixedPointIterationCounter run(Updaters updaters,
-                                 Preconditioner&& preconditioner,
                                  const std::vector<Matrix>& velocityMatrices,
                                  const std::vector<Vector>& velocityRHSs,
                                  std::vector<Vector>& velocityIterates);
 
 private:
   const Dune::ParameterTree& parset_;
-  NBodyAssembler& nBodyAssembler_;
+  const ContactNetwork& contactNetwork_;
   const IgnoreVector& ignoreNodes_;
 
   GlobalFriction& globalFriction_;

src/spatial-solving/fixedpointiterator_tmpl.cc

@@ -14,7 +14,6 @@
 #include "../time-stepping/state/stateupdater.hh"
 #include "../time-stepping/updaters.hh"
 
-using NBodyAssembler = typename MyContactNetwork::NBodyAssembler;
 using BoundaryNodes = typename MyContactNetwork::BoundaryNodes;
 using BoundaryFunctions = typename MyContactNetwork::BoundaryFunctions;
 
@@ -25,4 +24,4 @@ using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
 using ErrorNorms = typename MyContactNetwork::StateEnergyNorms;
 
 
-template class FixedPointIterator<MySolverFactory, NBodyAssembler, MyUpdaters, ErrorNorms>;
+template class FixedPointIterator<MySolverFactory, MyContactNetwork, MyUpdaters, ErrorNorms>;

src/spatial-solving/preconditioners/levelpatchpreconditioner.hh

@@ -13,6 +13,7 @@
 #include "../../data-structures/levelcontactnetwork.hh"
 
 #include "supportpatchfactory.hh"
+#include "../tnnmg/localproblem.hh"
 
 #include <dune/localfunctions/lagrange/pqkfactory.hh>
 
@@ -108,6 +109,8 @@ class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorTy
         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());
 
@@ -117,6 +120,27 @@ class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorTy
                     if (!vertexVisited[globalIdx][0]) {
                         vertexVisited[globalIdx][0] = true;
                         patchFactory_.build(bodyIdx, e, i, patches_[globalIdx], patchDepth_);
+
+                        /*print(patches_[globalIdx], "patch:");
+
+                        size_t c = 0;
+                        for (size_t j=0; j<levelContactNetwork_.nBodies(); j++) {
+                            const auto& gv = fineContactNetwork_.body(j)->gridView();
+
+                            printDofLocation(gv);
+
+                            Dune::BlockVector<Dune::FieldVector<ctype, 1>> patchVec(gv.size(dim));
+                            for (size_t l=0; l<patchVec.size(); l++) {
+                                if (patches_[globalIdx][c++][0]) {
+                                    patchVec[l][0] = 1;
+                                }
+                            }
+
+                            //print(patchVec, "patchVec");
+
+                            // output patch
+                            writeToVTK(gv, patchVec, "", "level_" + std::to_string(level_) + "_patch_" + std::to_string(globalIdx) + "_body_" + std::to_string(j));
+                        }*/
                     }
                 }
             }
@@ -160,6 +184,7 @@ class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorTy
         for (const auto& p : patches_) {
             x = 0;
 
+            /*
             auto& step = patchSolver_->getIterationStep();
             dynamic_cast<LinearIterationStep<MatrixType, VectorType>&>(step).setProblem(*this->mat_, x, *this->rhs_);
             step.setIgnore(p);
@@ -168,7 +193,29 @@ class LevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorTy
             patchSolver_->preprocess();
             patchSolver_->solve();
 
-            *(this->x_) += step.getSol();
+            *(this->x_) += step.getSol();*/
+
+            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;
+                }
+            }
+
+            LocalProblem<MatrixType, VectorType> localProblem(*this->mat_, *this->rhs_, ignore);
+            Vector newR;
+            localProblem.getLocalRhs(x, newR);
+
+           /* print(ignore, "ignore:");
+            print(*this->mat_, "Mat:");
+            print(localProblem.getMat(), "localMat:");*/
+
+            patchSolver_->setProblem(localProblem.getMat(), x, newR);
+            patchSolver_->preprocess();
+            patchSolver_->solve();
+
+            *(this->x_) += x;
         }
     }
 

src/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh

@@ -10,20 +10,27 @@
 
 #include <dune/solvers/norms/energynorm.hh>
 #include <dune/solvers/solvers/loopsolver.hh>
+#include <dune/solvers/solvers/umfpacksolver.hh>
 #include <dune/solvers/iterationsteps/blockgssteps.hh>
+//#include <dune/solvers/iterationsteps/blockgsstep.hh>
 #include <dune/solvers/iterationsteps/cgstep.hh>
 #include <dune/solvers/iterationsteps/lineariterationstep.hh>
 
 #include "nbodycontacttransfer.hh"
 #include "levelpatchpreconditioner.hh"
 
+#include "../../utils/debugutils.hh"
+
 template <class ContactNetwork, class MatrixType, class VectorType>
 class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, VectorType> {
 public:
     //struct Mode { enum ModeType {additive, multiplicative}; };
 
 private:
-    using PatchSolver = typename Dune::Solvers::LoopSolver<Vector, BitVector>;
+    using Base = LinearIterationStep<MatrixType, VectorType>;
+
+    //using PatchSolver = typename Dune::Solvers::LoopSolver<Vector, BitVector>;
+    using PatchSolver = Dune::Solvers::UMFPackSolver<MatrixType, VectorType>;
 
     using LevelContactNetwork = typename ContactNetwork::LevelContactNetwork;
     using LevelPatchPreconditioner = LevelPatchPreconditioner<LevelContactNetwork, PatchSolver, MatrixType, VectorType>;
@@ -44,10 +51,14 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
     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<BitVector> truncation_;
+    //std::vector<BitVector> truncation_;
     std::vector<std::shared_ptr<MGTransfer>> mgTransfer_; //TODO
 
 public:
@@ -84,29 +95,36 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
             }
         }
 
+        levelMat_.resize(size()-1);
         levelX_.resize(size());
         levelRhs_.resize(size());
 
         // init patch solvers
         levelSolvers_.resize(size());
-        auto coarseStep = Dune::Solvers::BlockGSStepFactory<MatrixType, VectorType>::create(Dune::Solvers::BlockGS::LocalSolvers::direct(0.0));
-        EnergyNorm<MatrixType, VectorType> coarseEnergyNorm(coarseStep);
-        levelSolvers_[0] = std::make_shared<PatchSolver>(coarseStep,
+        levelItSteps_.resize(size());
+        levelErrorNorms_.resize(size());
+        //auto coarseStep =
+        levelItSteps_[0] = Dune::Solvers::BlockGSStepFactory<MatrixType, VectorType>::createPtr(Dune::Solvers::BlockGS::LocalSolvers::direct(0.0));
+        levelErrorNorms_[0] = std::make_shared<EnergyNorm<MatrixType, VectorType>>(*levelItSteps_[0].get());
+        levelSolvers_[0] = /*std::make_shared<PatchSolver>(*levelItSteps_[0].get(),
                                    parset.get<size_t>("maximumIterations"),
                                    parset.get<double>("tolerance"),
-                                   coarseEnergyNorm,
+                                   *levelErrorNorms_[0].get(),
                                    parset.get<Solver::VerbosityMode>("verbosity"),
-                                   false); // absolute error
+                                   false); // absolute error */
+                std::make_shared<PatchSolver>();
 
         for (size_t i=1; i<levelSolvers_.size(); i++) {
-            auto gsStep = Dune::Solvers::BlockGSStepFactory<MatrixType, VectorType>::create(Dune::Solvers::BlockGS::LocalSolvers::direct(0.0));
-            EnergyNorm<MatrixType, VectorType> energyNorm(gsStep);
-            levelSolvers_[i] = std::make_shared<PatchSolver>(gsStep,
+            //auto gsStep = Dune::Solvers::BlockGSStepFactory<MatrixType, VectorType>::create(Dune::Solvers::BlockGS::LocalSolvers::direct(0.0));
+            levelItSteps_[i] = Dune::Solvers::BlockGSStepFactory<MatrixType, VectorType>::createPtr(Dune::Solvers::BlockGS::LocalSolvers::direct(0.0));
+            levelErrorNorms_[i] = std::make_shared<EnergyNorm<MatrixType, VectorType>>(*levelItSteps_[i].get());
+            levelSolvers_[i] = /*std::make_shared<PatchSolver>(*levelItSteps_[i].get(),
                                        parset.get<size_t>("maximumIterations"),
                                        parset.get<double>("tolerance"),
-                                       energyNorm,
+                                       *levelErrorNorms_[i].get(),
                                        parset.get<Solver::VerbosityMode>("verbosity"),
-                                       false); // absolute error
+                                       false); // absolute error */
+                    std::make_shared<PatchSolver>();
             levelPatchPreconditioners_[i]->setPatchSolver(levelSolvers_[i]);
         }
 
@@ -118,38 +136,91 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
             mgTransfer_[i]->setup(contactNetwork_, levelPatchPreconditioners_[i+1]->level(), levelPatchPreconditioners_[i+1]->fineLevel());
         }
 
-        //TODO
-        //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++)
             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_.resize(size());
         recompute_[recompute_.size()-1] = contactNetwork_.nBodyAssembler().totalHasObstacle_;
         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]);
+            print(recompute_[i], "recompute: ");
+        }*/
     }
+
+   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]);
        }
 
+      /* truncation_.resize(size());
+       truncation_[size()-1] = Base::ignore();
+       print(truncation_.back(), "truncation: ");
+       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+1]);
+           print(truncation_[i], "truncation: ");
+       }*/
+
+        size_t maxLevel = levelPatchPreconditioners_.size()-1;
+        levelPatchPreconditioners_[maxLevel]->setMatrix(mat);
+
+        //print(mat, "levelMat_" + std::to_string(maxLevel) + ":");
+
+        for (int i=levelPatchPreconditioners_.size()-2; i>0; i--) {
+            const auto& transfer = *mgTransfer_[i];
+
+            //print(*levelPatchPreconditioners_[i+1]->getMatrix(), "levelMat:");
+            //print(levelMat_[i], "coarseMat:");
+
+            transfer.galerkinRestrictSetOccupation(*levelPatchPreconditioners_[i+1]->getMatrix(), levelMat_[i]);
+            transfer.galerkinRestrict(*levelPatchPreconditioners_[i+1]->getMatrix(), levelMat_[i]);
+            levelPatchPreconditioners_[i]->setMatrix(levelMat_[i]);
 
-   void setIgnore(const BitVector& ignore) {
-       this->setIgnore(ignore);
+            //print(levelMat_[i], "levelMat_" + std::to_string(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]);
+
+        //print(levelMat_[0], "levelMat_0:");
+
+        /*dynamic_cast<Base&>(levelSolvers_[0]->getIterationStep()).setMatrix(levelMat_[0]); */
+
+        Dune::MatrixVector::resize(levelX_[0], levelMat_[0]);
+         //Dune::MatrixVector::resize(levelRhs_[0], levelMat_[0]);
+   }
+
+   void setTruncation(const BitVector& truncation) {
+       std::dynamic_pointer_cast<TruncatedMGTransfer<VectorType> >(mgTransfer_.back())->setCriticalBitField(&truncation);
+   }
+
+   /*void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs) override {
+       LinearIterationStep<MatrixType, VectorType>::setProblem(mat, x, rhs);
 
        truncation_.resize(size());
-       truncation_[size()-1] = ignore;
+       truncation_[size()-1] = this->ignore();
+       print(truncation_.back(), "truncation: ");
        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]);
-       }
+           print(truncation_[i], "truncation: ");
        }
 
-   void setProblem(const MatrixType& mat, VectorType& x, const VectorType& rhs) override {
-        this->setProblem(mat, x, rhs);
-
         size_t maxLevel = levelPatchPreconditioners_.size()-1;
         levelX_[maxLevel] = x;
         levelRhs_[maxLevel] = rhs;
@@ -176,10 +247,31 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
 
       //  EnergyNorm<MatrixType, VectorType> coarseErrorNorm(levelSolvers_[0]->getIterationStep());
       //  levelSolvers_[0]->setErrorNorm(coarseErrorNorm);
+    } */
+
+    void iterate() {
+        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];
+
+            /*print(levelX_[i+1], "levelX_:");
+            print(levelX_[i], "levelX_:");
+            print(levelRhs_[i+1], "levelRhs_:");
+            print(levelRhs_[i], "levelRhs_:");*/
+
+            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]);
+
 
-    void iterate() {
         if (mode_ == additive)
             iterateAdd();
         else
@@ -187,26 +279,46 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
     }
 
     void iterateAdd() {
+        //print(this->ignore(), "ignore: ");
+
         *(this->x_) = 0;	
 
         VectorType x;
 
+        /*
         // solve coarse global problem
         auto& coarseSolver = levelSolvers_[0];
-        coarseSolver->getIterationStep().setIgnore(truncation_[0]);
+        //coarseSolver->getIterationStep().setIgnore(ignoreHierarchy_[0]);
+
+        dynamic_cast<LinearIterationStep<MatrixType, VectorType>&>(levelSolvers_[0]->getIterationStep()).setProblem(levelMat_[0], levelX_[0], levelRhs_[0]);
 
         coarseSolver->check();
         coarseSolver->preprocess();
         coarseSolver->solve();
 
-        mgTransfer_[0]->prolong(coarseSolver->getIterationStep().getSol(), x);
+        mgTransfer_[0]->prolong(coarseSolver->getIterationStep().getSol(), x); */
 
-        for (size_t i=1; i<levelPatchPreconditioners_.size(); i++) {
+
+        BitVector emptyIgnore0(levelX_[0].size(), false);
+        LocalProblem<MatrixType, VectorType> localProblem(levelMat_[0], levelRhs_[0], ignoreHierarchy_[0]);
+        Vector newR;
+        localProblem.getLocalRhs(levelX_[0], newR);
+
+        auto& coarseSolver = levelSolvers_[0];
+        coarseSolver->setProblem(localProblem.getMat(), levelX_[0], newR);
+        coarseSolver->preprocess();
+        coarseSolver->solve();
+
+        mgTransfer_[0]->prolong(levelX_[0], x);
+
+        //print(x, "MultilevelPatchPreconditioner: x0");
+
+        for (size_t i=1; i<levelPatchPreconditioners_.size()-1; i++) {
             const auto& transfer = *mgTransfer_[i];
             auto& preconditioner = *levelPatchPreconditioners_[i];
 
-            preconditioner.setIgnore(truncation_[i]);
-            preconditioner.iterate();
+            preconditioner.setIgnore(ignoreHierarchy_[i]);
+            preconditioner.apply(levelX_[i], levelRhs_[i]);
 
             x += preconditioner.getSol();
 
@@ -214,8 +326,30 @@ class MultilevelPatchPreconditioner : public LinearIterationStep<MatrixType, Vec
             transfer.prolong(x, newX);
 
             x = newX;
+            //print(x, "MultilevelPatchPreconditioner: x" + std::to_string(i));
         }
 
+        auto& preconditioner = *levelPatchPreconditioners_.back();
+
+        preconditioner.setIgnore(ignoreHierarchy_.back());
+        preconditioner.apply(levelX_.back(), levelRhs_.back());
+
+        x += preconditioner.getSol();
+
+        //print(x, "MultilevelPatchPreconditioner: xmax");
+
+
+        BitVector emptyIgnore(this->x_->size(), false);
+        LocalProblem<MatrixType, VectorType> globalProblem(*this->mat_.get(), *this->rhs_, emptyIgnore);
+        globalProblem.getLocalRhs(*this->x_, newR);
+        VectorType totalX = *(this->x_);
+        auto& globalSolver = levelSolvers_.back();
+        globalSolver->setProblem(globalProblem.getMat(), totalX, newR);
+        globalSolver->preprocess();
+        globalSolver->solve();
+
+        //print(totalX, "MultilevelPatchPreconditioner: totalX");
+
         *(this->x_) = x;
     }
 

src/spatial-solving/preconditioners/supportpatchfactory.hh

@@ -538,7 +538,7 @@ class SupportPatchFactory
                     std::set<size_t> isDofs;
                     intersectionDofs(fineIndices_, is, bodyID, isDofs);
 
-                    addToPatchBoundary(isDofs, patchBoundary, patchDofs);
+                    addToPatch(isDofs, patchBoundary, patchDofs);
                 }
             }
         }

src/spatial-solving/solverfactory.cc

@@ -20,6 +20,9 @@ SolverFactory<Functional, BitVector>::SolverFactory(
     nonlinearSmoother_ = std::make_shared<NonlinearSmoother>(*J_, dummyIterate_, LocalSolver());
 
     auto linearSolver_ptr = Dune::Solvers::wrap_own_share<std::decay_t<LinearSolver>>(std::forward<LinearSolver>(linearSolver));
+
+    //tnnmgStep_ = std::make_shared<Step>(*J_, dummyIterate_, nonlinearSmoother_, linearSolver_ptr, 10, DefectProjection(), LineSearchSolver());
+
     tnnmgStep_ = std::make_shared<Step>(*J_, dummyIterate_, nonlinearSmoother_, linearSolver_ptr, DefectProjection(), LineSearchSolver());
     tnnmgStep_->setPreSmoothingSteps(parset.get<int>("main.pre"));
     tnnmgStep_->setIgnore(ignoreNodes);

src/spatial-solving/tnnmg/functional.hh

@@ -97,8 +97,9 @@ class DirectionalRestriction :
   {
     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;
+    //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_;
     }
@@ -401,7 +402,7 @@ class Functional : public Dune::TNNMG::BoxConstrainedQuadraticFunctional<M, V, L
   Range operator()(const Vector& v) const
   { 
     if (Dune::TNNMG::checkInsideBox(v, this->lower_.get(), this->upper_.get()))
-      return Base::Base::operator()(v) + phi_.get().operator()(v);
+      return Dune::TNNMG::QuadraticFunctional<M,V,R>::operator()(v) + phi_.get().operator()(v);
     return std::numeric_limits<Range>::max();
   }
 

src/spatial-solving/tnnmg/linearcorrection.hh

@@ -9,6 +9,10 @@
 #include <dune/solvers/solvers/linearsolver.hh>
 #include <dune/solvers/common/canignore.hh>
 
+#include <dune/solvers/common/resize.hh>
+#include <dune/solvers/solvers/umfpacksolver.hh>
+#include "localproblem.hh"
+
 /**
  * \brief linear correction step for use by \ref TNNMGStep
  *
@@ -44,7 +48,25 @@ LinearCorrection<Matrix, Vector>
 makeLinearCorrection(std::shared_ptr< Dune::Solvers::IterativeSolver<Vector> > iterativeSolver)
 {
   return [=](const Matrix& A, Vector& x, const Vector& b, const Dune::Solvers::DefaultBitVector_t<Vector>& ignore) {
-    using LinearIterationStep = Dune::Solvers::LinearIterationStep<Matrix, Vector>;
+
+      // 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)
@@ -53,7 +75,13 @@ makeLinearCorrection(std::shared_ptr< Dune::Solvers::IterativeSolver<Vector> > i
     linearIterationStep->setIgnore(ignore);
     linearIterationStep->setProblem(A, x, b);
     iterativeSolver->preprocess();
-    iterativeSolver->solve();
+    iterativeSolver->solve();*/
+
+    /*
+    const auto& norm = iterativeSolver->getErrorNorm();
+    auto error = norm.diff(linearIterationStep->getSol(), directX);
+
+    std::cout << "CG Solver iterations: " << iterativeSolver->getResult().iterations << " Error: " << error << std::endl;*/
   };
 }
 

src/spatial-solving/tnnmg/linearization.hh

@@ -145,10 +145,8 @@ class Linearization {
       negativeGradient_ *= -1;
 
       // truncate gradient and hessian
-      // truncateVector(negativeGradient_, truncationFlags_);
-      // truncateMatrix(hessian_, truncationFlags_, truncationFlags_);
-
-      // the actual truncation is done by the solver of linearized problem
+      truncateVector(negativeGradient_, truncationFlags_);
+      truncateMatrix(hessian_, truncationFlags_, truncationFlags_);
   }
 
 

src/spatial-solving/tnnmg/tnnmgstep.hh

@@ -125,6 +125,15 @@ class TNNMGStep :
     preSmoothingSteps_ = i;
   }
 
+  auto energy(const Vector& v) const {
+      Vector temp;
+      Dune::Solvers::resizeInitializeZero(temp, v);
+      f_->quadraticPart().umv(v, temp);
+      temp *= 0.5;
+      temp -= f_->linearPart();
+      return temp*v;
+  }
+
   /**
    * \brief Do one TNNMG step
    */
@@ -141,12 +150,14 @@ class TNNMGStep :
 
     //print(f.linearPart(), "f.linearPart():");
 
+    //std::cout << "-- energy before smoothing: " << energy(x) << std::endl;
+
     // Nonlinear presmoothing
     for (std::size_t i=0; i<preSmoothingSteps_; ++i)
         nonlinearSmoother_->iterate();
 
-    std::cout << "- nonlinear presmoothing: success" << std::endl;
-    print(x, "TNNMG iterate after smoothing:");
+    //std::cout << "- nonlinear presmoothing: success" << std::endl;
+    //print(x, "TNNMG iterate after smoothing:");
 
     // Compute constraint/truncated linearization
     if (not linearization_)
@@ -157,48 +168,46 @@ class TNNMGStep :
     auto&& A = linearization_->hessian();
     auto&& r = linearization_->negativeGradient();
 
-    //print(A, "TNNMG Linear problem matrix:");
-    //print(r, "TNNMG Linear problem rhs:");
+    /*print(A, "TNNMG Linear problem matrix:");
+    print(r, "TNNMG Linear problem rhs:");
+    print(ignore, "ignore:");
+    print(linearization_->truncated(), "truncation:");*/
 
     // Compute inexact solution of the linearized problem
     Solvers::resizeInitializeZero(correction_, x);
     Solvers::resizeInitializeZero(constrainedCorrection_, r);
 
-  //  std::cout << "- computing linear correction..." << std::endl;
+    //std::cout << "-- energy after presmoothing: " << energy(x) << std::endl;
 
-    linearCorrection_(A, constrainedCorrection_, r, linearization_->truncated());
+    //print(constrainedCorrection_compare, "direct solver solution:");
 
-    LocalProblem<ConstrainedMatrix, ConstrainedVector> localProblem(A, r, linearization_->truncated());
-    Vector newR, constrainedCorrection_compare;
-    Solvers::resizeInitializeZero(constrainedCorrection_compare, r);
-    localProblem.getLocalRhs(constrainedCorrection_compare, newR);
+    //std::cout << "- computing linear correction..." << std::endl;
 
-    /*print(*this->ignoreNodes_, "ignoreNodes:");
-    print(A, "A:");
-    print(localProblem.getMat(), "localMat:");*/
-    auto directSolver = std::make_shared<Dune::Solvers::UMFPackSolver<ConstrainedMatrix, ConstrainedVector>>();
+    linearCorrection_(A, constrainedCorrection_, r, ignore);
 
-    directSolver->setProblem(localProblem.getMat(), constrainedCorrection_compare, newR);
-    directSolver->preprocess();
-    directSolver->solve();
+    //DUNE_THROW(Dune::Exception, "TNNMGStep: Just need to terminate here!");
 
    // linearCorrection_(A, constrainedCorrection_, r);
 
-    std::cout << "- linear correction: success" << std::endl;
+    //std::cout << "- linear correction: success" << std::endl;
 
     //print(constrainedCorrection_, "contrained correction:");
 
     linearization_->extendCorrection(constrainedCorrection_, correction_);
 
-    std::cout << "- extended correction: success" << std::endl;
-    print(correction_, "correction:");
+    /*Vector h = x;
+    h += correction_;
+    std::cout << "-- energy after linear correction: " << energy(h) << std::endl;*/
+
+    //std::cout << "- extended correction: success" << std::endl;
+    //print(correction_, "correction:");
 
     // Project onto admissible set
     projection_(f, x, correction_);
 
-    std::cout << "- projection onto admissible set: success" << std::endl;
+    //std::cout << "- projection onto admissible set: success" << std::endl;
 
-    print(correction_, "correction:");
+    //print(correction_, "correction:");
 
     double const correctionNorm = correction_.two_norm();
     correction_ /= correctionNorm;
@@ -221,7 +230,11 @@ class TNNMGStep :
       dampingFactor_ = 0;
     correction_ *= dampingFactor_;
 
+    //std::cout << "damping factor: " << dampingFactor_ << std::endl;
+
     x += correction_;
+
+    //std::cout << "-- energy after correction: " << energy(x) << std::endl;
   }
 
   /**

src/time-stepping/adaptivetimestepper.cc

@@ -17,10 +17,10 @@ void IterationRegister::reset() {
   finalCount = FixedPointIterationCounter();
 }
 
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
-AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::AdaptiveTimeStepper(
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
+AdaptiveTimeStepper<Factory, ContactNetwork, Updaters, ErrorNorms>::AdaptiveTimeStepper(
         Dune::ParameterTree const &parset,
-        NBodyAssembler& nBodyAssembler,
+        const ContactNetwork& contactNetwork,
         const IgnoreVector& ignoreNodes,
         GlobalFriction& globalFriction,
         const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
@@ -34,7 +34,7 @@ AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::AdaptiveTime
       relativeTau_(relativeTau),
       finalTime_(parset.get<double>("problem.finalTime")),
       parset_(parset),
-      nBodyAssembler_(nBodyAssembler),
+      contactNetwork_(contactNetwork),
       ignoreNodes_(ignoreNodes),
       globalFriction_(globalFriction),
       bodywiseNonmortarBoundaries_(bodywiseNonmortarBoundaries),
@@ -44,14 +44,13 @@ AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::AdaptiveTime
       mustRefine_(mustRefine),
       errorNorms_(errorNorms) {}
 
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
-bool AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::reachedEnd() {
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
+bool AdaptiveTimeStepper<Factory, ContactNetwork, Updaters, ErrorNorms>::reachedEnd() {
   return relativeTime_ >= 1.0;
 }
 
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
-template <class Preconditioner>
-IterationRegister AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::advance(Preconditioner&& preconditioner) {
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
+IterationRegister AdaptiveTimeStepper<Factory, ContactNetwork, Updaters, ErrorNorms>::advance() {
   /*
     |     C     | We check here if making the step R1 of size tau is a
     |  R1 | R2  | good idea. To check if we can coarsen, we compare
@@ -64,7 +63,7 @@ IterationRegister AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNo
   std::cout << "AdaptiveTimeStepper::advance()" << std::endl;
 
   if (R1_.updaters == Updaters())
-    R1_ = step(current_, relativeTime_, relativeTau_, std::forward<Preconditioner>(preconditioner));
+    R1_ = step(current_, relativeTime_, relativeTau_);
 
   std::cout << "AdaptiveTimeStepper Step 1" << std::endl;
 
@@ -74,9 +73,9 @@ IterationRegister AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNo
   UpdatersWithCount R2;
   UpdatersWithCount C;
   while (relativeTime_ + relativeTau_ <= 1.0) {
-    R2 = step(R1_.updaters, relativeTime_ + relativeTau_, relativeTau_, std::forward<Preconditioner>(preconditioner));
+    R2 = step(R1_.updaters, relativeTime_ + relativeTau_, relativeTau_);
     std::cout << "AdaptiveTimeStepper R2 computed!" << std::endl << std::endl;
-    C = step(current_, relativeTime_, 2 * relativeTau_, std::forward<Preconditioner>(preconditioner));
+    C = step(current_, relativeTime_, 2 * relativeTau_);
     std::cout << "AdaptiveTimeStepper C computed!" << std::endl << std::endl;
     if (mustRefine_(R2.updaters, C.updaters))
       break;
@@ -91,10 +90,10 @@ IterationRegister AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNo
   UpdatersWithCount F2;
   if (!didCoarsen) {
     while (true) {
-      F1 = step(current_, relativeTime_, relativeTau_ / 2.0, std::forward<Preconditioner>(preconditioner));
+      F1 = step(current_, relativeTime_, relativeTau_ / 2.0);
       std::cout << "AdaptiveTimeStepper F1 computed!" << std::endl << std::endl;
       F2 = step(F1.updaters, relativeTime_ + relativeTau_ / 2.0,
-                relativeTau_ / 2.0, std::forward<Preconditioner>(preconditioner));
+                relativeTau_ / 2.0);
       std::cout << "AdaptiveTimeStepper F2 computed!" << std::endl << std::endl;
       if (!mustRefine_(F2.updaters, R1_.updaters)) {
         std::cout << "Sufficiently refined!" << std::endl;
@@ -122,17 +121,16 @@ IterationRegister AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNo
   return iterationRegister_;
 }
 
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
-template <class Preconditioner>
-typename AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::UpdatersWithCount
-AdaptiveTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::step(
-    Updaters const &oldUpdaters, double rTime, double rTau, Preconditioner&& preconditioner) {
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
+typename AdaptiveTimeStepper<Factory, ContactNetwork, Updaters, ErrorNorms>::UpdatersWithCount
+AdaptiveTimeStepper<Factory, ContactNetwork, Updaters, ErrorNorms>::step(
+    Updaters const &oldUpdaters, double rTime, double rTau) {
   UpdatersWithCount newUpdatersAndCount = {oldUpdaters.clone(), {}};
   newUpdatersAndCount.count =
-      MyCoupledTimeStepper(finalTime_, parset_, nBodyAssembler_, ignoreNodes_, globalFriction_, bodywiseNonmortarBoundaries_,
+      MyCoupledTimeStepper(finalTime_, parset_, contactNetwork_, ignoreNodes_, globalFriction_, bodywiseNonmortarBoundaries_,
                            newUpdatersAndCount.updaters, errorNorms_,
                            externalForces_)
-          .step(rTime, rTau, std::forward<Preconditioner>(preconditioner));
+          .step(rTime, rTau);
   iterationRegister_.registerCount(newUpdatersAndCount.count);
   return newUpdatersAndCount;
 }

src/time-stepping/adaptivetimestepper.hh

@@ -15,7 +15,7 @@ struct IterationRegister {
   FixedPointIterationCounter finalCount;
 };
 
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
 class AdaptiveTimeStepper {
   struct UpdatersWithCount {
     Updaters updaters;
@@ -27,7 +27,7 @@ class AdaptiveTimeStepper {
   //using ConvexProblem = typename Factory::ConvexProblem;
   //using Nonlinearity = typename Factory::Nonlinearity;
 
-  using MyCoupledTimeStepper = CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>;
+  using MyCoupledTimeStepper = CoupledTimeStepper<Factory, ContactNetwork, Updaters, ErrorNorms>;
 
   using GlobalFriction = typename MyCoupledTimeStepper::GlobalFriction;
   using BitVector = typename MyCoupledTimeStepper::BitVector;
@@ -36,7 +36,7 @@ class AdaptiveTimeStepper {
 public:
   AdaptiveTimeStepper(
                       Dune::ParameterTree const &parset,
-                      NBodyAssembler& nBodyAssembler,
+                      const ContactNetwork& contactNetwork,
                       const IgnoreVector& ignoreNodes,
                       GlobalFriction& globalFriction,
                       const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
@@ -49,20 +49,18 @@ class AdaptiveTimeStepper {
 
   bool reachedEnd();
 
-  template <class Preconditioner>
-  IterationRegister advance(Preconditioner&& preconditioner);
+  IterationRegister advance();
 
   double relativeTime_;
   double relativeTau_;
 
 private:
-  template <class Preconditioner>
   UpdatersWithCount step(Updaters const &oldUpdaters, double rTime,
-                         double rTau, Preconditioner&& preconditioner);
+                         double rTau);
 
   double finalTime_;
   Dune::ParameterTree const &parset_;
-  NBodyAssembler& nBodyAssembler_;
+  const ContactNetwork& contactNetwork_;
   const IgnoreVector& ignoreNodes_;
 
   GlobalFriction& globalFriction_;

src/time-stepping/adaptivetimestepper_tmpl.cc

@@ -14,7 +14,6 @@
 #include "state/stateupdater.hh"
 #include "updaters.hh"
 
-using NBodyAssembler = typename MyContactNetwork::NBodyAssembler;
 using BoundaryNodes = typename MyContactNetwork::BoundaryNodes;
 using BoundaryFunctions = typename MyContactNetwork::BoundaryFunctions;
 
@@ -24,6 +23,6 @@ using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
 
 using ErrorNorms = typename MyContactNetwork::StateEnergyNorms;
 
-using MyAdaptiveTimeStepper = AdaptiveTimeStepper<MySolverFactory, NBodyAssembler, MyUpdaters, ErrorNorms>;
+using MyAdaptiveTimeStepper = AdaptiveTimeStepper<MySolverFactory, MyContactNetwork, MyUpdaters, ErrorNorms>;
 
-template class AdaptiveTimeStepper<MySolverFactory, NBodyAssembler, MyUpdaters, ErrorNorms>;
+template class AdaptiveTimeStepper<MySolverFactory, MyContactNetwork, MyUpdaters, ErrorNorms>;

src/time-stepping/coupledtimestepper.cc

@@ -4,10 +4,10 @@
 
 #include "coupledtimestepper.hh"
 
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
-CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::CoupledTimeStepper(
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
+CoupledTimeStepper<Factory, ContactNetwork, Updaters, ErrorNorms>::CoupledTimeStepper(
     double finalTime, Dune::ParameterTree const &parset,
-    NBodyAssembler& nBodyAssembler,
+    const ContactNetwork& contactNetwork,
     const IgnoreVector& ignoreNodes,
     GlobalFriction& globalFriction,
     const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
@@ -16,7 +16,7 @@ CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::CoupledTimeSt
     ExternalForces& externalForces)
     : finalTime_(finalTime),
       parset_(parset),
-      nBodyAssembler_(nBodyAssembler),
+      contactNetwork_(contactNetwork),
       ignoreNodes_(ignoreNodes),
       globalFriction_(globalFriction),
       bodywiseNonmortarBoundaries_(bodywiseNonmortarBoundaries),
@@ -24,11 +24,10 @@ CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::CoupledTimeSt
       externalForces_(externalForces),
       errorNorms_(errorNorms) {}
 
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
-template <class Preconditioner>
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
 FixedPointIterationCounter
-CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::step(double relativeTime,
-                                                       double relativeTau, Preconditioner&& preconditioner) {
+CoupledTimeStepper<Factory, ContactNetwork, Updaters, ErrorNorms>::step(double relativeTime,
+                                                       double relativeTau) {
 
   std::cout << "CoupledTimeStepper::step()" << std::endl;
 
@@ -50,8 +49,8 @@ CoupledTimeStepper<Factory, NBodyAssembler, Updaters, ErrorNorms>::step(double r
   updaters_.rate_->setup(ell, tau, newRelativeTime, velocityRHS, velocityIterate, velocityMatrix);
 
   FixedPointIterator fixedPointIterator(
-      parset_, nBodyAssembler_, ignoreNodes_, globalFriction_, bodywiseNonmortarBoundaries_, errorNorms_);
-  auto const iterations = fixedPointIterator.run(updaters_, std::forward<Preconditioner>(preconditioner),
+      parset_, contactNetwork_, ignoreNodes_, globalFriction_, bodywiseNonmortarBoundaries_, errorNorms_);
+  auto const iterations = fixedPointIterator.run(updaters_,
                                                  velocityMatrix, velocityRHS, velocityIterate);
   return iterations;
 }

src/time-stepping/coupledtimestepper.hh

@@ -8,12 +8,12 @@
 
 #include "../spatial-solving/fixedpointiterator.hh"
 
-template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
 class CoupledTimeStepper {
   using Vector = typename Factory::Vector;
   using Matrix = typename Factory::Matrix;
   using IgnoreVector = typename Factory::BitVector;
-  using FixedPointIterator = FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>;
+  using FixedPointIterator = FixedPointIterator<Factory, ContactNetwork, Updaters, ErrorNorms>;
 
 public:
   using GlobalFriction = typename FixedPointIterator::GlobalFriction;
@@ -23,7 +23,7 @@ class CoupledTimeStepper {
 public:
   CoupledTimeStepper(double finalTime,
                      Dune::ParameterTree const &parset,
-                     NBodyAssembler& nBodyAssembler,
+                     const ContactNetwork& contactNetwork,
                      const IgnoreVector& ignoreNodes,
                      GlobalFriction& globalFriction,
                      const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
@@ -31,13 +31,12 @@ class CoupledTimeStepper {
                      const ErrorNorms& errorNorms,
                      ExternalForces& externalForces);
 
-  template <class Preconditioner>
-  FixedPointIterationCounter step(double relativeTime, double relativeTau, Preconditioner&& preconditioner);
+  FixedPointIterationCounter step(double relativeTime, double relativeTau);
 
 private:
   double finalTime_;
   Dune::ParameterTree const &parset_;
-  NBodyAssembler& nBodyAssembler_;
+  const ContactNetwork& contactNetwork_;
   const IgnoreVector& ignoreNodes_;
 
   GlobalFriction& globalFriction_;