Mercurial > repos > public > sbplib_julia

comparison src/Grids/mapped_grid.jl @ 1748:03894fd7b132 feature/grids/manifolds

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge feature/grids/curvilinear
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 11 Sep 2024 15:41:58 +0200
parents a4c52ae93b11 44faa334adc6
children e2abd72d7ce0
comparison
equal deleted inserted replaced
1695:a4c52ae93b11 1748:03894fd7b132
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: [`logicalgrid`](@ref), [`jacobian`](@ref), [`metric_tensor`](@ref).
9 """
1 struct MappedGrid{T,D, GT<:Grid{<:Any,D}, CT<:AbstractArray{T,D}, JT<:AbstractArray{<:AbstractArray{<:Any, 2}, D}} <: Grid{T,D} 10 struct MappedGrid{T,D, GT<:Grid{<:Any,D}, CT<:AbstractArray{T,D}, JT<:AbstractArray{<:AbstractArray{<:Any, 2}, D}} <: Grid{T,D}
2 logicalgrid::GT 11 logicalgrid::GT
3 physicalcoordinates::CT 12 physicalcoordinates::CT
4 jacobian::JT 13 jacobian::JT
5 end 14
6 15 """
16 MappedGrid(logicalgrid, physicalcoordinates, jacobian)
17
18 A MappedGrid with the given physical coordinates and jacobian.
19 """
20 function MappedGrid(logicalgrid::GT, physicalcoordinates::CT, jacobian::JT) where {T,D, GT<:Grid{<:Any,D}, CT<:AbstractArray{T,D}, JT<:AbstractArray{<:AbstractArray{<:Any, 2}, D}}
21 if !(size(logicalgrid) == size(physicalcoordinates) == size(jacobian))
22 throw(ArgumentError("Sizes must match"))
23 end
24
25 if size(first(jacobian)) != (length(first(physicalcoordinates)),D)
26 @show size(first(jacobian))
27 @show (length(first(physicalcoordinates)),D)
28 throw(ArgumentError("The size of the jacobian must match the dimensions of the grid and coordinates"))
29 end
30
31 return new{T,D,GT,CT,JT}(logicalgrid, physicalcoordinates, jacobian)
32 end
33 end
34
35 function Base.:(==)(a::MappedGrid, b::MappedGrid)
36 same_logicalgrid = logicalgrid(a) == logicalgrid(b)
37 same_coordinates = collect(a) == collect(b)
38 same_jacobian = jacobian(a) == jacobian(b)
39
40 return same_logicalgrid && same_coordinates && same_jacobian
41 end
42
43 """
44 logicalgrid(g::MappedGrid)
45
46 The logical grid of `g`.
47 """
48 logicalgrid(g::MappedGrid) = g.logicalgrid
49
50 """
51 jacobian(g::MappedGrid)
52
53 The Jacobian matrix of `g` as a grid function.
54 """
7 jacobian(g::MappedGrid) = g.jacobian 55 jacobian(g::MappedGrid) = g.jacobian
8 logicalgrid(g::MappedGrid) = g.logicalgrid
9 56
10 57
11 # Indexing interface 58 # Indexing interface
12 Base.getindex(g::MappedGrid, I::Vararg{Int}) = g.physicalcoordinates[I...] 59 Base.getindex(g::MappedGrid, I::Vararg{Int}) = g.physicalcoordinates[I...]
13 Base.eachindex(g::MappedGrid) = eachindex(g.logicalgrid) 60 Base.eachindex(g::MappedGrid) = eachindex(g.logicalgrid)
14 61
15 Base.firstindex(g::MappedGrid, d) = firstindex(g.logicalgrid, d) 62 Base.firstindex(g::MappedGrid, d) = firstindex(g.logicalgrid, d)
16 Base.lastindex(g::MappedGrid, d) = lastindex(g.logicalgrid, d) 63 Base.lastindex(g::MappedGrid, d) = lastindex(g.logicalgrid, d)
17 # TODO: axes(...)?
18 64
19 # Iteration interface 65 # Iteration interface
20
21 Base.iterate(g::MappedGrid) = iterate(g.physicalcoordinates) 66 Base.iterate(g::MappedGrid) = iterate(g.physicalcoordinates)
22 Base.iterate(g::MappedGrid, state) = iterate(g.physicalcoordinates, state) 67 Base.iterate(g::MappedGrid, state) = iterate(g.physicalcoordinates, state)
23 68
24 Base.IteratorSize(::Type{<:MappedGrid{<:Any, D}}) where D = Base.HasShape{D}() 69 Base.IteratorSize(::Type{<:MappedGrid{<:Any, D}}) where D = Base.HasShape{D}()
25 Base.length(g::MappedGrid) = length(g.logicalgrid) 70 Base.length(g::MappedGrid) = length(g.logicalgrid)
47 boundary_physicalcoordinates, 92 boundary_physicalcoordinates,
48 boundary_jacobian, 93 boundary_jacobian,
49 ) 94 )
50 end 95 end
51 96
52 # TBD: refine and coarsen could be implemented once we have a simple manifold implementation. 97
53 # Before we do, we should consider the overhead of including such a field in the mapped grid struct. 98 """
54 99 mapped_grid(x, J, size::Vararg{Int})
55 function mapped_grid(x, J, size...) 100
101 A `MappedGrid` with a default logical grid on a unit hyper box. `x` and `J`
102 are functions to be evaluated on the logical grid and `size` determines the
103 size of the logical grid.
104 """
105 function mapped_grid(x, J, size::Vararg{Int})
56 D = length(size) 106 D = length(size)
57 lg = equidistant_grid(ntuple(i->0., D), ntuple(i->1., D), size...) 107 lg = equidistant_grid(ntuple(i->0., D), ntuple(i->1., D), size...)
108 return mapped_grid(lg, x, J)
109 end
110
111 """
112 mapped_grid(lg::Grid, x, J)
113
114 A `MappedGrid` with logical grid `lg`. Physical coordinates and Jacobian are
115 determined by the functions `x` and `J`.
116 """
117 function mapped_grid(lg::Grid, x, J)
58 return MappedGrid( 118 return MappedGrid(
59 lg, 119 lg,
60 map(x,lg), 120 map(x,lg),
61 map(J,lg), 121 map(J,lg),
62 ) 122 )
63 end 123 end
64 # TODO: Delete this function 124
65 125 """
66 126 jacobian_determinant(g::MappedGrid)
127
128 The jacobian determinant of `g` as a grid function.
129 """
67 function jacobian_determinant(g::MappedGrid) 130 function jacobian_determinant(g::MappedGrid)
68 return map(jacobian(g)) do ∂x∂ξ 131 return map(jacobian(g)) do ∂x∂ξ
69 det(∂x∂ξ) 132 det(∂x∂ξ)
70 end 133 end
71 end 134 end
72 135 # TBD: Should this be changed to calculate sqrt(g) instead?
136 # This would make it well defined also for n-dim grids embedded in higher dimensions.
137 # TBD: Is there a better name? metric_determinant?
138 # TBD: Is the best option to delete it?
139
140 """
141 metric_tensor(g::MappedGrid)
142
143 The metric tensor of `g` as a grid function.
144 """
73 function metric_tensor(g::MappedGrid) 145 function metric_tensor(g::MappedGrid)
74 return map(jacobian(g)) do ∂x∂ξ 146 return map(jacobian(g)) do ∂x∂ξ
75 ∂x∂ξ'*∂x∂ξ 147 ∂x∂ξ'*∂x∂ξ
76 end 148 end
77 end 149 end
78 150
151 """
152 metric_tensor_inverse(g::MappedGrid)
153
154 The inverse of the metric tensor of `g` as a grid function.
155 """
79 function metric_tensor_inverse(g::MappedGrid) 156 function metric_tensor_inverse(g::MappedGrid)
80 return map(jacobian(g)) do ∂x∂ξ 157 return map(jacobian(g)) do ∂x∂ξ
81 inv(∂x∂ξ'*∂x∂ξ) 158 inv(∂x∂ξ'*∂x∂ξ)
82 end 159 end
83 end 160 end
117 end 194 end
118 195
119 """ 196 """
120 normal(g::MappedGrid, boundary) 197 normal(g::MappedGrid, boundary)
121 198
122 The outward pointing normal as a grid function on the boundary 199 The outward pointing normal as a grid function on the corresponding boundary grid.
123 """ 200 """
124 function normal(g::MappedGrid, boundary) 201 function normal(g::MappedGrid, boundary)
125 b_indices = boundary_indices(g, boundary) 202 b_indices = boundary_indices(g, boundary)
126 σ =_boundary_sign(component_type(g), boundary) 203 σ =_boundary_sign(component_type(g), boundary)
127 return map(jacobian(g)[b_indices...]) do ∂x∂ξ 204 return map(jacobian(g)[b_indices...]) do ∂x∂ξ
130 σ*∂ξ∂x[k,:]/norm(∂ξ∂x[k,:]) 207 σ*∂ξ∂x[k,:]/norm(∂ξ∂x[k,:])
131 end 208 end
132 end 209 end
133 210
134 function _boundary_sign(T, boundary) 211 function _boundary_sign(T, boundary)
135 if boundary_id(boundary) == Upper() 212 if boundary_id(boundary) == UpperBoundary()
136 return one(T) 213 return one(T)
137 elseif boundary_id(boundary) == Lower() 214 elseif boundary_id(boundary) == LowerBoundary()
138 return -one(T) 215 return -one(T)
139 else 216 else
140 throw(ArgumentError("The boundary identifier must be either `Lower()` or `Upper()`")) 217 throw(ArgumentError("The boundary identifier must be either `LowerBoundary()` or `UpperBoundary()`"))
141 end 218 end
142 end 219 end