Mercurial > repos > public > sbplib_julia

comparison src/Grids/manifolds.jl @ 1954:b0915f43b122 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 Sat, 08 Feb 2025 09:38:58 +0100
parents dd77b45ee1ac
children 3fb5b03162ee
comparison
equal deleted inserted replaced
1953:835b1dcee38e 1954:b0915f43b122
1 """
2 ParameterSpace{D}
3
4 A space of parameters of dimension `D`. Used with `Chart` to indicate which
5 parameters are valid for that chart.
6
7 Common parameter spaces are created using the functions unit sized spaces
8 * `unitinterval`
9 * `unitrectangle`
10 * `unitbox`
11 * `unittriangle`
12 * `unittetrahedron`
13 * `unithyperbox`
14 * `unitsimplex`
15
16 See also: [`Interval`](@ref), [`Rectangle`](@ref), [`Box`](@ref),
17 [`Triangle`](@ref), [`Tetrahedron`](@ref), [`HyperBox`](@ref),
18 [`Simplex`](@ref),
19 """
20 abstract type ParameterSpace{D} end
21 Base.ndims(::ParameterSpace{D}) where D = D
22
23 struct Interval{T} <: ParameterSpace{1}
24 a::T
25 b::T
26
27 function Interval(a,b)
28 a, b = promote(a, b)
29 new{typeof(a)}(a,b)
30 end
31 end
32
33 limits(i::Interval) = (i.a, i.b)
34
35 unitinterval(T=Float64) = Interval(zero(T), one(T))
36
37
38 struct HyperBox{T,D} <: ParameterSpace{D}
39 a::SVector{D,T}
40 b::SVector{D,T}
41 end
42
43 function HyperBox(a,b)
44 ET = promote_type(eltype(a),eltype(b))
45 T = SVector{length(a),ET}
46 HyperBox(convert(T,a), convert(T,b))
47 end
48
49 Rectangle{T} = HyperBox{T,2}
50 Box{T} = HyperBox{T,3}
51
52 limits(box::HyperBox, d) = (box.a[d], box.b[d])
53 limits(box::HyperBox) = (box.a, box.b)
54
55 unitsquare(T=Float64) = unithyperbox(T,2)
56 unitcube(T=Float64) = unithyperbox(T,3)
57 unithyperbox(T, D) = HyperBox((@SVector zeros(T,D)), (@SVector ones(T,D)))
58 unithyperbox(D) = unithyperbox(Float64,D)
59
60
61 struct Simplex{T,D,NV} <: ParameterSpace{D}
62 verticies::NTuple{NV,SVector{D,T}}
63 end
64
65 function Simplex(verticies::Vararg{AbstractArray})
66 ET = mapreduce(eltype,promote_type,verticies)
67 T = SVector{length(verticies[1]),ET}
68
69 return Simplex(Tuple(convert(T,v) for v ∈ verticies))
70 end
71
72 verticies(s::Simplex) = s.verticies
73
74 Triangle{T} = Simplex{T,2}
75 Tetrahedron{T} = Simplex{T,3}
76
77 unittriangle(T=Float64) = unitsimplex(T,2)
78 unittetrahedron(T=Float64) = unitsimplex(T,3)
79 function unitsimplex(T,D)
80 z = @SVector zeros(T,D)
81 unitelement = one(eltype(z))
82 verticies = ntuple(i->setindex(z, unitelement, i), D)
83 return Simplex((z,verticies...))
84 end
85 unitsimplex(D) = unitsimplex(Float64, D)
86
87 """ 1 """
88 Chart{D} 2 Chart{D}
89 3
90 A parametrized description of a manifold or part of a manifold. 4 A parametrized description of a manifold or part of a manifold.
91 """ 5 """
117 which will both allow calling `jacobian(c,ξ)`. 31 which will both allow calling `jacobian(c,ξ)`.
118 """ 32 """
119 jacobian(c::Chart, ξ) = jacobian(c.mapping, ξ) 33 jacobian(c::Chart, ξ) = jacobian(c.mapping, ξ)
120 # TBD: Can we register a error hint for when jacobian is called with a function that doesn't have a registered jacobian? 34 # TBD: Can we register a error hint for when jacobian is called with a function that doesn't have a registered jacobian?
121 35
36 boundary_identifiers(c::Chart) = boundary_identifiers(parameterspace(c))
122 37
123 # TBD: Should Charts, parameterspaces have boundary names?
124 38
125 """ 39 """
126 Atlas 40 Atlas
127 41
128 A collection of charts and their connections. 42 A collection of charts and their connections.
129 Should implement methods for `charts` and 43 Should implement methods for `charts` and `connections`.
130 """ 44 """
131 abstract type Atlas end 45 abstract type Atlas end
132 46
133 """ 47 """
134 charts(::Atlas) 48 charts(::Atlas)
136 The colloction of charts in the atlas. 50 The colloction of charts in the atlas.
137 """ 51 """
138 function charts end 52 function charts end
139 53
140 """ 54 """
141 connections 55 connections(::Atlas)
142 56
143 TBD: What exactly should this return? 57 Collection of pairs of multiblock boundary identifiers.
58 """
59 function connections end
60
144 61
145 """ 62 """
63 CartesianAtlas{D,C<:Chart,AT<:AbstractArray{C,D}} <: Atlas
146 64
147 struct CartesianAtlas <: Atlas 65 An atlas where the charts are arranged and connected like an array.
148 charts::Matrix{Chart} 66 """
67 struct CartesianAtlas{D,C<:Chart,AT<:AbstractArray{C,D}} <: Atlas
68 charts::AT
149 end 69 end
150 70
151 charts(a::CartesianAtlas) = a.charts 71 charts(a::CartesianAtlas) = a.charts
152 connections(a::CartesianAtlas) = nothing
153 72
154 struct UnstructuredAtlas <: Atlas 73 function connections(a::CartesianAtlas)
155 charts::Vector{Chart} 74 c = Tuple{MultiBlockBoundary, MultiBlockBoundary}[]
156 connections 75
76 for d ∈ 1:ndims(charts(a))
77 Is = eachslice(CartesianIndices(charts(a)); dims=d)
78 for i ∈ 1:length(Is)-1 # For each interface between slices
79 for jk ∈ eachindex(Is[i]) # For each block in slice
80 Iᵢⱼₖ = Tuple(Is[i][jk])
81 Iᵢ₊₁ⱼₖ = Tuple(Is[i+1][jk])
82 push!(c,
83 (
84 MultiBlockBoundary{Iᵢⱼₖ, CartesianBoundary{d,UpperBoundary}}(),
85 MultiBlockBoundary{Iᵢ₊₁ⱼₖ, CartesianBoundary{d,LowerBoundary}}(),
86 )
87 )
88 end
89 end
90 end
91
92 return c
93 end
94
95 """
96 boundary_identifiers(a::CartesianAtlas)
97
98 All non-connected boundaries of the charts of `a`.
99 """
100 function boundary_identifiers(a::CartesianAtlas)
101 bs = MultiBlockBoundary[]
102
103 for d ∈ 1:ndims(charts(a))
104 Is = eachslice(CartesianIndices(charts(a)); dims=d)
105
106 for (i,b) ∈ ((1,LowerBoundary),(length(Is),UpperBoundary)) # For first and last slice
107 for jk ∈ eachindex(Is[i]) # For each block in slice
108 Iᵢⱼₖ = Tuple(Is[i][jk])
109 push!(bs,
110 MultiBlockBoundary{Iᵢⱼₖ, CartesianBoundary{d,b}}(),
111 )
112 end
113 end
114 end
115
116 return bs
117 end
118
119
120 """
121 UnstructuredAtlas{C<:Chart, CN<:Tuple{MultiBlockBoundary,MultiBlockBoundary}, ...} <: Atlas
122
123 An atlas with connections determined by a vector `MultiBlockBoundary` pairs.
124 """
125 struct UnstructuredAtlas{C<:Chart, CN<:Tuple{MultiBlockBoundary,MultiBlockBoundary}, CV<:AbstractVector{C}, CNV<:AbstractVector{CN}} <: Atlas
126 charts::CV
127 connections::CNV
157 end 128 end
158 129
159 charts(a::UnstructuredAtlas) = a.charts 130 charts(a::UnstructuredAtlas) = a.charts
160 connections(a::UnstructuredAtlas) = nothing 131 connections(a::UnstructuredAtlas) = a.connections
161 132
133 """
134 boundary_identifiers(a::UnstructuredAtlas)
162 135
163 ### 136 All non-connected boundaries of the charts of `a`.
164 # Geometry 137 """
165 ### 138 function boundary_identifiers(a::UnstructuredAtlas)
139 bs = MultiBlockBoundary[]
166 140
167 abstract type Curve end 141 for (i,c) ∈ enumerate(charts(a))
168 abstract type Surface end 142 for b ∈ boundary_identifiers(c)
143 mbb = MultiBlockBoundary{i,typeof(b)}()
169 144
145 if !any(cn->mbb∈cn, connections(a))
146 push!(bs, mbb)
147 end
148 end
149 end
170 150
171 struct Line{PT} <: Curve 151 return bs
172 p::PT
173 tangent::PT
174 end 152 end
175
176 (c::Line)(s) = c.p + s*c.tangent
177
178
179 struct LineSegment{PT} <: Curve
180 a::PT
181 b::PT
182 end
183
184 (c::LineSegment)(s) = (1-s)*c.a + s*c.b
185
186
187 function linesegments(ps...)
188 return [LineSegment(ps[i], ps[i+1]) for i ∈ 1:length(ps)-1]
189 end
190
191
192 function polygon_edges(ps...)
193 n = length(ps)
194 return [LineSegment(ps[i], ps[mod1(i+1,n)]) for i ∈ eachindex(Ps)]
195 end
196
197 struct Circle{T,PT} <: Curve
198 c::PT
199 r::T
200 end
201
202 (c::Circle)(θ) = c.c + r*@SVector[cos(Θ), sin(Θ)]
203
204 struct TransfiniteInterpolationSurface{T1,T2,T3,T4} <: Surface
205 c₁::T1
206 c₂::T2
207 c₃::T3
208 c₄::T4
209 end
210
211 function (s::TransfiniteInterpolationSurface)(u,v)
212 c₁, c₂, c₃, c₄ = s.c₁, s.c₂, s.c₃, s.c₄
213 P₀₀ = c₁(0)
214 P₁₀ = c₂(0)
215 P₁₁ = c₃(0)
216 P₀₁ = c₄(0)
217 return (1-v)*c₁(u) + u*c₂(v) + v*c₃(1-u) + (1-u)*c₄(1-v) - (
218 (1-u)*(1-v)*P₀₀ + u*(1-v)*P₁₀ + u*v*P₁₁ + (1-u)*v*P₀₁
219 )
220 end
221
222 function (s::TransfiniteInterpolationSurface)(ξ̄::AbstractArray)
223 s(ξ̄...)
224 end
225
226 # TODO: Implement jacobian() for the different mapping helpers
227