dune/tectonic/quadraticenergy.hh

+#ifndef DUNE_TECTONIC_QUADRATICENERGY_HH
+#define DUNE_TECTONIC_QUADRATICENERGY_HH
+
+#include <memory>
+
+template <class NonlinearityTEMPLATE> class QuadraticEnergy {
+public:
+  using Nonlinearity = NonlinearityTEMPLATE;
+  using LocalVector = typename Nonlinearity::VectorType;
+
+  QuadraticEnergy(double alpha, LocalVector const &b, Nonlinearity const &phi)
+      : alpha(alpha), b(b), phi(phi) {}
+
+  double const alpha;
+  LocalVector const &b;
+  Nonlinearity const &phi;
+};
+#endif

dune/tectonic/tectonic.hh

-#ifndef DUNE_tectonic.hh
-#define DUNE_tectonic .hh
+#ifndef DUNE_TECTONIC_TECTONIC_HH
+#define DUNE_TECTONIC_TECTONIC_HH
 
 // add your classes here
-
-#endif // DUNE_tectonic.hh
+#endif

m4/Makefile.am

-M4FILES = dune-tectonic.m4
-
-aclocaldir = $(datadir)/aclocal
-aclocal_DATA = $(M4FILES)
-
-EXTRA_DIST = $(M4FILES)
-
-include $(top_srcdir)/am/global-rules

m4/dune-tectonic.m4

-dnl -*- autoconf -*-
-# Macros needed to find dune-tectonic and dependent libraries.  They are called by
-# the macros in ${top_src_dir}/dependencies.m4, which is generated by
-# "dunecontrol autogen"
-
-# Additional checks needed to build dune-tectonic
-# This macro should be invoked by every module which depends on dune-tectonic, as
-# well as by dune-tectonic itself
-AC_DEFUN([DUNE_TECTONIC_CHECKS])
-
-# Additional checks needed to find dune-tectonic
-# This macro should be invoked by every module which depends on dune-tectonic, but
-# not by dune-tectonic itself
-AC_DEFUN([DUNE_TECTONIC_CHECK_MODULE],
-[
-  DUNE_CHECK_MODULES([dune-tectonic],[tectonic/tectonic.hh])
-])

src/.gitignore

-timestepping.cc
-timestepping.hh
-timestepping_tmpl.cc
-one-body-sample.cc
\ No newline at end of file

src/CMakeLists.txt

+set(SW_SOURCE_FILES
+  assemblers.cc
+  enumparser.cc
+  hdf5/frictionalboundary-writer.cc
+  hdf5/iteration-writer.cc
+  hdf5/patchinfo-writer.cc
+  hdf5/restart-io.cc
+  hdf5/surface-writer.cc
+  hdf5/time-writer.cc
+  one-body-problem-data/mygeometry.cc
+  one-body-problem-data/mygrid.cc
+  one-body-problem.cc
+  spatial-solving/fixedpointiterator.cc
+  spatial-solving/solverfactory.cc
+  time-stepping/adaptivetimestepper.cc
+  time-stepping/coupledtimestepper.cc
+  time-stepping/rate.cc
+  time-stepping/rate/rateupdater.cc
+  time-stepping/state.cc
+  vtk.cc
+)
+set(UGW_SOURCE_FILES
+  assemblers.cc # FIXME
+  one-body-problem-data/mygrid.cc
+  uniform-grid-writer.cc
+  vtk.cc
+)
+
+foreach(_dim 2 3)
+  set(_sw_target one-body-problem-${_dim}D)
+  set(_ugw_target uniform-grid-writer-${_dim}D)
+
+  add_executable(${_sw_target} ${SW_SOURCE_FILES})
+  add_executable(${_ugw_target} ${UGW_SOURCE_FILES})
+  add_dune_ug_flags(${_sw_target})
+  add_dune_ug_flags(${_ugw_target})
+
+  add_dune_hdf5_flags(${_sw_target})
+
+  set_property(TARGET ${_sw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
+  set_property(TARGET ${_ugw_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
+endforeach()

src/Makefile.am

-SUBDIRS = tests
-
-bin_PROGRAMS = \
-	one-body-sample-2D \
-	one-body-sample-3D
-
-SOURCES = \
-	assemblers.cc \
-	compute_state_dieterich_euler.cc \
-	compute_state_dieterich_common.cc \
-	compute_state_ruina.cc \
-	mysolver.cc \
-	one-body-sample.cc \
-	timestepping.cc \
-	vtk.cc
-
-## 2D
-.PHONY: run-one-body-sample-2D
-run-one-body-sample-2D: one-body-sample-2D
-	libtool --mode execute ./one-body-sample-2D
-
-.PHONY: run-one-body-sample-2D-time
-run-one-body-sample-2D-time: one-body-sample-2D
-	libtool --mode execute time ./one-body-sample-2D >/dev/null
-
-.PHONY: run-one-body-sample-2D-valgrind
-run-one-body-sample-2D-valgrind: one-body-sample-2D
-	libtool --mode execute valgrind ./one-body-sample-2D
-
-.PHONY: run-one-body-sample-2D-gdb
-run-one-body-sample-2D-gdb: one-body-sample-2D
-	libtool --mode execute gdb ./one-body-sample-2D
-
-one_body_sample_2D_SOURCES = $(SOURCES)
-
-one_body_sample_2D_CPPFLAGS = \
-	$(AM_CPPFLAGS) -Dsrcdir=\"$(srcdir)\" -DDIM=2
-
-## 3D
-.PHONY: run-one-body-sample-3D
-run-one-body-sample-3D: one-body-sample-3D
-	libtool --mode execute ./one-body-sample-3D
-
-.PHONY: run-one-body-sample-3D-time
-run-one-body-sample-3D-time: one-body-sample-3D
-	libtool --mode execute time ./one-body-sample-3D >/dev/null
-
-.PHONY: run-one-body-sample-3D-valgrind
-run-one-body-sample-3D-valgrind: one-body-sample-3D
-	libtool --mode execute valgrind ./one-body-sample-3D
-
-.PHONY: run-one-body-sample-3D-gdb
-run-one-body-sample-3D-gdb: one-body-sample-3D
-	libtool --mode execute gdb ./one-body-sample-3D
-
-one_body_sample_3D_SOURCES = $(SOURCES)
-
-one_body_sample_3D_CPPFLAGS = \
-	$(AM_CPPFLAGS) -Dsrcdir=\"$(srcdir)\" -DDIM=3
-
-# Some are for clang, others are for gcc
-AM_CXXFLAGS = \
-	-Wall \
-	-Wextra \
-	-Wno-c++11-compat \
-	-Wno-c++11-extensions \
-	-Wno-deprecated-declarations \
-	-Wno-mismatched-tags \
-	-Wno-missing-declarations \
-	-Wno-overloaded-virtual \
-	-Wno-reorder \
-	-Wno-sign-compare \
-	-Wno-tautological-compare \
-	-Wno-type-limits \
-	-Wno-unused-parameter \
-	-Wno-unused-variable
-
-AM_CPPFLAGS = \
-	$(DUNE_CPPFLAGS) \
-	$(PYTHON_CPPFLAGS) \
-	$(ALUGRID_CPPFLAGS) \
-	-I$(top_srcdir)
-
-# The libraries have to be given in reverse order (most basic libraries
-# last).  Also, due to some misunderstanding, a lot of libraries include the
-# -L option in LDFLAGS instead of LIBS -- so we have to include the LDFLAGS
-# here as well.
-LDADD = \
-	$(DUNE_LDFLAGS) $(DUNE_LIBS) \
-	$(ALUGRID_LIBS) \
-	$(PYTHON_LIBS)
-AM_LDFLAGS = \
-	$(DUNE_LDFLAGS) \
-	$(ALUGRID_LDFLAGS) \
-	$(PYTHON_LDFLAGS)
-
-# pass most important options when "make distcheck" is used
-DISTCHECK_CONFIGURE_FLAGS = \
-	--with-dune-common=$(DUNE_COMMON_ROOT) \
-	--with-dune-fufem=$(DUNE_FUFEM_ROOT) \
-	--with-dune-tnnmg=$(DUNE_TNNMG_ROOT) \
-	CXX="$(CXX)" CC="$(CC)"
-
-include $(top_srcdir)/am/global-rules
-
-BUILT_SOURCES = timestepping.hh timestepping.cc
-# Make sure the two are not built in parallel
-$(srcdir)/timestepping.cc: $(srcdir)/timestepping.hh
-
-$(srcdir)/timestepping.hh: $(srcdir)/timestepping.org
-	emacs -Q --batch --eval \
-	  "(let (vc-handled-backends) \
-	      (org-babel-tangle-file \"$<\" nil 'c++))"
-
-$(srcdir)/one-body-sample.cc: $(srcdir)/one-body-sample.org
-	emacs -Q --batch --eval \
-	  "(let (vc-handled-backends) \
-	      (org-babel-tangle-file \"$<\" nil 'c++))"

src/TODO.org

-* consider truncating only in tangential direction

src/assemblers.cc

@@ -2,84 +2,181 @@
 #include "config.h"
 #endif
 
+#include <dune/istl/scaledidmatrix.hh>
+
+#include <dune/fufem/assemblers/localassemblers/boundarymassassembler.hh>
+#include <dune/fufem/assemblers/localassemblers/l2functionalassembler.hh>
+#include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
+#include <dune/fufem/assemblers/localassemblers/normalstressboundaryassembler.hh>
+#include <dune/fufem/assemblers/localassemblers/stvenantkirchhoffassembler.hh>
+#include <dune/fufem/assemblers/localassemblers/variablecoefficientviscosityassembler.hh>
+#include <dune/fufem/assemblers/localassemblers/vonmisesstressassembler.hh>
+#include <dune/fufem/assemblers/localassemblers/weightedmassassembler.hh>
 #include <dune/fufem/boundarypatch.hh>
+#include <dune/fufem/functions/basisgridfunction.hh>
 #include <dune/fufem/functions/constantfunction.hh>
-#include <dune/fufem/assemblers/localassemblers/neumannboundaryassembler.hh>
-#include <dune/fufem/assemblers/boundaryfunctionalassembler.hh>
+#include <dune/fufem/functiontools/p0p1interpolation.hh>
+#include <dune/fufem/quadraturerules/quadraturerulecache.hh>
 
-#include <dune/tectonic/globalruinanonlinearity.hh>
+#include <dune/tectonic/frictionpotential.hh>
+#include <dune/tectonic/globalratestatefriction.hh>
 
 #include "assemblers.hh"
 
-#include "enums.hh"
-#include "enum_parser.cc"
-
-// Assembles Neumann boundary term in f
-template <class GridType, class GridView, class LocalVectorType, class FEBasis>
-void assemble_neumann(GridView const &gridView, FEBasis const &feBasis,
-                      Dune::BitSetVector<1> const &neumannNodes,
-                      Dune::BlockVector<LocalVectorType> &f,
-                      Dune::VirtualFunction<double, double> const &neumann,
-                      double time) { // constant sample function on neumann
-                                     // boundary
-  BoundaryPatch<GridView> const neumannBoundary(gridView, neumannNodes);
-  LocalVectorType SampleVector(0);
-  neumann.evaluate(time, SampleVector[0]);
-  SampleVector[1] = 0;
-  ConstantFunction<LocalVectorType, LocalVectorType> const fNeumann(
-      SampleVector);
-  NeumannBoundaryAssembler<GridType, LocalVectorType> neumannBoundaryAssembler(
-      fNeumann);
-  BoundaryFunctionalAssembler<FEBasis>(feBasis, neumannBoundary)
-      .assemble(neumannBoundaryAssembler, f, true);
+template <class GridView, int dimension>
+MyAssembler<GridView, dimension>::MyAssembler(GridView const &_gridView)
+    : cellBasis(_gridView),
+      vertexBasis(_gridView),
+      gridView(_gridView),
+      cellAssembler(cellBasis, cellBasis),
+      vertexAssembler(vertexBasis, vertexBasis) {}
+
+template <class GridView, int dimension>
+void MyAssembler<GridView, dimension>::assembleFrictionalBoundaryMass(
+    BoundaryPatch<GridView> const &frictionalBoundary,
+    ScalarMatrix &frictionalBoundaryMass) const {
+  BoundaryMassAssembler<Grid, BoundaryPatch<GridView>, LocalVertexBasis,
+                        LocalVertexBasis, Dune::FieldMatrix<double, 1, 1>> const
+      frictionalBoundaryMassAssembler(frictionalBoundary);
+  vertexAssembler.assembleOperator(frictionalBoundaryMassAssembler,
+                                   frictionalBoundaryMass);
+}
+
+template <class GridView, int dimension>
+void MyAssembler<GridView, dimension>::assembleMass(
+    Dune::VirtualFunction<LocalVector, LocalScalarVector> const &
+        densityFunction,
+    Matrix &M) const {
+  // NOTE: We treat the weight as a constant function
+  QuadratureRuleKey quadKey(dimension, 0);
+
+  WeightedMassAssembler<Grid, LocalVertexBasis, LocalVertexBasis,
+                        Dune::VirtualFunction<LocalVector, LocalScalarVector>,
+                        Dune::ScaledIdentityMatrix<double, dimension>>
+      localWeightedMass(gridView.grid(), densityFunction, quadKey);
+  vertexAssembler.assembleOperator(localWeightedMass, M);
 }
 
-// Assembles constant 1-function on frictional boundary in nodalIntegrals
-template <class GridType, class GridView, class LocalVectorType, class FEBasis>
-Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>>
-assemble_frictional(GridView const &gridView, FEBasis const &feBasis,
-                    Dune::BitSetVector<1> const &frictionalNodes) {
-  typedef Dune::FieldVector<double, 1> Singleton;
-  BoundaryPatch<GridView> const frictionalBoundary(gridView, frictionalNodes);
-  ConstantFunction<LocalVectorType, Singleton> const constantOneFunction(1);
-  NeumannBoundaryAssembler<GridType, Singleton> frictionalBoundaryAssembler(
-      constantOneFunction);
-
-  auto const nodalIntegrals = Dune::make_shared<Dune::BlockVector<Singleton>>();
-  BoundaryFunctionalAssembler<FEBasis>(feBasis, frictionalBoundary)
-      .assemble(frictionalBoundaryAssembler, *nodalIntegrals, true);
-  return nodalIntegrals;
+template <class GridView, int dimension>
+void MyAssembler<GridView, dimension>::assembleElasticity(double E, double nu,
+                                                          Matrix &A) const {
+  StVenantKirchhoffAssembler<Grid, LocalVertexBasis, LocalVertexBasis> const
+      localStiffness(E, nu);
+  vertexAssembler.assembleOperator(localStiffness, A);
 }
 
-template <class MatrixType, class VectorType>
-Dune::shared_ptr<Dune::GlobalNonlinearity<MatrixType, VectorType> const>
-assemble_nonlinearity(
-    Dune::ParameterTree const &parset,
-    Dune::BlockVector<Dune::FieldVector<double, 1>> const &nodalIntegrals,
-    Dune::BlockVector<Dune::FieldVector<double, 1>> const &state,
-    Dune::BlockVector<Dune::FieldVector<double, 1>> const &normalStress) {
-  auto const size = nodalIntegrals.size();
+template <class GridView, int dimension>
+void MyAssembler<GridView, dimension>::assembleViscosity(
+    Dune::VirtualFunction<LocalVector, LocalScalarVector> const &shearViscosity,
+    Dune::VirtualFunction<LocalVector, LocalScalarVector> const &bulkViscosity,
+    Matrix &C) const {
+  // NOTE: We treat the weights as constant functions
+  QuadratureRuleKey shearViscosityKey(dimension, 0);
+  QuadratureRuleKey bulkViscosityKey(dimension, 0);
+  VariableCoefficientViscosityAssembler<
+      Grid, LocalVertexBasis, LocalVertexBasis,
+      Dune::VirtualFunction<LocalVector, LocalScalarVector>> const
+      localViscosity(gridView.grid(), shearViscosity, bulkViscosity,
+                     shearViscosityKey, bulkViscosityKey);
+  vertexAssembler.assembleOperator(localViscosity, C);
+}
+
+template <class GridView, int dimension>
+void MyAssembler<GridView, dimension>::assembleBodyForce(
+    Dune::VirtualFunction<LocalVector, LocalVector> const &gravityField,
+    Vector &f) const {
+  L2FunctionalAssembler<Grid, LocalVertexBasis, LocalVector>
+      gravityFunctionalAssembler(gravityField);
+  vertexAssembler.assembleFunctional(gravityFunctionalAssembler, f);
+}
 
-  typedef Dune::BlockVector<Dune::FieldVector<double, 1>> SingletonVectorType;
-  // {{{ Assemble terms for the nonlinearity
-  SingletonVectorType mu(size);
-  mu = parset.get<double>("mu");
+template <class GridView, int dimension>
+void MyAssembler<GridView, dimension>::assembleNeumann(
+    BoundaryPatch<GridView> const &neumannBoundary, Vector &f,
+    Dune::VirtualFunction<double, double> const &neumann,
+    double relativeTime) const {
+  LocalVector localNeumann(0);
+  neumann.evaluate(relativeTime, localNeumann[0]);
+  NeumannBoundaryAssembler<Grid, LocalVector> neumannBoundaryAssembler(
+      std::make_shared<ConstantFunction<LocalVector, LocalVector>>(
+          localNeumann));
+  vertexAssembler.assembleBoundaryFunctional(neumannBoundaryAssembler, f,
+                                             neumannBoundary);
+}
 
-  SingletonVectorType a(size);
-  a = parset.get<double>("a");
+template <class GridView, int dimension>
+void MyAssembler<GridView, dimension>::assembleWeightedNormalStress(
+    BoundaryPatch<GridView> const &frictionalBoundary,
+    ScalarVector &weightedNormalStress, double youngModulus,
+    double poissonRatio, Vector const &displacement) const {
 
-  SingletonVectorType V0(size);
-  V0 = parset.get<double>("V0");
+  BasisGridFunction<VertexBasis, Vector> displacementFunction(vertexBasis,
+                                                              displacement);
+  Vector traction(cellBasis.size());
+  NormalStressBoundaryAssembler<Grid> tractionBoundaryAssembler(
+      youngModulus, poissonRatio, &displacementFunction, 1);
+  cellAssembler.assembleBoundaryFunctional(tractionBoundaryAssembler, traction,
+                                           frictionalBoundary);
 
-  SingletonVectorType b(size);
-  b = parset.get<double>("b");
+  auto const nodalTractionAverage =
+      interpolateP0ToP1(frictionalBoundary, traction);
 
-  SingletonVectorType L(size);
-  L = parset.get<double>("L");
+  ScalarVector weights;
+  {
+    NeumannBoundaryAssembler<Grid, typename ScalarVector::block_type>
+        frictionalBoundaryAssembler(
+            std::make_shared<ConstantFunction<
+                LocalVector, typename ScalarVector::block_type>>(1));
+    vertexAssembler.assembleBoundaryFunctional(frictionalBoundaryAssembler,
+                                               weights, frictionalBoundary);
+  }
+  auto const normals = frictionalBoundary.getNormals();
+  for (size_t i = 0; i < vertexBasis.size(); ++i)
+    weightedNormalStress[i] =
+        std::fmin(normals[i] * nodalTractionAverage[i], 0) * weights[i];
+}
+
+template <class GridView, int dimension>
+auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
+    Config::FrictionModel frictionModel,
+    BoundaryPatch<GridView> const &frictionalBoundary,
+    GlobalFrictionData<dimension> const &frictionInfo,
+    ScalarVector const &weightedNormalStress) const
+    -> std::shared_ptr<GlobalFriction<Matrix, Vector>> {
+  // Lumping of the nonlinearity
+  ScalarVector weights;
+  {
+    NeumannBoundaryAssembler<Grid, typename ScalarVector::block_type>
+        frictionalBoundaryAssembler(std::make_shared<
+            ConstantFunction<LocalVector, typename ScalarVector::block_type>>(
+            1));
+    vertexAssembler.assembleBoundaryFunctional(frictionalBoundaryAssembler,
+                                               weights, frictionalBoundary);
+  }
+  switch (frictionModel) {
+    case Config::Truncated:
+      return std::make_shared<GlobalRateStateFriction<
+          Matrix, Vector, TruncatedRateState, GridView>>(
+          frictionalBoundary, frictionInfo, weights, weightedNormalStress);
+    case Config::Regularised:
+      return std::make_shared<GlobalRateStateFriction<
+          Matrix, Vector, RegularisedRateState, GridView>>(
+          frictionalBoundary, frictionInfo, weights, weightedNormalStress);
+    default:
+      assert(false);
+  }
+}
 
-  return Dune::make_shared<
-      Dune::GlobalRuinaNonlinearity<MatrixType, VectorType> const>(
-      nodalIntegrals, a, mu, V0, normalStress, b, state, L);
+template <class GridView, int dimension>
+void MyAssembler<GridView, dimension>::assembleVonMisesStress(
+    double youngModulus, double poissonRatio, Vector const &u,
+    ScalarVector &stress) const {
+  auto const gridDisplacement =
+      std::make_shared<BasisGridFunction<VertexBasis, Vector> const>(
+          vertexBasis, u);
+  VonMisesStressAssembler<Grid, LocalCellBasis> localStressAssembler(
+      youngModulus, poissonRatio, gridDisplacement);
+  cellAssembler.assembleFunctional(localStressAssembler, stress);
 }
 
 #include "assemblers_tmpl.cc"

src/assemblers.hh

-#ifndef ASSEMBLERS_HH
-#define ASSEMBLERS_HH
+#ifndef SRC_ASSEMBLERS_HH
+#define SRC_ASSEMBLERS_HH
 
 #include <dune/common/bitsetvector.hh>
 #include <dune/common/function.hh>
-#include <dune/common/fvector.hh>
-#include <dune/common/parametertree.hh>
-#include <dune/common/shared_ptr.hh>
+#include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/bvector.hh>
 
-#include <dune/tectonic/globalnonlinearity.hh>
-
-template <class GridType, class GridView, class LocalVectorType, class FEBasis>
-void assemble_neumann(GridView const &gridView, FEBasis const &feBasis,
-                      Dune::BitSetVector<1> const &neumannNodes,
-                      Dune::BlockVector<LocalVectorType> &f,
-                      Dune::VirtualFunction<double, double> const &neumann,
-                      double time);
-
-template <class GridType, class GridView, class LocalVectorType, class FEBasis>
-Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>>
-assemble_frictional(GridView const &gridView, FEBasis const &feBasis,
-                    Dune::BitSetVector<1> const &frictionalNodes);
-
-template <class MatrixType, class VectorType>
-Dune::shared_ptr<Dune::GlobalNonlinearity<MatrixType, VectorType> const>
-assemble_nonlinearity(
-    Dune::ParameterTree const &parset,
-    Dune::BlockVector<Dune::FieldVector<double, 1>> const &nodalIntegrals,
-    Dune::BlockVector<Dune::FieldVector<double, 1>> const &state,
-    Dune::BlockVector<Dune::FieldVector<double, 1>> const &normalStress);
+#include <dune/fufem/assemblers/assembler.hh>
+#pragma clang diagnostic push
+#pragma clang diagnostic ignored "-Wsign-compare"
+#include <dune/fufem/functionspacebases/p0basis.hh>
+#pragma clang diagnostic pop
+#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
+
+#include <dune/tectonic/globalfriction.hh>
+#include <dune/tectonic/globalfrictiondata.hh>
+
+#include "enums.hh"
+
+template <class GridView, int dimension> class MyAssembler {
+public:
+  using ScalarMatrix = Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>>;
+  using Matrix =
+      Dune::BCRSMatrix<Dune::FieldMatrix<double, dimension, dimension>>;
+  using ScalarVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
+  using Vector = Dune::BlockVector<Dune::FieldVector<double, dimension>>;
+
+  using CellBasis = P0Basis<GridView, double>;
+  using VertexBasis = P1NodalBasis<GridView, double>;
+
+  CellBasis const cellBasis;
+  VertexBasis const vertexBasis;
+  GridView const &gridView;
+
+private:
+  using Grid = typename GridView::Grid;
+  using LocalVector = typename Vector::block_type;
+  using LocalScalarVector = typename ScalarVector::block_type;
+
+  using LocalCellBasis = typename CellBasis::LocalFiniteElement;
+  using LocalVertexBasis = typename VertexBasis::LocalFiniteElement;
+
+  Assembler<CellBasis, CellBasis> cellAssembler;
+  Assembler<VertexBasis, VertexBasis> vertexAssembler;
+
+public:
+  MyAssembler(GridView const &gridView);
+
+  void assembleFrictionalBoundaryMass(
+      BoundaryPatch<GridView> const &frictionalBoundary,
+      ScalarMatrix &frictionalBoundaryMass) const;
+
+  void assembleMass(Dune::VirtualFunction<
+                        LocalVector, LocalScalarVector> const &densityFunction,
+                    Matrix &M) const;
+
+  void assembleElasticity(double E, double nu, Matrix &A) const;
+
+  void assembleViscosity(
+      Dune::VirtualFunction<LocalVector, LocalScalarVector> const &
+          shearViscosity,
+      Dune::VirtualFunction<LocalVector, LocalScalarVector> const &
+          bulkViscosity,
+      Matrix &C) const;
+
+  void assembleBodyForce(
+      Dune::VirtualFunction<LocalVector, LocalVector> const &gravityField,
+      Vector &f) const;
+
+  void assembleNeumann(BoundaryPatch<GridView> const &neumannBoundary,
+                       Vector &f,
+                       Dune::VirtualFunction<double, double> const &neumann,
+                       double relativeTime) const;
+
+  void assembleWeightedNormalStress(
+      BoundaryPatch<GridView> const &frictionalBoundary,
+      ScalarVector &weightedNormalStress, double youngModulus,
+      double poissonRatio, Vector const &displacement) const;
+
+  std::shared_ptr<GlobalFriction<Matrix, Vector>> assembleFrictionNonlinearity(
+      Config::FrictionModel frictionModel,
+      BoundaryPatch<GridView> const &frictionalBoundary,
+      GlobalFrictionData<dimension> const &frictionInfo,
+      ScalarVector const &weightedNormalStress) const;
+
+  void assembleVonMisesStress(double youngModulus, double poissonRatio,
+                              Vector const &u, ScalarVector &stress) const;
+};
 #endif

src/assemblers_tmpl.cc

-#ifndef DIM
-#error DIM unset
+#ifndef MY_DIM
+#error MY_DIM unset
 #endif
 
-#include <dune/common/fmatrix.hh>
-#include <dune/common/fvector.hh>
-#include <dune/grid/alugrid.hh>
-#include <dune/istl/bvector.hh>
-#include <dune/istl/bcrsmatrix.hh>
+#include "explicitgrid.hh"
 
-#include <dune/fufem/functionspacebases/p1nodalbasis.hh>
-
-typedef Dune::FieldVector<double, DIM> SmallVector;
-typedef Dune::FieldMatrix<double, DIM, DIM> SmallMatrix;
-typedef Dune::BCRSMatrix<SmallMatrix> MatrixType;
-typedef Dune::BlockVector<SmallVector> VectorType;
-
-typedef Dune::ALUGrid<DIM, DIM, Dune::simplex, Dune::nonconforming> GridType;
-typedef GridType::LeafGridView GridView;
-typedef P1NodalBasis<GridView, double> P1Basis;
-
-template void assemble_neumann<GridType, GridView, SmallVector, P1Basis>(
-    GridView const &gridView, P1Basis const &feBasis,
-    Dune::BitSetVector<1> const &neumannNodes,
-    Dune::BlockVector<SmallVector> &f,
-    Dune::VirtualFunction<double, double> const &neumann, double time);
-template Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>>
-assemble_frictional<GridType, GridView, SmallVector, P1Basis>(
-    GridView const &gridView, P1Basis const &feBasis,
-    Dune::BitSetVector<1> const &frictionalNodes);
-
-template Dune::shared_ptr<
-    Dune::GlobalNonlinearity<MatrixType, VectorType> const>
-assemble_nonlinearity<MatrixType, VectorType>(
-    Dune::ParameterTree const &parset,
-    Dune::BlockVector<Dune::FieldVector<double, 1>> const &nodalIntegrals,
-    Dune::BlockVector<Dune::FieldVector<double, 1>> const &state,
-    Dune::BlockVector<Dune::FieldVector<double, 1>> const &normalStress);
+template class MyAssembler<GridView, MY_DIM>;

src/commands.m

-## State evolution
-load ../data;
-graphics_toolkit fltk;
-path(path, '../../../');
-plot_state
-
-## Comparison of steady-state convergence rates, to and fro
-# 600 time steps
-load ../data;
-graphics_toolkit fltk;
-path(path, '../../../');
-plot(1:144,A(121:(120+144),1)-A(121,1),1:144,-interp1(1:240, A(361:600,1), linspace(1,240,144))+A(361,1));
-axis tight;
-legend('decline (first)','incline (second)');
-
-# 6000 time steps
-load ../data;
-graphics_toolkit fltk;
-path(path, '../../../');
-plot(1:1440,A(1201:(1200+1440),1)-A(1201,1),1:1440,-interp1(1:2400, A(3601:6000,1), linspace(1,2400,1440))+A(3601,1));
-axis tight;
-legend('decline (first)','incline (second)');

src/compute_state_dieterich_common.cc

-#include "compute_state_dieterich_common.hh"
-
-DieterichNonlinearity::DieterichNonlinearity(double _VoL) : VoL(_VoL) {}
-
-void DieterichNonlinearity::directionalSubDiff(VectorType const &u,
-                                               VectorType const &v,
-                                               Interval<double> &D) const {
-  D[0] = D[1] = v[0] * (VoL - std::exp(-u[0]));
-}
-
-void DieterichNonlinearity::directionalDomain(VectorType const &,
-                                              VectorType const &,
-                                              Interval<double> &dom) const {
-  dom[0] = -std::numeric_limits<double>::max();
-  dom[1] = std::numeric_limits<double>::max();
-}

src/compute_state_dieterich_common.hh

-#include <dune/common/fvector.hh>
-#include <dune/common/fmatrix.hh>
-#include <dune/fufem/interval.hh>
-
-class DieterichNonlinearity {
-public:
-  typedef Dune::FieldVector<double, 1> VectorType;
-  typedef Dune::FieldMatrix<double, 1, 1> MatrixType;
-
-  DieterichNonlinearity(double _VoL);
-  void directionalSubDiff(VectorType const &u, VectorType const &v,
-                          Interval<double> &D) const;
-  void directionalDomain(VectorType const &, VectorType const &,
-                         Interval<double> &dom) const;
-
-private:
-  double const VoL;
-};

src/compute_state_dieterich_euler.cc

-#include <cassert>
-
-#include <dune/tnnmg/problem-classes/bisection.hh>
-
-#include <dune/tectonic/mydirectionalconvexfunction.hh>
-
-#include "compute_state_dieterich_common.hh"
-#include "compute_state_dieterich_euler.hh"
-
-double state_update_dieterich_euler_bisection(double tau, double VoL,
-                                              double old_state) {
-  DieterichNonlinearity::VectorType const start(0);
-  DieterichNonlinearity::VectorType const direction(1);
-  DieterichNonlinearity const phi(VoL);
-  MyDirectionalConvexFunction<DieterichNonlinearity> const J(
-      1.0 / tau, old_state / tau, phi, start, direction);
-
-  int bisectionsteps = 0;
-  Bisection const bisection(0.0, 1.0, 1e-12, true, 0);    // TODO
-  return bisection.minimize(J, 0.0, 0.0, bisectionsteps); // TODO
-}
-
-double state_update_dieterich_euler(double tau, double VoL, double old_state) {
-  double const ret =
-      state_update_dieterich_euler_bisection(tau, VoL, old_state);
-  /* We have
-
-     ret - old_state + tau * VoL = tau*exp(-ret)
-
-     or
-
-     log((ret - old_state)/tau + VoL) = -ret
-  */
-  assert(std::min(std::abs(ret - old_state + tau * (VoL - std::exp(-ret))),
-                  std::abs(std::log((ret - old_state) / tau + VoL) + ret)) <
-         1e-12);
-  return ret;
-}

src/compute_state_dieterich_euler.hh

-#ifndef COMPUTE_STATE_DIETERICH_EULER_HH
-#define COMPUTE_STATE_DIETERICH_EULER_HH
-
-double state_update_dieterich_euler(double h, double VoL, double old_state);
-#endif

src/compute_state_ruina.cc

-#include <cassert>
-
-#include <cmath>
-
-#include "compute_state_ruina.hh"
-
-double state_update_ruina(double tau, double uol, double old_state) {
-  if (uol == 0)
-    return old_state;
-
-  double const ret = old_state - uol * std::log(uol / tau);
-  return ret / (1 + uol);
-}

src/compute_state_ruina.hh

-#ifndef COMPUTE_STATE_RUINA_HH
-#define COMPUTE_STATE_RUINA_HH
-
-double state_update_ruina(double h, double uol, double old_state);
-#endif

src/diameter.hh

+#ifndef SRC_DIAMETER_HH
+#define SRC_DIAMETER_HH
+
+template <class Geometry> double diameter(Geometry const &geometry) {
+  auto const numCorners = geometry.corners();
+  std::vector<typename Geometry::GlobalCoordinate> corners(numCorners);
+
+  double diameter = 0.0;
+  for (int i = 0; i < numCorners; ++i) {
+    corners[i] = geometry.corner(i);
+    for (int j = 0; j < i; ++j)
+      diameter = std::max(diameter, (corners[i] - corners[j]).two_norm());
+  }
+  return diameter;
+}
+#endif

src/dieterichstateupdater.hh

-#ifndef DIETERICH_STATE_UPDATER_HH
-#define DIETERICH_STATE_UPDATER_HH
-
-#include "compute_state_dieterich_euler.hh"
-#include "stateupdater.hh"
-
-template <class SingletonVectorType, class VectorType>
-class DieterichStateUpdater
-    : public StateUpdater<SingletonVectorType, VectorType> {
-public:
-  DieterichStateUpdater(SingletonVectorType _alpha_initial,
-                        Dune::BitSetVector<1> const &_nodes, double _L);
-
-  virtual void nextTimeStep();
-  virtual void setup(double _tau);
-  virtual void solve(VectorType const &velocity_field);
-  virtual void extractState(SingletonVectorType &state);
-
-private:
-  SingletonVectorType alpha_old;
-  SingletonVectorType alpha;
-  Dune::BitSetVector<1> const &nodes;
-  double const L;
-  double tau;
-};
-
-template <class SingletonVectorType, class VectorType>
-DieterichStateUpdater<SingletonVectorType, VectorType>::DieterichStateUpdater(
-    SingletonVectorType _alpha_initial, Dune::BitSetVector<1> const &_nodes,
-    double _L)
-    : alpha(_alpha_initial), nodes(_nodes), L(_L) {}
-
-template <class SingletonVectorType, class VectorType>
-void DieterichStateUpdater<SingletonVectorType, VectorType>::nextTimeStep() {
-  alpha_old = alpha;
-}
-
-template <class SingletonVectorType, class VectorType>
-void DieterichStateUpdater<SingletonVectorType, VectorType>::setup(
-    double _tau) {
-  tau = _tau;
-}
-
-template <class SingletonVectorType, class VectorType>
-void DieterichStateUpdater<SingletonVectorType, VectorType>::solve(
-    VectorType const &velocity_field) {
-  for (size_t i = 0; i < nodes.size(); ++i)
-    if (nodes[i][0]) {
-      double const V = velocity_field[i].two_norm();
-      alpha[i] = state_update_dieterich_euler(tau, V / L, alpha_old[i]);
-    }
-}
-
-template <class SingletonVectorType, class VectorType>
-void DieterichStateUpdater<SingletonVectorType, VectorType>::extractState(
-    SingletonVectorType &state) {
-  state = alpha;
-}
-
-#endif