From 4e03c456bf664b77d26955301cea276c350a1a33 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Thu, 12 Jan 2012 09:08:42 +0000
Subject: [PATCH] Reimplement method
derivativeOfDistanceSquaredWRTSecondArgument.
Now it properly handles the situation that p and q may be on
different sheets of the double cover.
[[Imported from SVN: r8367]]
---
dune/gfe/rotation.hh | 34 +++++++++++++++-------------------
1 file changed, 15 insertions(+), 19 deletions(-)
diff --git a/dune/gfe/rotation.hh b/dune/gfe/rotation.hh
index 6a52b7d7..e4f030a9 100644
--- a/dune/gfe/rotation.hh
+++ b/dune/gfe/rotation.hh
@@ -621,30 +621,26 @@ public:
}
+ /** \brief Compute the derivative of the squared distance function with respect to the second argument
+ *
+ * The squared distance function is
+ * \f[ 4 \arccos^2 (|<p,q>|) \f]
+ * Deriving this with respect to the second coordinate gives
+ * \f[ 4 (\arccos^2)'(x)|_{x = |<p,q>|} * P_qp \f]
+ * The whole thing has to be multiplied by -1 if \f$ <p,q> \f$ is negative
+ */
static EmbeddedTangentVector derivativeOfDistanceSquaredWRTSecondArgument(const Rotation<T,3>& p,
const Rotation<T,3>& q) {
+ T sp = p.globalCoordinates() * q.globalCoordinates();
- Rotation<T,3> pInv = p;
- pInv.invert();
-
- // the forth component of pInv times q
- T pInvq_4 = - pInv[0]*q[0] - pInv[1]*q[1] - pInv[2]*q[2] + pInv[3]*q[3];
-
- T arccosSquaredDer_pInvq_4 = derivativeOfArcCosSquared(pInvq_4);
+ // Compute the projection of p onto the tangent space of q
+ EmbeddedTangentVector result = p;
+ result.axpy(-1*(q*p), q);
- EmbeddedTangentVector result;
- result[0] = -4 * arccosSquaredDer_pInvq_4 * pInv[0];
- result[1] = -4 * arccosSquaredDer_pInvq_4 * pInv[1];
- result[2] = -4 * arccosSquaredDer_pInvq_4 * pInv[2];
- result[3] = 4 * arccosSquaredDer_pInvq_4 * pInv[3];
+ // The ternary operator comes from the derivative of the absolute value function
+ result *= 4 * derivativeOfArcCosSquared(std::fabs(sp)) * ( (sp<0) ? -1 : 1 );
- // project onto the tangent space at q
- EmbeddedTangentVector projectedResult = result;
- projectedResult.axpy(-1*(q*result), q);
-
- assert(std::fabs(projectedResult * q) < 1e-7);
-
- return projectedResult;
+ return result;
}
/** \brief Compute the Hessian of the squared distance function keeping the first argument fixed
--
GitLab