Mercurial > repos > public > sbplib_julia

comparison src/Grids/manifolds.jl @ 1558:81e97d3bec8c feature/grids/manifolds

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Start adding manifolds
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 24 Apr 2024 13:26:30 +0200
parents
children 35fe4375b35f
comparison
equal deleted inserted replaced
1556:ec5e7926c37b 1558:81e97d3bec8c
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
22 struct HyperBox{T,D} <: ParameterSpace{D}
23 a::SVector{D,T}
24 b::SVector{D,T}
25 end
26
27 function HyperBox(a,b)
28 T = SVector{length(a)}
29 HyperBox(convert(T,a), convert(T,b))
30 end
31
32 Interval{T} = HyperBox{T,1}
33 Rectangle{T} = HyperBox{T,2}
34 Box{T} = HyperBox{T,3}
35
36 limits(box::HyperBox, d) = (box.a[d], box.b[d])
37 limits(box::HyperBox) = (box.a, box.b)
38
39 unitinterval(T=Float64) = unithyperbox(T,1)
40 unitsquare(T=Float64) = unithyperbox(T,2)
41 unitcube(T=Float64) = unithyperbox(T,3)
42 unithyperbox(T, D) = HyperBox((@SVector zeros(T,D)), (@SVector ones(T,D)))
43 unithyperbox(D) = unithyperbox(Float64,D)
44
45
46 struct Simplex{T,D} <: ParameterSpace{D}
47 verticies::NTuple{D,SVector{D,T}}
48 end
49
50 Triangle{T} = Simplex{T,2}
51 Tetrahedron{T} = Simplex{T,3}
52
53 unittriangle(T) = unitsimplex(T,2)
54 unittetrahedron(T) = unitsimplex(T,3)
55 function unitsimplex(T,D)
56 z = @SVector zeros(T,D)
57 unitelement = one(eltype(z))
58 verticies = ntuple(i->setindex(z, unitelement, i), 4)
59 return Simplex(verticies)
60 end
61
62
63 """
64
65 A parametrized description of a manifold or part of a manifold.
66
67 Should implement a methods for
68 * `parameterspace`
69 * `(::Chart)(ξs...)`
70
71 There is a default implementation for `(::Chart{D})(::SVector{D})`
72 """
73 abstract type Chart{D} end
74 # abstract type Chart{D,R} end
75
76 domain_dim(::Chart{D}) where D = D
77 # range_dim(::Chart{D,R}) where {D,R} = R
78
79 """
80 The parameterspace of a chart
81 """
82 function parameterspace end
83
84 (c::Chart{D})(x̄::SVector{D}) where D = c(x̄...)
85
86
87 struct ConcereteChart{PST<:ParameterSpace, MT}
88 parameterspace::PST
89 mapping::MT
90 end
91
92 (c::Chart)(x̄) = c.mapping(x̄)
93
94
95 """
96 Atlas
97
98 A collection of charts and their connections.
99 Should implement methods for `charts` and
100 """
101 abstract type Atlas end
102
103 """
104 charts(::Atlas)
105
106 The colloction of charts in the atlas.
107 """
108 function charts end
109
110 """
111 connections
112
113 TBD: What exactly should this return?
114
115 """
116
117 struct CartesianAtlas <: Atlas
118 charts::Matrix{Chart}
119 end
120
121 charts(a::CartesianAtlas) = a.charts
122
123 struct UnstructuredAtlas <: Atlas
124 charts::Vector{Chart}
125 connections
126 end
127
128 charts(a::UnstructuredAtlas) = a.charts
129
130
131 ###
132 # Geometry
133 ###
134
135 abstract type Curve end
136 abstract type Surface end
137
138
139 struct Line{PT} <: Curve
140 p::PT
141 tangent::PT
142 end
143
144 (c::Line)(s) = c.p + s*c.tangent
145
146
147 struct LineSegment{PT} <: Curve
148 a::PT
149 b::PT
150 end
151
152 (c::LineSegment)(s) = (1-s)*c.a + s*c.b
153
154
155 struct Circle{T,PT} <: Curve
156 c::PT
157 r::T
158 end
159
160 (c::Circle)(θ) = c.c + r*@SVector[cos(Θ), sin(Θ)]
161
162 struct TransfiniteInterpolationSurface{T1,T2,T3,T4} <: Surface
163 c₁::T1
164 c₂::T2
165 c₃::T3
166 c₄::T4
167 end
168
169 function (s::TransfiniteInterpolationSurface)(u,v)
170 c₁, c₂, c₃, c₄ = s.c₁, s.c₂, s.c₃, s.c₄
171 P₀₀ = c₁(0)
172 P₁₀ = c₂(0)
173 P₁₁ = c₃(0)
174 P₀₁ = c₄(0)
175 return (1-v)*c₁(u) + u*c₂(v) + v*c₃(1-u) + (1-u)*c₄(1-v) - (
176 (1-u)*(1-v)*P₀₀ + u*(1-v)*P₁₀ + u*v*P₁₁ + (1-u)*v*P₀₁
177 )
178 end
179
180 function (s::TransfiniteInterpolationSurface)(ξ̄::AbstractArray)
181 s(ξ̄...)
182 end
183
184
185 function polygon_sides(Ps...)
186 n = length(Ps)
187 return [t->line(t,Ps[i],Ps[mod1(i+1,n)]) for i ∈ eachindex(Ps)]
188 end