static const Float R = 1;
static const Float r = 0.6;
Float phi (const point& x) {
 return sqr(sqrt(sqr(x[0])+sqr(x[1]))-sqr(R)) + sqr(x[2])-sqr(r);
}
void get_torus_coordinates (const point& x,
                            Float& rho, Float& theta, Float& phi) {
 static const Float pi = acos(Float(-1));
 rho = sqrt(sqr(x[2]) + sqr(sqrt(sqr(x[0]) + sqr(x[1])) - sqr(R)));
 phi = atan2(x[1], x[0]);
 theta = atan2(x[2], sqrt(sqr(x[0]) + sqr(x[1])) - R);
}
struct u_exact : field_functor<u_exact,Float> {
 Float operator() (const point& x) const {
  Float rho, theta, phi;
  get_torus_coordinates (x, rho, theta, phi);
  return sin(3*phi)*cos(3*theta+phi);
 }
 u_exact (size_t d=3) {}
};
struct f : field_functor<f,Float> {
 Float operator() (const point& x) const {
  Float rho, theta, phi;
  get_torus_coordinates (x, rho, theta, phi);
  Float fx = (9*sin(3*phi)*cos(3*theta+phi))/sqr(r)
    - (-10*sin(3*phi)*cos(3*theta+phi) - 6*cos(3*phi)*sin(3*theta+phi))
      /sqr(R + r*cos(theta))
    - (3*sin(theta)*sin(3*phi)*sin(3*theta+phi))
      /(r*(R + r*cos(theta)));
  return fx;
 }
 f (size_t d=3) {}
};
