ellipsoidgridfunction.hh

// -*- 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_ELLIPSOID_GRIDFUNCTION_HH
#define DUNE_CURVED_SURFACE_GRID_ELLIPSOID_GRIDFUNCTION_HH

#include <cmath>

#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/common/math.hh>
#include "analyticgridfunction.hh"

namespace Dune
{
  // Ellipsoid functor
  template< class T >
  class EllipsoidProjection
  {
    using Domain = FieldVector<T,3>;
    using Jacobian = FieldMatrix<T,3,3>;

    T a_;
    T b_;
    T c_;

  public:
    //! Constructor of ellipsoid by major axes
    EllipsoidProjection (T a, T b, T c)
      : a_(a)
      , b_(b)
      , c_(c)
    {}

    //! project the coordinate to the ellipsoid
    // NOTE: This is not a closest-point projection, but a spherical-coordinate projection
    Domain operator() (const Domain& X) const
    {
      using std::sqrt;
      T x = X[0], y = X[1], z = X[2];
      T nrm = sqrt(x*x + y*y + z*z);
      return { a_*x/nrm , b_*y/nrm , c_*z/nrm };
    }

    //! derivative of the projection
    friend auto derivative (const EllipsoidProjection& ellipsoid)
    {
      return [a=ellipsoid.a_,b=ellipsoid.b_,c=ellipsoid.c_](Domain const& X)
      {
        using std::sqrt;
        T x = X[0], y = X[1], z = X[2];
        T x2 = x*x, y2 = y*y, z2 = z*z;

        T nrm2 = x2 + y2 + z2;
        T nrm3 = sqrt(nrm2)*nrm2;
        return Jacobian{
          { a*(y2 + z2)/nrm3 ,  -b*y*x/nrm3 ,       -c*z*x/nrm3 },
          {-a*x*y/nrm3 ,         b*(x2 + z2)/nrm3 , -c*z*y/nrm3 },
          {-a*x*z/nrm3 ,        -b*y*z/nrm3 ,        c*(x2 + y2)/nrm3 }
        };
      };
    }

    //! Normal vector
    Domain normal (const Domain& X) const
    {
      using std::sqrt;
      T x = X[0], y = X[1], z = X[2];
      T a2 = a_*a_, b2 = b_*b_, c2 = c_*c_;

      auto div = sqrt(b2*b2*c2*c2*x*x + a2*a2*c2*c2*y*y + a2*a2*b2*b2*z*z);
      return { b2*c2*x/div , a2*c2*y/div , a2*b2*z/div };
    }

    //! Mean curvature
    T mean_curvature (const Domain& X) const
    {
      using std::sqrt; using std::abs;
      T x = X[0], y = X[1], z = X[2];
      T x2 = x*x, y2 = y*y, z2 = z*z;
      T a2 = a_*a_, b2 = b_*b_, c2 = c_*c_;

      auto div = 2*a2*b2*c2 * power(sqrt(x2/(a2*a2) + y2/(b2*b2) + z2/(c2*c2)), 3);
      return abs(x2 + y2 + z2 - a2 - b2 - c2)/div;
    }

    //! Gaussian curvature
    T gauss_curvature (const Domain& X) const
    {
      T x = X[0], y = X[1], z = X[2];
      T x2 = x*x, y2 = y*y, z2 = z*z;
      T a2 = a_*a_, b2 = b_*b_, c2 = c_*c_;

      auto div = a2*b2*c2 * power(x2/(a2*a2) + y2/(b2*b2) + z2/(c2*c2), 2);
      return T(1)/div;
    }

  private:
    std::array<T,2> angles (Domain x) const
    {
      using std::acos; using std::atan2;
      x /= x.two_norm();

      return { atan2(x[1], x[0]), acos(x[2]) };
    }
  };

  //! construct a grid function representing a sphere parametrization
  template< class Grid, class T >
  auto ellipsoidGridFunction (T a, T b, T c)
  {
    return analyticGridFunction<Grid>(EllipsoidProjection<T>{a,b,c});
  }

} // end namespace Dune

#endif // DUNE_CURVED_SURFACE_GRID_ELLIPSOID_GRIDFUNCTION_HH