[Algorit] Grid and patch construction overhaul

(1) Determine distance from patch properly (also in general polygon cases)
(2) Move refine()
(3) Rewrite the logic of the patch (allow it to be a polygon)

dune/tectonic/closestpointprojection.hh

-// Based on the class from dune-contact
-
-#ifndef CLOSEST_POINT_PROJECTION_HH
-#define CLOSEST_POINT_PROJECTION_HH
-
-template <class Coordinate> struct BoundarySegment {
-  int nVertices;
-  std::vector<Coordinate> vertices;
-};
-
-template <class Coordinate>
-Coordinate computeClosestPoint(const Coordinate &target,
-                               const BoundarySegment<Coordinate> &segment,
-                               int gaussSeidelIterations) {
-  /**
-   *  The functional to be minimized is: 0.5*norm(\sum_i lambda_i corner_i  -
-   *target)^2
-   *  where lambda_i are barycentric coordinates, i.e. 0\leq lambda_i \leq 1 and
-   *\sum_i lambda_i=1
-   *
-   *  The resulting quadratic minimization problem is given by 0.5
-   *Lambda^T*A*Lambda - l^T*Lambda
-   *  with
-   *      A_ij = (corner_i,corner_j)  and l_i = (corner_i, target)
-  */
-  int nCorners = segment.nVertices;
-  double A[nCorners][nCorners];
-
-  for (int i = 0; i < nCorners; i++)
-    for (int j = 0; j <= i; j++)
-      A[i][j] = segment.vertices[i] * segment.vertices[j];
-
-  // make symmetric
-  for (int i = 0; i < nCorners; i++)
-    for (int j = i + 1; j < nCorners; j++)
-      A[i][j] = A[j][i];
-
-  std::vector<double> l(nCorners);
-  for (int i = 0; i < nCorners; i++)
-    l[i] = target * segment.vertices[i];
-
-  // choose a feasible initial solution
-  std::vector<double> sol(nCorners);
-  for (int i = 0; i < nCorners; i++)
-    sol[i] = 1.0 / nCorners;
-
-  // use a Gauß-Seidel like method for possibly non-convex quadratic problems
-  // with simplex constraints.
-  for (int i = 0; i < gaussSeidelIterations; i++) {
-
-    // compute the residual
-    std::vector<double> rhs = l;
-
-    // compute rhs -= A*sol_
-    for (int k = 0; k < nCorners; k++)
-      for (int j = 0; j < nCorners; j++)
-        rhs[k] -= A[k][j] * sol[j];
-
-    // use the edge vectors as search directions
-    double alpha = 0.0;
-    for (int j1 = 0; j1 < (nCorners - 1); ++j1)
-      for (int j2 = j1 + 1; j2 < nCorners; ++j2) {
-        // compute matrix entry and rhs for edge direction
-        double a = A[j1][j1] - A[j1][j2] - A[j2][j1] + A[j2][j2];
-        double b = rhs[j1] - rhs[j2];
-
-        // compute minimizer for correction
-        if (a > 0) {
-          alpha = b / a;
-          // project alpha such that we stay positiv
-          if ((sol[j2] - alpha) < -1e-14) {
-            alpha = sol[j2];
-          } else if ((alpha < -1e-14) && ((sol[j1] + alpha) < -1e-14)) {
-            alpha = -sol[j1];
-          }
-        } else {
-          // if concave the minimum is achieved at one of the boundaries,
-          // i.e. at sol_[j2] or -sol_[j1]
-
-          double sum = sol[j1] + sol[j2];
-
-          double lValue = 0.5 * A[j1][j1] * sum * sum - rhs[j1] * sum;
-          double uValue = 0.5 * A[j2][j2] * sum * sum - rhs[j2] * sum;
-
-          alpha = (lValue < uValue) ? sol[j2] : -sol[j1];
-        }
-        // apply correction
-        sol[j1] += alpha;
-        sol[j2] -= alpha;
-
-        // update the local residual for corrections
-        for (int p = 0; p < nCorners; p++)
-          rhs[p] -= (A[p][j1] - A[p][j2]) * alpha;
-      }
-  }
-
-  double eps = 1e-5;
-  double sum = 0;
-  for (size_t i = 0; i < sol.size(); i++) {
-    if (sol[i] < -eps)
-      DUNE_THROW(Dune::Exception, "No barycentric coords " << sol[1] << " nr. "
-                                                           << i);
-    sum += sol[i];
-  }
-  if (sum > 1 + eps)
-    DUNE_THROW(Dune::Exception, "No barycentric coords, sum: " << sum);
-
-  Coordinate projected(0);
-  for (int i = 0; i < nCorners; i++)
-    projected.axpy(sol[i], segment.vertices[i]);
-
-  return projected;
-}
-#endif

src/distances.hh

+#ifndef SRC_DISTANCES_HH
+#define SRC_DISTANCES_HH
+
+#include "polyhedrondistance.hh"
+#include "sand-wedge-data/mygeometry.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;
+}
+
+template <class Geometry>
+double distanceToWeakeningRegion(
+    Geometry const &g,
+    ConvexPolyhedron<typename Geometry::GlobalCoordinate> const &weakPatch) {
+  using Coordinate = typename Geometry::GlobalCoordinate;
+
+  ConvexPolyhedron<Coordinate> bsg;
+  bsg.vertices.resize(g.corners());
+  for (size_t i = 0; i < bsg.vertices.size(); ++i)
+    bsg.vertices[i] = g.corner(i);
+
+  return distance(bsg, weakPatch, 1e-6 * MyGeometry::lengthScale);
+}
+#endif

src/polyhedrondistance.hh

+#ifndef SRC_POLYHEDRONDISTANCE_HH
+#define SRC_POLYHEDRONDISTANCE_HH
+
+// Based on the closest point projection from dune-contact
+
+#include <dune/common/dynmatrix.hh>
+#include <dune/solvers/common/arithmetic.hh>
+
+template <class Coordinate> struct ConvexPolyhedron {
+  std::vector<Coordinate> vertices;
+
+  template <class DynamicVector>
+  Coordinate barycentricToGlobal(DynamicVector const &b) const {
+    assert(b.size() == vertices.size());
+    Coordinate r(0);
+    for (size_t i = 0; i < b.size(); i++)
+      Arithmetic::addProduct(r, b[i], vertices[i]);
+    return r;
+  }
+
+  template <class DynamicVector>
+  void sanityCheck(DynamicVector const &b, double tol) const {
+    double sum = 0;
+    for (size_t i = 0; i < b.size(); i++) {
+      if (b[i] < -tol)
+        DUNE_THROW(Dune::Exception, "No barycentric coords " << b[1] << " nr. "
+                                                             << i);
+      sum += b[i];
+    }
+    if (sum > 1 + tol)
+      DUNE_THROW(Dune::Exception, "No barycentric coords, sum: " << sum);
+  }
+};
+
+template <class Coordinate, class Matrix>
+void populateMatrix(ConvexPolyhedron<Coordinate> const &segment,
+                    Matrix &matrix) {
+  size_t const nCorners = segment.vertices.size();
+  for (size_t i = 0; i < nCorners; i++)
+    for (size_t j = 0; j <= i; j++)
+      matrix[i][j] = segment.vertices[i] * segment.vertices[j];
+  // make symmetric
+  for (size_t i = 0; i < nCorners; i++)
+    for (size_t j = i + 1; j < nCorners; j++)
+      matrix[i][j] = matrix[j][i];
+}
+
+/**
+ * The functional to be minimized is
+ *
+ *   0.5*norm(\sum_i lambda_i corner_i - target)^2
+ *
+ * where lambda_i are barycentric coordinates, i.e.
+ *
+ *   0 \le lambda_i \le 1 and \sum_i lambda_i=1
+ *
+ * The resulting quadratic minimization problem is given by
+ *
+ *   0.5 lambda^T*A*lambda - l^T*lambda
+ *
+ * with A_ij = (corner_i,corner_j)  and l_i = (corner_i, target)
+ */
+template <class Coordinate>
+void approximate(Coordinate const &target, Dune::DynamicMatrix<double> const &A,
+                 Dune::DynamicVector<double> const &l,
+                 ConvexPolyhedron<Coordinate> const &segment,
+                 Dune::DynamicVector<double> &sol) {
+  size_t nCorners = segment.vertices.size();
+
+  // compute the residual
+  Dune::DynamicVector<double> rhs = l;
+  // compute rhs -= A*sol_
+  for (size_t k = 0; k < nCorners; k++)
+    for (size_t j = 0; j < nCorners; j++)
+      rhs[k] -= A[k][j] * sol[j];
+
+  // use the edge vectors as search directions
+  double alpha = 0.0;
+  for (size_t j1 = 0; j1 < nCorners - 1; ++j1) {
+    for (size_t j2 = j1 + 1; j2 < nCorners; ++j2) {
+      // compute matrix entry and rhs for edge direction
+      double a = A[j1][j1] - A[j1][j2] - A[j2][j1] + A[j2][j2];
+      double b = rhs[j1] - rhs[j2];
+
+      // compute minimizer for correction
+      if (a > 0) {
+        alpha = b / a;
+
+        double const largestAlpha = sol[j2];
+        double const smallestAlpha = -sol[j1];
+
+        // project alpha such that we stay positive
+        if (alpha > largestAlpha)
+          alpha = largestAlpha;
+        else if (alpha < smallestAlpha)
+          alpha = smallestAlpha;
+      } else {
+        // if the quadratic term vanishes, we have a linear function, which
+        // attains its extrema on the boundary
+        double sum = sol[j1] + sol[j2];
+
+        double lValue = 0.5 * A[j1][j1] * sum * sum - rhs[j1] * sum;
+        double uValue = 0.5 * A[j2][j2] * sum * sum - rhs[j2] * sum;
+
+        alpha = lValue < uValue ? sol[j2] : -sol[j1];
+      }
+      // apply correction
+      sol[j1] += alpha;
+      sol[j2] -= alpha;
+
+      // update the local residual for corrections
+      for (size_t p = 0; p < nCorners; p++)
+        rhs[p] -= (A[p][j1] - A[p][j2]) * alpha;
+    }
+  }
+}
+
+// Point-to-Polyhedron
+template <class Coordinate>
+double distance(const Coordinate &target,
+                const ConvexPolyhedron<Coordinate> &segment,
+                double correctionTolerance) {
+  size_t nCorners = segment.vertices.size();
+  double const barycentricTolerance = 1e-14;
+
+  Dune::DynamicMatrix<double> A(nCorners, nCorners);
+  populateMatrix(segment, A);
+
+  Dune::DynamicVector<double> l(nCorners);
+  for (size_t i = 0; i < nCorners; i++)
+    l[i] = target * segment.vertices[i];
+
+  // choose a feasible initial solution
+  Dune::DynamicVector<double> sol(nCorners);
+  for (size_t i = 0; i < nCorners; i++)
+    sol[i] = 1.0 / nCorners;
+
+  Coordinate oldCoordinates = segment.barycentricToGlobal(sol);
+  Coordinate newCoordinates;
+
+  size_t maxIterations = 1000;
+  size_t iterations;
+  for (iterations = 0; iterations < maxIterations; ++iterations) {
+    approximate(target, A, l, segment, sol);
+
+    newCoordinates = segment.barycentricToGlobal(sol);
+    if ((oldCoordinates - newCoordinates).two_norm() < correctionTolerance)
+      break;
+
+    oldCoordinates = newCoordinates;
+  }
+  assert(iterations < maxIterations);
+
+  segment.sanityCheck(sol, barycentricTolerance);
+  return (target - newCoordinates).two_norm();
+}
+
+// Polyhedron-to-Polyhedron
+template <class Coordinate>
+double distance(const ConvexPolyhedron<Coordinate> &s1,
+                const ConvexPolyhedron<Coordinate> &s2,
+                double valueCorrectionTolerance) {
+  size_t nCorners1 = s1.vertices.size();
+  size_t nCorners2 = s2.vertices.size();
+  double const barycentricTolerance = 1e-14;
+
+  // choose feasible initial solutions
+  Dune::DynamicVector<double> sol1(nCorners1);
+  for (size_t i = 0; i < s1.vertices.size(); i++)
+    sol1[i] = 1.0 / s1.vertices.size();
+  auto c1 = s1.barycentricToGlobal(sol1);
+
+  Dune::DynamicVector<double> sol2(nCorners2);
+  for (size_t i = 0; i < nCorners2; i++)
+    sol2[i] = 1.0 / nCorners2;
+  auto c2 = s2.barycentricToGlobal(sol2);
+
+  double oldDistance = (c2 - c1).two_norm();
+
+  Dune::DynamicMatrix<double> A1(nCorners1, nCorners1);
+  populateMatrix(s1, A1);
+
+  Dune::DynamicMatrix<double> A2(nCorners2, nCorners2);
+  populateMatrix(s2, A2);
+
+  size_t maxIterations = 1000;
+  size_t iterations;
+  for (iterations = 0; iterations < maxIterations; ++iterations) {
+    Dune::DynamicVector<double> l2(nCorners2);
+    for (size_t i = 0; i < nCorners2; i++)
+      l2[i] = c1 * s2.vertices[i];
+
+    approximate(c1, A2, l2, s2, sol2);
+    c2 = s2.barycentricToGlobal(sol2);
+
+    Dune::DynamicVector<double> l1(nCorners1);
+    for (size_t i = 0; i < nCorners1; i++)
+      l1[i] = c2 * s1.vertices[i];
+    approximate(c2, A1, l1, s1, sol1);
+    c1 = s1.barycentricToGlobal(sol1);
+
+    double const newDistance = (c2 - c1).two_norm();
+    assert(newDistance - oldDistance <= 1e-12); // debugging
+    if (oldDistance - newDistance <= valueCorrectionTolerance) {
+      s1.sanityCheck(sol1, barycentricTolerance);
+      s2.sanityCheck(sol2, barycentricTolerance);
+      return newDistance;
+    }
+
+    oldDistance = newDistance;
+  }
+  assert(false);
+}
+
+#endif

src/sand-wedge-data/myglobalfrictiondata.hh

@@ -2,27 +2,26 @@
 #define SRC_SAND_WEDGE_DATA_MYGLOBALFRICTIONDATA_HH
 
 #include <dune/common/function.hh>
-#include <dune/common/fvector.hh>
 
 #include <dune/tectonic/globalfrictiondata.hh>
 
-#include "mygeometry.hh"
 #include "patchfunction.hh"
 
-template <int dimension>
-class MyGlobalFrictionData : public GlobalFrictionData<dimension> {
+template <class LocalVector>
+class MyGlobalFrictionData : public GlobalFrictionData<LocalVector::dimension> {
 private:
-  using typename GlobalFrictionData<dimension>::VirtualFunction;
+  using typename GlobalFrictionData<LocalVector::dimension>::VirtualFunction;
 
 public:
-  MyGlobalFrictionData(Dune::ParameterTree const &parset)
+  MyGlobalFrictionData(Dune::ParameterTree const &parset,
+                       ConvexPolyhedron<LocalVector> const &segment)
       : C_(parset.get<double>("C")),
         L_(parset.get<double>("L")),
         V0_(parset.get<double>("V0")),
         a_(parset.get<double>("strengthening.a"),
-           parset.get<double>("weakening.a")),
+           parset.get<double>("weakening.a"), segment),
         b_(parset.get<double>("strengthening.b"),
-           parset.get<double>("weakening.b")),
+           parset.get<double>("weakening.b"), segment),
         mu0_(parset.get<double>("mu0")) {}
 
   double virtual const &C() const override { return C_; }

src/sand-wedge-data/mygrid.cc

@@ -4,6 +4,7 @@
 
 #include "mygrid.hh"
 #include "midpoint.hh"
+#include "../distances.hh"
 
 #if MY_DIM == 3
 SimplexManager::SimplexManager(unsigned int shift) : shift_(shift) {}
@@ -187,4 +188,38 @@ MyFaces<GridView>::MyFaces(GridView const &gridView)
   });
 }
 
+double computeAdmissibleDiameter(double distance, double smallestDiameter) {
+  return (distance / 0.0125 / MyGeometry::lengthScale + 1.0) * smallestDiameter;
+}
+
+template <class Grid, class LocalVector>
+void refine(Grid &grid, ConvexPolyhedron<LocalVector> const &weakPatch,
+            double smallestDiameter) {
+  bool needRefine = true;
+  while (true) {
+    needRefine = false;
+    for (auto it = grid.template leafbegin<0>();
+         it != grid.template leafend<0>(); ++it) {
+      auto const geometry = it->geometry();
+
+      auto const weakeningRegionDistance =
+          distanceToWeakeningRegion(geometry, weakPatch);
+      auto const admissibleDiameter =
+          computeAdmissibleDiameter(weakeningRegionDistance, smallestDiameter);
+
+      if (diameter(geometry) <= admissibleDiameter)
+        continue;
+
+      needRefine = true;
+      grid.mark(1, *it);
+    }
+    if (!needRefine)
+      break;
+
+    grid.preAdapt();
+    grid.adapt();
+    grid.postAdapt();
+  }
+}
+
 #include "mygrid_tmpl.cc"

src/sand-wedge-data/mygrid.hh

@@ -11,6 +11,8 @@
 #include <dune/fufem/boundarypatch.hh>
 #pragma clang diagnostic pop
 
+#include "../polyhedrondistance.hh"
+
 #include "mygeometry.hh"
 
 template <class GridView> struct MyFaces {
@@ -78,4 +80,11 @@ template <class Grid> class GridConstructor {
 private:
   Dune::GridFactory<Grid> gridFactory;
 };
+
+double computeAdmissibleDiameter(double distance, double smallestDiameter);
+
+template <class Grid, class LocalVector>
+void refine(Grid &grid, ConvexPolyhedron<LocalVector> const &weakPatch,
+            double smallestDiameter);
+
 #endif

src/sand-wedge-data/mygrid_tmpl.cc

@@ -3,6 +3,7 @@
 #endif
 
 #include "explicitgrid.hh"
+#include "explicitvectors.hh"
 
 template class GridConstructor<Grid>;
 
@@ -10,3 +11,7 @@ template struct MyFaces<GridView>;
 
 template MyFaces<GridView> GridConstructor<Grid>::constructFaces(
     GridView const &gridView);
+
+template void refine<Grid, LocalVector>(
+    Grid &grid, ConvexPolyhedron<LocalVector> const &weakPatch,
+    double smallestDiameter);

src/sand-wedge-data/patchfunction.hh

@@ -5,36 +5,27 @@
 #include <dune/common/fvector.hh>
 #include <dune/common/parametertree.hh>
 
+#include "../polyhedrondistance.hh"
+
 class PatchFunction
     : public Dune::VirtualFunction<Dune::FieldVector<double, MY_DIM>,
                                    Dune::FieldVector<double, 1>> {
 private:
-  bool static isBetween(double x, double lower, double upper) {
-    return lower <= x and x <= upper;
-  }
+  using Polyhedron = ConvexPolyhedron<Dune::FieldVector<double, MY_DIM>>;
 
-  bool static isClose(double a, double b) {
-    return std::abs(a - b) < MyGeometry::lengthScale * 1e-14;
-  };
-
-  bool insideRegion2(Dune::FieldVector<double, MY_DIM> const &z) const {
-    assert(isClose(0, z[1]));
-    return isBetween(z[0], _X[0], _Y[0]);
-  }
-
-  Dune::FieldVector<double, MY_DIM> const &_X;
-  Dune::FieldVector<double, MY_DIM> const &_Y;
-
-  double const _v1;
-  double const _v2;
+  double const v1_;
+  double const v2_;
+  Polyhedron const &segment_;
 
 public:
-  PatchFunction(double v1, double v2)
-      : _X(MyGeometry::X), _Y(MyGeometry::Y), _v1(v1), _v2(v2) {}
+  PatchFunction(double v1, double v2, Polyhedron const &segment)
+      : v1_(v1), v2_(v2), segment_(segment) {}
 
   void evaluate(Dune::FieldVector<double, MY_DIM> const &x,
                 Dune::FieldVector<double, 1> &y) const {
-    y = insideRegion2(x) ? _v2 : _v1;
+    y = distance(x, segment_, 1e-6 * MyGeometry::lengthScale) <= 1e-5 ? v2_
+                                                                      : v1_;
   }
 };
+
 #endif

src/sand-wedge-data/weakpatch.hh

+#ifndef SRC_SAND_WEDGE_DATA_WEAKPATCH_HH
+#define SRC_SAND_WEDGE_DATA_WEAKPATCH_HH
+
+template <class LocalVector> ConvexPolyhedron<LocalVector> getWeakPatch() {
+  ConvexPolyhedron<LocalVector> weakPatch;
+#if MY_DIM == 3
+  weakPatch.vertices.resize(4);
+  weakPatch.vertices[0] = weakPatch.vertices[2] = MyGeometry::X;
+  weakPatch.vertices[1] = weakPatch.vertices[3] = MyGeometry::Y;
+  for (size_t k = 0; k < 2; ++k) {
+    weakPatch.vertices[k][2] = -MyGeometry::depth / 2.0;
+    weakPatch.vertices[k + 2][2] = MyGeometry::depth / 2.0;
+  }
+#else
+  weakPatch.vertices.resize(2);
+  weakPatch.vertices[0] = MyGeometry::X;
+  weakPatch.vertices[1] = MyGeometry::Y;
+#endif
+  return weakPatch;
+};
+#endif

src/sand-wedge.cc

@@ -58,6 +58,7 @@
 #include "assemblers.hh"
 #include "tobool.hh"
 #include "coupledtimestepper.hh"
+#include "distances.hh"
 #include "enumparser.hh"
 #include "enums.hh"
 #include "fixedpointiterator.hh"
@@ -68,6 +69,7 @@
 #include "sand-wedge-data/myglobalfrictiondata.hh"
 #include "sand-wedge-data/mygrid.hh"
 #include "sand-wedge-data/special_writer.hh"
+#include "sand-wedge-data/weakpatch.hh"
 #include "solverfactory.hh"
 #include "state.hh"
 #include "timestepping.hh"
@@ -82,47 +84,6 @@ void initPython() {
   Python::run("sys.path.append('" datadir "')");
 }
 
-template <class Geometry> double diameter(Geometry const &geometry) {
-  auto const numCorners = geometry.corners();
-  std::vector<typename Geometry::GlobalCoordinate> corners(numCorners);
-  for (int i = 0; i < numCorners; ++i)
-    corners[i] = geometry.corner(i);
-
-  TwoNorm<typename Geometry::GlobalCoordinate> twoNorm;
-
-  double diameter = 0.0;
-  for (int i = 0; i < numCorners; ++i)
-    for (int j = 0; j < i; ++j)
-      diameter = std::max(diameter, twoNorm.diff(corners[i], corners[j]));
-  return diameter;
-}
-
-#include <dune/tectonic/closestpointprojection.hh>
-
-// we assume the weakening region to be a line segment for now
-template <class Geometry> double distanceToWeakeningRegion(Geometry const &g) {
-  TwoNorm<typename Geometry::GlobalCoordinate> twoNorm;
-
-  BoundarySegment<typename Geometry::GlobalCoordinate> bs;
-  bs.nVertices = 2;
-  bs.vertices.resize(bs.nVertices);
-  bs.vertices[0] = MyGeometry::X;
-  bs.vertices[1] = MyGeometry::Y;
-
-  double distance = std::numeric_limits<double>::infinity();
-  auto const numCorners = g.corners();
-  for (int i = 0; i < numCorners; ++i) {
-    auto const corner = g.corner(i);
-    auto const projection = computeClosestPoint(corner, bs, 10); // FIXME
-    distance = std::min(distance, twoNorm.diff(corner, projection));
-  }
-  return distance;
-}
-
-double computeAdmissibleDiameter(double distance, double smallestDiameter) {
-  return (distance / 0.0125 / MyGeometry::lengthScale + 1.0) * smallestDiameter;
-}
-
 int main(int argc, char *argv[]) {
   try {
     Dune::ParameterTree parset;
@@ -142,38 +103,15 @@ int main(int argc, char *argv[]) {
     using ScalarVector = MyAssembler::ScalarVector;
     using LocalScalarVector = ScalarVector::block_type;
 
+    auto weakPatch = getWeakPatch<LocalVector>();
+
     // {{{ Set up grid
     GridConstructor<Grid> gridConstructor;
     auto grid = gridConstructor.getGrid();
 
-    auto const smallestDiameter =
-        parset.get<double>("boundary.friction.smallestDiameter");
-    bool needRefine = true;
-
-    while (true) {
-      needRefine = false;
-      for (auto it = grid->template leafbegin<0>();
-           it != grid->template leafend<0>(); ++it) {
-        auto const geometry = it->geometry();
-
-        auto const weakeningRegionDistance =
-            distanceToWeakeningRegion(geometry);
-        auto const admissibleDiameter = computeAdmissibleDiameter(
-            weakeningRegionDistance, smallestDiameter);
-
-        if (diameter(geometry) <= admissibleDiameter)
-          continue;
+    refine(*grid, weakPatch,
+           parset.get<double>("boundary.friction.smallestDiameter"));
 
-        needRefine = true;
-        grid->mark(1, *it);
-      }
-      if (!needRefine)
-        break;
-
-      grid->preAdapt();
-      grid->adapt();
-      grid->postAdapt();
-    }
     auto const refinements = grid->maxLevel();
 
     double minDiameter = std::numeric_limits<double>::infinity();
@@ -318,7 +256,8 @@ int main(int argc, char *argv[]) {
     myAssembler.assembleNormalStress(frictionalBoundary, normalStress,
                                      body.getYoungModulus(),
                                      body.getPoissonRatio(), u_initial);
-    MyGlobalFrictionData<dims> frictionInfo(parset.sub("boundary.friction"));
+    MyGlobalFrictionData<LocalVector> frictionInfo(
+        parset.sub("boundary.friction"), weakPatch);
     auto myGlobalFriction = myAssembler.assembleFrictionNonlinearity(
         parset.get<Config::FrictionModel>("boundary.friction.frictionModel"),
         frictionalBoundary, frictionInfo, normalStress);