From d093bb18bf088faa6fc445a1b8bef8b5ae0cab7a Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Sun, 7 Mar 2010 17:55:27 +0000
Subject: [PATCH] factor out derivative of arccos squared into a separate
method
[[Imported from SVN: r5671]]
---
src/unitvector.hh | 20 ++++++++++++--------
1 file changed, 12 insertions(+), 8 deletions(-)
diff --git a/src/unitvector.hh b/src/unitvector.hh
index e4e63c08..446a6a30 100644
--- a/src/unitvector.hh
+++ b/src/unitvector.hh
@@ -11,6 +11,17 @@ class UnitVector
return (x < 1e-4) ? 1 + (x*x/6) : std::sin(x)/x;
}
+ /** \brief Compute the derivative of arccos^2 without getting unstable for x close to 1 */
+ static double derivativeOfArcCosSquared(const double& x) {
+ const double eps = 1e-12;
+ if (x > 1-eps) { // regular expression is unstable, use the series expansion instead
+ return -2 + 2*(x-1)/3 - 4/15*(x-1)*(x-1) + 4/35*(x-1)*(x-1)*(x-1);
+ } else if (x < -1+eps) { // The function is not differentiable
+ DUNE_THROW(Dune::Exception, "Distance is not differentiable for conjugate points!");
+ } else
+ return -2*std::acos(x) / std::sqrt(1-x*x);
+ }
+
public:
/** \brief The type used for coordinates */
@@ -58,17 +69,10 @@ public:
*/
static EmbeddedTangentVector derivativeOfDistanceSquaredWRTSecondArgument(const UnitVector& a, const UnitVector& b) {
double x = a.data_ * b.data_;
- const double eps = 1e-12;
EmbeddedTangentVector result = a.data_;
- if (x > 1-eps) { // regular expression is unstable, use the series expansion instead
- result *= -2 + 2*(x-1)/3 - 4/15*(x-1)*(x-1) + 4/35*(x-1)*(x-1)*(x-1);
- } else if (x < -1+eps) { // a and b are conjugate. The function is not differentiable
- DUNE_THROW(Dune::Exception, "Distance is not differentiable for conjugate points!");
- } else {
- result *= -2*std::acos(x) / std::sqrt(1-x*x);
- }
+ result *= derivativeOfArcCosSquared(x);
// Project gradient onto the tangent plane at b in order to obtain the surface gradient
result = b.projectOntoTangentSpace(result);
--
GitLab