Mercurial > repos > public > sbplib_julia

comparison src/Grids/grid.jl @ 1858:4a9be96f2569 feature/documenter_logo

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 Sun, 12 Jan 2025 21:18:44 +0100
parents 863385aae454
children 03894fd7b132
comparison
equal deleted inserted replaced
1857:ffde7dad9da5 1858:4a9be96f2569
1 """
2 Grid{T,D}
3
4 A grid with coordinates of type `T`, e.g. `SVector{3,Float64}`, and dimension
5 `D`. The grid can be embedded in a higher dimension in which case the number
6 of indices and the number of components of the coordinate vectors will be
7 different.
8
9 All grids are expected to behave as a grid function for the coordinates.
10
11 `Grids` is top level abstract type for grids. A grid should implement Julia's interfaces for
12 indexing and iteration.
13
14 ## Note
15
16 Importantly a grid does not have to be an `AbstractArray`. The reason is to
17 allow flexible handling of special types of grids like multi-block grids, or
18 grids with special indexing.
19 """
20 abstract type Grid{T,D} end
21
22 Base.ndims(::Grid{T,D}) where {T,D} = D
23 Base.eltype(::Type{<:Grid{T}}) where T = T
24
25 Base.getindex(g::Grid, I::CartesianIndex) = g[Tuple(I)...]
26
27 """
28 coordinate_size(g)
29
30 The length of the coordinate vector of `Grid` `g`.
31 """
32 coordinate_size(::Type{<:Grid{T}}) where T = _ncomponents(T)
33 coordinate_size(g::Grid) = coordinate_size(typeof(g)) # TBD: Name of this function?!
34
35 """
36 component_type(gf)
37
38 The type of the components of the elements of `gf`. I.e if `gf` is a vector
39 valued grid function, `component_view(gf)` is the element type of the vectors
40 at each grid point.
41
42 # Examples
43 ```julia-repl
44 julia> component_type([[1,2], [2,3], [3,4]])
45 Int64
46 ```
47 """
48 component_type(T::Type) = eltype(eltype(T))
49 component_type(t) = component_type(typeof(t))
50
51 """
52 componentview(gf, component_index...)
53
54 A view of `gf` with only the components specified by `component_index...`.
55
56 # Examples
57 ```julia-repl
58 julia> componentview([[1,2], [2,3], [3,4]],2)
59 3-element ArrayComponentView{Int64, Vector{Int64}, 1, Vector{Vector{Int64}}, Tuple{Int64}}:
60 2
61 3
62 4
63 ```
64 """
65 componentview(gf, component_index...) = ArrayComponentView(gf, component_index)
66
67 struct ArrayComponentView{CT,T,D,AT <: AbstractArray{T,D}, IT} <: AbstractArray{CT,D}
68 v::AT
69 component_index::IT
70
71 function ArrayComponentView(v, component_index)
72 CT = typeof(first(v)[component_index...])
73 return new{CT, eltype(v), ndims(v), typeof(v), typeof(component_index)}(v,component_index)
74 end
75 end
76
77 Base.size(cv::ArrayComponentView) = size(cv.v)
78 Base.getindex(cv::ArrayComponentView, i::Int) = cv.v[i][cv.component_index...]
79 Base.getindex(cv::ArrayComponentView, I::Vararg{Int}) = cv.v[I...][cv.component_index...]
80 IndexStyle(::Type{<:ArrayComponentView{<:Any,<:Any,AT}}) where AT = IndexStyle(AT)
81
82 # TODO: Implement `setindex!`?
83 # TODO: Implement a more general ComponentView that can handle non-AbstractArrays.
84
85
86 """
87 min_spacing(g::Grid)
88
89 The smallest distance between any pair of grid points in `g`.
90 """
91 function min_spacing end
92
93 """
94 refine(g::Grid, r)
95
96 The grid where `g` is refined by the factor `r`.
97
98 See also: [`coarsen`](@ref).
99 """
100 function refine end
101
102 """
103 coarsen(g::Grid, r)
104
105 The grid where `g` is coarsened by the factor `r`.
106
107 See also: [`refine`](@ref).
108 """
109 function coarsen end
110
111 """
112 BoundaryIdentifier
113
114 An identifier for a boundary of a grid.
115 """
116 abstract type BoundaryIdentifier end
117
118 """
119 boundary_identifiers(g::Grid)
120
121 Identifiers for all the boundaries of `g`.
122 """
123 function boundary_identifiers end
124
125 """
126 boundary_grid(g::Grid, id::BoundaryIdentifier)
127
128 The grid for the boundary specified by `id`.
129 """
130 function boundary_grid end
131 # TBD: Can we implement a version here that accepts multiple ids and grouped boundaries? Maybe we need multiblock stuff?
132
133 """
134 boundary_indices(g::Grid, id::BoundaryIdentifier)
135
136 A collection of indices corresponding to the boundary with given id. For grids
137 with Cartesian indexing these collections will be tuples with elements of type
138 ``Union{Int,Colon}``.
139
140 When implementing this method it is expected that the returned collection can
141 be used to index grid functions to obtain grid functions on the boundary grid.
142 """
143 function boundary_indices end
144
145 """
146 eval_on(g::Grid, f)
147
148 Lazy evaluation of `f` on the grid. `f` can either be on the form `f(x,y,...)`
149 with each coordinate as an argument, or on the form `f(x̄)` taking a
150 coordinate vector.
151
152 For concrete array grid functions `map(f,g)` can be used instead.
153 """
154 eval_on(g::Grid, f) = eval_on(g, f, Base.IteratorSize(g))
155 function eval_on(g::Grid, f, ::Base.HasShape)
156 if hasmethod(f, (Any,))
157 return LazyTensors.LazyFunctionArray((I...)->f(g[I...]), size(g))
158 else
159 # TBD This branch can be removed if we accept the trade off that we define f with the syntax f((x,y)) instead if we don't want to handle the vector in the body of f. (Add an example in the docs)
160 # Also see Notes.md
161 return LazyTensors.LazyFunctionArray((I...)->f(g[I...]...), size(g))
162 end
163 end
164
165 """
166 eval_on(g::Grid, f::Number)
167
168 Lazy evaluation of a scalar `f` on the grid.
169 """
170 eval_on(g::Grid, f::Number) = return LazyTensors.LazyConstantArray(f, size(g))
171
172 _ncomponents(::Type{<:Number}) = 1
173 _ncomponents(T::Type{<:SVector}) = length(T)
174
175