Newer
Older
1
2
3
4
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
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
#ifndef DUNE_GFE_TWISTEDSTRIP_GRIDFUNCTION_HH
#define DUNE_GFE_TWISTEDSTRIP_GRIDFUNCTION_HH
#include <cmath>
#include <dune/common/fmatrix.hh>
#include <dune/common/fvector.hh>
#include <dune/curvedgrid/gridfunctions/analyticgridfunction.hh>
namespace Dune {
/// \brief Functor representing a twisted strip in 3D, with length length_ and nTwists twists
template <class T = double>
class TwistedStripProjection
{
static const int dim = 3;
using Domain = FieldVector<T,dim>;
using Jacobian = FieldMatrix<T,dim,dim>;
T length_;
double nTwists_;
public:
TwistedStripProjection (T length, double nTwists)
: length_(length),
nTwists_(nTwists)
{}
/// \brief Project the coordinate to the twisted strip using closest point projection
Domain operator() (const Domain& x) const
{
// Angle for the x - position of the point
double alpha = std::acos(-1)*nTwists_*x[0]/length_;
double cosalpha = std::cos(alpha);
double sinalpha = std::sin(alpha);
// Add the line from middle to the new point - the direction is (0,cosalpha,sinalpha)
// We need the distance, calculated by (y^2 + z^2)^0.5
double dist = std::sqrt(x[1]*x[1] + x[2]*x[2]);
double y = std::abs(cosalpha*dist);
if (x[1] < 0 and std::abs(x[1]) > std::abs(x[2])) { // only trust the sign of x[1] if it is larger than x[2]
y *= (-1);
} else if (std::abs(x[1]) <= std::abs(x[2])) {
if (cosalpha * sinalpha >= 0 and x[2] < 0) // if cosalpha*sinalpha > 0, x[1] and x[2] must have the same sign
y *= (-1);
if (cosalpha * sinalpha <= 0 and x[2] > 0) // if cosalpha*sinalpha < 0, x[1] and x[2] must have opposite signs
y *= (-1);
}
double z = std::abs(sinalpha*dist);
if (x[2] < 0 and std::abs(x[2]) > std::abs(x[1])) {
z *= (-1);
} else if (std::abs(x[2]) <= std::abs(x[1])) {
if (cosalpha * sinalpha >= 0 and x[1] < 0)
z *= (-1);
if (cosalpha * sinalpha <= 0 and x[1] > 0)
z *= (-1);
}
return {x[0], y , z};
}
/// \brief derivative of the projection to the twistedstrip
friend auto derivative (const TwistedStripProjection& twistedStrip)
{
DUNE_THROW(NotImplemented,"The derivative of the projection to the twisted strip is not implemented yet!");
return [length = twistedStrip.length_, nTwists = twistedStrip.nTwists_](const Domain& x)
{
Jacobian out;
return out;
};
}
/// \brief Normal Vector of the twisted strip
Domain normal (const Domain& x) const
{
// Angle for the x - position of the point
double alpha = std::acos(-1)*nTwists_*x[0]/length_;
std::cout << "normal at " << x[0] << "is " << -std::sin(alpha) << ", " << std::cos(alpha) << std::endl;
return {0,-std::sin(alpha), std::cos(alpha)};
}
/// \brief The mean curvature of the twisted strip
T mean_curvature (const Domain& /*x*/) const
{
DUNE_THROW(NotImplemented,"The mean curvature of the twisted strip is not implemented yet!");
return 1;
}
/// \brief The area of the twisted strip
T area (T width) const
{
return length_*width;
}
};
/// \brief construct a grid function representing the parametrization of a twisted strip
template <class Grid, class T>
auto twistedStripGridFunction (T length, double nTwists)
{
return analyticGridFunction<Grid>(TwistedStripProjection<T>{length, nTwists});
}
} // end namespace Dune
#endif // DUNE_GFE_TWISTEDSTRIP_GRIDFUNCTION_HH