torusgridfunction.hh 3.5 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
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
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
50
51

#include <type_traits>

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

namespace Dune
{
  // torus functor
  template< class T >
  struct TorusProjection
  {
    T radius1_;
    T radius2_;

    TorusProjection (T radius1, T radius2)
      : radius1_(radius1)
      , radius2_(radius2)
    {}

    //! project the coordinate to the sphere at origin with `radius`
    template< class Domain >
    Domain operator() (const Domain& x) const
    {
      using std::sqrt;
      auto scale1 = radius1_ / sqrt(x[0]*x[0] + x[1]*x[1]);
      Domain center{x[0] * scale1, x[1] * scale1, 0};
      Domain out{x[0] - center[0], x[1] - center[1], x[2]};
      out *= radius2_ / out.two_norm();

      return out + center;
    }

    //! derivative of the projection
    friend auto derivative (const TorusProjection& t)
    {
      return [r1=t.radius1_,r2=t.radius2_](auto const& x)
      {
        using std::sqrt;
        using Domain = std::decay_t<decltype(x)>;
        using DerivativeTraits = Functions::DefaultDerivativeTraits<Domain(Domain)>;
        using Matrix = typename DerivativeTraits::Range;

        T x0 = x[0]*x[0], x1 = x[1]*x[1], x2 = x0+x1, x3 = sqrt(x2), x4 = T(1)/x3, x5 = r1*x4;
        T x6 = sqrt(x2*x2*x2), x7 = r1/x6, x8 = x0*x7, x9 = T(1)-x5, x10 = x8+x9;
        T x11 = x[2]*x[2], x12 = x9*x9, x13 = x0*x12, x14 = x1*x12, x15 = x11+x13+x14;
52
        T x16 = r2/sqrt(x15), x17 = x1*x7, x18 = r2*(x10 + x17)/sqrt(power(x15,3));
53
        T x19 = sqrt(power(x2,5)), x20 = r1*x19, x21 = x2*x2, x22 = r1*x3;
54
        T x23 = x21-r1*x6, x24 = r1-x3, x25 = x24*x24, x26 = x0*x25+x1*x25+x11*x2;
55
56
57
58
59
60
61
62
63
64
65
66
67
        T x27 = x26/x2, x28 = sqrt(x27), x29 = r2*x25*x28, x30 = x2*x2*x2*x29, x31 = x26*x26;
        T x32 = sqrt(x27*x27*x27), x33 = -r1*r2*sqrt(power(x2,13))*x32 + r1*sqrt(power(x2,9))*x31;
        T x34 = T(1)/x31, x35 = x[0]*x[1]*x34/power(x2,6), x36 = r2*x[2], x37 = x[0]*x36;
        T x38 = (x21*(x0*x24 + x1*x24 - x3*(x11 + x25)) + x26*x6)/(x19*x26*x28), x39 = x36*x[1];
        T x40 = x24*x4/x32;

        return Matrix{
          {x10*x16 - x13*x18 + x5 - x8, -x35*(x30*(x1*x20 + x21*(x0*x22 + x23)) + x33), x37*x38},
          {-x35*(x30*(x0*x20 + x21*(x1*x22 + x23)) + x33), -x14*x18 + x16*(x17 + x9) - x17 + x5, x38*x39},
          {x37*x40, x39*x40, x21*x29*x34}
        };
      };
    }
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91

    // phi = (sqrt(x**2 + y**2) - r2)**2 + z**2 - r1**2
    // N = grad(phi) / |grad(phi)|
    template< class Domain >
    Domain normal (const Domain& x)
    {
      using std::sqrt;
      T x0 = x[0]*x[0], x1 = x[1]*x[1], x2 = x0+x1, x3 = sqrt(x2), x4 = radius2_-x3;
      T x5 = x4*x4, x6 = sqrt((x0*x5 + x1*x5 + x2*x[2]*x[2])/x2), x7 = x4/(x6*x3);

      return {-x[0]*x7, -x[1]*x7, x[2]/x6};
    }

    template< class Domain >
    auto mean_curvature (const Domain& x)
    {
      return 0.0;
    }

    //! surface area of the torus = 4*pi^2*r1*r2
    T area () const
    {
      return 4*M_PI*M_PI*radius1_*radius2_;
    }
92
93
94
95
96
97
98
99
100
101
102
  };

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

} // end namespace Dune

103
#endif // DUNE_CURVED_SURFACE_GRID_TORUS_GRIDFUNCTION_HH