diff --git a/dune/tectonic/spatial-solving/tnnmg/linearization.hh b/dune/tectonic/spatial-solving/tnnmg/linearization.hh index 8abd074b094c20caf6070a44fa01d93c63ba7bbb..68cbd4176c419169f5617fe0cab1e460b3502fb1 100644 --- a/dune/tectonic/spatial-solving/tnnmg/linearization.hh +++ b/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_;