[Algorit] Use heterogeneous grid (requires UG)

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/gridselector.hh

@@ -5,7 +5,7 @@
 #define WANT_ALUGRID 0
 #define WANT_UG 1
 
-#define WANT_GRID WANT_ALUGRID
+#define WANT_GRID WANT_UG
 
 #if WANT_GRID == WANT_ALUGRID
 

src/sand-wedge-data/parset.cfg

@@ -21,6 +21,7 @@ shearViscosity  = 1e4   # [Pas]
 bulkViscosity   = 1e4   # [Pas]
 
 [boundary.friction]
+smallestDiameter= 1e-3  # [m]
 C               = 10    # [Pa]
 mu0             = 0.7   # [1]
 V0              = 5e-5  # [m/s]
@@ -40,9 +41,6 @@ refinementTolerance = 1e-5
 number = 100000
 scheme = newmark
 
-[grid]
-refinements = 4
-
 [u0.solver]
 tolerance         = 1e-10
 maximumIterations = 100000

src/sand-wedge.cc

@@ -97,6 +97,32 @@ template <class Geometry> double diameter(Geometry const &geometry) {
   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 + 1.0) * smallestDiameter;
+}
+
 int main(int argc, char *argv[]) {
   try {
     Dune::ParameterTree parset;
@@ -110,8 +136,35 @@ int main(int argc, char *argv[]) {
     GridConstructor<Grid> gridConstructor;
     auto grid = gridConstructor.getGrid();
 
-    auto const refinements = parset.get<size_t>("grid.refinements");
-    grid->globalRefine(refinements);
+    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;
+
+        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();
     double maxDiameter = 0.0;