From 032e25622713a7af6e0d975c2a0fb08112076130 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Mon, 30 May 2011 08:29:38 +0000
Subject: [PATCH] avoid the use of the pseudo inverse in
evaluateDerivativeOfValueWRTCoefficient
[[Imported from SVN: r7345]]
---
dune/gfe/localgeodesicfefunction.hh | 35 +++++++++++++++++++++++------
1 file changed, 28 insertions(+), 7 deletions(-)
diff --git a/dune/gfe/localgeodesicfefunction.hh b/dune/gfe/localgeodesicfefunction.hh
index fe4261b7..f987a13f 100644
--- a/dune/gfe/localgeodesicfefunction.hh
+++ b/dune/gfe/localgeodesicfefunction.hh
@@ -331,16 +331,37 @@ evaluateDerivativeOfValueWRTCoefficient(const Dune::FieldVector<ctype, dim>& loc
Dune::FieldMatrix<ctype,embeddedDim,embeddedDim> dFdq(0);
assembler.assembleEmbeddedHessian(q,dFdq);
- // dFDq is not invertible, if the target space is embedded into a higher-dimensional
- // Euclidean space. Therefore we use its pseudo inverse. I don't think that is the
- // best way, though.
- Dune::FieldMatrix<ctype,embeddedDim,embeddedDim> dFdqPseudoInv = pseudoInverse(dFdq);
+ const int shortDim = TargetSpace::TangentVector::size;
+
+ // the orthonormal frame
+ Dune::FieldMatrix<ctype,shortDim,embeddedDim> O = q.orthonormalFrame();
+
+ // compute A = O dFDq O^T
+ Dune::FieldMatrix<ctype,shortDim,shortDim> A;
+ for (int i=0; i<shortDim; i++)
+ for (int j=0; j<shortDim; j++) {
+ A[i][j] = 0;
+ for (int k=0; k<embeddedDim; k++)
+ for (int l=0; l<embeddedDim; l++)
+ A[i][j] += O[i][k]*dFdq[k][l]*O[j][l];
+ }
//
- result = TargetSpace::secondDerivativeOfDistanceSquaredWRTFirstAndSecondArgument(coefficients_[coefficient], q);
- result *= -w[coefficient];
+ Dune::FieldMatrix<double,embeddedDim,embeddedDim> rhs = TargetSpace::secondDerivativeOfDistanceSquaredWRTFirstAndSecondArgument(coefficients_[coefficient], q);
+ rhs *= -w[coefficient];
+
+ for (int i=0; i<embeddedDim; i++) {
+
+ Dune::FieldVector<ctype,shortDim> shortRhs;
+ O.mv(rhs[i],shortRhs);
- result = result * dFdqPseudoInv;
+ Dune::FieldVector<ctype,shortDim> shortX;
+ A.solve(shortX,shortRhs);
+
+ O.mtv(shortX,result[i]);
+
+ }
+
}
template <int dim, class ctype, class TargetSpace>
--
GitLab