From cdab8df6f98478051eb7fe0130881c92890bb01c Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Sun, 28 Mar 2010 20:14:53 +0000
Subject: [PATCH] bugfix: dFdq may not be invertible
[[Imported from SVN: r5804]]
---
src/localgeodesicfefunction.hh | 40 ++++++++++++++++++++++++++++++++--
1 file changed, 38 insertions(+), 2 deletions(-)
diff --git a/src/localgeodesicfefunction.hh b/src/localgeodesicfefunction.hh
index ea167299..b89f645a 100644
--- a/src/localgeodesicfefunction.hh
+++ b/src/localgeodesicfefunction.hh
@@ -9,6 +9,8 @@
#include <dune/src/averagedistanceassembler.hh>
#include <dune/src/targetspacertrsolver.hh>
+#include <dune/src/svd.hh>
+
/** \brief A geodesic function from the reference element to a manifold
\tparam dim Dimension of the reference element
@@ -20,6 +22,7 @@ class LocalGeodesicFEFunction
{
typedef typename TargetSpace::EmbeddedTangentVector EmbeddedTangentVector;
+ static const int targetDim = EmbeddedTangentVector::size;
public:
@@ -48,7 +51,34 @@ private:
}
return result;
}
+
+ static Dune::FieldMatrix<double,targetDim,targetDim> pseudoInverse(const Dune::FieldMatrix<double,targetDim,targetDim>& A)
+ {
+ Dune::FieldMatrix<double,targetDim,targetDim> U = A;
+ Dune::FieldVector<double,targetDim> W;
+ Dune::FieldMatrix<double,targetDim,targetDim> V;
+
+ svdcmp(U,W,V);
+
+ // pseudoInv = V W^{-1} U^T
+ Dune::FieldMatrix<double,targetDim,targetDim> UT;
+
+ for (int i=0; i<targetDim; i++)
+ for (int j=0; j<targetDim; j++)
+ UT[i][j] = U[j][i];
+
+ for (int i=0; i<targetDim; i++) {
+ if (std::abs(W[i]) > 1e-12) // Diagonal may be zero, that's why we're using the pseudo inverse
+ UT[i] /= W[i];
+ else
+ UT[i] = 0;
+ }
+
+ Dune::FieldMatrix<double,targetDim,targetDim> pseudoInv;
+ Dune::FMatrixHelp::multMatrix(V,UT,pseudoInv);
+ return pseudoInv;
+ }
/** \brief The coefficient vector */
std::vector<TargetSpace> coefficients_;
@@ -104,7 +134,7 @@ evaluateDerivative(const Dune::FieldVector<ctype, dim>& local)
// Hence we need to solve a small system of linear equations.
// ////////////////////////////////////////////////////////////////////////
- // the function value at the point where we are evaluation the derivative
+ // the function value at the point where we are evaluating the derivative
TargetSpace q = evaluate(local);
// the matrix that turns coordinates on the reference simplex into coordinates on the standard simplex
@@ -139,13 +169,19 @@ evaluateDerivative(const Dune::FieldVector<ctype, dim>& local)
// solve the system
// ////////////////////////////////////
+ // 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,targetDim,targetDim> dFdqPseudoInv = pseudoInverse(dFdq);
+
for (int i=0; i<dim; i++) {
Dune::FieldVector<ctype,targetDim> rhs, x;
for (int j=0; j<targetDim; j++)
rhs[j] = RHS[j][i];
- dFdq.solve(x, rhs);
+ //dFdq.solve(x, rhs);
+ dFdqPseudoInv.mv(rhs,x);
for (int j=0; j<targetDim; j++)
result[j][i] = x[j];
--
GitLab