CMakeLists.txt

@@ -18,8 +18,7 @@ include(DuneMacros)
 # start a dune project with information from dune.module
 dune_project()
 
-find_package(Boost REQUIRED system filesystem)
-find_package(HDF5 REQUIRED)
+find_package(HDF5 COMPONENTS C REQUIRED)
 
 add_subdirectory("src")
 add_subdirectory("dune")

dune-tectonic.pc.in

@@ -10,6 +10,6 @@ Name: @PACKAGE_NAME@
 Version: @VERSION@
 Description: dune-tectonic module
 URL: http://dune-project.org/
-Requires: dune-common dune-fufem dune-tnnmg
+Requires: dune-common dune-fufem dune-grid dune-istl dune-solvers dune-tnnmg dune-uggrid
 Libs: -L${libdir}
 Cflags: -I${includedir}

dune.module

@@ -3,6 +3,6 @@
 ################################
 
 Module: dune-tectonic
-Version: 2.4-1
+Version: 2.5-1
 Maintainer: elias.pipping@fu-berlin.de
-Depends: dune-common dune-fufem dune-geometry dune-grid  dune-istl dune-solvers dune-tnnmg
+Depends: dune-common dune-fufem dune-grid dune-istl dune-solvers dune-tnnmg dune-uggrid

dune/tectonic/CMakeLists.txt

@@ -10,7 +10,6 @@ install(FILES
   minimisation.hh
   myblockproblem.hh
   mydirectionalconvexfunction.hh
-  pointtractionboundaryassembler.hh
   quadraticenergy.hh
   tectonic.hh
   DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dune/tectonic)

dune/tectonic/frictionpotential.hh

@@ -33,32 +33,32 @@ class TruncatedRateState : public FrictionPotential {
                      FrictionData _fd)
       : fd(_fd), weight(_weight), weightedNormalStress(_weightedNormalStress) {}
 
-  double coefficientOfFriction(double V) const {
+  double coefficientOfFriction(double V) const override {
     if (V <= Vmin)
       return 0.0;
 
     return fd.a * std::log(V / Vmin);
   }
 
-  double differential(double V) const {
+  double differential(double V) const override {
     return weight * fd.C - weightedNormalStress * coefficientOfFriction(V);
   }
 
-  double second_deriv(double V) const {
+  double second_deriv(double V) const override {
     if (V <= Vmin)
       return 0;
 
     return -weightedNormalStress * (fd.a / V);
   }
 
-  double regularity(double V) const {
+  double regularity(double V) const override {
     if (std::abs(V - Vmin) < 1e-14) // TODO
       return std::numeric_limits<double>::infinity();
 
     return std::abs(second_deriv(V));
   }
 
-  void updateAlpha(double alpha) {
+  void updateAlpha(double alpha) override {
     double const logrest = (fd.mu0 + fd.b * alpha) / fd.a;
     Vmin = fd.V0 / std::exp(logrest);
   }
@@ -76,21 +76,23 @@ class RegularisedRateState : public FrictionPotential {
                        FrictionData _fd)
       : fd(_fd), weight(_weight), weightedNormalStress(_weightedNormalStress) {}
 
-  double coefficientOfFriction(double V) const {
+  double coefficientOfFriction(double V) const override {
     return fd.a * std::asinh(0.5 * V / Vmin);
   }
 
-  double differential(double V) const {
+  double differential(double V) const override {
     return weight * fd.C - weightedNormalStress * coefficientOfFriction(V);
   }
 
-  double second_deriv(double V) const {
+  double second_deriv(double V) const override {
     return -weightedNormalStress * fd.a / std::hypot(2.0 * Vmin, V);
   }
 
-  double regularity(double V) const { return std::abs(second_deriv(V)); }
+  double regularity(double V) const override {
+    return std::abs(second_deriv(V));
+  }
 
-  void updateAlpha(double alpha) {
+  void updateAlpha(double alpha) override {
     double const logrest = (fd.mu0 + fd.b * alpha) / fd.a;
     Vmin = fd.V0 / std::exp(logrest);
   }
@@ -104,16 +106,16 @@ class RegularisedRateState : public FrictionPotential {
 
 class ZeroFunction : public FrictionPotential {
 public:
-  double evaluate(double) const { return 0; }
+  double evaluate(double) const override { return 0; }
 
-  double coefficientOfFriction(double s) const { return 0; }
+  double coefficientOfFriction(double s) const override { return 0; }
 
-  double differential(double) const { return 0; }
+  double differential(double) const override { return 0; }
 
-  double second_deriv(double) const { return 0; }
+  double second_deriv(double) const override { return 0; }
 
-  double regularity(double) const { return 0; }
+  double regularity(double) const override { return 0; }
 
-  void updateAlpha(double) {}
+  void updateAlpha(double) override {}
 };
 #endif

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,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 &,
@@ -58,8 +54,7 @@ template <class Matrix, class Vector> class GlobalFriction {
     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,18 @@ 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);
   }
 
   ScalarVector coefficientOfFriction(Vector const &x) const {
     ScalarVector ret(x.size());
     for (size_t i = 0; i < x.size(); ++i)
-      ret[i] = restriction(i)->coefficientOfFriction(x[i]);
+      ret[i] = restriction(i).coefficientOfFriction(x[i]);
     return ret;
   }
 

dune/tectonic/globalratestatefriction.hh

@@ -10,57 +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,
                           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.
   */
-  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>
 
@@ -256,7 +255,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_);
     }

dune/tectonic/pointtractionboundaryassembler.hh

-#ifndef DUNE_TECTONIC_POINTTRACTIONBOUNDARYASSEMBLER_HH
-#define DUNE_TECTONIC_POINTTRACTIONBOUNDARYASSEMBLER_HH
-
-// Based on
-// dune/fufem/assemblers/localassemblers/normalstressboundaryassembler.hh
-
-#include <dune/common/fvector.hh>
-#include <dune/geometry/quadraturerules.hh>
-#include <dune/istl/bvector.hh>
-
-#include <dune/fufem/assemblers/localboundaryassembler.hh>
-#include <dune/fufem/symmetrictensor.hh>
-#include <dune/fufem/functions/virtualgridfunction.hh>
-#include <dune/fufem/quadraturerules/quadraturerulecache.hh>
-
-/** \brief Computes the surface traction of a linear elastic material,
-    pointwise. This is done in such a manner that the sum over all
-    vertices equals the (approximate) integral */
-template <class GridType>
-class PointTractionBoundaryAssembler
-    : public LocalBoundaryAssembler<
-          GridType,
-          Dune::FieldVector<typename GridType::ctype, GridType::dimension>> {
-  static const int dim = GridType::dimension;
-  typedef typename GridType::ctype ctype;
-  typedef typename Dune::FieldVector<ctype, dim> T;
-
-public:
-  typedef typename LocalBoundaryAssembler<GridType, T>::LocalVector LocalVector;
-  typedef VirtualGridFunction<
-      GridType, Dune::FieldVector<ctype, GridType::dimension>> GridFunction;
-
-  //! Create assembler for grid
-  PointTractionBoundaryAssembler(double E, double nu,
-                                 const GridFunction* displacement, int order)
-      : E_(E), nu_(nu), displacement_(displacement), order_(order) {}
-
-  /** \brief Assemble the normal stress at a boundary face.*/
-  template <class TrialLocalFE, class BoundaryIterator>
-  void assemble(const BoundaryIterator& it, LocalVector& localVector,
-                const TrialLocalFE& tFE) {
-
-    localVector = 0.0;
-
-    // geometry of the boundary face
-    const typename BoundaryIterator::Intersection::Geometry segmentGeometry =
-        it->geometry();
-
-    // get quadrature rule
-    const int order = (it->inside()->type().isSimplex())
-                          ? 2 * (order_ - 1)
-                          : 2 * (order_ * dim - 1);
-
-    // get quadrature rule
-    const Dune::QuadratureRule<ctype, dim - 1>& quad =
-        QuadratureRuleCache<ctype, dim - 1>::rule(segmentGeometry.type(), order,
-                                                  false);
-
-    // loop over quadrature points
-    for (size_t pt = 0; pt < quad.size(); ++pt) {
-      // get quadrature point
-      const Dune::FieldVector<ctype, dim - 1>& quadPos = quad[pt].position();
-
-      // get integration factor
-      const ctype integrationElement =
-          segmentGeometry.integrationElement(quadPos);
-
-      // get outer unit normal at quad pos
-      Dune::FieldVector<ctype, dim> normalAtQuadPos =
-          it->unitOuterNormal(quadPos);
-
-      // position of the quadrature point within the element
-      const Dune::FieldVector<ctype, dim> elementQuadPos =
-          it->geometryInInside().global(quadPos);
-
-      // evaluate the displacement gradient at quad point of the element
-      typename GridFunction::DerivativeType localDispGradient;
-      if (displacement_->isDefinedOn(*it->inside()))
-        displacement_->evaluateDerivativeLocal(*it->inside(), elementQuadPos,
-                                               localDispGradient);
-      else
-        displacement_->evaluateDerivative(segmentGeometry.global(quadPos),
-                                          localDispGradient);
-
-      // Compute strain tensor from the displacement gradient
-      SymmetricTensor<dim> strain;
-      computeStrain(localDispGradient, strain);
-
-      // and the stress
-      SymmetricTensor<dim> stress = hookeTimesStrain(strain);
-
-      normalAtQuadPos *= quad[pt].weight() * integrationElement;
-      typename LocalVector::block_type tractionAtQuadPos;
-      stress.umv(normalAtQuadPos, tractionAtQuadPos);
-
-      // Divide the traction
-      tractionAtQuadPos /= segmentGeometry.corners();
-
-      for (size_t k = 0; k < localVector.size(); k++)
-        localVector[k] += tractionAtQuadPos;
-    }
-  }
-
-  ~PointTractionBoundaryAssembler() {}
-
-private:
-  /** \brief Young's modulus */
-  double E_;
-
-  /** \brief The Poisson ratio */
-  double nu_;
-
-  /** \brief The displacement.*/
-  const GridFunction* displacement_;
-
-  /** \brief The gridfunction is usually corresponding to some finite element
-   * solution. This is its order.*/
-  int order_;
-
-  /** \brief Compute linear strain tensor from the deformation gradient. */
-  void computeStrain(const Dune::FieldMatrix<double, dim, dim>& grad,
-                     SymmetricTensor<dim>& strain) const {
-    for (int i = 0; i < dim; ++i) {
-      strain(i, i) = grad[i][i];
-      for (int j = i + 1; j < dim; ++j)
-        strain(i, j) = 0.5 * (grad[i][j] + grad[j][i]);
-    }
-  }
-
-  /** \brief Compute the stress tensor from strain tensor. */
-  SymmetricTensor<dim> hookeTimesStrain(
-      const SymmetricTensor<dim>& strain) const {
-
-    SymmetricTensor<dim> stress = strain;
-    stress *= E_ / (1 + nu_);
-
-    double f = E_ * nu_ / ((1 + nu_) * (1 - 2 * nu_)) * strain.trace();
-
-    for (int i = 0; i < dim; i++)
-      stress[i] += f;
-
-    return stress;
-  }
-};
-
-#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

src/CMakeLists.txt

-set(SOURCE_FILES
+set(SW_SOURCE_FILES
   assemblers.cc
   enumparser.cc
   hdf5/frictionalboundary-writer.cc
@@ -7,9 +7,9 @@ set(SOURCE_FILES
   hdf5/restart-io.cc
   hdf5/surface-writer.cc
   hdf5/time-writer.cc
-  sand-wedge-data/mygeometry.cc
-  sand-wedge-data/mygrid.cc
-  sand-wedge.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
@@ -19,30 +19,24 @@ set(SOURCE_FILES
   time-stepping/state.cc
   vtk.cc
 )
-
-file(MAKE_DIRECTORY "${CMAKE_CURRENT_BINARY_DIR}/sand-wedge-data")
-dune_symlink_to_source_files("sand-wedge-data/boundaryconditions.py")
-dune_symlink_to_source_files("sand-wedge-data/parset.cfg")
-dune_symlink_to_source_files("sand-wedge-data/parset-2D.cfg")
-dune_symlink_to_source_files("sand-wedge-data/parset-3D.cfg")
+set(UGW_SOURCE_FILES
+  assemblers.cc # FIXME
+  one-body-problem-data/mygrid.cc
+  uniform-grid-writer.cc
+  vtk.cc
+)
 
 foreach(_dim 2 3)
-  set(_target sand-wedge-${_dim}D)
-  add_executable(${_target} ${SOURCE_FILES})
-  add_dune_pythonlibs_flags(${_target})
-  add_dune_ug_flags(${_target})
+  set(_sw_target one-body-problem-${_dim}D)
+  set(_ugw_target uniform-grid-writer-${_dim}D)
 
-  set_property(TARGET ${_target} APPEND PROPERTY INCLUDE_DIRECTORIES ${Boost_INCLUDE_DIRS})
+  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})
 
-  target_link_libraries(${_target} ${Boost_FILESYSTEM_LIBRARY})
-  target_link_libraries(${_target} ${Boost_SYSTEM_LIBRARY})
+  add_dune_hdf5_flags(${_sw_target})
 
-  set_property(TARGET ${_target} APPEND PROPERTY INCLUDE_DIRECTORIES ${HDF5_INCLUDE_DIR})
-  target_link_libraries(${_target} ${HDF5_LIBRARIES})
-  target_link_libraries(${_target} "dl") # FIXME: missing from HDF5_LIBRARIES
-
-  set_property(TARGET ${_target} APPEND PROPERTY COMPILE_DEFINITIONS "MY_DIM=${_dim}")
+  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()
-
-# Mark UG as required
-find_package(UG REQUIRED)

src/assemblers.cc

@@ -7,6 +7,7 @@
 #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>
@@ -14,11 +15,11 @@
 #include <dune/fufem/boundarypatch.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>
 #include <dune/tectonic/globalratestatefriction.hh>
-#include <dune/tectonic/pointtractionboundaryassembler.hh>
 
 #include "assemblers.hh"
 
@@ -111,24 +112,28 @@ void MyAssembler<GridView, dimension>::assembleWeightedNormalStress(
 
   BasisGridFunction<VertexBasis, Vector> displacementFunction(vertexBasis,
                                                               displacement);
-  Vector traction;
-  // Note: We assume v = 0 here so that there is no viscous
-  // contribution to the stress.
-  PointTractionBoundaryAssembler<Grid> weightedTractionBoundaryAssembler(
+  Vector traction(cellBasis.size());
+  NormalStressBoundaryAssembler<Grid> tractionBoundaryAssembler(
       youngModulus, poissonRatio, &displacementFunction, 1);
-  vertexAssembler.assembleBoundaryFunctional(weightedTractionBoundaryAssembler,
-                                             traction, frictionalBoundary);
-
-  std::vector<typename Vector::block_type> normals;
-  frictionalBoundary.getNormals(normals);
-  for (size_t i = 0; i < traction.size(); ++i) {
-    weightedNormalStress[i] = normals[i] * traction[i];
-    if (weightedNormalStress[i] > 0.0) {
-      weightedNormalStress[i] = 0.0;
-      std::cout << "Warning: Manually reducing positive normal stress to zero."
-                << std::endl;
-    }
+  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>
@@ -152,13 +157,11 @@ auto MyAssembler<GridView, dimension>::assembleFrictionNonlinearity(
     case Config::Truncated:
       return std::make_shared<GlobalRateStateFriction<
           Matrix, Vector, TruncatedRateState, GridView>>(
-          frictionalBoundary, gridView, frictionInfo, weights,
-          weightedNormalStress);
+          frictionalBoundary, frictionInfo, weights, weightedNormalStress);
     case Config::Regularised:
       return std::make_shared<GlobalRateStateFriction<
           Matrix, Vector, RegularisedRateState, GridView>>(
-          frictionalBoundary, gridView, frictionInfo, weights,
-          weightedNormalStress);
+          frictionalBoundary, frictionInfo, weights, weightedNormalStress);
     default:
       assert(false);
   }

src/enumparser.cc

@@ -9,8 +9,8 @@
 #include "enumparser.hh"
 
 template <class Enum>
-typename Dune::enable_if<
-    !Dune::IsBaseOf<Dune::NotImplemented, StringToEnum<Enum>>::value,
+typename std::enable_if<
+    !std::is_base_of<Dune::NotImplemented, StringToEnum<Enum>>::value,
     std::istream &>::type
 operator>>(std::istream &lhs, Enum &e) {
   std::string s;

src/enumparser.hh

@@ -3,7 +3,7 @@
 
 // Copyright Carsten Graeser 2012
 
-#include <dune/common/typetraits.hh>
+#include <type_traits>
 
 #include <dune/solvers/solvers/solver.hh>
 
@@ -28,8 +28,8 @@ template <> struct StringToEnum<Config::PatchType> {
 };
 
 template <class Enum>
-typename Dune::enable_if<
-    !Dune::IsBaseOf<Dune::NotImplemented, StringToEnum<Enum>>::value,
+typename std::enable_if<
+    !std::is_base_of<Dune::NotImplemented, StringToEnum<Enum>>::value,
     std::istream &>::type
 operator>>(std::istream &lhs, Enum &e);
 #endif

src/gridselector.hh

@@ -12,10 +12,7 @@
 #if !HAVE_ALUGRID
 #error ALUGRID was requested but not found
 #endif
-#pragma clang diagnostic push
-#pragma clang diagnostic ignored "-Wignored-qualifiers"
 #include <dune/grid/alugrid.hh>
-#pragma clang diagnostic pop
 using Grid = Dune::ALUGrid<MY_DIM, MY_DIM, Dune::simplex, Dune::nonconforming>;
 
 #elif WANT_GRID == WANT_UG

src/hdf5-writer.hh

@@ -3,12 +3,11 @@
 
 #include <dune/fufem/functions/basisgridfunction.hh>
 #include <dune/fufem/geometry/convexpolyhedron.hh>
-#include <dune/fufem/hdf5/hdf5file.hh>
+#include <dune/fufem/hdf5/file.hh>
 
 #include "hdf5/frictionalboundary-writer.hh"
 #include "hdf5/iteration-writer.hh"
 #include "hdf5/patchinfo-writer.hh"
-#include "hdf5/restart-io.hh"
 #include "hdf5/surface-writer.hh"
 #include "hdf5/time-writer.hh"
 
@@ -21,7 +20,7 @@ class HDF5Writer {
   using LocalVector = typename Vector::block_type;
 
 public:
-  HDF5Writer(HDF5Grouplike &file, Vector const &vertexCoordinates,
+  HDF5Writer(HDF5::Grouplike &file, Vector const &vertexCoordinates,
              VertexBasis const &vertexBasis, Patch const &surface,
              Patch const &frictionalBoundary,
              ConvexPolyhedron<LocalVector> const &weakPatch)
@@ -32,8 +31,8 @@ class HDF5Writer {
         patchInfoWriter_(file_, vertexBasis, frictionalBoundary, weakPatch),
 #endif
         surfaceWriter_(file_, vertexCoordinates, surface),
-        frictionalBoundaryWriter_(file_, vertexCoordinates, frictionalBoundary),
-        restartWriter_(file_, vertexCoordinates.size()) {
+        frictionalBoundaryWriter_(file_, vertexCoordinates,
+                                  frictionalBoundary) {
   }
 
   template <class Friction>
@@ -53,12 +52,8 @@ class HDF5Writer {
     iterationWriter_.write(programState.timeStep, iterationCount);
   }
 
-  void writeRestart(ProgramState const &programState) {
-    restartWriter_.write(programState);
-  }
-
 private:
-  HDF5Grouplike &file_;
+  HDF5::Grouplike &file_;
 
   IterationWriter iterationWriter_;
   TimeWriter<ProgramState> timeWriter_;
@@ -67,7 +62,5 @@ class HDF5Writer {
 #endif
   SurfaceWriter<ProgramState, GridView> surfaceWriter_;
   FrictionalBoundaryWriter<ProgramState, GridView> frictionalBoundaryWriter_;
-
-  RestartIO<ProgramState> restartWriter_;
 };
 #endif