From db42041f0ca3a9871bcaf52d4f8e392ce842cc0e Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Wed, 22 Oct 2014 15:33:27 +0200
Subject: [PATCH] [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 | 114 -----------
src/distances.hh | 33 +++
src/polyhedrondistance.hh | 215 ++++++++++++++++++++
src/sand-wedge-data/myglobalfrictiondata.hh | 15 +-
src/sand-wedge-data/mygrid.cc | 35 ++++
src/sand-wedge-data/mygrid.hh | 9 +
src/sand-wedge-data/mygrid_tmpl.cc | 5 +
src/sand-wedge-data/patchfunction.hh | 31 +--
src/sand-wedge-data/weakpatch.hh | 21 ++
src/sand-wedge.cc | 77 +------
10 files changed, 344 insertions(+), 211 deletions(-)
delete mode 100644 dune/tectonic/closestpointprojection.hh
create mode 100644 src/distances.hh
create mode 100644 src/polyhedrondistance.hh
create mode 100644 src/sand-wedge-data/weakpatch.hh
diff --git a/dune/tectonic/closestpointprojection.hh b/dune/tectonic/closestpointprojection.hh
deleted file mode 100644
index 75a513e9..00000000
--- a/dune/tectonic/closestpointprojection.hh
+++ /dev/null
@@ -1,114 +0,0 @@
-// 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/distances.hh b/src/distances.hh
new file mode 100644
index 00000000..2f01e5ab
--- /dev/null
+++ b/src/distances.hh
@@ -0,0 +1,33 @@
+#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
diff --git a/src/polyhedrondistance.hh b/src/polyhedrondistance.hh
new file mode 100644
index 00000000..885735f6
--- /dev/null
+++ b/src/polyhedrondistance.hh
@@ -0,0 +1,215 @@
+#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
diff --git a/src/sand-wedge-data/myglobalfrictiondata.hh b/src/sand-wedge-data/myglobalfrictiondata.hh
index 64db619c..2b14abce 100644
--- a/src/sand-wedge-data/myglobalfrictiondata.hh
+++ b/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_; }
diff --git a/src/sand-wedge-data/mygrid.cc b/src/sand-wedge-data/mygrid.cc
index 604d2db6..aa4175cb 100644
--- a/src/sand-wedge-data/mygrid.cc
+++ b/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"
diff --git a/src/sand-wedge-data/mygrid.hh b/src/sand-wedge-data/mygrid.hh
index 2a208eca..01b17053 100644
--- a/src/sand-wedge-data/mygrid.hh
+++ b/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
diff --git a/src/sand-wedge-data/mygrid_tmpl.cc b/src/sand-wedge-data/mygrid_tmpl.cc
index 0ffd11dc..c076d610 100644
--- a/src/sand-wedge-data/mygrid_tmpl.cc
+++ b/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);
diff --git a/src/sand-wedge-data/patchfunction.hh b/src/sand-wedge-data/patchfunction.hh
index 864ee8c8..25099d52 100644
--- a/src/sand-wedge-data/patchfunction.hh
+++ b/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
diff --git a/src/sand-wedge-data/weakpatch.hh b/src/sand-wedge-data/weakpatch.hh
new file mode 100644
index 00000000..a6162ebe
--- /dev/null
+++ b/src/sand-wedge-data/weakpatch.hh
@@ -0,0 +1,21 @@
+#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
diff --git a/src/sand-wedge.cc b/src/sand-wedge.cc
index fcb31bdd..042ad384 100644
--- a/src/sand-wedge.cc
+++ b/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);
--
GitLab