src/sand-wedge-data/patchfunction.hh → src/one-body-problem-data/patchfunction.hh

-#ifndef SRC_SAND_WEDGE_DATA_PATCHFUNCTION_HH
-#define SRC_SAND_WEDGE_DATA_PATCHFUNCTION_HH
+#ifndef SRC_ONE_BODY_PROBLEM_DATA_PATCHFUNCTION_HH
+#define SRC_ONE_BODY_PROBLEM_DATA_PATCHFUNCTION_HH
 
 #include <dune/common/function.hh>
 #include <dune/common/fvector.hh>
 #include <dune/common/parametertree.hh>
 
+#include <dune/fufem/geometry/polyhedrondistance.hh>
+
 class PatchFunction
     : public Dune::VirtualFunction<Dune::FieldVector<double, MY_DIM>,
                                    Dune::FieldVector<double, 1>> {
 private:
-  bool static isBetween(double x, double lower, double upper) {
-    return lower <= x and x <= upper;
-  }
+  using Polyhedron = ConvexPolyhedron<Dune::FieldVector<double, MY_DIM>>;
 
-  bool static isClose(double a, double b) {
-    return std::abs(a - b) < 1e-14;
-  };
-
-  bool insideRegion2(Dune::FieldVector<double, MY_DIM> const &z) const {
-    assert(isClose(0, z[1]));
-    return isBetween(z[0], _X[0], _Y[0]);
-  }
-
-  Dune::FieldVector<double, MY_DIM> const &_X;
-  Dune::FieldVector<double, MY_DIM> const &_Y;
-
-  double const _v1;
-  double const _v2;
+  double const v1_;
+  double const v2_;
+  Polyhedron const &segment_;
 
 public:
-  PatchFunction(double v1, double v2)
-      : _X(MyGeometry::X), _Y(MyGeometry::Y), _v1(v1), _v2(v2) {}
+  PatchFunction(double v1, double v2, Polyhedron const &segment)
+      : v1_(v1), v2_(v2), segment_(segment) {}
 
   void evaluate(Dune::FieldVector<double, MY_DIM> const &x,
                 Dune::FieldVector<double, 1> &y) const {
-    y = insideRegion2(x) ? _v2 : _v1;
+    y = distance(x, segment_, 1e-6 * MyGeometry::lengthScale) <= 1e-5 ? v2_
+                                                                      : v1_;
   }
 };
+
 #endif

src/sand-wedge-data/twopiece.hh → src/one-body-problem-data/segmented-function.hh

-#ifndef SRC_SAND_WEDGE_DATA_TWOPIECE_HH
-#define SRC_SAND_WEDGE_DATA_TWOPIECE_HH
+#ifndef SRC_ONE_BODY_PROBLEM_DATA_SEGMENTED_FUNCTION_HH
+#define SRC_ONE_BODY_PROBLEM_DATA_SEGMENTED_FUNCTION_HH
 
 #include <dune/common/function.hh>
 #include <dune/common/fvector.hh>
@@ -7,28 +7,24 @@
 
 #include "mygeometry.hh"
 
-class TwoPieceFunction
+class SegmentedFunction
     : public Dune::VirtualFunction<Dune::FieldVector<double, MY_DIM>,
                                    Dune::FieldVector<double, 1>> {
 private:
   bool liesBelow(Dune::FieldVector<double, MY_DIM> const &x,
                  Dune::FieldVector<double, MY_DIM> const &y,
                  Dune::FieldVector<double, MY_DIM> const &z) const {
-    return (z[0] - x[0]) * (y[1] - x[1]) / (y[0] - x[0]) >= z[1] - x[1];
+    return x[1] + (z[0] - x[0]) * (y[1] - x[1]) / (y[0] - x[0]) >= z[1];
   };
   bool insideRegion2(Dune::FieldVector<double, MY_DIM> const &z) const {
-    return liesBelow(_K, _M, z);
+    return liesBelow(MyGeometry::K, MyGeometry::M, z);
   };
 
-  Dune::FieldVector<double, MY_DIM> const &_K;
-  Dune::FieldVector<double, MY_DIM> const &_M;
-
   double const _v1;
   double const _v2;
 
 public:
-  TwoPieceFunction(double v1, double v2)
-      : _K(MyGeometry::K), _M(MyGeometry::M), _v1(v1), _v2(v2) {}
+  SegmentedFunction(double v1, double v2) : _v1(v1), _v2(v2) {}
 
   void evaluate(Dune::FieldVector<double, MY_DIM> const &x,
                 Dune::FieldVector<double, 1> &y) const {

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/sand-wedge-data/parset.cfg → src/one-body-problem.cfg

@@ -2,68 +2,64 @@
 gravity         = 9.81  # [m/s^2]
 
 [io]
+data.write      = true
 printProgress   = false
-writeVTK        = false
+restarts.first  = 0
+restarts.spacing= 20
+restarts.write  = true
+vtk.write       = false
 
 [problem]
-finalTime       = 1800  # [s]
+finalTime       = 1000  # [s]
 
 [body]
-bulkModulus     = 1e5   # [Pa]
-poissonRatio    = 0.3   # [1]      0.2 - 0.3
+bulkModulus     = 0.5e5 # [Pa]
+poissonRatio    = 0.3   # [1]
 [body.elastic]
 density         = 900   # [kg/m^3]
-shearViscosity  = 1e3   # [Pas]    0
-bulkViscosity   = 1e3   # [Pas]    0
+shearViscosity  = 1e3   # [Pas]
+bulkViscosity   = 1e3   # [Pas]
 [body.viscoelastic]
 density         = 1000  # [kg/m^3]
 shearViscosity  = 1e4   # [Pas]
 bulkViscosity   = 1e4   # [Pas]
 
 [boundary.friction]
-smallestDiameter= 1e-3  # [m]
 C               = 10    # [Pa]
-mu0             = 0.7   # [1]
+mu0             = 0.7   # [ ]
 V0              = 5e-5  # [m/s]
-L               = 2e-5  # [m]      ?
-initialAlpha    = 0.916290731874155 # [ ]      ?
+L               = 2.25e-5 # [m]
+initialAlpha    = 0     # [ ]
 stateModel      = AgeingLaw
-frictionModel   = Truncated
+frictionModel   = Regularised
 [boundary.friction.weakening]
-a               = 0.015 # [1]      ?
-b               = 0.030 # [1]      ?
+a               = 0.002 # [ ]
+b               = 0.017 # [ ]
 [boundary.friction.strengthening]
-a               = 0.030 # [1]      ?
-b               = 0.015 # [1]      ?
+a               = 0.020 # [ ]
+b               = 0.005 # [ ]
 
 [timeSteps]
-refinementTolerance = 1e-5
-number = 100000
 scheme = newmark
 
 [u0.solver]
-tolerance         = 1e-10
 maximumIterations = 100000
 verbosity         = quiet
 
 [a0.solver]
-tolerance         = 1e-10
 maximumIterations = 100000
 verbosity         = quiet
 
 [v.solver]
-tolerance         = 1e-10
 maximumIterations = 100000
 verbosity         = quiet
 
 [v.fpi]
-tolerance         = 1e-10
 maximumIterations = 10000
 lambda            = 0.5
 
 [solver.tnnmg.linear]
 maxiumumIterations = 100000
-tolerance          = 1e-10
 pre                = 3
 cycle              = 1  # 1 = V, 2 = W, etc.
 post               = 3

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/sand-wedge-data/boundaryconditions.py

-import math
-
-class neumannCondition:
-    def __call__(self, relativeTime):
-        return 0
-
-class velocityDirichletCondition:
-    def __call__(self, relativeTime):
-        finalVelocity = -5e-5
-        if relativeTime <= 0.1:
-            return finalVelocity * ( 1.0 - math.cos(relativeTime * math.pi / 0.1) ) / 2.0
-        return finalVelocity
-
-Functions = {
-    'neumannCondition' : neumannCondition(),
-    'velocityDirichletCondition' : velocityDirichletCondition()
-}

src/sand-wedge-data/special_writer.hh

-#ifndef SRC_SPECIAL_WRITER_HH
-#define SRC_SPECIAL_WRITER_HH
-
-#include <fstream>
-#include <utility>
-
-#include <dune/common/fvector.hh>
-#include <dune/grid/utility/hierarchicsearch.hh>
-#include <dune/fufem/functions/virtualgridfunction.hh>
-
-#include "mygeometry.hh"
-
-template <class GridView, int dimension> class SpecialWriter {
-  using LocalVector = Dune::FieldVector<double, dimension>;
-  using Element = typename GridView::Grid::template Codim<0>::Entity;
-  using ElementPointer =
-      typename GridView::Grid::template Codim<0>::EntityPointer;
-
-  void writeHorizontal(LocalVector const &v) {
-    writer_ << MyGeometry::horizontalProjection(v) << " ";
-  }
-  void writeVertical(LocalVector const &v) {
-    writer_ << MyGeometry::verticalProjection(v) << " ";
-  }
-
-  std::pair<ElementPointer, LocalVector> globalToLocal(LocalVector const &x)
-      const {
-    auto const element = hsearch_.findEntity(x);
-    return std::make_pair(element, element->geometry().local(x));
-  }
-
-  std::fstream writer_;
-
-  Dune::HierarchicSearch<typename GridView::Grid, GridView> const hsearch_;
-
-  std::pair<ElementPointer, LocalVector> const G;
-  std::pair<ElementPointer, LocalVector> const H;
-  std::pair<ElementPointer, LocalVector> const J;
-  std::pair<ElementPointer, LocalVector> const I;
-  std::pair<ElementPointer, LocalVector> const U;
-  std::pair<ElementPointer, LocalVector> const Z;
-
-public:
-  SpecialWriter(std::string filename, GridView const &gridView)
-      : writer_(filename, std::fstream::out),
-        hsearch_(gridView.grid(), gridView),
-        G(globalToLocal(MyGeometry::G)),
-        H(globalToLocal(MyGeometry::H)),
-        J(globalToLocal(MyGeometry::J)),
-        I(globalToLocal(MyGeometry::I)),
-        U(globalToLocal(MyGeometry::U)),
-        Z(globalToLocal(MyGeometry::Z)) {
-    writer_ << "Gh Hh Jh Ih Uv Uh Zv Zh" << std::endl;
-  }
-
-  void write(VirtualGridFunction<typename GridView::Grid, LocalVector> const &
-                 specialField) {
-    LocalVector value;
-
-    specialField.evaluateLocal(*G.first, G.second, value);
-    writeHorizontal(value);
-
-    specialField.evaluateLocal(*H.first, H.second, value);
-    writeHorizontal(value);
-
-    specialField.evaluateLocal(*J.first, J.second, value);
-    writeHorizontal(value);
-
-    specialField.evaluateLocal(*I.first, I.second, value);
-    writeHorizontal(value);
-
-    specialField.evaluateLocal(*U.first, U.second, value);
-    writeVertical(value);
-    writeHorizontal(value);
-
-    specialField.evaluateLocal(*Z.first, Z.second, value);
-    writeVertical(value);
-    writeHorizontal(value);
-
-    writer_ << std::endl;
-  }
-};
-#endif

src/sand-wedge.cc

-#include <Python.h>
-
-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#ifndef HAVE_PYTHON
-#error Python is required
-#endif
-
-#ifdef HAVE_IPOPT
-#undef HAVE_IPOPT
-#endif
-
-#ifndef datadir
-#error datadir unset
-#endif
-
-#include <cmath>
-#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/parametertree.hh>
-#include <dune/common/parametertreeparser.hh>
-
-#include <dune/grid/common/mcmgmapper.hh>
-#include <dune/istl/bcrsmatrix.hh>
-#include <dune/istl/bvector.hh>
-
-#pragma clang diagnostic push
-#pragma clang diagnostic ignored "-Wdeprecated-declarations"
-#include <dune/fufem/boundarypatch.hh>
-#pragma clang diagnostic pop
-
-#include <dune/fufem/dunepython.hh>
-#include <dune/fufem/functions/basisgridfunction.hh>
-#include <dune/fufem/sharedpointermap.hh>
-#include <dune/solvers/norms/energynorm.hh>
-#include <dune/solvers/norms/sumnorm.hh>
-#include <dune/solvers/norms/twonorm.hh>
-#include <dune/solvers/solvers/loopsolver.hh>
-#include <dune/solvers/solvers/solver.hh>
-#include <dune/tnnmg/nonlinearities/zerononlinearity.hh>
-#include <dune/tnnmg/problem-classes/blocknonlineartnnmgproblem.hh>
-#include <dune/tnnmg/problem-classes/convexproblem.hh>
-
-#include <dune/tectonic/myblockproblem.hh>
-#include <dune/tectonic/globalfriction.hh>
-
-#include "adaptivetimestepper.hh" // FIXME
-#include "assemblers.hh"
-#include "tobool.hh"
-#include "coupledtimestepper.hh"
-#include "enumparser.hh"
-#include "enums.hh"
-#include "fixedpointiterator.hh"
-#include "friction_writer.hh"
-#include "gridselector.hh"
-#include "sand-wedge-data/mybody.hh"
-#include "sand-wedge-data/mygeometry.hh"
-#include "sand-wedge-data/myglobalfrictiondata.hh"
-#include "sand-wedge-data/mygrid.hh"
-#include "sand-wedge-data/special_writer.hh"
-#include "solverfactory.hh"
-#include "state.hh"
-#include "timestepping.hh"
-#include "vtk.hh"
-
-size_t const dims = MY_DIM;
-
-void initPython() {
-  Python::start();
-
-  Python::run("import sys");
-  Python::run("sys.path.append('" datadir "')");
-}
-
-template <class Geometry> double diameter(Geometry const &geometry) {
-  auto const numCorners = geometry.corners();
-  std::vector<typename Geometry::GlobalCoordinate> corners(numCorners);
-  for (int i = 0; i < numCorners; ++i)
-    corners[i] = geometry.corner(i);
-
-  TwoNorm<typename Geometry::GlobalCoordinate> twoNorm;
-
-  double diameter = 0.0;
-  for (int i = 0; i < numCorners; ++i)
-    for (int j = 0; j < i; ++j)
-      diameter = std::max(diameter, twoNorm.diff(corners[i], corners[j]));
-  return diameter;
-}
-
-#include <dune/tectonic/closestpointprojection.hh>
-
-// we assume the weakening region to be a line segment for now
-template <class Geometry> double distanceToWeakeningRegion(Geometry const &g) {
-  TwoNorm<typename Geometry::GlobalCoordinate> twoNorm;
-
-  BoundarySegment<typename Geometry::GlobalCoordinate> bs;
-  bs.nVertices = 2;
-  bs.vertices.resize(bs.nVertices);
-  bs.vertices[0] = MyGeometry::X;
-  bs.vertices[1] = MyGeometry::Y;
-
-  double distance = std::numeric_limits<double>::infinity();
-  auto const numCorners = g.corners();
-  for (int i = 0; i < numCorners; ++i) {
-    auto const corner = g.corner(i);
-    auto const projection = computeClosestPoint(corner, bs, 10); // FIXME
-    distance = std::min(distance, twoNorm.diff(corner, projection));
-  }
-  return distance;
-}
-
-double computeAdmissibleDiameter(double distance, double smallestDiameter) {
-  return (distance / 0.0125 + 1.0) * smallestDiameter;
-}
-
-int main(int argc, char *argv[]) {
-  try {
-    Dune::ParameterTree parset;
-    Dune::ParameterTreeParser::readINITree(datadir "/parset.cfg", parset);
-    Dune::ParameterTreeParser::readOptions(argc, argv, parset);
-
-    MyGeometry::render();
-    MyGeometry::write();
-
-    // {{{ Set up grid
-    GridConstructor<Grid> gridConstructor;
-    auto grid = gridConstructor.getGrid();
-
-    auto const smallestDiameter =
-        parset.get<double>("boundary.friction.smallestDiameter");
-    bool needRefine = true;
-
-    while (true) {
-      needRefine = false;
-      for (auto it = grid->template leafbegin<0>();
-           it != grid->template leafend<0>(); ++it) {
-        auto const geometry = it->geometry();
-
-        auto const weakeningRegionDistance =
-            distanceToWeakeningRegion(geometry);
-        auto const admissibleDiameter = computeAdmissibleDiameter(
-            weakeningRegionDistance, smallestDiameter);
-
-        if (diameter(geometry) <= admissibleDiameter)
-          continue;
-
-        needRefine = true;
-        grid->mark(1, *it);
-      }
-      if (!needRefine)
-        break;
-
-      grid->preAdapt();
-      grid->adapt();
-      grid->postAdapt();
-    }
-    auto const refinements = grid->maxLevel();
-
-    double minDiameter = std::numeric_limits<double>::infinity();
-    double maxDiameter = 0.0;
-    for (auto it = grid->template leafbegin<0>();
-         it != grid->template leafend<0>(); ++it) {
-      auto const geometry = it->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;
-
-    using GridView = Grid::LeafGridView;
-    GridView const leafView = grid->leafGridView();
-    size_t const leafVertexCount = leafView.size(dims);
-
-    std::cout << "Number of DOFs: " << leafVertexCount << std::endl;
-
-    auto myFaces = gridConstructor.constructFaces(leafView);
-
-    // Neumann boundary
-    BoundaryPatch<GridView> const neumannBoundary(leafView);
-
-    // Frictional Boundary
-    BoundaryPatch<GridView> const &frictionalBoundary = myFaces.lower;
-    Dune::BitSetVector<1> frictionalNodes(leafVertexCount);
-    frictionalBoundary.getVertices(frictionalNodes);
-
-    // Surface
-    BoundaryPatch<GridView> const &surface = myFaces.upper;
-    Dune::BitSetVector<1> surfaceNodes(leafVertexCount);
-    surface.getVertices(surfaceNodes);
-
-    // 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>;
-    using FunctionMap = SharedPointerMap<std::string, Function>;
-    FunctionMap functions;
-    {
-      initPython();
-      Python::import("boundaryconditions")
-          .get("Functions")
-          .toC<typename FunctionMap::Base>(functions);
-    }
-    auto const &velocityDirichletFunction =
-                   functions.get("velocityDirichletCondition"),
-               &neumannFunction = functions.get("neumannCondition");
-
-    std::fstream timeStepWriter("timeSteps", std::fstream::out);
-    timeStepWriter << std::fixed << std::setprecision(10);
-    auto const reportTimeStep = [&](double _relativeTime, double _relativeTau) {
-      timeStepWriter << _relativeTime << " " << _relativeTau << std::endl;
-    };
-
-    using MyAssembler = MyAssembler<GridView, dims>;
-    using Matrix = MyAssembler::Matrix;
-    using LocalMatrix = Matrix::block_type;
-    using Vector = MyAssembler::Vector;
-    using LocalVector = Vector::block_type;
-    using ScalarMatrix = MyAssembler::ScalarMatrix;
-    using ScalarVector = MyAssembler::ScalarVector;
-    using LocalScalarVector = ScalarVector::block_type;
-
-    MyAssembler myAssembler(leafView);
-
-    MyBody<dims> const body(parset);
-
-    Matrix A, C, M;
-    myAssembler.assembleElasticity(body.getYoungModulus(),
-                                   body.getPoissonRatio(), A);
-    myAssembler.assembleViscosity(body.getShearViscosityField(),
-                                  body.getBulkViscosityField(), C);
-    myAssembler.assembleMass(body.getDensityField(), M);
-    EnergyNorm<Matrix, Vector> const ANorm(A), MNorm(M);
-
-    ScalarMatrix frictionalBoundaryMass;
-    myAssembler.assembleFrictionalBoundaryMass(frictionalBoundary,
-                                               frictionalBoundaryMass);
-    EnergyNorm<ScalarMatrix, ScalarVector> const stateEnergyNorm(
-        frictionalBoundaryMass);
-
-    // 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;
-    };
-    Vector ell0(leafVertexCount);
-    computeExternalForces(0.0, ell0);
-
-    // {{{ Initial conditions
-    using LinearFactory = SolverFactory<
-        dims, BlockNonlinearTNNMGProblem<ConvexProblem<
-                  ZeroNonlinearity<LocalVector, LocalMatrix>, Matrix>>,
-        Grid>;
-    ZeroNonlinearity<LocalVector, LocalMatrix> zeroNonlinearity;
-
-    auto const solveLinearProblem = [&](
-        Dune::BitSetVector<dims> const &_dirichletNodes, Matrix const &_matrix,
-        Vector const &_rhs, Vector &_x, EnergyNorm<Matrix, Vector> const &_norm,
-        Dune::ParameterTree const &_localParset) {
-      LinearFactory factory(parset.sub("solver.tnnmg"), // FIXME
-                            refinements, *grid, _dirichletNodes);
-
-      typename LinearFactory::ConvexProblem convexProblem(
-          1.0, _matrix, zeroNonlinearity, _rhs, _x);
-      typename LinearFactory::BlockProblem problem(parset, convexProblem);
-
-      auto multigridStep = factory.getSolver();
-      multigridStep->setProblem(_x, problem);
-      LoopSolver<Vector> solver(
-          multigridStep, _localParset.get<size_t>("maximumIterations"),
-          _localParset.get<double>("tolerance"), &_norm,
-          _localParset.get<Solver::VerbosityMode>("verbosity"),
-          false); // absolute error
-
-      solver.preprocess();
-      solver.solve();
-    };
-
-    // Solve the stationary problem
-    Vector u_initial(leafVertexCount);
-    u_initial = 0.0;
-    solveLinearProblem(dirichletNodes, A, ell0, u_initial, ANorm,
-                       parset.sub("u0.solver"));
-    Vector ur_initial(leafVertexCount);
-    ur_initial = 0.0;
-
-    ScalarVector alpha_initial(leafVertexCount);
-    alpha_initial = parset.get<double>("boundary.friction.initialAlpha");
-    ScalarVector normalStress(leafVertexCount);
-    myAssembler.assembleNormalStress(frictionalBoundary, normalStress,
-                                     body.getYoungModulus(),
-                                     body.getPoissonRatio(), u_initial);
-    MyGlobalFrictionData<dims> frictionInfo(parset.sub("boundary.friction"));
-    auto myGlobalFriction = myAssembler.assembleFrictionNonlinearity(
-        parset.get<Config::FrictionModel>("boundary.friction.frictionModel"),
-        frictionalBoundary, frictionInfo, normalStress);
-    myGlobalFriction->updateAlpha(alpha_initial);
-
-    Vector v_initial(leafVertexCount);
-    v_initial = 0.0;
-    {
-      double v_initial_const;
-      velocityDirichletFunction.evaluate(0.0, v_initial_const);
-      assert(std::abs(v_initial_const) < 1e-14);
-    }
-    Vector vr_initial(leafVertexCount);
-    vr_initial = 0.0;
-
-    Vector a_initial(leafVertexCount);
-    a_initial = 0.0;
-    {
-      // We solve Ma = ell0 - [Au + Cv + DPsi(v)]
-      Vector accelerationRHS(leafVertexCount);
-      {
-        accelerationRHS = 0.0;
-        Arithmetic::addProduct(accelerationRHS, A, u_initial);
-        Arithmetic::addProduct(accelerationRHS, C, v_initial);
-        // NOTE: We assume differentiability of Psi at 0 here!
-        myGlobalFriction->addGradient(v_initial, accelerationRHS);
-        accelerationRHS *= -1.0;
-        accelerationRHS += ell0;
-      }
-      solveLinearProblem(noNodes, M, accelerationRHS, a_initial, MNorm,
-                         parset.sub("a0.solver"));
-    }
-    // }}}
-
-    Vector vertexCoordinates(leafVertexCount);
-    {
-      Dune::MultipleCodimMultipleGeomTypeMapper<
-          GridView, Dune::MCMGVertexLayout> const vertexMapper(leafView);
-      for (auto it = leafView.begin<dims>(); it != leafView.end<dims>(); ++it) {
-        auto const geometry = it->geometry();
-        assert(geometry.corners() == 1);
-        vertexCoordinates[vertexMapper.map(*it)] = geometry.corner(0);
-      }
-    }
-
-    FrictionWriter<ScalarVector, Vector> frictionWriter(
-        vertexCoordinates, frictionalNodes, "friction",
-        MyGeometry::horizontalProjection);
-    BoundaryWriter<ScalarVector, Vector> verticalSurfaceWriter(
-        vertexCoordinates, surfaceNodes, "verticalSurface",
-        MyGeometry::verticalProjection);
-    BoundaryWriter<ScalarVector, Vector> horizontalSurfaceWriter(
-        vertexCoordinates, surfaceNodes, "horizontalSurface",
-        MyGeometry::horizontalProjection);
-
-    auto const report = [&](Vector const &_u, Vector const &_v,
-                            ScalarVector const &_alpha) {
-      horizontalSurfaceWriter.writeKinetics(_u, _v);
-      verticalSurfaceWriter.writeKinetics(_u, _v);
-
-      ScalarVector c;
-      myGlobalFriction->coefficientOfFriction(_v, c);
-      frictionWriter.writeKinetics(_u, _v);
-      frictionWriter.writeOther(c, _alpha);
-    };
-    report(ur_initial, vr_initial, alpha_initial);
-
-    MyVTKWriter<typename MyAssembler::VertexBasis,
-                typename MyAssembler::CellBasis> const
-    vtkWriter(myAssembler.cellBasis, myAssembler.vertexBasis, "obs");
-
-    if (parset.get<bool>("io.writeVTK")) {
-      ScalarVector stress;
-      myAssembler.assembleVonMisesStress(
-          body.getYoungModulus(), body.getPoissonRatio(), u_initial, stress);
-      vtkWriter.write(0, ur_initial, vr_initial, alpha_initial, stress);
-    }
-
-    SpecialWriter<GridView, dims> specialVelocityWriter("specialVelocities",
-                                                        leafView);
-    SpecialWriter<GridView, dims> specialDisplacementWriter(
-        "specialDisplacements", leafView);
-
-    // Set up TNNMG solver
-    using NonlinearFactory = SolverFactory<
-        dims,
-        MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
-        Grid>;
-    NonlinearFactory factory(parset.sub("solver.tnnmg"), refinements, *grid,
-                             dirichletNodes);
-
-    using UpdaterPair = std::pair<
-        std::shared_ptr<StateUpdater<ScalarVector, Vector>>,
-        std::shared_ptr<TimeSteppingScheme<Vector, Matrix, Function, dims>>>;
-    UpdaterPair current(
-        initStateUpdater<ScalarVector, Vector>(
-            parset.get<Config::stateModel>("boundary.friction.stateModel"),
-            alpha_initial, frictionalNodes,
-            parset.get<double>("boundary.friction.L"),
-            parset.get<double>("boundary.friction.V0")),
-        initTimeStepper(parset.get<Config::scheme>("timeSteps.scheme"),
-                        velocityDirichletFunction, dirichletNodes, M, A, C,
-                        u_initial, ur_initial, v_initial, vr_initial,
-                        a_initial));
-
-    auto const refinementTolerance =
-        parset.get<double>("timeSteps.refinementTolerance");
-    auto const mustRefine = [&](UpdaterPair &coarseUpdater,
-                                UpdaterPair &fineUpdater) {
-      ScalarVector coarseAlpha;
-      coarseUpdater.first->extractAlpha(coarseAlpha);
-
-      ScalarVector fineAlpha;
-      fineUpdater.first->extractAlpha(fineAlpha);
-
-      return stateEnergyNorm.diff(fineAlpha, coarseAlpha) > refinementTolerance;
-    };
-
-    size_t timeStep = 1;
-
-    AdaptiveTimeStepper<NonlinearFactory, UpdaterPair> adaptiveTimeStepper(
-        factory, parset, myGlobalFriction, current, computeExternalForces,
-        mustRefine);
-    while (!adaptiveTimeStepper.reachedEnd()) {
-      adaptiveTimeStepper.advance();
-      reportTimeStep(adaptiveTimeStepper.getRelativeTime(),
-                     adaptiveTimeStepper.getRelativeTau());
-
-      Vector u, ur, vr;
-      ScalarVector alpha;
-      current.second->extractDisplacement(u);
-      current.second->extractRelativeDisplacement(ur);
-      current.second->extractRelativeVelocity(vr);
-      current.first->extractAlpha(alpha);
-
-      report(ur, vr, alpha);
-      {
-        BasisGridFunction<typename MyAssembler::VertexBasis, Vector>
-        relativeVelocity(myAssembler.vertexBasis, vr);
-        BasisGridFunction<typename MyAssembler::VertexBasis, Vector>
-        relativeDisplacement(myAssembler.vertexBasis, ur);
-        specialVelocityWriter.write(relativeVelocity);
-        specialDisplacementWriter.write(relativeDisplacement);
-      }
-
-      if (parset.get<bool>("io.writeVTK")) {
-        ScalarVector stress;
-        myAssembler.assembleVonMisesStress(body.getYoungModulus(),
-                                           body.getPoissonRatio(), u, stress);
-        vtkWriter.write(timeStep, ur, vr, alpha, stress);
-      }
-      timeStep++;
-    }
-    timeStepWriter.close();
-    Python::stop();
-  }
-  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/fixedpointiterator.cc → src/spatial-solving/fixedpointiterator.cc

@@ -5,18 +5,25 @@
 #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 "../enums.hh"
+#include "../enumparser.hh"
 
 #include "fixedpointiterator.hh"
 
-template <class Factory, class StateUpdater, class VelocityUpdater>
-FixedPointIterator<Factory, StateUpdater, VelocityUpdater>::FixedPointIterator(
+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)
-    : factory_(factory),
+    std::shared_ptr<Nonlinearity> globalFriction, ErrorNorm const &errorNorm)
+    : step_(factory.getStep()),
       parset_(parset),
       globalFriction_(globalFriction),
       fixedPointMaxIterations_(parset.get<size_t>("v.fpi.maximumIterations")),
@@ -24,57 +31,70 @@ FixedPointIterator<Factory, StateUpdater, VelocityUpdater>::FixedPointIterator(
       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")) {}
-
-template <class Factory, class StateUpdater, class VelocityUpdater>
-int FixedPointIterator<Factory, StateUpdater, VelocityUpdater>::run(
-    std::shared_ptr<StateUpdater> stateUpdater,
-    std::shared_ptr<VelocityUpdater> velocityUpdater,
-    Matrix const &velocityMatrix, Norm<Vector> const &velocityMatrixNorm,
-    Vector const &velocityRHS, Vector &velocityIterate) {
-  auto multigridStep = factory_.getSolver();
-
-  LoopSolver<Vector> velocityProblemSolver(
-      multigridStep, velocityMaxIterations_, velocityTolerance_,
-      &velocityMatrixNorm, verbosity_, false); // absolute error
-
-  Vector previousVelocityIterate = velocityIterate;
+      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;
-  for (fixedPointIteration = 1; fixedPointIteration <= fixedPointMaxIterations_;
+  size_t multigridIterations = 0;
+  ScalarVector alpha;
+  updaters.state_->extractAlpha(alpha);
+  for (fixedPointIteration = 0; fixedPointIteration < fixedPointMaxIterations_;
        ++fixedPointIteration) {
-    Vector v_m;
-    velocityUpdater->extractOldVelocity(v_m);
-    v_m *= 1.0 - lambda_;
-    Arithmetic::addProduct(v_m, lambda_, velocityIterate);
-
-    // solve a state problem
-    stateUpdater->solve(v_m);
-    ScalarVector alpha;
-    stateUpdater->extractAlpha(alpha);
-
     // solve a velocity problem
     globalFriction_->updateAlpha(alpha);
     ConvexProblem convexProblem(1.0, velocityMatrix, *globalFriction_,
                                 velocityRHS, velocityIterate);
     BlockProblem velocityProblem(parset_, convexProblem);
-    multigridStep->setProblem(velocityIterate, velocityProblem);
+    step_->setProblem(velocityIterate, velocityProblem);
     velocityProblemSolver.preprocess();
     velocityProblemSolver.solve();
 
-    if (velocityMatrixNorm.diff(previousVelocityIterate, velocityIterate) <
-        fixedPointTolerance_)
-      break;
+    multigridIterations += velocityProblemSolver.getResult().iterations;
+
+    Vector v_m;
+    updaters.rate_->extractOldVelocity(v_m);
+    v_m *= 1.0 - lambda_;
+    Arithmetic::addProduct(v_m, lambda_, velocityIterate);
 
-    previousVelocityIterate = 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");
 
-  velocityUpdater->postProcess(velocityIterate);
-  velocityUpdater->postProcessRelativeQuantities();
+  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;
+}
 
-  return fixedPointIteration;
+std::ostream &operator<<(std::ostream &stream,
+                         FixedPointIterationCounter const &fpic) {
+  return stream << "(" << fpic.iterations << "," << fpic.multigridIterations
+                << ")";
 }
 
 #include "fixedpointiterator_tmpl.cc"

src/fixedpointiterator.hh → src/spatial-solving/fixedpointiterator.hh

-#ifndef SRC_FIXEDPOINTITERATOR_HH
-#define SRC_FIXEDPOINTITERATOR_HH
+#ifndef SRC_SPATIAL_SOLVING_FIXEDPOINTITERATOR_HH
+#define SRC_SPATIAL_SOLVING_FIXEDPOINTITERATOR_HH
 
 #include <memory>
 
@@ -8,9 +8,19 @@
 #include <dune/solvers/norms/norm.hh>
 #include <dune/solvers/solvers/solver.hh>
 
-template <class Factory, class StateUpdater, class VelocityUpdater>
+struct FixedPointIterationCounter {
+  size_t iterations = 0;
+  size_t multigridIterations = 0;
+
+  void operator+=(FixedPointIterationCounter const &other);
+};
+
+std::ostream &operator<<(std::ostream &stream,
+                         FixedPointIterationCounter const &fpic);
+
+template <class Factory, class Updaters, class ErrorNorm>
 class FixedPointIterator {
-  using ScalarVector = typename StateUpdater::ScalarVector;
+  using ScalarVector = typename Updaters::StateUpdater::ScalarVector;
   using Vector = typename Factory::Vector;
   using Matrix = typename Factory::Matrix;
   using ConvexProblem = typename Factory::ConvexProblem;
@@ -19,15 +29,16 @@ class FixedPointIterator {
 
 public:
   FixedPointIterator(Factory &factory, Dune::ParameterTree const &parset,
-                     std::shared_ptr<Nonlinearity> globalFriction);
+                     std::shared_ptr<Nonlinearity> globalFriction,
+                     ErrorNorm const &errorNorm_);
 
-  int run(std::shared_ptr<StateUpdater> stateUpdater,
-          std::shared_ptr<VelocityUpdater> velocityUpdater,
-          Matrix const &velocityMatrix, Norm<Vector> const &velocityMatrixNorm,
-          Vector const &velocityRHS, Vector &velocityIterate);
+  FixedPointIterationCounter run(Updaters updaters,
+                                 Matrix const &velocityMatrix,
+                                 Vector const &velocityRHS,
+                                 Vector &velocityIterate);
 
 private:
-  Factory &factory_;
+  std::shared_ptr<typename Factory::Step> step_;
   Dune::ParameterTree const &parset_;
   std::shared_ptr<Nonlinearity> globalFriction_;
 
@@ -37,5 +48,6 @@ class FixedPointIterator {
   size_t velocityMaxIterations_;
   double velocityTolerance_;
   Solver::VerbosityMode verbosity_;
+  ErrorNorm const &errorNorm_;
 };
 #endif

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/solverfactory.cc → src/spatial-solving/solverfactory.cc

@@ -13,15 +13,16 @@
 
 template <size_t dim, class BlockProblem, class Grid>
 SolverFactory<dim, BlockProblem, Grid>::SolverFactory(
-    Dune::ParameterTree const &parset, size_t refinements, Grid const &grid,
+    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(refinements),
-      multigridStep(new Solver(linearIterationStep, nonlinearSmoother)) {
+      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"),
@@ -47,12 +48,11 @@ template <size_t dim, class BlockProblem, class Grid>
 SolverFactory<dim, BlockProblem, Grid>::~SolverFactory() {
   for (auto &&x : transferOperators)
     delete x;
-
-  delete multigridStep;
 }
 
 template <size_t dim, class BlockProblem, class Grid>
-auto SolverFactory<dim, BlockProblem, Grid>::getSolver() -> Solver *{
+auto SolverFactory<dim, BlockProblem, Grid>::getStep()
+    -> std::shared_ptr<Step> {
   return multigridStep;
 }
 

src/solverfactory.hh → src/spatial-solving/solverfactory.hh

-#ifndef SRC_SOLVERFACTORY_HH
-#define SRC_SOLVERFACTORY_HH
+#ifndef SRC_SPATIAL_SOLVING_SOLVERFACTORY_HH
+#define SRC_SPATIAL_SOLVING_SOLVERFACTORY_HH
 
 #include <dune/common/bitsetvector.hh>
 #include <dune/common/parametertree.hh>
 
-#pragma clang diagnostic push
-#pragma clang diagnostic ignored "-Wsign-compare"
 #include <dune/solvers/iterationsteps/multigridstep.hh>
-#pragma clang diagnostic pop
-
 #include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
 #include <dune/solvers/norms/energynorm.hh>
 #include <dune/solvers/solvers/loopsolver.hh>
@@ -26,16 +22,17 @@ class SolverFactory {
 
 private:
   using NonlinearSmoother = GenericNonlinearGS<BlockProblem>;
-  using Solver =
-      TruncatedNonsmoothNewtonMultigrid<BlockProblem, NonlinearSmoother>;
 
 public:
-  SolverFactory(Dune::ParameterTree const &parset, size_t refinements,
-                Grid const &grid, Dune::BitSetVector<dim> const &ignoreNodes);
+  using Step =
+      TruncatedNonsmoothNewtonMultigrid<BlockProblem, NonlinearSmoother>;
+
+  SolverFactory(Dune::ParameterTree const &parset, Grid const &grid,
+                Dune::BitSetVector<dim> const &ignoreNodes);
 
   ~SolverFactory();
 
-  Solver *getSolver();
+  std::shared_ptr<Step> getStep();
 
 private:
   TruncatedBlockGSStep<Matrix, Vector> linearBaseSolverStep;
@@ -46,6 +43,6 @@ class SolverFactory {
   MultigridStep<Matrix, Vector> linearIterationStep;
   std::vector<CompressedMultigridTransfer<Vector> *> transferOperators;
   NonlinearSmoother nonlinearSmoother;
-  Solver *multigridStep;
+  std::shared_ptr<Step> multigridStep;
 };
 #endif

src/solverfactory_tmpl.cc → src/spatial-solving/solverfactory_tmpl.cc

@@ -2,8 +2,8 @@
 #error MY_DIM unset
 #endif
 
-#include "explicitgrid.hh"
-#include "explicitvectors.hh"
+#include "../explicitgrid.hh"
+#include "../explicitvectors.hh"
 
 #include <dune/tnnmg/nonlinearities/zerononlinearity.hh>
 #include <dune/tnnmg/problem-classes/blocknonlineartnnmgproblem.hh>
@@ -17,6 +17,6 @@ template class SolverFactory<
     MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
     Grid>;
 template class SolverFactory<
-    MY_DIM, BlockNonlinearTNNMGProblem<
-             ConvexProblem<ZeroNonlinearity<LocalVector, LocalMatrix>, Matrix>>,
+    MY_DIM, BlockNonlinearTNNMGProblem<ConvexProblem<
+                ZeroNonlinearity<LocalVector, LocalMatrix>, Matrix>>,
     Grid>;

src/state_tmpl.cc

-#include "explicitvectors.hh"
-
-template std::shared_ptr<StateUpdater<ScalarVector, Vector>> initStateUpdater<
-    ScalarVector, Vector>(Config::stateModel model,
-                          ScalarVector const &alpha_initial,
-                          Dune::BitSetVector<1> const &frictionalNodes,
-                          double L, double V0);

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.hh

+#ifndef SRC_TIME_STEPPING_ADAPTIVETIMESTEPPER_HH
+#define SRC_TIME_STEPPING_ADAPTIVETIMESTEPPER_HH
+
+#include <fstream>
+
+#include "coupledtimestepper.hh"
+
+struct IterationRegister {
+  void registerCount(FixedPointIterationCounter count);
+  void registerFinalCount(FixedPointIterationCounter count);
+
+  void reset();
+
+  FixedPointIterationCounter totalCount;
+  FixedPointIterationCounter finalCount;
+};
+
+template <class Factory, class Updaters, class ErrorNorm>
+class AdaptiveTimeStepper {
+  struct UpdatersWithCount {
+    Updaters updaters;
+    FixedPointIterationCounter count;
+  };
+
+  using Vector = typename Factory::Vector;
+  using ConvexProblem = typename Factory::ConvexProblem;
+  using Nonlinearity = typename ConvexProblem::NonlinearityType;
+
+  using MyCoupledTimeStepper = CoupledTimeStepper<Factory, Updaters, ErrorNorm>;
+
+public:
+  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);
+
+  bool reachedEnd();
+  IterationRegister advance();
+
+  double relativeTime_;
+  double relativeTau_;
+
+private:
+  UpdatersWithCount step(Updaters const &oldUpdaters, double rTime,
+                         double rTau);
+
+  double finalTime_;
+  Factory &factory_;
+  Dune::ParameterTree const &parset_;
+  std::shared_ptr<Nonlinearity> globalFriction_;
+  Updaters &current_;
+  UpdatersWithCount R1_;
+  std::function<void(double, Vector &)> externalForces_;
+  std::function<bool(Updaters &, Updaters &)> mustRefine_;
+  ErrorNorm const &errorNorm_;
+
+  IterationRegister iterationRegister_;
+};
+#endif

src/fixedpointiterator_tmpl.cc → src/time-stepping/adaptivetimestepper_tmpl.cc

@@ -2,25 +2,31 @@
 #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 "explicitgrid.hh"
-#include "explicitvectors.hh"
-
-#include "solverfactory.hh"
+#include "../spatial-solving/solverfactory.hh"
+#include "rate/rateupdater.hh"
 #include "state/stateupdater.hh"
-#include "timestepping.hh"
+#include "updaters.hh"
 
 using Function = Dune::VirtualFunction<double, double>;
 using Factory = SolverFactory<
-    MY_DIM, MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
+    MY_DIM,
+    MyBlockProblem<ConvexProblem<GlobalFriction<Matrix, Vector>, Matrix>>,
     Grid>;
 using MyStateUpdater = StateUpdater<ScalarVector, Vector>;
-using MyVelocityUpdater = TimeSteppingScheme<Vector, Matrix, Function, MY_DIM>;
+using MyRateUpdater = RateUpdater<Vector, Matrix, Function, MY_DIM>;
+using MyUpdaters = Updaters<MyRateUpdater, MyStateUpdater>;
+
+using ErrorNorm = EnergyNorm<ScalarMatrix, ScalarVector>;
 
-template class FixedPointIterator<Factory, MyStateUpdater, MyVelocityUpdater>;
+template class AdaptiveTimeStepper<Factory, MyUpdaters, ErrorNorm>;

src/coupledtimestepper.cc → src/time-stepping/coupledtimestepper.cc

@@ -2,31 +2,28 @@
 #include "config.h"
 #endif
 
-#include <dune/solvers/norms/energynorm.hh>
-
 #include "coupledtimestepper.hh"
-#include "fixedpointiterator.hh"
 
-template <class Factory, class StateUpdater, class VelocityUpdater>
-CoupledTimeStepper<Factory, StateUpdater, VelocityUpdater>::CoupledTimeStepper(
+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,
-    std::shared_ptr<StateUpdater> stateUpdater,
-    std::shared_ptr<VelocityUpdater> velocityUpdater,
+    std::shared_ptr<Nonlinearity> globalFriction, Updaters updaters,
+    ErrorNorm const &errorNorm,
     std::function<void(double, Vector &)> externalForces)
     : finalTime_(finalTime),
       factory_(factory),
       parset_(parset),
       globalFriction_(globalFriction),
-      stateUpdater_(stateUpdater),
-      velocityUpdater_(velocityUpdater),
-      externalForces_(externalForces) {}
+      updaters_(updaters),
+      externalForces_(externalForces),
+      errorNorm_(errorNorm) {}
 
-template <class Factory, class StateUpdater, class VelocityUpdater>
-int CoupledTimeStepper<Factory, StateUpdater, VelocityUpdater>::step(
-    double relativeTime, double relativeTau) {
-  stateUpdater_->nextTimeStep();
-  velocityUpdater_->nextTimeStep();
+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;
@@ -37,16 +34,13 @@ int CoupledTimeStepper<Factory, StateUpdater, VelocityUpdater>::step(
   Vector velocityIterate;
 
   auto const tau = relativeTau * finalTime_;
-  stateUpdater_->setup(tau);
-  velocityUpdater_->setup(ell, tau, newRelativeTime, velocityRHS,
-                          velocityIterate, velocityMatrix);
-  EnergyNorm<Matrix, Vector> const velocityMatrixNorm(velocityMatrix);
-
-  FixedPointIterator<Factory, StateUpdater, VelocityUpdater> fixedPointIterator(
-      factory_, parset_, globalFriction_);
-  auto const iterations =
-      fixedPointIterator.run(stateUpdater_, velocityUpdater_, velocityMatrix,
-                             velocityMatrixNorm, velocityRHS, velocityIterate);
+  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;
 }