Fix up CircularCurvedFunction::cartesian()

dune/tectonic/circularconvexfunction.hh

@@ -13,7 +13,13 @@ template <class NonlinearityType> class CircularConvexFunction {
   CircularConvexFunction(MatrixType const &A, VectorType const &b,
                          NonlinearityType const &phi, VectorType const &x,
                          VectorType const &dir)
-      : A(A), b(b), phi(phi), x(x), dir(dir) {}
+      : A(A),
+        b(b),
+        phi(phi),
+        x(x),
+        dir(dir),
+        xnorm(x.two_norm()),
+        dnorm(dir.two_norm()) {}
 
   double quadraticPart() const { return 0; }
 
@@ -43,7 +49,7 @@ template <class NonlinearityType> class CircularConvexFunction {
   void cartesian(double m, VectorType &y) const {
     y = 0;
     y.axpy(std::cos(m), x);
-    y.axpy(std::sin(m), dir);
+    y.axpy(std::sin(m) * xnorm / dnorm, dir);
   }
 
 private:
@@ -53,6 +59,9 @@ template <class NonlinearityType> class CircularConvexFunction {
   VectorType const &x;
   VectorType const &dir;
 
+  double const dnorm;
+  double const xnorm;
+
   /* If x and d were normalised, we would have
 
        cartesian(a) =  cos(a) * x + sin(a) * d and
@@ -63,8 +72,8 @@ template <class NonlinearityType> class CircularConvexFunction {
   */
   void tangent(double m, VectorType &y) const {
     y = 0;
-    y.axpy(-std::sin(m) * dir.two_norm2(), x);
-    y.axpy(std::cos(m) * x.two_norm2(), dir);
+    y.axpy(-std::sin(m) * dnorm, x);
+    y.axpy(std::cos(m) * xnorm, dir);
   }
 };
 }

src/test-gradient-sample-steep.cc

@@ -36,18 +36,18 @@ int main() {
   start[0] = 0;
   start[1] = 1;
 
-  double const ret1 = functionTester(J, start, 5);
+  double const ret1 = functionTester(J, start, 3);
 
   // Something random
   start[0] = 279;
   start[1] = -96;
 
-  double const ret2 = functionTester(J, start, 5);
-  assert(std::abs(ret1 - ret2) < 1e-7);
+  double const ret2 = functionTester(J, start, 3);
+  assert(std::abs(ret1 - ret2) < 1e-5);
 
   start[0] = 0;
   start[1] = 0;
 
   double const ret3 = functionTester(J, start, 1);
-  assert(std::abs(ret1 - ret3) < 1e-7);
+  assert(std::abs(ret1 - ret3) < 1e-5);
 }

src/test-gradient-sample-steep2.cc

@@ -42,12 +42,12 @@ int main() {
   start[0] = 279;
   start[1] = -96;
 
-  double const ret2 = functionTester(J, start, 5);
-  assert(std::abs(ret1 - ret2) < 1e-7);
+  double const ret2 = functionTester(J, start, 3);
+  assert(std::abs(ret1 - ret2) < 1e-5);
 
   start[0] = 0;
   start[1] = 0;
 
-  double const ret3 = functionTester(J, start, 2);
-  assert(std::abs(ret1 - ret3) < 1e-7);
+  double const ret3 = functionTester(J, start, 1);
+  assert(std::abs(ret1 - ret3) < 1e-5);
 }

src/test-gradient-sample-verysteep.cc

@@ -42,12 +42,12 @@ int main() {
   start[0] = 279;
   start[1] = -96;
 
-  double const ret2 = functionTester(J, start, 5);
-  assert(std::abs(ret1 - ret2) < 1e-7);
+  double const ret2 = functionTester(J, start, 4);
+  assert(std::abs(ret1 - ret2) < 1e-8);
 
   start[0] = 0;
   start[1] = 0;
 
-  double const ret3 = functionTester(J, start, 2);
-  assert(std::abs(ret1 - ret3) < 1e-7);
+  double const ret3 = functionTester(J, start, 1);
+  assert(std::abs(ret1 - ret3) < 1e-5);
 }