Mercurial > repos > public > sbplib_julia

comparison src/Grids/equidistant_grid.jl @ 1222:5f677cd6f0b6 refactor/grids

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Start refactoring
author Jonatan Werpers <jonatan@werpers.com>
date Sat, 18 Feb 2023 11:37:35 +0100
parents 50b008d2e937
children 2abec782cf5b
comparison
equal deleted inserted replaced
1221:b3b4d29b46c3 1222:5f677cd6f0b6
1 struct EquidistantGrid{T,R<:AbstractRange{T}} <: Grid{T,1,1}
2 points::R
3 end
4
5 Base.eltype(g::EquidistantGrid{T}) where T = T
6 Base.getindex(g::EquidistantGrid, i) = g.points[i]
7 Base.size(g::EquidistantGrid) = size(g.points)
8 Base.length(g::EquidistantGrid) = length(g.points)
9 Base.eachindex(g::EquidistantGrid) = eachindex(g.points)
10
11 # TODO: Make sure collect works!
12
13
1 """ 14 """
2 EquidistantGrid{Dim,T<:Real} <: Grid 15 spacing(grid::EquidistantGrid)
3 16
4 `Dim`-dimensional equidistant grid with coordinates of type `T`. 17 The spacing between grid points.
5 """ 18 """
6 struct EquidistantGrid{Dim,T<:Real} <: Grid 19 spacing(g::EquidistantGrid) = step(g.points)
7 size::NTuple{Dim, Int}
8 limit_lower::NTuple{Dim, T}
9 limit_upper::NTuple{Dim, T}
10 20
11 function EquidistantGrid{Dim,T}(size::NTuple{Dim, Int}, limit_lower::NTuple{Dim, T}, limit_upper::NTuple{Dim, T}) where {Dim,T} 21
12 if any(size .<= 0) 22 """
13 throw(DomainError("all components of size must be postive")) 23 inverse_spacing(grid::EquidistantGrid)
14 end 24
15 if any(limit_upper.-limit_lower .<= 0) 25 The reciprocal of the spacing between grid points.
16 throw(DomainError("all side lengths must be postive")) 26 """
17 end 27 inverse_spacing(g::EquidistantGrid) = 1/step(g.points)
18 return new{Dim,T}(size, limit_lower, limit_upper) 28
19 end 29
30 boundary_identifiers(::EquidistantGrid) = (Lower(), Upper())
31 boundary_grid(g::EquidistantGrid, id::Lower) = ZeroDimGrid(g[begin])
32 boundary_grid(g::EquidistantGrid, id::Upper) = ZeroDimGrid(g[end])
33
34
35 """
36 refine(g::EquidistantGrid, r::Int)
37
38 Refines `grid` by a factor `r`. The factor is applied to the number of
39 intervals which is 1 less than the size of the grid.
40
41 See also: [`coarsen`](@ref)
42 """
43 function refine(g::EquidistantGrid, r::Int)
44 new_sz = (length(g) - 1)*r + 1
45 return EquidistantGrid(change_length(g.points, new_sz))
20 end 46 end
21 47
22 """ 48 """
23 EquidistantGrid(size, limit_lower, limit_upper) 49 coarsen(grid::EquidistantGrid, r::Int)
50
51 Coarsens `grid` by a factor `r`. The factor is applied to the number of
52 intervals which is 1 less than the size of the grid. If the number of
53 intervals are not divisible by `r` an error is raised.
54
55 See also: [`refine`](@ref)
56 """
57 function coarsen(g::EquidistantGrid, r::Int)
58 if (length(g)-1)%r != 0
59 throw(DomainError(r, "Size minus 1 must be divisible by the ratio."))
60 end
61
62 new_sz = (length(g) - 1)÷r + 1
63
64 return EquidistantGrid(change_length(g.points), new_sz)
65 end
66
67
68
69
70
71
72
73 """
74 equidistant_grid(size::Dims, limit_lower, limit_upper)
24 75
25 Construct an equidistant grid with corners at the coordinates `limit_lower` and 76 Construct an equidistant grid with corners at the coordinates `limit_lower` and
26 `limit_upper`. 77 `limit_upper`.
27 78
28 The length of the domain sides are given by the components of 79 The length of the domain sides are given by the components of
30 as `(-1,1)x(0,2)`. The side lengths of the grid are not allowed to be negative. 81 as `(-1,1)x(0,2)`. The side lengths of the grid are not allowed to be negative.
31 82
32 The number of equidistantly spaced points in each coordinate direction are given 83 The number of equidistantly spaced points in each coordinate direction are given
33 by the tuple `size`. 84 by the tuple `size`.
34 """ 85 """
35 function EquidistantGrid(size, limit_lower, limit_upper) 86 function equidistant_grid(size::Dims, limit_lower, limit_upper)
36 return EquidistantGrid{length(size), eltype(limit_lower)}(size, limit_lower, limit_upper) 87 gs = map(size, limit_lower, limit_upper) do s,l,u
37 end 88 EquidistantGrid(range(l, u, length=s)) # TBD: Should it use LinRange instead?
89 end
38 90
39 """ 91 return TensorGrid(gs...)
40 EquidistantGrid{T}()
41
42 Constructs a 0-dimensional grid.
43 """
44 EquidistantGrid{T}() where T = EquidistantGrid{0,T}((),(),()) # Convenience constructor for 0-dim grid
45
46
47 """
48 EquidistantGrid(size::Int, limit_lower::T, limit_upper::T)
49
50 Convenience constructor for 1D grids.
51 """
52 function EquidistantGrid(size::Int, limit_lower::T, limit_upper::T) where T
53 return EquidistantGrid((size,),(limit_lower,),(limit_upper,))
54 end
55
56 Base.eltype(grid::EquidistantGrid{Dim,T}) where {Dim,T} = T
57
58 Base.eachindex(grid::EquidistantGrid) = CartesianIndices(grid.size)
59
60 Base.size(g::EquidistantGrid) = g.size
61
62 Base.ndims(::EquidistantGrid{Dim}) where Dim = Dim
63
64 function Base.getindex(g::EquidistantGrid, I::Vararg{Int})
65 h = spacing(g)
66 return g.limit_lower .+ (I.-1).*h
67 end
68
69 Base.getindex(g::EquidistantGrid, I::CartesianIndex) = g[Tuple(I)...]
70
71 # Review:
72 # Is it not strange that evalOn(::Grid) is non-lazy while evalOn(::EquidistantGrid) is?
73 # Also: Change name to evalon or eval_on!!!!!!
74 function evalOn(grid::EquidistantGrid, f::Function)
75 F(I...) = f(grid[I...]...)
76
77 return LazyFunctionArray(F, size(grid))
78 end
79
80 """
81 spacing(grid::EquidistantGrid)
82
83 The spacing between grid points.
84 """
85 spacing(grid::EquidistantGrid) = (grid.limit_upper.-grid.limit_lower)./(grid.size.-1)
86
87
88 """
89 inverse_spacing(grid::EquidistantGrid)
90
91 The reciprocal of the spacing between grid points.
92 """
93 inverse_spacing(grid::EquidistantGrid) = 1 ./ spacing(grid)
94
95
96 """
97 points(grid::EquidistantGrid)
98
99 The point of the grid as an array of tuples with the same dimension as the grid.
100 The points are stored as [(x1,y1), (x1,y2), … (x1,yn);
101 (x2,y1), (x2,y2), … (x2,yn);
102 ⋮ ⋮ ⋮
103 (xm,y1), (xm,y2), … (xm,yn)]
104 """
105 function points(grid::EquidistantGrid)
106 indices = Tuple.(CartesianIndices(grid.size))
107 h = spacing(grid)
108 return broadcast(I -> grid.limit_lower .+ (I.-1).*h, indices)
109 end
110
111 """
112 restrict(::EquidistantGrid, dim)
113
114 Pick out given dimensions from the grid and return a grid for them.
115 """
116 function restrict(grid::EquidistantGrid, dim)
117 size = grid.size[dim]
118 limit_lower = grid.limit_lower[dim]
119 limit_upper = grid.limit_upper[dim]
120
121 return EquidistantGrid(size, limit_lower, limit_upper)
122 end 92 end
123 93
124 94
125 """ 95 """
126 orthogonal_dims(grid::EquidistantGrid,dim) 96 equidistant_grid(size::Int, limit_lower::T, limit_upper::T)
127 97
128 Returns the dimensions of grid orthogonal to that of dim. 98 Constructs a 1D equidistant grid.
129 """ 99 """
130 function orthogonal_dims(grid::EquidistantGrid, dim) 100 function equidistant_grid(size::Int, limit_lower::T, limit_upper::T) where T
131 orth_dims = filter(i -> i != dim, dims(grid)) 101 return equidistant_grid((size,),(limit_lower,),(limit_upper,))
132 if orth_dims == dims(grid)
133 throw(DomainError(string("dimension ",string(dim)," not matching grid")))
134 end
135 return orth_dims
136 end 102 end
137 103
138 104
139 """
140 boundary_identifiers(::EquidistantGrid)
141
142 Returns a tuple containing the boundary identifiers for the grid, stored as
143 (CartesianBoundary(1,Lower),
144 CartesianBoundary(1,Upper),
145 CartesianBoundary(2,Lower),
146 ...)
147 """
148 boundary_identifiers(g::EquidistantGrid) = (((ntuple(i->(CartesianBoundary{i,Lower}(),CartesianBoundary{i,Upper}()),ndims(g)))...)...,)
149
150 105
151 """ 106 """
152 boundary_grid(grid::EquidistantGrid, id::CartesianBoundary) 107 change_length(::AbstractRange, n)
153 108
154 Creates the lower-dimensional restriciton of `grid` spanned by the dimensions 109 Change the length of a range to `n`, keeping the same start and stop.
155 orthogonal to the boundary specified by `id`. The boundary grid of a 1-dimensional
156 grid is a zero-dimensional grid.
157 """ 110 """
158 function boundary_grid(grid::EquidistantGrid, id::CartesianBoundary) 111 function change_length(::AbstractRange, n) end
159 orth_dims = orthogonal_dims(grid, dim(id))
160 return restrict(grid, orth_dims)
161 end
162 boundary_grid(::EquidistantGrid{1,T},::CartesianBoundary{1}) where T = EquidistantGrid{T}()
163 112
164 113 change_length(r::LinRange, n) = LinRange(r[begin], r[end], n)
165 """ 114 change_length(r::StepRangeLen, n) = range(r[begin], r[end], n)
166 refine(grid::EquidistantGrid, r::Int) 115 # TODO: Test the above
167
168 Refines `grid` by a factor `r`. The factor is applied to the number of
169 intervals which is 1 less than the size of the grid.
170
171 See also: [`coarsen`](@ref)
172 """
173 function refine(grid::EquidistantGrid, r::Int)
174 sz = size(grid)
175 new_sz = (sz .- 1).*r .+ 1
176 return EquidistantGrid{ndims(grid), eltype(grid)}(new_sz, grid.limit_lower, grid.limit_upper)
177 end
178
179
180 """
181 coarsen(grid::EquidistantGrid, r::Int)
182
183 Coarsens `grid` by a factor `r`. The factor is applied to the number of
184 intervals which is 1 less than the size of the grid. If the number of
185 intervals are not divisible by `r` an error is raised.
186
187 See also: [`refine`](@ref)
188 """
189 function coarsen(grid::EquidistantGrid, r::Int)
190 sz = size(grid)
191
192 if !all(n -> (n % r == 0), sz.-1)
193 throw(DomainError(r, "Size minus 1 must be divisible by the ratio."))
194 end
195
196 new_sz = (sz .- 1).÷r .+ 1
197
198 return EquidistantGrid{ndims(grid), eltype(grid)}(new_sz, grid.limit_lower, grid.limit_upper)
199 end