From 9cfd884337e8ff575617173f07b48136a38b4a14 Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Fri, 31 May 2013 17:33:05 +0200
Subject: [PATCH] SmoothesNorm -> SmallestPositivePoint
---
dune/tectonic/frictionpotential.hh | 14 ++++++++------
dune/tectonic/localfriction.hh | 24 +++++++++++++-----------
2 files changed, 21 insertions(+), 17 deletions(-)
diff --git a/dune/tectonic/frictionpotential.hh b/dune/tectonic/frictionpotential.hh
index 1f150890..a46b1ba5 100644
--- a/dune/tectonic/frictionpotential.hh
+++ b/dune/tectonic/frictionpotential.hh
@@ -18,15 +18,15 @@ class FrictionPotentialWrapper {
double virtual differential(double s) const = 0;
double virtual second_deriv(double x) const = 0;
-
double virtual regularity(double s) const = 0;
- // Whether H(|.|) is smooth at zero
- bool virtual smoothesNorm() const { return false; }
-
double virtual evaluate(double x) const {
DUNE_THROW(NotImplemented, "evaluation not implemented");
}
+
+ // Between 0 and this point, the function is constantly zero (and
+ // thus so are all derivatives)
+ double virtual smallestPositivePoint() const = 0;
};
// phi(V) = V log(V/V_m) - V + V_m if V >= V_m
@@ -87,7 +87,7 @@ class FrictionPotential : public FrictionPotentialWrapper {
return std::abs(second_deriv(V));
}
- bool virtual smoothesNorm() const { return true; }
+ double virtual smallestPositivePoint() const { return V_m; }
private:
double const coefficientProduct;
@@ -104,7 +104,9 @@ class TrivialFunction : public FrictionPotentialWrapper {
double virtual regularity(double) const { return 0; }
- bool virtual smoothesNorm() const { return true; }
+ double virtual smallestPositivePoint() const {
+ return std::numeric_limits<double>::infinity();
+ }
};
}
#endif
diff --git a/dune/tectonic/localfriction.hh b/dune/tectonic/localfriction.hh
index 1791d5a3..02b5b2e2 100644
--- a/dune/tectonic/localfriction.hh
+++ b/dune/tectonic/localfriction.hh
@@ -18,16 +18,16 @@ template <int dimension> class LocalFriction {
using VectorType = FieldVector<double, dimension>;
using MatrixType = FieldMatrix<double, dimension, dimension>;
- LocalFriction(shared_ptr<FrictionPotentialWrapper const> func) : func(func) {}
+ LocalFriction(shared_ptr<FrictionPotentialWrapper const> func)
+ : func(func), smp(func->smallestPositivePoint()) {}
double operator()(VectorType const &x) const {
return func->evaluate(x.two_norm());
}
double regularity(VectorType const &x) const {
- if (!func->smoothesNorm() && x.two_norm() < 1e-14) // TODO: Make this
- // controllable
- // (truncationRadius?)
+ double const xnorm = x.two_norm();
+ if (xnorm < 1e-14 && xnorm >= smp)
return std::numeric_limits<double>::infinity();
return func->regularity(x.two_norm());
@@ -37,12 +37,11 @@ template <int dimension> class LocalFriction {
// u and v are assumed to be non-zero
void directionalSubDiff(VectorType const &x, VectorType const &v,
Interval<double> &D) const {
- if (x.two_norm() == 0) {
- D[0] = D[1] = func->differential(0) * v.two_norm();
- return;
- }
- double const un = x.two_norm();
- D[0] = D[1] = func->differential(un) * (x * v) / un;
+ double const xnorm = x.two_norm();
+ if (xnorm < smp)
+ D[0] = D[1] = 0;
+ else
+ D[0] = D[1] = func->differential(xnorm) / xnorm * (x * v);
}
/** Formula for the derivative:
@@ -64,6 +63,8 @@ template <int dimension> class LocalFriction {
void addHessian(VectorType const &x, MatrixType &A) const {
double const xnorm2 = x.two_norm2();
double const xnorm = std::sqrt(xnorm2);
+ if (xnorm < smp)
+ return;
double const xnorm3 = xnorm * xnorm2;
double const H1 = func->differential(xnorm);
@@ -110,7 +111,7 @@ template <int dimension> class LocalFriction {
void addGradient(VectorType const &x, VectorType &y) const {
double const xnorm = x.two_norm();
- if (xnorm == 0 or std::isinf(func->regularity(xnorm)))
+ if (std::isinf(func->regularity(xnorm)))
return;
y.axpy(func->differential(xnorm) / xnorm, x);
@@ -129,6 +130,7 @@ template <int dimension> class LocalFriction {
private:
shared_ptr<FrictionPotentialWrapper const> const func;
+ double const smp;
};
}
#endif
--
GitLab