[Kill] Remove the sliding block problem

src/Makefile.am

-bin_PROGRAMS = sand-wedge sliding-block-2D
+bin_PROGRAMS = sand-wedge
 
 common_sources = \
 	assemblers.cc \
@@ -12,10 +12,6 @@ sand_wedge_SOURCES = $(common_sources) sand-wedge.cc
 sand_wedge_CPPFLAGS = \
 	$(AM_CPPFLAGS) $(BOOST_CPPFLAGS) \
 	-Ddatadir=\"$(abs_srcdir)/sand-wedge-data/\" -DDIM=2
-sliding_block_2D_SOURCES = $(common_sources) sliding-block.cc
-sliding_block_2D_CPPFLAGS = \
-	$(AM_CPPFLAGS) \
-	-Ddatadir=\"$(abs_srcdir)/sliding-block-data/\" -DDIM=2
 
 # Some are for clang, others are for gcc
 AM_CXXFLAGS = \

src/sliding-block-data/boundaryconditions.py

-class neumannCondition:
-    def __call__(self, relativeTime):
-        return 0
-
-class velocityDirichletCondition:
-    def __call__(self, relativeTime):
-        if relativeTime <= 0.25:
-            factor = 4*relativeTime;
-        else:
-            factor = 1
-        return 2e-4 * factor
-
-Functions = {
-    'neumannCondition' : neumannCondition(),
-    'velocityDirichletCondition' : velocityDirichletCondition()
-}

src/sliding-block-data/mybody.hh

-#ifndef SRC_SLIDING_BLOCK_DATA_MYBODY_HH
-#define SRC_SLIDING_BLOCK_DATA_MYBODY_HH
-
-#include <dune/common/fvector.hh>
-
-#include <dune/fufem/functions/constantfunction.hh>
-
-#include <dune/tectonic/body.hh>
-#include <dune/tectonic/gravity.hh>
-
-#include "mygeometry.hh"
-
-template <int dimension> class MyBody : public Body<dimension> {
-  using typename Body<dimension>::ScalarFunction;
-  using typename Body<dimension>::VectorField;
-  using MyConstantFunction = ConstantFunction<
-      Dune::FieldVector<double, dimension>, Dune::FieldVector<double, 1>>;
-
-public:
-  MyBody(Dune::ParameterTree const &parset)
-      : poissonRatio_(parset.get<double>("body.poissonRatio")),
-        youngModulus_(parset.get<double>("body.youngModulus")),
-        shearViscosityField_(parset.get<double>("body.shearViscosity")),
-        bulkViscosityField_(parset.get<double>("body.bulkViscosity")),
-        densityField_(parset.get<double>("body.density")),
-        gravityField_(densityField_, MyGeometry::zenith,
-                      parset.get<double>("gravity")) {}
-
-  double getPoissonRatio() const override { return poissonRatio_; }
-
-  double getYoungModulus() const override { return youngModulus_; }
-
-  ScalarFunction const &getShearViscosityField() const override {
-    return shearViscosityField_;
-  }
-  ScalarFunction const &getBulkViscosityField() const override {
-    return bulkViscosityField_;
-  }
-
-  ScalarFunction const &getDensityField() const override {
-    return densityField_;
-  }
-  VectorField const &getGravityField() const override { return gravityField_; }
-
-private:
-  double const poissonRatio_;
-  double const youngModulus_;
-  MyConstantFunction const shearViscosityField_;
-  MyConstantFunction const bulkViscosityField_;
-  MyConstantFunction const densityField_;
-  Gravity<dimension> const gravityField_;
-};
-
-#endif

src/sliding-block-data/mygeometry.hh

-#ifndef SRC_SLIDING_BLOCK_DATA_MYGEOMETRY_HH
-#define SRC_SLIDING_BLOCK_DATA_MYGEOMETRY_HH
-
-#include <dune/common/fassign.hh>
-#include <dune/common/fvector.hh>
-
-namespace MyGeometry {
-namespace {
-  using LocalVector = Dune::FieldVector<double, 2>;
-
-  // kludge because fieldvectors have no initialiser_list constructor, see
-  // https://dune-project.org/flyspray/index.php?do=details&task_id=1166
-  LocalVector generateVector(double x, double y) {
-    LocalVector tmp;
-    tmp <<= x, y;
-    return tmp;
-  }
-}
-
-LocalVector const A = generateVector(0, 0);
-LocalVector const B = generateVector(5, 0);
-LocalVector const C = generateVector(5, 1);
-LocalVector const D = generateVector(0, 1);
-
-LocalVector const zenith = generateVector(0, 1);
-}
-#endif

src/sliding-block-data/myglobalfrictiondata.hh

-#ifndef SRC_SLIDING_BLOCK_DATA_MYGLOBALFRICTIONDATA_HH
-#define SRC_SLIDING_BLOCK_DATA_MYGLOBALFRICTIONDATA_HH
-
-#include <dune/common/function.hh>
-#include <dune/common/fvector.hh>
-
-#include <dune/fufem/functions/constantfunction.hh>
-
-#include <dune/tectonic/globalfrictiondata.hh>
-
-#include "mygeometry.hh"
-
-template <int dimension>
-class MyGlobalFrictionData : public GlobalFrictionData<dimension> {
-private:
-  using typename GlobalFrictionData<dimension>::VirtualFunction;
-
-public:
-  MyGlobalFrictionData(Dune::ParameterTree const &parset)
-      : C_(parset.get<double>("C")),
-        L_(parset.get<double>("L")),
-        V0_(parset.get<double>("V0")),
-        a_(parset.get<double>("a")),
-        b_(parset.get<double>("b")),
-        mu0_(parset.get<double>("mu0")) {}
-
-  double virtual const &C() const override { return C_; }
-  double virtual const &L() const override { return L_; }
-  double virtual const &V0() const override { return V0_; }
-  VirtualFunction virtual const &a() const override { return a_; }
-  VirtualFunction virtual const &b() const override { return b_; }
-  double virtual const &mu0() const override { return mu0_; }
-
-private:
-  using MyConstantFunction = ConstantFunction<
-      Dune::FieldVector<double, dimension>, Dune::FieldVector<double, 1>>;
-
-  double const C_;
-  double const L_;
-  double const V0_;
-  MyConstantFunction const a_;
-  MyConstantFunction const b_;
-  double const mu0_;
-};
-#endif

src/sliding-block-data/mygrid.hh

-#ifndef SRC_SLIDING_BLOCK_DATA_MYGRID_HH
-#define SRC_SLIDING_BLOCK_DATA_MYGRID_HH
-
-#include <dune/grid/common/gridfactory.hh>
-
-#include <dune/fufem/boundarypatch.hh>
-
-#include "mygeometry.hh"
-
-template <class Grid> std::shared_ptr<Grid> constructGrid() {
-  std::array<unsigned int, 2> elements = { { 5, 1 } };
-  return Dune::StructuredGridFactory<Grid>::createSimplexGrid(
-      MyGeometry::A, MyGeometry::C, elements);
-}
-
-template <class GridView> class MyFaces {
-private:
-  bool isClose(double a, double b) {
-    return std::abs(a - b) < 1e-14;
-  };
-
-public:
-  MyFaces(GridView const &gridView) : lower(gridView), upper(gridView) {
-    lower.insertFacesByProperty([&](typename GridView::Intersection const &in) {
-      return isClose(MyGeometry::A[1], in.geometry().center()[1]);
-    });
-
-    upper.insertFacesByProperty([&](typename GridView::Intersection const &in) {
-      return isClose(MyGeometry::C[1], in.geometry().center()[1]);
-    });
-  }
-
-  BoundaryPatch<GridView> lower;
-  BoundaryPatch<GridView> upper;
-};
-#endif

src/sliding-block-data/parset.cfg

-# -*- mode:conf -*-
-gravity         = 9.81
-
-[io]
-printProgress = false
-writeVTK      = false
-
-[problem]
-finalTime = 15
-
-[body]
-youngModulus = 5e7
-poissonRatio = 0.3
-density = 5000
-shearViscosity = 0
-bulkViscosity  = 0
-
-[boundary.friction]
-C            = 0.0
-mu0          = 0.6
-a            = 0.010
-b            = 0.015
-V0           = 1e-6
-L            = 1e-5
-initialState = 4.54e-05 # = exp(-10) = theta, so that alpha = -10
-stateModel   = Dieterich
-
-[timeSteps]
-number = 10000
-scheme = newmark
-
-[grid]
-refinements = 4
-
-[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
-relaxation        = 0.5
-requiredReduction = 0.5
-
-[solver.tnnmg.linear]
-maxiumumIterations = 100000
-tolerance          = 1e-10
-pre                = 3
-cycle              = 1  # 1 = V, 2 = W, etc.
-post               = 3
-
-[solver.tnnmg.main]
-pre   = 1
-multi = 5 # number of multigrid steps
-post  = 0
-
-[localsolver]
-steps = 1

src/sliding-block.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
-
-#ifndef DIM
-#error DIM unset
-#endif
-
-#if !HAVE_ALUGRID
-#error ALUGRID is required
-#endif
-
-#include <cmath>
-#include <exception>
-#include <fstream>
-#include <iostream>
-
-#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>
-
-#pragma clang diagnostic push
-#pragma clang diagnostic ignored "-Wignored-qualifiers"
-#include <dune/grid/alugrid.hh>
-#pragma clang diagnostic pop
-
-#include <dune/grid/common/mcmgmapper.hh>
-#include <dune/grid/utility/structuredgridfactory.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/functionspacebases/p1nodalbasis.hh>
-#include <dune/fufem/sharedpointermap.hh>
-#include <dune/solvers/norms/energynorm.hh>
-#include <dune/solvers/norms/sumnorm.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/globalnonlinearity.hh>
-
-#include "assemblers.hh"
-#include "enum_parser.cc"
-#include "enum_scheme.cc"
-#include "enum_state_model.cc"
-#include "enum_verbosity.cc"
-#include "enums.hh"
-#include "friction_writer.hh"
-#include "sliding-block-data/mybody.hh"
-#include "sliding-block-data/mygeometry.hh"
-#include "sliding-block-data/myglobalfrictiondata.hh"
-#include "sliding-block-data/mygrid.hh"
-#include "solverfactory.hh"
-#include "state.hh"
-#include "timestepping.hh"
-#include "vtk.hh"
-size_t const dims = DIM;
-
-void initPython() {
-  Python::start();
-
-  Python::run("import sys");
-  Python::run("sys.path.append('" datadir "')");
-}
-
-int main(int argc, char *argv[]) {
-  try {
-    Dune::ParameterTree parset;
-    Dune::ParameterTreeParser::readINITree(datadir "/parset.cfg", parset);
-    Dune::ParameterTreeParser::readOptions(argc, argv, parset);
-
-    // {{{ Set up grid
-    using Grid = Dune::ALUGrid<dims, dims, Dune::simplex, Dune::nonconforming>;
-    auto grid = constructGrid<Grid>();
-
-    auto const refinements = parset.get<size_t>("grid.refinements");
-    grid->globalRefine(refinements);
-    size_t const fineVertexCount = grid->size(grid->maxLevel(), dims);
-
-    using GridView = Grid::LeafGridView;
-    GridView const leafView = grid->leafView();
-    // }}}
-
-    // Set up myFaces
-    MyFaces<GridView> myFaces(leafView);
-
-    // Neumann boundary
-    BoundaryPatch<GridView> const neumannBoundary(leafView);
-
-    // Frictional Boundary
-    BoundaryPatch<GridView> const &frictionalBoundary = myFaces.lower;
-    Dune::BitSetVector<1> frictionalNodes(fineVertexCount);
-    frictionalBoundary.getVertices(frictionalNodes);
-
-    // Dirichlet Boundary
-    Dune::BitSetVector<dims> noNodes(fineVertexCount);
-    Dune::BitSetVector<dims> dirichletNodes(fineVertexCount);
-    for (size_t i = 0; i < fineVertexCount; ++i) {
-      if (myFaces.lower.containsVertex(i))
-        dirichletNodes[i][1] = true;
-
-      if (myFaces.upper.containsVertex(i))
-        dirichletNodes[i] = true;
-    }
-
-    // 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");
-
-    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;
-
-    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);
-    // Q: Does it make sense to weigh them in this manner?
-    SumNorm<Vector> const AMNorm(1.0, ANorm, 1.0, MNorm);
-
-    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
-    auto const computeExternalForces = [&](double _relativeTime, Vector &_ell) {
-      myAssembler.assembleNeumann(neumannBoundary, _ell, neumannFunction,
-                                  _relativeTime);
-      _ell += gravityFunctional;
-    };
-    Vector ell(fineVertexCount);
-    computeExternalForces(0.0, ell);
-
-    // {{{ Initial conditions
-    using LinearFactory = SolverFactory<
-        dims, BlockNonlinearTNNMGProblem<ConvexProblem<
-                  ZeroNonlinearity<LocalVector, LocalMatrix>, Matrix>>,
-        Grid>;
-    ZeroNonlinearity<LocalVector, LocalMatrix> zeroNonlinearity;
-
-    // TODO: clean up once generic lambdas arrive
-    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(fineVertexCount);
-    u_initial = 0.0;
-    solveLinearProblem(dirichletNodes, A, ell, u_initial, ANorm,
-                       parset.sub("u0.solver"));
-    Vector ur_initial(fineVertexCount);
-    ur_initial = 0.0;
-
-    ScalarVector alpha_initial(fineVertexCount);
-    alpha_initial =
-        std::log(parset.get<double>("boundary.friction.initialState"));
-    ScalarVector normalStress(fineVertexCount);
-    myAssembler.assembleNormalStress(frictionalBoundary, normalStress,
-                                     body.getYoungModulus(),
-                                     body.getPoissonRatio(), u_initial);
-    MyGlobalFrictionData<dims> frictionInfo(parset.sub("boundary.friction"));
-    auto myGlobalNonlinearity = myAssembler.assembleFrictionNonlinearity(
-        frictionalBoundary, frictionInfo, normalStress);
-    myGlobalNonlinearity->updateLogState(alpha_initial);
-
-    Vector v_initial(fineVertexCount);
-    v_initial = 0.0;
-    {
-      // Prescribe a homogeneous velocity field in the x-direction
-      // This way, the initial condition for v and the Dirichlet
-      // condition match up at t=0
-      double v_initial_const;
-      velocityDirichletFunction.evaluate(0.0, v_initial_const);
-      for (auto &x : v_initial)
-        x[0] = v_initial_const;
-    }
-    Vector vr_initial(fineVertexCount);
-    vr_initial = 0.0;
-
-    Vector a_initial(fineVertexCount);
-    a_initial = 0.0;
-    {
-      // We solve Ma = ell - [Au + Cv + Psi(v)]
-      Vector accelerationRHS(fineVertexCount);
-      {
-        accelerationRHS = 0.0;
-        Arithmetic::addProduct(accelerationRHS, A, u_initial);
-        Arithmetic::addProduct(accelerationRHS, C, v_initial);
-        // NOTE: We assume differentiability of Psi at v0 here!
-        myGlobalNonlinearity->addGradient(v_initial, accelerationRHS);
-        accelerationRHS *= -1.0;
-        accelerationRHS += ell;
-      }
-      solveLinearProblem(noNodes, M, accelerationRHS, a_initial, MNorm,
-                         parset.sub("a0.solver"));
-    }
-    // }}}
-
-    Vector vertexCoordinates(fineVertexCount);
-    {
-      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",
-        [](LocalVector const &x) { return x[0]; });
-    auto const reportFriction = [&](Vector const &_u, Vector const &_v,
-                                    ScalarVector const &_alpha) {
-      ScalarVector c;
-      myGlobalNonlinearity->coefficientOfFriction(_v, c);
-      frictionWriter.writeKinetics(_u, _v);
-      frictionWriter.writeOther(c, _alpha);
-    };
-    reportFriction(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);
-    }
-
-    // Set up TNNMG solver
-    using NonlinearFactory =
-        SolverFactory<dims, MyBlockProblem<ConvexProblem<
-                                GlobalNonlinearity<Matrix, Vector>, Matrix>>,
-                      Grid>;
-    NonlinearFactory factory(parset.sub("solver.tnnmg"), refinements, *grid,
-                             dirichletNodes);
-    auto multigridStep = factory.getSolver();
-
-    std::fstream iterationWriter("iterations", std::fstream::out),
-        relaxationWriter("relaxation", std::fstream::out);
-
-    auto timeSteppingScheme = initTimeStepper(
-        parset.get<Config::scheme>("timeSteps.scheme"),
-        velocityDirichletFunction, dirichletNodes, M, A, C, u_initial,
-        ur_initial, v_initial, vr_initial, a_initial);
-    auto stateUpdater = initStateUpdater<ScalarVector, Vector>(
-        parset.get<Config::stateModel>("boundary.friction.stateModel"),
-        alpha_initial, frictionalNodes,
-        parset.get<double>("boundary.friction.L"));
-
-    Vector v = v_initial;
-    Vector v_m(fineVertexCount);
-    ScalarVector alpha(fineVertexCount);
-
-    auto const timeSteps = parset.get<size_t>("timeSteps.number"),
-               maximumStateFPI = parset.get<size_t>("v.fpi.maximumIterations"),
-               maximumIterations =
-                   parset.get<size_t>("v.solver.maximumIterations");
-    auto const tau = parset.get<double>("problem.finalTime") / timeSteps,
-               tolerance = parset.get<double>("v.solver.tolerance"),
-               fixedPointTolerance = parset.get<double>("v.fpi.tolerance"),
-               relaxation = parset.get<double>("v.fpi.relaxation"),
-               requiredReduction =
-                   parset.get<double>("v.fpi.requiredReduction");
-    auto const printProgress = parset.get<bool>("io.printProgress");
-    auto const verbosity =
-        parset.get<Solver::VerbosityMode>("v.solver.verbosity");
-    for (size_t timeStep = 1; timeStep <= timeSteps; ++timeStep) {
-      if (printProgress)
-        std::cout << std::setw(7) << timeStep << " " << std::flush;
-
-      stateUpdater->nextTimeStep();
-      timeSteppingScheme->nextTimeStep();
-
-      auto const relativeTime = double(timeStep) / double(timeSteps);
-      computeExternalForces(relativeTime, ell);
-
-      Matrix velocityMatrix;
-      Vector velocityRHS(fineVertexCount);
-      Vector velocityIterate(fineVertexCount);
-
-      stateUpdater->setup(tau);
-      timeSteppingScheme->setup(ell, tau, relativeTime, velocityRHS,
-                                velocityIterate, velocityMatrix);
-
-      LoopSolver<Vector> velocityProblemSolver(multigridStep, maximumIterations,
-                                               tolerance, &AMNorm, verbosity,
-                                               false); // absolute error
-
-      size_t iterationCounter;
-      auto solveVelocityProblem = [&](Vector &_velocityIterate,
-                                      ScalarVector const &_alpha) {
-        myGlobalNonlinearity->updateLogState(_alpha);
-
-        // NIT: Do we really need to pass u here?
-        typename NonlinearFactory::ConvexProblem convexProblem(
-            1.0, velocityMatrix, *myGlobalNonlinearity, velocityRHS,
-            _velocityIterate);
-        typename NonlinearFactory::BlockProblem velocityProblem(parset,
-                                                                convexProblem);
-        multigridStep->setProblem(_velocityIterate, velocityProblem);
-
-        velocityProblemSolver.preprocess();
-        velocityProblemSolver.solve();
-        iterationCounter = velocityProblemSolver.getResult().iterations;
-      };
-
-      Vector u;
-      Vector v_saved;
-      ScalarVector alpha_saved;
-      double lastStateCorrection;
-      for (size_t stateFPI = 1; stateFPI <= maximumStateFPI; ++stateFPI) {
-        timeSteppingScheme->extractOldVelocity(v_m);
-        v_m *= 0.5;
-        Arithmetic::addProduct(v_m, 0.5, v);
-
-        stateUpdater->solve(v_m);
-        stateUpdater->extractLogState(alpha);
-
-        if (stateFPI == 1)
-          relaxationWriter << "N ";
-        else {
-          double const stateCorrection =
-              stateEnergyNorm.diff(alpha, alpha_saved);
-          if (stateFPI <= 2 // lastStateCorrection is only set for stateFPI > 2
-              or stateCorrection < requiredReduction * lastStateCorrection)
-            relaxationWriter << "N ";
-          else {
-            alpha *= (1.0 - relaxation);
-            Arithmetic::addProduct(alpha, relaxation, alpha_saved);
-            relaxationWriter << "Y ";
-          }
-          lastStateCorrection = stateCorrection;
-        }
-
-        solveVelocityProblem(velocityIterate, alpha);
-        timeSteppingScheme->postProcess(velocityIterate);
-        timeSteppingScheme->extractDisplacement(u);
-        timeSteppingScheme->extractVelocity(v);
-
-        iterationWriter << iterationCounter << " ";
-        if (printProgress)
-          std::cout << '.' << std::flush;
-
-        if (stateFPI > 1) {
-          double const velocityCorrection = AMNorm.diff(v_saved, v);
-          if (velocityCorrection < fixedPointTolerance)
-            break;
-        }
-        if (stateFPI == maximumStateFPI)
-          DUNE_THROW(Dune::Exception, "FPI failed to converge");
-
-        alpha_saved = alpha;
-        v_saved = v;
-      }
-      if (printProgress)
-        std::cout << std::endl;
-
-      Vector ur;
-      Vector vr;
-      timeSteppingScheme->postProcessRelativeQuantities();
-      timeSteppingScheme->extractRelativeDisplacement(ur);
-      timeSteppingScheme->extractRelativeVelocity(vr);
-
-      reportFriction(ur, vr, alpha);
-      iterationWriter << std::endl;
-      relaxationWriter << std::endl;
-
-      if (parset.get<bool>("io.writeVTK")) {
-        ScalarVector stress;
-        myAssembler.assembleVonMisesStress(body.getYoungModulus(),
-                                           body.getPoissonRatio(), u, stress);
-        vtkWriter.write(timeStep, ur, vr, alpha, stress);
-      }
-    }
-    iterationWriter.close();
-    relaxationWriter.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;
-  }
-}