From 118743c65fd9c947ca5d497cd74594d74cf95ae3 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Tue, 7 Oct 2014 00:15:43 +0200
Subject: [PATCH] [Algorit] Use heterogeneous grid (requires UG)
---
dune/tectonic/closestpointprojection.hh | 114 ++++++++++++++++++++++++
src/gridselector.hh | 2 +-
src/sand-wedge-data/parset.cfg | 4 +-
src/sand-wedge.cc | 57 +++++++++++-
4 files changed, 171 insertions(+), 6 deletions(-)
create mode 100644 dune/tectonic/closestpointprojection.hh
diff --git a/dune/tectonic/closestpointprojection.hh b/dune/tectonic/closestpointprojection.hh
new file mode 100644
index 00000000..75a513e9
--- /dev/null
+++ b/dune/tectonic/closestpointprojection.hh
@@ -0,0 +1,114 @@
+// 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
diff --git a/src/gridselector.hh b/src/gridselector.hh
index 0e27a54f..a07b0610 100644
--- a/src/gridselector.hh
+++ b/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
diff --git a/src/sand-wedge-data/parset.cfg b/src/sand-wedge-data/parset.cfg
index 5d254440..ab6670d0 100644
--- a/src/sand-wedge-data/parset.cfg
+++ b/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
diff --git a/src/sand-wedge.cc b/src/sand-wedge.cc
index cf831898..799b6509 100644
--- a/src/sand-wedge.cc
+++ b/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;
--
GitLab