[Cleanup] Store only the necessary DOFs

dune/tectonic/geocoordinate.hh

+#ifndef DUNE_TECTONIC_GEOCOORDINATE_HH
+#define DUNE_TECTONIC_GEOCOORDINATE_HH
+
+// tiny helper to make a common piece of code pleasanter to read
+
+template <class Geometry>
+typename Geometry::GlobalCoordinate geoToPoint(Geometry geo) {
+  assert(geo.corners() == 1);
+  return geo.corner(0);
+}
+
+#endif

dune/tectonic/globalfriction.hh

@@ -22,7 +22,7 @@ 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;
@@ -35,7 +35,7 @@ template <class Matrix, class Vector> class GlobalFriction {
   /*
     Return a restriction of the outer function to the i'th node.
   */
-  Friction const 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)

dune/tectonic/globalratestatefriction.hh

@@ -10,14 +10,16 @@
 #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;
@@ -27,40 +29,43 @@ class GlobalRateStateFriction : public GlobalFriction<Matrix, Vector> {
                           GlobalFrictionData<block_size> const &frictionInfo,
                           ScalarVector const &weights,
                           ScalarVector const &weightedNormalStress)
-      : restrictions(weightedNormalStress.size()) {
-    auto zeroNonlinearity =
-        std::make_shared<Friction>(std::make_shared<ZeroFunction>());
+      : 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.index(*it);
-      auto const coordinate = it->geometry().corner(0);
-      if (not frictionalBoundary.containsVertex(i)) {
-        restrictions[i] = zeroNonlinearity;
+
+      if (not frictionalBoundary.containsVertex(i))
         continue;
-      }
-      auto const fp = std::make_shared<ScalarFriction>(
-          weights[i], weightedNormalStress[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.
   */
-  Friction const &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/myblockproblem.hh

@@ -256,7 +256,8 @@ class MyBlockProblem<ConvexProblem>::IterateObject {
         else
           ui[j] = 0;
 
-      QuadraticEnergy<typename ConvexProblem::NonlinearityType::Friction>
+      QuadraticEnergy<
+          typename ConvexProblem::NonlinearityType::LocalNonlinearity>
           localJ(maxEigenvalues_[m], localb, problem.phi.restriction(m));
       minimise(localJ, ui, bisection_);
     }

src/sand-wedge.cc

@@ -46,6 +46,7 @@
 #include <dune/solvers/solvers/solver.hh>
 #include <dune/tnnmg/problem-classes/convexproblem.hh>
 
+#include <dune/tectonic/geocoordinate.hh>
 #include <dune/tectonic/myblockproblem.hh>
 #include <dune/tectonic/globalfriction.hh>
 #include <dune/fufem/hdf5/hdf5file.hh>
@@ -236,11 +237,8 @@ int main(int argc, char *argv[]) {
     {
       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.index(*it)] = geometry.corner(0);
-      }
+      for (auto it = leafView.begin<dims>(); it != leafView.end<dims>(); ++it)
+        vertexCoordinates[vertexMapper.index(*it)] = geoToPoint(it->geometry());
     }
 
     HDF5Writer<ProgramState<Vector, ScalarVector>,