Mercurial > repos > public > sbplib_julia

comparison src/Grids/mapped_grid.jl @ 1777:819ab806960f 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 Sun, 15 Sep 2024 23:00:15 +0200
parents e2abd72d7ce0 265a740a49da
children 2b5f81e288f1
comparison
equal deleted inserted replaced
1750:e2abd72d7ce0 1777:819ab806960f
3 3
4 A grid defined by a coordinate mapping from a logical grid to some physical 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 5 coordinates. The physical coordinates and the Jacobian are stored as grid
6 functions corresponding to the logical grid. 6 functions corresponding to the logical grid.
7 7
8 See also: [`logicalgrid`](@ref), [`jacobian`](@ref), [`metric_tensor`](@ref). 8 See also: [`logical_grid`](@ref), [`jacobian`](@ref), [`metric_tensor`](@ref).
9 """ 9 """
10 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{<:AbstractMatrix{<:Any}, D}} <: Grid{T,D}
11 logicalgrid::GT 11 logical_grid::GT
12 physicalcoordinates::CT 12 physicalcoordinates::CT
13 jacobian::JT 13 jacobian::JT
14 14
15 """ 15 """
16 MappedGrid(logicalgrid, physicalcoordinates, jacobian) 16 MappedGrid(logical_grid, physicalcoordinates, jacobian)
17 17
18 A MappedGrid with the given physical coordinates and jacobian. 18 A MappedGrid with the given physical coordinates and jacobian.
19 """ 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}} 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(logicalgrid) == size(physicalcoordinates) == size(jacobian)) 21 if !(size(logical_grid) == size(physicalcoordinates) == size(jacobian))
22 throw(ArgumentError("Sizes must match")) 22 throw(ArgumentError("Sizes must match"))
23 end 23 end
24 24
25 if size(first(jacobian)) != (length(first(physicalcoordinates)),D) 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")) 26 throw(ArgumentError("The size of the jacobian must match the dimensions of the grid and coordinates"))
27 end 27 end
28 28
29 return new{T,D,GT,CT,JT}(logicalgrid, physicalcoordinates, jacobian) 29 return new{T,D,GT,CT,JT}(logical_grid, physicalcoordinates, jacobian)
30 end 30 end
31 end 31 end
32 32
33 function Base.:(==)(a::MappedGrid, b::MappedGrid) 33 function Base.:(==)(a::MappedGrid, b::MappedGrid)
34 same_logicalgrid = logicalgrid(a) == logicalgrid(b) 34 same_logical_grid = logical_grid(a) == logical_grid(b)
35 same_coordinates = collect(a) == collect(b) 35 same_coordinates = collect(a) == collect(b)
36 same_jacobian = jacobian(a) == jacobian(b) 36 same_jacobian = jacobian(a) == jacobian(b)
37 37
38 return same_logicalgrid && same_coordinates && same_jacobian 38 return same_logical_grid && same_coordinates && same_jacobian
39 end 39 end
40 40
41 """ 41 """
42 logicalgrid(g::MappedGrid) 42 logical_grid(g::MappedGrid)
43 43
44 The logical grid of `g`. 44 The logical grid of `g`.
45 """ 45 """
46 logicalgrid(g::MappedGrid) = g.logicalgrid 46 logical_grid(g::MappedGrid) = g.logical_grid
47 47
48 """ 48 """
49 jacobian(g::MappedGrid) 49 jacobian(g::MappedGrid)
50 50
51 The Jacobian matrix of `g` as a grid function. 51 The Jacobian matrix of `g` as a grid function.
53 jacobian(g::MappedGrid) = g.jacobian 53 jacobian(g::MappedGrid) = g.jacobian
54 54
55 55
56 # Indexing interface 56 # Indexing interface
57 Base.getindex(g::MappedGrid, I::Vararg{Int}) = g.physicalcoordinates[I...] 57 Base.getindex(g::MappedGrid, I::Vararg{Int}) = g.physicalcoordinates[I...]
58 Base.eachindex(g::MappedGrid) = eachindex(g.logicalgrid) 58 Base.eachindex(g::MappedGrid) = eachindex(g.logical_grid)
59 59
60 Base.firstindex(g::MappedGrid, d) = firstindex(g.logicalgrid, d) 60 Base.firstindex(g::MappedGrid, d) = firstindex(g.logical_grid, d)
61 Base.lastindex(g::MappedGrid, d) = lastindex(g.logicalgrid, d) 61 Base.lastindex(g::MappedGrid, d) = lastindex(g.logical_grid, d)
62 62
63 # Iteration interface 63 # Iteration interface
64 Base.iterate(g::MappedGrid) = iterate(g.physicalcoordinates) 64 Base.iterate(g::MappedGrid) = iterate(g.physicalcoordinates)
65 Base.iterate(g::MappedGrid, state) = iterate(g.physicalcoordinates, state) 65 Base.iterate(g::MappedGrid, state) = iterate(g.physicalcoordinates, state)
66 66
67 Base.IteratorSize(::Type{<:MappedGrid{<:Any, D}}) where D = Base.HasShape{D}() 67 Base.IteratorSize(::Type{<:MappedGrid{<:Any, D}}) where D = Base.HasShape{D}()
68 Base.length(g::MappedGrid) = length(g.logicalgrid) 68 Base.length(g::MappedGrid) = length(g.logical_grid)
69 Base.size(g::MappedGrid) = size(g.logicalgrid) 69 Base.size(g::MappedGrid) = size(g.logical_grid)
70 Base.size(g::MappedGrid, d) = size(g.logicalgrid, d) 70 Base.size(g::MappedGrid, d) = size(g.logical_grid, d)
71 71
72 boundary_identifiers(g::MappedGrid) = boundary_identifiers(g.logicalgrid) 72 boundary_identifiers(g::MappedGrid) = boundary_identifiers(g.logical_grid)
73 boundary_indices(g::MappedGrid, id::TensorGridBoundary) = boundary_indices(g.logicalgrid, 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)))
75 function boundary_grid(g::MappedGrid, id::TensorGridBoundary) 77 function boundary_grid(g::MappedGrid, id::TensorGridBoundary)
76 b_indices = boundary_indices(g.logicalgrid, id) 78 b_indices = boundary_indices(g.logical_grid, id)
77 79
78 # Calculate indices of needed jacobian components 80 # Calculate indices of needed jacobian components
79 D = ndims(g) 81 D = ndims(g)
80 all_indices = SVector{D}(1:D) 82 all_indices = SVector{D}(1:D)
81 free_variable_indices = deleteat(all_indices, grid_id(id)) 83 free_variable_indices = deleteat(all_indices, grid_id(id))
84 # Create grid function for boundary grid jacobian 86 # Create grid function for boundary grid jacobian
85 boundary_jacobian = componentview((@view g.jacobian[b_indices...]) , jacobian_components...) 87 boundary_jacobian = componentview((@view g.jacobian[b_indices...]) , jacobian_components...)
86 boundary_physicalcoordinates = @view g.physicalcoordinates[b_indices...] 88 boundary_physicalcoordinates = @view g.physicalcoordinates[b_indices...]
87 89
88 return MappedGrid( 90 return MappedGrid(
89 boundary_grid(g.logicalgrid, id), 91 boundary_grid(g.logical_grid, id),
90 boundary_physicalcoordinates, 92 boundary_physicalcoordinates,
91 boundary_jacobian, 93 boundary_jacobian,
92 ) 94 )
93 end 95 end
94 96
117 lg, 119 lg,
118 map(x,lg), 120 map(x,lg),
119 map(J,lg), 121 map(J,lg),
120 ) 122 )
121 end 123 end
122
123 """
124 jacobian_determinant(g::MappedGrid)
125
126 The jacobian determinant of `g` as a grid function.
127 """
128 function jacobian_determinant(g::MappedGrid)
129 return map(jacobian(g)) do ∂x∂ξ
130 det(∂x∂ξ)
131 end
132 end
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 124
138 """ 125 """
139 metric_tensor(g::MappedGrid) 126 metric_tensor(g::MappedGrid)
140 127
141 The metric tensor of `g` as a grid function. 128 The metric tensor of `g` as a grid function.