field elasticity_criterion (Float lambda, const field& uh) {
  string grad_approx = (uh.get_approx() == "P2") ? "P1d" : "P0";
  space Xh  (uh.get_geo(), grad_approx);
  if (grad_approx == "P0") return interpolate (Xh, norm(uh));
  space T0h  (uh.get_geo(), grad_approx);
  size_t d = uh.get_geo().dimension();
  tensor I = tensor::eye (d);
  return interpolate (T0h, sqrt(2*norm2(D(uh)) + lambda*sqr(div(uh))));
}
