directionalSubDiff works on BlockVectors!

src/bisection-example-new.cc

@@ -62,22 +62,28 @@ class SampleFunctional {
        since A is symmetric.
      */
 
-    SmallVector tmp2;
-    A_.mv(descDir, tmp2);
-    double const rest_A = tmp2 * descDir;
-
-    SmallVector tmp3;
-    A_.mv(x, tmp3);
-    double const rest_b = (b_ - tmp3) * descDir;
-
-    Dune::BlockVector<SmallVector> xx; // FIXME: we actually want x here
-    Dune::BlockVector<SmallVector> dd; // FIXME: we actually want descDir here
-
-    typedef MyNonlinearity<dimension, SampleFunction> MyNonlinearityType;
-    MyNonlinearityType phi;
-    typedef DirectionalConvexFunction<MyNonlinearityType>
-    MyDirectionalConvexFunctionType;
-    MyDirectionalConvexFunctionType rest(rest_A, rest_b, phi, xx, dd);
+    {
+      SmallVector tmp2;
+      A_.mv(descDir, tmp2);
+      double const rest_A = tmp2 * descDir;
+
+      SmallVector tmp3;
+      A_.mv(x, tmp3);
+      double const rest_b = (b_ - tmp3) * descDir;
+
+      Dune::BlockVector<SmallVector> xx; // FIXME: we actually want x here
+      Dune::BlockVector<SmallVector> dd; // FIXME: we actually want descDir here
+
+      typedef MyNonlinearity<dimension, SampleFunction> MyNonlinearityType;
+      MyNonlinearityType phi;
+      typedef DirectionalConvexFunction<MyNonlinearityType>
+      MyDirectionalConvexFunctionType;
+      MyDirectionalConvexFunctionType rest(rest_A, rest_b, phi, xx, dd);
+
+      // Experiment a bit
+      Interval<double> D;
+      rest.subDiff(0, D);
+    }
 
     if (descDir == SmallVector(0.0))
       return SmallVector(0.0);

src/mynonlinearity.cc

@@ -24,16 +24,16 @@ class MyNonlinearity {
   typedef Dune::BlockVector<SmallVector> VectorType;
   typedef Dune::BlockVector<SmallMatrix> MatrixType;
 
-  void directionalSubDiff(SmallVector u, SmallVector v, Interval<double>& D) {
-    if (u == SmallVector(0.0)) {
-      D[0] = D[1] = func_.rightDifferential(0);
-    } else if (u * v > 0) {
-      D[1] = (v * u) * func_.rightDifferential(u.two_norm()) / u.two_norm();
-      D[0] = (v * u) * func_.leftDifferential(u.two_norm()) / u.two_norm();
-    } else {
-      D[1] = (v * u) * func_.leftDifferential(u.two_norm()) / u.two_norm();
-      D[0] = (v * u) * func_.rightDifferential(u.two_norm()) / u.two_norm();
-    }
+  void directionalSubDiff(VectorType u, VectorType v, Interval<double>& D) {
+    // if (u == SmallVector(0.0)) {
+    //   D[0] = D[1] = func_.rightDifferential(0);
+    // } else if (u * v > 0) {
+    //   D[1] = (v * u) * func_.rightDifferential(u.two_norm())/u.two_norm();
+    //   D[0] = (v * u) * func_.leftDifferential(u.two_norm())/u.two_norm();
+    // } else {
+    //   D[1] = (v * u) * func_.leftDifferential(u.two_norm())/u.two_norm();
+    //   D[0] = (v * u) * func_.rightDifferential(u.two_norm())/u.two_norm();
+    // }
   }
 
   void directionalDomain(const VectorType&, const VectorType&,