changeset 1669:a750992d870a feature/grids/min_spacing

Merge default
author Jonatan Werpers <jonatan@werpers.com>
date Sat, 29 Jun 2024 17:06:14 +0200
parents 5f348cc5598e (current diff) bcc2be934326 (diff)
children 1212ba23ee32
files
diffstat 7 files changed, 122 insertions(+), 32 deletions(-) [+]
line wrap: on
line diff
--- a/Project.toml	Wed Jun 26 11:30:38 2024 +0200
+++ b/Project.toml	Sat Jun 29 17:06:14 2024 +0200
@@ -8,5 +8,11 @@
 TOML = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
 TiledIteration = "06e1c1a7-607b-532d-9fad-de7d9aa2abac"
 
+[weakdeps]
+Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+
+[extensions]
+SbplibMakieExt = "Makie"
+
 [compat]
 julia = "1.5"
--- a/docs/src/grids_and_grid_functions.md	Wed Jun 26 11:30:38 2024 +0200
+++ b/docs/src/grids_and_grid_functions.md	Sat Jun 29 17:06:14 2024 +0200
@@ -7,6 +7,24 @@
 ## Interface for grids
 All grids are expected to work as a grid function for the coordinate function, and thus implements Julia's Indexing- and Iteration-interfaces. Notably they are *not* abstract arrays because that inteface is too restrictive for the types of grids we wish to implement.
 
+
+## Plotting
+Plotting of grids and grid functions is supported through a package extension with Makie.jl.
+
+For grids we have:
+* `plot(::Grid{<:Any,2})` (same as `lines`)
+* `lines(::Grid{<:Any,2})`
+* `scatter(::Grid{<:Any,2})`
+
+For 1D grid functions we have:
+* `plot(::Grid{<:Any,1}, ::AbstractVector)` (same as `lines`)
+* `lines(::Grid{<:Any,1}, ::AbstractVector)`
+* `scatter(::Grid{<:Any,1}, ::AbstractVector)`
+
+For 2D grid functions we have:
+* `plot(::Grid{<:Any,2}, ::AbstractArray{<:Any,2})` (constructs a 2d mesh)
+* `surface(::Grid{<:Any,2}, ::AbstractArray{<:Any,2})`
+
 ## To write about
 <!-- # TODO: -->
 * Grid functions
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/ext/SbplibMakieExt.jl	Sat Jun 29 17:06:14 2024 +0200
@@ -0,0 +1,88 @@
+module SbplibMakieExt
+
+using Sbplib.Grids
+using Makie
+using StaticArrays
+
+
+function verticies_and_faces_and_values(g::Grid{<:Any,2}, gf::AbstractArray{<:Any, 2})
+    ps = map(Tuple, g)[:]
+    values = gf[:]
+    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
+
+
+## Grids
+
+Makie.convert_arguments(::Type{<:Scatter}, g::Grid) = (reshape(map(Point,g),:),) # (map(Point,collect(g)[:]),)
+function Makie.convert_arguments(::Type{<:Lines}, g::Grid{<:Any,2})
+    M = collect(g)
+
+    function cat_with_NaN(a,b)
+        vcat(a,[@SVector[NaN,NaN]],b)
+    end
+
+    xlines = reduce(cat_with_NaN, eachrow(M))
+    ylines = reduce(cat_with_NaN, eachcol(M))
+
+    return (cat_with_NaN(xlines,ylines),)
+end
+
+Makie.plot!(plot::Plot(Grid{<:Any,2})) = lines!(plot, plot.attributes, plot[1])
+
+
+## Grid functions
+
+### 1D
+function Makie.convert_arguments(::Type{<:Lines}, g::Grid{<:Any,1}, gf::AbstractArray{<:Any, 1})
+    (collect(g), gf)
+end
+
+function Makie.convert_arguments(::Type{<:Scatter}, g::Grid{<:Any,1}, gf::AbstractArray{<:Any, 1})
+    (collect(g), gf)
+end
+
+Makie.plot!(plot::Plot(Grid{<:Any,1}, AbstractArray{<:Any,1})) = lines!(plot, plot.attributes, plot[1], plot[2])
+
+### 2D
+function Makie.convert_arguments(::Type{<:Surface}, g::Grid{<:Any,2}, gf::AbstractArray{<:Any, 2})
+    (getindex.(g,1), getindex.(g,2), gf)
+end
+
+function Makie.plot!(plot::Plot(Grid{<:Any,2},AbstractArray{<:Any, 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, plot.attributes, v, f;
+        color=c,
+        shading = NoShading,
+    )
+end
+# TBD: Can we define `mesh` instead of the above function and then forward plot! to that?
+
+function Makie.convert_arguments(::Type{<:Scatter}, g::Grid{<:Any,2}, gf::AbstractArray{<:Any, 2})
+    ps = map(g,gf) do (x,y), z
+        @SVector[x,y,z]
+    end
+    (reshape(ps,:),)
+end
+
+end
--- a/src/RegionIndices/RegionIndices.jl	Wed Jun 26 11:30:38 2024 +0200
+++ b/src/RegionIndices/RegionIndices.jl	Sat Jun 29 17:06:14 2024 +0200
@@ -13,7 +13,7 @@
     Index{R,T}(i::T) where {R<:Region,T<:Integer} = new{R,T}(i)
     Index{R}(i::T) where {R<:Region,T<:Integer} = new{R,T}(i)
     Index(i::T, ::Type{R}) where {R<:Region,T<:Integer} = Index{R,T}(i)
-    Index(t::Tuple{T, DataType}) where {R<:Region,T<:Integer} = Index{t[2],T}(t[1]) # TBD: This is not very specific in what types are allowed in t[2]. Can this be fixed?
+    Index(t::Tuple{T, DataType}) where T<:Integer = Index{t[2],T}(t[1]) # TBD: This is not very specific in what types are allowed in t[2]. Can this be fixed?
 end
 
 export Index
--- a/test/LazyTensors/lazy_tensor_operations_test.jl	Wed Jun 26 11:30:38 2024 +0200
+++ b/test/LazyTensors/lazy_tensor_operations_test.jl	Sat Jun 29 17:06:14 2024 +0200
@@ -4,13 +4,13 @@
 
 using Tullio
 
-struct DummyMapping{T,R,D} <: LazyTensor{T,R,D} end
+struct TransposableDummyMapping{T,R,D} <: LazyTensor{T,R,D} end
 
-LazyTensors.apply(m::DummyMapping{T,R}, v, I::Vararg{Any,R}) where {T,R} = :apply
-LazyTensors.apply_transpose(m::DummyMapping{T,R,D}, v, I::Vararg{Any,D}) where {T,R,D} = :apply_transpose
+LazyTensors.apply(m::TransposableDummyMapping{T,R}, v, I::Vararg{Any,R}) where {T,R} = :apply
+LazyTensors.apply_transpose(m::TransposableDummyMapping{T,R,D}, v, I::Vararg{Any,D}) where {T,R,D} = :apply_transpose
 
-LazyTensors.range_size(m::DummyMapping) = :range_size
-LazyTensors.domain_size(m::DummyMapping) = :domain_size
+LazyTensors.range_size(m::TransposableDummyMapping) = :range_size
+LazyTensors.domain_size(m::TransposableDummyMapping) = :domain_size
 
 
 struct SizeDoublingMapping{T,R,D} <: LazyTensor{T,R,D}
@@ -24,7 +24,7 @@
 
 
 @testset "Mapping transpose" begin
-    m = DummyMapping{Float64,2,3}()
+    m = TransposableDummyMapping{Float64,2,3}()
     @test m' isa LazyTensor{Float64, 3,2}
     @test m'' == m
     @test apply(m',zeros(Float64,(0,0)), 0, 0, 0) == :apply_transpose
--- a/test/SbpOperators/volumeops/derivatives/dissipation_test.jl	Wed Jun 26 11:30:38 2024 +0200
+++ b/test/SbpOperators/volumeops/derivatives/dissipation_test.jl	Sat Jun 29 17:06:14 2024 +0200
@@ -13,20 +13,9 @@
 using Sbplib.SbpOperators: dissipation_lower_closure_stencils,dissipation_upper_closure_stencils
 using Sbplib.SbpOperators: dissipation_transpose_lower_closure_stencils, dissipation_transpose_upper_closure_stencils
 
-"""
-    monomial(x,k)
-
-Evaluates ``x^k/k!` with the convetion that it is ``0`` for all ``k<0``.
-Has the property that ``d/dx monomial(x,k) = monomial(x,k-1)``
-"""
-function monomial(x,k)
-    if k < 0
-        return zero(x)
-    end
-    x^k/factorial(k)
-end
 
 @testset "undivided_skewed04" begin
+    monomial(x,k) = k < 0 ? zero(x) : x^k/factorial(k)
     g = equidistant_grid(0., 11., 20)
     D,Dįµ€ = undivided_skewed04(g, 1)
 
--- a/test/SbpOperators/volumeops/derivatives/first_derivative_test.jl	Wed Jun 26 11:30:38 2024 +0200
+++ b/test/SbpOperators/volumeops/derivatives/first_derivative_test.jl	Sat Jun 29 17:06:14 2024 +0200
@@ -7,19 +7,6 @@
 
 using Sbplib.SbpOperators: closure_size, Stencil, VolumeOperator
 
-"""
-    monomial(x,k)
-
-Evaluates ``x^k/k!` with the convetion that it is ``0`` for all ``k<0``.
-Has the property that ``d/dx monomial(x,k) = monomial(x,k-1)``
-"""
-function monomial(x,k)
-    if k < 0
-        return zero(x)
-    end
-    x^k/factorial(k)
-end
-
 @testset "first_derivative" begin
     @testset "Constructors" begin
         stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order=2)
@@ -39,6 +26,8 @@
     @testset "Accuracy conditions" begin
         N = 20
         g = equidistant_grid(0//1, 2//1, N)
+
+        monomial(x,k) = k < 0 ? zero(x) : x^k/factorial(k)
         @testset for order ∈ [2,4]
             stencil_set = read_stencil_set(sbp_operators_path()*"standard_diagonal.toml"; order)
             D₁ = first_derivative(g, stencil_set)