a default implementation of the energy functional gradient using finite differences

[[Imported from SVN: r4021]]

src/localgeodesicfestiffness.hh

@@ -87,13 +87,37 @@ public:
     
 };
 
-template <class GridType, class TargetSpace>
-void LocalGeodesicFEStiffness<GridType, TargetSpace>::
+template <class GridView, class TargetSpace>
+void LocalGeodesicFEStiffness<GridView, TargetSpace>::
 assembleGradient(const Entity& element,
-                 const std::vector<TargetSpace>& solution,
-                 Dune::array<Dune::FieldVector<double,6>, 2>& gradient) const
+                 const std::vector<TargetSpace>& localSolution,
+                 Dune::array<Dune::FieldVector<double,6>, 2>& localGradient) const
 {
- 
+    // ///////////////////////////////////////////////////////////
+    //   Compute gradient by finite-difference approximation
+    // ///////////////////////////////////////////////////////////
+
+    double eps = 1e-6;
+
+    std::vector<TargetSpace> forwardSolution = localSolution;
+    std::vector<TargetSpace> backwardSolution = localSolution;
+
+    for (size_t i=0; i<localSolution.size(); i++) {
+        
+        for (int j=0; j<6; j++) {
+            
+            infinitesimalVariation(forwardSolution[i],   eps, j);
+            infinitesimalVariation(backwardSolution[i], -eps, j);
+            
+            localGradient[i][j] = (energy(element,forwardSolution) - energy(element,backwardSolution))
+                / (2*eps);
+            
+            forwardSolution[i]  = localSolution[i];
+            backwardSolution[i] = localSolution[i];
+        }
+        
+    }
+
 }
 
 

src/rodlocalstiffness.hh

@@ -30,19 +30,6 @@ class RodLocalStiffness
     enum {bendingQuadOrder = 2};
 
 public:
-    /** \brief For the fd approximations 
-        \todo This is public because RodAssembler uses it
-    */
-    static void infinitesimalVariation(RigidBodyMotion<3>& c, double eps, int i)
-    {
-        if (i<3)
-            c.r[i] += eps;
-        else
-            c.q = c.q.mult(Rotation<3,double>::exp((i==3)*eps, 
-                                                   (i==4)*eps, 
-                                                   (i==5)*eps));
-    }
-    
     std::vector<RigidBodyMotion<3> > localReferenceConfiguration_;
 
 public: