fix regularity truncation, only truncates tangential dofs

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

@@ -34,6 +34,9 @@ class Linearization {
   void determineRegularityTruncation(const Vector& x, BitVector& truncationFlags, const T& truncationTolerance)
   {
     //std::cout << "x: " << x << std::endl;
+    size_t blocksize = truncationFlags[0].size();
+
+    size_t count = 0;
 
     for (size_t i = 0; i < x.size(); ++i) {
       //if (truncationFlags[i].all())
@@ -42,9 +45,14 @@ class Linearization {
       //std::cout << f_.phi().regularity(i, x[i]) << " xnorm: " << x[i].two_norm() << std::endl;
 
       if (f_.phi().regularity(i, x[i]) > truncationTolerance) {
-        truncationFlags[i] = true;
+        for (size_t j=1; j<blocksize; j++) {
+            truncationFlags[i][j] = true;
+        }
+        count++;
       }
     }
+
+     std::cout << "regularityTruncation: " << count << std::endl;
   }
 
   template<class NV, class NBV, class T>
@@ -52,8 +60,9 @@ class Linearization {
   {
     namespace H = Dune::Hybrid;
     H::ifElse(Dune::IsNumber<NV>(), [&](auto id){
-      if ((id(x) <= id(lower)+truncationTolerance) || (id(x) >= id(upper) - truncationTolerance))
+      if ((id(x) <= id(lower)+truncationTolerance) || (id(x) >= id(upper) - truncationTolerance)) {
         id(truncationFlags) = true;
+      }
     }, [&](auto id){
       H::forEach(H::integralRange(H::size(id(x))), [&](auto&& i) {
         This::determineTruncation(x[i], lower[i], upper[i], truncationFlags[i], truncationTolerance);
@@ -160,7 +169,9 @@ class Linearization {
       //print(hessian_, "Hessian: ");
       //std::cout << "--------" << (double) hessian_[21][21] << "------------" << std::endl;
 
-      //auto B = hessian_;
+      //std::cout << hessian_[0][0] << std::endl;
+
+      auto B = hessian_;
 
       //print(x, "x: ");
 
@@ -168,7 +179,9 @@ class Linearization {
       phi.addHessian(x, hessian_);
       //truncateMatrix(hessian_, regularityTruncationFlags, regularityTruncationFlags);
 
-      //B -= hessian_;
+      //std::cout << x[0] << " " << hessian_[0][0] << std::endl;
+
+      B -= hessian_;
       //print(B, "phi hessian: ");
 
       // compute gradient
@@ -206,9 +219,9 @@ class Linearization {
     truncateVector(v, truncationFlags_);
   }
 
-  const BitVector& truncated() const {
+  /*const BitVector& truncated() const {
     return truncationFlags_;
-  }
+  }*/
 
   const auto& negativeGradient() const {
     return negativeGradient_;