Mercurial > repos > public > sbplib_julia
changeset 1943:48c49c04f3b2 feature/grids/manifolds
Merge feature/grids/parameter_spaces
| author | Jonatan Werpers <jonatan@werpers.com> |
|---|---|
| date | Fri, 07 Feb 2025 15:29:40 +0100 |
| parents | b8395f69ad80 (diff) 15be190a40cd (current diff) |
| children | 2a923e673cfc |
| files | |
| diffstat | 13 files changed, 498 insertions(+), 1 deletions(-) [+] |
line wrap: on
line diff
--- a/Manifest.toml Fri Feb 07 15:29:05 2025 +0100 +++ b/Manifest.toml Fri Feb 07 15:29:40 2025 +0100 @@ -2,7 +2,7 @@ julia_version = "1.11.3" manifest_format = "2.0" -project_hash = "383b10bc09a08d519764abf167e40dae29b16b2e" +project_hash = "51374910f6b7831392791bcca0079860f1023f47" [[deps.Artifacts]] uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33"
--- a/Project.toml Fri Feb 07 15:29:05 2025 +0100 +++ b/Project.toml Fri Feb 07 15:29:40 2025 +0100 @@ -9,6 +9,7 @@ TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76" [weakdeps] +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" SparseArrayKit = "a9a3c162-d163-4c15-8926-b8794fbefed2" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" @@ -16,6 +17,7 @@ [extensions] DiffinitiveMakieExt = "Makie" +DiffinitivePlotsExt = "Plots" DiffinitiveSparseArrayKitExt = ["SparseArrayKit", "Tokens"] DiffinitiveSparseArraysExt = ["SparseArrays", "Tokens"] @@ -24,6 +26,7 @@ StaticArrays = "1.0" TOML = "1.0" Makie = "0.21" +Plots = "1.40" SparseArrayKit = "0.4" Tokens = "0.1.1" SparseArrays = "1.10"
--- a/docs/make.jl Fri Feb 07 15:29:05 2025 +0100 +++ b/docs/make.jl Fri Feb 07 15:29:40 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 15:29:40 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? --/>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/ext/DiffinitivePlotsExt.jl Fri Feb 07 15:29:40 2025 +0100 @@ -0,0 +1,72 @@ +module DiffinitivePlotsExt + +using Diffinitive.Grids +using Plots + +@recipe f(::Type{<:Grid}, g::Grid) = map(Tuple,g)[:] + +@recipe function f(c::Chart{2,<:Rectangle}, n=5, m=n; draw_border=true, bordercolor=1) + Ξ = parameterspace(c) + ξs = range(limits(Ξ,1)..., n) + ηs = range(limits(Ξ,2)..., m) + + label := false + seriescolor --> 2 + for ξ ∈ ξs + @series adapted_curve_grid(η->c((ξ,η)),limits(Ξ,1)) + end + + for η ∈ ηs + @series adapted_curve_grid(ξ->c((ξ,η)),limits(Ξ,2)) + end + + if ~draw_border + return + end + + for ξ ∈ limits(Ξ,1) + @series begin + linewidth --> 3 + seriescolor := bordercolor + adapted_curve_grid(η->c((ξ,η)),limits(Ξ,1)) + end + end + + for η ∈ limits(Ξ,2) + @series begin + linewidth --> 3 + seriescolor := bordercolor + adapted_curve_grid(ξ->c((ξ,η)),limits(Ξ,2)) + end + end +end + +function adapted_curve_grid(g, minmax) + t1, _ = PlotUtils.adapted_grid(t->g(t)[1], minmax) + t2, _ = PlotUtils.adapted_grid(t->g(t)[2], minmax) + + ts = sort(vcat(t1,t2)) + + x = map(ts) do t + g(t)[1] + end + y = map(ts) do t + g(t)[2] + end + + return x, y +end + +# get_axis_limits(plt, :x) + + +# ReicpesPipline/src/user_recipe.jl +# @recipe function f(f::FuncOrFuncs{F}) where {F<:Function} + +# @recipe function f(f::Function, xmin::Number, xmax::Number) + +# _scaled_adapted_grid(f, xscale, yscale, xmin, xmax) + +end + +
--- a/src/Grids/Grids.jl Fri Feb 07 15:29:05 2025 +0100 +++ b/src/Grids/Grids.jl Fri Feb 07 15:29:40 2025 +0100 @@ -24,6 +24,17 @@ export unittetrahedron export unitsimplex +export Chart + +export Atlas +export charts +export connections +export boundaries +export CartesianAtlas +export UnstructuredAtlas + +export parameterspace + # Grid export Grid export coordinate_size @@ -55,6 +66,8 @@ export spacing export equidistant_grid +export MultiBlockBoundary + # MappedGrid export MappedGrid @@ -65,6 +78,8 @@ include("parameter_space.jl") include("grid.jl") +include("multiblockgrids.jl") +include("manifolds.jl") include("tensor_grid.jl") include("equidistant_grid.jl") include("zero_dim_grid.jl")
--- a/src/Grids/equidistant_grid.jl Fri Feb 07 15:29:05 2025 +0100 +++ b/src/Grids/equidistant_grid.jl Fri Feb 07 15:29:40 2025 +0100 @@ -154,6 +154,10 @@ equidistant_grid(d::Interval, size::Int) = equidistant_grid(limits(d)..., size) equidistant_grid(hb::HyperBox, dims::Vararg{Int}) = equidistant_grid(limits(hb)..., dims...) +function equidistant_grid(c::Chart, dims::Vararg{Int}) + mapped_grid(c, ξ->jacobian(c,ξ), parameterspace(c), dims...) +end + CartesianBoundary{D,BID} = TensorGridBoundary{D,BID} # TBD: What should we do about the naming of this boundary?
--- a/src/Grids/grid.jl Fri Feb 07 15:29:05 2025 +0100 +++ b/src/Grids/grid.jl Fri Feb 07 15:29:40 2025 +0100 @@ -122,6 +122,8 @@ """ function boundary_identifiers end +# TBD: Boundary identifiers for charts and atlases? + """ boundary_grid(g::Grid, id::BoundaryIdentifier)
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/Grids/manifolds.jl Fri Feb 07 15:29:40 2025 +0100 @@ -0,0 +1,127 @@ +""" + Chart{D} + +A parametrized description of a manifold or part of a manifold. +""" +struct Chart{D, PST<:ParameterSpace{D}, MT} + mapping::MT + parameterspace::PST +end + +Base.ndims(::Chart{D}) where D = D +(c::Chart)(ξ) = c.mapping(ξ) +parameterspace(c::Chart) = c.parameterspace + +""" + jacobian(c::Chart, ξ) + +The jacobian of the mapping evaluated at `ξ`. This defers to the +implementation of `jacobian` for the mapping itself. If no implementation is +available one can easily be specified for either the mapping function or the +chart itself. +```julia +c = Chart(f, ps) +jacobian(f::typeof(f), ξ) = f′(ξ) +``` +or +```julia +c = Chart(f, ps) +jacobian(c::typeof(c),ξ) = f′(ξ) +``` +which will both allow calling `jacobian(c,ξ)`. +""" +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? + +boundaries(c::Chart) = boundaries(parameterspace(c)) + + +""" + Atlas + +A collection of charts and their connections. +Should implement methods for `charts` and `connections`. +""" +abstract type Atlas end + +""" + charts(::Atlas) + +The colloction of charts in the atlas. +""" +function charts end + +""" + connections(::Atlas) + +TBD: What exactly should this return? +""" +function connections end + +struct CartesianAtlas{D,C<:Chart,AT<:AbstractArray{C,D}} <: Atlas + charts::AT +end + +charts(a::CartesianAtlas) = a.charts + +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 + +function boundaries(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 + + +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) = a.connections + +function boundaries(a::UnstructuredAtlas) + bs = MultiBlockBoundary[] + + for c ∈ charts(a) + for b ∈ boundaries(c) + if !any(cn->b∈cn, connections(a)) + push!(bs, b) + end + end + end +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/Grids/multiblockgrids.jl Fri Feb 07 15:29:40 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/test/Grids/equidistant_grid_test.jl Fri Feb 07 15:29:05 2025 +0100 +++ b/test/Grids/equidistant_grid_test.jl Fri Feb 07 15:29:40 2025 +0100 @@ -163,6 +163,16 @@ @test equidistant_grid(unitinterval(),3) == equidistant_grid(0,1,3) @test equidistant_grid(HyperBox((0,),(2,)),4) == equidistant_grid(@SVector[0], @SVector[2], 4) end + + + @testset "equidistant_grid(::Chart)" begin + c = Chart(unitsquare()) do (ξ,η) + @SVector[2ξ, 3η] + end + Grids.jacobian(c::typeof(c), ξ̄) = @SMatrix[2 0; 0 3] + + @test equidistant_grid(c, 5, 4) isa Grid + end end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test/Grids/manifolds_test.jl Fri Feb 07 15:29:40 2025 +0100 @@ -0,0 +1,208 @@ +using Test + +using Diffinitive.Grids +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 + 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 jacobian(c, [3,2]) == [2,2] + + @test_broken Set(boundaries(X,unitsquare())) == Set([east,west,south,north]) +end + +@testset "CartesianAtlas" begin + @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}()), + + (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 "boundaries" begin + # 2D + a = CartesianAtlas(fill(Chart(identity, unitcube()), 2,3)) + @test Set(boundaries(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(boundaries(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 + @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 "boundaries" 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_broken Set(boundaries(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 15:29:40 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