rescale DirectionalRestriction to direction.two_norm = 1.0

dune/tectonic/spatial-solving/tnnmg/functional.hh

@@ -162,11 +162,12 @@ class DirectionalRestriction :
   using Matrix = typename Base::Matrix;
   using Vector = typename Base::Vector;
 
-  DirectionalRestriction(const GlobalMatrix& matrix, const GlobalVector& linearTerm, const Nonlinearity& phi, const GlobalLowerObstacle& lower, const GlobalLowerObstacle& upper, const GlobalVector& origin, const GlobalVector& direction) :
+  DirectionalRestriction(const GlobalMatrix& matrix, const GlobalVector& linearTerm, const Nonlinearity& phi, const GlobalLowerObstacle& lower, const GlobalLowerObstacle& upper, const GlobalVector& origin, const GlobalVector& direction, const double scaling) :
     Base(matrix, linearTerm, lower, upper, origin, direction),
     origin_(origin),
     direction_(direction),
-    phi_(phi)
+    phi_(phi),
+    scaling_(scaling)
   {
     phi_.directionalDomain(origin_, direction_, domain_);
 
@@ -207,12 +208,17 @@ class DirectionalRestriction :
         return direction_;
     }
 
+    double scaling() const {
+        return scaling_;
+    }
 protected:
   const GlobalVector& origin_;
-  const GlobalVector& direction_;
+  GlobalVector direction_;
 
   const Nonlinearity& phi_;
   Interval domain_;
+
+  const double scaling_;
 };
 
 
@@ -474,7 +480,18 @@ class Functional : public Dune::TNNMG::BoxConstrainedQuadraticFunctional<M, V, L
   friend auto directionalRestriction(const Functional& f, const Vector& origin, const Vector& direction)
     -> DirectionalRestriction<Matrix, Vector, Nonlinearity, LowerObstacle, UpperObstacle, Range>
   {
-    return DirectionalRestriction<Matrix, Vector, Nonlinearity, LowerObstacle, UpperObstacle, Range>(f.quadraticPart(), f.linearPart(), f.phi(), f.lowerObstacle(), f.upperObstacle(), origin, direction);
+    // rescale direction for numerical stability
+    auto dir = direction;
+    auto dirNorm = dir.two_norm();
+    if (dirNorm > 0.0)
+        dir /= dirNorm;
+    else {
+        dirNorm = 0.0;
+        dir = 0.0;
+    }
+
+
+    return DirectionalRestriction<Matrix, Vector, Nonlinearity, LowerObstacle, UpperObstacle, Range>(f.quadraticPart(), f.linearPart(), f.phi(), f.lowerObstacle(), f.upperObstacle(), origin, dir, dirNorm);
   }
 
   friend auto shift(const Functional& f, const Vector& origin)