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

@@ -11,24 +11,6 @@ namespace {
   using LocalMatrix2D = Dune::FieldMatrix<double, 2, 2>;
 
   using LocalVector = Dune::FieldVector<double, MY_DIM>;
-
-  // kludge because fieldvectors have no initialiser_list constructor, see
-  // https://dune-project.org/flyspray/index.php?do=details&task_id=1166
-  LocalVector2D createVector(double x, double y) {
-    LocalVector2D ret(0);
-    ret[0] = x;
-    ret[1] = y;
-    return ret;
-  }
-
-  LocalMatrix2D createMatrix(double a11, double a12, double a21, double a22) {
-    LocalMatrix2D ret(0);
-    ret[0][0] = a11;
-    ret[0][1] = a12;
-    ret[1][0] = a21;
-    ret[1][1] = a22;
-    return ret;
-  }
 }
 
 namespace reference {
@@ -42,34 +24,31 @@ namespace reference {
 
   double const zDistance = 0.35;
 
-  LocalVector2D const A = createVector(0, 0);
-  LocalVector2D const B = createVector(leftLeg, -rightLeg);
-  LocalVector2D const C = createVector(leftLeg, 0);
+  LocalVector2D const A = {0, 0};
+  LocalVector2D const B = {leftLeg, -rightLeg};
+  LocalVector2D const C = {leftLeg, 0};
 
-  LocalVector2D const Z = createVector(zDistance * s, 0);
-  LocalVector2D const Y =
-      createVector(zDistance * s, -zDistance *s / leftLeg * rightLeg);
-  LocalVector2D const X = createVector(Y[0] - weakLen * std::cos(leftAngle),
-                                       Y[1] + weakLen * std::sin(leftAngle));
+  LocalVector2D const Z = {zDistance * s, 0};
+  LocalVector2D const Y = {zDistance * s, -zDistance *s / leftLeg *rightLeg};
+  LocalVector2D const X = {Y[0] - weakLen * std::cos(leftAngle),
+                           Y[1] + weakLen *std::sin(leftAngle)};
 
-  LocalVector2D const U = createVector(X[0], 0);
+  LocalVector2D const U = {X[0], 0};
 
-  LocalVector2D const K =
-      createVector(B[0] - leftLeg * viscoHeight / rightLeg, B[1] + viscoHeight);
-  LocalVector2D const M = createVector(B[0], B[1] + viscoHeight);
+  LocalVector2D const K = {B[0] - leftLeg * viscoHeight / rightLeg,
+                           B[1] + viscoHeight};
+  LocalVector2D const M = {B[0], B[1] + viscoHeight};
 
   LocalVector2D const G = midPoint(A, X);
   LocalVector2D const H = midPoint(X, Y);
   LocalVector2D const J = midPoint(Y, B);
 
-  LocalVector2D const I = createVector(Y[0] + G[0], Y[1] + G[1]);
+  LocalVector2D const I = {Y[0] + G[0], Y[1] + G[1]};
 
-  LocalVector2D const zenith = createVector(0, 1);
-  LocalVector2D const orthogonalZenith = createVector(-1, 0);
+  LocalVector2D const zenith = {0, 1};
 
-  LocalMatrix2D const rotation =
-      createMatrix(std::cos(leftAngle), -std::sin(leftAngle),
-                   std::sin(leftAngle), std::cos(leftAngle));
+  LocalMatrix2D const rotation = {{std::cos(leftAngle), -std::sin(leftAngle)},
+                                  {std::sin(leftAngle), std::cos(leftAngle)}};
 }
 
 namespace {
@@ -102,10 +81,6 @@ LocalVector const Y = rotate(reference::Y);
 LocalVector const Z = rotate(reference::Z);
 
 LocalVector const zenith = rotate(reference::zenith);
-LocalVector const orthogonalZenith = rotate(reference::orthogonalZenith);
-
-double verticalProjection(LocalVector const &);
-double horizontalProjection(LocalVector const &);
 
 void write();
 

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

-#ifndef SRC_SAND_WEDGE_DATA_MYGLOBALFRICTIONDATA_HH
-#define SRC_SAND_WEDGE_DATA_MYGLOBALFRICTIONDATA_HH
+#ifndef SRC_ONE_BODY_PROBLEM_DATA_MYGLOBALFRICTIONDATA_HH
+#define SRC_ONE_BODY_PROBLEM_DATA_MYGLOBALFRICTIONDATA_HH
 
 #include <dune/common/function.hh>
 

src/sand-wedge-data/mygrid.cc → src/one-body-problem-data/mygrid.cc

@@ -2,9 +2,11 @@
 #include "config.h"
 #endif
 
+#include <dune/fufem/geometry/polyhedrondistance.hh>
+
 #include "mygrid.hh"
 #include "midpoint.hh"
-#include "../distances.hh"
+#include "../diameter.hh"
 
 #if MY_DIM == 3
 SimplexManager::SimplexManager(unsigned int shift) : shift_(shift) {}
@@ -165,7 +167,6 @@ MyFaces<GridView>::MyFaces(GridView const &gridView)
       upper(gridView)
 #endif
 {
-  using Grid = typename GridView::Grid;
   assert(isClose(MyGeometry::A[1], 0));
   assert(isClose(MyGeometry::B[1], 0));
   lower.insertFacesByProperty([&](typename GridView::Intersection const &in) {
@@ -200,12 +201,11 @@ void refine(Grid &grid, ConvexPolyhedron<LocalVector> const &weakPatch,
   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();
+    for (auto &&e : elements(grid.leafGridView())) {
+      auto const geometry = e.geometry();
 
       auto const weakeningRegionDistance =
-          distanceToWeakeningRegion(geometry, weakPatch);
+          distance(weakPatch, geometry, 1e-6 * MyGeometry::lengthScale);
       auto const admissibleDiameter =
           computeAdmissibleDiameter(weakeningRegionDistance, smallestDiameter);
 
@@ -213,7 +213,7 @@ void refine(Grid &grid, ConvexPolyhedron<LocalVector> const &weakPatch,
         continue;
 
       needRefine = true;
-      grid.mark(1, *it);
+      grid.mark(1, e);
     }
     if (!needRefine)
       break;

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

-#ifndef SRC_SAND_WEDGE_DATA_MYGRID_HH
-#define SRC_SAND_WEDGE_DATA_MYGRID_HH
-
-#include <boost/range/irange.hpp>
+#ifndef SRC_ONE_BODY_PROBLEM_DATA_MYGRID_HH
+#define SRC_ONE_BODY_PROBLEM_DATA_MYGRID_HH
 
 #include <dune/common/fmatrix.hh>
 #include <dune/grid/common/gridfactory.hh>
 
-#pragma clang diagnostic push
-#pragma clang diagnostic ignored "-Wdeprecated-declarations"
 #include <dune/fufem/boundarypatch.hh>
-#pragma clang diagnostic pop
-
-#include <dune/tectonic/polyhedrondistance.hh>
+#include <dune/fufem/geometry/convexpolyhedron.hh>
 
 #include "mygeometry.hh"
 

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/tectonic/polyhedrondistance.hh>
+#include <dune/fufem/geometry/polyhedrondistance.hh>
 
 class PatchFunction
     : public Dune::VirtualFunction<Dune::FieldVector<double, MY_DIM>,

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,7 +7,7 @@
 
 #include "mygeometry.hh"
 
-class TwoPieceFunction
+class SegmentedFunction
     : public Dune::VirtualFunction<Dune::FieldVector<double, MY_DIM>,
                                    Dune::FieldVector<double, 1>> {
 private:
@@ -24,7 +24,7 @@ class TwoPieceFunction
   double const _v2;
 
 public:
-  TwoPieceFunction(double v1, double v2) : _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/sand-wedge-data/weakpatch.hh → src/one-body-problem-data/weakpatch.hh

-#ifndef SRC_SAND_WEDGE_DATA_WEAKPATCH_HH
-#define SRC_SAND_WEDGE_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() {
+template <class LocalVector>
+ConvexPolyhedron<LocalVector> getWeakPatch(Dune::ParameterTree const &parset) {
   ConvexPolyhedron<LocalVector> weakPatch;
 #if MY_DIM == 3
   weakPatch.vertices.resize(4);
@@ -11,6 +12,16 @@ template <class LocalVector> ConvexPolyhedron<LocalVector> getWeakPatch() {
     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;

src/sand-wedge-data/parset.cfg → src/one-body-problem.cfg

@@ -2,11 +2,15 @@
 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     = 0.5e5 # [Pa]
@@ -21,49 +25,41 @@ shearViscosity  = 1e4   # [Pas]
 bulkViscosity   = 1e4   # [Pas]
 
 [boundary.friction]
-smallestDiameter= 2e-3  # [m]
 C               = 10    # [Pa]
 mu0             = 0.7   # [ ]
 V0              = 5e-5  # [m/s]
-L               = 2.5e-5# [m]
+L               = 2.25e-5 # [m]
 initialAlpha    = 0     # [ ]
 stateModel      = AgeingLaw
-frictionModel   = Truncated
+frictionModel   = Regularised
 [boundary.friction.weakening]
 a               = 0.002 # [ ]
-b               = 0.014 # [ ]
+b               = 0.017 # [ ]
 [boundary.friction.strengthening]
-a               = 0.025 # [ ]
+a               = 0.020 # [ ]
 b               = 0.005 # [ ]
 
 [timeSteps]
-refinementTolerance = 1e-5
-number = 100000
 scheme = newmark
 
 [u0.solver]
-tolerance         = 1e-8
 maximumIterations = 100000
 verbosity         = quiet
 
 [a0.solver]
-tolerance         = 1e-8
 maximumIterations = 100000
 verbosity         = quiet
 
 [v.solver]
-tolerance         = 1e-8
 maximumIterations = 100000
 verbosity         = quiet
 
 [v.fpi]
-tolerance         = 1e-5
 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/fixedpointiterator.cc → src/spatial-solving/fixedpointiterator.cc

@@ -8,18 +8,22 @@
 #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,
-          class ErrorNorm>
-FixedPointIterator<Factory, StateUpdater, VelocityUpdater, ErrorNorm>::
-    FixedPointIterator(Factory &factory, Dune::ParameterTree const &parset,
-                       std::shared_ptr<Nonlinearity> globalFriction,
-                       ErrorNorm const &errorNorm)
-    : factory_(factory),
+void FixedPointIterationCounter::operator+=(
+    FixedPointIterationCounter const &other) {
+  iterations += other.iterations;
+  multigridIterations += other.multigridIterations;
+}
+
+template <class Factory, class Updaters, class ErrorNorm>
+FixedPointIterator<Factory, Updaters, ErrorNorm>::FixedPointIterator(
+    Factory &factory, Dune::ParameterTree const &parset,
+    std::shared_ptr<Nonlinearity> globalFriction, ErrorNorm const &errorNorm)
+    : step_(factory.getStep()),
       parset_(parset),
       globalFriction_(globalFriction),
       fixedPointMaxIterations_(parset.get<size_t>("v.fpi.maximumIterations")),
@@ -30,25 +34,20 @@ FixedPointIterator<Factory, StateUpdater, VelocityUpdater, ErrorNorm>::
       verbosity_(parset.get<Solver::VerbosityMode>("v.solver.verbosity")),
       errorNorm_(errorNorm) {}
 
-template <class Factory, class StateUpdater, class VelocityUpdater,
-          class ErrorNorm>
+template <class Factory, class Updaters, class ErrorNorm>
 FixedPointIterationCounter
-FixedPointIterator<Factory, StateUpdater, VelocityUpdater, ErrorNorm>::run(
-    std::shared_ptr<StateUpdater> stateUpdater,
-    std::shared_ptr<VelocityUpdater> velocityUpdater,
-    Matrix const &velocityMatrix, Vector const &velocityRHS,
+FixedPointIterator<Factory, Updaters, ErrorNorm>::run(
+    Updaters updaters, Matrix const &velocityMatrix, Vector const &velocityRHS,
     Vector &velocityIterate) {
-  auto multigridStep = factory_.getSolver();
-
   EnergyNorm<Matrix, Vector> energyNorm(velocityMatrix);
-  LoopSolver<Vector> velocityProblemSolver(
-      multigridStep, velocityMaxIterations_, velocityTolerance_, &energyNorm,
-      verbosity_, false); // absolute error
+  LoopSolver<Vector> velocityProblemSolver(step_.get(), velocityMaxIterations_,
+                                           velocityTolerance_, &energyNorm,
+                                           verbosity_, false); // absolute error
 
   size_t fixedPointIteration;
   size_t multigridIterations = 0;
   ScalarVector alpha;
-  stateUpdater->extractAlpha(alpha);
+  updaters.state_->extractAlpha(alpha);
   for (fixedPointIteration = 0; fixedPointIteration < fixedPointMaxIterations_;
        ++fixedPointIteration) {
     // solve a velocity problem
@@ -56,21 +55,21 @@ FixedPointIterator<Factory, StateUpdater, VelocityUpdater, ErrorNorm>::run(
     ConvexProblem convexProblem(1.0, velocityMatrix, *globalFriction_,
                                 velocityRHS, velocityIterate);
     BlockProblem velocityProblem(parset_, convexProblem);
-    multigridStep->setProblem(velocityIterate, velocityProblem);
+    step_->setProblem(velocityIterate, velocityProblem);
     velocityProblemSolver.preprocess();
     velocityProblemSolver.solve();
 
     multigridIterations += velocityProblemSolver.getResult().iterations;
 
     Vector v_m;
-    velocityUpdater->extractOldVelocity(v_m);
+    updaters.rate_->extractOldVelocity(v_m);
     v_m *= 1.0 - lambda_;
     Arithmetic::addProduct(v_m, lambda_, velocityIterate);
 
     // solve a state problem
-    stateUpdater->solve(v_m);
+    updaters.state_->solve(v_m);
     ScalarVector newAlpha;
-    stateUpdater->extractAlpha(newAlpha);
+    updaters.state_->extractAlpha(newAlpha);
 
     if (lambda_ < 1e-12 or
         errorNorm_.diff(alpha, newAlpha) < fixedPointTolerance_) {
@@ -82,9 +81,14 @@ FixedPointIterator<Factory, StateUpdater, VelocityUpdater, ErrorNorm>::run(
   if (fixedPointIteration == fixedPointMaxIterations_)
     DUNE_THROW(Dune::Exception, "FPI failed to converge");
 
-  velocityUpdater->postProcess(velocityIterate);
+  updaters.rate_->postProcess(velocityIterate);
 
-  return { fixedPointIteration, multigridIterations };
+  // Cannot use return { fixedPointIteration, multigridIterations };
+  // with gcc 4.9.2, see also http://stackoverflow.com/a/37777814/179927
+  FixedPointIterationCounter ret;
+  ret.iterations = fixedPointIteration;
+  ret.multigridIterations = multigridIterations;
+  return ret;
 }
 
 std::ostream &operator<<(std::ostream &stream,

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>
 
@@ -9,17 +9,18 @@
 #include <dune/solvers/solvers/solver.hh>
 
 struct FixedPointIterationCounter {
-  size_t iterations;
-  size_t multigridIterations;
+  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 StateUpdater, class VelocityUpdater,
-          class ErrorNorm>
+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;
@@ -31,14 +32,13 @@ class FixedPointIterator {
                      std::shared_ptr<Nonlinearity> globalFriction,
                      ErrorNorm const &errorNorm_);
 
-  FixedPointIterationCounter run(
-      std::shared_ptr<StateUpdater> stateUpdater,
-      std::shared_ptr<VelocityUpdater> velocityUpdater,
-      Matrix const &velocityMatrix, 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_;
 

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>

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"