From 71025c09c61a6715adde8603e2471c4721fdb95f Mon Sep 17 00:00:00 2001
From: Elias Pipping <elias.pipping@fu-berlin.de>
Date: Fri, 19 Jul 2013 14:47:09 +0200
Subject: [PATCH] [Algorit] Do away with the horrible near-zero mess
x/|x| isn't that big a deal and the logic was a gigantic mess
---
dune/tectonic/frictionpotential.hh | 9 -----
dune/tectonic/localfriction.hh | 65 +++++++++++++-----------------
2 files changed, 27 insertions(+), 47 deletions(-)
diff --git a/dune/tectonic/frictionpotential.hh b/dune/tectonic/frictionpotential.hh
index c88afb32..c3542ff0 100644
--- a/dune/tectonic/frictionpotential.hh
+++ b/dune/tectonic/frictionpotential.hh
@@ -24,9 +24,6 @@ class FrictionPotentialWrapper {
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;
void virtual updateState(double) = 0;
};
@@ -76,8 +73,6 @@ class FrictionPotential : public FrictionPotentialWrapper {
return std::abs(second_deriv(V));
}
- double smallestPositivePoint() const { return V_m; }
-
// V_m = V_0 exp(-K/a),
// with K = -a log(V_m / V_0) = mu_0 + b [ alpha + log(V_0 / L) ]
void updateState(double state) {
@@ -102,10 +97,6 @@ class TrivialFunction : public FrictionPotentialWrapper {
double regularity(double) const { return 0; }
- double smallestPositivePoint() const {
- return std::numeric_limits<double>::infinity();
- }
-
void updateState(double) {}
};
}
diff --git a/dune/tectonic/localfriction.hh b/dune/tectonic/localfriction.hh
index 870b9e8c..490d663f 100644
--- a/dune/tectonic/localfriction.hh
+++ b/dune/tectonic/localfriction.hh
@@ -13,6 +13,23 @@
#include "frictionpotential.hh"
+// In order to compute (x * y) / |x|, compute x/|x| first
+template <class VectorType>
+double dotFirstNormalised(VectorType const &x, VectorType const &y) {
+ double const xnorm = x.two_norm();
+ if (xnorm <= 0.0)
+ return 0.0;
+
+ size_t const xsize = x.size();
+ assert(xsize == y.size());
+
+ double sum = 0;
+ for (size_t i = 0; i < xsize; ++i)
+ sum += x[i] / xnorm * y[i];
+
+ return sum;
+}
+
namespace Dune {
template <int dimension> class LocalFriction {
public:
@@ -20,20 +37,17 @@ template <int dimension> class LocalFriction {
using MatrixType = FieldMatrix<double, dimension, dimension>;
explicit LocalFriction(shared_ptr<FrictionPotentialWrapper> func)
- : func(func), smp(func->smallestPositivePoint()) {}
+ : func(func) {}
double operator()(VectorType const &x) const {
return func->evaluate(x.two_norm());
}
- void updateState(double state) {
- func->updateState(state);
- smp = func->smallestPositivePoint();
- }
+ void updateState(double state) { func->updateState(state); }
double regularity(VectorType const &x) const {
double const xnorm = x.two_norm();
- if (xnorm < 1e-14 && xnorm >= smp)
+ if (xnorm <= 0.0)
return std::numeric_limits<double>::infinity();
return func->regularity(xnorm);
@@ -44,10 +58,10 @@ template <int dimension> class LocalFriction {
void directionalSubDiff(VectorType const &x, VectorType const &v,
Interval<double> &D) const {
double const xnorm = x.two_norm();
- if (xnorm < smp)
- D[0] = D[1] = 0;
+ if (xnorm <= 0.0)
+ D[0] = D[1] = 0.0;
else
- D[0] = D[1] = func->differential(xnorm) / xnorm * (x * v);
+ D[0] = D[1] = func->differential(xnorm) * dotFirstNormalised(x, v);
}
/** Formula for the derivative:
@@ -69,38 +83,14 @@ 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)
+ if (xnorm2 <= 0.0)
return;
- double const xnorm3 = xnorm * xnorm2;
double const H1 = func->differential(xnorm);
double const H2 = func->second_deriv(xnorm);
- // TODO: potential optimisation: factor out (H1 / xnorm), get rid of xnorm3
- double const weight1 = H2 / xnorm2;
- double const weight2 = -H1 / xnorm3;
- double const weight3 = H1 / xnorm;
-
- // {{{ In what follows, we handle the case 0 * (1/x) = 0 with x
- // close to 0.
- double tensorweight = 0;
-
- if (std::isnan(weight1))
- assert(H2 == 0);
- else
- tensorweight += weight1;
-
- if (std::isnan(weight2))
- assert(H1 == 0);
- else
- tensorweight += weight2;
-
- double idweight = 0;
- if (std::isnan(weight3))
- assert(H1 == 0);
- else
- idweight += weight3;
- // }}}
+ double const tensorweight = (H2 - H1 / xnorm) / xnorm2;
+ double const idweight = H1 / xnorm;
for (int i = 0; i < dimension; ++i)
for (int j = 0; j < i; ++j) {
@@ -120,7 +110,7 @@ template <int dimension> class LocalFriction {
if (std::isinf(func->regularity(xnorm)))
return;
- if (xnorm >= smp)
+ if (xnorm > 0.0)
Arithmetic::addProduct(y, func->differential(xnorm) / xnorm, x);
}
@@ -132,7 +122,6 @@ template <int dimension> class LocalFriction {
private:
shared_ptr<FrictionPotentialWrapper> const func;
- double smp;
};
}
#endif
--
GitLab