changeset 1500:0d630d4e1d21 feature/grids/componentview

Merge default
author Jonatan Werpers <jonatan@werpers.com>
date Thu, 25 Jan 2024 10:54:00 +0100
parents 06dec7b68f82 (diff) 8ef207e4bc87 (current diff)
children de7b76b61ecd
files
diffstat 3 files changed, 79 insertions(+), 8 deletions(-) [+]
line wrap: on
line diff
--- a/src/Grids/Grids.jl	Tue Jan 23 20:39:06 2024 +0100
+++ b/src/Grids/Grids.jl	Thu Jan 25 10:54:00 2024 +0100
@@ -13,10 +13,11 @@
 export boundary_indices
 export boundary_identifiers
 export boundary_grid
+export coarsen
+export refine
 export eval_on
-export coarsen
-export getcomponent
-export refine
+export componentview
+export ArrayComponentView
 
 export BoundaryIdentifier
 export TensorGridBoundary
--- a/src/Grids/grid.jl	Tue Jan 23 20:39:06 2024 +0100
+++ b/src/Grids/grid.jl	Thu Jan 25 10:54:00 2024 +0100
@@ -33,12 +33,48 @@
 coordinate_size(g::Grid) = coordinate_size(typeof(g)) # TBD: Name of this function?!
 
 """
-    component_type(g)
+    component_type(gf)
+
+The type of the components of `gf`.
+"""
+component_type(T::Type) = eltype(eltype(T))
+component_type(t) = component_type(typeof(t))
+
+componentview(gf, component_index...) = ArrayComponentView(gf, component_index)
+
+# REVIEW: Should this only be defined for vector-valued component types of the same dimension?
+# Now one can for instance do:  v = [rand(2,2),rand(2,2), rand(2,1)] and cv = componentview(v,1,2)
+# resulting in #undef in the third component of cv.
+struct ArrayComponentView{CT,T,D,AT <: AbstractArray{T,D}, IT} <: AbstractArray{CT,D}
+    v::AT
+    component_index::IT
 
-The type of the components of the coordinate vector of `Grid` `g`.
-"""
-component_type(::Type{<:Grid{T}}) where T = eltype(T)
-component_type(g::Grid) = component_type(typeof(g))
+    function ArrayComponentView(v, component_index)
+        CT = typeof(first(v)[component_index...])
+        return new{CT, eltype(v), ndims(v), typeof(v), typeof(component_index)}(v,component_index)
+    end
+end
+
+Base.size(cv) = size(cv.v)
+Base.getindex(cv::ArrayComponentView, i::Int) = cv.v[i][cv.component_index...]
+Base.getindex(cv::ArrayComponentView, I::Vararg{Int}) = cv.v[I...][cv.component_index...] #REVIEW: Will this allocate if I... slices v? if so, we should probably use a view on v?
+IndexStyle(::Type{<:ArrayComponentView{<:Any,<:Any,AT}}) where AT = IndexStyle(AT)
+
+# TODO: Implement the remaining optional methods from the array interface
+# setindex!(A, v, i::Int)
+# setindex!(A, v, I::Vararg{Int, N})
+# iterate
+# length(A)
+# similar(A)
+# similar(A, ::Type{S})
+# similar(A, dims::Dims)
+# similar(A, ::Type{S}, dims::Dims)
+# # Non-traditional indices
+# axes(A)
+# similar(A, ::Type{S}, inds)
+# similar(T::Union{Type,Function}, inds)
+
+# TODO: Implement a more general ComponentView that can handle non-AbstractArrays.
 
 """
     refine(g::Grid, r)
--- a/test/Grids/grid_test.jl	Tue Jan 23 20:39:06 2024 +0100
+++ b/test/Grids/grid_test.jl	Thu Jan 25 10:54:00 2024 +0100
@@ -63,6 +63,40 @@
     @test eval_on(g, f) == map(x̄->f(x̄...), g)
 end
 
+@testset "componentview" begin
+    # REVIEW: I think we can reduce the index ranges.
+    v = [@SMatrix[1 3; 2 4] .+ 100*i .+ 10*j for i ∈ 1:3, j∈ 1:4]
+
+    @test componentview(v, 1, 1) == [1 .+ 100*i .+ 10*j for i ∈ 1:3, j∈ 1:4]
+    @test componentview(v, 1, 2) == [3 .+ 100*i .+ 10*j for i ∈ 1:3, j∈ 1:4]
+    @test componentview(v, 2, 1) == [2 .+ 100*i .+ 10*j for i ∈ 1:3, j∈ 1:4]
+
+    @test componentview(v, 1, :) == [@SVector[1,3] .+ 100*i .+ 10*j for i ∈ 1:3, j∈ 1:4]
+    @test componentview(v, 2, :) == [@SVector[2,4] .+ 100*i .+ 10*j for i ∈ 1:3, j∈ 1:4]
+    @test componentview(v, :, 1) == [@SVector[1,2] .+ 100*i .+ 10*j for i ∈ 1:3, j∈ 1:4]
+    @test componentview(v, :, 2) == [@SVector[3,4] .+ 100*i .+ 10*j for i ∈ 1:3, j∈ 1:4]
+
+
+    A = @SMatrix[
+        1 4 7;
+        2 5 8;
+        3 6 9;
+    ]
+    v = [A .+ 100*i .+ 10*j for i ∈ 1:3, j∈ 1:4]
+    @test componentview(v, 1:2, 1:2) == [@SMatrix[1 4;2 5] .+ 100*i .+ 10*j for i ∈ 1:3, j∈ 1:4]
+    @test componentview(v, 2:3, 1:2) == [@SMatrix[2 5;3 6] .+ 100*i .+ 10*j for i ∈ 1:3, j∈ 1:4]
+end
+
+@testset "ArrayComponentView" begin
+    v = [@SMatrix[1 3; 2 4] .+ 100*i .+ 10*j for i ∈ 1:3, j∈ 1:4]
+
+    @test ArrayComponentView(v, (1,1)) == ArrayComponentView(v, (1,1)) # REVIEW: Does not test anything?
+    @test ArrayComponentView(v, (1,1)) == ArrayComponentView(copy(v), (1,1))
+    # REVIEW: The two below are equivalent?
+    @test ArrayComponentView(v, (1,1)) == [1 .+ 100*i .+ 10*j for i ∈ 1:3, j∈ 1:4]
+    @test [1 .+ 100*i .+ 10*j for i ∈ 1:3, j∈ 1:4] == ArrayComponentView(v, (1,1))
+end
+
 @testset "_ncomponents" begin
     @test Grids._ncomponents(Int) == 1
     @test Grids._ncomponents(Float64) == 1