diff --git a/dune/contact/solvers/filtermultigridsolver.hh b/dune/contact/solvers/filtermultigridsolver.hh index 6cf919026ac0c85bf733bba1ece771fd58a71c5d..3aba5cd28618fa8c8e993703de7ca557a41c60be 100644 --- a/dune/contact/solvers/filtermultigridsolver.hh +++ b/dune/contact/solvers/filtermultigridsolver.hh @@ -3,7 +3,6 @@ #ifndef DUNE_CONTACT_SOLVERS_FILTER_MULTIGRID_CONTACT_SOLVER_HH #define DUNE_CONTACT_SOLVERS_FILTER_MULTIGRID_CONTACT_SOLVER_HH - #include <dune/common/parametertree.hh> #include <dune/common/bitsetvector.hh> @@ -182,6 +181,29 @@ protected: //! Compute criticality measure for the current sub-problem field_type criticality() const { + + if (dirichletValues_.infinity_norm()>1e-10) + return dirichletValues_.infinity_norm(); + + const auto& obstacles = problem_->constraints(); + const auto& f = problem_->f(); + + VectorType x(f.size()); + for (size_t i=0; i<f.size();i++) + for (int k=0; k<dim; k++) + if (dirichletNodes_[i][k]) + x[i][k] = 0; + else if (f[i][k]>0) + x[i][k] = std::min(1.0, obstacles[i].upper(k)); + else if (f[i][k]<0) + x[i][k] = -1; + return std::fabs(f*x); + + } + + //! Compute criticality measure for the current sub-problem + field_type criticality2() const { + BoxConstraint<field_type,1> oneConstraint; oneConstraint[0] = {-1, 1}; // the criticality measure is given by computing the projected gradient path @@ -196,7 +218,6 @@ protected: gradient[j][k] = -correction_[j][k]; } } - return gradient.two_norm(); }