diff --git a/dune/tectonic/myblockproblem.hh b/dune/tectonic/myblockproblem.hh index 2d7028796c94832b9dee18e913bb8c2ef14e23ba..de1d51304a7515b90f8c66bd33bbaa874e7560b1 100644 --- a/dune/tectonic/myblockproblem.hh +++ b/dune/tectonic/myblockproblem.hh @@ -187,6 +187,25 @@ template <class MyConvexProblemTypeTEMPLATE> class MyBlockProblem { // -grad is needed for Newton step linearization.b *= -1.0; + // b should be a descent direction + { + VectorType const direction = linearization.b; + VectorType tmp = linearization.b; // b + linearization.A.mmv(u, tmp); // b-Au + double const localA = tmp * direction; // <b-Au,v> + + linearization.A.mv(direction, tmp); // Av + double const localb = tmp * direction; // <Av,v> + + MyDirectionalConvexFunction<Dune::GlobalNonlinearity<block_size>> const + psi(localA, localb, problem.phi, u, direction); + + Interval<double> D; + psi.subDiff(0, D); + if (!isnan(D[1])) + assert(D[1] <= 0); + } + // apply truncation to system for (size_t row = 0; row < linearization.A.N(); ++row) { auto const end = linearization.A[row].end();