Mercurial > repos > public > sbplib_julia
changeset 1951:6c1bb9bdb092 feature/grids/geometry_functions
Merge feature/grids/manifolds
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Fri, 07 Feb 2025 22:38:39 +0100 |
| parents | 6859089b361e (current diff) dd77b45ee1ac (diff) |
| children | b0915f43b122 f46b4f1fa82c |
| files | src/Grids/Grids.jl |
| diffstat | 10 files changed, 369 insertions(+), 37 deletions(-) [+] |
line wrap: on
line diff
--- a/docs/make.jl Mon Feb 03 15:44:07 2025 +0100 +++ b/docs/make.jl Fri Feb 07 22:38:39 2025 +0100 @@ -30,6 +30,7 @@ "operator_file_format.md", "grids_and_grid_functions.md", "matrix_and_tensor_representations.md", + "manifolds_charts_atlases.md" "Submodules" => [ "submodules/grids.md", "submodules/lazy_tensors.md",
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/docs/src/manifolds_charts_atlases.md Fri Feb 07 22:38:39 2025 +0100 @@ -0,0 +1,36 @@ +# Manifolds, Charts, and Atlases + +To construct grids on more complicated geometries we use manifolds described +by one or more charts. The charts describe a mapping from some parameter space +to the geometry that we are interested in. If there are more than one chart +for a given geometry this collection of charts and their connection is +described by and atlas. + +For the construction of differential and difference operators on a manifold +with a chart the library needs to know the Jacobian of the mapping as a +function of coordinates in the logical parameter space. Internally, +Diffinitive.jl uses a local Jacobian function, `Grids.jacobian(f, ξ)`. For +geometry objects provided by the library this function should have fast and +efficient implementations. If you are creating your own mapping functions you +can implement `Grids.jacobian` for your function or type, for example + +```julia +f(x) = 2x +Grids.jacobian(::typeof(f), x) = fill(2, length(x)) +``` + +```julia +struct F end +(::F)(x) = 2x +Grids.jacobian(::F, x) = fill(2,length(x)) +``` + +You can also provide a fallback function using one of the many automatic +differentiation packages, for example + +```julia +using ForwardDiff +Grids.jacobian(f,x) = ForwardDiff.jacobian(f,x) +``` + +<!-- What more needs to be said here? --/>
--- a/src/Grids/Grids.jl Mon Feb 03 15:44:07 2025 +0100 +++ b/src/Grids/Grids.jl Fri Feb 07 22:38:39 2025 +0100 @@ -30,6 +30,7 @@ export charts export connections export CartesianAtlas +export UnstructuredAtlas export parameterspace @@ -64,6 +65,8 @@ export spacing export equidistant_grid +export MultiBlockBoundary + # MappedGrid export MappedGrid @@ -73,8 +76,9 @@ export metric_tensor include("parameter_space.jl") +include("grid.jl") +include("multiblockgrids.jl") include("manifolds.jl") -include("grid.jl") include("tensor_grid.jl") include("equidistant_grid.jl") include("zero_dim_grid.jl")
--- a/src/Grids/grid.jl Mon Feb 03 15:44:07 2025 +0100 +++ b/src/Grids/grid.jl Fri Feb 07 22:38:39 2025 +0100 @@ -122,8 +122,6 @@ """ function boundary_identifiers end -# TBD: Boundary identifiers for charts and atlases? - """ boundary_grid(g::Grid, id::BoundaryIdentifier)
--- a/src/Grids/manifolds.jl Mon Feb 03 15:44:07 2025 +0100 +++ b/src/Grids/manifolds.jl Fri Feb 07 22:38:39 2025 +0100 @@ -33,8 +33,8 @@ jacobian(c::Chart, ξ) = jacobian(c.mapping, ξ) # TBD: Can we register a error hint for when jacobian is called with a function that doesn't have a registered jacobian? +boundary_identifiers(c::Chart) = boundary_identifiers(parameterspace(c)) -# TBD: Should Charts, parameterspaces, Atlases, have boundary names? """ Atlas @@ -54,21 +54,99 @@ """ connections(::Atlas) -TBD: What exactly should this return? +Collection of pairs of multiblock boundary identifiers. """ function connections end -struct CartesianAtlas <: Atlas - charts::Matrix{Chart} + +""" + CartesianAtlas{D,C<:Chart,AT<:AbstractArray{C,D}} <: Atlas + +An atlas where the charts are arranged and connected like an array. +""" +struct CartesianAtlas{D,C<:Chart,AT<:AbstractArray{C,D}} <: Atlas + charts::AT end charts(a::CartesianAtlas) = a.charts -connections(a::CartesianAtlas) = nothing + +function connections(a::CartesianAtlas) + c = Tuple{MultiBlockBoundary, MultiBlockBoundary}[] + + for d ∈ 1:ndims(charts(a)) + Is = eachslice(CartesianIndices(charts(a)); dims=d) + for i ∈ 1:length(Is)-1 # For each interface between slices + for jk ∈ eachindex(Is[i]) # For each block in slice + Iᵢⱼₖ = Tuple(Is[i][jk]) + Iᵢ₊₁ⱼₖ = Tuple(Is[i+1][jk]) + push!(c, + ( + MultiBlockBoundary{Iᵢⱼₖ, CartesianBoundary{d,UpperBoundary}}(), + MultiBlockBoundary{Iᵢ₊₁ⱼₖ, CartesianBoundary{d,LowerBoundary}}(), + ) + ) + end + end + end + + return c +end + +""" + boundary_identifiers(a::CartesianAtlas) -struct UnstructuredAtlas <: Atlas - charts::Vector{Chart} - connections +All non-connected boundaries of the charts of `a`. +""" +function boundary_identifiers(a::CartesianAtlas) + bs = MultiBlockBoundary[] + + for d ∈ 1:ndims(charts(a)) + Is = eachslice(CartesianIndices(charts(a)); dims=d) + + for (i,b) ∈ ((1,LowerBoundary),(length(Is),UpperBoundary)) # For first and last slice + for jk ∈ eachindex(Is[i]) # For each block in slice + Iᵢⱼₖ = Tuple(Is[i][jk]) + push!(bs, + MultiBlockBoundary{Iᵢⱼₖ, CartesianBoundary{d,b}}(), + ) + end + end + end + + return bs +end + + +""" + UnstructuredAtlas{C<:Chart, CN<:Tuple{MultiBlockBoundary,MultiBlockBoundary}, ...} <: Atlas + +An atlas with connections determined by a vector `MultiBlockBoundary` pairs. +""" +struct UnstructuredAtlas{C<:Chart, CN<:Tuple{MultiBlockBoundary,MultiBlockBoundary}, CV<:AbstractVector{C}, CNV<:AbstractVector{CN}} <: Atlas + charts::CV + connections::CNV end charts(a::UnstructuredAtlas) = a.charts -connections(a::UnstructuredAtlas) = nothing +connections(a::UnstructuredAtlas) = a.connections + +""" + boundary_identifiers(a::UnstructuredAtlas) + +All non-connected boundaries of the charts of `a`. +""" +function boundary_identifiers(a::UnstructuredAtlas) + bs = MultiBlockBoundary[] + + for (i,c) ∈ enumerate(charts(a)) + for b ∈ boundary_identifiers(c) + mbb = MultiBlockBoundary{i,typeof(b)}() + + if !any(cn->mbb∈cn, connections(a)) + push!(bs, mbb) + end + end + end + + return bs +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/Grids/multiblockgrids.jl Fri Feb 07 22:38:39 2025 +0100 @@ -0,0 +1,9 @@ +""" + MultiBlockBoundary{N, BID} <: BoundaryIdentifier + +A boundary identifier for a multiblock grids. `N` Specifies which grid and +`BID` which boundary on that grid. +""" +struct MultiBlockBoundary{N, BID} <: BoundaryIdentifier end +grid_id(::MultiBlockBoundary{N, BID}) where {N, BID} = N +boundary_id(::MultiBlockBoundary{N, BID}) where {N, BID} = BID()
--- a/src/Grids/parameter_space.jl Mon Feb 03 15:44:07 2025 +0100 +++ b/src/Grids/parameter_space.jl Fri Feb 07 22:38:39 2025 +0100 @@ -45,6 +45,8 @@ """ limits(i::Interval) = (i.a, i.b) +boundary_identifiers(::Interval) = (LowerBoundary(), UpperBoundary()) + """ unitinterval(T=Float64) @@ -91,6 +93,16 @@ """ limits(box::HyperBox) = (box.a, box.b) +function boundary_identifiers(box::HyperBox) + mapreduce(vcat, 1:ndims(box)) do d + [ + CartesianBoundary{d, LowerBoundary}(), + CartesianBoundary{d, UpperBoundary}(), + ] + end +end + + """ unitsquare(T=Float64)
--- a/test/Grids/manifolds_test.jl Mon Feb 03 15:44:07 2025 +0100 +++ b/test/Grids/manifolds_test.jl Fri Feb 07 22:38:39 2025 +0100 @@ -4,42 +4,205 @@ using Diffinitive.RegionIndices using Diffinitive.LazyTensors +using StaticArrays + +west = CartesianBoundary{1,LowerBoundary} +east = CartesianBoundary{1,UpperBoundary} +south = CartesianBoundary{2,LowerBoundary} +north = CartesianBoundary{2,UpperBoundary} +bottom = CartesianBoundary{3, LowerBoundary} +top = CartesianBoundary{3, UpperBoundary} + @testset "Chart" begin - c = Chart(x->2x, unitsquare()) + X(ξ) = 2ξ + Grids.jacobian(::typeof(X), ξ) = @SVector[2,2] + c = Chart(X, unitsquare()) @test c isa Chart{2} @test c([3,2]) == [6,4] @test parameterspace(c) == unitsquare() @test ndims(c) == 2 - @test_broken jacobian(c, [3,2]) + @test jacobian(c, [3,2]) == [2,2] + + @test Set(boundary_identifiers(Chart(X,unitsquare()))) == Set([east(),west(),south(),north()]) end @testset "CartesianAtlas" begin - c = Chart(identity, unitsquare()) + @testset "Constructors" begin + c = Chart(identity, unitsquare()) + @test CartesianAtlas([c c; c c]) isa Atlas + + c2 = Chart(x->2x, unitsquare()) + @test CartesianAtlas([c c2; c2 c]) isa CartesianAtlas + @test CartesianAtlas(@SMatrix[c c; c c]) isa CartesianAtlas + @test CartesianAtlas(@SMatrix[c c2; c2 c]) isa CartesianAtlas + end + + @testset "Getters" begin + c = Chart(identity, unitsquare()) + a = CartesianAtlas([c c; c c]) + @test charts(a) == [c c; c c] + end + + @testset "connections" begin + # 2D + a = CartesianAtlas(fill(Chart(identity, unitsquare()), 2,3)) + + @test Set(connections(a)) == Set([ + (MultiBlockBoundary{(1,1), east}(), MultiBlockBoundary{(2,1), west}()), + (MultiBlockBoundary{(1,1), north}(), MultiBlockBoundary{(1,2), south}()), + (MultiBlockBoundary{(2,1), north}(), MultiBlockBoundary{(2,2), south}()), + (MultiBlockBoundary{(1,2), east}(), MultiBlockBoundary{(2,2), west}()), + (MultiBlockBoundary{(1,2), north}(), MultiBlockBoundary{(1,3), south}()), + (MultiBlockBoundary{(2,2), north}(), MultiBlockBoundary{(2,3), south}()), + (MultiBlockBoundary{(1,3), east}(), MultiBlockBoundary{(2,3), west}()), + ]) + + # 3D + a = CartesianAtlas(fill(Chart(identity, unitcube()), 2,2,3)) + @test Set(connections(a)) == Set([ + (MultiBlockBoundary{(1,1,1), east}(), MultiBlockBoundary{(2,1,1), west}()), + (MultiBlockBoundary{(1,1,1), north}(), MultiBlockBoundary{(1,2,1), south}()), + (MultiBlockBoundary{(2,1,1), north}(), MultiBlockBoundary{(2,2,1), south}()), + (MultiBlockBoundary{(1,2,1), east}(), MultiBlockBoundary{(2,2,1), west}()), + + (MultiBlockBoundary{(1,1,2), east}(), MultiBlockBoundary{(2,1,2), west}()), + (MultiBlockBoundary{(1,1,2), north}(), MultiBlockBoundary{(1,2,2), south}()), + (MultiBlockBoundary{(2,1,2), north}(), MultiBlockBoundary{(2,2,2), south}()), + (MultiBlockBoundary{(1,2,2), east}(), MultiBlockBoundary{(2,2,2), west}()), + + (MultiBlockBoundary{(1,1,3), east}(), MultiBlockBoundary{(2,1,3), west}()), + (MultiBlockBoundary{(1,1,3), north}(), MultiBlockBoundary{(1,2,3), south}()), + (MultiBlockBoundary{(2,1,3), north}(), MultiBlockBoundary{(2,2,3), south}()), + (MultiBlockBoundary{(1,2,3), east}(), MultiBlockBoundary{(2,2,3), west}()), + + (MultiBlockBoundary{(1,1,1), top}(), MultiBlockBoundary{(1,1,2), bottom}()), + (MultiBlockBoundary{(2,1,1), top}(), MultiBlockBoundary{(2,1,2), bottom}()), + (MultiBlockBoundary{(1,2,1), top}(), MultiBlockBoundary{(1,2,2), bottom}()), + (MultiBlockBoundary{(2,2,1), top}(), MultiBlockBoundary{(2,2,2), bottom}()), - a = CartesianAtlas([c c; c c]) - @test a isa Atlas - @test charts(a) == [c c; c c] - @test_broken connections(a) == [ - ( - ((1,1), CartesianBoundary{1,UpperBoundary}()), - ((1,2), CartesianBoundary{1,LowerBoundary}()), - ), - ( - ((1,1), CartesianBoundary{2,LowerBoundary}()), - ((2,1), CartesianBoundary{2,UpperBoundary}()), - ), - ( - ((2,2), CartesianBoundary{1,LowerBoundary}()), - ((2,1), CartesianBoundary{1,UpperBoundary}()), - ), - ( - ((2,2), CartesianBoundary{2,UpperBoundary}()), - ((1,2), CartesianBoundary{2,LowerBoundary}()), - ) - ] + (MultiBlockBoundary{(1,1,2), top}(), MultiBlockBoundary{(1,1,3), bottom}()), + (MultiBlockBoundary{(2,1,2), top}(), MultiBlockBoundary{(2,1,3), bottom}()), + (MultiBlockBoundary{(1,2,2), top}(), MultiBlockBoundary{(1,2,3), bottom}()), + (MultiBlockBoundary{(2,2,2), top}(), MultiBlockBoundary{(2,2,3), bottom}()), + ]) + end + + @testset "boundary_identifiers" begin + # 2D + a = CartesianAtlas(fill(Chart(identity, unitcube()), 2,3)) + @test Set(boundary_identifiers(a)) == Set([ + MultiBlockBoundary{(1,1), south}(), + MultiBlockBoundary{(2,1), south}(), + MultiBlockBoundary{(2,1), east}(), + MultiBlockBoundary{(2,2), east}(), + MultiBlockBoundary{(2,3), east}(), + MultiBlockBoundary{(1,3), north}(), + MultiBlockBoundary{(2,3), north}(), + MultiBlockBoundary{(1,1), west}(), + MultiBlockBoundary{(1,2), west}(), + MultiBlockBoundary{(1,3), west}(), + ]) + + # 3D + a = CartesianAtlas(fill(Chart(identity, unitsquare()), 2,2,3)) + @test Set(boundary_identifiers(a)) == Set([ + MultiBlockBoundary{(1,1,1), bottom}(), + MultiBlockBoundary{(2,1,1), bottom}(), + MultiBlockBoundary{(1,2,1), bottom}(), + MultiBlockBoundary{(2,2,1), bottom}(), + + MultiBlockBoundary{(1,1,3), top}(), + MultiBlockBoundary{(2,1,3), top}(), + MultiBlockBoundary{(1,2,3), top}(), + MultiBlockBoundary{(2,2,3), top}(), + + MultiBlockBoundary{(1,1,1), west}(), + MultiBlockBoundary{(1,2,1), west}(), + MultiBlockBoundary{(1,1,2), west}(), + MultiBlockBoundary{(1,2,2), west}(), + MultiBlockBoundary{(1,1,3), west}(), + MultiBlockBoundary{(1,2,3), west}(), + + MultiBlockBoundary{(2,1,1), east}(), + MultiBlockBoundary{(2,2,1), east}(), + MultiBlockBoundary{(2,1,2), east}(), + MultiBlockBoundary{(2,2,2), east}(), + MultiBlockBoundary{(2,1,3), east}(), + MultiBlockBoundary{(2,2,3), east}(), + + MultiBlockBoundary{(1,1,1), south}(), + MultiBlockBoundary{(2,1,1), south}(), + MultiBlockBoundary{(1,1,2), south}(), + MultiBlockBoundary{(2,1,2), south}(), + MultiBlockBoundary{(1,1,3), south}(), + MultiBlockBoundary{(2,1,3), south}(), + + MultiBlockBoundary{(1,2,1), north}(), + MultiBlockBoundary{(2,2,1), north}(), + MultiBlockBoundary{(1,2,2), north}(), + MultiBlockBoundary{(2,2,2), north}(), + MultiBlockBoundary{(1,2,3), north}(), + MultiBlockBoundary{(2,2,3), north}(), + ]) + end end @testset "UnstructuredAtlas" begin - @test_broken false + @testset "Constructors" begin + c1 = Chart(identity, unitsquare()) + c2 = Chart(x->2x, unitsquare()) + cn = [ + (MultiBlockBoundary{1, east}(), MultiBlockBoundary{2, west}()), + (MultiBlockBoundary{1, north}(), MultiBlockBoundary{3, west}()), + (MultiBlockBoundary{2, north}(), MultiBlockBoundary{3, south}()), + ] + + @test UnstructuredAtlas([c1, c1, c1], cn) isa UnstructuredAtlas + @test UnstructuredAtlas([c1, c2, c1, c2], cn) isa UnstructuredAtlas + + + cn = @SVector[ + (MultiBlockBoundary{1, east}(), MultiBlockBoundary{2, west}()), + (MultiBlockBoundary{1, north}(), MultiBlockBoundary{3, west}()), + (MultiBlockBoundary{2, north}(), MultiBlockBoundary{3, south}()), + ] + @test UnstructuredAtlas(@SVector[c1, c1, c1], cn) isa UnstructuredAtlas + @test UnstructuredAtlas(@SVector[c1, c2, c1, c2], cn) isa UnstructuredAtlas + end + + @testset "Getters" begin + c = Chart(identity, unitsquare()) + cn = [ + (MultiBlockBoundary{1, east}(), MultiBlockBoundary{2, west}()), + (MultiBlockBoundary{1, north}(), MultiBlockBoundary{3, west}()), + (MultiBlockBoundary{2, north}(), MultiBlockBoundary{3, south}()), + ] + + a = UnstructuredAtlas([c, c, c], cn) + + @test charts(a) == [c,c,c] + @test connections(a) == cn + end + + @testset "boundary_identifiers" begin + c = Chart(identity, unitsquare()) + cn = [ + (MultiBlockBoundary{1, east}(), MultiBlockBoundary{2, west}()), + (MultiBlockBoundary{1, north}(), MultiBlockBoundary{3, west}()), + (MultiBlockBoundary{2, north}(), MultiBlockBoundary{3, south}()), + ] + + a = UnstructuredAtlas([c, c, c], cn) + + @test Set(boundary_identifiers(a)) == Set([ + MultiBlockBoundary{1, west}(), + MultiBlockBoundary{1, south}(), + MultiBlockBoundary{2, south}(), + MultiBlockBoundary{2, east}(), + MultiBlockBoundary{3, north}(), + MultiBlockBoundary{3, east}(), + ]) + + end end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test/Grids/multiblockgrids_test.jl Fri Feb 07 22:38:39 2025 +0100 @@ -0,0 +1,10 @@ +using Diffinitive.Grids + +@testset "MultiBlockBoundary" begin + @test MultiBlockBoundary{1,UpperBoundary}() isa BoundaryIdentifier + + @test grid_id(MultiBlockBoundary{1,UpperBoundary}()) == 1 + + @test boundary_id(MultiBlockBoundary{1,UpperBoundary}()) == UpperBoundary() + +end
--- a/test/Grids/parameter_space_test.jl Mon Feb 03 15:44:07 2025 +0100 +++ b/test/Grids/parameter_space_test.jl Fri Feb 07 22:38:39 2025 +0100 @@ -1,5 +1,7 @@ using Test +using Diffinitive.Grids + @testset "ParameterSpace" begin @test ndims(HyperBox([1,1], [2,2])) == 2 @test ndims(unittetrahedron()) == 3 @@ -18,6 +20,8 @@ @test unitinterval(Int) isa Interval{Int} @test unitinterval(Int) == Interval(0,1) @test limits(unitinterval(Int)) == (0,1) + + @test boundary_identifiers(unitinterval()) == (LowerBoundary(), UpperBoundary()) end @testset "HyperBox" begin @@ -38,6 +42,23 @@ @test unithyperbox(4) isa HyperBox{Float64,4} @test limits(unithyperbox(4)) == ([0,0,0,0],[1,1,1,1]) + + + @test boundary_identifiers(unitsquare()) == [ + CartesianBoundary{1,LowerBoundary}(), + CartesianBoundary{1,UpperBoundary}(), + CartesianBoundary{2,LowerBoundary}(), + CartesianBoundary{2,UpperBoundary}(), + ] + + @test boundary_identifiers(unitcube()) == [ + CartesianBoundary{1,LowerBoundary}(), + CartesianBoundary{1,UpperBoundary}(), + CartesianBoundary{2,LowerBoundary}(), + CartesianBoundary{2,UpperBoundary}(), + CartesianBoundary{3,LowerBoundary}(), + CartesianBoundary{3,UpperBoundary}(), + ] end @testset "Simplex" begin