Mercurial > repos > public > sbplib_julia

comparison src/Grids/mapped_grid.jl @ 1835:a6f28a8b8f3f refactor/lazy_tensors/elementwise_ops

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge default
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 09 Jan 2025 12:40:49 +0100
parents f21bfc5f21aa
children d91a9f47380f 85f8855473ab 1c58005429fd
comparison
equal deleted inserted replaced
1789:48eaa973159a 1835:a6f28a8b8f3f
1 """
2 MappedGrid{T,D} <: Grid{T,D}
3
4 A grid defined by a coordinate mapping from a logical grid to some physical
5 coordinates. The physical coordinates and the Jacobian are stored as grid
6 functions corresponding to the logical grid.
7
8 See also: [`logical_grid`](@ref), [`jacobian`](@ref), [`metric_tensor`](@ref).
9 """
10 struct MappedGrid{T,D, GT<:Grid{<:Any,D}, CT<:AbstractArray{T,D}, JT<:AbstractArray{<:AbstractMatrix{<:Any}, D}} <: Grid{T,D}
11 logical_grid::GT
12 physicalcoordinates::CT
13 jacobian::JT
14
15 """
16 MappedGrid(logical_grid, physicalcoordinates, jacobian)
17
18 A MappedGrid with the given physical coordinates and jacobian.
19 """
20 function MappedGrid(logical_grid::GT, physicalcoordinates::CT, jacobian::JT) where {T,D, GT<:Grid{<:Any,D}, CT<:AbstractArray{T,D}, JT<:AbstractArray{<:AbstractMatrix{<:Any}, D}}
21 if !(size(logical_grid) == size(physicalcoordinates) == size(jacobian))
22 throw(ArgumentError("Sizes must match"))
23 end
24
25 if size(first(jacobian)) != (length(first(physicalcoordinates)),D)
26 throw(ArgumentError("The size of the jacobian must match the dimensions of the grid and coordinates"))
27 end
28
29 return new{T,D,GT,CT,JT}(logical_grid, physicalcoordinates, jacobian)
30 end
31 end
32
33 function Base.:(==)(a::MappedGrid, b::MappedGrid)
34 same_logical_grid = logical_grid(a) == logical_grid(b)
35 same_coordinates = collect(a) == collect(b)
36 same_jacobian = jacobian(a) == jacobian(b)
37
38 return same_logical_grid && same_coordinates && same_jacobian
39 end
40
41 """
42 logical_grid(g::MappedGrid)
43
44 The logical grid of `g`.
45 """
46 logical_grid(g::MappedGrid) = g.logical_grid
47
48 """
49 jacobian(g::MappedGrid)
50
51 The Jacobian matrix of `g` as a grid function.
52 """
53 jacobian(g::MappedGrid) = g.jacobian
54
55
56 # Indexing interface
57 Base.getindex(g::MappedGrid, I::Vararg{Int}) = g.physicalcoordinates[I...]
58 Base.eachindex(g::MappedGrid) = eachindex(g.logical_grid)
59
60 Base.firstindex(g::MappedGrid, d) = firstindex(g.logical_grid, d)
61 Base.lastindex(g::MappedGrid, d) = lastindex(g.logical_grid, d)
62
63 # Iteration interface
64 Base.iterate(g::MappedGrid) = iterate(g.physicalcoordinates)
65 Base.iterate(g::MappedGrid, state) = iterate(g.physicalcoordinates, state)
66
67 Base.IteratorSize(::Type{<:MappedGrid{<:Any, D}}) where D = Base.HasShape{D}()
68 Base.length(g::MappedGrid) = length(g.logical_grid)
69 Base.size(g::MappedGrid) = size(g.logical_grid)
70 Base.size(g::MappedGrid, d) = size(g.logical_grid, d)
71
72 boundary_identifiers(g::MappedGrid) = boundary_identifiers(g.logical_grid)
73 boundary_indices(g::MappedGrid, id::TensorGridBoundary) = boundary_indices(g.logical_grid, id)
74
75 function boundary_grid(g::MappedGrid, id::TensorGridBoundary)
76 b_indices = boundary_indices(g.logical_grid, id)
77
78 # Calculate indices of needed jacobian components
79 D = ndims(g)
80 all_indices = SVector{D}(1:D)
81 free_variable_indices = deleteat(all_indices, grid_id(id))
82 jacobian_components = (:, free_variable_indices)
83
84 # Create grid function for boundary grid jacobian
85 boundary_jacobian = componentview((@view g.jacobian[b_indices...]) , jacobian_components...)
86 boundary_physicalcoordinates = @view g.physicalcoordinates[b_indices...]
87
88 return MappedGrid(
89 boundary_grid(g.logical_grid, id),
90 boundary_physicalcoordinates,
91 boundary_jacobian,
92 )
93 end
94
95
96 """
97 mapped_grid(x, J, size::Vararg{Int})
98
99 A `MappedGrid` with a default logical grid on the D-dimensional unit hyper
100 box [0,1]ᴰ. `x` and `J`are functions to be evaluated on the logical grid
101 and `size` determines the size of the logical grid.
102 """
103 function mapped_grid(x, J, size::Vararg{Int})
104 D = length(size)
105 lg = equidistant_grid(ntuple(i->0., D), ntuple(i->1., D), size...)
106 return mapped_grid(lg, x, J)
107 end
108
109 """
110 mapped_grid(lg::Grid, x, J)
111
112 A `MappedGrid` with logical grid `lg`. Physical coordinates and Jacobian are
113 determined by the functions `x` and `J`.
114 """
115 function mapped_grid(lg::Grid, x, J)
116 return MappedGrid(
117 lg,
118 map(x,lg),
119 map(J,lg),
120 )
121 end
122
123 """
124 metric_tensor(g::MappedGrid)
125
126 The metric tensor of `g` as a grid function.
127 """
128 function metric_tensor(g::MappedGrid)
129 return map(jacobian(g)) do ∂x∂ξ
130 ∂x∂ξ'*∂x∂ξ
131 end
132 end
133
134 function min_spacing(g::MappedGrid{T,1} where T)
135 n, = size(g)
136
137 ms = Inf
138 for i ∈ 1:n-1
139 ms = min(ms, norm(g[i+1]-g[i]))
140 end
141
142 return ms
143 end
144
145 function min_spacing(g::MappedGrid{T,2} where T)
146 n, m = size(g)
147
148 ms = Inf
149 for i ∈ 1:n-1, j ∈ 1:m-1 # loop over each cell of the grid
150
151 ms = min(
152 ms,
153 norm(g[i+1,j]-g[i,j]),
154 norm(g[i,j+1]-g[i,j]),
155
156 norm(g[i+1,j]-g[i+1,j+1]),
157 norm(g[i,j+1]-g[i+1,j+1]),
158
159 norm(g[i+1,j+1]-g[i,j]),
160 norm(g[i+1,j]-g[i,j+1]),
161 )
162 # NOTE: This could be optimized to avoid checking all interior edges twice.
163 end
164
165 return ms
166 end
167
168 """
169 normal(g::MappedGrid, boundary)
170
171 The outward pointing normal as a grid function on the corresponding boundary grid.
172 """
173 function normal(g::MappedGrid, boundary)
174 b_indices = boundary_indices(g, boundary)
175 σ = _boundary_sign(component_type(g), boundary)
176
177 # TODO: Refactor this when `boundary_indices(g, ...)` has been made iterable.
178 return map(jacobian(g)[b_indices...]) do ∂x∂ξ
179 ∂ξ∂x = inv(∂x∂ξ)
180 k = grid_id(boundary)
181 σ*∂ξ∂x[k,:]/norm(∂ξ∂x[k,:])
182 end
183 end
184
185 """
186 normal(g::MappedGrid, boundary, i...)
187
188 The outward pointing normal to the specified boundary in grid point `i`.
189 """
190 function normal(g::MappedGrid, boundary, i...)
191 σ = _boundary_sign(component_type(g), boundary)
192 ∂ξ∂x = inv(jacobian(g)[i...])
193
194 k = grid_id(boundary)
195 return σ*∂ξ∂x[k,:]/norm(∂ξ∂x[k,:])
196 end
197
198
199 function _boundary_sign(T, boundary)
200 if boundary_id(boundary) == UpperBoundary()
201 return one(T)
202 elseif boundary_id(boundary) == LowerBoundary()
203 return -one(T)
204 else
205 throw(ArgumentError("The boundary identifier must be either `LowerBoundary()` or `UpperBoundary()`"))
206 end
207 end