Mercurial > repos > public > sbplib_julia

comparison src/Grids/geometry.jl @ 2003:524a52f190d7 feature/sbp_operators/laplace_curvilinear

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge feature/grids/geometry_functions
author Jonatan Werpers <jonatan@werpers.com>
date Tue, 29 Apr 2025 09:03:05 +0200
parents 730c9848ad0b
children d0b6c63c506e
comparison
equal deleted inserted replaced
1962:496b8d3903c6 2003:524a52f190d7
1 struct Line{PT} 1 struct Line{PT}
2 p::PT 2 p::PT
3 tangent::PT 3 tangent::PT
4
5 Line{PT}(p::PT, tangent::PT) where PT = new{PT}(p,tangent)
6 end
7
8 """
9 Line(p,t)
10
11 A line, as a callable object, starting at `p` with tangent `t`.
12
13 # Example
14 ```julia-repl
15 julia> l = Grids.Line([1,1],[2,1])
16 Diffinitive.Grids.Line{StaticArraysCore.SVector{2, Int64}}([1, 1], [2, 1])
17
18 julia> l(0)
19 2-element StaticArraysCore.SVector{2, Int64} with indices SOneTo(2):
20 1
21 1
22
23 julia> l(1)
24 2-element StaticArraysCore.SVector{2, Int64} with indices SOneTo(2):
25 3
26 2
27 ```
28
29 See also: [`LineSegment`](@ref).
30 """
31 function Line(p, t)
32 p = SVector{length(p)}(p)
33 t = SVector{length(t)}(t)
34 p, t = promote(p, t)
35
36 return Line{typeof(p)}(p,t)
4 end 37 end
5 38
6 (c::Line)(s) = c.p + s*c.tangent 39 (c::Line)(s) = c.p + s*c.tangent
7 40
41 Grids.jacobian(l::Line, t) = l.tangent
8 42
9 struct LineSegment{PT} 43 struct LineSegment{PT}
10 a::PT 44 a::PT
11 b::PT 45 b::PT
46
47 LineSegment{PT}(p::PT, tangent::PT) where PT = new{PT}(p,tangent)
48 end
49
50 """
51 LineSegment(a,b)
52
53 A line segment, as a callable object, from `a` to `b`.
54
55 # Example
56 ```julia-repl
57 julia> l = Grids.LineSegment([1,1],[2,1])
58 Diffinitive.Grids.LineSegment{StaticArraysCore.SVector{2, Int64}}([1, 1], [2, 1])
59
60 julia> l(0)
61 2-element StaticArraysCore.SVector{2, Int64} with indices SOneTo(2):
62 1
63 1
64
65 julia> l(0.5)
66 2-element StaticArraysCore.SVector{2, Float64} with indices SOneTo(2):
67 1.5
68 1.0
69
70 julia> l(1)
71 2-element StaticArraysCore.SVector{2, Int64} with indices SOneTo(2):
72 2
73 1
74 ```
75
76 See also: [`Line`](@ref).
77 """
78 function LineSegment(a, b)
79 a = SVector{length(a)}(a)
80 b = SVector{length(b)}(b)
81 a, b = promote(a, b)
82
83 return LineSegment{typeof(a)}(a,b)
12 end 84 end
13 85
14 (c::LineSegment)(s) = (1-s)*c.a + s*c.b 86 (c::LineSegment)(s) = (1-s)*c.a + s*c.b
15 87
16 88 Grids.jacobian(c::LineSegment, s) = c.b - c.a
89
90 """
91 linesegments(ps...)
92
93 An array of line segments between the points `ps[1]`, `ps[2]`, and so on.
94
95 See also: [`polygon_edges`](@ref).
96 """
17 function linesegments(ps...) 97 function linesegments(ps...)
18 return [LineSegment(ps[i], ps[i+1]) for i ∈ 1:length(ps)-1] 98 return [LineSegment(ps[i], ps[i+1]) for i ∈ 1:length(ps)-1]
19 end 99 end
20 100
21 101
102 """
103 polygon_edges(ps...)
104
105 An array of line segments between the points `ps[1]`, `ps[2]`, and so on
106 including the segment between `ps[end]` and `ps[1]`.
107
108 See also: [`linesegments`](@ref).
109 """
22 function polygon_edges(ps...) 110 function polygon_edges(ps...)
23 n = length(ps) 111 n = length(ps)
24 return [LineSegment(ps[i], ps[mod1(i+1,n)]) for i ∈ eachindex(ps)] 112 return [LineSegment(ps[i], ps[mod1(i+1,n)]) for i ∈ eachindex(ps)]
25 end 113 end
26 114
27 struct Circle{T,PT} 115 struct Circle{PT,T}
28 c::PT 116 c::PT
29 r::T 117 r::T
118
119 Circle{PT,T}(c,r) where {PT,T} = new{PT,T}(c,r)
120 end
121
122 """
123 Circle(c,r)
124
125 A circle with center `c` and radius `r` paramatrized with the angle to the x-axis.
126
127 # Example
128 ```julia-repl
129 julia> c = Grids.Circle([1,1], 2)
130 Diffinitive.Grids.Circle{StaticArraysCore.SVector{2, Int64}, Int64}([1, 1], 2)
131
132 julia> c(0)
133 2-element StaticArraysCore.SVector{2, Float64} with indices SOneTo(2):
134 3.0
135 1.0
136
137 julia> c(π/2)
138 2-element StaticArraysCore.SVector{2, Float64} with indices SOneTo(2):
139 1.0000000000000002
140 3.0
141 ```
142 """
143 function Circle(c,r)
144 c = SVector{2}(c)
145 return Circle{typeof(c), typeof(r)}(c,r)
30 end 146 end
31 147
32 function (C::Circle)(θ) 148 function (C::Circle)(θ)
33 (;c, r) = C 149 (;c, r) = C
34 c + r*@SVector[cos(θ), sin(θ)] 150 c + r*@SVector[cos(θ), sin(θ)]
35 end 151 end
36 152
153 function Grids.jacobian(C::Circle, θ)
154 (;r) = C
155 r*@SVector[-sin(θ), cos(θ)]
156 end
157
158 """
159 TransfiniteInterpolationSurface(c₁, c₂, c₃, c₄)
160
161 A surface defined by the transfinite interpolation of the curves `c₁`, `c₂`, `c₃`, and `c₄`.
162 """
37 struct TransfiniteInterpolationSurface{T1,T2,T3,T4} 163 struct TransfiniteInterpolationSurface{T1,T2,T3,T4}
38 c₁::T1 164 c₁::T1
39 c₂::T2 165 c₂::T2
40 c₃::T3 166 c₃::T3
41 c₄::T4 167 c₄::T4
54 180
55 function (s::TransfiniteInterpolationSurface)(ξ̄::AbstractArray) 181 function (s::TransfiniteInterpolationSurface)(ξ̄::AbstractArray)
56 s(ξ̄...) 182 s(ξ̄...)
57 end 183 end
58 184
59 # TODO: Implement jacobian() for the different mapping helpers 185 """
186 check_transfiniteinterpolation(s::TransfiniteInterpolationSurface)
187
188 Throw an error if the ends of the curves in the transfinite interpolation do not match.
189 """
190 function check_transfiniteinterpolation(s::TransfiniteInterpolationSurface)
191 if check_transfiniteinterpolation(Bool, s)
192 return nothing
193 else
194 throw(ArgumentError("The end of each curve in the transfinite interpolation should be the same as the beginning of the next curve."))
195 end
196 end
197
198 """
199 check_transfiniteinterpolation(Bool, s::TransfiniteInterpolationSurface)
200
201 Return true if the ends of the curves in the transfinite interpolation match.
202 """
203 function check_transfiniteinterpolation(::Type{Bool}, s::TransfiniteInterpolationSurface)
204 if !isapprox(s.c₁(1), s.c₂(0))
205 return false
206 end
207
208 if !isapprox(s.c₂(1), s.c₃(0))
209 return false
210 end
211
212 if !isapprox(s.c₃(1), s.c₄(0))
213 return false
214 end
215
216 if !isapprox(s.c₄(1), s.c₁(0))
217 return false
218 end
219
220 return true
221 end
222
223 function Grids.jacobian(s::TransfiniteInterpolationSurface, ξ̄)
224 u, v = ξ̄
225
226 c₁, c₂, c₃, c₄ = s.c₁, s.c₂, s.c₃, s.c₄
227 P₀₀ = c₁(0)
228 P₁₀ = c₂(0)
229 P₁₁ = c₃(0)
230 P₀₁ = c₄(0)
231
232 ∂x̄∂ξ₁ = (1-v)*jacobian(c₁,u) + c₂(v) - v*jacobian(c₃,1-u) -c₄(1-v) - (
233 -(1-v)*P₀₀ + (1-v)*P₁₀ + v*P₁₁ - v*P₀₁
234 )
235
236 ∂x̄∂ξ₂ = -c₁(u) + u*jacobian(c₂,v) + c₃(1-u) - (1-u)*jacobian(c₄,1-v) - (
237 -(1-u)*P₀₀ - u*P₁₀ + u*P₁₁ + (1-u)*P₀₁
238 )
239
240 return [∂x̄∂ξ₁ ∂x̄∂ξ₂]
241 end