[cleanup] Move the FD approximation of the GFE function gradient from the GFE...

[cleanup] Move the FD approximation of the GFE function gradient from the GFE function into the test

It is a part of the test harness, and should not clutter the implementation of the
GFE functions themselves.

dune/gfe/localgeodesicfefunction.hh

@@ -84,9 +84,6 @@ public:
     DerivativeType evaluateDerivative(const Dune::FieldVector<ctype, dim>& local,
                                       const TargetSpace& q) const;
 
-    /** \brief For debugging: Evaluate the derivative of the function using a finite-difference approximation*/
-    DerivativeType evaluateDerivativeFD(const Dune::FieldVector<ctype, dim>& local) const;
-
     /** \brief Evaluate the derivative of the function value with respect to a coefficient */
     void evaluateDerivativeOfValueWRTCoefficient(const Dune::FieldVector<ctype, dim>& local,
                                                  int coefficient,
@@ -301,34 +298,6 @@ evaluateDerivative(const Dune::FieldVector<ctype, dim>& local, const TargetSpace
     return result;
 }
 
-template <int dim, class ctype, class LocalFiniteElement, class TargetSpace>
-typename LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::DerivativeType
-LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::
-evaluateDerivativeFD(const Dune::FieldVector<ctype, dim>& local) const
-{
-    double eps = 1e-6;
-
-    Dune::FieldMatrix<ctype, embeddedDim, 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).globalCoordinates() - evaluate(backward).globalCoordinates();
-        fdDer /= 2*eps;
-
-        for (int j=0; j<embeddedDim; j++)
-            result[j][i] = fdDer[j];
-
-    }
-
-    return result;
-}
-
 template <int dim, class ctype, class LocalFiniteElement, class TargetSpace>
 void LocalGeodesicFEFunction<dim,ctype,LocalFiniteElement,TargetSpace>::
 evaluateDerivativeOfValueWRTCoefficient(const Dune::FieldVector<ctype, dim>& local,
@@ -699,28 +668,6 @@ public:
         return result;
     }
 
-    /** \brief For debugging: Evaluate the derivative of the function using a finite-difference approximation*/
-    DerivativeType evaluateDerivativeFD(const Dune::FieldVector<ctype, dim>& local) const
-    {
-        Dune::FieldMatrix<ctype, embeddedDim, dim> result(0);
-
-        // get translation part
-        std::vector<Dune::FieldMatrix<ctype,1,dim> > sfDer(translationCoefficients_.size());
-        localFiniteElement_.localBasis().evaluateJacobian(local, sfDer);
-
-        for (size_t i=0; i<translationCoefficients_.size(); i++)
-            for (int j=0; j<3; j++)
-                result[j].axpy(translationCoefficients_[i][j], sfDer[i][0]);
-
-        // get orientation part
-        Dune::FieldMatrix<ctype,4,dim> qResult = orientationFEFunction_->evaluateDerivativeFD(local);
-        for (int i=0; i<4; i++)
-            for (int j=0; j<dim; j++)
-                result[3+i][j] = qResult[i][j];
-
-        return result;
-    }
-
     /** \brief Evaluate the derivative of the function value with respect to a coefficient */
     void evaluateDerivativeOfValueWRTCoefficient(const Dune::FieldVector<ctype, dim>& local,
                                                  int coefficient,

test/localgeodesicfefunctiontest.cc

@@ -34,6 +34,34 @@ double diameter(const std::vector<TargetSpace>& v)
     return d;
 }
 
+template <int dim, class ctype, class LocalFunction>
+auto
+evaluateDerivativeFD(const LocalFunction& f, const Dune::FieldVector<ctype, dim>& local)
+-> decltype(f.evaluateDerivative(local))
+{
+    double eps = 1e-8;
+    //static const int embeddedDim = LocalFunction::TargetSpace::embeddedDim;
+    decltype(f.evaluateDerivative(local)) 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;
+
+        auto fdDer = f.evaluate(forward).globalCoordinates() - f.evaluate(backward).globalCoordinates();
+        fdDer /= 2*eps;
+
+        for (size_t j=0; j<result.N(); j++)
+            result[j][i] = fdDer[j];
+
+    }
+
+    return result;
+}
+
 
 template <int domainDim>
 void testDerivativeTangentiality(const RealTuple<double,1>& x,
@@ -155,7 +183,7 @@ void testDerivative(const LocalGeodesicFEFunction<domainDim,double,typename PQkL
         Dune::FieldMatrix<double, embeddedDim, domainDim> derivative = f.evaluateDerivative(quadPos);
 
         // evaluate fd approximation of derivative
-        Dune::FieldMatrix<double, embeddedDim, domainDim> fdDerivative = f.evaluateDerivativeFD(quadPos);
+        Dune::FieldMatrix<double, embeddedDim, domainDim> fdDerivative = evaluateDerivativeFD(f,quadPos);
 
         Dune::FieldMatrix<double, embeddedDim, domainDim> diff = derivative;
         diff -= fdDerivative;