dune/tectonic/time-stepping/uniformtimestepper.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/iterationsteps/multigridstep.hh>
+#include <dune/solvers/iterationsteps/cgstep.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+
+#include "../spatial-solving/preconditioners/multilevelpatchpreconditioner.hh"
+
+#include <dune/tectonic/utils/reductionfactors.hh>
+
+#include "uniformtimestepper.hh"
+#include "step.hh"
+
+
+/*
+ * Implementation: AdaptiveTimeStepper
+ */
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
+UniformTimeStepper<Factory, ContactNetwork, Updaters, ErrorNorms>::UniformTimeStepper(
+        const StepBase& stepBase,
+        ContactNetwork& contactNetwork,
+        Updaters &current,
+        double relativeTime,
+        double relativeTau)
+    : relativeTime_(relativeTime),
+      relativeTau_(relativeTau),
+      stepBase_(stepBase),
+      contactNetwork_(contactNetwork),
+      current_(current) {}
+
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
+bool UniformTimeStepper<Factory, ContactNetwork, Updaters, ErrorNorms>::reachedEnd() {
+  return relativeTime_ >= 1.0;
+}
+
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
+IterationRegister UniformTimeStepper<Factory, ContactNetwork, Updaters, ErrorNorms>::advance() {
+  //std::cout << "AdaptiveTimeStepper::advance()" << std::endl;
+
+  using Step = Step<Factory, ContactNetwork, Updaters, ErrorNorms>;
+  iterationRegister_.reset();
+
+  UpdatersWithCount N;
+
+  auto step = Step(stepBase_, current_, contactNetwork_.nBodyAssembler(), relativeTime_, relativeTau_, iterationRegister_);
+  step.run(Step::Mode::sameThread);
+  N = step.get();
+
+  current_ = N.updaters;
+  iterationRegister_.registerFinalCount(N.count);
+  relativeTime_ += relativeTau_;
+
+  return iterationRegister_;
+}
+
+#include "uniformtimestepper_tmpl.cc"

dune/tectonic/time-stepping/uniformtimestepper.hh

+#ifndef SRC_TIME_STEPPING_UNIFORMTIMESTEPPER_HH
+#define SRC_TIME_STEPPING_UNIFORMTIMESTEPPER_HH
+
+#include "../spatial-solving/fixedpointiterator.hh"
+#include "adaptivetimestepper.hh"
+#include "stepbase.hh"
+
+template <class Factory, class ContactNetwork, class Updaters, class ErrorNorms>
+class UniformTimeStepper {
+
+  using UpdatersWithCount = UpdatersWithCount<Updaters>;
+  using StepBase = StepBase<Factory, ContactNetwork, Updaters, ErrorNorms>;
+
+public:
+  UniformTimeStepper(const StepBase& stepBase,
+                      ContactNetwork& contactNetwork,
+                      Updaters &current,
+                      double relativeTime,
+                      double relativeTau);
+
+  bool reachedEnd();
+
+  IterationRegister advance();
+
+  double relativeTime_;
+  const double relativeTau_;
+
+private:
+  const StepBase& stepBase_;
+  ContactNetwork& contactNetwork_;
+
+  Updaters &current_;
+
+  IterationRegister iterationRegister_;
+};
+
+#endif

dune/tectonic/time-stepping/uniformtimestepper_tmpl.cc

+#ifndef MY_DIM
+#error MY_DIM unset
+#endif
+
+#include <future>
+#include <thread>
+
+#include "../explicitgrid.hh"
+#include "../explicitvectors.hh"
+
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+
+#include "../spatial-solving/tnnmg/functional.hh"
+#include "../spatial-solving/solverfactory.hh"
+
+#include "../data-structures/network/contactnetwork.hh"
+#include "../data-structures/friction/globalfriction.hh"
+#include "../spatial-solving/tnnmg/zerononlinearity.hh"
+
+#include "rate/rateupdater.hh"
+#include "state/stateupdater.hh"
+#include "updaters.hh"
+
+using MyContactNetwork = ContactNetwork<Grid, Vector>;
+
+using BoundaryNodes = typename MyContactNetwork::BoundaryNodes;
+using BoundaryFunctions = typename MyContactNetwork::BoundaryFunctions;
+
+using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
+using MyRateUpdater = RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>;
+using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
+
+using LinearSolver = Dune::Solvers::LoopSolver<Vector>;
+using MyNBodyAssembler = typename MyContactNetwork::NBodyAssembler;
+using ErrorNorms = typename MyContactNetwork::StateEnergyNorms;
+
+using MyGlobalFriction = GlobalFriction<Matrix, Vector>;
+using MyFunctional = Functional<Matrix&, Vector&, MyGlobalFriction&, Vector&, Vector&, double>;
+using MySolverFactory = SolverFactory<MyFunctional, BitVector>;
+
+template class UniformTimeStepper<MySolverFactory, MyContactNetwork, MyUpdaters, ErrorNorms>;
+
+
+using NoFriction = ZeroNonlinearity;
+using NoFrictionFunctional = Functional<Matrix&, Vector&, NoFriction&, Vector&, Vector&, double>;
+using NoFrictionSolverFactory = SolverFactory<NoFrictionFunctional, BitVector>;
+
+template class UniformTimeStepper<NoFrictionSolverFactory, MyContactNetwork, MyUpdaters, ErrorNorms>;
+
+/*
+template std::packaged_task<typename AdaptiveTimeStepper<MySolverFactory, MyContactNetwork, MyUpdaters, ErrorNorms>::UpdatersWithCount()>
+AdaptiveTimeStepper<MySolverFactory, MyContactNetwork, MyUpdaters, ErrorNorms>::step<LinearSolver>(
+        const MyUpdaters&, const MyNBodyAssembler&, std::shared_ptr<LinearSolver>&, double, double); */

dune/tectonic/utils/CMakeLists.txt

+add_custom_target(tectonic_dune_utils SOURCES
+  almostequal.hh
+  debugutils.hh 
+  diameter.hh
+  geocoordinate.hh
+  index-in-sorted-range.hh
+  maxeigenvalue.hh
+  tobool.hh
+)
+
+#install headers
+install(FILES
+  almostequal.hh
+  debugutils.hh 
+  diameter.hh
+  geocoordinate.hh
+  index-in-sorted-range.hh
+  maxeigenvalue.hh
+  tobool.hh
+DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)

dune/tectonic/utils/almostequal.hh

+#ifndef SRC_ALMOSTEQUAL_HH
+#define SRC_ALMOSTEQUAL_HH
+
+#include <type_traits>
+#include <limits>
+#include <math.h>
+
+template <typename ctype>
+typename std::enable_if<!std::numeric_limits<ctype>::is_integer, bool>::type almost_equal(ctype x, ctype y, int ulp) {
+    return std::abs(x-y) < std::numeric_limits<ctype>::epsilon() * std::abs(x+y) * ulp || std::abs(x-y) < std::numeric_limits<ctype>::min();
+}
+
+#endif

dune/tectonic/utils/debugutils.hh

+#ifndef DEBUG_UTILS_
+#define DEBUG_UTILS_
+
+#include <string>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/bitsetvector.hh>
+
+#include <dune/istl/matrix.hh>
+#include <dune/istl/bcrsmatrix.hh>
+
+#include <dune/grid/io/file/vtk/vtkwriter.hh>
+
+//using namespace std;
+
+    template <class Friction, typename VectorType>
+    void printRegularityTruncation(const Friction& friction, const VectorType& x, double truncationTolerance = 1e8) {
+        using BitVector = Dune::BitSetVector<MY_DIM>;
+
+        BitVector truncationFlags(x.size());
+        truncationFlags.unsetAll();
+
+        size_t count = 0;
+        size_t vsmall = 0;
+
+
+        for (size_t i = 0; i < x.size(); ++i) {
+          //std::cout << f_.phi().regularity(i, x[i]) << " xnorm: " << x[i].two_norm() << std::endl;
+          auto tangential_x = x[i];
+          tangential_x[0] = 0.0;
+
+          if (tangential_x.two_norm()<1e-14) {
+              vsmall++;
+          }
+
+          if (friction.regularity(i, x[i]) > truncationTolerance) {
+            count++;
+          }
+        }
+
+         std::cout << "V<1e-14: " << vsmall << " regularityTruncation: " << count << std::endl;
+    }
+
+    template <typename VectorType>
+    auto average(const VectorType& vec) {
+        using BlockType = typename VectorType::block_type;
+
+        BlockType res(0.0);
+        auto len = vec.size();
+
+        if (len>0) {
+            for (size_t i=0; i<vec.size(); i++) {
+                res += vec[i];
+            }
+            res *= 1.0/len;
+        }
+
+        return res;
+    }
+
+    template <typename VectorType>
+    auto minmax(const VectorType& vec) {
+        using BlockType = typename VectorType::block_type;
+
+        BlockType min(0.0), max(0.0);
+        auto len = vec.size();
+
+        if (len>0) {
+            for (size_t i=0; i<vec.size(); i++) {
+                if (vec[i] > max) {
+                    max = vec[i];
+                }
+                if (vec[i] < min) {
+                    min = vec[i];
+                }
+            }
+        }
+        return std::array{min, max};
+    }
+
+    template <int s>
+    void print(const Dune::BitSetVector<s>& x, std::string message){
+        std::cout << message << std::endl;
+        for (size_t i=0; i<x.size(); i++) {
+            std::cout << "block " << i << ": ";
+            for (size_t j=0; j<s; j++) {
+                 std::cout << x[i][j] << " ";
+            }
+            std::cout << std::endl;
+        }
+        std::cout << std::endl << std::endl;
+    }
+
+    template <int dim, typename ctype=double>
+    void print(const Dune::FieldVector<ctype, dim>& x, std::string message){
+       if (message != "")
+        std::cout << message << std::endl;
+
+       for (size_t i=0; i<x.size(); i++) {
+           std::cout << x[i] << ";" << std::endl;
+       }
+
+   }
+
+    template <int dim, typename ctype=double>
+    void print(const Dune::BlockVector<Dune::FieldVector<ctype, dim>>& x, std::string message){
+       std::cout << message  << " " << "size " << x.size() << std::endl;
+       for (size_t i=0; i<x.size(); i++) {
+           print(x[i], "");
+       }
+       std::cout << std::endl << std::endl;
+   }
+
+   template <int dim, typename ctype=double>
+   void print(const Dune::FieldMatrix<ctype, dim, dim>& mat, std::string message, bool endOfLine = true){
+        if (message != "")
+            std::cout << message << std::endl;
+
+        for (size_t i=0; i<mat.N(); i++) {
+            for (size_t j=0; j<mat.M(); j++) {
+                if (mat.exists(i,j))
+                    std::cout << mat[i][j] << " ";
+                else
+                    std::cout << "0 ";
+            }
+
+            if (endOfLine)
+                std::cout << ";"<< std::endl;
+        }
+   }
+
+    template <int dim, typename ctype=double>
+    void print(const Dune::FieldMatrix<ctype, dim, dim>& mat, std::string message, const size_t line, bool endOfLine = true){
+        if (message != "")
+            std::cout << message << std::endl;
+
+        assert(line<mat.N());
+        for (size_t j=0; j<mat.M(); j++) {
+            if (mat.exists(line,j))
+                    std::cout << mat[line][j] << " ";
+                else
+                    std::cout << "0 ";
+            }
+
+            if (endOfLine)
+                std::cout << ";"<< std::endl;
+   }
+
+    template <int dim, typename ctype=double>
+    void print(const Dune::BCRSMatrix<Dune::FieldMatrix<ctype, dim, dim>>& mat, std::string message){
+        std::cout << message << " " << "size (" << mat.N() << ", " << mat.M() << ")" <<std::endl;
+
+        Dune::FieldMatrix<ctype, dim, dim> zeroEntry = 0;
+
+        for (size_t i=0; i<mat.N(); i++) {
+            for (size_t d=0; d<dim; d++) {
+                for (size_t j=0; j<mat.M(); j++) {
+                    if (mat.exists(i,j))
+                        print(mat[i][j], "", d, j==mat.M()-1);
+                    else
+                        print(zeroEntry, "", d, j==mat.M()-1);
+                }
+            }
+       }
+       std::cout << std::endl;
+   }
+
+   template <class T=Dune::FieldMatrix<double,1,1>>
+   void print(const Dune::Matrix<T>& mat, std::string message){
+       std::cout << message << std::endl;
+       for (size_t i=0; i<mat.N(); i++) {
+           for (size_t j=0; j<mat.M(); j++) {
+                   std::cout << mat[i][j] << " ";
+           }
+           std::cout << std::endl;
+       }
+       std::cout << std::endl;
+   }
+
+   template <class T>
+   void print(const std::vector<T>& x, std::string message){
+      std::cout << message << std::endl;
+      for (size_t i=0; i<x.size(); i++) {
+          std::cout << x[i] << std::endl;
+          //print(x[i], "");
+      }
+      std::cout << std::endl << std::endl;
+   }
+
+ /*  void print(const std::vector<Dune::FieldVector<double,1>>& x, std::string message){
+      std::cout << message << std::endl;
+      for (size_t i=0; i<x.size(); i++) {
+          std::cout << x[i] << " ";
+      }
+      std::cout << std::endl << std::endl;
+   }*/
+
+   template <class T>
+   void print(const std::set<T>& x, std::string message){
+      std::cout << message << std::endl;
+
+      for (const auto& c : x) {
+          std::cout << c << " ";
+      }
+      std::cout << std::endl << std::endl;
+   }
+/*
+   void print(const std::set<int>& x, std::string message){
+      std::cout << message << std::endl;
+      std::set<int>::iterator dofIt = x.begin();
+      std::set<int>::iterator dofEndIt = x.end();
+      for (; dofIt != dofEndIt; ++dofIt) {
+          std::cout << *dofIt << " ";
+      }
+      std::cout << std::endl << std::endl;
+   }
+
+   void step(const double stepNumber, std::string message=""){
+      std::cout << message << " Step " << stepNumber << "!" << std::endl;
+   }*/
+
+   template <class GridView, class VectorType>
+   void writeToVTK(const GridView& gridView, const VectorType& x, const std::string path, const std::string name) {
+       Dune::VTKWriter<GridView> vtkwriter(gridView);
+
+       std::ofstream lStream( "garbage.txt" );
+       std::streambuf* lBufferOld = std::cout.rdbuf();
+       std::cout.rdbuf(lStream.rdbuf());
+
+       vtkwriter.addVertexData(x,"data");
+       vtkwriter.pwrite(name, path, "");
+
+       std::cout.rdbuf( lBufferOld );
+   }
+
+   template <class GridView>
+   void writeToVTK(const GridView& gridView, const std::string path, const std::string name) {
+       Dune::VTKWriter<GridView> vtkwriter(gridView);
+
+       std::ofstream lStream( "garbage.txt" );
+       std::streambuf* lBufferOld = std::cout.rdbuf();
+       std::cout.rdbuf(lStream.rdbuf());
+
+       vtkwriter.pwrite(name, path, "");
+
+       std::cout.rdbuf( lBufferOld );
+   }
+
+   /*
+   template <class VectorType, class DGBasisType>
+   void writeToVTK(const DGBasisType& basis, const VectorType& x, const std::string path, const std::string name) {
+       Dune::VTKWriter<typename DGBasisType::GridView> vtkwriter(basis.getGridView());
+       VectorType toBePlotted(x);
+
+       toBePlotted.resize(basis.getGridView().indexSet().size(DGBasisType::GridView::Grid::dimension));
+
+       std::ofstream lStream( "garbage.txt" );
+       std::streambuf* lBufferOld = std::cout.rdbuf();
+       std::cout.rdbuf( lStream.rdbuf() );
+
+       vtkwriter.addVertexData(toBePlotted,"data");
+       vtkwriter.pwrite(name, path, path);
+
+       std::cout.rdbuf( lBufferOld );
+   }*/
+
+   template <class Intersection, class Basis>
+   void printIntersection(const Intersection& it, const Basis& basis, std::string message = "") {
+       if (message != "") {
+           std::cout << message << std::endl;
+       }
+
+       const auto& inside = it.inside();
+       const auto& localCoefficients = basis.getLocalFiniteElement(inside).localCoefficients();
+
+       std::cout << "dofs: ";
+       for (size_t i=0; i<localCoefficients.size(); i++) {
+           unsigned int entity = localCoefficients.localKey(i).subEntity();
+           unsigned int codim  = localCoefficients.localKey(i).codim();
+
+           if (it.containsInsideSubentity(entity, codim)) {
+               int dofIndex = basis.index(it.inside(), i);
+               std::cout << dofIndex << " ";
+               break;
+           }
+       }
+       std::cout << std::endl;
+   }
+
+   template <class BoundaryPatch, class Basis>
+   void printBoundary(const BoundaryPatch& patch, const Basis& basis, std::string message = "") {
+       if (message != "") {
+           std::cout << message << std::endl;
+       }
+
+       for (auto bIt = patch.begin(); bIt != patch.end(); ++bIt) {
+            printIntersection(*bIt, basis, "");
+       }
+
+       std::cout << std::endl;
+   }
+
+   template <class BasisType, typename ctype=double>
+   void printBasisDofLocation(const BasisType& basis) {
+     /*  typedef typename BasisType::GridView GridView;
+
+       const int dim = GridView::dimension;
+
+       std::map<int, int> indexTransformation;
+       std::map<int, Dune::FieldVector<ctype, dim>> indexCoords;
+
+
+       const GridView& gridView = basis.getGridView();
+       const int gridN = std::pow(gridView.indexSet().size(dim), 1.0/dim)-1;
+
+       for (auto& it : elements(gridView)) {
+           const typename BasisType::LocalFiniteElement& tFE = basis.getLocalFiniteElement(it);
+           const auto& geometry = it.geometry();
+
+           for(int i=0; i<geometry.corners(); ++i) {
+               const auto& vertex = geometry.corner(i);
+               const auto& local = geometry.local(vertex);
+
+               int vertexIndex = vertex[0]*gridN;
+               for (int j=1; j<dim; ++j){
+                   vertexIndex += vertex[j]*gridN*std::pow(gridN+1, j);
+               }
+
+               const int localBasisSize = tFE.localBasis().size();
+               std::vector<Dune::FieldVector<double, 1>, std::allocator<Dune::FieldVector<double, 1> > > localBasisRep(localBasisSize);
+               tFE.localBasis().evaluateFunction(local, localBasisRep);
+
+               for(int k=0; k<localBasisSize; ++k) {
+                   if (localBasisRep[k]==1) {
+                       int dofIndex = basis.index(it, k);
+
+                       if (indexTransformation.count(dofIndex)==0) {
+                           indexTransformation[dofIndex] = vertexIndex;
+                           indexCoords[dofIndex] = vertex;
+                       }
+
+                       break;
+                   }
+               }
+           }
+       }
+
+       std::cout << std::endl;
+       std::cout << "Coarse level: Geometric vertex indices vs. associated basis dof indices " << indexTransformation.size() << std::endl;
+       std::map<int,int>::iterator mapIt = indexTransformation.begin();
+       std::map<int,int>::iterator mapEndIt = indexTransformation.end();
+       for (; mapIt!=mapEndIt; ++mapIt) {
+           std::cout << mapIt->second << " => " << mapIt->first << " Coords: ";
+
+           const Dune::FieldVector<ctype, dim>& coords = indexCoords[mapIt->first];
+           for (size_t i=0; i<coords.size(); i++) {
+               std::cout << coords[i] << " ";
+           }
+           std::cout << std::endl;
+       }
+       std::cout << std::endl;*/
+
+       const int dim = BasisType::GridView::dimension;
+
+       std::map<int, Dune::FieldVector<ctype, dim>> indexCoords;
+
+
+       const auto& gridView = basis.getGridView();
+
+       for (auto& it : elements(gridView)) {
+            const typename BasisType::LocalFiniteElement& tFE = basis.getLocalFiniteElement(it);
+            const auto& geometry = it.geometry();
+
+            for(int i=0; i<geometry.corners(); ++i) {
+                const auto& vertex = geometry.corner(i);
+                const auto& local = geometry.local(vertex);
+
+                const int localBasisSize = tFE.localBasis().size();
+                std::vector<Dune::FieldVector<double, 1>, std::allocator<Dune::FieldVector<double, 1> > > localBasisRep(localBasisSize);
+                tFE.localBasis().evaluateFunction(local, localBasisRep);
+
+                for(int k=0; k<localBasisSize; ++k) {
+                    if (1.0 - localBasisRep[k]< 1e-14) {
+                        int dofIndex = basis.index(it, k);
+                        indexCoords[dofIndex] = vertex;
+                        break;
+                    }
+                }
+              }
+       }
+
+       std::cout << "Coords of basis dofs: " << indexCoords.size() << std::endl;
+       auto&& mapIt = indexCoords.begin();
+       const auto& mapEndIt = indexCoords.end();
+       for (; mapIt!=mapEndIt; ++mapIt) {
+            std::cout << mapIt->first << " Coords: ";
+
+            const auto& coords = mapIt->second;
+            for (size_t i=0; i<coords.size(); i++) {
+                std::cout << coords[i] << " ";
+            }
+            std::cout << std::endl;
+       }
+
+       std::cout << std::endl;
+   }
+
+   template <class GridView>
+   void printDofLocation(const GridView& gridView) {
+
+       std::cout << "-- GridView vertices --" << std::endl;
+       std::cout << "Index   Coords         " << std::endl;
+       std::cout << "-----------------------" << std::endl;
+
+       Dune::MultipleCodimMultipleGeomTypeMapper<
+           GridView, Dune::MCMGVertexLayout> const vertexMapper(gridView, Dune::mcmgVertexLayout());
+
+       for (auto&& it : vertices(gridView)) {
+            const auto i = vertexMapper.index(it);
+            const auto& geometry = it.geometry();
+            const auto& corner = geometry.corner(0);
+
+            std::cout << std::setw(5) << i << "   ";
+            for(size_t i=0; i<corner.size(); ++i) {
+                std::cout << std::setw(5) << std::setprecision(3) << corner[i] << " ";
+            }
+            std::cout << std::endl;
+       }
+       std::cout << std::endl;
+   }
+
+   template <class Vector, class NBodyAssembler>
+   void printMortarBasis(const NBodyAssembler& nBodyAssembler) {
+       std::cout << "-- Mortar basis in canonical coords --" << std::endl;
+       std::cout << "--------------------------------------" << std::endl;
+
+       const auto& dim = nBodyAssembler.getTransformationMatrix().N();
+
+       Vector mortarBasisVector(dim);
+       std::vector<Vector> canonicalBasisVector(nBodyAssembler.nGrids());
+       for (size_t i=0; i<dim; i++) {
+            mortarBasisVector = 0;
+            mortarBasisVector[i] = 1;
+
+            nBodyAssembler.postprocess(mortarBasisVector, canonicalBasisVector);
+            print(canonicalBasisVector, "canonicalBasisVector " + std::to_string(i));
+       }
+       std::cout << std::endl;
+   }
+#endif

dune/tectonic/utils/maxeigenvalue.hh

+// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
+// vi: set et ts=4 sw=2 sts=2:
+#ifndef DUNE_TECTONIC_UTILS_MAXEIGENVALUES_HH
+#define DUNE_TECTONIC_UTILS_MAXEIGENVALUES_HH
+
+/** \file
+ * \brief max eigenvalue computations for the FieldMatrix class
+ */
+
+#include <iostream>
+#include <cmath>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/fmatrixev.hh>
+
+template <typename K>
+static K maxEigenvalue(const Dune::FieldMatrix<K, 1, 1>& matrix) {
+    return matrix[0][0];
+}
+
+/** \brief calculates the eigenvalues of a symmetric field matrix
+    \param[in]  matrix matrix eigenvalues are calculated for
+    \param[out] max eigenvalue
+*/
+template <typename K>
+static K maxEigenvalue(const Dune::FieldMatrix<K, 2, 2>& matrix) {
+    const K p = 0.5 * (matrix[0][0] + matrix [1][1]);
+    const K p2 = p - matrix[1][1];
+    K q = p2 * p2 + matrix[1][0] * matrix[0][1];
+    if( q < 0 && q > -1e-14 ) q = 0;
+    if (q < 0) {
+        std::cout << matrix << std::endl;
+        // Complex eigenvalues are either caused by non-symmetric matrices or by round-off errors
+        DUNE_THROW(Dune::MathError, "Complex eigenvalue detected (which this implementation cannot handle).");
+    }
+
+    // get square root
+    q = std::sqrt(q);
+
+    return p + q;
+}
+
+/** \brief Calculates the eigenvalues of a symmetric 3x3 field matrix
+    \param[in]  matrix eigenvalues are calculated for
+    \param[out] max eigenvalue
+
+    \note If the input matrix is not symmetric the behavior of this method is undefined.
+ */
+template <typename K>
+static K maxEigenvalue(const Dune::FieldMatrix<K, 3, 3>& matrix) {
+    Dune::FieldVector<K, 3> eigenvalues;
+    Dune::FMatrixHelp::eigenValues(matrix, eigenvalues);
+    return eigenvalues[0];
+}
+
+  /** @} end documentation */
+
+#endif

dune/tectonic/utils/reductionfactors.hh

+extern std::vector<std::vector<double>> allReductionFactors;

src/tobool.hh → dune/tectonic/utils/tobool.hh

@@ -4,7 +4,7 @@
 #include <dune/common/bitsetvector.hh>
 
 template <class Alloc>
-bool toBool(Dune::BitSetVectorConstReference<1, Alloc> x) {
+bool toBool(const Dune::BitSetVectorConstReference<1, Alloc> x) {
   return x[0];
 }
 #endif

src/CMakeLists.txt

-set(SW_SOURCE_FILES
-  assemblers.cc
-  enumparser.cc
-  hdf5/frictionalboundary-writer.cc
-  hdf5/iteration-writer.cc
-  hdf5/patchinfo-writer.cc
-  hdf5/restart-io.cc
-  hdf5/surface-writer.cc
-  hdf5/time-writer.cc
-  one-body-problem-data/mygeometry.cc
-  one-body-problem-data/mygrid.cc
-  one-body-problem.cc
-  spatial-solving/fixedpointiterator.cc
-  spatial-solving/solverfactory.cc
-  time-stepping/adaptivetimestepper.cc
-  time-stepping/coupledtimestepper.cc
-  time-stepping/rate.cc
-  time-stepping/rate/rateupdater.cc
-  time-stepping/state.cc
-  vtk.cc
-)
+add_subdirectory("foam")
+add_subdirectory("multi-body-problem")
+add_subdirectory("strikeslip")
+
 set(UGW_SOURCE_FILES
-  assemblers.cc # FIXME
-  one-body-problem-data/mygrid.cc
-  uniform-grid-writer.cc
-  vtk.cc
+  ../dune/tectonic/assemblers.cc # FIXME
+  #../dune/tectonic/io/uniform-grid-writer.cc
+  ../dune/tectonic/io/vtk.cc
+  ../dune/tectonic/problem-data/grid/mygrids.cc
 )
 
 foreach(_dim 2 3)
-  set(_sw_target one-body-problem-${_dim}D)
-  set(_ugw_target uniform-grid-writer-${_dim}D)
-
-  add_executable(${_sw_target} ${SW_SOURCE_FILES})
-  add_executable(${_ugw_target} ${UGW_SOURCE_FILES})
-  add_dune_ug_flags(${_sw_target})
-  add_dune_ug_flags(${_ugw_target})
+  #set(_ugw_target uniform-grid-writer-${_dim}D)
 
-  add_dune_hdf5_flags(${_sw_target})
+  #add_executable(${_ugw_target} ${UGW_SOURCE_FILES})
+  
+  #add_dune_ug_flags(${_ugw_target})
 
-  set_property(TARGET ${_sw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
-  set_property(TARGET ${_ugw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
+  #set_property(TARGET ${_ugw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
 endforeach()

src/assemblers.hh

-#ifndef SRC_ASSEMBLERS_HH
-#define SRC_ASSEMBLERS_HH
-
-#include <dune/common/bitsetvector.hh>
-#include <dune/common/function.hh>
-#include <dune/istl/bcrsmatrix.hh>
-#include <dune/istl/bvector.hh>
-
-#include <dune/fufem/assemblers/assembler.hh>
-#pragma clang diagnostic push
-#pragma clang diagnostic ignored "-Wsign-compare"
-#include <dune/fufem/functionspacebases/p0basis.hh>
-#pragma clang diagnostic pop
-#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
-
-#include <dune/tectonic/globalfriction.hh>
-#include <dune/tectonic/globalfrictiondata.hh>
-
-#include "enums.hh"
-
-template <class GridView, int dimension> class MyAssembler {
-public:
-  using ScalarMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
-  using Matrix =
-      Dune::BCRSMatrix<Dune::FieldMatrix<double, dimension, dimension>>;
-  using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
-  using Vector = Dune::BlockVector<Dune::FieldVector<double, dimension>>;
-
-  using CellBasis = P0Basis<GridView, double>;
-  using VertexBasis = P1NodalBasis<GridView, double>;
-
-  CellBasis const cellBasis;
-  VertexBasis const vertexBasis;
-  GridView const &gridView;
-
-private:
-  using Grid = typename GridView::Grid;
-  using LocalVector = typename Vector::block_type;
-  using LocalScalarVector = typename ScalarVector::block_type;
-
-  using LocalCellBasis = typename CellBasis::LocalFiniteElement;
-  using LocalVertexBasis = typename VertexBasis::LocalFiniteElement;
-
-  Assembler<CellBasis, CellBasis> cellAssembler;
-  Assembler<VertexBasis, VertexBasis> vertexAssembler;
-
-public:
-  MyAssembler(GridView const &gridView);
-
-  void assembleFrictionalBoundaryMass(
-      BoundaryPatch<GridView> const &frictionalBoundary,
-      ScalarMatrix &frictionalBoundaryMass) const;
-
-  void assembleMass(Dune::VirtualFunction<
-                        LocalVector, LocalScalarVector> const &densityFunction,
-                    Matrix &M) const;
-
-  void assembleElasticity(double E, double nu, Matrix &A) const;
-
-  void assembleViscosity(
-      Dune::VirtualFunction<LocalVector, LocalScalarVector> const &
-          shearViscosity,
-      Dune::VirtualFunction<LocalVector, LocalScalarVector> const &
-          bulkViscosity,
-      Matrix &C) const;
-
-  void assembleBodyForce(
-      Dune::VirtualFunction<LocalVector, LocalVector> const &gravityField,
-      Vector &f) const;
-
-  void assembleNeumann(BoundaryPatch<GridView> const &neumannBoundary,
-                       Vector &f,
-                       Dune::VirtualFunction<double, double> const &neumann,
-                       double relativeTime) const;
-
-  void assembleWeightedNormalStress(
-      BoundaryPatch<GridView> const &frictionalBoundary,
-      ScalarVector &weightedNormalStress, double youngModulus,
-      double poissonRatio, Vector const &displacement) const;
-
-  std::shared_ptr<GlobalFriction<Matrix, Vector>> assembleFrictionNonlinearity(
-      Config::FrictionModel frictionModel,
-      BoundaryPatch<GridView> const &frictionalBoundary,
-      GlobalFrictionData<dimension> const &frictionInfo,
-      ScalarVector const &weightedNormalStress) const;
-
-  void assembleVonMisesStress(double youngModulus, double poissonRatio,
-                              Vector const &u, ScalarVector &stress) const;
-};
-#endif

src/explicitgrid.hh

-#ifndef SRC_EXPLICITGRID_HH
-#define SRC_EXPLICITGRID_HH
-
-#include "gridselector.hh"
-
-using GridView = Grid::LeafGridView;
-#endif

src/foam/CMakeLists.txt

+add_custom_target(tectonic_src_foam SOURCES
+  foam.cfg
+  foam-2D.cfg
+) 
+
+set(FOAM_SOURCE_FILES
+  ../../dune/tectonic/assemblers.cc
+  ../../dune/tectonic/data-structures/body/body.cc
+  ../../dune/tectonic/data-structures/network/levelcontactnetwork.cc
+  ../../dune/tectonic/data-structures/network/contactnetwork.cc
+  ../../dune/tectonic/data-structures/enumparser.cc
+  ../../dune/tectonic/factories/twoblocksfactory.cc
+  #../../dune/tectonic/io/vtk.cc
+  #../../dune/tectonic/io/hdf5/frictionalboundary-writer.cc
+  #../../dune/tectonic/io/hdf5/iteration-writer.cc
+  #../../dune/tectonic/io/hdf5/patchinfo-writer.cc
+  #../../dune/tectonic/io/hdf5/restart-io.cc
+  #../../dune/tectonic/io/hdf5/surface-writer.cc
+  #../../dune/tectonic/io/hdf5/time-writer.cc
+  ../../dune/tectonic/problem-data/grid/cuboidgeometry.cc
+  ../../dune/tectonic/problem-data/grid/mygrids.cc
+  ../../dune/tectonic/problem-data/grid/simplexmanager.cc
+  ../../dune/tectonic/spatial-solving/solverfactory.cc
+  ../../dune/tectonic/spatial-solving/fixedpointiterator.cc
+  ../../dune/tectonic/time-stepping/coupledtimestepper.cc
+  ../../dune/tectonic/time-stepping/adaptivetimestepper.cc
+  ../../dune/tectonic/time-stepping/rate.cc
+  ../../dune/tectonic/time-stepping/rate/rateupdater.cc
+  ../../dune/tectonic/time-stepping/state.cc
+  ../../dune/tectonic/time-stepping/uniformtimestepper.cc
+  foam.cc
+)
+
+foreach(_dim 2 3)
+  set(_foam_target foam-${_dim}D)
+
+  add_executable(${_foam_target} ${FOAM_SOURCE_FILES})
+
+  add_dune_ug_flags(${_foam_target})
+  add_dune_hdf5_flags(${_foam_target})
+
+  set_property(TARGET ${_foam_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
+  #set_property(TARGET ${_foam_target} APPEND PROPERTY COMPILE_DEFINITIONS "NEW_TNNMG_COMPUTE_ITERATES_DIRECTLY=1")
+endforeach()

src/foam/foam-2D.cfg

+# -*- mode:conf -*-
+[body0]
+smallestDiameter = 6.25e-2  # 2e-3 [m]
+
+[body1]
+smallestDiameter = 6.25e-2 # 2e-3 [m]
+
+[timeSteps]
+refinementTolerance = 1e-5 # 1e-5
+
+[u0.solver]
+tolerance         = 1e-8
+
+[a0.solver]
+tolerance         = 1e-8
+
+[v.solver]
+tolerance         = 1e-8
+
+[v.fpi]
+tolerance         = 1e-10
+
+[solver.tnnmg.preconditioner.basesolver]
+tolerance          = 1e-10
+
+[solver.tnnmg.preconditioner.patchsolver]
+tolerance          = 1e-10

src/foam/foam.cc

+#ifdef HAVE_CONFIG_H
+#include "config.h"
+#endif
+
+#ifdef HAVE_IPOPT
+#undef HAVE_IPOPT
+#endif
+
+#include <atomic>
+#include <cmath>
+#include <csignal>
+#include <exception>
+#include <fstream>
+#include <iostream>
+#include <iomanip>
+#include <memory>
+
+#include <dune/common/bitsetvector.hh>
+#include <dune/common/exceptions.hh>
+#include <dune/common/fmatrix.hh>
+#include <dune/common/function.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/parallel/mpihelper.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+
+#include <dune/grid/common/mcmgmapper.hh>
+#include <dune/istl/bcrsmatrix.hh>
+#include <dune/istl/bvector.hh>
+
+#include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/geometry/convexpolyhedron.hh>
+#include <dune/fufem/formatstring.hh>
+#include <dune/fufem/hdf5/file.hh>
+#include <dune/fufem/assemblers/transferoperatorassembler.hh>
+
+#include <dune/solvers/norms/energynorm.hh>
+#include <dune/solvers/solvers/loopsolver.hh>
+#include <dune/solvers/iterationsteps/blockgssteps.hh>
+
+#include <dune/contact/common/deformedcontinuacomplex.hh>
+#include <dune/contact/common/couplingpair.hh>
+#include <dune/contact/projections/normalprojection.hh>
+
+#include <dune/tnnmg/functionals/quadraticfunctional.hh>
+
+
+#include <dune/tectonic/assemblers.hh>
+#include <dune/tectonic/gridselector.hh>
+#include <dune/tectonic/explicitgrid.hh>
+#include <dune/tectonic/explicitvectors.hh>
+
+#include <dune/tectonic/data-structures/enumparser.hh>
+#include <dune/tectonic/data-structures/enums.hh>
+#include <dune/tectonic/data-structures/network/contactnetwork.hh>
+#include <dune/tectonic/data-structures/matrices.hh>
+#include <dune/tectonic/data-structures/program_state.hh>
+#include <dune/tectonic/data-structures/friction/globalfriction.hh>
+
+#include <dune/tectonic/factories/twoblocksfactory.hh>
+
+#include <dune/tectonic/io/io-handler.hh>
+#include <dune/tectonic/io/hdf5-writer.hh>
+#include <dune/tectonic/io/hdf5/restart-io.hh>
+#include <dune/tectonic/io/vtk.hh>
+
+#include <dune/tectonic/problem-data/bc.hh>
+#include <dune/tectonic/problem-data/mybody.hh>
+#include <dune/tectonic/problem-data/grid/mygrids.hh>
+
+#include <dune/tectonic/spatial-solving/tnnmg/functional.hh>
+//#include <dune/tectonic/spatial-solving/preconditioners/multilevelpatchpreconditioner.hh>
+#include <dune/tectonic/spatial-solving/tnnmg/localbisectionsolver.hh>
+#include <dune/tectonic/spatial-solving/solverfactory.hh>
+#include <dune/tectonic/spatial-solving/makelinearsolver.hh>
+
+#include <dune/tectonic/time-stepping/uniformtimestepper.hh>
+#include <dune/tectonic/time-stepping/rate.hh>
+#include <dune/tectonic/time-stepping/state.hh>
+#include <dune/tectonic/time-stepping/updaters.hh>
+
+#include <dune/tectonic/utils/debugutils.hh>
+#include <dune/tectonic/utils/diameter.hh>
+#include <dune/tectonic/utils/geocoordinate.hh>
+
+#include <dune/matrix-vector/axpy.hh>
+
+#include <dune/contact/assemblers/nbodyassembler.hh>
+#include "dune/tectonic/spatial-solving/tnnmg/zerononlinearity.hh"
+
+
+
+// for getcwd
+#include <unistd.h>
+
+//#include <tbb/tbb.h> //TODO multi threading preconditioner?
+//#include <pthread.h>
+
+#include <dune/tectonic/utils/reductionfactors.hh>
+std::vector<std::vector<double>> allReductionFactors;
+
+size_t const dims = MY_DIM;
+const std::string sourcePath = "/home/joscha/software/dune/dune-tectonic/src/foam/";
+
+Dune::ParameterTree getParameters(int argc, char *argv[]) {
+  Dune::ParameterTree parset;
+  Dune::ParameterTreeParser::readINITree(sourcePath + "foam.cfg", parset);
+  Dune::ParameterTreeParser::readINITree(
+      Dune::Fufem::formatString(sourcePath + "foam-%dD.cfg", dims), parset);
+  Dune::ParameterTreeParser::readOptions(argc, argv, parset);
+  return parset;
+}
+
+static std::atomic<bool> terminationRequested(false);
+void handleSignal(int signum) { terminationRequested = true; }
+
+int main(int argc, char *argv[]) {
+  using BlocksFactory = TwoBlocksFactory<Grid, Vector>;
+  using ContactNetwork = typename BlocksFactory::ContactNetwork;
+  using MyProgramState = ProgramState<Vector, ScalarVector>;
+  using Assembler = MyAssembler<DefLeafGridView, dims>;
+
+  using IOHandler = IOHandler<Assembler, ContactNetwork, MyProgramState>;
+  std::unique_ptr<IOHandler> ioHandler;
+
+  try {
+    Dune::MPIHelper::instance(argc, argv);
+
+    char buffer[256];
+    char *val = getcwd(buffer, sizeof(buffer));
+    if (val) {
+        std::cout << buffer << std::endl;
+        std::cout << argv[0] << std::endl;
+    }
+
+    auto const parset = getParameters(argc, argv);
+
+    auto outPath = std::filesystem::current_path();
+    outPath +=  "/output/" + parset.get<std::string>("general.outPath");
+    if (!std::filesystem::is_directory(outPath))
+        std::filesystem::create_directories(outPath);
+
+    const auto copyOptions = std::filesystem::copy_options::overwrite_existing;
+    std::filesystem::copy(sourcePath + "foam.cfg", outPath, copyOptions);
+    std::filesystem::copy(Dune::Fufem::formatString(sourcePath + "foam-%dD.cfg", dims), outPath, copyOptions);
+    std::filesystem::current_path(outPath);
+
+    std::ofstream out("foam.log");
+    std::streambuf *coutbuf = std::cout.rdbuf(); //save old buffer
+    std::cout.rdbuf(out.rdbuf()); //redirect std::cout to log.txt
+
+
+    using field_type = Matrix::field_type;
+
+    // ----------------------
+    // set up contact network
+    // ----------------------
+    BlocksFactory blocksFactory(parset);
+    blocksFactory.build();
+
+    ContactNetwork& contactNetwork = blocksFactory.contactNetwork();
+
+    const size_t bodyCount = contactNetwork.nBodies();
+
+    for (size_t i=0; i<contactNetwork.nLevels(); i++) {
+        // printDofLocation(contactNetwork.body(i)->gridView());
+
+        //Vector def(contactNetwork.deformedGrids()[i]->size(dims));
+        //def = 1;
+        //deformedGridComplex.setDeformation(def, i);
+
+        const auto& level = *contactNetwork.level(i);
+
+        //for (size_t j=0; j<level.nBodies(); j++) {
+            //writeToVTK(level.body(j)->gridView(), "../debug_print/bodies/", "body_" + std::to_string(j) + "_level_" + std::to_string(i));
+        //}
+    }
+
+    //for (size_t i=0; i<bodyCount; i++) {
+        //writeToVTK(contactNetwork.body(i)->gridView(), "../debug_print/bodies/", "body_" + std::to_string(i) + "_leaf");
+    //}
+
+    // ----------------------------
+    // assemble contactNetwork
+    // ----------------------------
+    contactNetwork.assemble();
+
+    //printMortarBasis<Vector>(contactNetwork.nBodyAssembler());
+
+    // -----------------
+    // init input/output
+    // -----------------
+    std::vector<size_t> nVertices(bodyCount);
+    for (size_t i=0; i<bodyCount; i++) {
+        nVertices[i] = contactNetwork.body(i)->nVertices();
+    }
+    print(nVertices, "#dofs: ");
+
+    MyProgramState programState(nVertices);
+
+    ioHandler = std::make_unique<IOHandler>(parset.sub("io"), contactNetwork);
+
+    bool restartRead = ioHandler->read(programState);
+    if (!restartRead) {
+      programState.setupInitialConditions(parset, contactNetwork);
+    }
+
+
+    const auto& nBodyAssembler = contactNetwork.nBodyAssembler();
+   /* auto linearSolver = makeLinearSolver<ContactNetwork, Vector>(parset, contactNetwork);
+    auto nonlinearity = ZeroNonlinearity();
+
+    // Solving a linear problem with a multigrid solver
+    auto const solveLinearProblem = [&](
+        const BitVector& _dirichletNodes, const std::vector<std::shared_ptr<Matrix>>& _matrices,
+        const std::vector<Vector>& _rhs, std::vector<Vector>& _x,
+        Dune::ParameterTree const &_localParset) {
+
+        Vector totalX;
+        nBodyAssembler.nodalToTransformed(_x, totalX);
+
+        FunctionalFactory<std::decay_t<decltype(nBodyAssembler)>, decltype(nonlinearity), Matrix, Vector> functionalFactory(nBodyAssembler, nonlinearity, _matrices, _rhs);
+        functionalFactory.build();
+        auto functional = functionalFactory.functional();
+
+        NonlinearSolver<std::remove_reference_t<decltype(*functional)>, BitVector> solver(parset.sub("solver.tnnmg"), linearSolver, functional, _dirichletNodes);
+        solver.solve(_localParset, totalX);
+
+        nBodyAssembler.postprocess(totalX, _x);
+    };
+
+    programState.timeStep = parset.get<size_t>("initialTime.timeStep");
+    programState.relativeTime = parset.get<double>("initialTime.relativeTime");
+    programState.relativeTau = parset.get<double>("initialTime.relativeTau");
+
+    std::vector<Vector> ell0(bodyCount);
+    for (size_t i=0; i<bodyCount; i++) {
+      programState.u[i] = 0.0;
+      programState.v[i] = 0.0;
+      programState.a[i] = 0.0;
+
+      ell0[i].resize(programState.u[i].size());
+      ell0[i] = 0.0;
+
+      contactNetwork.body(i)->externalForce()(programState.relativeTime, ell0[i]);
+    }
+    */
+    // 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
+    BitVector totalDirichletNodes;
+    contactNetwork.totalNodes("dirichlet", totalDirichletNodes);
+
+    using BoundaryNodes = typename ContactNetwork::BoundaryNodes;
+    BoundaryNodes dirichletNodes;
+    contactNetwork.boundaryNodes("dirichlet", dirichletNodes);
+
+    size_t dof = 0;
+    for (size_t i=0; i<bodyCount; i++) {
+        const auto& body = *contactNetwork.body(i);
+
+        if (body.data()->getYoungModulus() == 0.0) {
+            for (size_t j=0; j<body.nVertices(); j++) {
+                totalDirichletNodes[dof] = true;
+                dof++;
+            }
+        } else {
+            dof += body.nVertices();
+        }
+    }
+
+
+
+    std::vector<const Dune::BitSetVector<1>*> frictionNodes;
+    contactNetwork.frictionNodes(frictionNodes);
+
+    /*
+    std::cout << "solving linear problem for u..." << std::endl;
+
+    {
+        int nVertices = contactNetwork.body(1)->nVertices();
+        BitVector uDirichletNodes(nVertices, false);
+
+        int idx=0;
+        const auto& body = contactNetwork.body(1);
+        using LeafBody = typename ContactNetwork::LeafBody;
+        std::vector<std::shared_ptr<typename LeafBody::BoundaryCondition>> boundaryConditions;
+        body->boundaryConditions("dirichlet", boundaryConditions);
+
+        if (boundaryConditions.size()>0) {
+            const int idxBackup = idx;
+            for (size_t bc = 0; bc<boundaryConditions.size(); bc++) {
+                const auto& nodes = boundaryConditions[bc]->boundaryNodes();
+                for (size_t j=0; j<nodes->size(); j++, idx++)
+                    for (int k=0; k<dims; k++)
+                        uDirichletNodes[idx][k] = uDirichletNodes[idx][k] or (*nodes)[j][k];
+                idx = (bc==boundaryConditions.size()-1 ? idx : idxBackup);
+            }
+        }
+
+        for (auto i=0; i<nVertices; i++) {
+            if ((*frictionNodes[1])[i][0]) {
+                uDirichletNodes[i][1] = true;
+            }
+        }
+
+        using Norm = EnergyNorm<Matrix, Vector>;
+        using LinearSolver = typename Dune::Solvers::LoopSolver<Vector>;
+
+        // transfer operators need to be recomputed on change due to setDeformation()
+        using TransferOperator = CompressedMultigridTransfer<Vector>;
+        using TransferOperators = std::vector<std::shared_ptr<TransferOperator>>;
+
+        const auto& grid = contactNetwork.body(1)->grid();
+        TransferOperators transfer(grid->maxLevel());
+        for (size_t i = 0; i < transfer.size(); ++i)
+        {
+            // create transfer operators from level i to i+1 (note that this will only work for either uniformly refined grids or adaptive grids with RefinementType=COPY)
+            auto t = std::make_shared<TransferOperator>();
+            t->setup(*grid, i, i+1);
+            transfer[i] = t;
+        }
+
+        auto linearMultigridStep = std::make_shared<Dune::Solvers::MultigridStep<Matrix, Vector> >(
+                    *contactNetwork.matrices().elasticity[1],
+                    programState.u[1],
+                    ell0[1]);
+        linearMultigridStep->setIgnore(uDirichletNodes);
+        linearMultigridStep->setMGType(1, 3, 3);
+        linearMultigridStep->setSmoother(TruncatedBlockGSStep<Matrix, Vector>());
+        linearMultigridStep->setTransferOperators(transfer);
+
+        const auto& solverParset = parset.sub("u0.solver");
+        Dune::Solvers::LoopSolver<Vector> solver(linearMultigridStep, solverParset.get<size_t>("maximumIterations"),
+            solverParset.get<double>("tolerance"), Norm(*linearMultigridStep),
+            solverParset.get<Solver::VerbosityMode>("verbosity")); // absolute error
+
+
+        solver.preprocess();
+        solver.solve();
+    }
+
+    //print(u, "initial u:");
+
+    // Initial acceleration: Computed in agreement with Ma = ell0 - Au
+    // (without Dirichlet constraints), again assuming dPhi(v = 0) = 0
+    std::vector<Vector> accelerationRHS = ell0;
+    for (size_t i=0; i<bodyCount; i++) {
+      const auto& body = contactNetwork.body(i);
+      Dune::MatrixVector::subtractProduct(accelerationRHS[i], *body->matrices().elasticity, programState.u[i]);
+    }
+
+    std::cout << "solving linear problem for a..." << std::endl;
+
+    BitVector noNodes(totalDirichletNodes.size(), false);
+    solveLinearProblem(noNodes, contactNetwork.matrices().mass, accelerationRHS, programState.a,
+                       parset.sub("a0.solver"));
+
+    // setup initial conditions in program state
+    programState.alpha[0] = 0;
+    programState.alpha[1] = parset.get<double>("boundary.friction.initialAlpha");
+
+    for (size_t i=0; i<contactNetwork.nCouplings(); i++) {
+      const auto& coupling = contactNetwork.coupling(i);
+      const auto& contactCoupling = contactNetwork.nBodyAssembler().getContactCouplings()[i];
+
+      const auto nonmortarIdx = coupling->gridIdx_[0];
+      const auto& body = contactNetwork.body(nonmortarIdx);
+
+      ScalarVector frictionBoundaryStress(programState.weightedNormalStress[nonmortarIdx].size());
+
+      body->assembler()->assembleWeightedNormalStress(
+        contactCoupling->nonmortarBoundary(), frictionBoundaryStress, body->data()->getYoungModulus(),
+        body->data()->getPoissonRatio(), programState.u[nonmortarIdx]);
+
+      programState.weightedNormalStress[nonmortarIdx] += frictionBoundaryStress;
+    } */
+
+    for (size_t i=0; i<bodyCount; i++) {
+      contactNetwork.body(i)->setDeformation(programState.u[i]);
+    }
+    nBodyAssembler.assembleTransferOperator();
+    nBodyAssembler.assembleObstacle();
+
+    // ------------------------
+    // assemble global friction
+    // ------------------------
+    //print(programState.weightedNormalStress, "weightedNormalStress");
+
+    contactNetwork.assembleFriction(parset.get<Config::FrictionModel>("boundary.friction.frictionModel"), programState.weightedNormalStress);
+
+    auto& globalFriction = contactNetwork.globalFriction();
+    globalFriction.updateAlpha(programState.alpha);
+
+    IterationRegister iterationCount;
+
+    ioHandler->write(programState, contactNetwork, globalFriction, iterationCount, true);
+
+    // -------------------
+    // Set up TNNMG solver
+    // -------------------
+
+    //using Functional = Functional<Matrix&, Vector&, ZeroNonlinearity&, Vector&, Vector&, field_type>;
+    using Functional = Functional<Matrix&, Vector&, GlobalFriction<Matrix, Vector>&, Vector&, Vector&, field_type>;
+    using NonlinearFactory = SolverFactory<Functional, BitVector>;
+
+    using BoundaryFunctions = typename ContactNetwork::BoundaryFunctions;
+    using Updaters = Updaters<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>,
+                               StateUpdater<ScalarVector, Vector>>;
+
+    BoundaryFunctions velocityDirichletFunctions;
+    contactNetwork.boundaryFunctions("dirichlet", velocityDirichletFunctions);
+
+    /*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));
+        }
+    }*/
+
+
+    /*for (size_t i=0; i<frictionNodes.size(); i++) {
+        print(*frictionNodes[i], "frictionNodes_body_" + std::to_string(i));
+    }*/
+
+    //DUNE_THROW(Dune::Exception, "Just need to stop here!");
+
+    Updaters current(
+        initRateUpdater(
+            parset.get<Config::scheme>("timeSteps.scheme"),
+            velocityDirichletFunctions,
+            dirichletNodes,
+            contactNetwork.matrices(),
+            programState.u,
+            programState.v,
+            programState.a),
+        initStateUpdater<ScalarVector, Vector>(
+            parset.get<Config::stateModel>("boundary.friction.stateModel"),
+            programState.alpha,
+            nBodyAssembler.getContactCouplings(),
+            contactNetwork.couplings())
+            );
+
+
+    auto const refinementTolerance = parset.get<double>("timeSteps.refinementTolerance");
+
+    const auto& stateEnergyNorms = contactNetwork.stateEnergyNorms();
+
+    std::signal(SIGXCPU, handleSignal);
+    std::signal(SIGINT, handleSignal);
+    std::signal(SIGTERM, handleSignal);
+
+/*
+    // 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<Preconditioner::Mode>("mode"));
+    preconditioner.setPatchSolver(patchSolver);
+    preconditioner.setPatchDepth(preconditionerParset.get<size_t>("patchDepth"));
+*/
+    // set adaptive time stepper
+    typename ContactNetwork::ExternalForces externalForces;
+    contactNetwork.externalForces(externalForces);
+
+    auto&& noFriction = ZeroNonlinearity();
+
+    StepBase<NonlinearFactory, std::decay_t<decltype(contactNetwork)>, Updaters, std::decay_t<decltype(stateEnergyNorms)>>
+        stepBase(parset, contactNetwork, totalDirichletNodes, globalFriction, frictionNodes,
+                 externalForces, stateEnergyNorms);
+
+    UniformTimeStepper<NonlinearFactory, std::decay_t<decltype(contactNetwork)>, Updaters, std::decay_t<decltype(stateEnergyNorms)>>
+        timeStepper(stepBase, contactNetwork, current,
+                            programState.relativeTime, programState.relativeTau);
+
+    size_t timeSteps = std::round(parset.get<double>("timeSteps.timeSteps"));
+
+    while (!timeStepper.reachedEnd()) {
+      programState.timeStep++;
+
+      //preconditioner.build();
+      iterationCount = timeStepper.advance();
+
+      programState.relativeTime = timeStepper.relativeTime_;
+      programState.relativeTau = timeStepper.relativeTau_;
+      current.rate_->extractDisplacement(programState.u);
+      current.rate_->extractVelocity(programState.v);
+      current.rate_->extractAcceleration(programState.a);
+      current.state_->extractAlpha(programState.alpha);
+
+      globalFriction.updateAlpha(programState.alpha);
+
+      /*print(programState.u, "current u:");
+      print(programState.v, "current v:");
+      print(programState.a, "current a:");
+      print(programState.alpha, "current alpha:");*/
+
+      //using Vector = typename ProgramState::Vector;
+      /*Vector mortarV;
+      contactNetwork.nBodyAssembler().nodalToTransformed(programState.v, mortarV);
+      printRegularityTruncation(globalFriction, mortarV);*/
+
+      contactNetwork.setDeformation(programState.u);
+
+      ioHandler->write(programState, contactNetwork, globalFriction, iterationCount, false);
+
+
+
+      if (programState.timeStep==timeSteps) {
+        std::cout << "limit of timeSteps reached!" << std::endl;
+        break; // TODO remove after debugging
+      }
+
+      if (terminationRequested) {
+        std::cerr << "Terminating prematurely" << std::endl;
+        break;
+      }
+
+
+    }
+
+
+    std::cout.rdbuf(coutbuf); //reset to standard output again
+
+  } catch (Dune::Exception &e) {
+    Dune::derr << "Dune reported error: " << e << std::endl;
+  } catch (std::exception &e) {
+    std::cerr << "Standard exception: " << e.what() << std::endl;
+  } catch (...) {
+
+  }
+
+  ioHandler.reset();
+}

src/foam/foam.cfg

+# -*- mode:conf -*-
+[general]
+outPath = friction-0.6 # output written to ./output/outPath
+gravity         = 9.81     # [m/s^2]
+
+[body0]
+length          = 6.0      # [m]
+height          = 1.0     # [m]
+bulkModulus     = 0.0 #4.12e9   # [Pa] #2190
+poissonRatio    = 0.0 #0.3     # [1]  #0.11
+[body0.elastic]
+density         = 5e3      # [kg/m^3] #750
+shearViscosity  = 0.0     # [Pas]
+bulkViscosity   = 0.0     # [Pas]
+[body0.viscoelastic]
+density         = 5e3      # [kg/m^3] #750
+shearViscosity  = 0.0     # [Pas]
+bulkViscosity   = 0.0     # [Pas]
+
+[body1]
+length          = 5.00     # [m]
+height          = 1.00     # [m]
+bulkModulus     = 4.12e7    # [Pa]
+poissonRatio    = 0.3     # [1]
+[body1.elastic]
+density         = 5e3      # [kg/m^3]
+shearViscosity  = 0.0     # [Pas]
+bulkViscosity   = 0.0     # [Pas]
+[body1.viscoelastic]
+density         = 5e3      # [kg/m^3]
+shearViscosity  = 0.0     # [Pas]
+bulkViscosity   = 0.0     # [Pas]
+
+
+[boundary.friction]
+C               = 0.0       # [Pa]
+mu0             = 0.6      # [ ]
+V0              = 1e-6     # [m/s]
+L               = 1e-5  # [m]
+initialAlpha    = -10.0  # [ ]
+stateModel      = AgeingLaw
+frictionModel   = Truncated #Tresca #None #Truncated #Regularised
+[boundary.friction.weakening]
+a               = 0.010    # [ ]
+b               = 0.015    # [ ]
+[boundary.friction.strengthening]
+a               = 0.010    # [ ]
+b               = 0.015    # [ ]
+
+
+[boundary.neumann]
+sigmaN          = 0.0      # [Pa]
+
+[boundary.dirichlet]
+finalVelocity   = 2e-3 #2e-4     # [m/s]
+
+[io]
+data.write      = true
+printProgress   = true
+restarts.first  = 0
+restarts.spacing= 50
+restarts.write  = true #true
+vtk.write       = true
+
+[problem]
+finalTime       = 15     # [s] #1000
+bodyCount       = 2
+
+[initialTime]
+timeStep = 0
+relativeTime = 0.0
+relativeTau = 5e-4 # 1e-6
+
+[timeSteps]
+scheme = newmark
+timeSteps = 2e3
+
+[u0.solver]
+maximumIterations = 100
+verbosity         = full
+
+[a0.solver]
+maximumIterations = 100
+verbosity         = full
+
+[v.solver]
+maximumIterations = 100
+verbosity         = quiet
+
+[v.fpi]
+maximumIterations = 10000
+lambda            = 0.5
+
+[solver.tnnmg.preconditioner]
+mode         = additive
+patchDepth   = 1
+maximumIterations = 2
+verbosity         = quiet
+[solver.tnnmg.preconditioner.patchsolver]
+maximumIterations = 100
+verbosity         = quiet
+[solver.tnnmg.preconditioner.basesolver]
+maximumIterations = 10000
+verbosity         = quiet
+
+[solver.tnnmg.main]
+pre   = 1
+multi = 5 # number of multigrid steps
+post  = 0
+
+
+