///
/// This file is part of Rheolef.
///
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
///
/// Rheolef is free software; you can redistribute it and/or modify
/// it under the terms of the GNU General Public License as published by
/// the Free Software Foundation; either version 2 of the License, or
/// (at your option) any later version.
///
/// Rheolef is distributed in the hope that it will be useful,
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
/// GNU General Public License for more details.
///
/// You should have received a copy of the GNU General Public License
/// along with Rheolef; if not, write to the Free Software
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
/// 
/// =========================================================================
//! @examplefile torus.icc The torus benchmark -- level set, right-hand-side and exact solution
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 {
 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 {
 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) {}
};
