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]
+refinements = 4  # 2e-3 [m]
+initialXElements = 5
+initialYElements = 1
+
+[body1]
+refinements = 4 # 2e-3 [m]
+initialXElements = 5
+initialYElements = 1
+
+[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();
+
+    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);
+
+    const auto minTau = parset.get<double>("initialTime.minRelativeTau");
+    AdaptiveTimeStepper<NonlinearFactory, std::decay_t<decltype(contactNetwork)>, Updaters, std::decay_t<decltype(stateEnergyNorms)>>
+        timeStepper(stepBase, contactNetwork, current,
+                            programState.relativeTime, programState.relativeTau, minTau,
+                            mustRefine);
+
+    /*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 = pipping-2013-newmark-double-adaptive  # output written to ./output/outPath
+gravity         = 9.81     # [m/s^2]
+
+[body0]
+length          = 5.0      # [m]
+height          = 1.0     # [m]
+bulkModulus     = 4.12e7 #4.12e9   # [Pa] #2190
+poissonRatio    = 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 #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-4     # [m/s]
+
+[io]
+data.write      = true
+printProgress   = true
+restarts.first  = 0
+restarts.spacing= 50
+restarts.write  = true #true
+vtk.write       = false
+
+[problem]
+finalTime       = 50     # [s] #1000
+bodyCount       = 2
+
+[initialTime]
+timeStep = 0
+relativeTime = 0.0
+relativeTau = 1e-3 # 1e-6
+minRelativeTau = 1e-6
+
+[timeSteps]
+scheme = newmark # newmark, backwardEuler
+timeSteps = 1e6
+
+[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/hdf5-writer.hh

-#ifndef SRC_HDF5_WRITER_HH
-#define SRC_HDF5_WRITER_HH
-
-#include <dune/fufem/functions/basisgridfunction.hh>
-#include <dune/fufem/geometry/convexpolyhedron.hh>
-#include <dune/fufem/hdf5/file.hh>
-
-#include "hdf5/frictionalboundary-writer.hh"
-#include "hdf5/iteration-writer.hh"
-#include "hdf5/patchinfo-writer.hh"
-#include "hdf5/surface-writer.hh"
-#include "hdf5/time-writer.hh"
-
-template <class ProgramState, class VertexBasis, class GridView>
-class HDF5Writer {
-private:
-  using Vector = typename ProgramState::Vector;
-  using Patch = BoundaryPatch<GridView>;
-
-  using LocalVector = typename Vector::block_type;
-
-public:
-  HDF5Writer(HDF5::Grouplike &file, Vector const &vertexCoordinates,
-             VertexBasis const &vertexBasis, Patch const &surface,
-             Patch const &frictionalBoundary,
-             ConvexPolyhedron<LocalVector> const &weakPatch)
-      : file_(file),
-        iterationWriter_(file_),
-        timeWriter_(file_),
-#if MY_DIM == 3
-        patchInfoWriter_(file_, vertexBasis, frictionalBoundary, weakPatch),
-#endif
-        surfaceWriter_(file_, vertexCoordinates, surface),
-        frictionalBoundaryWriter_(file_, vertexCoordinates,
-                                  frictionalBoundary) {
-  }
-
-  template <class Friction>
-  void reportSolution(ProgramState const &programState,
-                      // for the friction coefficient
-                      Friction &friction) {
-    timeWriter_.write(programState);
-#if MY_DIM == 3
-    patchInfoWriter_.write(programState);
-#endif
-    surfaceWriter_.write(programState);
-    frictionalBoundaryWriter_.write(programState, friction);
-  }
-
-  void reportIterations(ProgramState const &programState,
-                        IterationRegister const &iterationCount) {
-    iterationWriter_.write(programState.timeStep, iterationCount);
-  }
-
-private:
-  HDF5::Grouplike &file_;
-
-  IterationWriter iterationWriter_;
-  TimeWriter<ProgramState> timeWriter_;
-#if MY_DIM == 3
-  PatchInfoWriter<ProgramState, VertexBasis, GridView> patchInfoWriter_;
-#endif
-  SurfaceWriter<ProgramState, GridView> surfaceWriter_;
-  FrictionalBoundaryWriter<ProgramState, GridView> frictionalBoundaryWriter_;
-};
-#endif

src/hdf5/frictionalboundary-writer.cc

-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include <dune/fufem/hdf5/singletonwriter.hh>
-
-#include "frictionalboundary-writer.hh"
-#include "restrict.hh"
-
-template <class ProgramState, class GridView>
-FrictionalBoundaryWriter<ProgramState, GridView>::FrictionalBoundaryWriter(
-    HDF5::Grouplike &file, Vector const &vertexCoordinates,
-    Patch const &frictionalBoundary)
-    : group_(file, "frictionalBoundary"),
-      frictionalBoundary_(frictionalBoundary),
-      frictionalBoundaryDisplacementWriter_(group_, "displacement",
-                                            frictionalBoundary.numVertices(),
-                                            Vector::block_type::dimension),
-      frictionalBoundaryVelocityWriter_(group_, "velocity",
-                                        frictionalBoundary.numVertices(),
-                                        Vector::block_type::dimension),
-      frictionalBoundaryStateWriter_(group_, "state",
-                                     frictionalBoundary.numVertices()),
-      frictionalBoundaryCoefficientWriter_(group_, "coefficient",
-                                           frictionalBoundary.numVertices()) {
-  auto const frictionalBoundaryCoordinates =
-      restrictToSurface(vertexCoordinates, frictionalBoundary);
-  HDF5::SingletonWriter<2> frictionalBoundaryCoordinateWriter(
-      group_, "coordinates", frictionalBoundaryCoordinates.size(),
-      Vector::block_type::dimension);
-  setEntry(frictionalBoundaryCoordinateWriter, frictionalBoundaryCoordinates);
-}
-
-template <class ProgramState, class GridView>
-template <class Friction>
-void FrictionalBoundaryWriter<ProgramState, GridView>::write(
-    ProgramState const &programState, Friction &friction) {
-  auto const frictionalBoundaryDisplacements =
-      restrictToSurface(programState.u, frictionalBoundary_);
-  addEntry(frictionalBoundaryDisplacementWriter_, programState.timeStep,
-           frictionalBoundaryDisplacements);
-
-  auto const frictionalBoundaryVelocities =
-      restrictToSurface(programState.v, frictionalBoundary_);
-  addEntry(frictionalBoundaryVelocityWriter_, programState.timeStep,
-           frictionalBoundaryVelocities);
-
-  auto const frictionalBoundaryStates =
-      restrictToSurface(programState.alpha, frictionalBoundary_);
-  addEntry(frictionalBoundaryStateWriter_, programState.timeStep,
-           frictionalBoundaryStates);
-
-  friction.updateAlpha(programState.alpha);
-  auto const c = friction.coefficientOfFriction(programState.v);
-  auto const frictionalBoundaryCoefficient =
-      restrictToSurface(c, frictionalBoundary_);
-  addEntry(frictionalBoundaryCoefficientWriter_, programState.timeStep,
-           frictionalBoundaryCoefficient);
-}
-
-#include "frictionalboundary-writer_tmpl.cc"

src/hdf5/frictionalboundary-writer.hh

-#ifndef SRC_HDF_FRICTIONALBOUNDARY_WRITER_HH
-#define SRC_HDF_FRICTIONALBOUNDARY_WRITER_HH
-
-#include <dune/fufem/boundarypatch.hh>
-
-#include <dune/fufem/hdf5/sequenceio.hh>
-
-template <class ProgramState, class GridView> class FrictionalBoundaryWriter {
-  using ScalarVector = typename ProgramState::ScalarVector;
-  using Vector = typename ProgramState::Vector;
-  using Patch = BoundaryPatch<GridView>;
-
-public:
-  FrictionalBoundaryWriter(HDF5::Grouplike &file, Vector const &vertexCoordinates,
-                           Patch const &frictionalBoundary);
-
-  template <class Friction>
-  void write(ProgramState const &programState, Friction &friction);
-
-private:
-  HDF5::Group group_;
-
-  Patch const &frictionalBoundary_;
-
-  HDF5::SequenceIO<2> frictionalBoundaryDisplacementWriter_;
-  HDF5::SequenceIO<2> frictionalBoundaryVelocityWriter_;
-  HDF5::SequenceIO<1> frictionalBoundaryStateWriter_;
-  HDF5::SequenceIO<1> frictionalBoundaryCoefficientWriter_;
-};
-#endif

src/hdf5/frictionalboundary-writer_tmpl.cc

-#include "../explicitvectors.hh"
-#include "../explicitgrid.hh"
-
-#include "../program_state.hh"
-
-using MyProgramState = ProgramState<Vector, ScalarVector>;
-using MyFriction = GlobalFriction<Matrix, Vector>;
-
-template class FrictionalBoundaryWriter<MyProgramState, GridView>;
-template void FrictionalBoundaryWriter<MyProgramState, GridView>::write(
-    MyProgramState const &programState, MyFriction &friction);

src/hdf5/iteration-writer.cc

-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include "iteration-writer.hh"
-
-IterationWriter::IterationWriter(HDF5::Grouplike &file)
-    : group_(file, "iterations"),
-      fpiSubGroup_(group_, "fixedPoint"),
-      mgSubGroup_(group_, "multiGrid"),
-      finalMGIterationWriter_(mgSubGroup_, "final"),
-      finalFPIIterationWriter_(fpiSubGroup_, "final"),
-      totalMGIterationWriter_(mgSubGroup_, "total"),
-      totalFPIIterationWriter_(fpiSubGroup_, "total") {}
-
-void IterationWriter::write(size_t timeStep,
-                            IterationRegister const &iterationCount) {
-  addEntry(finalMGIterationWriter_, timeStep,
-           iterationCount.finalCount.multigridIterations);
-  addEntry(finalFPIIterationWriter_, timeStep,
-           iterationCount.finalCount.iterations);
-  addEntry(totalMGIterationWriter_, timeStep,
-           iterationCount.totalCount.multigridIterations);
-  addEntry(totalFPIIterationWriter_, timeStep,
-           iterationCount.totalCount.iterations);
-}

src/hdf5/iteration-writer.hh

-#ifndef SRC_HDF_ITERATION_WRITER_HH
-#define SRC_HDF_ITERATION_WRITER_HH
-
-#include <dune/fufem/hdf5/sequenceio.hh>
-
-#include "../time-stepping/adaptivetimestepper.hh"
-
-class IterationWriter {
-public:
-  IterationWriter(HDF5::Grouplike &file);
-
-  void write(size_t timeStep, IterationRegister const &iterationCount);
-
-private:
-  HDF5::Group group_;
-  HDF5::Group fpiSubGroup_;
-  HDF5::Group mgSubGroup_;
-
-  HDF5::SequenceIO<0, size_t> finalMGIterationWriter_;
-  HDF5::SequenceIO<0, size_t> finalFPIIterationWriter_;
-  HDF5::SequenceIO<0, size_t> totalMGIterationWriter_;
-  HDF5::SequenceIO<0, size_t> totalFPIIterationWriter_;
-};
-#endif

src/hdf5/patchinfo-writer.cc

-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include <dune/fufem/grid/hierarchic-approximation.hh>
-#include <dune/fufem/hdf5/singletonwriter.hh>
-
-#include "patchinfo-writer.hh"
-
-template <class LocalVector, class GridView>
-GridEvaluator<LocalVector, GridView>::GridEvaluator(
-    ConvexPolyhedron<LocalVector> const &weakPatch, GridView const &gridView) {
-  double const bufferWidth = 0.05 * MyGeometry::lengthScale;
-  auto const xminmax = std::minmax_element(
-      weakPatch.vertices.begin(), weakPatch.vertices.end(),
-      [](LocalVector const &a, LocalVector const &b) { return a[0] < b[0]; });
-  double const xmin = (*xminmax.first)[0] - bufferWidth;
-  double const xspan = (*xminmax.second)[0] + bufferWidth - xmin;
-
-  double const zmin = -MyGeometry::depth / 2.0;
-  double const zspan = MyGeometry::depth;
-
-  double spatialResolution = 0.01 * MyGeometry::lengthScale;
-  size_t const xsteps = std::round(xspan / spatialResolution);
-  size_t const zsteps = std::round(zspan / spatialResolution);
-
-  xCoordinates.resize(xsteps + 1);
-  for (size_t xi = 0; xi <= xsteps; xi++)
-    xCoordinates[xi] = xmin + xi * xspan / xsteps;
-  zCoordinates.resize(zsteps + 1);
-  for (size_t zi = 0; zi <= zsteps; zi++)
-    zCoordinates[zi] = zmin + zi * zspan / zsteps;
-
-  HierarchicApproximation<typename GridView::Grid, GridView> const
-      hApproximation(gridView.grid(), gridView, 1e-6 * MyGeometry::lengthScale);
-
-  LocalVector global(0);
-  localInfo.resize(xsteps + 1);
-  for (size_t xi = 0; xi < xCoordinates.size(); ++xi) {
-    localInfo[xi].resize(zsteps + 1);
-    for (size_t zi = 0; zi < zCoordinates.size(); ++zi) {
-      global[0] = xCoordinates[xi];
-      global[2] = zCoordinates[zi];
-      auto const element = hApproximation.findEntity(global);
-      localInfo[xi][zi] =
-          std::make_pair(element, element.geometry().local(global));
-    }
-  }
-}
-
-template <class LocalVector, class GridView>
-template <class Function>
-Dune::Matrix<typename Function::RangeType>
-GridEvaluator<LocalVector, GridView>::evaluate(Function const &f) const {
-  Dune::Matrix<typename Function::RangeType> ret(xCoordinates.size(),
-                                                 zCoordinates.size());
-  for (size_t xi = 0; xi < localInfo.size(); ++xi) {
-    auto const &localInfoX = localInfo[xi];
-    for (size_t zi = 0; zi < localInfoX.size(); ++zi) {
-      auto const &localInfoXZ = localInfoX[zi];
-      f.evaluateLocal(localInfoXZ.first, localInfoXZ.second, ret[xi][zi]);
-    }
-  }
-  return ret;
-}
-
-template <class ProgramState, class VertexBasis, class GridView>
-PatchInfoWriter<ProgramState, VertexBasis, GridView>::PatchInfoWriter(
-    HDF5::Grouplike &file, VertexBasis const &vertexBasis,
-    Patch const &frictionalBoundary,
-    ConvexPolyhedron<LocalVector> const &weakPatch)
-    : group_(file, "weakPatchGrid"),
-      vertexBasis_(vertexBasis),
-      gridEvaluator_(weakPatch, frictionalBoundary.gridView()),
-      weakPatchGridVelocityWriter_(
-          group_, "velocity", gridEvaluator_.xCoordinates.size(),
-          gridEvaluator_.zCoordinates.size(), Vector::block_type::dimension) {
-  HDF5::SingletonWriter<1> weakPatchGridXCoordinateWriter(
-      group_, "xCoordinates", gridEvaluator_.xCoordinates.size());
-  setEntry(weakPatchGridXCoordinateWriter, gridEvaluator_.xCoordinates);
-
-  HDF5::SingletonWriter<1> weakPatchGridZCoordinateWriter(
-      group_, "zCoordinates", gridEvaluator_.zCoordinates.size());
-  setEntry(weakPatchGridZCoordinateWriter, gridEvaluator_.zCoordinates);
-}
-
-template <class ProgramState, class VertexBasis, class GridView>
-void PatchInfoWriter<ProgramState, VertexBasis, GridView>::write(
-    ProgramState const &programState) {
-  BasisGridFunction<VertexBasis, Vector> velocity(vertexBasis_, programState.v);
-  auto const gridVelocity = gridEvaluator_.evaluate(velocity);
-  addEntry(weakPatchGridVelocityWriter_, programState.timeStep, gridVelocity);
-}
-
-#include "patchinfo-writer_tmpl.cc"

src/hdf5/patchinfo-writer.hh

-#ifndef SRC_HDF_PATCHINFO_WRITER_HH
-#define SRC_HDF_PATCHINFO_WRITER_HH
-
-#include <dune/istl/matrix.hh>
-
-#include <dune/fufem/boundarypatch.hh>
-#include <dune/fufem/functions/basisgridfunction.hh>
-#include <dune/fufem/geometry/convexpolyhedron.hh>
-#include <dune/fufem/hdf5/sequenceio.hh>
-
-#include "../one-body-problem-data/mygeometry.hh"
-
-template <class LocalVector, class GridView> class GridEvaluator {
-  using Element = typename GridView::Grid::template Codim<0>::Entity;
-
-public:
-  GridEvaluator(ConvexPolyhedron<LocalVector> const &weakPatch,
-                GridView const &gridView);
-
-  template <class Function>
-  Dune::Matrix<typename Function::RangeType> evaluate(Function const &f) const;
-
-  Dune::BlockVector<Dune::FieldVector<double, 1>> xCoordinates;
-  Dune::BlockVector<Dune::FieldVector<double, 1>> zCoordinates;
-
-private:
-  std::vector<std::vector<std::pair<Element, LocalVector>>> localInfo;
-};
-
-template <class ProgramState, class VertexBasis, class GridView>
-class PatchInfoWriter {
-  using Vector = typename ProgramState::Vector;
-  using LocalVector = typename Vector::block_type;
-  using Patch = BoundaryPatch<GridView>;
-
-public:
-  PatchInfoWriter(HDF5::Grouplike &file, VertexBasis const &vertexBasis,
-                  Patch const &frictionalBoundary,
-                  ConvexPolyhedron<LocalVector> const &weakPatch);
-
-  void write(ProgramState const &programState);
-
-private:
-  HDF5::Group group_;
-
-  VertexBasis const &vertexBasis_;
-
-  GridEvaluator<LocalVector, GridView> const gridEvaluator_;
-  HDF5::SequenceIO<3> weakPatchGridVelocityWriter_;
-};
-#endif

src/hdf5/patchinfo-writer_tmpl.cc

-#include "../explicitvectors.hh"
-#include "../explicitgrid.hh"
-
-#include "../program_state.hh"
-
-using MyProgramState = ProgramState<Vector, ScalarVector>;
-using P1Basis = P1NodalBasis<GridView, double>;
-using MyFunction = BasisGridFunction<P1Basis, Vector>;
-
-template class GridEvaluator<LocalVector, GridView>;
-template Dune::Matrix<typename MyFunction::RangeType>
-GridEvaluator<LocalVector, GridView>::evaluate(MyFunction const &f) const;
-
-template class PatchInfoWriter<MyProgramState, P1Basis, GridView>;

src/hdf5/restart-io.cc

-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include "restart-io.hh"
-
-template <class ProgramState>
-RestartIO<ProgramState>::RestartIO(HDF5::Grouplike &group, size_t vertexCount)
-    : displacementWriter_(group, "displacement", vertexCount,
-                          Vector::block_type::dimension),
-      velocityWriter_(group, "velocity", vertexCount,
-                      Vector::block_type::dimension),
-      accelerationWriter_(group, "acceleration", vertexCount,
-                          Vector::block_type::dimension),
-      stateWriter_(group, "state", vertexCount),
-      weightedNormalStressWriter_(group, "weightedNormalStress", vertexCount),
-      relativeTimeWriter_(group, "relativeTime"),
-      relativeTimeIncrementWriter_(group, "relativeTimeIncrement") {}
-
-template <class ProgramState>
-void RestartIO<ProgramState>::write(ProgramState const &programState) {
-  addEntry(displacementWriter_, programState.timeStep, programState.u);
-  addEntry(velocityWriter_, programState.timeStep, programState.v);
-  addEntry(accelerationWriter_, programState.timeStep, programState.a);
-  addEntry(stateWriter_, programState.timeStep, programState.alpha);
-  addEntry(weightedNormalStressWriter_, programState.timeStep,
-           programState.weightedNormalStress);
-  addEntry(relativeTimeWriter_, programState.timeStep,
-           programState.relativeTime);
-  addEntry(relativeTimeIncrementWriter_, programState.timeStep,
-           programState.relativeTau);
-}
-
-template <class ProgramState>
-void RestartIO<ProgramState>::read(size_t timeStep,
-                                   ProgramState &programState) {
-  programState.timeStep = timeStep;
-  readEntry(displacementWriter_, timeStep, programState.u);
-  readEntry(velocityWriter_, timeStep, programState.v);
-  readEntry(accelerationWriter_, timeStep, programState.a);
-  readEntry(stateWriter_, timeStep, programState.alpha);
-  readEntry(weightedNormalStressWriter_, timeStep,
-            programState.weightedNormalStress);
-  readEntry(relativeTimeWriter_, timeStep, programState.relativeTime);
-  readEntry(relativeTimeIncrementWriter_, timeStep, programState.relativeTau);
-}
-
-#include "restart-io_tmpl.cc"

src/hdf5/restart-io.hh

-#ifndef SRC_HDF_RESTART_HDF_HH
-#define SRC_HDF_RESTART_HDF_HH
-
-#include <dune/fufem/hdf5/file.hh>
-#include <dune/fufem/hdf5/sequenceio.hh>
-
-template <class ProgramState> class RestartIO {
-  using ScalarVector = typename ProgramState::ScalarVector;
-  using Vector = typename ProgramState::Vector;
-
-public:
-  RestartIO(HDF5::Grouplike &group, size_t vertexCount);
-
-  void write(ProgramState const &programState);
-
-  void read(size_t timeStep, ProgramState &programState);
-
-private:
-  HDF5::SequenceIO<2> displacementWriter_;
-  HDF5::SequenceIO<2> velocityWriter_;
-  HDF5::SequenceIO<2> accelerationWriter_;
-  HDF5::SequenceIO<1> stateWriter_;
-  HDF5::SequenceIO<1> weightedNormalStressWriter_;
-  HDF5::SequenceIO<0> relativeTimeWriter_;
-  HDF5::SequenceIO<0> relativeTimeIncrementWriter_;
-};
-#endif

src/hdf5/restart-io_tmpl.cc

-#include "../explicitvectors.hh"
-
-#include "../program_state.hh"
-
-using MyProgramState = ProgramState<Vector, ScalarVector>;
-
-template class RestartIO<MyProgramState>;

src/hdf5/restrict.hh

-#ifndef SRC_HDF5_RESTRICT_HH
-#define SRC_HDF5_RESTRICT_HH
-
-#include <cassert>
-
-#include <dune/common/bitsetvector.hh>
-
-#include "../tobool.hh"
-
-template <class Vector, class Patch>
-Vector restrictToSurface(Vector const &v1, Patch const &patch) {
-  auto const &vertices = *patch.getVertices();
-  assert(vertices.size() == v1.size());
-
-  Vector ret(vertices.count());
-  auto target = ret.begin();
-  for (size_t i = 0; i < v1.size(); ++i)
-    if (toBool(vertices[i]))
-      *(target++) = v1[i];
-  assert(target == ret.end());
-  return ret;
-}
-#endif

src/hdf5/surface-writer.hh

-#ifndef SRC_HDF_SURFACE_WRITER_HH
-#define SRC_HDF_SURFACE_WRITER_HH
-
-#include <dune/fufem/boundarypatch.hh>
-
-#include <dune/fufem/hdf5/sequenceio.hh>
-#include <dune/fufem/hdf5/singletonwriter.hh>
-
-template <class ProgramState, class GridView> class SurfaceWriter {
-  using Vector = typename ProgramState::Vector;
-  using Patch = BoundaryPatch<GridView>;
-
-public:
-  SurfaceWriter(HDF5::Grouplike &file, Vector const &vertexCoordinates,
-                Patch const &surface);
-
-  void write(ProgramState const &programState);
-
-private:
-  HDF5::Group group_;
-
-  Patch const &surface_;
-
-  HDF5::SequenceIO<2> surfaceDisplacementWriter_;
-  HDF5::SequenceIO<2> surfaceVelocityWriter_;
-};
-#endif

src/hdf5/time-writer.cc

-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include "time-writer.hh"
-
-template <class ProgramState>
-TimeWriter<ProgramState>::TimeWriter(HDF5::Grouplike &file)
-    : file_(file),
-      relativeTimeWriter_(file_, "relativeTime"),
-      relativeTimeIncrementWriter_(file_, "relativeTimeIncrement") {}
-
-template <class ProgramState>
-void TimeWriter<ProgramState>::write(ProgramState const &programState) {
-  addEntry(relativeTimeWriter_, programState.timeStep,
-           programState.relativeTime);
-  addEntry(relativeTimeIncrementWriter_, programState.timeStep,
-           programState.relativeTau);
-}
-
-#include "time-writer_tmpl.cc"

src/hdf5/time-writer.hh

-#ifndef SRC_HDF_TIME_WRITER_HH
-#define SRC_HDF_TIME_WRITER_HH
-
-#include <dune/fufem/hdf5/file.hh>
-#include <dune/fufem/hdf5/sequenceio.hh>
-
-template <class ProgramState> class TimeWriter {
-public:
-  TimeWriter(HDF5::Grouplike &file);
-  void write(ProgramState const &programState);
-
-private:
-  HDF5::Grouplike &file_;
-  HDF5::SequenceIO<0> relativeTimeWriter_;
-  HDF5::SequenceIO<0> relativeTimeIncrementWriter_;
-};
-#endif