Mercurial > repos > public > sbplib_julia
changeset 1903:04c251bccbd4 feature/grids/manifolds
Merge feature/grids/parameter_spaces
author | Jonatan Werpers <jonatan@werpers.com> |
---|---|
date | Sat, 01 Feb 2025 22:17:39 +0100 |
parents | edee7d677efb (diff) f93ba5832146 (current diff) |
children | e97f4352b8d0 |
files | src/Grids/Grids.jl src/Grids/equidistant_grid.jl src/Grids/manifolds.jl src/Grids/mapped_grid.jl test/Grids/equidistant_grid_test.jl test/Grids/manifolds_test.jl test/Grids/mapped_grid_test.jl |
diffstat | 9 files changed, 315 insertions(+), 1 deletions(-) [+] |
line wrap: on
line diff
diff -r f93ba5832146 -r 04c251bccbd4 Manifest.toml --- a/Manifest.toml Sat Feb 01 22:12:53 2025 +0100 +++ b/Manifest.toml Sat Feb 01 22:17:39 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"
diff -r f93ba5832146 -r 04c251bccbd4 Project.toml --- a/Project.toml Sat Feb 01 22:12:53 2025 +0100 +++ b/Project.toml Sat Feb 01 22:17:39 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"
diff -r f93ba5832146 -r 04c251bccbd4 ext/DiffinitivePlotsExt.jl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/ext/DiffinitivePlotsExt.jl Sat Feb 01 22:17:39 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 + +
diff -r f93ba5832146 -r 04c251bccbd4 src/Grids/Grids.jl --- a/src/Grids/Grids.jl Sat Feb 01 22:12:53 2025 +0100 +++ b/src/Grids/Grids.jl Sat Feb 01 22:17:39 2025 +0100 @@ -24,6 +24,15 @@ export unittetrahedron export unitsimplex +export Chart + +export Atlas +export charts +export connections +export CartesianAtlas + +export parameterspace + # Grid export Grid export coordinate_size @@ -64,6 +73,7 @@ export metric_tensor include("parameter_space.jl") +include("manifolds.jl") include("grid.jl") include("tensor_grid.jl") include("equidistant_grid.jl")
diff -r f93ba5832146 -r 04c251bccbd4 src/Grids/equidistant_grid.jl --- a/src/Grids/equidistant_grid.jl Sat Feb 01 22:12:53 2025 +0100 +++ b/src/Grids/equidistant_grid.jl Sat Feb 01 22:17:39 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?
diff -r f93ba5832146 -r 04c251bccbd4 src/Grids/grid.jl --- a/src/Grids/grid.jl Sat Feb 01 22:12:53 2025 +0100 +++ b/src/Grids/grid.jl Sat Feb 01 22:17:39 2025 +0100 @@ -122,6 +122,8 @@ """ function boundary_identifiers end +# TBD: Boundary identifiers for charts and atlases? + """ boundary_grid(g::Grid, id::BoundaryIdentifier)
diff -r f93ba5832146 -r 04c251bccbd4 src/Grids/manifolds.jl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/src/Grids/manifolds.jl Sat Feb 01 22:17:39 2025 +0100 @@ -0,0 +1,144 @@ +""" + 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? + + +# TBD: Should Charts, parameterspaces, Atlases, have boundary names? + +""" + 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 <: Atlas + charts::Matrix{Chart} +end + +charts(a::CartesianAtlas) = a.charts +connections(a::CartesianAtlas) = nothing + +struct UnstructuredAtlas <: Atlas + charts::Vector{Chart} + connections +end + +charts(a::UnstructuredAtlas) = a.charts +connections(a::UnstructuredAtlas) = nothing + + +### +# Geometry +### + +abstract type Curve end +abstract type Surface end + + +struct Line{PT} <: Curve + p::PT + tangent::PT +end + +(c::Line)(s) = c.p + s*c.tangent + + +struct LineSegment{PT} <: Curve + a::PT + b::PT +end + +(c::LineSegment)(s) = (1-s)*c.a + s*c.b + + +function linesegments(ps...) + return [LineSegment(ps[i], ps[i+1]) for i ∈ 1:length(ps)-1] +end + + +function polygon_edges(ps...) + n = length(ps) + return [LineSegment(ps[i], ps[mod1(i+1,n)]) for i ∈ eachindex(ps)] +end + +struct Circle{T,PT} <: Curve + c::PT + r::T +end + +function (C::Circle)(θ) + (;c, r) = C + c + r*@SVector[cos(θ), sin(θ)] +end + +struct TransfiniteInterpolationSurface{T1,T2,T3,T4} <: Surface + c₁::T1 + c₂::T2 + c₃::T3 + c₄::T4 +end + +function (s::TransfiniteInterpolationSurface)(u,v) + c₁, c₂, c₃, c₄ = s.c₁, s.c₂, s.c₃, s.c₄ + P₀₀ = c₁(0) + P₁₀ = c₂(0) + P₁₁ = c₃(0) + P₀₁ = c₄(0) + return (1-v)*c₁(u) + u*c₂(v) + v*c₃(1-u) + (1-u)*c₄(1-v) - ( + (1-u)*(1-v)*P₀₀ + u*(1-v)*P₁₀ + u*v*P₁₁ + (1-u)*v*P₀₁ + ) +end + +function (s::TransfiniteInterpolationSurface)(ξ̄::AbstractArray) + s(ξ̄...) +end + +# TODO: Implement jacobian() for the different mapping helpers +
diff -r f93ba5832146 -r 04c251bccbd4 test/Grids/equidistant_grid_test.jl --- a/test/Grids/equidistant_grid_test.jl Sat Feb 01 22:12:53 2025 +0100 +++ b/test/Grids/equidistant_grid_test.jl Sat Feb 01 22:17:39 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
diff -r f93ba5832146 -r 04c251bccbd4 test/Grids/manifolds_test.jl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/test/Grids/manifolds_test.jl Sat Feb 01 22:17:39 2025 +0100 @@ -0,0 +1,69 @@ +using Test + +using Diffinitive.Grids +using Diffinitive.RegionIndices +using Diffinitive.LazyTensors + +@testset "Chart" begin + c = Chart(x->2x, 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]) +end + +@testset "CartesianAtlas" begin + c = Chart(identity, unitsquare()) + + 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}()), + ) + ] +end + +@testset "UnstructuredAtlas" begin + @test_broken false +end + +@testset "Line" begin + @test_broken false +end + +@testset "LineSegment" begin + @test_broken false +end + +@testset "linesegments" begin + @test_broken false +end + +@testset "polygon_edges" begin + @test_broken false +end + +@testset "Circle" begin + @test_broken false +end + +@testset "TransfiniteInterpolationSurface" begin + @test_broken false +end