diff --git a/dune/gfe/localgeodesicfefunction.hh b/dune/gfe/localgeodesicfefunction.hh index fe4261b7f2fd5cee4f3e439ba1073abcb51b3295..f987a13f4a3241f8ce8c6e0bc52a71a753a1ea50 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>