From 9ea062b0d4120306e12a400ccabba31deb2eb565 Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Thu, 18 Feb 2010 17:17:06 +0000
Subject: [PATCH] implement fd approximation of the gradient for
 unit-vector-valued functions

[[Imported from SVN: r5590]]
---
 src/localgeodesicfestiffness.hh | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/src/localgeodesicfestiffness.hh b/src/localgeodesicfestiffness.hh
index b1862fa3..e2b16b7d 100644
--- a/src/localgeodesicfestiffness.hh
+++ b/src/localgeodesicfestiffness.hh
@@ -304,7 +304,7 @@ class LocalGeodesicFEStiffness <GridView,UnitVector<dim> >
 
     /** \brief For the fd approximations 
     */
-    static Dune::FieldVector<double,dim> infinitesimalVariation(UnitVector<dim>& c, double eps, int i)
+    static Dune::FieldVector<double,dim> infinitesimalVariation(const UnitVector<dim>& c, double eps, int i)
     {
         Dune::FieldVector<double,dim> result = c.globalCoordinates();
         result[i] += eps;
@@ -386,12 +386,13 @@ assembleGradient(const Entity& element,
         
         for (int j=0; j<blocksize; j++) {
             
-            infinitesimalVariation(forwardSolution[i],   eps, j);
-            infinitesimalVariation(backwardSolution[i], -eps, j);
-            
+            // Brute force: the return value does not have unit norm.  Stuff it in there anyways
+            forwardSolution[i].data_  = infinitesimalVariation(localSolution[i],  eps, j);
+            backwardSolution[i].data_ = infinitesimalVariation(localSolution[i], -eps, j);
+
             localGradient[i][j] = (energy(element,forwardSolution) - energy(element,backwardSolution))
                 / (2*eps);
-            
+
             forwardSolution[i]  = localSolution[i];
             backwardSolution[i] = localSolution[i];
         }
-- 
GitLab