torusgridfunction.hh 2.26 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_TORUS_GRIDFUNCTION_HH
#define DUNE_CURVED_SURFACE_GRID_TORUS_GRIDFUNCTION_HH
5

6
7
8
#include <cassert>
#include <cmath>
#include <limits>
9
10
11
#include <type_traits>

#include <dune/common/math.hh>
12

13
#include "analyticgridfunction.hh"
14
15
16
17
18

namespace Dune
{
  // torus functor
  template< class T >
19
  class TorusProjection
20
  {
21
22
23
    using Domain = FieldVector<T,3>;
    using Jacobian = FieldMatrix<T,3,3>;

24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
    // Implicit function representation of the torus surface
    struct Phi
    {
      T R_, r_;

      // phi(x,y,z) = (sqrt(x^2 + y^2) - R)+2 + z^2 = r^2
      T operator() (const Domain& x) const
      {
        T phi0 = sqrt(x[0]*x[0] + x[1]*x[1]);
        return (phi0 - R_)*(phi0 - R_) + x[2]*x[2] - r_*r_;
      }

      // grad(phi)
      friend auto derivative (Phi phi)
      {
        return [R=phi.R_,r=phi.r_](const Domain& x) -> Domain
        {
          T phi0 = sqrt(x[0]*x[0] + x[1]*x[1]);
          return {
            -2*R*x[0]/phi0 + 2*x[0] ,
            -2*R*x[1]/phi0 + 2*x[1] ,
             2*x[2]
          };
        };
      }
    };
50

51
  public:
52
    //! Construction of the torus function by major and minor radius
53
54
55
    TorusProjection (T R, T r)
      : R_(R)
      , r_(r)
56
      , implicit_(Phi{R,r},100)
57
58
    {}

59
60
    //! Closest point projection
    Domain operator() (const Domain& x) const
61
    {
62
      return implicit_(x);
63
    }
64

65
66
    //! Normal vector
    Domain normal (const Domain& x) const
67
    {
68
      return implicit_.normal(x);
69
70
    }

71
    //! Mean curvature
Praetorius, Simon's avatar
Praetorius, Simon committed
72
    T mean_curvature (FieldVector<T,3> X) const
73
    {
74
75
76
77
      using std::sqrt;
      X = (*this)(X);

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

      return -(3 - 2/x2y2_1_2)/2;
81
82
    }

83
    //! surface area of the torus = 4*pi^2*R*r
84
85
    T area () const
    {
86
      return 4*M_PI*M_PI*R_*r_;
87
    }
88
89
90
91

  private:
    T R_, r_;
    ImplicitSurfaceProjection<Phi> implicit_;
92
93
94
95
  };

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

} // end namespace Dune

103
#endif // DUNE_CURVED_SURFACE_GRID_TORUS_GRIDFUNCTION_HH