fd approximation of the gradient in a dimension-independent way

[[Imported from SVN: r4072]]

src/localgeodesicfefunction.hh

@@ -259,46 +259,26 @@ template <int dim, class ctype, class TargetSpace>
 Dune::FieldMatrix<ctype, TargetSpace::EmbeddedTangentVector::size, dim> LocalGeodesicFEFunction<dim,ctype,TargetSpace>::
 evaluateDerivativeFD(const Dune::FieldVector<ctype, dim>& local)
 {
-    Dune::FieldMatrix<ctype, EmbeddedTangentVector::size, dim> result;
-
-    if (dim==1) {
-
-        EmbeddedTangentVector tmp = TargetSpace::interpolateDerivative(coefficients_[0], coefficients_[1], local[0]);
-
-        for (int i=0; i<EmbeddedTangentVector::size; i++)
-            result[i][0] = tmp[i];
-
-    }
-
-    if (dim==2) {
-
-        double eps = 1e-6;
+    double eps = 1e-6;
 
-        Dune::FieldVector<ctype, dim> fX = local;
-        Dune::FieldVector<ctype, dim> bX = local;
-        Dune::FieldVector<ctype, dim> fY = local;
-        Dune::FieldVector<ctype, dim> bY = local;
-
-        fX[0] += eps;
-        bX[0] -= eps;
-        fY[1] += eps;
-        bY[1] -= eps;
-
-        Quaternion<double> dX = (evaluate(fX) - evaluate(bX));
-        Quaternion<double> dY = (evaluate(fY) - evaluate(bY));
-
-        dX /= 2*eps;
-        dY /= 2*eps;
-
-        for (int i=0; i<EmbeddedTangentVector::size; i++){
-            result[i][0] = dX[i];
-            result[i][1] = dY[i];
-        }
+    Dune::FieldMatrix<ctype, EmbeddedTangentVector::size, dim> result;
 
+    for (int i=0; i<dim; i++) {
+        
+        Dune::FieldVector<ctype, dim> forward  = local;
+        Dune::FieldVector<ctype, dim> backward = local;
+        
+        forward[i]  += eps;
+        backward[i] -= eps;
+        
+        EmbeddedTangentVector fdDer = evaluate(forward) - evaluate(backward);
+        fdDer /= 2*eps;
+        
+        for (int j=0; j<EmbeddedTangentVector::size; j++)
+            result[j][i] = fdDer[j];
+        
     }
 
-    assert(dim==1 || dim==2);
-
     return result;
 }