src/one-body-problem-data/mygrid_tmpl.cc

-#ifndef MY_DIM
-#error MY_DIM unset
-#endif
-
-#include "../explicitgrid.hh"
-#include "../explicitvectors.hh"
-
-template class GridConstructor<Grid>;
-
-template struct MyFaces<GridView>;
-
-template MyFaces<GridView> GridConstructor<Grid>::constructFaces(
-    GridView const &gridView);
-
-template void refine<Grid, LocalVector>(
-    Grid &grid, ConvexPolyhedron<LocalVector> const &weakPatch,
-    double smallestDiameter);

src/one-body-problem-data/weakpatch.hh

-#ifndef SRC_ONE_BODY_PROBLEM_DATA_WEAKPATCH_HH
-#define SRC_ONE_BODY_PROBLEM_DATA_WEAKPATCH_HH
-
-template <class LocalVector>
-ConvexPolyhedron<LocalVector> getWeakPatch(Dune::ParameterTree const &parset) {
-  ConvexPolyhedron<LocalVector> weakPatch;
-#if MY_DIM == 3
-  weakPatch.vertices.resize(4);
-  weakPatch.vertices[0] = weakPatch.vertices[2] = MyGeometry::X;
-  weakPatch.vertices[1] = weakPatch.vertices[3] = MyGeometry::Y;
-  for (size_t k = 0; k < 2; ++k) {
-    weakPatch.vertices[k][2] = -MyGeometry::depth / 2.0;
-    weakPatch.vertices[k + 2][2] = MyGeometry::depth / 2.0;
-  }
-  switch (parset.get<Config::PatchType>("patchType")) {
-    case Config::Rectangular:
-      break;
-    case Config::Trapezoidal:
-      weakPatch.vertices[1][0] += 0.05 * MyGeometry::lengthScale;
-      weakPatch.vertices[3][0] -= 0.05 * MyGeometry::lengthScale;
-      break;
-    default:
-      assert(false);
-  }
-#else
-  weakPatch.vertices.resize(2);
-  weakPatch.vertices[0] = MyGeometry::X;
-  weakPatch.vertices[1] = MyGeometry::Y;
-#endif
-  return weakPatch;
-};
-#endif

src/one-body-problem.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 <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/formatstring.hh>
-
-#include <dune/solvers/norms/energynorm.hh>
-#include <dune/solvers/solvers/loopsolver.hh>
-#include <dune/solvers/solvers/solver.hh>
-#include <dune/tnnmg/problem-classes/convexproblem.hh>
-
-#include <dune/tectonic/geocoordinate.hh>
-#include <dune/tectonic/myblockproblem.hh>
-#include <dune/tectonic/globalfriction.hh>
-#include <dune/fufem/hdf5/file.hh>
-
-#include "assemblers.hh"
-#include "diameter.hh"
-#include "enumparser.hh"
-#include "enums.hh"
-#include "gridselector.hh"
-#include "hdf5-writer.hh"
-#include "hdf5/restart-io.hh"
-#include "matrices.hh"
-#include "program_state.hh"
-#include "one-body-problem-data/bc.hh"
-#include "one-body-problem-data/mybody.hh"
-#include "one-body-problem-data/mygeometry.hh"
-#include "one-body-problem-data/myglobalfrictiondata.hh"
-#include "one-body-problem-data/mygrid.hh"
-#include "one-body-problem-data/weakpatch.hh"
-#include "spatial-solving/solverfactory.hh"
-#include "time-stepping/adaptivetimestepper.hh"
-#include "time-stepping/rate.hh"
-#include "time-stepping/state.hh"
-#include "time-stepping/updaters.hh"
-#include "vtk.hh"
-
-size_t const dims = MY_DIM;
-
-Dune::ParameterTree getParameters(int argc, char *argv[]) {
-  Dune::ParameterTree parset;
-  Dune::ParameterTreeParser::readINITree("one-body-problem.cfg", parset);
-  Dune::ParameterTreeParser::readINITree(
-      Dune::Fufem::formatString("one-body-problem-%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[]) {
-  try {
-    Dune::MPIHelper::instance(argc, argv);
-    auto const parset = getParameters(argc, argv);
-
-    MyGeometry::render();
-    MyGeometry::write();
-
-    using GridView = Grid::LeafGridView;
-    using MyAssembler = MyAssembler<GridView, dims>;
-    using Matrix = MyAssembler::Matrix;
-    using Vector = MyAssembler::Vector;
-    using LocalVector = Vector::block_type;
-    using ScalarMatrix = MyAssembler::ScalarMatrix;
-    using ScalarVector = MyAssembler::ScalarVector;
-
-    auto const weakPatch =
-        getWeakPatch<LocalVector>(parset.sub("boundary.friction.weakening"));
-
-    // {{{ Set up grid
-    GridConstructor<Grid> gridConstructor;
-    auto grid = gridConstructor.getGrid();
-
-    refine(*grid, weakPatch,
-           parset.get<double>("boundary.friction.smallestDiameter"));
-
-    double minDiameter = std::numeric_limits<double>::infinity();
-    double maxDiameter = 0.0;
-    for (auto &&e : elements(grid->leafGridView())) {
-      auto const geometry = e.geometry();
-      auto const diam = diameter(geometry);
-      minDiameter = std::min(minDiameter, diam);
-      maxDiameter = std::max(maxDiameter, diam);
-    }
-    std::cout << "min diameter: " << minDiameter << std::endl;
-    std::cout << "max diameter: " << maxDiameter << std::endl;
-
-    auto const leafView = grid->leafGridView();
-    auto const leafVertexCount = leafView.size(dims);
-
-    std::cout << "Number of DOFs: " << leafVertexCount << std::endl;
-
-    auto myFaces = gridConstructor.constructFaces(leafView);
-
-    BoundaryPatch<GridView> const neumannBoundary(leafView);
-    BoundaryPatch<GridView> const &frictionalBoundary = myFaces.lower;
-    BoundaryPatch<GridView> const &surface = myFaces.upper;
-
-    // Dirichlet Boundary
-    Dune::BitSetVector<dims> noNodes(leafVertexCount);
-    Dune::BitSetVector<dims> dirichletNodes(leafVertexCount);
-    for (size_t i = 0; i < leafVertexCount; ++i) {
-      if (myFaces.right.containsVertex(i))
-        dirichletNodes[i][0] = true;
-
-      if (myFaces.lower.containsVertex(i))
-        dirichletNodes[i][1] = true;
-
-#if MY_DIM == 3
-      if (myFaces.front.containsVertex(i) || myFaces.back.containsVertex(i))
-        dirichletNodes[i][2] = true;
-#endif
-    }
-
-    // Set up functions for time-dependent boundary conditions
-    using Function = Dune::VirtualFunction<double, double>;
-    Function const &velocityDirichletFunction = VelocityDirichletCondition();
-    Function const &neumannFunction = NeumannCondition();
-
-    MyAssembler const myAssembler(leafView);
-
-    MyBody<dims> const body(parset);
-
-    Matrices<Matrix> matrices;
-    myAssembler.assembleElasticity(body.getYoungModulus(),
-                                   body.getPoissonRatio(), matrices.elasticity);
-    myAssembler.assembleViscosity(body.getShearViscosityField(),
-                                  body.getBulkViscosityField(),
-                                  matrices.damping);
-    myAssembler.assembleMass(body.getDensityField(), matrices.mass);
-
-    ScalarMatrix relativeFrictionalBoundaryMass;
-    myAssembler.assembleFrictionalBoundaryMass(frictionalBoundary,
-                                               relativeFrictionalBoundaryMass);
-    relativeFrictionalBoundaryMass /= frictionalBoundary.area();
-    EnergyNorm<ScalarMatrix, ScalarVector> const stateEnergyNorm(
-        relativeFrictionalBoundaryMass);
-
-    // Assemble forces
-    Vector gravityFunctional;
-    myAssembler.assembleBodyForce(body.getGravityField(), gravityFunctional);
-
-    // Problem formulation: right-hand side
-    std::function<void(double, Vector &)> computeExternalForces =
-        [&](double _relativeTime, Vector &_ell) {
-          myAssembler.assembleNeumann(neumannBoundary, _ell, neumannFunction,
-                                      _relativeTime);
-          _ell += gravityFunctional;
-        };
-
-    using MyProgramState = ProgramState<Vector, ScalarVector>;
-    MyProgramState programState(leafVertexCount);
-    auto const firstRestart = parset.get<size_t>("io.restarts.first");
-    auto const restartSpacing = parset.get<size_t>("io.restarts.spacing");
-    auto const writeRestarts = parset.get<bool>("io.restarts.write");
-    auto const writeData = parset.get<bool>("io.data.write");
-    bool const handleRestarts = writeRestarts or firstRestart > 0;
-
-    auto dataFile =
-        writeData ? std::make_unique<HDF5::File>("output.h5") : nullptr;
-    auto restartFile = handleRestarts
-                           ? std::make_unique<HDF5::File>(
-                                 "restarts.h5",
-                                 writeRestarts ? HDF5::Access::READWRITE
-                                               : HDF5::Access::READONLY)
-                           : nullptr;
-    auto restartIO = handleRestarts
-                         ? std::make_unique<RestartIO<MyProgramState>>(
-                               *restartFile, leafVertexCount)
-                         : nullptr;
-
-    if (firstRestart > 0) // automatically adjusts the time and timestep
-      restartIO->read(firstRestart, programState);
-    else
-      programState.setupInitialConditions(parset, computeExternalForces,
-                                          matrices, myAssembler, dirichletNodes,
-                                          noNodes, frictionalBoundary, body);
-
-    MyGlobalFrictionData<LocalVector> frictionInfo(
-        parset.sub("boundary.friction"), weakPatch);
-    auto myGlobalFriction = myAssembler.assembleFrictionNonlinearity(
-        parset.get<Config::FrictionModel>("boundary.friction.frictionModel"),
-        frictionalBoundary, frictionInfo, programState.weightedNormalStress);
-    myGlobalFriction->updateAlpha(programState.alpha);
-
-    Vector vertexCoordinates(leafVertexCount);
-    {
-      Dune::MultipleCodimMultipleGeomTypeMapper<
-          GridView, Dune::MCMGVertexLayout> const vertexMapper(leafView);
-      for (auto &&v : vertices(leafView))
-        vertexCoordinates[vertexMapper.index(v)] = geoToPoint(v.geometry());
-    }
-
-    using MyVertexBasis = typename MyAssembler::VertexBasis;
-    auto dataWriter =
-        writeData ? std::make_unique<
-                        HDF5Writer<MyProgramState, MyVertexBasis, GridView>>(
-                        *dataFile, vertexCoordinates, myAssembler.vertexBasis,
-                        surface, frictionalBoundary, weakPatch)
-                  : nullptr;
-    MyVTKWriter<MyVertexBasis, typename MyAssembler::CellBasis> const vtkWriter(
-        myAssembler.cellBasis, myAssembler.vertexBasis, "obs");
-
-    IterationRegister iterationCount;
-    auto const report = [&](bool initial = false) {
-      if (writeData) {
-        dataWriter->reportSolution(programState, *myGlobalFriction);
-        if (!initial)
-          dataWriter->reportIterations(programState, iterationCount);
-        dataFile->flush();
-      }
-
-      if (writeRestarts and !initial and
-          programState.timeStep % restartSpacing == 0) {
-        restartIO->write(programState);
-        restartFile->flush();
-      }
-
-      if (parset.get<bool>("io.printProgress"))
-        std::cout << "timeStep = " << std::setw(6) << programState.timeStep
-                  << ", time = " << std::setw(12) << programState.relativeTime
-                  << ", tau = " << std::setw(12) << programState.relativeTau
-                  << std::endl;
-
-      if (parset.get<bool>("io.vtk.write")) {
-        ScalarVector stress;
-        myAssembler.assembleVonMisesStress(body.getYoungModulus(),
-                                           body.getPoissonRatio(),
-                                           programState.u, stress);
-        vtkWriter.write(programState.timeStep, programState.u, programState.v,
-                        programState.alpha, stress);
-      }
-    };
-    report(true);
-
-    // Set up TNNMG solver
-    using NonlinearFactory = SolverFactory<
-        dims,
-        MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
-        Grid>;
-    NonlinearFactory factory(parset.sub("solver.tnnmg"), *grid, dirichletNodes);
-
-    using MyUpdater = Updaters<RateUpdater<Vector, Matrix, Function, dims>,
-                               StateUpdater<ScalarVector, Vector>>;
-    MyUpdater current(
-        initRateUpdater(parset.get<Config::scheme>("timeSteps.scheme"),
-                        velocityDirichletFunction, dirichletNodes, matrices,
-                        programState.u, programState.v, programState.a),
-        initStateUpdater<ScalarVector, Vector>(
-            parset.get<Config::stateModel>("boundary.friction.stateModel"),
-            programState.alpha, *frictionalBoundary.getVertices(),
-            parset.get<double>("boundary.friction.L"),
-            parset.get<double>("boundary.friction.V0")));
-
-    auto const refinementTolerance =
-        parset.get<double>("timeSteps.refinementTolerance");
-    auto const mustRefine = [&](MyUpdater &coarseUpdater,
-                                MyUpdater &fineUpdater) {
-      ScalarVector coarseAlpha;
-      coarseUpdater.state_->extractAlpha(coarseAlpha);
-
-      ScalarVector fineAlpha;
-      fineUpdater.state_->extractAlpha(fineAlpha);
-
-      return stateEnergyNorm.diff(fineAlpha, coarseAlpha) > refinementTolerance;
-    };
-
-    std::signal(SIGXCPU, handleSignal);
-    std::signal(SIGINT, handleSignal);
-    std::signal(SIGTERM, handleSignal);
-
-    AdaptiveTimeStepper<NonlinearFactory, MyUpdater,
-                        EnergyNorm<ScalarMatrix, ScalarVector>>
-        adaptiveTimeStepper(factory, parset, myGlobalFriction, current,
-                            programState.relativeTime, programState.relativeTau,
-                            computeExternalForces, stateEnergyNorm, mustRefine);
-    while (!adaptiveTimeStepper.reachedEnd()) {
-      programState.timeStep++;
-
-      iterationCount = adaptiveTimeStepper.advance();
-
-      programState.relativeTime = adaptiveTimeStepper.relativeTime_;
-      programState.relativeTau = adaptiveTimeStepper.relativeTau_;
-      current.rate_->extractDisplacement(programState.u);
-      current.rate_->extractVelocity(programState.v);
-      current.rate_->extractAcceleration(programState.a);
-      current.state_->extractAlpha(programState.alpha);
-
-      report();
-
-      if (terminationRequested) {
-        std::cerr << "Terminating prematurely" << std::endl;
-        break;
-      }
-    }
-  } catch (Dune::Exception &e) {
-    Dune::derr << "Dune reported error: " << e << std::endl;
-  } catch (std::exception &e) {
-    std::cerr << "Standard exception: " << e.what() << std::endl;
-  }
-}

src/program_state.hh

-#ifndef SRC_PROGRAM_STATE_HH
-#define SRC_PROGRAM_STATE_HH
-
-#include <dune/common/parametertree.hh>
-
-#include <dune/fufem/boundarypatch.hh>
-#include <dune/tnnmg/nonlinearities/zerononlinearity.hh>
-#include <dune/tnnmg/problem-classes/blocknonlineartnnmgproblem.hh>
-
-#include <dune/tectonic/body.hh>
-
-#include "assemblers.hh"
-#include "matrices.hh"
-#include "spatial-solving/solverfactory.hh"
-
-template <class VectorTEMPLATE, class ScalarVectorTEMPLATE> class ProgramState {
-public:
-  using Vector = VectorTEMPLATE;
-  using ScalarVector = ScalarVectorTEMPLATE;
-
-  ProgramState(int leafVertexCount)
-      : u(leafVertexCount),
-        v(leafVertexCount),
-        a(leafVertexCount),
-        alpha(leafVertexCount),
-        weightedNormalStress(leafVertexCount) {}
-
-  // Set up initial conditions
-  template <class Matrix, class GridView>
-  void setupInitialConditions(
-      Dune::ParameterTree const &parset,
-      std::function<void(double, Vector &)> externalForces,
-      Matrices<Matrix> const matrices,
-      MyAssembler<GridView, Vector::block_type::dimension> const &myAssembler,
-      Dune::BitSetVector<Vector::block_type::dimension> const &dirichletNodes,
-      Dune::BitSetVector<Vector::block_type::dimension> const &noNodes,
-      BoundaryPatch<GridView> const &frictionalBoundary,
-      Body<Vector::block_type::dimension> const &body) {
-
-    using LocalVector = typename Vector::block_type;
-    using LocalMatrix = typename Matrix::block_type;
-    auto constexpr dims = LocalVector::dimension;
-
-    // Solving a linear problem with a multigrid solver
-    auto const solveLinearProblem = [&](
-        Dune::BitSetVector<dims> const &_dirichletNodes, Matrix const &_matrix,
-        Vector const &_rhs, Vector &_x,
-        Dune::ParameterTree const &_localParset) {
-
-      using LinearFactory = SolverFactory<
-          dims, BlockNonlinearTNNMGProblem<ConvexProblem<
-                    ZeroNonlinearity<LocalVector, LocalMatrix>, Matrix>>,
-          typename GridView::Grid>;
-      ZeroNonlinearity<LocalVector, LocalMatrix> zeroNonlinearity;
-
-      LinearFactory factory(parset.sub("solver.tnnmg"), // FIXME
-                            myAssembler.gridView.grid(), _dirichletNodes);
-
-      typename LinearFactory::ConvexProblem convexProblem(
-          1.0, _matrix, zeroNonlinearity, _rhs, _x);
-      typename LinearFactory::BlockProblem problem(parset, convexProblem);
-
-      auto multigridStep = factory.getStep();
-      multigridStep->setProblem(_x, problem);
-      EnergyNorm<Matrix, Vector> const norm(_matrix);
-      LoopSolver<Vector> solver(
-          multigridStep.get(), _localParset.get<size_t>("maximumIterations"),
-          _localParset.get<double>("tolerance"), &norm,
-          _localParset.get<Solver::VerbosityMode>("verbosity"),
-          false); // absolute error
-
-      solver.preprocess();
-      solver.solve();
-    };
-
-    timeStep = 0;
-    relativeTime = 0.0;
-    relativeTau = 1e-6;
-
-    Vector ell0(u.size());
-    externalForces(relativeTime, ell0);
-
-    // Initial velocity
-    v = 0.0;
-
-    // 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
-    solveLinearProblem(dirichletNodes, matrices.elasticity, ell0, u,
-                       parset.sub("u0.solver"));
-
-    // Initial acceleration: Computed in agreement with Ma = ell0 - Au
-    // (without Dirichlet constraints), again assuming dPhi(v = 0) = 0
-    Vector accelerationRHS = ell0;
-    Arithmetic::subtractProduct(accelerationRHS, matrices.elasticity, u);
-    solveLinearProblem(noNodes, matrices.mass, accelerationRHS, a,
-                       parset.sub("a0.solver"));
-
-    // Initial state
-    alpha = parset.get<double>("boundary.friction.initialAlpha");
-
-    // Initial normal stress
-    myAssembler.assembleWeightedNormalStress(
-        frictionalBoundary, weightedNormalStress, body.getYoungModulus(),
-        body.getPoissonRatio(), u);
-  }
-
-public:
-  Vector u;
-  Vector v;
-  Vector a;
-  ScalarVector alpha;
-  ScalarVector weightedNormalStress;
-  double relativeTime;
-  double relativeTau;
-  size_t timeStep;
-};
-
-#endif

src/spatial-solving/fixedpointiterator.cc

-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include <dune/common/exceptions.hh>
-
-#include <dune/solvers/common/arithmetic.hh>
-#include <dune/solvers/norms/energynorm.hh>
-#include <dune/solvers/solvers/loopsolver.hh>
-
-#include "../enums.hh"
-#include "../enumparser.hh"
-
-#include "fixedpointiterator.hh"
-
-void FixedPointIterationCounter::operator+=(
-    FixedPointIterationCounter const &other) {
-  iterations += other.iterations;
-  multigridIterations += other.multigridIterations;
-}
-
-template <class Factory, class Updaters, class ErrorNorm>
-FixedPointIterator<Factory, Updaters, ErrorNorm>::FixedPointIterator(
-    Factory &factory, Dune::ParameterTree const &parset,
-    std::shared_ptr<Nonlinearity> globalFriction, ErrorNorm const &errorNorm)
-    : step_(factory.getStep()),
-      parset_(parset),
-      globalFriction_(globalFriction),
-      fixedPointMaxIterations_(parset.get<size_t>("v.fpi.maximumIterations")),
-      fixedPointTolerance_(parset.get<double>("v.fpi.tolerance")),
-      lambda_(parset.get<double>("v.fpi.lambda")),
-      velocityMaxIterations_(parset.get<size_t>("v.solver.maximumIterations")),
-      velocityTolerance_(parset.get<double>("v.solver.tolerance")),
-      verbosity_(parset.get<Solver::VerbosityMode>("v.solver.verbosity")),
-      errorNorm_(errorNorm) {}
-
-template <class Factory, class Updaters, class ErrorNorm>
-FixedPointIterationCounter
-FixedPointIterator<Factory, Updaters, ErrorNorm>::run(
-    Updaters updaters, Matrix const &velocityMatrix, Vector const &velocityRHS,
-    Vector &velocityIterate) {
-  EnergyNorm<Matrix, Vector> energyNorm(velocityMatrix);
-  LoopSolver<Vector> velocityProblemSolver(step_.get(), velocityMaxIterations_,
-                                           velocityTolerance_, &energyNorm,
-                                           verbosity_, false); // absolute error
-
-  size_t fixedPointIteration;
-  size_t multigridIterations = 0;
-  ScalarVector alpha;
-  updaters.state_->extractAlpha(alpha);
-  for (fixedPointIteration = 0; fixedPointIteration < fixedPointMaxIterations_;
-       ++fixedPointIteration) {
-    // solve a velocity problem
-    globalFriction_->updateAlpha(alpha);
-    ConvexProblem convexProblem(1.0, velocityMatrix, *globalFriction_,
-                                velocityRHS, velocityIterate);
-    BlockProblem velocityProblem(parset_, convexProblem);
-    step_->setProblem(velocityIterate, velocityProblem);
-    velocityProblemSolver.preprocess();
-    velocityProblemSolver.solve();
-
-    multigridIterations += velocityProblemSolver.getResult().iterations;
-
-    Vector v_m;
-    updaters.rate_->extractOldVelocity(v_m);
-    v_m *= 1.0 - lambda_;
-    Arithmetic::addProduct(v_m, lambda_, velocityIterate);
-
-    // solve a state problem
-    updaters.state_->solve(v_m);
-    ScalarVector newAlpha;
-    updaters.state_->extractAlpha(newAlpha);
-
-    if (lambda_ < 1e-12 or
-        errorNorm_.diff(alpha, newAlpha) < fixedPointTolerance_) {
-      fixedPointIteration++;
-      break;
-    }
-    alpha = newAlpha;
-  }
-  if (fixedPointIteration == fixedPointMaxIterations_)
-    DUNE_THROW(Dune::Exception, "FPI failed to converge");
-
-  updaters.rate_->postProcess(velocityIterate);
-
-  // Cannot use return { fixedPointIteration, multigridIterations };
-  // with gcc 4.9.2, see also http://stackoverflow.com/a/37777814/179927
-  FixedPointIterationCounter ret;
-  ret.iterations = fixedPointIteration;
-  ret.multigridIterations = multigridIterations;
-  return ret;
-}
-
-std::ostream &operator<<(std::ostream &stream,
-                         FixedPointIterationCounter const &fpic) {
-  return stream << "(" << fpic.iterations << "," << fpic.multigridIterations
-                << ")";
-}
-
-#include "fixedpointiterator_tmpl.cc"

src/spatial-solving/fixedpointiterator_tmpl.cc

-#ifndef MY_DIM
-#error MY_DIM unset
-#endif
-
-#include "../explicitgrid.hh"
-#include "../explicitvectors.hh"
-
-#include <dune/common/function.hh>
-
-#include <dune/solvers/norms/energynorm.hh>
-#include <dune/tnnmg/problem-classes/convexproblem.hh>
-
-#include <dune/tectonic/globalfriction.hh>
-#include <dune/tectonic/myblockproblem.hh>
-
-#include "../time-stepping/rate/rateupdater.hh"
-#include "../time-stepping/state/stateupdater.hh"
-#include "../time-stepping/updaters.hh"
-#include "solverfactory.hh"
-
-using Function = Dune::VirtualFunction<double, double>;
-using Factory = SolverFactory<
-    MY_DIM,
-    MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
-    Grid>;
-using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
-using MyRateUpdater = RateUpdater<Vector, Matrix, Function, MY_DIM>;
-using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
-
-using ErrorNorm = EnergyNorm<ScalarMatrix, ScalarVector>;
-
-template class FixedPointIterator<Factory, MyUpdaters, ErrorNorm>;

src/spatial-solving/solverfactory.cc

-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#ifdef HAVE_IPOPT
-#undef HAVE_IPOPT
-#endif
-
-#include <dune/fufem/assemblers/transferoperatorassembler.hh>
-#include <dune/solvers/solvers/solver.hh>
-
-#include "solverfactory.hh"
-
-template <size_t dim, class BlockProblem, class Grid>
-SolverFactory<dim, BlockProblem, Grid>::SolverFactory(
-    Dune::ParameterTree const &parset, Grid const &grid,
-    Dune::BitSetVector<dim> const &ignoreNodes)
-    : baseEnergyNorm(linearBaseSolverStep),
-      linearBaseSolver(&linearBaseSolverStep,
-                       parset.get<size_t>("linear.maxiumumIterations"),
-                       parset.get<double>("linear.tolerance"), &baseEnergyNorm,
-                       Solver::QUIET),
-      transferOperators(grid.maxLevel()),
-      multigridStep(
-          std::make_shared<Step>(linearIterationStep, nonlinearSmoother)) {
-  // linear iteration step
-  linearIterationStep.setMGType(parset.get<int>("linear.cycle"),
-                                parset.get<int>("linear.pre"),
-                                parset.get<int>("linear.post"));
-  linearIterationStep.basesolver_ = &linearBaseSolver;
-  linearIterationStep.setSmoother(&linearPresmoother, &linearPostsmoother);
-
-  // transfer operators
-  for (auto &&x : transferOperators)
-    x = new CompressedMultigridTransfer<Vector>;
-  TransferOperatorAssembler<Grid>(grid)
-      .assembleOperatorPointerHierarchy(transferOperators);
-  linearIterationStep.setTransferOperators(transferOperators);
-
-  // tnnmg iteration step
-  multigridStep->setSmoothingSteps(parset.get<int>("main.pre"),
-                                   parset.get<int>("main.multi"),
-                                   parset.get<int>("main.post"));
-  multigridStep->ignoreNodes_ = &ignoreNodes;
-}
-
-template <size_t dim, class BlockProblem, class Grid>
-SolverFactory<dim, BlockProblem, Grid>::~SolverFactory() {
-  for (auto &&x : transferOperators)
-    delete x;
-}
-
-template <size_t dim, class BlockProblem, class Grid>
-auto SolverFactory<dim, BlockProblem, Grid>::getStep()
-    -> std::shared_ptr<Step> {
-  return multigridStep;
-}
-
-#include "solverfactory_tmpl.cc"

src/spatial-solving/solverfactory.hh

-#ifndef SRC_SPATIAL_SOLVING_SOLVERFACTORY_HH
-#define SRC_SPATIAL_SOLVING_SOLVERFACTORY_HH
-
-#include <dune/common/bitsetvector.hh>
-#include <dune/common/parametertree.hh>
-
-#include <dune/solvers/iterationsteps/multigridstep.hh>
-#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
-#include <dune/solvers/norms/energynorm.hh>
-#include <dune/solvers/solvers/loopsolver.hh>
-#include <dune/solvers/transferoperators/compressedmultigridtransfer.hh>
-#include <dune/tnnmg/iterationsteps/genericnonlineargs.hh>
-#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
-
-template <size_t dim, class BlockProblemTEMPLATE, class Grid>
-class SolverFactory {
-public:
-  using BlockProblem = BlockProblemTEMPLATE;
-  using ConvexProblem = typename BlockProblem::ConvexProblemType;
-  using Vector = typename BlockProblem::VectorType;
-  using Matrix = typename BlockProblem::MatrixType;
-
-private:
-  using NonlinearSmoother = GenericNonlinearGS<BlockProblem>;
-
-public:
-  using Step =
-      TruncatedNonsmoothNewtonMultigrid<BlockProblem, NonlinearSmoother>;
-
-  SolverFactory(Dune::ParameterTree const &parset, Grid const &grid,
-                Dune::BitSetVector<dim> const &ignoreNodes);
-
-  ~SolverFactory();
-
-  std::shared_ptr<Step> getStep();
-
-private:
-  TruncatedBlockGSStep<Matrix, Vector> linearBaseSolverStep;
-  EnergyNorm<Matrix, Vector> baseEnergyNorm;
-  LoopSolver<Vector> linearBaseSolver;
-  TruncatedBlockGSStep<Matrix, Vector> linearPresmoother;
-  TruncatedBlockGSStep<Matrix, Vector> linearPostsmoother;
-  MultigridStep<Matrix, Vector> linearIterationStep;
-  std::vector<CompressedMultigridTransfer<Vector> *> transferOperators;
-  NonlinearSmoother nonlinearSmoother;
-  std::shared_ptr<Step> multigridStep;
-};
-#endif

src/spatial-solving/solverfactory_tmpl.cc

-#ifndef MY_DIM
-#error MY_DIM unset
-#endif
-
-#include "../explicitgrid.hh"
-#include "../explicitvectors.hh"
-
-#include <dune/tnnmg/nonlinearities/zerononlinearity.hh>
-#include <dune/tnnmg/problem-classes/blocknonlineartnnmgproblem.hh>
-#include <dune/tnnmg/problem-classes/convexproblem.hh>
-
-#include <dune/tectonic/globalfriction.hh>
-#include <dune/tectonic/myblockproblem.hh>
-
-template class SolverFactory<
-    MY_DIM,
-    MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
-    Grid>;
-template class SolverFactory<
-    MY_DIM, BlockNonlinearTNNMGProblem<ConvexProblem<
-                ZeroNonlinearity<LocalVector, LocalMatrix>, Matrix>>,
-    Grid>;

src/strikeslip/CMakeLists.txt

+add_custom_target(tectonic_src_strikeslip SOURCES
+  strikeslip.cfg
+  strikeslip-2D.cfg
+) 
+
+set(STRIKESLIP_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/strikeslipfactory.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/simplexmanager.cc
+  ../../dune/tectonic/problem-data/grid/trianglegeometry.cc
+  ../../dune/tectonic/problem-data/grid/trianglegrids.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
+  strikeslip.cc
+)
+
+foreach(_dim 2 3)
+  set(_strikeslip_target strikeslip-${_dim}D)
+
+  add_executable(${_strikeslip_target} ${STRIKESLIP_SOURCE_FILES})
+
+  add_dune_ug_flags(${_strikeslip_target})
+  add_dune_hdf5_flags(${_strikeslip_target})
+
+  set_property(TARGET ${_strikeslip_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
+  #set_property(TARGET ${_strikeslip_target} APPEND PROPERTY COMPILE_DEFINITIONS "NEW_TNNMG_COMPUTE_ITERATES_DIRECTLY=1")
+endforeach()

src/strikeslip/strikeslip-2D.cfg

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

src/strikeslip/strikeslip.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 <filesystem>
+
+#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/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/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/strikeslipfactory.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/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/time-stepping/adaptivetimestepper.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>
+
+// 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;
+
+const size_t dims = MY_DIM;
+const std::string sourcePath = "/home/joscha/software/dune/dune-tectonic/src/strikeslip/";
+
+Dune::ParameterTree getParameters(int argc, char *argv[]) {
+  Dune::ParameterTree parset;
+  Dune::ParameterTreeParser::readINITree(sourcePath + "strikeslip.cfg", parset);
+  Dune::ParameterTreeParser::readINITree(
+      Dune::Fufem::formatString(sourcePath + "strikeslip-%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[]) {
+  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>("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 + "strikeslip.cfg", outPath, copyOptions);
+    std::filesystem::copy(Dune::Fufem::formatString(sourcePath + "strikeslip-%dD.cfg", dims), outPath, copyOptions);
+    std::filesystem::current_path(outPath);
+
+    std::ofstream out("strikeslip.log");
+    std::streambuf *coutbuf = std::cout.rdbuf(); //save old buffer
+    std::cout.rdbuf(out.rdbuf()); //redirect std::cout to log.txt
+
+    using Assembler = MyAssembler<DefLeafGridView, dims>;
+    using field_type = Matrix::field_type;
+
+    // ----------------------
+    // set up contact network
+    // ----------------------
+    using BlocksFactory = StrikeSlipFactory<Grid, Vector>;
+    BlocksFactory blocksFactory(parset);
+    using ContactNetwork = typename BlocksFactory::ContactNetwork;
+    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();
+    }
+
+    using MyProgramState = ProgramState<Vector, ScalarVector>;
+    MyProgramState programState(nVertices);
+
+    IOHandler<Assembler, ContactNetwork> ioHandler(parset.sub("io"), contactNetwork);
+
+    bool restartRead = ioHandler.read(programState);
+    if (!restartRead) {
+      programState.setupInitialConditions(parset, contactNetwork);
+    }
+
+    auto& nBodyAssembler = contactNetwork.nBodyAssembler();
+    for (size_t i=0; i<bodyCount; i++) {
+      contactNetwork.body(i)->setDeformation(programState.u[i]);
+    }
+    nBodyAssembler.assembleTransferOperator();
+    nBodyAssembler.assembleObstacle();
+
+    // ------------------------
+    // assemble global friction
+    // ------------------------
+    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
+    // -------------------
+
+    BitVector totalDirichletNodes;
+    contactNetwork.totalNodes("dirichlet", totalDirichletNodes);
+
+    /*for (size_t i=0; i<totalDirichletNodes.size(); i++) {
+        bool val = false;
+        for (size_t d=0; d<dims; d++) {
+            val = val || totalDirichletNodes[i][d];
+        }
+
+        totalDirichletNodes[i] = val;
+        for (size_t d=0; d<dims; d++) {
+            totalDirichletNodes[i][d] = val;
+        }
+    }*/
+
+    //print(totalDirichletNodes, "totalDirichletNodes:");
+
+    //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 BoundaryNodes = typename ContactNetwork::BoundaryNodes;
+    using Updaters = Updaters<RateUpdater<Vector, Matrix, BoundaryFunctions, BoundaryNodes>,
+                               StateUpdater<ScalarVector, Vector>>;
+
+    BoundaryFunctions velocityDirichletFunctions;
+    contactNetwork.boundaryFunctions("dirichlet", velocityDirichletFunctions);
+
+    BoundaryNodes dirichletNodes;
+    contactNetwork.boundaryNodes("dirichlet", dirichletNodes);
+
+    /*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++) {
+        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();
+
+    auto const mustRefine = [&](Updaters &coarseUpdater,
+                                Updaters &fineUpdater) {
+
+        //return false;
+      //std::cout << "Step 1" << std::endl;
+
+      std::vector<ScalarVector> coarseAlpha;
+      coarseAlpha.resize(bodyCount);
+      coarseUpdater.state_->extractAlpha(coarseAlpha);
+
+      //print(coarseAlpha, "coarseAlpha:");
+
+      std::vector<ScalarVector> fineAlpha;
+      fineAlpha.resize(bodyCount);
+      fineUpdater.state_->extractAlpha(fineAlpha);
+
+      //print(fineAlpha, "fineAlpha:");
+
+      //std::cout << "Step 3" << std::endl;
+
+      ScalarVector::field_type energyNorm = 0;
+      for (size_t i=0; i<stateEnergyNorms.size(); i++) {
+          //std::cout << "for " << i << std::endl;
+
+          //std::cout << not stateEnergyNorms[i] << std::endl;
+
+          if (coarseAlpha[i].size()==0 || fineAlpha[i].size()==0)
+              continue;
+
+          energyNorm += stateEnergyNorms[i]->diff(fineAlpha[i], coarseAlpha[i]);
+      }
+      //std::cout << "energy norm: " << energyNorm << " tol: " << refinementTolerance <<  std::endl;
+      //std::cout << "must refine: " << (energyNorm > refinementTolerance) <<  std::endl;
+      return energyNorm > refinementTolerance;
+    };
+
+
+    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);
+
+    StepBase<NonlinearFactory, std::decay_t<decltype(contactNetwork)>, Updaters, std::decay_t<decltype(stateEnergyNorms)>>
+        stepBase(parset, contactNetwork, totalDirichletNodes, globalFriction, frictionNodes,
+                 externalForces, stateEnergyNorms);
+
+    AdaptiveTimeStepper<NonlinearFactory, std::decay_t<decltype(contactNetwork)>, Updaters, std::decay_t<decltype(stateEnergyNorms)>>
+        adaptiveTimeStepper(stepBase, contactNetwork, current,
+                            programState.relativeTime, programState.relativeTau,
+                            mustRefine);
+
+    size_t timeSteps = parset.get<size_t>("timeSteps.timeSteps");
+
+    while (!adaptiveTimeStepper.reachedEnd()) {
+      programState.timeStep++;
+
+      //preconditioner.build();
+      iterationCount = adaptiveTimeStepper.advance();
+
+      programState.relativeTime = adaptiveTimeStepper.relativeTime_;
+      programState.relativeTau = adaptiveTimeStepper.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:");*/
+
+      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;
+  }
+}

src/strikeslip/strikeslip.cfg

+outPath = test  # output written to ./output/strikeslip/outPath
+
+# -*- mode:conf -*-
+gravity         = 0.0     # [m/s^2]
+
+[body0]
+length          = 0.5      # [m]
+height          = 0.5     # [m]
+depth           = 0.12     # [m]
+bulkModulus     = 0.5e5    # [Pa] #2190
+poissonRatio    = 0.3     # [1]  #0.11
+[body0.elastic]
+distFromDiag    = 0.2
+density         = 900      # [kg/m^3] #750
+shearViscosity  = 1e3     # [Pas]
+bulkViscosity   = 1e3     # [Pas]
+[body0.viscoelastic]
+density         = 1000     # [kg/m^3]
+shearViscosity  = 1e4     # [Pas]
+bulkViscosity   = 1e4     # [Pas]
+
+[body1]
+length          = 0.5     # [m]
+height          = 0.5     # [m]
+depth           = 0.12     # [m]
+bulkModulus     = 1e5    # [Pa]
+poissonRatio    = 0.0     # [1]
+[body1.elastic]
+distFromDiag    = 0.2
+density         = 900      # [kg/m^3]
+shearViscosity  = 0.0    # [Pas]
+bulkViscosity   = 0.0     # [Pas]
+[body1.viscoelastic]
+density         = 1000     # [kg/m^3]
+shearViscosity  = 0.0     # [Pas]
+bulkViscosity   = 0.0     # [Pas]
+
+[boundary.friction]
+C               = 10       # [Pa]
+mu0             = 0.7      # [ ]
+V0              = 5e-5     # [m/s]
+L               = 2.25e-5  # [m]
+initialAlpha    = 0        # [ ]
+stateModel      = AgeingLaw
+frictionModel   = Truncated #Regularised
+[boundary.friction.weakening]
+a               = 0.002    # [ ]
+b               = 0.017    # [ ]
+[boundary.friction.strengthening]
+a               = 0.020    # [ ]
+b               = 0.005    # [ ]
+
+[boundary.neumann]
+sigmaN          = 200.0      # [Pa]
+
+[boundary.dirichlet]
+finalVelocity   = 1e-4     # [m/s]
+
+[io]
+data.write      = true
+printProgress   = true
+restarts.first  = 0
+restarts.spacing= 50
+restarts.write  = true #true
+vtk.write       = true
+
+[problem]
+finalTime       = 100     # [s] #1000
+bodyCount       = 2
+
+[initialTime]
+timeStep = 0
+relativeTime = 0.0
+relativeTau = 2e-2 # 1e-6
+
+[timeSteps]
+scheme = newmark
+timeSteps = 10
+
+[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
+
+
+

src/time-stepping/adaptivetimestepper.cc

-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include "adaptivetimestepper.hh"
-
-void IterationRegister::registerCount(FixedPointIterationCounter count) {
-  totalCount += count;
-}
-
-void IterationRegister::registerFinalCount(FixedPointIterationCounter count) {
-  finalCount = count;
-}
-
-void IterationRegister::reset() {
-  totalCount = FixedPointIterationCounter();
-  finalCount = FixedPointIterationCounter();
-}
-
-template <class Factory, class Updaters, class ErrorNorm>
-AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::AdaptiveTimeStepper(
-    Factory &factory, Dune::ParameterTree const &parset,
-    std::shared_ptr<Nonlinearity> globalFriction, Updaters &current,
-    double relativeTime, double relativeTau,
-    std::function<void(double, Vector &)> externalForces,
-    ErrorNorm const &errorNorm,
-    std::function<bool(Updaters &, Updaters &)> mustRefine)
-    : relativeTime_(relativeTime),
-      relativeTau_(relativeTau),
-      finalTime_(parset.get<double>("problem.finalTime")),
-      factory_(factory),
-      parset_(parset),
-      globalFriction_(globalFriction),
-      current_(current),
-      R1_(),
-      externalForces_(externalForces),
-      mustRefine_(mustRefine),
-      errorNorm_(errorNorm) {}
-
-template <class Factory, class Updaters, class ErrorNorm>
-bool AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::reachedEnd() {
-  return relativeTime_ >= 1.0;
-}
-
-template <class Factory, class Updaters, class ErrorNorm>
-IterationRegister AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::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
-    |F1|F2|     | the result of (R1+R2) with C, i.e. two steps of size
-                  tau with one of size 2*tau. To check if we need to
-    refine, we compare the result of (F1+F2) with R1, i.e. two steps
-    of size tau/2 with one of size tau. The method makes multiple
-    coarsening/refining attempts, with coarsening coming first. */
-  if (R1_.updaters == Updaters())
-    R1_ = step(current_, relativeTime_, relativeTau_);
-
-  bool didCoarsen = false;
-  iterationRegister_.reset();
-  UpdatersWithCount R2;
-  UpdatersWithCount C;
-  while (relativeTime_ + relativeTau_ <= 1.0) {
-    R2 = step(R1_.updaters, relativeTime_ + relativeTau_, relativeTau_);
-    C = step(current_, relativeTime_, 2 * relativeTau_);
-    if (mustRefine_(R2.updaters, C.updaters))
-      break;
-
-    didCoarsen = true;
-    relativeTau_ *= 2;
-    R1_ = C;
-  }
-  UpdatersWithCount F1;
-  UpdatersWithCount F2;
-  if (!didCoarsen) {
-    while (true) {
-      F1 = step(current_, relativeTime_, relativeTau_ / 2.0);
-      F2 = step(F1.updaters, relativeTime_ + relativeTau_ / 2.0,
-                relativeTau_ / 2.0);
-      if (!mustRefine_(F2.updaters, R1_.updaters))
-        break;
-
-      relativeTau_ /= 2.0;
-      R1_ = F1;
-      R2 = F2;
-    }
-  }
-
-  iterationRegister_.registerFinalCount(R1_.count);
-  relativeTime_ += relativeTau_;
-  current_ = R1_.updaters;
-  R1_ = R2;
-
-  return iterationRegister_;
-}
-
-template <class Factory, class Updaters, class ErrorNorm>
-typename AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::UpdatersWithCount
-AdaptiveTimeStepper<Factory, Updaters, ErrorNorm>::step(
-    Updaters const &oldUpdaters, double rTime, double rTau) {
-  UpdatersWithCount newUpdatersAndCount = {oldUpdaters.clone(), {}};
-  newUpdatersAndCount.count =
-      MyCoupledTimeStepper(finalTime_, factory_, parset_, globalFriction_,
-                           newUpdatersAndCount.updaters, errorNorm_,
-                           externalForces_)
-          .step(rTime, rTau);
-  iterationRegister_.registerCount(newUpdatersAndCount.count);
-  return newUpdatersAndCount;
-}
-
-#include "adaptivetimestepper_tmpl.cc"

src/time-stepping/adaptivetimestepper_tmpl.cc

-#ifndef MY_DIM
-#error MY_DIM unset
-#endif
-
-#include "../explicitgrid.hh"
-#include "../explicitvectors.hh"
-
-#include <dune/common/function.hh>
-
-#include <dune/solvers/norms/energynorm.hh>
-#include <dune/tnnmg/problem-classes/convexproblem.hh>
-
-#include <dune/tectonic/globalfriction.hh>
-#include <dune/tectonic/myblockproblem.hh>
-
-#include "../spatial-solving/solverfactory.hh"
-#include "rate/rateupdater.hh"
-#include "state/stateupdater.hh"
-#include "updaters.hh"
-
-using Function = Dune::VirtualFunction<double, double>;
-using Factory = SolverFactory<
-    MY_DIM,
-    MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
-    Grid>;
-using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
-using MyRateUpdater = RateUpdater<Vector, Matrix, Function, MY_DIM>;
-using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
-
-using ErrorNorm = EnergyNorm<ScalarMatrix, ScalarVector>;
-
-template class AdaptiveTimeStepper<Factory, MyUpdaters, ErrorNorm>;

src/time-stepping/coupledtimestepper.cc

-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include "coupledtimestepper.hh"
-
-template <class Factory, class Updaters, class ErrorNorm>
-CoupledTimeStepper<Factory, Updaters, ErrorNorm>::CoupledTimeStepper(
-    double finalTime, Factory &factory, Dune::ParameterTree const &parset,
-    std::shared_ptr<Nonlinearity> globalFriction, Updaters updaters,
-    ErrorNorm const &errorNorm,
-    std::function<void(double, Vector &)> externalForces)
-    : finalTime_(finalTime),
-      factory_(factory),
-      parset_(parset),
-      globalFriction_(globalFriction),
-      updaters_(updaters),
-      externalForces_(externalForces),
-      errorNorm_(errorNorm) {}
-
-template <class Factory, class Updaters, class ErrorNorm>
-FixedPointIterationCounter
-CoupledTimeStepper<Factory, Updaters, ErrorNorm>::step(double relativeTime,
-                                                       double relativeTau) {
-  updaters_.state_->nextTimeStep();
-  updaters_.rate_->nextTimeStep();
-
-  auto const newRelativeTime = relativeTime + relativeTau;
-  Vector ell;
-  externalForces_(newRelativeTime, ell);
-
-  Matrix velocityMatrix;
-  Vector velocityRHS;
-  Vector velocityIterate;
-
-  auto const tau = relativeTau * finalTime_;
-  updaters_.state_->setup(tau);
-  updaters_.rate_->setup(ell, tau, newRelativeTime, velocityRHS,
-                         velocityIterate, velocityMatrix);
-  FixedPointIterator<Factory, Updaters, ErrorNorm> fixedPointIterator(
-      factory_, parset_, globalFriction_, errorNorm_);
-  auto const iterations = fixedPointIterator.run(updaters_, velocityMatrix,
-                                                 velocityRHS, velocityIterate);
-  return iterations;
-}
-
-#include "coupledtimestepper_tmpl.cc"

src/time-stepping/coupledtimestepper.hh

-#ifndef SRC_TIME_STEPPING_COUPLEDTIMESTEPPER_HH
-#define SRC_TIME_STEPPING_COUPLEDTIMESTEPPER_HH
-
-#include <functional>
-#include <memory>
-
-#include <dune/common/parametertree.hh>
-
-#include "../spatial-solving/fixedpointiterator.hh"
-
-template <class Factory, class Updaters, class ErrorNorm>
-class CoupledTimeStepper {
-  using Vector = typename Factory::Vector;
-  using Matrix = typename Factory::Matrix;
-  using ConvexProblem = typename Factory::ConvexProblem;
-  using Nonlinearity = typename ConvexProblem::NonlinearityType;
-
-public:
-  CoupledTimeStepper(double finalTime, Factory &factory,
-                     Dune::ParameterTree const &parset,
-                     std::shared_ptr<Nonlinearity> globalFriction,
-                     Updaters updaters, ErrorNorm const &errorNorm,
-                     std::function<void(double, Vector &)> externalForces);
-
-  FixedPointIterationCounter step(double relativeTime, double relativeTau);
-
-private:
-  double finalTime_;
-  Factory &factory_;
-  Dune::ParameterTree const &parset_;
-  std::shared_ptr<Nonlinearity> globalFriction_;
-  Updaters updaters_;
-  std::function<void(double, Vector &)> externalForces_;
-  ErrorNorm const &errorNorm_;
-};
-#endif

src/time-stepping/coupledtimestepper_tmpl.cc

-#ifndef MY_DIM
-#error MY_DIM unset
-#endif
-
-#include "../explicitgrid.hh"
-#include "../explicitvectors.hh"
-
-#include <dune/common/function.hh>
-
-#include <dune/solvers/norms/energynorm.hh>
-#include <dune/tnnmg/problem-classes/convexproblem.hh>
-
-#include <dune/tectonic/globalfriction.hh>
-#include <dune/tectonic/myblockproblem.hh>
-
-#include "../spatial-solving/solverfactory.hh"
-#include "rate/rateupdater.hh"
-#include "state/stateupdater.hh"
-#include "updaters.hh"
-
-using Function = Dune::VirtualFunction<double, double>;
-using Factory = SolverFactory<
-    MY_DIM,
-    MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
-    Grid>;
-using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
-using MyRateUpdater = RateUpdater<Vector, Matrix, Function, MY_DIM>;
-using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
-
-using ErrorNorm = EnergyNorm<ScalarMatrix, ScalarVector>;
-
-template class CoupledTimeStepper<Factory, MyUpdaters, ErrorNorm>;

src/time-stepping/rate.hh

-#ifndef SRC_TIME_STEPPING_RATE_HH
-#define SRC_TIME_STEPPING_RATE_HH
-
-#include <memory>
-
-#include "../enums.hh"
-#include "rate/rateupdater.hh"
-
-template <class Vector, class Matrix, class Function, int dimension>
-std::shared_ptr<RateUpdater<Vector, Matrix, Function, dimension>>
-initRateUpdater(Config::scheme scheme,
-                Function const &velocityDirichletFunction,
-                Dune::BitSetVector<dimension> const &velocityDirichletNodes,
-                Matrices<Matrix> const &matrices, Vector const &u_initial,
-                Vector const &v_initial, Vector const &a_initial);
-#endif

src/time-stepping/rate/backward_euler.cc

-#include <dune/solvers/common/arithmetic.hh>
-
-#include "backward_euler.hh"
-
-template <class Vector, class Matrix, class Function, size_t dim>
-BackwardEuler<Vector, Matrix, Function, dim>::BackwardEuler(
-    Matrices<Matrix> const &_matrices, Vector const &_u_initial,
-    Vector const &_v_initial, Vector const &_a_initial,
-    Dune::BitSetVector<dim> const &_dirichletNodes,
-    Function const &_dirichletFunction)
-    : RateUpdater<Vector, Matrix, Function, dim>(
-          _matrices, _u_initial, _v_initial, _a_initial, _dirichletNodes,
-          _dirichletFunction) {}
-
-template <class Vector, class Matrix, class Function, size_t dim>
-void BackwardEuler<Vector, Matrix, Function, dim>::setup(
-    Vector const &ell, double _tau, double relativeTime, Vector &rhs,
-    Vector &iterate, Matrix &AM) {
-  this->dirichletFunction.evaluate(relativeTime, this->dirichletValue);
-  this->tau = _tau;
-
-  /* We start out with the formulation
-
-       M a + C v + A u = ell
-
-     Backward Euler means
-
-       a1 = 1.0/tau ( v1 - v0 )
-       u1 = tau v1 + u0
-
-     in summary, we get at time t=1
-
-       M [1.0/tau ( v1 - v0 )] + C v1
-       + A [tau v1 + u0] = ell
-
-     or
-
-       1.0/tau M v1 + C v1 + tau A v1
-       = [1.0/tau M + C + tau A] v1
-       = ell + 1.0/tau M v0 - A u0
-  */
-
-  // set up LHS (for fixed tau, we'd only really have to do this once)
-  {
-    Dune::MatrixIndexSet indices(this->matrices.elasticity.N(),
-                                 this->matrices.elasticity.M());
-    indices.import(this->matrices.elasticity);
-    indices.import(this->matrices.mass);
-    indices.import(this->matrices.damping);
-    indices.exportIdx(AM);
-  }
-  AM = 0.0;
-  Arithmetic::addProduct(AM, 1.0 / this->tau, this->matrices.mass);
-  Arithmetic::addProduct(AM, 1.0, this->matrices.damping);
-  Arithmetic::addProduct(AM, this->tau, this->matrices.elasticity);
-
-  // set up RHS
-  {
-    rhs = ell;
-    Arithmetic::addProduct(rhs, 1.0 / this->tau, this->matrices.mass,
-                           this->v_o);
-    Arithmetic::subtractProduct(rhs, this->matrices.elasticity, this->u_o);
-  }
-
-  iterate = this->v_o;
-
-  for (size_t i = 0; i < this->dirichletNodes.size(); ++i)
-    for (size_t j = 0; j < dim; ++j)
-      if (this->dirichletNodes[i][j])
-        iterate[i][j] = (j == 0) ? this->dirichletValue : 0;
-}
-
-template <class Vector, class Matrix, class Function, size_t dim>
-void BackwardEuler<Vector, Matrix, Function, dim>::postProcess(
-    Vector const &iterate) {
-  this->postProcessCalled = true;
-
-  this->v = iterate;
-
-  this->u = this->u_o;
-  Arithmetic::addProduct(this->u, this->tau, this->v);
-
-  this->a = this->v;
-  this->a -= this->v_o;
-  this->a /= this->tau;
-}
-
-template <class Vector, class Matrix, class Function, size_t dim>
-std::shared_ptr<RateUpdater<Vector, Matrix, Function, dim>>
-BackwardEuler<Vector, Matrix, Function, dim>::clone() const {
-  return std::make_shared<BackwardEuler<Vector, Matrix, Function, dim>>(*this);
-}