From 94265a24f7c4cb4a14714bee8c75f5127e31e12a Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Fri, 1 Apr 2011 10:10:17 +0000
Subject: [PATCH] bugfix: properly implement
secondDerivativeOfDistanceSquaredWRTFirstAndSecondArgument
[[Imported from SVN: r7048]]
---
dune/gfe/unitvector.hh | 27 +++++++++++++--------------
1 file changed, 13 insertions(+), 14 deletions(-)
diff --git a/dune/gfe/unitvector.hh b/dune/gfe/unitvector.hh
index 4251dae8..3e038e44 100644
--- a/dune/gfe/unitvector.hh
+++ b/dune/gfe/unitvector.hh
@@ -199,16 +199,13 @@ public:
*/
static Dune::FieldMatrix<double,N,N> secondDerivativeOfDistanceSquaredWRTFirstAndSecondArgument(const UnitVector& a, const UnitVector& b) {
- Dune::FieldMatrix<double,N,N> result;
-
double sp = a.data_ * b.data_;
// Compute vector A (see notes)
Dune::FieldMatrix<double,1,N> row;
- row[0] = b.globalCoordinates();
- row *= secondDerivativeOfArcCosSquared(sp);
+ row[0] = b.projectOntoTangentSpace(a.globalCoordinates());
- Dune::FieldVector<double,N> tmp = b.projectOntoTangentSpace(a.globalCoordinates());
+ Dune::FieldVector<double,N> tmp = a.projectOntoTangentSpace(b.globalCoordinates());
Dune::FieldMatrix<double,N,1> column;
for (int i=0; i<N; i++) // turn row vector into column vector
column[i] = tmp[i];
@@ -216,21 +213,23 @@ public:
Dune::FieldMatrix<double,N,N> A;
// A = row * column
Dune::FMatrixHelp::multMatrix(column,row,A);
+ A *= secondDerivativeOfArcCosSquared(sp);
// Compute matrix B (see notes)
- Dune::FieldMatrix<double,N,N> B;
+ Dune::FieldMatrix<double,N,N> Pp, Pq;
for (int i=0; i<N; i++)
- for (int j=0; j<N; j++)
- B[i][j] = (i==j) - b.data_[i]*b.data_[j];
+ for (int j=0; j<N; j++) {
+ Pp[i][j] = (i==j) - a.data_[i]*a.data_[j];
+ Pq[i][j] = (i==j) - b.data_[i]*b.data_[j];
+ }
+
+ Dune::FieldMatrix<double,N,N> B;
+ Dune::FMatrixHelp::multMatrix(Pp,Pq,B);
// Bring it all together
- result = A;
- result.axpy(derivativeOfArcCosSquared(sp), B);
+ A.axpy(derivativeOfArcCosSquared(sp), B);
- for (int i=0; i<N; i++)
- result[i] = a.projectOntoTangentSpace(result[i]);
-
- return result;
+ return A;
}
--
GitLab