Mercurial > repos > public > sbplib_julia

comparison src/Grids/mapped_grid.jl @ 1801:2b5f81e288f1 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, 02 Oct 2024 08:51:37 +0200
parents 819ab806960f 8583f6379bd8
children 1c58005429fd
comparison
equal deleted inserted replaced
1784:c5070edd0ebb 1801:2b5f81e288f1
70 Base.size(g::MappedGrid, d) = size(g.logical_grid, d) 70 Base.size(g::MappedGrid, d) = size(g.logical_grid, d)
71 71
72 boundary_identifiers(g::MappedGrid) = boundary_identifiers(g.logical_grid) 72 boundary_identifiers(g::MappedGrid) = boundary_identifiers(g.logical_grid)
73 boundary_indices(g::MappedGrid, id::TensorGridBoundary) = boundary_indices(g.logical_grid, id) 73 boundary_indices(g::MappedGrid, id::TensorGridBoundary) = boundary_indices(g.logical_grid, id)
74 74
75 # Review: Error when calling plot(boundary_grid(g, id))
76 # Currently need to collect first, i.e., plot(collect(boundary_grid(g, id)))
77 function boundary_grid(g::MappedGrid, id::TensorGridBoundary) 75 function boundary_grid(g::MappedGrid, id::TensorGridBoundary)
78 b_indices = boundary_indices(g.logical_grid, id) 76 b_indices = boundary_indices(g.logical_grid, id)
79 77
80 # Calculate indices of needed jacobian components 78 # Calculate indices of needed jacobian components
81 D = ndims(g) 79 D = ndims(g)
96 94
97 95
98 """ 96 """
99 mapped_grid(x, J, size::Vararg{Int}) 97 mapped_grid(x, J, size::Vararg{Int})
100 98
101 A `MappedGrid` with a default logical grid on a unit hyper box. `x` and `J` 99 A `MappedGrid` with a default logical grid on the D-dimensional unit hyper
102 are functions to be evaluated on the logical grid and `size` determines the 100 box [0,1]ᴰ. `x` and `J`are functions to be evaluated on the logical grid
103 size of the logical grid. 101 and `size` determines the size of the logical grid.
104 """ 102 """
105 function mapped_grid(x, J, size::Vararg{Int}) 103 function mapped_grid(x, J, size::Vararg{Int})
106 D = length(size) 104 D = length(size)
107 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...)
108 return mapped_grid(lg, x, J) 106 return mapped_grid(lg, x, J)
131 return map(jacobian(g)) do ∂x∂ξ 129 return map(jacobian(g)) do ∂x∂ξ
132 ∂x∂ξ'*∂x∂ξ 130 ∂x∂ξ'*∂x∂ξ
133 end 131 end
134 end 132 end
135 133
136 """
137 metric_tensor_inverse(g::MappedGrid)
138
139 The inverse of the metric tensor of `g` as a grid function.
140 """
141 function metric_tensor_inverse(g::MappedGrid)
142 return map(jacobian(g)) do ∂x∂ξ
143 inv(∂x∂ξ'*∂x∂ξ)
144 end
145 end
146
147 function min_spacing(g::MappedGrid{T,1} where T) 134 function min_spacing(g::MappedGrid{T,1} where T)
148 n, = size(g) 135 n, = size(g)
149 136
150 ms = Inf 137 ms = Inf
151 for i ∈ 1:n-1 138 for i ∈ 1:n-1
183 170
184 The outward pointing normal as a grid function on the corresponding boundary grid. 171 The outward pointing normal as a grid function on the corresponding boundary grid.
185 """ 172 """
186 function normal(g::MappedGrid, boundary) 173 function normal(g::MappedGrid, boundary)
187 b_indices = boundary_indices(g, boundary) 174 b_indices = boundary_indices(g, boundary)
188 σ =_boundary_sign(component_type(g), boundary) 175 σ = _boundary_sign(component_type(g), boundary)
189 return map(jacobian(g)[b_indices...]) do ∂x∂ξ 176 return map(jacobian(g)[b_indices...]) do ∂x∂ξ
190 ∂ξ∂x = inv(∂x∂ξ) 177 ∂ξ∂x = inv(∂x∂ξ)
191 k = grid_id(boundary) 178 k = grid_id(boundary)
192 σ*∂ξ∂x[k,:]/norm(∂ξ∂x[k,:]) 179 σ*∂ξ∂x[k,:]/norm(∂ξ∂x[k,:])
193 end 180 end
194 end 181 end
182
183 function normal(g::MappedGrid, boundary, i...)
184 σ = _boundary_sign(component_type(g), boundary)
185 ∂ξ∂x = inv(jacobian(g)[i...])
186
187 k = grid_id(boundary)
188 return σ*∂ξ∂x[k,:]/norm(∂ξ∂x[k,:])
189 end
190
195 191
196 function _boundary_sign(T, boundary) 192 function _boundary_sign(T, boundary)
197 if boundary_id(boundary) == UpperBoundary() 193 if boundary_id(boundary) == UpperBoundary()
198 return one(T) 194 return one(T)
199 elseif boundary_id(boundary) == LowerBoundary() 195 elseif boundary_id(boundary) == LowerBoundary()