From 6f3b968988350b1f08befcfff226f698a56e751b Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Sat, 30 Jan 2016 06:56:05 +0100
Subject: [PATCH] Implement new method distanceSquared, between a double vector
and an adouble one
This is needed for the new gradientflow application. The standard 'distance' method
doesn't cut it, because it is not differentiable near zero. Therefore, even the
differentiable expression 'distance*distance' will fail to be differentiable for
an automatic-differentiation system. I think in the long run we should replace
'distance' by 'distanceSquared' everywhere it is used.
Unfortunately, this patch is hacky: it only adds the method for the double/adouble
combination.
---
dune/gfe/unitvector.hh | 26 ++++++++++++++++++++++++++
1 file changed, 26 insertions(+)
diff --git a/dune/gfe/unitvector.hh b/dune/gfe/unitvector.hh
index 6621a674..168d3153 100644
--- a/dune/gfe/unitvector.hh
+++ b/dune/gfe/unitvector.hh
@@ -26,6 +26,15 @@ class UnitVector
return (x < 1e-4) ? 1 - (x*x/6) : std::sin(x)/x;
}
+ /** \brief Compute arccos^2 without using the (at 1) nondifferentiable function acos x close to 1 */
+ static T arcCosSquared(const T& x) {
+ const T eps = 1e-4;
+ if (x > 1-eps) { // acos is not differentiable, use the series expansion instead
+ return -2 * (x-1) + 1.0/3 * (x-1)*(x-1) - 4.0/45 * (x-1)*(x-1)*(x-1);
+ } else
+ return Dune::Power<2>::eval(std::acos(x));
+ }
+
/** \brief Compute the derivative of arccos^2 without getting unstable for x close to 1 */
static T derivativeOfArcCosSquared(const T& x) {
const T eps = 1e-4;
@@ -177,6 +186,23 @@ public:
return std::acos(x);
}
+#if ADOLC_ADOUBLE_H
+ /** \brief Squared length of the great arc connecting the two points */
+ static adouble distanceSquared(const UnitVector<double,N>& a, const UnitVector<adouble,N>& b)
+ {
+ // Not nice: we are in a class for unit vectors, but the class is actually
+ // supposed to handle perturbations of unit vectors as well. Therefore
+ // we normalize here.
+ adouble x = a.data_ * b.data_ / (a.data_.two_norm()*b.data_.two_norm());
+
+ // paranoia: if the argument is just eps larger than 1 acos returns NaN
+ x = std::min(x,1.0);
+
+ // Special implementation that remains AD-differentiable near x==1
+ return arcCosSquared(x);
+ }
+#endif
+
/** \brief Compute the gradient of the squared distance function keeping the first argument fixed
Unlike the distance itself the squared distance is differentiable at zero
--
GitLab