From f2726eb998cfccd8d5162b1bd8cb9da0d4c11c21 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Tue, 29 Mar 2011 20:46:20 +0000
Subject: [PATCH] Use coordinates in the tangent spaces to evaluate gfe
 functions

... instead of coordinates in the embedding space.  The usual problem
is that when working in the embedding space, a projection onto the
trust-region makes the iterate leave the tangent space.

[[Imported from SVN: r7035]]
---
 dune/gfe/averagedistanceassembler.hh | 23 +++++++++++++++++++++--
 dune/gfe/targetspacertrsolver.hh     |  2 +-
 2 files changed, 22 insertions(+), 3 deletions(-)

diff --git a/dune/gfe/averagedistanceassembler.hh b/dune/gfe/averagedistanceassembler.hh
index ff09c61e..92914068 100644
--- a/dune/gfe/averagedistanceassembler.hh
+++ b/dune/gfe/averagedistanceassembler.hh
@@ -8,7 +8,8 @@
 template <class TargetSpace>
 class AverageDistanceAssembler
 {
-    static const int size = TargetSpace::EmbeddedTangentVector::size;
+    static const int size         = TargetSpace::TangentVector::size;
+    static const int embeddedSize = TargetSpace::EmbeddedTangentVector::size;
 
 public:
 
@@ -38,6 +39,24 @@ public:
                           TargetSpace::derivativeOfDistanceSquaredWRTSecondArgument(coefficients_[i], x));
     }
 
+    void assembleGradient(const TargetSpace& x,
+                          typename TargetSpace::TangentVector& gradient) const
+    {
+        typename TargetSpace::EmbeddedTangentVector embeddedGradient;
+        assembleGradient(x,embeddedGradient);
+        
+        Dune::FieldMatrix<double,size,embeddedSize> orthonormalFrame = x.orthonormalFrame();
+        orthonormalFrame.mv(embeddedGradient,gradient);
+    }
+
+    void assembleHessianApproximation(const TargetSpace& x,
+                                      Dune::FieldMatrix<double,embeddedSize,embeddedSize>& matrix) const
+    {
+        for (int i=0; i<embeddedSize; i++)
+            for (int j=0; j<embeddedSize; j++)
+                matrix[i][j] = (i==j);
+    }
+
     void assembleHessianApproximation(const TargetSpace& x,
                                       Dune::FieldMatrix<double,size,size>& matrix) const
     {
@@ -47,7 +66,7 @@ public:
     }
 
     void assembleHessian(const TargetSpace& x,
-                         Dune::FieldMatrix<double,size,size>& matrix) const
+                         Dune::FieldMatrix<double,embeddedSize,embeddedSize>& matrix) const
     {
         matrix = 0;
         for (size_t i=0; i<coefficients_.size(); i++)
diff --git a/dune/gfe/targetspacertrsolver.hh b/dune/gfe/targetspacertrsolver.hh
index c05aadea..cb6928b4 100644
--- a/dune/gfe/targetspacertrsolver.hh
+++ b/dune/gfe/targetspacertrsolver.hh
@@ -13,7 +13,7 @@ class TargetSpaceRiemannianTRSolver
 //     : public IterativeSolver<std::vector<TargetSpace>,
 //                             Dune::BitSetVector<TargetSpace::TangentVector::size> >
 { 
-    const static int blocksize = TargetSpace::EmbeddedTangentVector::size;
+    const static int blocksize = TargetSpace::TangentVector::size;
 
     // Centralize the field type here
     typedef double field_type;
-- 
GitLab