ellipsoidgridfunction.hh 3.11 KB
Newer Older
1
2
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
3
4
#ifndef DUNE_CURVED_SURFACE_GRID_ELLIPSOID_GRIDFUNCTION_HH
#define DUNE_CURVED_SURFACE_GRID_ELLIPSOID_GRIDFUNCTION_HH
5

6
#include <cmath>
7

8
9
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
Praetorius, Simon's avatar
Praetorius, Simon committed
10
11
12
#include <dune/common/math.hh>
#include "analyticgridfunction.hh"

13
14
15
16
17
18
namespace Dune
{
  // Ellipsoid functor
  template< class T >
  class EllipsoidProjection
  {
19
20
21
    using Domain = FieldVector<T,3>;
    using Jacobian = FieldMatrix<T,3,3>;

22
23
24
25
26
27
28
29
30
31
32
33
34
    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
35
    // NOTE: This is not a closest-point projection, but a spherical-coordinate projection
36
37
    Domain operator() (const Domain& X) const
    {
38
39
40
41
      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 };
42
43
44
    }

    //! derivative of the projection
Praetorius, Simon's avatar
Praetorius, Simon committed
45
    friend auto derivative (const EllipsoidProjection& ellipsoid)
46
    {
47
      return [a=ellipsoid.a_,b=ellipsoid.b_,c=ellipsoid.c_](Domain const& X)
48
      {
Praetorius, Simon's avatar
Praetorius, Simon committed
49
        using std::sqrt;
50
51
        T x = X[0], y = X[1], z = X[2];
        T x2 = x*x, y2 = y*y, z2 = z*z;
52
53
54
55
56
57
58

        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 }
59
60
61
62
63
64
65
66
67
68
69
70
        };
      };
    }

    //! 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);
71
      return { b2*c2*x/div , a2*c2*y/div , a2*b2*z/div };
72
    }
73

74
75
76
77
78
    //! 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];
79
      T x2 = x*x, y2 = y*y, z2 = z*z;
80
81
      T a2 = a_*a_, b2 = b_*b_, c2 = c_*c_;

82
83
      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;
84
85
86
    }

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

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

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

103
      return { atan2(x[1], x[0]), acos(x[2]) };
104
105
106
107
108
109
110
111
112
113
114
115
    }
  };

  //! 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

116
#endif // DUNE_CURVED_SURFACE_GRID_ELLIPSOID_GRIDFUNCTION_HH