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 src/Grids/mapped_grid.jl
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
diff -r f93ba5832146 -r 04c251bccbd4 test/Grids/mapped_grid_test.jl