Mercurial > repos > public > sbplib_julia

comparison src/Grids/mapped_grid.jl @ 1751:f3d7e2d7a43f feature/sbp_operators/laplace_curvilinear

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