Implement test for the directionalSubDiff method

dune/tnnmg/functionals/test/functionaltest.hh

@@ -245,7 +245,39 @@ template <class Functional>
 void testDirectionalSubdifferential(const Functional& functional,
                                     const std::vector<typename Functional::VectorType>& testPoints)
 {
-  // TODO: Implement me!
+  // Step size. Best value: square root of the machine precision
+  const double eps = std::sqrt(std::numeric_limits<double>::epsilon());
+
+  // Loop over all test points
+  for (auto&& testPoint : testPoints)
+  {
+    // Abuse test points also as test directions
+    for (auto&& testDirection : testPoints)
+    {
+      Solvers::Interval<double> subDifferential;
+      functional.directionalSubDiff(testPoint, testDirection, subDifferential);
+
+      auto forwardPoint = testPoint;
+      forwardPoint.axpy(eps,testDirection);
+      auto backwardPoint = testPoint;
+      backwardPoint.axpy(-eps,testDirection);
+
+      auto forwardFDGradient  = (functional(forwardPoint) - functional(testPoint)) / eps;
+      auto backwardFDGradient = (functional(testPoint) - functional(backwardPoint)) / eps;
+
+      if (std::abs(backwardFDGradient - subDifferential[0]) > 100*eps
+          or std::abs(forwardFDGradient - subDifferential[1]) > 100*eps)
+      {
+        std::cout << "Bug in directionalSubDiff for functional type " << Dune::className<Functional>() << ":" << std::endl;
+        std::cerr << "Subdifferential doesn't match FD approximation" << std::endl;
+        std::cerr << "SubDifferential: " << subDifferential
+                  << ",   FD: [" << backwardFDGradient << ", " << forwardFDGradient << "]" << std::endl;
+        std::cerr << "Test point:" << std::endl << testPoint << std::endl;
+        std::cerr << "Test direction:" << std::endl << testDirection << std::endl;
+        DUNE_THROW(MathError,"");
+      }
+    }
+  }
 }