field elasticity_solve (const space& Xh, Float lambda) {
  size_t d = Xh.get_geo().dimension();
  point f (0,0,0);
  f[d-1] = -1;
  trial u (Xh); test v (Xh);
  field lh = integrate (dot(f,v));
  form  a  = integrate (lambda*div(u)*div(v) + 2*ddot(D(u),D(v)));
  solver sa (a.uu());
  field uh (Xh, 0);
  uh.set_u() = sa.solve (lh.u() - a.ub()*uh.b());
  return uh;
}
