coordinateRestriction now yields ShiftedFunctional

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

@@ -116,6 +116,22 @@ class ShiftedFunctional {
       return phi_;
   }
 
+  friend auto directionalRestriction(const ShiftedFunctional& f, const Vector& origin, const Vector& direction)
+    -> DirectionalRestriction<Matrix, Vector, Nonlinearity, LowerObstacle, UpperObstacle, Range>
+  {
+    // 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.originalLinearPart(), f.phi(), f.originalLowerObstacle(), f.originalUpperObstacle(), origin, dir, dirNorm);
+  }
+
 protected:
   Dune::Solvers::ConstCopyOrReference<M> quadraticPart_;
 
@@ -267,9 +283,9 @@ auto coordinateRestriction(const GlobalShiftedFunctional& f, const Index& i)
   dli -= origini;
   dui -= origini;
 
-  return Functional<decltype(Aii), LocalVector, decltype(phii), LocalLowerObstacle, LocalUpperObstacle, Range>(Aii, std::move(ri), phii, std::move(dli), std::move(dui));
+  //return Functional<decltype(Aii), LocalVector, decltype(phii), LocalLowerObstacle, LocalUpperObstacle, Range>(Aii, std::move(ri), phii, std::move(dli), std::move(dui));
 
-  //return ShiftedFunctional<decltype(Aii), LocalVector, decltype(phii), LocalLowerObstacle, LocalUpperObstacle, LocalVector, Range>(Aii, std::move(ri), phii, std::move(dli), std::move(dui), origini);
+  return ShiftedFunctional<decltype(Aii), LocalVector, decltype(phii), LocalLowerObstacle, LocalUpperObstacle, LocalVector, Range>(Aii, std::move(ri), phii, std::move(dli), std::move(dui), origini);
 }
 
 

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

@@ -26,9 +26,8 @@ class LocalBisectionSolver
           return;
       }
 
-      auto old_x = x;
-      auto origin = x;
-      auto linearPart = f.linearPart();
+      auto origin = f.origin();
+      auto linearPart = f.originalLinearPart();
 
 
       //std::cout << x << std::endl;
@@ -36,11 +35,11 @@ class LocalBisectionSolver
       //std::cout << "ignore: " << ignore << std::endl;
       //std::cout << "orig linear part: " << linearPart << std::endl;
 
-      auto lower = f.lowerObstacle();
-      auto upper = f.upperObstacle();
+      auto lower = f.originalLowerObstacle();
+      auto upper = f.originalUpperObstacle();
       for (size_t j = 0; j < lower.size(); ++j) {
         if (std::abs(upper[j]-lower[j])<= safety) {
-          //linearPart[j] = 0;
+          linearPart[j] = 0;
         }
       }
 
@@ -54,13 +53,14 @@ class LocalBisectionSolver
         }
       }
 
+      /*
       auto linearPartNorm = linearPart.two_norm();
       if (linearPartNorm > 0.0)
           linearPart /= linearPartNorm;
       else {
           linearPartNorm = 0.0;
           linearPart = 0.0;
-      }
+      }*/
 
       //using DirectionalRestriction = DirectionalRestriction<Matrix, Vector, Nonlinearity, LowerObstacle, UpperObstacle, Range>;
       //DirectionalRestriction fv(f.quadraticPart(), f.linearPart(), f.phi(), f.lowerObstacle(), f.upperObstacle(), f.origin(), linearPart, linearPartNorm);
@@ -70,10 +70,10 @@ class LocalBisectionSolver
       //std::cout << "origin: " << origin << std::endl;
       //std::cout << " linearPart: " << linearPart << std::endl;
 
-      std::cout << "maxEigenvalue " << f.quadraticPart() << std::endl;
-      print(linearPart, "linearPart:");
-      std::cout << "linearPartNorm " << linearPartNorm << std::endl;
-      print(origin, "origin:");
+      //std::cout << "maxEigenvalue " << f.quadraticPart() << std::endl;
+      //print(linearPart, "linearPart:");
+      //std::cout << "linearPartNorm " << linearPartNorm << std::endl;
+      //print(origin, "origin:");
 
 
 
@@ -115,11 +115,11 @@ class LocalBisectionSolver
 
       if (fv.scaling()<=0.0) {
         alpha = 0.0;
-      } /* else {
+      } else {
         LineSearchSolver lineSearch;
         lineSearch(alpha, fv, false);
-      }*/
-
+      }
+/*
       const auto& domain = fv.domain();
       if (std::abs(domain[0]-domain[1])>safety) {
         int bisectionsteps = 0;
@@ -128,15 +128,16 @@ class LocalBisectionSolver
         alpha = globalBisection.minimize(fv, 0.0, 0.0, bisectionsteps);
       } else {
           alpha = domain[0];
-      }
+      }*/
 
 
-      std::cout << "alpha: " << alpha << std::endl;
+      //std::cout << "alpha: " << alpha << std::endl;
 
       linearPart *= alpha;
 
-      std::cout << linearPart << std::endl;
+      //std::cout << linearPart << std::endl;
 
+      //std::cout << "old_x: " << old_x << std::endl;
 #ifdef NEW_TNNMG_COMPUTE_ITERATES_DIRECTLY
       if (std::abs(alpha)> safety) {
         x = origin;
@@ -145,14 +146,14 @@ class LocalBisectionSolver
         x = origin;
       }
 #else
-      std::cout << "use old" << std::endl;
+      //std::cout << "use old" << std::endl;
       if (std::abs(alpha)> safety) {
         x = origin;
         x += linearPart;
-        x -= old_x;//f.origin();
+        x -= f.origin();
       }
 #endif
-      std::cout << "new x: "<< x << std::endl << std::endl;
+      //std::cout << "new x: "<< x << std::endl << std::endl;
   }
 };