Mercurial > repos > public > sbplib_julia

changeset 1587:aef3827ef522 feature/sbp_operators/laplace_curvilinear

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Merge feature/grids/manifolds
author Jonatan Werpers <jonatan@werpers.com>
date Fri, 26 Apr 2024 22:49:36 +0200
parents d359d0d7fb12 (current diff) d4a6f9effcdd (diff)
children f6774e98d223
files src/Grids/Grids.jl
diffstat 9 files changed, 173 insertions(+), 20 deletions(-) [+]
line wrap: on
line diff
--- a/Project.toml	Fri Apr 26 08:53:32 2024 +0200
+++ b/Project.toml	Fri Apr 26 22:49:36 2024 +0200
@@ -8,5 +8,13 @@
 TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 TiledIteration = "06e1c1a7-607b-532d-9fad-de7d9aa2abac"
 
+[weakdeps]
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+
+[extensions]
+SbplibPlotsExt = "Plots"
+SbplibMakieExt = "Makie"
+
 [compat]
 julia = "1.5"
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/ext/SbplibMakieExt.jl	Fri Apr 26 22:49:36 2024 +0200
@@ -0,0 +1,71 @@
+module SbplibMakieExt
+
+using Sbplib.Grids
+using Makie
+
+
+function verticies_and_faces_and_values(g::Grid{<:Any,2}, gf::AbstractArray{<:Any, 2})
+    ps = map(Tuple, g)[:]
+    values = gf[:]
+
+    N = length(ps)
+
+    faces = Vector{NTuple{3,Int}}()
+
+    n,m = size(g)
+    Li = LinearIndices((1:n, 1:m))
+    for i ∈ 1:n-1, j = 1:m-1
+
+        # Add point in the middle of the patch to preserve symmetries
+        push!(ps, Tuple((g[i,j] + g[i+1,j] + g[i+1,j+1] + g[i,j+1])/4))
+        push!(values, (gf[i,j] + gf[i+1,j] + gf[i+1,j+1] + gf[i,j+1])/4)
+
+        push!(faces, (Li[i,j],     Li[i+1,j],   length(ps)))
+        push!(faces, (Li[i+1,j],   Li[i+1,j+1], length(ps)))
+        push!(faces, (Li[i+1,j+1], Li[i,j+1],   length(ps)))
+        push!(faces, (Li[i,j+1],   Li[i,j],     length(ps)))
+    end
+
+    verticies = permutedims(reinterpret(reshape,eltype(eltype(ps)), ps))
+    faces = permutedims(reinterpret(reshape,Int, faces))
+
+    return verticies, faces, values
+end
+
+
+function make_plot(g,gf)
+    v,f,c = verticies_and_faces_and_values(g,gf)
+    mesh(v,f,color=c,
+        shading = NoShading,
+    )
+end
+
+# scatter(collect(g)[:])
+
+function Makie.surface(g::Grid{<:Any,2}, gf::AbstractArray{<:Any, 2}; kwargs...)
+    surface(getindex.(g,1), getindex.(g,2), gf;
+        shading = NoShading,
+        kwargs...,
+    )
+end
+
+function Makie.mesh(g::Grid{<:Any,2}, gf::AbstractArray{<:Any, 2}; kwargs...)
+    v,f,c = verticies_and_faces_and_values(g, gf)
+    mesh(v,f,color=c,
+        shading = NoShading,
+        kwargs...,
+    )
+end
+
+function Makie.plot!(plot::Plot(Grid{<:Any,2},AbstractArray{<:Any, 2}))
+    # TODO: How to handle kwargs?
+    # v,f,c = verticies_and_faces_and_values(plot[1], plot[2])
+    r = @lift verticies_and_faces_and_values($(plot[1]), $(plot[2]))
+    v,f,c = (@lift $r[1]), (@lift $r[2]), (@lift $r[3])
+    mesh!(plot, v, f, color=c,
+        shading = NoShading,
+    )
+end
+
+Makie.convert_arguments(::Type{<:Scatter}, g::Grid) = (map(Tuple,collect(g)[:]),)
+end
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/ext/SbplibPlotsExt.jl	Fri Apr 26 22:49:36 2024 +0200
@@ -0,0 +1,72 @@
+module SbplibPlotsExt
+
+using Sbplib.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 Apr 26 08:53:32 2024 +0200
+++ b/src/Grids/Grids.jl	Fri Apr 26 22:49:36 2024 +0200
@@ -5,7 +5,7 @@
 using Sbplib.LazyTensors
 using StaticArrays
 
-
+export ParameterSpace
 export HyperBox
 export Simplex
 export Interval
--- a/src/Grids/equidistant_grid.jl	Fri Apr 26 08:53:32 2024 +0200
+++ b/src/Grids/equidistant_grid.jl	Fri Apr 26 22:49:36 2024 +0200
@@ -130,6 +130,15 @@
 
 equidistant_grid(hb::HyperBox, dims::Vararg{Int}) = equidistant_grid(hb.a, hb.b, dims...)
 
+function equidistant_grid(c::Chart, dims::Vararg{Int})
+    lg = equidistant_grid(parameterspace(c), dims...)
+    return MappedGrid(
+        lg,
+        map(c,lg),
+        map(ξ->jacobian(c, ξ), lg),
+    )
+end
+
 
 CartesianBoundary{D,BID} = TensorGridBoundary{D,BID} # TBD: What should we do about the naming of this boundary?
 
--- a/src/Grids/mapped_grid.jl	Fri Apr 26 08:53:32 2024 +0200
+++ b/src/Grids/mapped_grid.jl	Fri Apr 26 22:49:36 2024 +0200
@@ -61,14 +61,6 @@
     )
 end
 
-function mapped_grid(c::Chart, size...)
-    lg = equidistant_grid(parameterspace(c), size...)
-    return MappedGrid(
-        lg,
-        map(c,lg),
-        map(ξ->jacobian(c, ξ), lg),
-    )
-end
 
 function jacobian_determinant(g::MappedGrid)
     return map(jacobian(g)) do ∂x∂ξ
--- a/test/Grids/equidistant_grid_test.jl	Fri Apr 26 08:53:32 2024 +0200
+++ b/test/Grids/equidistant_grid_test.jl	Fri Apr 26 22:49:36 2024 +0200
@@ -151,6 +151,16 @@
 
         @test equidistant_grid(ps, 3,4) == equidistant_grid((0,0), (2,1), 3,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
 
 
--- a/test/Grids/manifolds_test.jl	Fri Apr 26 08:53:32 2024 +0200
+++ b/test/Grids/manifolds_test.jl	Fri Apr 26 22:49:36 2024 +0200
@@ -12,6 +12,7 @@
 end
 
 @testset "HyperBox" begin
+    @test HyperBox{<:Any, 2} <: ParameterSpace{2}
     @test HyperBox([1,1], [2,2]) isa HyperBox{Int, 2}
 
     @test limits(HyperBox([1,2], [3,4])) == ([1,2], [3,4])
@@ -35,6 +36,7 @@
 end
 
 @testset "Simplex" begin
+    @test Simplex{<:Any, 3} <: ParameterSpace{3}
     @test Simplex([1,2], [3,4]) isa Simplex{Int, 2}
     @test Simplex([1,2,3], [4,5,6],[1,1,1]) isa Simplex{Int, 3}
 
--- a/test/Grids/mapped_grid_test.jl	Fri Apr 26 08:53:32 2024 +0200
+++ b/test/Grids/mapped_grid_test.jl	Fri Apr 26 22:49:36 2024 +0200
@@ -183,15 +183,4 @@
     lg = equidistant_grid((0,0), (1,1), 10, 11)
     @test logicalgrid(mg) == lg
     @test collect(mg) == map(x̄, lg)
-
-
-    @testset "mapped_grid(::Chart)" begin
-        c = Chart(unitsquare()) do (ξ,η)
-            @SVector[2ξ, 3η]
-        end
-        Grids.jacobian(c::typeof(c), ξ̄) = @SMatrix[2 0; 0 3]
-
-        @test mapped_grid(c, 5, 4) isa Grid
-    end
-
 end