diff --git a/dune/tectonic/globalfriction.hh b/dune/tectonic/globalfriction.hh
index db9cec0eabb9a89015c49a11db8319efbfefe7e9..3fb7c39999db0cc34b21e08e1034c82918156faf 100644
--- a/dune/tectonic/globalfriction.hh
+++ b/dune/tectonic/globalfriction.hh
@@ -82,7 +82,7 @@ template <class Matrix, class Vector> class GlobalFriction {
     return res->regularity(x);
   }
 
-  void coefficientOfFriction(Vector const &x, ScalarVector &coefficient) {
+  void coefficientOfFriction(Vector const &x, ScalarVector &coefficient) const {
     coefficient.resize(x.size());
     for (size_t i = 0; i < x.size(); ++i)
       coefficient[i] = restriction(i)->coefficientOfFriction(x[i]);
diff --git a/src/assemblers.cc b/src/assemblers.cc
index 64a999043f1ffce5151bc87683a91889e7941b6d..29ceb758f9c7eecda7539426aee415102e9113db 100644
--- a/src/assemblers.cc
+++ b/src/assemblers.cc
@@ -33,7 +33,7 @@ MyAssembler<GridView, dimension>::MyAssembler(GridView const &_gridView)
 template <class GridView, int dimension>
 void MyAssembler<GridView, dimension>::assembleFrictionalBoundaryMass(
     BoundaryPatch<GridView> const &frictionalBoundary,
-    ScalarMatrix &frictionalBoundaryMass) {
+    ScalarMatrix &frictionalBoundaryMass) const {
   BoundaryMassAssembler<Grid, BoundaryPatch<GridView>, LocalVertexBasis,
                         LocalVertexBasis, Dune::FieldMatrix<double, 1, 1>> const
       frictionalBoundaryMassAssembler(frictionalBoundary);
@@ -45,7 +45,7 @@ template <class GridView, int dimension>
 void MyAssembler<GridView, dimension>::assembleMass(
     Dune::VirtualFunction<LocalVector, LocalScalarVector> const &
         densityFunction,
-    Matrix &M) {
+    Matrix &M) const {
   // NOTE: We treat the weight as a constant function
   QuadratureRuleKey quadKey(dimension, 0);
 
@@ -58,7 +58,7 @@ void MyAssembler<GridView, dimension>::assembleMass(
 
 template <class GridView, int dimension>
 void MyAssembler<GridView, dimension>::assembleElasticity(double E, double nu,
-                                                          Matrix &A) {
+                                                          Matrix &A) const {
   StVenantKirchhoffAssembler<Grid, LocalVertexBasis, LocalVertexBasis> const
       localStiffness(E, nu);
   vertexAssembler.assembleOperator(localStiffness, A);
@@ -68,7 +68,7 @@ template <class GridView, int dimension>
 void MyAssembler<GridView, dimension>::assembleViscosity(
     Dune::VirtualFunction<LocalVector, LocalScalarVector> const &shearViscosity,
     Dune::VirtualFunction<LocalVector, LocalScalarVector> const &bulkViscosity,
-    Matrix &C) {
+    Matrix &C) const {
   // NOTE: We treat the weights as constant functions
   QuadratureRuleKey shearViscosityKey(dimension, 0);
   QuadratureRuleKey bulkViscosityKey(dimension, 0);
@@ -83,7 +83,7 @@ void MyAssembler<GridView, dimension>::assembleViscosity(
 template <class GridView, int dimension>
 void MyAssembler<GridView, dimension>::assembleBodyForce(
     Dune::VirtualFunction<LocalVector, LocalVector> const &gravityField,
-    Vector &f) {
+    Vector &f) const {
   L2FunctionalAssembler<Grid, LocalVertexBasis, LocalVector>
       gravityFunctionalAssembler(gravityField);
   vertexAssembler.assembleFunctional(gravityFunctionalAssembler, f);
@@ -92,7 +92,8 @@ void MyAssembler<GridView, dimension>::assembleBodyForce(
 template <class GridView, int dimension>
 void MyAssembler<GridView, dimension>::assembleNeumann(
     BoundaryPatch<GridView> const &neumannBoundary, Vector &f,
-    Dune::VirtualFunction<double, double> const &neumann, double relativeTime) {
+    Dune::VirtualFunction<double, double> const &neumann,
+    double relativeTime) const {
   LocalVector localNeumann(0);
   neumann.evaluate(relativeTime, localNeumann[0]);
   NeumannBoundaryAssembler<Grid, LocalVector> neumannBoundaryAssembler(
@@ -106,7 +107,7 @@ template <class GridView, int dimension>
 void MyAssembler<GridView, dimension>::assembleNormalStress(
     BoundaryPatch<GridView> const &frictionalBoundary,
     ScalarVector &normalStress, double youngModulus, double poissonRatio,
-    Vector const &displacement) {
+    Vector const &displacement) const {
 
   Vector traction;
   Stress<GridView>::getElasticSurfaceNormalStress // misnomer(!)
@@ -129,7 +130,7 @@ auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
     Config::FrictionModel frictionModel,
     BoundaryPatch<GridView> const &frictionalBoundary,
     GlobalFrictionData<dimension> const &frictionInfo,
-    ScalarVector const &normalStress)
+    ScalarVector const &normalStress) const
     -> std::shared_ptr<GlobalFriction<Matrix, Vector>> {
   // Lump negative normal stress (kludge)
   ScalarVector weights;
@@ -158,7 +159,7 @@ auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
 template <class GridView, int dimension>
 void MyAssembler<GridView, dimension>::assembleVonMisesStress(
     double youngModulus, double poissonRatio, Vector const &u,
-    ScalarVector &stress) {
+    ScalarVector &stress) const {
   auto const gridDisplacement =
       std::make_shared<BasisGridFunction<VertexBasis, Vector> const>(
           vertexBasis, u);
diff --git a/src/assemblers.hh b/src/assemblers.hh
index 207bbde220327e999d418f2906de7e09f467a99f..14abad33f9bee7587f0f1182ac5fb9452f657c0e 100644
--- a/src/assemblers.hh
+++ b/src/assemblers.hh
@@ -29,8 +29,8 @@ template <class GridView, int dimension> class MyAssembler {
   using CellBasis = P0Basis<GridView, double>;
   using VertexBasis = P1NodalBasis<GridView, double>;
 
-  CellBasis cellBasis;
-  VertexBasis vertexBasis;
+  CellBasis const cellBasis;
+  VertexBasis const vertexBasis;
 
 private:
   using Grid = typename GridView::Grid;
@@ -49,41 +49,42 @@ template <class GridView, int dimension> class MyAssembler {
 
   void assembleFrictionalBoundaryMass(
       BoundaryPatch<GridView> const &frictionalBoundary,
-      ScalarMatrix &frictionalBoundaryMass);
+      ScalarMatrix &frictionalBoundaryMass) const;
 
   void assembleMass(Dune::VirtualFunction<
                         LocalVector, LocalScalarVector> const &densityFunction,
-                    Matrix &M);
+                    Matrix &M) const;
 
-  void assembleElasticity(double E, double nu, Matrix &A);
+  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);
+      Matrix &C) const;
 
   void assembleBodyForce(
       Dune::VirtualFunction<LocalVector, LocalVector> const &gravityField,
-      Vector &f);
+      Vector &f) const;
 
   void assembleNeumann(BoundaryPatch<GridView> const &neumannBoundary,
                        Vector &f,
                        Dune::VirtualFunction<double, double> const &neumann,
-                       double relativeTime);
+                       double relativeTime) const;
 
   void assembleNormalStress(BoundaryPatch<GridView> const &frictionalBoundary,
                             ScalarVector &normalStress, double youngModulus,
-                            double poissonRatio, Vector const &displacement);
+                            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 &normalStress);
+      ScalarVector const &normalStress) const;
 
   void assembleVonMisesStress(double youngModulus, double poissonRatio,
-                              Vector const &u, ScalarVector &stress);
+                              Vector const &u, ScalarVector &stress) const;
 };
 #endif
diff --git a/src/sand-wedge.cc b/src/sand-wedge.cc
index 78ef4edeebd324738d9b67ec7594779a66414ad5..af0fe30b561556be77a833600d57787f28068c14 100644
--- a/src/sand-wedge.cc
+++ b/src/sand-wedge.cc
@@ -183,7 +183,7 @@ int main(int argc, char *argv[]) {
       timeStepWriter << _relativeTime << " " << _relativeTau << std::endl;
     };
 
-    MyAssembler myAssembler(leafView);
+    MyAssembler const myAssembler(leafView);
 
     MyBody<dims> const body(parset);