torusgridfunction.hh 3.15 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 <type_traits>

#include <dune/common/math.hh>
9
#include <dune/curvedsurfacegrid/gridfunctions/analyticgridfunction.hh>
10
11
12
13
14
15
16
17
#include <dune/functions/common/defaultderivativetraits.hh>

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

21
22
23
    TorusProjection (T R, T r)
      : R_(R)
      , r_(r)
24
25
    {}

26
27
    // closest point projection
    FieldVector<T,3> operator()(FieldVector<T,3> x) const
28
    {
29
30
31
32
33
34
35
36
37
38
39
40
      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);
41

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

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

53
    FieldVector<T,3> normal (FieldVector<T,3> X) const
54
55
    {
      using std::sqrt;
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
      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
      };
75
76
77

    }

78
    T mean_curvature (FieldVector<T,3> X)
79
    {
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
      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) ;
97
98
99
100
101
    }

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

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

} // end namespace Dune

115
#endif // DUNE_CURVED_SURFACE_GRID_TORUS_GRIDFUNCTION_HH