dune/tectonic/globalfriction.hh

@@ -22,13 +22,12 @@ template <class Matrix, class Vector> class GlobalFriction {
   using LocalMatrix = typename Matrix::block_type;
   using LocalVectorType = typename Vector::block_type;
   size_t static const block_size = LocalVectorType::dimension;
-  using Friction = LocalFriction<block_size>;
+  using LocalNonlinearity = LocalFriction<block_size>;
 
   double operator()(Vector const &x) const {
     double tmp = 0;
     for (size_t i = 0; i < x.size(); ++i) {
-      auto const res = restriction(i);
-      tmp += (*res)(x[i]);
+      tmp += restriction(i)(x[i]);
     }
     return tmp;
   }
@@ -36,14 +35,11 @@ template <class Matrix, class Vector> class GlobalFriction {
   /*
     Return a restriction of the outer function to the i'th node.
   */
-  std::shared_ptr<LocalFriction<block_size>> virtual restriction(size_t i)
-      const = 0;
+  LocalNonlinearity const virtual &restriction(size_t i) const = 0;
 
   void addHessian(Vector const &v, Matrix &hessian) const {
-    for (size_t i = 0; i < v.size(); ++i) {
-      auto const res = restriction(i);
-      res->addHessian(v[i], hessian[i][i]);
-    }
+    for (size_t i = 0; i < v.size(); ++i)
+      restriction(i).addHessian(v[i], hessian[i][i]);
   }
 
   void directionalDomain(Vector const &, Vector const &,
@@ -52,14 +48,13 @@ template <class Matrix, class Vector> class GlobalFriction {
     dom[1] = std::numeric_limits<double>::max();
   }
 
-  void directionalSubDiff(Vector const &u, Vector const &v,
-                          Dune::Solvers::Interval<double> &subdifferential)
-      const {
+  void directionalSubDiff(
+      Vector const &u, Vector const &v,
+      Dune::Solvers::Interval<double> &subdifferential) const {
     subdifferential[0] = subdifferential[1] = 0;
     for (size_t i = 0; i < u.size(); ++i) {
       Dune::Solvers::Interval<double> D;
-      auto const res = restriction(i);
-      res->directionalSubDiff(u[i], v[i], D);
+      restriction(i).directionalSubDiff(u[i], v[i], D);
       subdifferential[0] += D[0];
       subdifferential[1] += D[1];
     }
@@ -71,21 +66,19 @@ template <class Matrix, class Vector> class GlobalFriction {
   }
 
   void addGradient(Vector const &v, Vector &gradient) const {
-    for (size_t i = 0; i < v.size(); ++i) {
-      auto const res = restriction(i);
-      res->addGradient(v[i], gradient[i]);
-    }
+    for (size_t i = 0; i < v.size(); ++i)
+      restriction(i).addGradient(v[i], gradient[i]);
   }
 
   double regularity(size_t i, typename Vector::block_type const &x) const {
-    auto const res = restriction(i);
-    return res->regularity(x);
+    return restriction(i).regularity(x);
   }
 
-  void coefficientOfFriction(Vector const &x, ScalarVector &coefficient) {
-    coefficient.resize(x.size());
+  ScalarVector coefficientOfFriction(Vector const &x) const {
+    ScalarVector ret(x.size());
     for (size_t i = 0; i < x.size(); ++i)
-      coefficient[i] = restriction(i)->coefficientOfFriction(x[i]);
+      ret[i] = restriction(i).coefficientOfFriction(x[i]);
+    return ret;
   }
 
   void virtual updateAlpha(ScalarVector const &alpha) = 0;

dune/tectonic/globalfrictiondata.hh

@@ -23,8 +23,9 @@ template <int dimension> class GlobalFrictionData {
   }
 
 protected:
-  using VirtualFunction = Dune::VirtualFunction<
-      Dune::FieldVector<double, dimension>, Dune::FieldVector<double, 1>>;
+  using VirtualFunction =
+      Dune::VirtualFunction<Dune::FieldVector<double, dimension>,
+                            Dune::FieldVector<double, 1>>;
 
   double virtual const &C() const = 0;
   double virtual const &L() const = 0;

dune/tectonic/globalratestatefriction.hh

@@ -10,58 +10,62 @@
 #include <dune/istl/bcrsmatrix.hh>
 #include <dune/istl/bvector.hh>
 
+#include <dune/tectonic/geocoordinate.hh>
 #include <dune/tectonic/globalfrictiondata.hh>
 #include <dune/tectonic/globalfriction.hh>
+#include <dune/tectonic/index-in-sorted-range.hh>
 
 template <class Matrix, class Vector, class ScalarFriction, class GridView>
 class GlobalRateStateFriction : public GlobalFriction<Matrix, Vector> {
 public:
   using GlobalFriction<Matrix, Vector>::block_size;
-  using typename GlobalFriction<Matrix, Vector>::Friction;
+  using typename GlobalFriction<Matrix, Vector>::LocalNonlinearity;
 
 private:
   using typename GlobalFriction<Matrix, Vector>::ScalarVector;
 
 public:
   GlobalRateStateFriction(BoundaryPatch<GridView> const &frictionalBoundary,
-                          GridView const &gridView,
                           GlobalFrictionData<block_size> const &frictionInfo,
-                          // Note: passing the following two makes no sense
                           ScalarVector const &weights,
-                          ScalarVector const &normalStress)
-      : restrictions(normalStress.size()) {
-    auto trivialNonlinearity =
-        std::make_shared<Friction>(std::make_shared<TrivialFunction>());
-
+                          ScalarVector const &weightedNormalStress)
+      : restrictions(), localToGlobal(), zeroFriction() {
+    auto const gridView = frictionalBoundary.gridView();
     Dune::MultipleCodimMultipleGeomTypeMapper<
         GridView, Dune::MCMGVertexLayout> const vertexMapper(gridView);
     for (auto it = gridView.template begin<block_size>();
          it != gridView.template end<block_size>(); ++it) {
-      auto const i = vertexMapper.map(*it);
-      auto const coordinate = it->geometry().corner(0);
-      if (not frictionalBoundary.containsVertex(i)) {
-        restrictions[i] = trivialNonlinearity;
+      auto const i = vertexMapper.index(*it);
+
+      if (not frictionalBoundary.containsVertex(i))
         continue;
-      }
-      auto const fp = std::make_shared<ScalarFriction>(
-          weights[i], normalStress[i], frictionInfo(coordinate));
-      restrictions[i] = std::make_shared<Friction>(fp);
+
+      localToGlobal.emplace_back(i);
+      restrictions.emplace_back(weights[i], weightedNormalStress[i],
+                                frictionInfo(geoToPoint(it->geometry())));
     }
+    assert(restrictions.size() == frictionalBoundary.numVertices());
+    assert(localToGlobal.size() == frictionalBoundary.numVertices());
   }
 
   void updateAlpha(ScalarVector const &alpha) override {
-    for (size_t i = 0; i < restrictions.size(); ++i)
-      restrictions[i]->updateAlpha(alpha[i]);
+    for (size_t j = 0; j < restrictions.size(); ++j)
+      restrictions[j].updateAlpha(alpha[localToGlobal[j]]);
   }
 
   /*
     Return a restriction of the outer function to the i'th node.
   */
-  std::shared_ptr<Friction> restriction(size_t i) const override {
-    return restrictions[i];
+  LocalNonlinearity const &restriction(size_t i) const override {
+    auto const index = indexInSortedRange(localToGlobal, i);
+    if (index == localToGlobal.size())
+      return zeroFriction;
+    return restrictions[index];
   }
 
 private:
-  std::vector<std::shared_ptr<Friction>> restrictions;
+  std::vector<WrappedScalarFriction<block_size, ScalarFriction>> restrictions;
+  std::vector<size_t> localToGlobal;
+  WrappedScalarFriction<block_size, ZeroFunction> const zeroFriction;
 };
 #endif

dune/tectonic/index-in-sorted-range.hh

+#ifndef DUNE_TECTONIC_INDEX_IN_SORTED_RANGE_HH
+#define DUNE_TECTONIC_INDEX_IN_SORTED_RANGE_HH
+
+#include <algorithm>
+
+// returns v.size() if value does not exist
+template <typename T>
+size_t indexInSortedRange(std::vector<T> const &v, T value) {
+  size_t const specialReturnValue = v.size();
+
+  auto const b = std::begin(v);
+  auto const e = std::end(v);
+  auto const lb = std::lower_bound(b, e, value);
+
+  if (lb == e) // all elements are strictly smaller
+    return specialReturnValue;
+
+  if (value < *lb) // value falls between to elements
+    return specialReturnValue;
+
+  return std::distance(b, lb);
+}
+#endif

dune/tectonic/localfriction.hh

@@ -14,39 +14,59 @@
 
 template <size_t dimension> class LocalFriction {
 public:
+  virtual ~LocalFriction() {}
+
   using VectorType = Dune::FieldVector<double, dimension>;
   using MatrixType = Dune::FieldMatrix<double, dimension, dimension>;
 
-  explicit LocalFriction(std::shared_ptr<FrictionPotential> func)
-      : func(func) {}
+  void virtual updateAlpha(double alpha) = 0;
+  double virtual regularity(VectorType const &x) const = 0;
+  double virtual coefficientOfFriction(VectorType const &x) const = 0;
+  void virtual directionalSubDiff(VectorType const &x, VectorType const &v,
+                                  Dune::Solvers::Interval<double> &D) const = 0;
 
-  double operator()(VectorType const &x) const {
-    return func->evaluate(x.two_norm());
-  }
+  void virtual addHessian(VectorType const &x, MatrixType &A) const = 0;
+
+  void virtual addGradient(VectorType const &x, VectorType &y) const = 0;
+
+  void virtual directionalDomain(
+      VectorType const &, VectorType const &,
+      Dune::Solvers::Interval<double> &dom) const = 0;
+};
+
+template <size_t dimension, class ScalarFriction>
+class WrappedScalarFriction : public LocalFriction<dimension> {
+  using VectorType = typename LocalFriction<dimension>::VectorType;
+  using MatrixType = typename LocalFriction<dimension>::MatrixType;
+
+public:
+  template <typename... Args>
+  WrappedScalarFriction(Args... args)
+      : func_(args...) {}
 
-  void updateAlpha(double alpha) { func->updateAlpha(alpha); }
+  void updateAlpha(double alpha) override { func_.updateAlpha(alpha); }
 
-  double regularity(VectorType const &x) const {
+  double regularity(VectorType const &x) const override {
     double const xnorm = x.two_norm();
     if (xnorm <= 0.0)
       return std::numeric_limits<double>::infinity();
 
-    return func->regularity(xnorm);
+    return func_.regularity(xnorm);
   }
 
-  double coefficientOfFriction(VectorType const &x) const {
-    return func->coefficientOfFriction(x.two_norm());
+  double coefficientOfFriction(VectorType const &x) const override {
+    return func_.coefficientOfFriction(x.two_norm());
   }
 
   // directional subdifferential: at u on the line u + t*v
   // u and v are assumed to be non-zero
   void directionalSubDiff(VectorType const &x, VectorType const &v,
-                          Dune::Solvers::Interval<double> &D) const {
+                          Dune::Solvers::Interval<double> &D) const override {
     double const xnorm = x.two_norm();
     if (xnorm <= 0.0)
-      D[0] = D[1] = func->differential(0.0) * v.two_norm();
+      D[0] = D[1] = func_.differential(0.0) * v.two_norm();
     else
-      D[0] = D[1] = func->differential(xnorm) * (x * v) / xnorm;
+      D[0] = D[1] = func_.differential(xnorm) * (x * v) / xnorm;
   }
 
   /** Formula for the derivative:
@@ -65,14 +85,14 @@ template <size_t dimension> class LocalFriction {
       + \frac {H'(|z|)}{|z|} \operatorname{id}
       \f}
   */
-  void addHessian(VectorType const &x, MatrixType &A) const {
+  void addHessian(VectorType const &x, MatrixType &A) const override {
     double const xnorm2 = x.two_norm2();
     double const xnorm = std::sqrt(xnorm2);
     if (xnorm2 <= 0.0)
       return;
 
-    double const H1 = func->differential(xnorm);
-    double const H2 = func->second_deriv(xnorm);
+    double const H1 = func_.differential(xnorm);
+    double const H2 = func_.second_deriv(xnorm);
 
     double const tensorweight = (H2 - H1 / xnorm) / xnorm2;
     double const idweight = H1 / xnorm;
@@ -90,22 +110,22 @@ template <size_t dimension> class LocalFriction {
     }
   }
 
-  void addGradient(VectorType const &x, VectorType &y) const {
+  void addGradient(VectorType const &x, VectorType &y) const override {
     double const xnorm = x.two_norm();
-    if (std::isinf(func->regularity(xnorm)))
+    if (std::isinf(func_.regularity(xnorm)))
       return;
 
     if (xnorm > 0.0)
-      Arithmetic::addProduct(y, func->differential(xnorm) / xnorm, x);
+      Arithmetic::addProduct(y, func_.differential(xnorm) / xnorm, x);
   }
 
   void directionalDomain(VectorType const &, VectorType const &,
-                         Dune::Solvers::Interval<double> &dom) const {
+                         Dune::Solvers::Interval<double> &dom) const override {
     dom[0] = -std::numeric_limits<double>::max();
     dom[1] = std::numeric_limits<double>::max();
   }
 
 private:
-  std::shared_ptr<FrictionPotential> const func;
+  ScalarFriction func_;
 };
 #endif

dune/tectonic/minimisation.hh

@@ -18,7 +18,7 @@ double lineSearch(Functional const &J,
                   typename Functional::LocalVector const &v,
                   Bisection const &bisection) {
   MyDirectionalConvexFunction<typename Functional::Nonlinearity> const JRest(
-      J.alpha * v.two_norm2(), J.b * v, *J.phi, x, v);
+      J.alpha * v.two_norm2(), J.b * v, J.phi, x, v);
   int count;
   return bisection.minimize(JRest, 0.0, 0.0, count);
 }

dune/tectonic/myblockproblem.hh

@@ -4,7 +4,6 @@
 // Based on dune/tnnmg/problem-classes/blocknonlineartnnmgproblem.hh
 
 #include <dune/common/bitsetvector.hh>
-#include <dune/common/nullptr.hh>
 #include <dune/common/parametertree.hh>
 #include <dune/common/fmatrixev.hh>
 
@@ -22,16 +21,16 @@
 /** \brief Base class for problems where each block can be solved with a
  * modified gradient method */
 template <class ConvexProblem>
-class MyBlockProblem : /* NOT PUBLIC */ BlockNonlinearGSProblem<ConvexProblem> {
+class MyBlockProblem : /* not public */ BlockNonlinearGSProblem<ConvexProblem> {
 private:
-  typedef BlockNonlinearGSProblem<ConvexProblem> BNGSP;
+  typedef BlockNonlinearGSProblem<ConvexProblem> Base;
 
 public:
-  using typename BNGSP::ConvexProblemType;
-  using typename BNGSP::LocalMatrixType;
-  using typename BNGSP::LocalVectorType;
-  using typename BNGSP::MatrixType;
-  using typename BNGSP::VectorType;
+  using typename Base::ConvexProblemType;
+  using typename Base::LocalMatrixType;
+  using typename Base::LocalVectorType;
+  using typename Base::MatrixType;
+  using typename Base::VectorType;
 
   size_t static const block_size = ConvexProblem::block_size;
   size_t static const coarse_block_size = block_size;
@@ -56,7 +55,7 @@ class MyBlockProblem : /* NOT PUBLIC */ BlockNonlinearGSProblem<ConvexProblem> {
   };
 
   MyBlockProblem(Dune::ParameterTree const &parset, ConvexProblem &problem)
-      : BNGSP(parset, problem),
+      : Base(parset, problem),
         maxEigenvalues_(problem.f.size()),
         localBisection(0.0, 1.0, 0.0, true, 0.0) {
     for (size_t i = 0; i < problem.f.size(); ++i) {
@@ -104,14 +103,10 @@ class MyBlockProblem : /* NOT PUBLIC */ BlockNonlinearGSProblem<ConvexProblem> {
 
     v /= vnorm; // Rescale for numerical stability
 
-    MyDirectionalConvexFunction<GlobalFriction<MatrixType, VectorType>> const
-    psi(computeDirectionalA(problem_.A, v),
-        computeDirectionalb(problem_.A, problem_.f, u, v), problem_.phi, u, v);
-
+    auto const psi = restrict(problem_.A, problem_.f, u, v, problem_.phi);
     Dune::Solvers::Interval<double> D;
     psi.subDiff(0, D);
-    // NOTE: Numerical instability can actually get us here
-    if (D[1] > 0)
+    if (D[1] > 0) // NOTE: Numerical instability can actually get us here
       return 0;
 
     int bisectionsteps = 0;
@@ -178,17 +173,15 @@ class MyBlockProblem : /* NOT PUBLIC */ BlockNonlinearGSProblem<ConvexProblem> {
           auto const blockEnd = std::end((*col_it)[i]);
           for (auto blockIt = std::begin((*col_it)[i]); blockIt != blockEnd;
                ++blockIt)
-            if (linearization.truncation[row][i] ||
+            if (linearization.truncation[row][i] or
                 linearization.truncation[col][blockIt.index()])
               *blockIt = 0.0;
         }
       }
-
       for (size_t j = 0; j < block_size; ++j)
         if (linearization.truncation[row][j])
           linearization.b[row][j] = 0.0;
     }
-
     for (size_t j = 0; j < block_size; ++j)
       outStream << std::setw(9) << linearization.truncation.countmasked(j);
   }
@@ -202,7 +195,7 @@ class MyBlockProblem : /* NOT PUBLIC */ BlockNonlinearGSProblem<ConvexProblem> {
   std::vector<double> maxEigenvalues_;
 
   // problem data
-  using BNGSP::problem_;
+  using Base::problem_;
 
   Bisection const localBisection;
 
@@ -262,20 +255,17 @@ class MyBlockProblem<ConvexProblem>::IterateObject {
         else
           ui[j] = 0;
 
-      QuadraticEnergy<typename ConvexProblem::NonlinearityType::Friction>
-      localJ(maxEigenvalues_[m], localb, problem.phi.restriction(m));
+      QuadraticEnergy<
+          typename ConvexProblem::NonlinearityType::LocalNonlinearity>
+          localJ(maxEigenvalues_[m], localb, problem.phi.restriction(m));
       minimise(localJ, ui, bisection_);
     }
   }
 
 private:
-  // problem data
   ConvexProblem const &problem;
-
   std::vector<double> maxEigenvalues_;
-
   Bisection const bisection_;
-
   // state data for smoothing procedure used by:
   // setIterate, updateIterate, solveLocalProblem
   VectorType u;

dune/tectonic/mydirectionalconvexfunction.hh

@@ -7,25 +7,6 @@
 #include <dune/fufem/arithmetic.hh>
 #include <dune/solvers/common/interval.hh>
 
-/*
-  1/2 <A(u + hv),u + hv> - <b, u + hv>
-  = 1/2 <Av,v> h^2 - <b - Au, v> h + const.
-
-  localA = <Av,v>
-  localb = <b - Au, v>
-*/
-
-template <class Matrix, class Vector>
-double computeDirectionalA(Matrix const &A, Vector const &v) {
-  return Arithmetic::Axy(A, v, v);
-}
-
-template <class Matrix, class Vector>
-double computeDirectionalb(Matrix const &A, Vector const &b, Vector const &u,
-                           Vector const &v) {
-  return Arithmetic::bmAxy(A, b, u, v);
-}
-
 template <class Nonlinearity> class MyDirectionalConvexFunction {
 public:
   using Vector = typename Nonlinearity::VectorType;
@@ -65,4 +46,23 @@ template <class Nonlinearity> class MyDirectionalConvexFunction {
 
   Dune::Solvers::Interval<double> dom;
 };
+
+/*
+  1/2 <A(u + hv),u + hv> - <b, u + hv>
+  = 1/2 <Av,v> h^2 - <b - Au, v> h + const.
+
+  localA = <Av,v>
+  localb = <b - Au, v>
+*/
+
+template <class Matrix, class Vector, class Nonlinearity>
+MyDirectionalConvexFunction<Nonlinearity> restrict(Matrix const &A,
+                                                   Vector const &b,
+                                                   Vector const &u,
+                                                   Vector const &v,
+                                                   Nonlinearity const &phi) {
+  return MyDirectionalConvexFunction<Nonlinearity>(
+      Arithmetic::Axy(A, v, v), Arithmetic::bmAxy(A, b, u, v), phi, u, v);
+}
+
 #endif

dune/tectonic/quadraticenergy.hh

@@ -8,12 +8,11 @@ template <class NonlinearityTEMPLATE> class QuadraticEnergy {
   using Nonlinearity = NonlinearityTEMPLATE;
   using LocalVector = typename Nonlinearity::VectorType;
 
-  QuadraticEnergy(double alpha, LocalVector const &b,
-                  std::shared_ptr<Nonlinearity const> phi)
+  QuadraticEnergy(double alpha, LocalVector const &b, Nonlinearity const &phi)
       : alpha(alpha), b(b), phi(phi) {}
 
   double const alpha;
   LocalVector const &b;
-  std::shared_ptr<Nonlinearity const> const phi;
+  Nonlinearity const &phi;
 };
 #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],[
-  AC_REQUIRE([AX_BOOST_BASE])
-])
-
-# 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/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

-bin_PROGRAMS = sand-wedge-2D sand-wedge-3D
-
-common_sources = \
-	assemblers.cc \
-	boundary_writer.cc \
-	coupledtimestepper.cc \
-	enumparser.cc \
-	fixedpointiterator.cc \
-	friction_writer.cc \
-	sand-wedge-data/mygeometry.cc \
-	sand-wedge-data/mygrid.cc \
-	solverfactory.cc \
-	state.cc \
-	timestepping.cc \
-	vtk.cc
-
-sand_wedge_2D_SOURCES = $(common_sources) sand-wedge.cc
-sand_wedge_2D_CPPFLAGS = \
-	$(AM_CPPFLAGS) $(BOOST_CPPFLAGS) \
-	-Ddatadir=\"$(abs_srcdir)/sand-wedge-data/\" -DMY_DIM=2
-sand_wedge_3D_SOURCES = $(common_sources) sand-wedge.cc
-sand_wedge_3D_CPPFLAGS = \
-	$(AM_CPPFLAGS) $(BOOST_CPPFLAGS) \
-	-Ddatadir=\"$(abs_srcdir)/sand-wedge-data/\" -DMY_DIM=3
-
-# Some are for clang, others are for gcc
-AM_CXXFLAGS = \
-	-Wall \
-	-Wextra \
-	-Wno-unused-parameter \
-	-Wno-overloaded-virtual \
-	-Wno-new-returns-null \
-	-Wno-unknown-warning-option \
-	-Wno-unknown-pragmas \
-	-Wno-deprecated
-
-AM_CPPFLAGS = \
-	-DDUNE_FMatrix_WITH_CHECKING \
-	$(DUNE_CPPFLAGS) \
-	$(PYTHON_CPPFLAGS) \
-	$(ALUGRID_CPPFLAGS) \
-	$(UG_CPPFLAGS) \
-	-I$(top_srcdir)
-
-# The libraries have to be given in reverse order (most basic libraries
-# last).
-LDADD = \
-	$(DUNE_LDFLAGS) $(DUNE_LIBS) \
-	$(ALUGRID_LIBS) \
-	$(UG_LIBS) \
-	$(PYTHON_LIBS)
-AM_LDFLAGS = \
-	$(DUNE_LDFLAGS) \
-	$(ALUGRID_LDFLAGS) \
-	$(UG_LDFLAGS) \
-	$(PYTHON_LDFLAGS)
-
-if CAIROMM
-AM_CPPFLAGS += $(CAIROMM_CFLAGS)
-LDADD += $(CAIROMM_LIBS)
-endif
-
-include $(top_srcdir)/am/global-rules

src/adaptivetimestepper.hh

-#include "coupledtimestepper.hh"
-
-template <typename T1, typename T2>
-std::pair<T1, T2> clonePair(std::pair<T1, T2> in) {
-  return { in.first->clone(), in.second->clone() };
-}
-
-template <class Factory, class UpdaterPair> class AdaptiveTimeStepper {
-  using StateUpdater = typename UpdaterPair::first_type::element_type;
-  using VelocityUpdater = typename UpdaterPair::second_type::element_type;
-  using Vector = typename Factory::Vector;
-  using ConvexProblem = typename Factory::ConvexProblem;
-  using Nonlinearity = typename ConvexProblem::NonlinearityType;
-
-  using MyCoupledTimeStepper =
-      CoupledTimeStepper<Factory, StateUpdater, VelocityUpdater>;
-
-public:
-  AdaptiveTimeStepper(
-      Factory &factory, Dune::ParameterTree const &parset,
-      std::shared_ptr<Nonlinearity> globalFriction, UpdaterPair &current,
-      std::function<void(double, Vector &)> externalForces,
-      std::function<bool(UpdaterPair &, UpdaterPair &)> mustRefine)
-      : finalTime_(parset.get<double>("problem.finalTime")),
-        relativeTime_(0.0),
-        relativeTau_(1e-6), // FIXME (not really important, though)
-        factory_(factory),
-        parset_(parset),
-        globalFriction_(globalFriction),
-        current_(current),
-        R1_(clonePair(current_)),
-        externalForces_(externalForces),
-        mustRefine_(mustRefine),
-        iterationWriter_("iterations", std::fstream::out) {
-    MyCoupledTimeStepper coupledTimeStepper(finalTime_, factory_, parset_,
-                                            globalFriction_, R1_.first,
-                                            R1_.second, externalForces_);
-    stepAndReport("R1", coupledTimeStepper, relativeTime_, relativeTau_);
-    iterationWriter_ << std::endl;
-  }
-
-  bool reachedEnd() { return relativeTime_ >= 1.0; }
-
-  bool coarsen() {
-    bool didCoarsen = false;
-
-    while (relativeTime_ + relativeTau_ < 1.0) {
-      R2_ = clonePair(R1_);
-      {
-        MyCoupledTimeStepper coupledTimeStepper(finalTime_, factory_, parset_,
-                                                globalFriction_, R2_.first,
-                                                R2_.second, externalForces_);
-        stepAndReport("R2", coupledTimeStepper, relativeTime_ + relativeTau_,
-                      relativeTau_);
-      }
-
-      UpdaterPair C = clonePair(current_);
-      {
-        MyCoupledTimeStepper coupledTimeStepper(finalTime_, factory_, parset_,
-                                                globalFriction_, C.first,
-                                                C.second, externalForces_);
-
-        stepAndReport("C", coupledTimeStepper, relativeTime_,
-                      2.0 * relativeTau_);
-      }
-
-      if (!mustRefine_(C, R2_)) {
-        R2_ = { nullptr, nullptr };
-        R1_ = C;
-        relativeTau_ *= 2.0;
-        didCoarsen = true;
-      } else {
-        break;
-      }
-    }
-    return didCoarsen;
-  }
-
-  void refine() {
-    while (true) {
-      UpdaterPair F2 = clonePair(current_);
-      UpdaterPair F1;
-      {
-        MyCoupledTimeStepper coupledTimeStepper(finalTime_, factory_, parset_,
-                                                globalFriction_, F2.first,
-                                                F2.second, externalForces_);
-        stepAndReport("F1", coupledTimeStepper, relativeTime_,
-                      relativeTau_ / 2.0);
-
-        F1 = clonePair(F2);
-        stepAndReport("F2", coupledTimeStepper,
-                      relativeTime_ + relativeTau_ / 2.0, relativeTau_ / 2.0);
-      }
-
-      if (!mustRefine_(R1_, F2)) {
-        break;
-      } else {
-        R1_ = F1;
-        R2_ = F2;
-        relativeTau_ /= 2.0;
-      }
-    }
-  }
-
-  void advance() {
-    if (!coarsen())
-      refine();
-
-    iterationWriter_ << std::endl;
-
-    current_ = R1_;
-    R1_ = R2_;
-    relativeTime_ += relativeTau_;
-  }
-
-  double getRelativeTime() { return relativeTime_; }
-  double getRelativeTau() { return relativeTau_; }
-
-private:
-  void stepAndReport(std::string type, MyCoupledTimeStepper &stepper,
-                     double rTime, double rTau) {
-    iterationWriter_ << type << " " << stepper.step(rTime, rTau) << " "
-                     << std::flush;
-  }
-
-  double finalTime_;
-  double relativeTime_;
-  double relativeTau_;
-  Factory &factory_;
-  Dune::ParameterTree const &parset_;
-  std::shared_ptr<Nonlinearity> globalFriction_;
-  UpdaterPair &current_;
-  UpdaterPair R1_;
-  UpdaterPair R2_;
-  std::function<void(double, Vector &)> externalForces_;
-  std::function<bool(UpdaterPair &, UpdaterPair &)> mustRefine_;
-  std::fstream iterationWriter_;
-};

src/assemblers.cc

@@ -7,14 +7,15 @@
 #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/computestress.hh>
 #include <dune/fufem/functions/basisgridfunction.hh>
 #include <dune/fufem/functions/constantfunction.hh>
+#include <dune/fufem/functiontools/p0p1interpolation.hh>
 #include <dune/fufem/quadraturerules/quadraturerulecache.hh>
 
 #include <dune/tectonic/frictionpotential.hh>
@@ -33,10 +34,10 @@ 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);
+      frictionalBoundaryMassAssembler(frictionalBoundary);
   vertexAssembler.assembleOperator(frictionalBoundaryMassAssembler,
                                    frictionalBoundaryMass);
 }
@@ -45,22 +46,22 @@ 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);
 
   WeightedMassAssembler<Grid, LocalVertexBasis, LocalVertexBasis,
                         Dune::VirtualFunction<LocalVector, LocalScalarVector>,
                         Dune::ScaledIdentityMatrix<double, dimension>>
-  localWeightedMass(gridView.grid(), densityFunction, quadKey);
+      localWeightedMass(gridView.grid(), densityFunction, quadKey);
   vertexAssembler.assembleOperator(localWeightedMass, M);
 }
 
 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);
+      localStiffness(E, nu);
   vertexAssembler.assembleOperator(localStiffness, A);
 }
 
@@ -68,60 +69,71 @@ 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);
   VariableCoefficientViscosityAssembler<
       Grid, LocalVertexBasis, LocalVertexBasis,
       Dune::VirtualFunction<LocalVector, LocalScalarVector>> const
-  localViscosity(gridView.grid(), shearViscosity, bulkViscosity,
-                 shearViscosityKey, bulkViscosityKey);
+      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) {
+    Vector &f) const {
   L2FunctionalAssembler<Grid, LocalVertexBasis, LocalVector>
-  gravityFunctionalAssembler(gravityField);
+      gravityFunctionalAssembler(gravityField);
   vertexAssembler.assembleFunctional(gravityFunctionalAssembler, f);
 }
 
 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]);
-  ConstantFunction<LocalVector, LocalVector> const fNeumann(localNeumann);
   NeumannBoundaryAssembler<Grid, LocalVector> neumannBoundaryAssembler(
-      fNeumann);
+      std::make_shared<ConstantFunction<LocalVector, LocalVector>>(
+          localNeumann));
   vertexAssembler.assembleBoundaryFunctional(neumannBoundaryAssembler, f,
                                              neumannBoundary);
 }
 
 template <class GridView, int dimension>
-void MyAssembler<GridView, dimension>::assembleNormalStress(
+void MyAssembler<GridView, dimension>::assembleWeightedNormalStress(
     BoundaryPatch<GridView> const &frictionalBoundary,
-    ScalarVector &normalStress, double youngModulus, double poissonRatio,
-    Vector const &displacement) {
-
-  Vector traction;
-  Stress<GridView>::getElasticSurfaceNormalStress // misnomer(!)
-      (frictionalBoundary, displacement, traction, youngModulus, poissonRatio);
-
-  std::vector<typename Vector::block_type> normals;
-  frictionalBoundary.getNormals(normals);
-  for (size_t i = 0; i < traction.size(); ++i) {
-    normalStress[i] = normals[i] * traction[i];
-    if (normalStress[i] > 0.0) {
-      normalStress[i] = 0.0;
-      std::cout << "Warning: Manually reducing positive normal stress to zero."
-                << std::endl;
-    }
+    ScalarVector &weightedNormalStress, double youngModulus,
+    double poissonRatio, Vector const &displacement) const {
+
+  BasisGridFunction<VertexBasis, Vector> displacementFunction(vertexBasis,
+                                                              displacement);
+  Vector traction(cellBasis.size());
+  NormalStressBoundaryAssembler<Grid> tractionBoundaryAssembler(
+      youngModulus, poissonRatio, &displacementFunction, 1);
+  cellAssembler.assembleBoundaryFunctional(tractionBoundaryAssembler, traction,
+                                           frictionalBoundary);
+
+  auto const nodalTractionAverage =
+      interpolateP0ToP1(frictionalBoundary, traction);
+
+  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>
@@ -129,15 +141,15 @@ auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
     Config::FrictionModel frictionModel,
     BoundaryPatch<GridView> const &frictionalBoundary,
     GlobalFrictionData<dimension> const &frictionInfo,
-    ScalarVector const &normalStress)
+    ScalarVector const &weightedNormalStress) const
     -> std::shared_ptr<GlobalFriction<Matrix, Vector>> {
-  // Lump negative normal stress (kludge)
+  // Lumping of the nonlinearity
   ScalarVector weights;
   {
-    ConstantFunction<LocalVector, typename ScalarVector::block_type> const
-    constantOneFunction(1);
     NeumannBoundaryAssembler<Grid, typename ScalarVector::block_type>
-    frictionalBoundaryAssembler(constantOneFunction);
+        frictionalBoundaryAssembler(std::make_shared<
+            ConstantFunction<LocalVector, typename ScalarVector::block_type>>(
+            1));
     vertexAssembler.assembleBoundaryFunctional(frictionalBoundaryAssembler,
                                                weights, frictionalBoundary);
   }
@@ -145,11 +157,11 @@ auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
     case Config::Truncated:
       return std::make_shared<GlobalRateStateFriction<
           Matrix, Vector, TruncatedRateState, GridView>>(
-          frictionalBoundary, gridView, frictionInfo, weights, normalStress);
+          frictionalBoundary, frictionInfo, weights, weightedNormalStress);
     case Config::Regularised:
       return std::make_shared<GlobalRateStateFriction<
           Matrix, Vector, RegularisedRateState, GridView>>(
-          frictionalBoundary, gridView, frictionInfo, weights, normalStress);
+          frictionalBoundary, frictionInfo, weights, weightedNormalStress);
     default:
       assert(false);
   }
@@ -158,7 +170,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);

src/assemblers.hh

@@ -29,8 +29,9 @@ 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;
+  GridView const &gridView;
 
 private:
   using Grid = typename GridView::Grid;
@@ -40,7 +41,6 @@ template <class GridView, int dimension> class MyAssembler {
   using LocalCellBasis = typename CellBasis::LocalFiniteElement;
   using LocalVertexBasis = typename VertexBasis::LocalFiniteElement;
 
-  GridView const &gridView;
   Assembler<CellBasis, CellBasis> cellAssembler;
   Assembler<VertexBasis, VertexBasis> vertexAssembler;
 
@@ -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);
+  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 &normalStress);
+      ScalarVector const &weightedNormalStress) const;
 
   void assembleVonMisesStress(double youngModulus, double poissonRatio,
-                              Vector const &u, ScalarVector &stress);
+                              Vector const &u, ScalarVector &stress) const;
 };
 #endif

src/boundary_writer.cc

-#ifdef HAVE_CONFIG_H
-#include "config.h"
-#endif
-
-#include "boundary_writer.hh"
-#include "tobool.hh"
-
-template <class ScalarVector, class Vector>
-BoundaryWriter<ScalarVector, Vector>::BoundaryWriter(
-    Vector const &vertexCoordinates,
-    Dune::BitSetVector<1> const &_boundaryNodes, std::string prefix,
-    Projector projector)
-    : displacementWriter(prefix + "Displacements", std::fstream::out),
-      velocityWriter(prefix + "Velocities", std::fstream::out),
-      boundaryNodes(_boundaryNodes),
-      projector_(projector) {
-  std::fstream vertexCoordinateWriter(prefix + "Coordinates",
-                                      std::fstream::out);
-  for (size_t i = 0; i < boundaryNodes.size(); ++i)
-    if (toBool(boundaryNodes[i]))
-      vertexCoordinateWriter << vertexCoordinates[i] << std::endl;
-  vertexCoordinateWriter.close();
-}
-
-template <class ScalarVector, class Vector>
-BoundaryWriter<ScalarVector, Vector>::~BoundaryWriter() {
-  displacementWriter.close();
-  velocityWriter.close();
-}
-
-template <class ScalarVector, class Vector>
-void BoundaryWriter<ScalarVector, Vector>::writeKinetics(Vector const &u,
-                                                         Vector const &v) {
-  for (size_t i = 0; i < boundaryNodes.size(); ++i) {
-    if (!toBool(boundaryNodes[i]))
-      continue;
-
-    displacementWriter << projector_(u[i]) << " ";
-    velocityWriter << projector_(v[i]) << " ";
-  }
-
-  displacementWriter << std::endl;
-  velocityWriter << std::endl;
-}
-
-#include "boundary_writer_tmpl.cc"

src/boundary_writer.hh

-#ifndef SRC_BOUNDARY_WRITER_HH
-#define SRC_BOUNDARY_WRITER_HH
-
-#include <fstream>
-#include <string>
-
-#include <dune/common/bitsetvector.hh>
-
-template <class ScalarVector, class Vector> class BoundaryWriter {
-protected:
-  using Projector = std::function<double(typename Vector::block_type const &)>;
-
-public:
-  BoundaryWriter(Vector const &vertexCoordinates,
-                 Dune::BitSetVector<1> const &_boundaryNodes,
-                 std::string prefix, Projector projector);
-
-  virtual ~BoundaryWriter();
-
-  void writeKinetics(Vector const &u, Vector const &v);
-
-protected:
-  std::fstream displacementWriter;
-  std::fstream velocityWriter;
-  Dune::BitSetVector<1> const &boundaryNodes;
-
-  Projector projector_;
-};
-#endif

src/boundary_writer_tmpl.cc

-#ifndef MY_DIM
-#error MY_DIM unset
-#endif
-
-#include "explicitvectors.hh"
-
-template class BoundaryWriter<ScalarVector, Vector>;

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