Mercurial > repos > public > sbplib_julia

comparison src/Grids/mapped_grid.jl @ 2057:8a2a0d678d6f feature/lazy_tensors/pretty_printing

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 Tue, 10 Feb 2026 22:41:19 +0100
parents e97f4352b8d0
children 83759a0d0f28
comparison
equal deleted inserted replaced
1110:c0bff9f6e0fb 2057:8a2a0d678d6f
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 mapped_grid(x, J, size...)
97
98 A `MappedGrid` with a default logical grid on the D-dimensional unit hyper
99 box [0,1]ᴰ. `x` and `J` are functions to be evaluated on the logical grid
100 and `size` determines the size of the logical grid.
101 """
102 function mapped_grid(x, J, size::Vararg{Int})
103 D = length(size)
104 return mapped_grid(x, J, unithyperbox(D), size...)
105 end
106
107 """
108 mapped_grid(x, J, lg::Grid)
109
110 A `MappedGrid` with logical grid `lg`. Physical coordinates and Jacobian are
111 determined by the functions `x` and `J`.
112 """
113 function mapped_grid(x, J, lg::Grid)
114 return MappedGrid(
115 lg,
116 map(x,lg),
117 map(J,lg),
118 )
119 end
120
121 """
122 mapped_grid(x, J, ps::ParameterSpace, size...)
123
124 A `MappedGrid` with logical grid `lg`. Physical coordinates and Jacobian are
125 determined by the functions `x` and `J`.
126 """
127 function mapped_grid(x, J, ps::ParameterSpace, size::Vararg{Int})
128 lg = equidistant_grid(ps, size...)
129 return mapped_grid(x, J, lg)
130 end
131
132 """
133 metric_tensor(g::MappedGrid)
134
135 The metric tensor of `g` as a grid function.
136 """
137 function metric_tensor(g::MappedGrid)
138 return map(jacobian(g)) do ∂x∂ξ
139 ∂x∂ξ'*∂x∂ξ
140 end
141 end
142
143 function min_spacing(g::MappedGrid{T,1} where T)
144 n, = size(g)
145
146 ms = Inf
147 for i ∈ 1:n-1
148 ms = min(ms, norm(g[i+1]-g[i]))
149 end
150
151 return ms
152 end
153
154 function min_spacing(g::MappedGrid{T,2} where T)
155 n, m = size(g)
156
157 ms = Inf
158 for i ∈ 1:n-1, j ∈ 1:m-1 # loop over each cell of the grid
159
160 ms = min(
161 ms,
162 norm(g[i+1,j]-g[i,j]),
163 norm(g[i,j+1]-g[i,j]),
164
165 norm(g[i+1,j]-g[i+1,j+1]),
166 norm(g[i,j+1]-g[i+1,j+1]),
167
168 norm(g[i+1,j+1]-g[i,j]),
169 norm(g[i+1,j]-g[i,j+1]),
170 )
171 # NOTE: This could be optimized to avoid checking all interior edges twice.
172 end
173
174 return ms
175 end
176
177 """
178 normal(g::MappedGrid, boundary)
179
180 The outward pointing normal as a grid function on the corresponding boundary grid.
181 """
182 function normal(g::MappedGrid, boundary)
183 return map(boundary_indices(g, boundary)) do I
184 normal(g, boundary, Tuple(I)...)
185 end
186 end
187
188 """
189 normal(g::MappedGrid, boundary, i...)
190
191 The outward pointing normal to the specified boundary in grid point `i`.
192 """
193 function normal(g::MappedGrid, boundary, i...)
194 σ = _boundary_sign(component_type(g), boundary)
195 ∂ξ∂x = inv(jacobian(g)[i...])
196
197 k = grid_id(boundary)
198 return σ*∂ξ∂x[k,:]/norm(∂ξ∂x[k,:])
199 end
200
201
202 function _boundary_sign(T, boundary)
203 if boundary_id(boundary) == UpperBoundary()
204 return one(T)
205 elseif boundary_id(boundary) == LowerBoundary()
206 return -one(T)
207 else
208 throw(ArgumentError("The boundary identifier must be either `LowerBoundary()` or `UpperBoundary()`"))
209 end
210 end