void streamf_cavity_solver (const space& Xh, solver& sa, form& a, form& b) {
  string valued =  (Xh.size() == 3) ? "vector" : "scalar";
  space Ph (Xh.get_geo(), "P2", valued);
  Ph.block("top");  Ph.block("bottom");
  size_t d = Xh.get_geo().dimension();
  if (d == 3) {
    Ph.block("back"); Ph.block("front");
  } else {
    Ph.block("left"); Ph.block("right");
  }
  trial u (Xh), psi (Ph); test phi (Ph);
  a = (d == 3) ? integrate (ddot(grad(psi),grad(phi)))
               : integrate ( dot(grad(psi),grad(phi)));
  b = (d == 3) ? integrate (dot(curl(u),phi))
               : integrate (curl(u)*phi);
  sa = solver (a.uu());
}
field streamf (const field& uh, const solver& sa, const form& a, const form& b) {
  field lh = b*uh;
  field psi_h (b.get_second_space(), 0);
  psi_h.set_u() = sa.solve (lh.u() - a.ub()*psi_h.b());
  return psi_h;
}
