From c218cd9e9b0ba7f49b2c894ae150e1cfc43e56cc Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Fri, 13 Dec 2013 09:50:31 +0000
Subject: [PATCH] Use SymmetricMatrix for Hessian of the distance functional
exclusively
[[Imported from SVN: r9588]]
---
dune/gfe/averagedistanceassembler.hh | 19 +++++++++++--------
dune/gfe/localgeodesicfefunction.hh | 5 +++--
dune/gfe/realtuple.hh | 9 +++++----
dune/gfe/rotation.hh | 15 ++++++++-------
dune/gfe/unitvector.hh | 15 ++++++++-------
5 files changed, 35 insertions(+), 28 deletions(-)
diff --git a/dune/gfe/averagedistanceassembler.hh b/dune/gfe/averagedistanceassembler.hh
index c426a6ba..b4def6ba 100644
--- a/dune/gfe/averagedistanceassembler.hh
+++ b/dune/gfe/averagedistanceassembler.hh
@@ -73,20 +73,23 @@ public:
Dune::FieldMatrix<ctype,embeddedSize,embeddedSize>& matrix) const
{
matrix = 0;
- for (size_t i=0; i<coefficients_.size(); i++)
- matrix.axpy(weights_[i], TargetSpace::secondDerivativeOfDistanceSquaredWRTSecondArgument(coefficients_[i], x));
+ for (size_t i=0; i<coefficients_.size(); i++) {
+ Dune::SymmetricMatrix<ctype,embeddedSize> tmp = TargetSpace::secondDerivativeOfDistanceSquaredWRTSecondArgument(coefficients_[i], x);
+ for (int j=0; j<embeddedSize; j++)
+ for (int k=0; k<=j; k++) {
+ matrix[j][k] += weights_[i] * tmp(j,k);
+ if (j!=k)
+ matrix[k][j] += weights_[i] * tmp(j,k);
+ }
+ }
}
void assembleEmbeddedHessian(const TargetSpace& x,
Dune::SymmetricMatrix<ctype,embeddedSize>& matrix) const
{
matrix = 0;
- for (size_t i=0; i<coefficients_.size(); i++) {
- Dune::FieldMatrix<ctype,embeddedSize,embeddedSize> tmp = TargetSpace::secondDerivativeOfDistanceSquaredWRTSecondArgument(coefficients_[i], x);
- for (size_t j=0; j<embeddedSize; j++)
- for (size_t k=0; k<=j; k++)
- matrix(j,k) += weights_[i] * tmp[j][k];
- }
+ for (size_t i=0; i<coefficients_.size(); i++)
+ matrix.axpy(weights_[i], TargetSpace::secondDerivativeOfDistanceSquaredWRTSecondArgument(coefficients_[i], x));
}
void assembleHessian(const TargetSpace& x,
diff --git a/dune/gfe/localgeodesicfefunction.hh b/dune/gfe/localgeodesicfefunction.hh
index 0b8b5434..9eadf45e 100644
--- a/dune/gfe/localgeodesicfefunction.hh
+++ b/dune/gfe/localgeodesicfefunction.hh
@@ -497,10 +497,11 @@ evaluateDerivativeOfGradientWRTCoefficient(const Dune::FieldVector<ctype, dim>&
TensorSSD<RT, embeddedDim,embeddedDim> dqdwF(coefficients_.size());
for (size_t k=0; k<coefficients_.size(); k++) {
- Dune::FieldMatrix<RT,embeddedDim,embeddedDim> hesse = TargetSpace::secondDerivativeOfDistanceSquaredWRTSecondArgument(coefficients_[k], q);
+ Dune::SymmetricMatrix<RT,embeddedDim> hesse = TargetSpace::secondDerivativeOfDistanceSquaredWRTSecondArgument(coefficients_[k], q);
for (int i=0; i<embeddedDim; i++)
for (int j=0; j<embeddedDim; j++)
- dqdwF(i, j, k) = hesse[i][j];
+ dqdwF(i, j, k) = dqdwF(j, i, k) = hesse(i,j);
+
}
TensorSSD<RT, embeddedDim,embeddedDim> dqdwF_times_dvq(coefficients_.size());
diff --git a/dune/gfe/realtuple.hh b/dune/gfe/realtuple.hh
index afd7396a..9fecb7e9 100644
--- a/dune/gfe/realtuple.hh
+++ b/dune/gfe/realtuple.hh
@@ -4,6 +4,7 @@
#include <dune/common/array.hh>
#include <dune/common/fvector.hh>
#include <dune/gfe/tensor3.hh>
+#include <dune/gfe/symmetricmatrix.hh>
/** \brief Implement a tuple of real numbers as a Riemannian manifold
@@ -101,12 +102,12 @@ public:
Unlike the distance itself the squared distance is differentiable at zero
*/
- static Dune::FieldMatrix<T,N,N> secondDerivativeOfDistanceSquaredWRTSecondArgument(const RealTuple& a, const RealTuple& b) {
+ static Dune::SymmetricMatrix<T,N> secondDerivativeOfDistanceSquaredWRTSecondArgument(const RealTuple& a, const RealTuple& b) {
- Dune::FieldMatrix<T,N,N> result;
+ Dune::SymmetricMatrix<T,N> result;
for (int i=0; i<N; i++)
- for (int j=0; j<N; j++)
- result[i][j] = 2*(i==j);
+ for (int j=0; j<=i; j++)
+ result(i,j) = 2*(i==j);
return result;
}
diff --git a/dune/gfe/rotation.hh b/dune/gfe/rotation.hh
index 24536405..3fc8f2d7 100644
--- a/dune/gfe/rotation.hh
+++ b/dune/gfe/rotation.hh
@@ -14,6 +14,7 @@
#include <dune/gfe/tensor3.hh>
#include <dune/gfe/unitvector.hh>
#include <dune/gfe/skewmatrix.hh>
+#include <dune/gfe/symmetricmatrix.hh>
template <class T, int dim>
class Rotation
@@ -684,24 +685,24 @@ public:
}
/** \brief Compute the Hessian of the squared distance function keeping the first argument fixed */
- static Dune::FieldMatrix<T,4,4> secondDerivativeOfDistanceSquaredWRTSecondArgument(const Rotation<T,3>& p, const Rotation<T,3>& q) {
+ static Dune::SymmetricMatrix<T,4> secondDerivativeOfDistanceSquaredWRTSecondArgument(const Rotation<T,3>& p, const Rotation<T,3>& q) {
T sp = p.globalCoordinates() * q.globalCoordinates();
EmbeddedTangentVector pProjected = q.projectOntoTangentSpace(p.globalCoordinates());
- Dune::FieldMatrix<T,4,4> A;
+ Dune::SymmetricMatrix<T,4> A;
for (int i=0; i<4; i++)
- for (int j=0; j<4; j++)
- A[i][j] = pProjected[i]*pProjected[j];
+ for (int j=0; j<=i; j++)
+ A(i,j) = pProjected[i]*pProjected[j];
A *= 4*UnitVector<T,4>::secondDerivativeOfArcCosSquared(std::abs(sp));
// Compute matrix B (see notes)
- Dune::FieldMatrix<T,4,4> Pq;
+ Dune::SymmetricMatrix<T,4> Pq;
for (int i=0; i<4; i++)
- for (int j=0; j<4; j++)
- Pq[i][j] = (i==j) - q.globalCoordinates()[i]*q.globalCoordinates()[j];
+ for (int j=0; j<=i; j++)
+ Pq(i,j) = (i==j) - q.globalCoordinates()[i]*q.globalCoordinates()[j];
// Bring it all together
A.axpy(-4* ((sp<0) ? -1 : 1) * UnitVector<T,4>::derivativeOfArcCosSquared(std::abs(sp))*sp, Pq);
diff --git a/dune/gfe/unitvector.hh b/dune/gfe/unitvector.hh
index a218ee40..f37fa4d9 100644
--- a/dune/gfe/unitvector.hh
+++ b/dune/gfe/unitvector.hh
@@ -5,6 +5,7 @@
#include <dune/common/fmatrix.hh>
#include <dune/gfe/tensor3.hh>
+#include <dune/gfe/symmetricmatrix.hh>
template <class T, int N>
class Rotation;
@@ -189,24 +190,24 @@ public:
Unlike the distance itself the squared distance is differentiable at zero
*/
- static Dune::FieldMatrix<T,N,N> secondDerivativeOfDistanceSquaredWRTSecondArgument(const UnitVector& p, const UnitVector& q) {
+ static Dune::SymmetricMatrix<T,N> secondDerivativeOfDistanceSquaredWRTSecondArgument(const UnitVector& p, const UnitVector& q) {
T sp = p.data_ * q.data_;
Dune::FieldVector<T,N> pProjected = q.projectOntoTangentSpace(p.globalCoordinates());
- Dune::FieldMatrix<T,N,N> A;
+ Dune::SymmetricMatrix<T,N> A;
for (int i=0; i<N; i++)
- for (int j=0; j<N; j++)
- A[i][j] = pProjected[i]*pProjected[j];
+ for (int j=0; j<=i; j++)
+ A(i,j) = pProjected[i]*pProjected[j];
A *= secondDerivativeOfArcCosSquared(sp);
// Compute matrix B (see notes)
- Dune::FieldMatrix<T,N,N> Pq;
+ Dune::SymmetricMatrix<T,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];
+ for (int j=0; j<=i; j++)
+ Pq(i,j) = (i==j) - q.data_[i]*q.data_[j];
// Bring it all together
A.axpy(-1*derivativeOfArcCosSquared(sp)*sp, Pq);
--
GitLab