// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_CURVED_SURFACE_GRID_TORUS_GRIDFUNCTION_HH
#define DUNE_CURVED_SURFACE_GRID_TORUS_GRIDFUNCTION_HH

#include <type_traits>

#include <dune/common/math.hh>
#include <dune/curvedsurfacegrid/gridfunctions/analyticgridfunction.hh>
#include <dune/functions/common/defaultderivativetraits.hh>

namespace Dune
{
  // torus functor
  template< class T >
  struct TorusProjection
  {
    T R_;
    T r_;

    TorusProjection (T R, T r)
      : R_(R)
      , r_(r)
    {}

    // closest point projection
    FieldVector<T,3> operator()(FieldVector<T,3> x) const
    {
      using std::abs;

      FieldVector<T,3> x2{x[0], x[1], T(0)};
      T norm2 = x2.two_norm();
      assert(norm2 > std::numeric_limits<T>::epsilon());

      FieldVector<T,3> c = x2 * (R_ / norm2);
      x2 = x - c;

      norm2 = x2.two_norm();
      assert(norm2 > std::numeric_limits<T>::epsilon());
      x = c + x2 * (r_ / norm2);

      assert(abs(phi(x)) < 10*std::numeric_limits<T>::epsilon());
      return x;
    }

    T phi(FieldVector<T,3> const& x) const
    {
      using std::sqrt;
      T phi0 = sqrt(x[0]*x[0] + x[1]*x[1]) - R_;
      return phi0*phi0 + x[2]*x[2] - r_*r_;
    }

    FieldVector<T,3> normal (FieldVector<T,3> X) const
    {
      using std::sqrt;
      X = (*this)(X);

      T x = X[0], y = X[1], z = X[2];
      // T x2y2_1_2 = sqrt(x*x + y*y);
      // assert(x2y2_1_2 > std::numeric_limits<T>::epsilon());

      // return { x - 2*x/x2y2_1_2, y - 2*y/x2y2_1_2, z };

      T x2 = x*x, y2 = y*y, z2 = z*z;
      auto factor1 = x2 + y2 + z2 - 5;
      auto factor2 = x2 + y2 + z2 + 3;
      auto factor1_2 = factor1*factor1;
      auto factor2_2 = factor2*factor2;
      auto div = sqrt(x2*factor1_2 + y2*factor1_2 + z2*factor2_2);
      return {
        x*factor1/div,
        y*factor1/div,
        z*factor2/div
      };

    }

    T mean_curvature (FieldVector<T,3> X) const
    {
      using std::sqrt;
      X = (*this)(X);

      T x = X[0], y = X[1], z = X[2];
      // T x2y2_1_2 = sqrt(x*x + y*y);

      // return -(3 - 2/x2y2_1_2)/2;

      T x2 = x*x, y2 = y*y, z2 = z*z;
      T x3 = x*x2,y3 = y*y2,z3 = z*z2;
      auto factor1 = x2 + y2 + z2 - 5;
      auto factor2 = x2 + y2 + z2 + 3;
      auto factor1_2 = factor1*factor1;
      auto factor2_2 = factor2*factor2;
      auto div = sqrt(x2*factor1_2 + y2*factor1_2 + z2*factor2_2);
      auto div2= power(div, 3);
      return -(2*x2/div + x*factor1*(-2*x3*factor1 - 2*x*y2*factor1 - 2*x*z2*factor2 - x*factor1_2)/div2 + 2*y2/div + y*factor1*(-2*x2*y*factor1 - 2*y3*factor1 - 2*y*z2*factor2 - y*factor1_2)/div2 + 2*z2/div + z*factor2*(-2*x2*z*factor1 - 2*y2*z*factor1 - 2*z3*factor2 - z*factor2_2)/div2 + 2*factor1/div + factor2/div) ;
    }

    //! surface area of the torus = 4*pi^2*r1*r2
    T area () const
    {
      return 4*M_PI*M_PI*R_*r_;
    }
  };

  //! construct a grid function representing a torus parametrization
  template< class Grid, class T >
  auto torusGridFunction (T R, T r)
  {
    return analyticGridFunction<Grid>(TorusProjection<T>{R,r});
  }

} // end namespace Dune

#endif // DUNE_CURVED_SURFACE_GRID_TORUS_GRIDFUNCTION_HH