Mercurial > repos > public > sbplib_julia
changeset 1506:535f32316637 feature/grids/curvilinear
Rename from curvilinear to mapped
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Fri, 16 Feb 2024 15:28:19 +0100 |
| parents | 63101a1cd0e6 |
| children | 4df668d00d03 |
| files | src/Grids/Grids.jl src/Grids/curvilinear_grid.jl src/Grids/mapped_grid.jl test/Grids/curvilinear_grid_test.jl test/Grids/mapped_grid_test.jl |
| diffstat | 5 files changed, 256 insertions(+), 257 deletions(-) [+] |
line wrap: on
line diff
--- a/src/Grids/Grids.jl Fri Feb 16 14:33:13 2024 +0100 +++ b/src/Grids/Grids.jl Fri Feb 16 15:28:19 2024 +0100 @@ -33,11 +33,11 @@ export equidistant_grid -# CurvilinearGrid -export CurvilinearGrid +# MappedGrid +export MappedGrid export jacobian export logicalgrid -export curvilinear_grid +export mapped_grid abstract type BoundaryIdentifier end @@ -45,6 +45,6 @@ include("tensor_grid.jl") include("equidistant_grid.jl") include("zero_dim_grid.jl") -include("curvilinear_grid.jl") +include("mapped_grid.jl") end # module
--- a/src/Grids/curvilinear_grid.jl Fri Feb 16 14:33:13 2024 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,60 +0,0 @@ -# TBD: Rename to MappedGrid? -struct CurvilinearGrid{T,D, GT<:Grid{<:Any,D}, CT<:AbstractArray{T,D}, JT<:AbstractArray{<:AbstractArray{<:Any, 2}, D}} <: Grid{T,D} - logicalgrid::GT - physicalcoordinates::CT - jacobian::JT -end - -jacobian(g::CurvilinearGrid) = g.jacobian -logicalgrid(g::CurvilinearGrid) = g.logicalgrid - - -# Indexing interface -Base.getindex(g::CurvilinearGrid, I::Vararg{Int}) = g.physicalcoordinates[I...] -Base.eachindex(g::CurvilinearGrid) = eachindex(g.logicalgrid) - -Base.firstindex(g::CurvilinearGrid, d) = firstindex(g.logicalgrid, d) -Base.lastindex(g::CurvilinearGrid, d) = lastindex(g.logicalgrid, d) - -# Iteration interface - -Base.iterate(g::CurvilinearGrid) = iterate(g.physicalcoordinates) -Base.iterate(g::CurvilinearGrid, state) = iterate(g.physicalcoordinates, state) - -Base.IteratorSize(::Type{<:CurvilinearGrid{<:Any, D}}) where D = Base.HasShape{D}() -Base.length(g::CurvilinearGrid) = length(g.logicalgrid) -Base.size(g::CurvilinearGrid) = size(g.logicalgrid) -Base.size(g::CurvilinearGrid, d) = size(g.logicalgrid, d) - -boundary_identifiers(g::CurvilinearGrid) = boundary_identifiers(g.logicalgrid) -boundary_indices(g::CurvilinearGrid, id::TensorGridBoundary) = boundary_indices(g.logicalgrid, id) - -function boundary_grid(g::CurvilinearGrid, id::TensorGridBoundary) - b_indices = boundary_indices(g.logicalgrid, id) - - # Calculate indices of needed jacobian components - D = ndims(g) - all_indices = SVector{D}(1:D) - free_variable_indices = deleteat(all_indices, grid_id(id)) - jacobian_components = (:, free_variable_indices) - - # Create grid function for boundary grid jacobian - boundary_jacobian = componentview((@view g.jacobian[b_indices...]) , jacobian_components...) - boundary_physicalcoordinates = @view g.physicalcoordinates[b_indices...] - - return CurvilinearGrid( - boundary_grid(g.logicalgrid, id), - boundary_physicalcoordinates, - boundary_jacobian, - ) -end - -function curvilinear_grid(x, J, size...) - D = length(size) - lg = equidistant_grid(size, ntuple(i->0., D), ntuple(i->1., D)) - return CurvilinearGrid( - lg, - map(x,lg), - map(J,lg), - ) -end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/Grids/mapped_grid.jl Fri Feb 16 15:28:19 2024 +0100 @@ -0,0 +1,59 @@ +struct MappedGrid{T,D, GT<:Grid{<:Any,D}, CT<:AbstractArray{T,D}, JT<:AbstractArray{<:AbstractArray{<:Any, 2}, D}} <: Grid{T,D} + logicalgrid::GT + physicalcoordinates::CT + jacobian::JT +end + +jacobian(g::MappedGrid) = g.jacobian +logicalgrid(g::MappedGrid) = g.logicalgrid + + +# Indexing interface +Base.getindex(g::MappedGrid, I::Vararg{Int}) = g.physicalcoordinates[I...] +Base.eachindex(g::MappedGrid) = eachindex(g.logicalgrid) + +Base.firstindex(g::MappedGrid, d) = firstindex(g.logicalgrid, d) +Base.lastindex(g::MappedGrid, d) = lastindex(g.logicalgrid, d) + +# Iteration interface + +Base.iterate(g::MappedGrid) = iterate(g.physicalcoordinates) +Base.iterate(g::MappedGrid, state) = iterate(g.physicalcoordinates, state) + +Base.IteratorSize(::Type{<:MappedGrid{<:Any, D}}) where D = Base.HasShape{D}() +Base.length(g::MappedGrid) = length(g.logicalgrid) +Base.size(g::MappedGrid) = size(g.logicalgrid) +Base.size(g::MappedGrid, d) = size(g.logicalgrid, d) + +boundary_identifiers(g::MappedGrid) = boundary_identifiers(g.logicalgrid) +boundary_indices(g::MappedGrid, id::TensorGridBoundary) = boundary_indices(g.logicalgrid, id) + +function boundary_grid(g::MappedGrid, id::TensorGridBoundary) + b_indices = boundary_indices(g.logicalgrid, id) + + # Calculate indices of needed jacobian components + D = ndims(g) + all_indices = SVector{D}(1:D) + free_variable_indices = deleteat(all_indices, grid_id(id)) + jacobian_components = (:, free_variable_indices) + + # Create grid function for boundary grid jacobian + boundary_jacobian = componentview((@view g.jacobian[b_indices...]) , jacobian_components...) + boundary_physicalcoordinates = @view g.physicalcoordinates[b_indices...] + + return MappedGrid( + boundary_grid(g.logicalgrid, id), + boundary_physicalcoordinates, + boundary_jacobian, + ) +end + +function mapped_grid(x, J, size...) + D = length(size) + lg = equidistant_grid(size, ntuple(i->0., D), ntuple(i->1., D)) + return MappedGrid( + lg, + map(x,lg), + map(J,lg), + ) +end
--- a/test/Grids/curvilinear_grid_test.jl Fri Feb 16 14:33:13 2024 +0100 +++ /dev/null Thu Jan 01 00:00:00 1970 +0000 @@ -1,193 +0,0 @@ -using Sbplib.Grids -using Sbplib.RegionIndices -using Test -using StaticArrays - -@testset "CurvilinearGrid" begin - lg = equidistant_grid((11,11), (0,0), (1,1)) # TODO: Change dims of the grid to be different - x̄ = map(ξ̄ -> 2ξ̄, lg) - J = map(ξ̄ -> @SArray(fill(2., 2, 2)), lg) - cg = CurvilinearGrid(lg, x̄, J) - - # TODO: Test constructor for different dims of range and domain for the coordinates - # TODO: Test constructor with different type than TensorGrid. a dummy type? - - @test_broken false # @test_throws ArgumentError("Sizes must match") CurvilinearGrid(lg, map(ξ̄ -> @SArray[ξ̄[1], ξ̄[2], -ξ̄[1]], lg), rand(SMatrix{2,3,Float64},15,11)) - - - @test cg isa Grid{SVector{2, Float64},2} - - @test jacobian(cg) isa Array{<:AbstractMatrix} - @test logicalgrid(cg) isa Grid - - @testset "Indexing Interface" begin - cg = CurvilinearGrid(lg, x̄, J) - @test cg[1,1] == [0.0, 0.0] - @test cg[4,2] == [0.6, 0.2] - @test cg[6,10] == [1., 1.8] - - @test cg[begin, begin] == [0.0, 0.0] - @test cg[end,end] == [2.0, 2.0] - @test cg[begin,end] == [0., 2.] - - @test eachindex(cg) == CartesianIndices((11,11)) - - @testset "cartesian indexing" begin - cases = [ - (1,1) , - (3,5) , - (10,6), - (1,1) , - (3,2) , - ] - - @testset "i = $is" for (lg, is) ∈ cases - @test cg[CartesianIndex(is...)] == cg[is...] - end - end - - @testset "eachindex" begin - @test eachindex(cg) == CartesianIndices((11,11)) - end - - @testset "firstindex" begin - @test firstindex(cg, 1) == 1 - @test firstindex(cg, 2) == 1 - end - - @testset "lastindex" begin - @test lastindex(cg, 1) == 11 - @test lastindex(cg, 2) == 11 - end - end - # TODO: Test with different types of logical grids - - @testset "Iterator interface" begin - sg = CurvilinearGrid( - equidistant_grid((15,11), (0,0), (1,1)), - map(ξ̄ -> @SArray[ξ̄[1], ξ̄[2], -ξ̄[1]], lg), rand(SMatrix{2,3,Float64},15,11) - ) - - @test eltype(cg) == SVector{2,Float64} - @test eltype(sg) == SVector{3,Float64} - - @test eltype(typeof(cg)) == SVector{2,Float64} - @test eltype(typeof(sg)) == SVector{3,Float64} - - @test size(cg) == (11,11) - @test size(sg) == (15,11) - - @test size(cg,2) == 11 - @test size(sg,2) == 11 - - @test length(cg) == 121 - @test length(sg) == 165 - - @test Base.IteratorSize(cg) == Base.HasShape{2}() - @test Base.IteratorSize(typeof(cg)) == Base.HasShape{2}() - - @test Base.IteratorSize(sg) == Base.HasShape{2}() - @test Base.IteratorSize(typeof(sg)) == Base.HasShape{2}() - - element, state = iterate(cg) - @test element == lg[1,1].*2 - element, _ = iterate(cg, state) - @test element == lg[2,1].*2 - - element, state = iterate(sg) - @test element == sg.physicalcoordinates[1,1] - element, _ = iterate(sg, state) - @test element == sg.physicalcoordinates[2,1] - - @test collect(cg) == 2 .* lg - end - - @testset "Base" begin - @test ndims(cg) == 2 - end - - @testset "boundary_identifiers" begin - @test boundary_identifiers(cg) == boundary_identifiers(lg) - end - - @testset "boundary_indices" begin - @test boundary_indices(cg, CartesianBoundary{1,Lower}()) == boundary_indices(lg,CartesianBoundary{1,Lower}()) - @test boundary_indices(cg, CartesianBoundary{2,Lower}()) == boundary_indices(lg,CartesianBoundary{2,Lower}()) - @test boundary_indices(cg, CartesianBoundary{1,Upper}()) == boundary_indices(lg,CartesianBoundary{1,Upper}()) - end - - @testset "boundary_grid" begin - x̄((ξ, η)) = @SVector[ξ, η*(1+ξ*(ξ-1))] - J((ξ, η)) = @SMatrix[ - 1 0; - η*(2ξ-1) 1+ξ*(ξ-1); - ] - - cg = curvilinear_grid(x̄, J, 10, 11) - J1((ξ, η)) = @SMatrix[ - 1 ; - η*(2ξ-1); - ] - J2((ξ, η)) = @SMatrix[ - 0; - 1+ξ*(ξ-1); - ] - - function test_boundary_grid(cg, bId, Jb) - bg = boundary_grid(cg, bId) - - lg = logicalgrid(cg) - expected_bg = CurvilinearGrid( - boundary_grid(lg, bId), - map(x̄, boundary_grid(lg, bId)), - map(Jb, boundary_grid(lg, bId)), - ) - - @testset let bId=bId, bg=bg, expected_bg=expected_bg - @test collect(bg) == collect(expected_bg) - @test logicalgrid(bg) == logicalgrid(expected_bg) - @test jacobian(bg) == jacobian(expected_bg) - # TODO: Implement equality of a curvilinear grid and simlify the above - end - end - - @testset test_boundary_grid(cg, TensorGridBoundary{1, Lower}(), J2) - @testset test_boundary_grid(cg, TensorGridBoundary{1, Upper}(), J2) - @testset test_boundary_grid(cg, TensorGridBoundary{2, Lower}(), J1) - @testset test_boundary_grid(cg, TensorGridBoundary{2, Upper}(), J1) - end - - # TBD: Should curvilinear grid support refining and coarsening? - # This would require keeping the coordinate mapping around which seems burdensome, and might increase compilation time? - @testset "refine" begin - @test_broken refine(cg, 1) == cg - @test_broken refine(cg, 2) == CurvilinearGrid(refine(lg,2), x̄, J) - @test_broken refine(cg, 3) == CurvilinearGrid(refine(lg,3), x̄, J) - end - - @testset "coarsen" begin - lg = equidistant_grid((11,11), (0,0), (1,1)) # TODO: Change dims of the grid to be different - x̄ = map(ξ̄ -> 2ξ̄, lg) - J = map(ξ̄ -> @SArray(fill(2., 2, 2)), lg) - cg = CurvilinearGrid(lg, x̄, J) - - @test_broken coarsen(cg, 1) == cg - @test_broken coarsen(cg, 2) == CurvilinearGrid(coarsen(lg,2), x̄, J) - - @test_broken false # @test_throws DomainError(3, "Size minus 1 must be divisible by the ratio.") coarsen(cg, 3) - end -end - -@testset "curvilinear_grid" begin - x̄((ξ, η)) = @SVector[ξ, η*(1+ξ*(ξ-1))] - J((ξ, η)) = @SMatrix[ - 1 0; - η*(2ξ-1) 1+ξ*(ξ-1); - ] - cg = curvilinear_grid(x̄, J, 10, 11) - @test cg isa CurvilinearGrid{SVector{2,Float64}, 2} - - lg = equidistant_grid((10,11), (0,0), (1,1)) - @test logicalgrid(cg) == lg - @test collect(cg) == map(x̄, lg) -end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test/Grids/mapped_grid_test.jl Fri Feb 16 15:28:19 2024 +0100 @@ -0,0 +1,193 @@ +using Sbplib.Grids +using Sbplib.RegionIndices +using Test +using StaticArrays + +@testset "MappedGrid" begin + lg = equidistant_grid((11,11), (0,0), (1,1)) # TODO: Change dims of the grid to be different + x̄ = map(ξ̄ -> 2ξ̄, lg) + J = map(ξ̄ -> @SArray(fill(2., 2, 2)), lg) + mg = MappedGrid(lg, x̄, J) + + # TODO: Test constructor for different dims of range and domain for the coordinates + # TODO: Test constructor with different type than TensorGrid. a dummy type? + + @test_broken false # @test_throws ArgumentError("Sizes must match") MappedGrid(lg, map(ξ̄ -> @SArray[ξ̄[1], ξ̄[2], -ξ̄[1]], lg), rand(SMatrix{2,3,Float64},15,11)) + + + @test mg isa Grid{SVector{2, Float64},2} + + @test jacobian(mg) isa Array{<:AbstractMatrix} + @test logicalgrid(mg) isa Grid + + @testset "Indexing Interface" begin + mg = MappedGrid(lg, x̄, J) + @test mg[1,1] == [0.0, 0.0] + @test mg[4,2] == [0.6, 0.2] + @test mg[6,10] == [1., 1.8] + + @test mg[begin, begin] == [0.0, 0.0] + @test mg[end,end] == [2.0, 2.0] + @test mg[begin,end] == [0., 2.] + + @test eachindex(mg) == CartesianIndices((11,11)) + + @testset "cartesian indexing" begin + cases = [ + (1,1) , + (3,5) , + (10,6), + (1,1) , + (3,2) , + ] + + @testset "i = $is" for (lg, is) ∈ cases + @test mg[CartesianIndex(is...)] == mg[is...] + end + end + + @testset "eachindex" begin + @test eachindex(mg) == CartesianIndices((11,11)) + end + + @testset "firstindex" begin + @test firstindex(mg, 1) == 1 + @test firstindex(mg, 2) == 1 + end + + @testset "lastindex" begin + @test lastindex(mg, 1) == 11 + @test lastindex(mg, 2) == 11 + end + end + # TODO: Test with different types of logical grids + + @testset "Iterator interface" begin + sg = MappedGrid( + equidistant_grid((15,11), (0,0), (1,1)), + map(ξ̄ -> @SArray[ξ̄[1], ξ̄[2], -ξ̄[1]], lg), rand(SMatrix{2,3,Float64},15,11) + ) + + @test eltype(mg) == SVector{2,Float64} + @test eltype(sg) == SVector{3,Float64} + + @test eltype(typeof(mg)) == SVector{2,Float64} + @test eltype(typeof(sg)) == SVector{3,Float64} + + @test size(mg) == (11,11) + @test size(sg) == (15,11) + + @test size(mg,2) == 11 + @test size(sg,2) == 11 + + @test length(mg) == 121 + @test length(sg) == 165 + + @test Base.IteratorSize(mg) == Base.HasShape{2}() + @test Base.IteratorSize(typeof(mg)) == Base.HasShape{2}() + + @test Base.IteratorSize(sg) == Base.HasShape{2}() + @test Base.IteratorSize(typeof(sg)) == Base.HasShape{2}() + + element, state = iterate(mg) + @test element == lg[1,1].*2 + element, _ = iterate(mg, state) + @test element == lg[2,1].*2 + + element, state = iterate(sg) + @test element == sg.physicalcoordinates[1,1] + element, _ = iterate(sg, state) + @test element == sg.physicalcoordinates[2,1] + + @test collect(mg) == 2 .* lg + end + + @testset "Base" begin + @test ndims(mg) == 2 + end + + @testset "boundary_identifiers" begin + @test boundary_identifiers(mg) == boundary_identifiers(lg) + end + + @testset "boundary_indices" begin + @test boundary_indices(mg, CartesianBoundary{1,Lower}()) == boundary_indices(lg,CartesianBoundary{1,Lower}()) + @test boundary_indices(mg, CartesianBoundary{2,Lower}()) == boundary_indices(lg,CartesianBoundary{2,Lower}()) + @test boundary_indices(mg, CartesianBoundary{1,Upper}()) == boundary_indices(lg,CartesianBoundary{1,Upper}()) + end + + @testset "boundary_grid" begin + x̄((ξ, η)) = @SVector[ξ, η*(1+ξ*(ξ-1))] + J((ξ, η)) = @SMatrix[ + 1 0; + η*(2ξ-1) 1+ξ*(ξ-1); + ] + + mg = mapped_grid(x̄, J, 10, 11) + J1((ξ, η)) = @SMatrix[ + 1 ; + η*(2ξ-1); + ] + J2((ξ, η)) = @SMatrix[ + 0; + 1+ξ*(ξ-1); + ] + + function test_boundary_grid(mg, bId, Jb) + bg = boundary_grid(mg, bId) + + lg = logicalgrid(mg) + expected_bg = MappedGrid( + boundary_grid(lg, bId), + map(x̄, boundary_grid(lg, bId)), + map(Jb, boundary_grid(lg, bId)), + ) + + @testset let bId=bId, bg=bg, expected_bg=expected_bg + @test collect(bg) == collect(expected_bg) + @test logicalgrid(bg) == logicalgrid(expected_bg) + @test jacobian(bg) == jacobian(expected_bg) + # TODO: Implement equality of a curvilinear grid and simlify the above + end + end + + @testset test_boundary_grid(mg, TensorGridBoundary{1, Lower}(), J2) + @testset test_boundary_grid(mg, TensorGridBoundary{1, Upper}(), J2) + @testset test_boundary_grid(mg, TensorGridBoundary{2, Lower}(), J1) + @testset test_boundary_grid(mg, TensorGridBoundary{2, Upper}(), J1) + end + + # TBD: Should curvilinear grid support refining and coarsening? + # This would require keeping the coordinate mapping around which seems burdensome, and might increase compilation time? + @testset "refine" begin + @test_broken refine(mg, 1) == mg + @test_broken refine(mg, 2) == MappedGrid(refine(lg,2), x̄, J) + @test_broken refine(mg, 3) == MappedGrid(refine(lg,3), x̄, J) + end + + @testset "coarsen" begin + lg = equidistant_grid((11,11), (0,0), (1,1)) # TODO: Change dims of the grid to be different + x̄ = map(ξ̄ -> 2ξ̄, lg) + J = map(ξ̄ -> @SArray(fill(2., 2, 2)), lg) + mg = MappedGrid(lg, x̄, J) + + @test_broken coarsen(mg, 1) == mg + @test_broken coarsen(mg, 2) == MappedGrid(coarsen(lg,2), x̄, J) + + @test_broken false # @test_throws DomainError(3, "Size minus 1 must be divisible by the ratio.") coarsen(mg, 3) + end +end + +@testset "mapped_grid" begin + x̄((ξ, η)) = @SVector[ξ, η*(1+ξ*(ξ-1))] + J((ξ, η)) = @SMatrix[ + 1 0; + η*(2ξ-1) 1+ξ*(ξ-1); + ] + mg = mapped_grid(x̄, J, 10, 11) + @test mg isa MappedGrid{SVector{2,Float64}, 2} + + lg = equidistant_grid((10,11), (0,0), (1,1)) + @test logicalgrid(mg) == lg + @test collect(mg) == map(x̄, lg) +end