Mercurial > repos > public > sbplib_julia

comparison src/Grids/equidistant_grid.jl @ 1360:f59228534d3a tooling/benchmarks

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge default
author Jonatan Werpers <jonatan@werpers.com>
date Sat, 20 May 2023 15:15:22 +0200
parents 08f06bfacd5c
children 4684c7f1c4cb
comparison
equal deleted inserted replaced
1321:42738616422e 1360:f59228534d3a
1 """
2 EquidistantGrid{T,R<:AbstractRange{T}} <: Grid{T,1}
1 3
4 A one-dimensional equidistant grid. Most users are expected to use
5 [`equidistant_grid`](@ref) for constructing equidistant grids.
6
7 See also: [`equidistant_grid`](@ref)
8
9
10 ## Note
11 The type of range used for the points can likely impact performance.
2 """ 12 """
3 EquidistantGrid{Dim,T<:Real} <: Grid 13 struct EquidistantGrid{T,R<:AbstractRange{T}} <: Grid{T,1}
4 14 points::R
5 `Dim`-dimensional equidistant grid with coordinates of type `T`.
6 """
7 struct EquidistantGrid{Dim,T<:Real} <: Grid
8 size::NTuple{Dim, Int}
9 limit_lower::NTuple{Dim, T}
10 limit_upper::NTuple{Dim, T}
11
12 function EquidistantGrid{Dim,T}(size::NTuple{Dim, Int}, limit_lower::NTuple{Dim, T}, limit_upper::NTuple{Dim, T}) where {Dim,T}
13 if any(size .<= 0)
14 throw(DomainError("all components of size must be postive"))
15 end
16 if any(limit_upper.-limit_lower .<= 0)
17 throw(DomainError("all side lengths must be postive"))
18 end
19 return new{Dim,T}(size, limit_lower, limit_upper)
20 end
21 end 15 end
22 16
17 # Indexing interface
18 Base.getindex(g::EquidistantGrid, i) = g.points[i]
19 Base.eachindex(g::EquidistantGrid) = eachindex(g.points)
20 Base.firstindex(g::EquidistantGrid) = firstindex(g.points)
21 Base.lastindex(g::EquidistantGrid) = lastindex(g.points)
23 22
24 """ 23 # Iteration interface
25 EquidistantGrid(size, limit_lower, limit_upper) 24 Base.iterate(g::EquidistantGrid) = iterate(g.points)
25 Base.iterate(g::EquidistantGrid, state) = iterate(g.points, state)
26 26
27 Construct an equidistant grid with corners at the coordinates `limit_lower` and 27 Base.IteratorSize(::Type{<:EquidistantGrid}) = Base.HasShape{1}()
28 `limit_upper`. 28 Base.length(g::EquidistantGrid) = length(g.points)
29 29 Base.size(g::EquidistantGrid) = size(g.points)
30 The length of the domain sides are given by the components of
31 `limit_upper-limit_lower`. E.g for a 2D grid with `limit_lower=(-1,0)` and `limit_upper=(1,2)` the domain is defined
32 as `(-1,1)x(0,2)`. The side lengths of the grid are not allowed to be negative.
33
34 The number of equidistantly spaced points in each coordinate direction are given
35 by the tuple `size`.
36 """
37 function EquidistantGrid(size, limit_lower, limit_upper)
38 return EquidistantGrid{length(size), eltype(limit_lower)}(size, limit_lower, limit_upper)
39 end
40
41
42 """
43 EquidistantGrid{T}()
44
45 Constructs a 0-dimensional grid.
46 """
47 EquidistantGrid{T}() where T = EquidistantGrid{0,T}((),(),()) # Convenience constructor for 0-dim grid
48
49
50 """
51 EquidistantGrid(size::Int, limit_lower::T, limit_upper::T)
52
53 Convenience constructor for 1D grids.
54 """
55 function EquidistantGrid(size::Int, limit_lower::T, limit_upper::T) where T
56 return EquidistantGrid((size,),(limit_lower,),(limit_upper,))
57 end
58
59 Base.eltype(grid::EquidistantGrid{Dim,T}) where {Dim,T} = T
60
61 Base.eachindex(grid::EquidistantGrid) = CartesianIndices(grid.size)
62
63 Base.size(g::EquidistantGrid) = g.size
64
65 Base.ndims(::EquidistantGrid{Dim}) where Dim = Dim
66
67
68
69 30
70 31
71 """ 32 """
72 spacing(grid::EquidistantGrid) 33 spacing(grid::EquidistantGrid)
73 34
74 The spacing between grid points. 35 The spacing between grid points.
75 """ 36 """
76 spacing(grid::EquidistantGrid) = (grid.limit_upper.-grid.limit_lower)./(grid.size.-1) 37 spacing(g::EquidistantGrid) = step(g.points)
77 38
78 39
79 """ 40 """
80 inverse_spacing(grid::EquidistantGrid) 41 inverse_spacing(grid::EquidistantGrid)
81 42
82 The reciprocal of the spacing between grid points. 43 The reciprocal of the spacing between grid points.
83 """ 44 """
84 inverse_spacing(grid::EquidistantGrid) = 1 ./ spacing(grid) 45 inverse_spacing(g::EquidistantGrid) = 1/step(g.points)
46
47
48 boundary_identifiers(::EquidistantGrid) = (Lower(), Upper())
49 boundary_grid(g::EquidistantGrid, id::Lower) = ZeroDimGrid(g[begin])
50 boundary_grid(g::EquidistantGrid, id::Upper) = ZeroDimGrid(g[end])
85 51
86 52
87 """ 53 """
88 points(grid::EquidistantGrid) 54 refine(g::EquidistantGrid, r::Int)
89 55
90 The point of the grid as an array of tuples with the same dimension as the grid. 56 The grid where `g` is refined by the factor `r`. The factor is applied to the number of
91 The points are stored as [(x1,y1), (x1,y2), … (x1,yn); 57 intervals, i.e., 1 less than the size of `g`.
92 (x2,y1), (x2,y2), … (x2,yn); 58
93 ⋮ ⋮ ⋮ 59 See also: [`coarsen`](@ref)
94 (xm,y1), (xm,y2), … (xm,yn)]
95 """ 60 """
96 function points(grid::EquidistantGrid) 61 function refine(g::EquidistantGrid, r::Int)
97 indices = Tuple.(CartesianIndices(grid.size)) 62 new_sz = (length(g) - 1)*r + 1
98 h = spacing(grid) 63 return EquidistantGrid(change_length(g.points, new_sz))
99 return broadcast(I -> grid.limit_lower .+ (I.-1).*h, indices) 64 end
65
66 """
67 coarsen(g::EquidistantGrid, r::Int)
68
69 The grid where `g` is coarsened by the factor `r`. The factor is applied to the number of
70 intervals, i.e., 1 less than the size of `g`. If the number of
71 intervals are not divisible by `r` an error is raised.
72
73 See also: [`refine`](@ref)
74 """
75 function coarsen(g::EquidistantGrid, r::Int)
76 if (length(g)-1)%r != 0
77 throw(DomainError(r, "Size minus 1 must be divisible by the ratio."))
78 end
79
80 new_sz = (length(g) - 1)÷r + 1
81
82 return EquidistantGrid(change_length(g.points, new_sz))
100 end 83 end
101 84
102 85
103 """ 86 """
104 restrict(::EquidistantGrid, dim) 87 equidistant_grid(size::Dims, limit_lower, limit_upper)
105 88
106 Pick out given dimensions from the grid and return a grid for them. 89 Construct an equidistant grid with corners at the coordinates `limit_lower` and
90 `limit_upper`.
91
92 The length of the domain sides are given by the components of
93 `limit_upper-limit_lower`. E.g for a 2D grid with `limit_lower=(-1,0)` and
94 `limit_upper=(1,2)` the domain is defined as `(-1,1)x(0,2)`. The side lengths
95 of the grid are not allowed to be negative.
96
97 The number of equispaced points in each coordinate direction are given
98 by the tuple `size`.
99
100 Note: If `limit_lower` and `limit_upper` are integers and `size` would allow a
101 completely integer grid, `equidistant_grid` will still return a floating point
102 grid. This simlifies the implementation and avoids certain surprise
103 behaviours.
107 """ 104 """
108 function restrict(grid::EquidistantGrid, dim) 105 function equidistant_grid(size::Dims, limit_lower, limit_upper)
109 size = grid.size[dim] 106 gs = map(equidistant_grid, size, limit_lower, limit_upper)
110 limit_lower = grid.limit_lower[dim] 107 return TensorGrid(gs...)
111 limit_upper = grid.limit_upper[dim] 108 end
112 109
113 return EquidistantGrid(size, limit_lower, limit_upper) 110 """
111 equidistant_grid(size::Int, limit_lower::T, limit_upper::T)
112
113 Constructs a 1D equidistant grid.
114 """
115 function equidistant_grid(size::Int, limit_lower::T, limit_upper::T) where T
116 if any(size .<= 0)
117 throw(DomainError("size must be postive"))
118 end
119
120 if any(limit_upper.-limit_lower .<= 0)
121 throw(DomainError("side length must be postive"))
122 end
123 return EquidistantGrid(range(limit_lower, limit_upper, length=size)) # TBD: Should it use LinRange instead?
114 end 124 end
125
126 CartesianBoundary{D,BID} = TensorGridBoundary{D,BID} # TBD: What should we do about the naming of this boundary?
115 127
116 128
117 """ 129 """
118 orthogonal_dims(grid::EquidistantGrid,dim) 130 change_length(r::AbstractRange, n)
119 131
120 Returns the dimensions of grid orthogonal to that of dim. 132 Change the length of `r` to `n`, keeping the same start and stop.
121 """ 133 """
122 function orthogonal_dims(grid::EquidistantGrid, dim) 134 function change_length end
123 orth_dims = filter(i -> i != dim, dims(grid))
124 if orth_dims == dims(grid)
125 throw(DomainError(string("dimension ",string(dim)," not matching grid")))
126 end
127 return orth_dims
128 end
129 135
130 136 change_length(r::UnitRange, n) = StepRange{Int,Int}(range(r[begin], r[end], n))
131 """ 137 change_length(r::StepRange, n) = StepRange{Int,Int}(range(r[begin], r[end], n))
132 boundary_identifiers(::EquidistantGrid) 138 change_length(r::StepRangeLen, n) = range(r[begin], r[end], n)
133 139 change_length(r::LinRange, n) = LinRange(r[begin], r[end], n)
134 Returns a tuple containing the boundary identifiers for the grid, stored as
135 (CartesianBoundary(1,Lower),
136 CartesianBoundary(1,Upper),
137 CartesianBoundary(2,Lower),
138 ...)
139 """
140 boundary_identifiers(g::EquidistantGrid) = (((ntuple(i->(CartesianBoundary{i,Lower}(),CartesianBoundary{i,Upper}()),ndims(g)))...)...,)
141
142
143 """
144 boundary_grid(grid::EquidistantGrid, id::CartesianBoundary)
145
146 Creates the lower-dimensional restriciton of `grid` spanned by the dimensions
147 orthogonal to the boundary specified by `id`. The boundary grid of a 1-dimensional
148 grid is a zero-dimensional grid.
149 """
150 function boundary_grid(grid::EquidistantGrid, id::CartesianBoundary)
151 orth_dims = orthogonal_dims(grid, dim(id))
152 return restrict(grid, orth_dims)
153 end
154 boundary_grid(::EquidistantGrid{1,T},::CartesianBoundary{1}) where T = EquidistantGrid{T}()
155
156
157 """
158 refine(grid::EquidistantGrid, r::Int)
159
160 Refines `grid` by a factor `r`. The factor is applied to the number of
161 intervals which is 1 less than the size of the grid.
162
163 See also: [`coarsen`](@ref)
164 """
165 function refine(grid::EquidistantGrid, r::Int)
166 sz = size(grid)
167 new_sz = (sz .- 1).*r .+ 1
168 return EquidistantGrid{ndims(grid), eltype(grid)}(new_sz, grid.limit_lower, grid.limit_upper)
169 end
170
171
172 """
173 coarsen(grid::EquidistantGrid, r::Int)
174
175 Coarsens `grid` by a factor `r`. The factor is applied to the number of
176 intervals which is 1 less than the size of the grid. If the number of
177 intervals are not divisible by `r` an error is raised.
178
179 See also: [`refine`](@ref)
180 """
181 function coarsen(grid::EquidistantGrid, r::Int)
182 sz = size(grid)
183
184 if !all(n -> (n % r == 0), sz.-1)
185 throw(DomainError(r, "Size minus 1 must be divisible by the ratio."))
186 end
187
188 new_sz = (sz .- 1).÷r .+ 1
189
190 return EquidistantGrid{ndims(grid), eltype(grid)}(new_sz, grid.limit_lower, grid.limit_upper)
191 end