Mercurial > repos > public > sbplib_julia
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)