Mercurial > repos > public > sbplib_julia

comparison src/Grids/equidistant_grid.jl @ 1365:4684c7f1c4cb feature/variable_derivatives

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge with default
author Vidar Stiernström <vidar.stiernstrom@it.uu.se>
date Sun, 21 May 2023 21:55:14 +0200
parents 102ebdaf7c11 08f06bfacd5c
children d2219cc8316b
comparison
equal deleted inserted replaced
1358:e7861cfb6ede 1365:4684c7f1c4cb
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 Base.size(g::EquidistantGrid) = size(g.points)
29 30
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 function Base.getindex(g::EquidistantGrid, I::Vararg{Int})
66 h = spacing(g)
67 return g.limit_lower .+ (I.-1).*h
68 end
69
70 Base.getindex(g::EquidistantGrid, I::CartesianIndex) = g[Tuple(I)...]
71 # TBD: Can this method be removed if `EquidistantGrid` is an AbstractArray?
72
73 Base.ndims(::EquidistantGrid{Dim}) where Dim = Dim
74 31
75 """ 32 """
76 spacing(grid::EquidistantGrid) 33 spacing(grid::EquidistantGrid)
77 34
78 The spacing between grid points. 35 The spacing between grid points.
79 """ 36 """
80 spacing(grid::EquidistantGrid) = (grid.limit_upper.-grid.limit_lower)./(grid.size.-1) 37 spacing(g::EquidistantGrid) = step(g.points)
81 38
82 39
83 """ 40 """
84 inverse_spacing(grid::EquidistantGrid) 41 inverse_spacing(grid::EquidistantGrid)
85 42
86 The reciprocal of the spacing between grid points. 43 The reciprocal of the spacing between grid points.
87 """ 44 """
88 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])
89 51
90 52
91 """ 53 """
92 points(grid::EquidistantGrid) 54 refine(g::EquidistantGrid, r::Int)
93 55
94 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
95 The points are stored as [(x1,y1), (x1,y2), … (x1,yn); 57 intervals, i.e., 1 less than the size of `g`.
96 (x2,y1), (x2,y2), … (x2,yn); 58
97 ⋮ ⋮ ⋮ 59 See also: [`coarsen`](@ref)
98 (xm,y1), (xm,y2), … (xm,yn)]
99 """ 60 """
100 function points(grid::EquidistantGrid) 61 function refine(g::EquidistantGrid, r::Int)
101 indices = Tuple.(CartesianIndices(grid.size)) 62 new_sz = (length(g) - 1)*r + 1
102 h = spacing(grid) 63 return EquidistantGrid(change_length(g.points, new_sz))
103 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))
104 end 83 end
105 84
106 85
107 """ 86 """
108 restrict(::EquidistantGrid, dim) 87 equidistant_grid(size::Dims, limit_lower, limit_upper)
109 88
110 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.
111 """ 104 """
112 function restrict(grid::EquidistantGrid, dim) 105 function equidistant_grid(size::Dims, limit_lower, limit_upper)
113 size = grid.size[dim] 106 gs = map(equidistant_grid, size, limit_lower, limit_upper)
114 limit_lower = grid.limit_lower[dim] 107 return TensorGrid(gs...)
115 limit_upper = grid.limit_upper[dim] 108 end
116 109
117 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?
118 end 124 end
125
126 CartesianBoundary{D,BID} = TensorGridBoundary{D,BID} # TBD: What should we do about the naming of this boundary?
119 127
120 128
121 """ 129 """
122 orthogonal_dims(grid::EquidistantGrid,dim) 130 change_length(r::AbstractRange, n)
123 131
124 Returns the dimensions of grid orthogonal to that of dim. 132 Change the length of `r` to `n`, keeping the same start and stop.
125 """ 133 """
126 function orthogonal_dims(grid::EquidistantGrid, dim) 134 function change_length end
127 orth_dims = filter(i -> i != dim, dims(grid))
128 if orth_dims == dims(grid)
129 throw(DomainError(string("dimension ",string(dim)," not matching grid")))
130 end
131 return orth_dims
132 end
133 135
134 136 change_length(r::UnitRange, n) = StepRange{Int,Int}(range(r[begin], r[end], n))
135 """ 137 change_length(r::StepRange, n) = StepRange{Int,Int}(range(r[begin], r[end], n))
136 boundary_identifiers(::EquidistantGrid) 138 change_length(r::StepRangeLen, n) = range(r[begin], r[end], n)
137 139 change_length(r::LinRange, n) = LinRange(r[begin], r[end], n)
138 Returns a tuple containing the boundary identifiers for the grid, stored as
139 (CartesianBoundary(1,Lower),
140 CartesianBoundary(1,Upper),
141 CartesianBoundary(2,Lower),
142 ...)
143 """
144 boundary_identifiers(g::EquidistantGrid) = (((ntuple(i->(CartesianBoundary{i,Lower}(),CartesianBoundary{i,Upper}()),ndims(g)))...)...,)
145
146
147 """
148 boundary_grid(grid::EquidistantGrid, id::CartesianBoundary)
149
150 Creates the lower-dimensional restriciton of `grid` spanned by the dimensions
151 orthogonal to the boundary specified by `id`. The boundary grid of a 1-dimensional
152 grid is a zero-dimensional grid.
153 """
154 function boundary_grid(grid::EquidistantGrid, id::CartesianBoundary)
155 orth_dims = orthogonal_dims(grid, dim(id))
156 return restrict(grid, orth_dims)
157 end
158 boundary_grid(::EquidistantGrid{1,T},::CartesianBoundary{1}) where T = EquidistantGrid{T}()
159
160
161 """
162 refine(grid::EquidistantGrid, r::Int)
163
164 Refines `grid` by a factor `r`. The factor is applied to the number of
165 intervals which is 1 less than the size of the grid.
166
167 See also: [`coarsen`](@ref)
168 """
169 function refine(grid::EquidistantGrid, r::Int)
170 sz = size(grid)
171 new_sz = (sz .- 1).*r .+ 1
172 return EquidistantGrid{ndims(grid), eltype(grid)}(new_sz, grid.limit_lower, grid.limit_upper)
173 end
174
175
176 """
177 coarsen(grid::EquidistantGrid, r::Int)
178
179 Coarsens `grid` by a factor `r`. The factor is applied to the number of
180 intervals which is 1 less than the size of the grid. If the number of
181 intervals are not divisible by `r` an error is raised.
182
183 See also: [`refine`](@ref)
184 """
185 function coarsen(grid::EquidistantGrid, r::Int)
186 sz = size(grid)
187
188 if !all(n -> (n % r == 0), sz.-1)
189 throw(DomainError(r, "Size minus 1 must be divisible by the ratio."))
190 end
191
192 new_sz = (sz .- 1).÷r .+ 1
193
194 return EquidistantGrid{ndims(grid), eltype(grid)}(new_sz, grid.limit_lower, grid.limit_upper)
195 end