auto && tau_star{.5*Hooke::evaluate_stress(this->lambda, this->mu, ln_be_star)};
// deviatoric part of Kirchhoff stress
// TODO: There seems to be disagreement regarding whether the deviator in 2D also needs to have 3 times the trace subtracted or just 2 times. Have an eye on this!
auto && tau_d_star{tau_star - tau_star.trace()/3*tau_star.Identity()};
auto && tau_eq_star{std::sqrt(3*.5*(tau_d_star.array()*
tau_d_star.transpose().array()).sum())};
auto && N_star{3*.5*tau_d_star/tau_eq_star};
// this is eq (27), and the std::min enforces the Kuhn-Tucker relation (16)
auto && phi_star{std::min(tau_eq_star - this->tau_y0 - this->H * eps_p.old(), 0.)};
// return mapping
auto && Del_gamma{phi_star/(this->H + 3 * this->mu)};
auto && tau{tau_star - 2*Del_gamma*this->mu*N_star};