From b61bc8d92d5ab916a661e07473fe48af2dd12d28 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Tue, 1 Jan 2013 17:28:57 +0000
Subject: [PATCH] Implement method
'secondDerivativeOfDistanceSquaredWRTSecondArgument'
Not tested yet.
[[Imported from SVN: r9082]]
---
dune/gfe/hyperbolichalfspacepoint.hh | 76 +++++++++++++++++++++-------
1 file changed, 59 insertions(+), 17 deletions(-)
diff --git a/dune/gfe/hyperbolichalfspacepoint.hh b/dune/gfe/hyperbolichalfspacepoint.hh
index e1ebf4db..04baba2c 100644
--- a/dune/gfe/hyperbolichalfspacepoint.hh
+++ b/dune/gfe/hyperbolichalfspacepoint.hh
@@ -54,6 +54,15 @@ class HyperbolicHalfspacePoint
return 2/(1-x*x) - 2*x*std::acos(x) / std::pow(1-x*x,1.5);
}
+ /** \brief Compute the second derivative of arccosh^2 without getting unstable for x close to 1 */
+ static T secondDerivativeOfArcCosHSquared(const T& x) {
+ const T eps = 1e-4;
+ if (x < 1+eps) { // regular expression is unstable, use the series expansion instead
+ return -2.0/3 + 8*(x-1)/15;
+ } else
+ return 2/(x*x-1) - 2*x*std::acosh(x) / std::pow(x*x-1,1.5);
+ }
+
/** \brief Compute the third derivative of arccos^2 without getting unstable for x close to 1 */
static T thirdDerivativeOfArcCosSquared(const T& x) {
const T eps = 1e-4;
@@ -183,29 +192,62 @@ public:
Unlike the distance itself the squared distance is differentiable at zero
*/
- static Dune::FieldMatrix<T,N,N> secondDerivativeOfDistanceSquaredWRTSecondArgument(const HyperbolicHalfspacePoint& p, const HyperbolicHalfspacePoint& q) {
+ static Dune::FieldMatrix<T,N,N> secondDerivativeOfDistanceSquaredWRTSecondArgument(const HyperbolicHalfspacePoint& a, const HyperbolicHalfspacePoint& b) {
- T sp = p.data_ * q.data_;
+ // abbreviate notation
+ const Dune::FieldVector<T,N>& p = a.data_;
+ const Dune::FieldVector<T,N>& q = b.data_;
- Dune::FieldVector<T,N> pProjected = q.projectOntoTangentSpace(p.globalCoordinates());
+ T diffNormSquared = (p-q).two_norm2();
- Dune::FieldMatrix<T,N,N> A;
- for (int i=0; i<N; i++)
- for (int j=0; j<N; j++)
- A[i][j] = pProjected[i]*pProjected[j];
+ // Compute first derivative of F
+ Dune::FieldVector<T,N> dFdq;
+ for (size_t i=0; i<N-1; i++)
+ dFdq[i] = ( b.data_[i] - a.data_[i] ) / (a.data_[N-1] * b.data_[N-1]);
- A *= secondDerivativeOfArcCosSquared(sp);
+ dFdq[N-1] = - diffNormSquared / (2*a.data_[N-1]*b.data_[N-1]*b.data_[N-1]) - (a.data_[N-1] - b.data_[N-1]) / (a.data_[N-1]*b.data_[N-1]);
- // Compute matrix B (see notes)
- Dune::FieldMatrix<T,N,N> Pq;
- for (int i=0; i<N; i++)
- for (int j=0; j<N; j++)
- Pq[i][j] = (i==j) - q.data_[i]*q.data_[j];
-
- // Bring it all together
- A.axpy(-1*derivativeOfArcCosSquared(sp)*sp, Pq);
+ // Compute second derivatives of F
+ Dune::FieldMatrix<T,N,N> dFdqdq;
+
+ for (size_t i=0; i<N; i++) {
+
+ for (size_t j=0; j<N; j++) {
+
+ if (i!=N-1 and j!=N-1) {
+
+ dFdqdq[i][j] = (i==j) / (p[N-1]*q[N-1]);
+
+ } else if (i!=N-1 and j==N-1) {
+
+ dFdqdq[i][j] = (p[i] - q[i]) / (p[N-1]*q[N-1]*q[N-1]);
+
+ } else if (i!=N-1 and j==N-1) {
+
+ dFdqdq[i][j] = (p[j] - q[j]) / (p[N-1]*q[N-1]*q[N-1]);
+
+ } else if (i==N-1 and j==N-1) {
+
+ dFdqdq[i][j] = 1/(q[N-1]*q[N-1]) + (p[N-1]*q[N-1]) / (p[N-1]*q[N-1]*q[N-1]) + diffNormSquared / (p[N-1]*q[N-1]*q[N-1]*q[N-1]);
+
+ }
+
+ }
+
+ }
+
+ //
+ T x = 1 + diffNormSquared/ (2*p[N-1]*q[N-1]);
+ T alphaPrime = derivativeOfArcCosHSquared(x);
+ T alphaPrimePrime = secondDerivativeOfArcCosHSquared(x);
- return A;
+ // Sum it all together
+ Dune::FieldMatrix<T,N,N> result;
+ for (size_t i=0; i<N; i++)
+ for (size_t j=0; j<N; j++)
+ result[i][j] = alphaPrimePrime * dFdq[i] * dFdq[j] + alphaPrime * dFdqdq[i][j];
+
+ return result;
}
/** \brief Compute the mixed second derivate \partial d^2 / \partial da db
--
GitLab