Mercurial > repos > public > sbplib_julia

changeset 1480:4550beef9694 feature/boundary_conditions

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 Sat, 23 Dec 2023 23:03:13 +0100
parents b96858a50e35 (current diff) 869a40cda18c (diff)
children ee242c3fe4af
files Notes.md
diffstat 9 files changed, 114 insertions(+), 25 deletions(-) [+]
line wrap: on
line diff
--- a/Notes.md	Sat Dec 23 23:01:28 2023 +0100
+++ b/Notes.md	Sat Dec 23 23:03:13 2023 +0100
@@ -331,3 +331,15 @@
 Could the implementation of LazyOuterProduct be simplified by making it a
 struct containing two or more LazyTensors? (using split_tuple in a similar way
 as TensorGrid)
+
+## Implementation of boundary_indices for more complex grids
+To represent boundaries of for example tet-elements we can use a type `IndexCollection` to index a grid function directly.
+
+```julia
+I = IndexCollection(...)
+v[I]
+```
+
+* This would impact how tensor grid works.
+* To make things homogenous maybe these index collections should be used for the more simple grids too.
+* The function `to_indices` from Base could be useful to implement for `IndexCollection`
--- a/src/Grids/Grids.jl	Sat Dec 23 23:01:28 2023 +0100
+++ b/src/Grids/Grids.jl	Sat Dec 23 23:03:13 2023 +0100
@@ -8,33 +8,29 @@
 export Grid
 export coordinate_size
 export component_type
+export grid_id
+export boundary_id
+export boundary_indices
+export boundary_identifiers
+export boundary_grid
+export eval_on
+export coarsen
+export getcomponent
+export refine
+
+export BoundaryIdentifier
+export TensorGridBoundary
+export CartesianBoundary
 
 export TensorGrid
 export ZeroDimGrid
 
-export TensorGridBoundary
-
-export grid_id
-export boundary_id
-
-export eval_on
-export getcomponent
-
-# BoundaryIdentifier
-export BoundaryIdentifier
+export EquidistantGrid
+export inverse_spacing
+export spacing
+export equidistant_grid
 
 
-# EquidistantGrid
-export EquidistantGrid
-export spacing
-export inverse_spacing
-export boundary_identifiers
-export boundary_grid
-export refine
-export coarsen
-export equidistant_grid
-export CartesianBoundary
-
 abstract type BoundaryIdentifier end
 
 include("grid.jl")
--- a/src/Grids/equidistant_grid.jl	Sat Dec 23 23:01:28 2023 +0100
+++ b/src/Grids/equidistant_grid.jl	Sat Dec 23 23:03:13 2023 +0100
@@ -29,6 +29,7 @@
 Base.IteratorSize(::Type{<:EquidistantGrid}) = Base.HasShape{1}()
 Base.length(g::EquidistantGrid) = length(g.points)
 Base.size(g::EquidistantGrid) = size(g.points)
+Base.size(g::EquidistantGrid, d) = size(g.points)[d]
 
 
 """
@@ -50,7 +51,8 @@
 boundary_identifiers(::EquidistantGrid) = (Lower(), Upper())
 boundary_grid(g::EquidistantGrid, id::Lower) = ZeroDimGrid(g[begin])
 boundary_grid(g::EquidistantGrid, id::Upper) = ZeroDimGrid(g[end])
-
+boundary_indices(g::EquidistantGrid, id::Lower) = (1,)
+boundary_indices(g::EquidistantGrid, id::Upper) = (length(g),)
 
 """
     refine(g::EquidistantGrid, r::Int)
--- a/src/Grids/grid.jl	Sat Dec 23 23:01:28 2023 +0100
+++ b/src/Grids/grid.jl	Sat Dec 23 23:03:13 2023 +0100
@@ -74,6 +74,18 @@
 # TBD: Can we implement a version here that accepts multiple ids and grouped boundaries? Maybe we need multiblock stuff?
 
 """
+    boundary_indices(g::Grid, id::BoundaryIdentifier)
+
+A collection of indices corresponding to the boundary with given id. For grids
+with Cartesian indexing these collections will be tuples with elements of type
+``Union{Int,Colon}``.
+
+When implementing this method it is expected that the returned collection can
+be used to index grid functions to obtain grid functions on the boundary grid.
+"""
+function boundary_indices end
+
+"""
     eval_on(g::Grid, f)
 
 Lazy evaluation of `f` on the grid. `f` can either be on the form `f(x,y,...)`
@@ -100,3 +112,5 @@
 
 _ncomponents(::Type{<:Number}) = 1
 _ncomponents(T::Type{<:SVector}) = length(T)
+
+
--- a/src/Grids/tensor_grid.jl	Sat Dec 23 23:01:28 2023 +0100
+++ b/src/Grids/tensor_grid.jl	Sat Dec 23 23:03:13 2023 +0100
@@ -43,9 +43,9 @@
 _iterate_combine_coords((next,state)) = combine_coordinates(next...), state
 
 Base.IteratorSize(::Type{<:TensorGrid{<:Any, D}}) where D = Base.HasShape{D}()
-Base.eltype(::Type{<:TensorGrid{T}}) where T = T
-Base.length(g::TensorGrid) = sum(length, g.grids)
+Base.length(g::TensorGrid) = prod(length, g.grids)
 Base.size(g::TensorGrid) = LazyTensors.concatenate_tuples(size.(g.grids)...)
+Base.size(g::TensorGrid, d) = size(g)[d]
 
 
 refine(g::TensorGrid, r::Int) = mapreduce(g->refine(g,r), TensorGrid, g.grids)
@@ -73,7 +73,6 @@
     return LazyTensors.concatenate_tuples(per_grid...)
 end
 
-
 """
     boundary_grid(g::TensorGrid, id::TensorGridBoundary)
 
@@ -85,6 +84,16 @@
     return TensorGrid(new_grids...)
 end
 
+function boundary_indices(g::TensorGrid, id::TensorGridBoundary)
+    per_grid_ind = map(g.grids) do g
+        ntuple(i->:, ndims(g))
+    end
+
+    local_b_ind = boundary_indices(g.grids[grid_id(id)], boundary_id(id))
+    b_ind = Base.setindex(per_grid_ind, local_b_ind, grid_id(id))
+
+    return LazyTensors.concatenate_tuples(b_ind...)
+end
 
 function combined_coordinate_vector_type(coordinate_types...)
     combined_coord_length = mapreduce(_ncomponents, +, coordinate_types)
--- a/src/Grids/zero_dim_grid.jl	Sat Dec 23 23:01:28 2023 +0100
+++ b/src/Grids/zero_dim_grid.jl	Sat Dec 23 23:03:13 2023 +0100
@@ -25,3 +25,4 @@
 
 boundary_identifiers(g::ZeroDimGrid) = ()
 boundary_grid(g::ZeroDimGrid, ::Any) = throw(ArgumentError("ZeroDimGrid has no boundaries"))
+boundary_indices(g::ZeroDimGrid, ::Any) = throw(ArgumentError("ZeroDimGrid has no boundaries"))
--- a/test/Grids/equidistant_grid_test.jl	Sat Dec 23 23:01:28 2023 +0100
+++ b/test/Grids/equidistant_grid_test.jl	Sat Dec 23 23:03:13 2023 +0100
@@ -31,6 +31,8 @@
         @test size(EquidistantGrid(0:10)) == (11,)
         @test size(EquidistantGrid(0:0.1:10)) == (101,)
 
+        @test size(EquidistantGrid(0:0.1:10),1) == 101
+
         @test collect(EquidistantGrid(0:0.1:0.5)) == [0.0, 0.1, 0.2, 0.3, 0.4, 0.5]
 
         @test Base.IteratorSize(EquidistantGrid{Float64, StepRange{Float64}}) == Base.HasShape{1}()
@@ -66,6 +68,17 @@
         @test boundary_grid(g, Upper()) == ZeroDimGrid(1.0)
     end
 
+    @testset "boundary_indices" begin
+        g = EquidistantGrid(0:0.1:1)
+        @test boundary_indices(g, Lower()) == (1,)
+        @test boundary_indices(g, Upper()) == (11,)
+
+        g = EquidistantGrid(2:0.1:10)
+        @test boundary_indices(g, Lower()) == (1,)
+        @test boundary_indices(g, Upper()) == (81,)
+
+    end
+
     @testset "refine" begin
         g = EquidistantGrid(0:0.1:1)
         @test refine(g, 1) == g
@@ -107,8 +120,13 @@
     @testset "Base" begin
         @test eltype(equidistant_grid(4,0.0,1.0)) == Float64
         @test eltype(equidistant_grid((4,3),(0,0),(1,3))) <: AbstractVector{Float64}
+
         @test size(equidistant_grid(4,0.0,1.0)) == (4,)
         @test size(equidistant_grid((5,3), (0.0,0.0), (2.0,1.0))) == (5,3)
+
+        @test size(equidistant_grid((5,3), (0.0,0.0), (2.0,1.0)),1) == 5
+        @test size(equidistant_grid((5,3), (0.0,0.0), (2.0,1.0)),2) == 3
+
         @test ndims(equidistant_grid(4,0.0,1.0)) == 1
         @test ndims(equidistant_grid((5,3), (0.0,0.0), (2.0,1.0))) == 2
     end
--- a/test/Grids/tensor_grid_test.jl	Sat Dec 23 23:01:28 2023 +0100
+++ b/test/Grids/tensor_grid_test.jl	Sat Dec 23 23:03:13 2023 +0100
@@ -85,12 +85,30 @@
         @test eltype(TensorGrid(g₁, g₄)) == SVector{3,Float64}
         @test eltype(TensorGrid(g₁, g₄, g₂)) == SVector{4,Float64}
 
+        @test eltype(typeof(TensorGrid(g₁, g₂))) == SVector{2,Float64}
+        @test eltype(typeof(TensorGrid(g₁, g₃))) == SVector{2,Float64}
+        @test eltype(typeof(TensorGrid(g₁, g₂, g₃))) == SVector{3,Float64}
+        @test eltype(typeof(TensorGrid(g₁, g₄))) == SVector{3,Float64}
+        @test eltype(typeof(TensorGrid(g₁, g₄, g₂))) == SVector{4,Float64}
+
         @test size(TensorGrid(g₁, g₂)) == (11,6)
         @test size(TensorGrid(g₁, g₃)) == (11,10)
         @test size(TensorGrid(g₁, g₂, g₃)) == (11,6,10)
         @test size(TensorGrid(g₁, g₄)) == (11,)
         @test size(TensorGrid(g₁, g₄, g₂)) == (11,6)
 
+        @test size(TensorGrid(g₁, g₂, g₃),1) == 11
+        @test size(TensorGrid(g₁, g₂, g₃),2) == 6
+        @test size(TensorGrid(g₁, g₂, g₃),3) == 10
+        @test size(TensorGrid(g₁, g₄, g₂),1) == 11
+        @test size(TensorGrid(g₁, g₄, g₂),2) == 6
+
+        @test length(TensorGrid(g₁, g₂)) == 66
+        @test length(TensorGrid(g₁, g₃)) == 110
+        @test length(TensorGrid(g₁, g₂, g₃)) == 660
+        @test length(TensorGrid(g₁, g₄)) == 11
+        @test length(TensorGrid(g₁, g₄, g₂)) == 66
+
         @test Base.IteratorSize(TensorGrid(g₁, g₂)) == Base.HasShape{2}()
         @test Base.IteratorSize(TensorGrid(g₁, g₃)) == Base.HasShape{2}()
         @test Base.IteratorSize(TensorGrid(g₁, g₂, g₃)) == Base.HasShape{3}()
@@ -151,6 +169,21 @@
         @test boundary_grid(TensorGrid(g₁, g₂), TensorGridBoundary{1, Upper}()) == TensorGrid(ZeroDimGrid(g₁[end]), g₂)
         @test boundary_grid(TensorGrid(g₁, g₄), TensorGridBoundary{1, Upper}()) == TensorGrid(ZeroDimGrid(g₁[end]), g₄)
     end
+
+    @testset "boundary_indices" begin
+        g₁ = EquidistantGrid(range(0,1,length=11))
+        g₂ = EquidistantGrid(range(2,3,length=6))
+        g₄ = ZeroDimGrid(@SVector[1,2])
+
+        @test boundary_indices(TensorGrid(g₁, g₂), TensorGridBoundary{1, Lower}()) == (1,:)
+        @test boundary_indices(TensorGrid(g₁, g₂), TensorGridBoundary{1, Upper}()) == (11,:)
+        @test boundary_indices(TensorGrid(g₁, g₂), TensorGridBoundary{2, Lower}()) == (:,1)
+        @test boundary_indices(TensorGrid(g₁, g₂), TensorGridBoundary{2, Upper}()) == (:,6)
+        @test boundary_indices(TensorGrid(g₁, g₄), TensorGridBoundary{1, Lower}()) == (1,)
+        @test boundary_indices(TensorGrid(g₁, g₄), TensorGridBoundary{1, Upper}()) == (11,)
+        @test boundary_indices(TensorGrid(g₄,g₁), TensorGridBoundary{2, Lower}()) == (1,)
+        @test boundary_indices(TensorGrid(g₄,g₁), TensorGridBoundary{2, Upper}()) == (11,)
+    end
 end
 
 @testset "combined_coordinate_vector_type" begin
--- a/test/Grids/zero_dim_grid_test.jl	Sat Dec 23 23:01:28 2023 +0100
+++ b/test/Grids/zero_dim_grid_test.jl	Sat Dec 23 23:03:13 2023 +0100
@@ -41,4 +41,8 @@
     @testset "boundary_grid" begin
         @test_throws ArgumentError("ZeroDimGrid has no boundaries") boundary_grid(ZeroDimGrid(@SVector[1.0,2.0]), :bid)
     end
+
+    @testset "boundary_indices" begin
+        @test_throws ArgumentError("ZeroDimGrid has no boundaries") boundary_indices(ZeroDimGrid(@SVector[1.0,2.0]), :bid)
+    end
 end