comparison src/Grids/grid.jl @ 1531:9da4ab4fb85e bugfix/sbp_operators/stencil_return_type

Merge default
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 11 Apr 2024 22:50:42 +0200
parents d641798539c2
children 9113f437431d 611ae2308aa1
comparison
equal deleted inserted replaced
1456:f13857f37b8f 1531:9da4ab4fb85e
25 Base.getindex(g::Grid, I::CartesianIndex) = g[Tuple(I)...] 25 Base.getindex(g::Grid, I::CartesianIndex) = g[Tuple(I)...]
26 26
27 """ 27 """
28 coordinate_size(g) 28 coordinate_size(g)
29 29
30 The lenght of the coordinate vector of `Grid` `g`. 30 The length of the coordinate vector of `Grid` `g`.
31 """ 31 """
32 coordinate_size(::Type{<:Grid{T}}) where T = _ncomponents(T) 32 coordinate_size(::Type{<:Grid{T}}) where T = _ncomponents(T)
33 coordinate_size(g::Grid) = coordinate_size(typeof(g)) # TBD: Name of this function?! 33 coordinate_size(g::Grid) = coordinate_size(typeof(g)) # TBD: Name of this function?!
34 34
35 """ 35 """
36 component_type(g) 36 component_type(gf)
37 37
38 The type of the components of the coordinate vector of `Grid` `g`. 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 ```
39 """ 47 """
40 component_type(::Type{<:Grid{T}}) where T = eltype(T) 48 component_type(T::Type) = eltype(eltype(T))
41 component_type(g::Grid) = component_type(typeof(g)) 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) = 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.
42 84
43 """ 85 """
44 refine(g::Grid, r) 86 refine(g::Grid, r)
45 87
46 The grid where `g` is refined by the factor `r`. 88 The grid where `g` is refined by the factor `r`.
72 """ 114 """
73 function boundary_grid end 115 function boundary_grid end
74 # TBD: Can we implement a version here that accepts multiple ids and grouped boundaries? Maybe we need multiblock stuff? 116 # TBD: Can we implement a version here that accepts multiple ids and grouped boundaries? Maybe we need multiblock stuff?
75 117
76 """ 118 """
119 boundary_indices(g::Grid, id::BoundaryIdentifier)
120
121 A collection of indices corresponding to the boundary with given id. For grids
122 with Cartesian indexing these collections will be tuples with elements of type
123 ``Union{Int,Colon}``.
124
125 When implementing this method it is expected that the returned collection can
126 be used to index grid functions to obtain grid functions on the boundary grid.
127 """
128 function boundary_indices end
129
130 """
77 eval_on(g::Grid, f) 131 eval_on(g::Grid, f)
78 132
79 Lazy evaluation of `f` on the grid. `f` can either be on the form `f(x,y,...)` 133 Lazy evaluation of `f` on the grid. `f` can either be on the form `f(x,y,...)`
80 with each coordinate as an argument, or on the form `f(x̄)` taking a 134 with each coordinate as an argument, or on the form `f(x̄)` taking a
81 coordinate vector. 135 coordinate vector.
85 eval_on(g::Grid, f) = eval_on(g, f, Base.IteratorSize(g)) 139 eval_on(g::Grid, f) = eval_on(g, f, Base.IteratorSize(g))
86 function eval_on(g::Grid, f, ::Base.HasShape) 140 function eval_on(g::Grid, f, ::Base.HasShape)
87 if hasmethod(f, (Any,)) 141 if hasmethod(f, (Any,))
88 return LazyTensors.LazyFunctionArray((I...)->f(g[I...]), size(g)) 142 return LazyTensors.LazyFunctionArray((I...)->f(g[I...]), size(g))
89 else 143 else
144 # 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)
145 # Also see Notes.md
90 return LazyTensors.LazyFunctionArray((I...)->f(g[I...]...), size(g)) 146 return LazyTensors.LazyFunctionArray((I...)->f(g[I...]...), size(g))
91 end 147 end
92 end 148 end
93 149
94 """ 150 """
98 """ 154 """
99 eval_on(g::Grid, f::Number) = return LazyTensors.LazyConstantArray(f, size(g)) 155 eval_on(g::Grid, f::Number) = return LazyTensors.LazyConstantArray(f, size(g))
100 156
101 _ncomponents(::Type{<:Number}) = 1 157 _ncomponents(::Type{<:Number}) = 1
102 _ncomponents(T::Type{<:SVector}) = length(T) 158 _ncomponents(T::Type{<:SVector}) = length(T)
159
160