-
Notifications
You must be signed in to change notification settings - Fork 79
/
trajecs.jl
67 lines (52 loc) · 1.92 KB
/
trajecs.jl
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
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------
"""
CylindricalTrajectory(centroids, radius)
Trajectory of cylinders of given `radius` positioned at the `centroids`.
"""
struct CylindricalTrajectory{T} <: Domain{3,T}
centroids::Vector{Point{3,T}}
radius::T
end
CylindricalTrajectory(centroids::AbstractVector{Point{3,T}}, radius) where {T} =
CylindricalTrajectory(centroids, T(radius))
CylindricalTrajectory(centroids) = CylindricalTrajectory(centroids, 1)
topology(t::CylindricalTrajectory) = GridTopology(length(t.centroids))
function element(t::CylindricalTrajectory{T}, ind::Int) where {T}
c = t.centroids
r = t.radius
n = length(c)
if n == 1 # single vertical cylinder
p₁ = c[1] - Vec{3,T}(0, 0, 0.5)
p₂ = c[1] + Vec{3,T}(0, 0, 0.5)
return Cylinder(p₁, p₂, r)
end
if ind == 1 # head of trajectory
# points at cylinder planes
p₂ = center(Segment(c[ind], c[ind + 1]))
p₁ = p₂ - 2 * (p₂ - c[ind])
# normals to cylinder planes
n₂ = c[ind + 1] - c[ind]
n₁ = n₂
elseif ind == n # tail of trajectory
# points at cylinder planes
p₁ = center(Segment(c[ind - 1], c[ind]))
p₂ = p₁ + 2 * (c[ind] - p₁)
# normals to cylinder planes
n₁ = c[ind] - c[ind - 1]
n₂ = n₁
else # middle of trajectory
# points at cylinder planes
p₁ = center(Segment(c[ind - 1], c[ind]))
p₂ = center(Segment(c[ind], c[ind + 1]))
# normals to cylinder planes
n₁ = c[ind] - c[ind - 1]
n₂ = c[ind + 1] - c[ind]
end
# cylinder with given radius and planes
Cylinder(Plane(p₁, n₁), Plane(p₂, n₂), r)
end
nelements(t::CylindricalTrajectory) = length(t.centroids)
Base.eltype(t::CylindricalTrajectory) = typeof(first(t))
radius(t::CylindricalTrajectory) = t.radius