diff --git a/dune/tectonic/geocoordinate.hh b/dune/tectonic/geocoordinate.hh
new file mode 100644
index 0000000000000000000000000000000000000000..641bd50fee29129ecbf2c1db90c8d8f00e5336f9
--- /dev/null
+++ b/dune/tectonic/geocoordinate.hh
@@ -0,0 +1,12 @@
+#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
diff --git a/dune/tectonic/globalfriction.hh b/dune/tectonic/globalfriction.hh
index 538eec785150700e6351672ceaf2d6c1611789a0..41ddd18f892ea5b1130a63d9620ab5502d4f7421 100644
--- a/dune/tectonic/globalfriction.hh
+++ b/dune/tectonic/globalfriction.hh
@@ -22,7 +22,7 @@ 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;
@@ -35,7 +35,7 @@ template <class Matrix, class Vector> class GlobalFriction {
/*
Return a restriction of the outer function to the i'th node.
*/
- Friction const 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)
diff --git a/dune/tectonic/globalratestatefriction.hh b/dune/tectonic/globalratestatefriction.hh
index 64684f2c6a18702e4bae7daa69b12be7a0dd294c..d6c61470cd76fb0219ab0c72bdaa9ab8fb37a42f 100644
--- a/dune/tectonic/globalratestatefriction.hh
+++ b/dune/tectonic/globalratestatefriction.hh
@@ -10,14 +10,16 @@
#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;
@@ -27,40 +29,43 @@ class GlobalRateStateFriction : public GlobalFriction<Matrix, Vector> {
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.
*/
- Friction const &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
diff --git a/dune/tectonic/index-in-sorted-range.hh b/dune/tectonic/index-in-sorted-range.hh
new file mode 100644
index 0000000000000000000000000000000000000000..ecc4038f335f66a1ed7caad80f354ab4f10c96f2
--- /dev/null
+++ b/dune/tectonic/index-in-sorted-range.hh
@@ -0,0 +1,23 @@
+#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
diff --git a/dune/tectonic/localfriction.hh b/dune/tectonic/localfriction.hh
index 5cda30a6bca5ad82baf149dacb42782d45728bc5..4bdd054e17438d0c5bcb89dfc4d1c0778b9018e5 100644
--- a/dune/tectonic/localfriction.hh
+++ b/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
diff --git a/dune/tectonic/myblockproblem.hh b/dune/tectonic/myblockproblem.hh
index 50fbf2e7864b97cec97ee0d1562e7f12cf623ef2..b0de9b122469cfafc54776fca5a12b332185bb4e 100644
--- a/dune/tectonic/myblockproblem.hh
+++ b/dune/tectonic/myblockproblem.hh
@@ -256,7 +256,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_);
}
diff --git a/src/sand-wedge.cc b/src/sand-wedge.cc
index 692397ea88e9ba3ba6e4863bae1869ce9fee4baf..e89adce903dd295018e282e935a51c34cd4e9040 100644
--- a/src/sand-wedge.cc
+++ b/src/sand-wedge.cc
@@ -46,6 +46,7 @@
#include <dune/solvers/solvers/solver.hh>
#include <dune/tnnmg/problem-classes/convexproblem.hh>
+#include <dune/tectonic/geocoordinate.hh>
#include <dune/tectonic/myblockproblem.hh>
#include <dune/tectonic/globalfriction.hh>
#include <dune/fufem/hdf5/hdf5file.hh>
@@ -236,11 +237,8 @@ int main(int argc, char *argv[]) {
{
Dune::MultipleCodimMultipleGeomTypeMapper<
GridView, Dune::MCMGVertexLayout> const vertexMapper(leafView);
- for (auto it = leafView.begin<dims>(); it != leafView.end<dims>(); ++it) {
- auto const geometry = it->geometry();
- assert(geometry.corners() == 1);
- vertexCoordinates[vertexMapper.index(*it)] = geometry.corner(0);
- }
+ for (auto it = leafView.begin<dims>(); it != leafView.end<dims>(); ++it)
+ vertexCoordinates[vertexMapper.index(*it)] = geoToPoint(it->geometry());
}
HDF5Writer<ProgramState<Vector, ScalarVector>,