From 6608115aba6fcb7e682e84c2d5d67f9b76c012d5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Wed, 24 Jan 2024 18:19:06 +0100 Subject: [PATCH 1/9] Start with iterator based interface --- Project.toml | 2 ++ src/GridVisualizeTools.jl | 8 +++++ src/linearsimplex.jl | 75 +++++++++++++++++++++++++++++++++++++++ src/marching_iterators.jl | 59 ++++++++++++++++++++++++++++++ 4 files changed, 144 insertions(+) create mode 100644 src/linearsimplex.jl create mode 100644 src/marching_iterators.jl diff --git a/Project.toml b/Project.toml index 88dee08..eff6939 100644 --- a/Project.toml +++ b/Project.toml @@ -7,11 +7,13 @@ version = "1.1.0" ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4" Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" [compat] ColorSchemes = "3" Colors = "0.12" DocStringExtensions = "0.8, 0.9" +GeometryBasics = "0.4" StaticArraysCore = "1.4" julia = "1.6" diff --git a/src/GridVisualizeTools.jl b/src/GridVisualizeTools.jl index ac91403..d95dc73 100644 --- a/src/GridVisualizeTools.jl +++ b/src/GridVisualizeTools.jl @@ -5,6 +5,7 @@ import ColorSchemes using DocStringExtensions: SIGNATURES, TYPEDEF, TYPEDSIGNATURES using StaticArraysCore: SVector +using GeometryBasics: Point include("colors.jl") export region_cmap, bregion_cmap, rgbtuple @@ -12,9 +13,16 @@ export region_cmap, bregion_cmap, rgbtuple include("extraction.jl") export extract_visible_cells3D, extract_visible_bfaces3D +include("linearsimplex.jl") +export LinearSimplex, LinearSimplexIterator +export LinearEdge, LinearTriangle, LinearTetrahedron +export testloop + include("marching.jl") +include("marching_iterators.jl") export marching_tetrahedra, marching_triangles + include("markerpoints.jl") export markerpoints diff --git a/src/linearsimplex.jl b/src/linearsimplex.jl new file mode 100644 index 0000000..52ba6e0 --- /dev/null +++ b/src/linearsimplex.jl @@ -0,0 +1,75 @@ +struct LinearSimplex{D,N,Tv} + points::SVector{N,Point{D,Tv}} + values::SVector{N,Tv} +end + +values(s::LinearSimplex)=s.values +points(s::LinearSimplex)=s.points + + +function LinearSimplex(::Type{Val{D}},points::Union{Vector,Tuple}, values::Union{Vector,Tuple}) where {D} + spoints=SVector{D+1,Point{D,Float32}}(points...) + svalues=SVector{D+1,Float32}(values...) + LinearSimplex(spoints,svalues) +end + +function LinearSimplex(::Type{Val{1}},points::AbstractMatrix, values::AbstractVector) + @views spoints=SVector{2,Point{1,Float32}}(points[:,1],points[:,2]) + svalues=SVector{2,Float32}(values[1],values[2]) + LinearSimplex(spoints,svalues) +end + +function LinearSimplex(::Type{Val{2}},points::AbstractMatrix, values::AbstractVector) + @views spoints=SVector{3,Point{2,Float32}}(points[:,1],points[:,2],points[:,3]) + svalues=SVector{3,Float32}(values[1],values[2],values[3]) + LinearSimplex(spoints,svalues) +end + +function LinearSimplex(::Type{Val{3}},points::AbstractMatrix, values::AbstractVector) + @views spoints=SVector{4,Point{3,Float32}}(points[:,1],points[:,2],points[:,3],points[:,4]) + svalues=SVector{4,Float32}(values[1],values[2],values[3],values[4]) + LinearSimplex(spoints,svalues) +end + +LinearEdge(points,values)=LinearSimplex(Val{1},points,values) +LinearTriangle(points,values)=LinearSimplex(Val{2},points,values) +LinearTetrahedron(points,values)=LinearSimplex(Val{3},points,values) + + +""" + abstract type LinearSimplexIterator{D} + +Iterator over D-dimensional linear simplices. + +Any subtype `TSub` should comply with the [iteration interface](https://docs.julialang.org/en/v1/manual/interfaces/#man-interface-iteration) +and implement +```julia + Base.iterate(vs::TSub{D}, state):: (LinearSimplex{D}, state) where D +``` +The optional methods (in particular, size, length) are not needed. +""" +abstract type LinearSimplexIterator{D} end + +""" + testloop(iterators::AbstractVector{T}) where T<:LinearSimplexIterator + +Sum up all values passed via the iterator. Designed for testing iterators. + +Useful for: +- Consistency test - e.g. when passing constant ones, the sum can be calculated by + other means and compared to the return value of testloop +- Allocation test. `testloop` itself does not allocate, so any allocations from a + call to `testloop` hint at possibilities to improve an iterator implementation. +""" +function testloop(iterators::AbstractVector{T}) where T<:LinearSimplexIterator + threads=map(iterators) do iterator + begin + local x=0.0 + for vt in iterator + x+=sum(vt.values) + end + x + end + end + sum(fetch.(threads)) +end diff --git a/src/marching_iterators.jl b/src/marching_iterators.jl new file mode 100644 index 0000000..338fc44 --- /dev/null +++ b/src/marching_iterators.jl @@ -0,0 +1,59 @@ +function pushintersection!(intersection_points,triangle,levels) + f = values(triangle) + coord=points(triangle) + + (n1, n2, n3) = (1, 2, 3) + + f[1] <= f[2] ? (n1, n2) = (1, 2) : (n1, n2) = (2, 1) + f[n2] <= f[3] ? n3 = 3 : (n2, n3) = (3, n2) + f[n1] > f[n2] ? (n1, n2) = (n2, n1) : nothing + + dx31 = coord[n3][1] - coord[n1][1] + dx21 = coord[n2][1] - coord[n1][1] + dx32 = coord[n3][1] - coord[n2][1] + + dy31 = coord[n3][2] - coord[n1][2] + dy21 = coord[n2][2] - coord[n1][2] + dy32 = coord[n3][2] - coord[n2][2] + + df31 = f[n3] != f[n1] ? 1 / (f[n3] - f[n1]) : 0.0 + df21 = f[n2] != f[n1] ? 1 / (f[n2] - f[n1]) : 0.0 + df32 = f[n3] != f[n2] ? 1 / (f[n3] - f[n2]) : 0.0 + + for level ∈ levels + if (f[n1] <= level) && (level < f[n3]) + α = (level - f[n1]) * df31 + x1 = coord[n1][1] + α * dx31 + y1 = coord[n1][2] + α * dy31 + + if (level < f[n2]) + α = (level - f[n1]) * df21 + x2 = coord[n1][1] + α * dx21 + y2 = coord[n1][2] + α * dy21 + else + α = (level - f[n2]) * df32 + x2 = coord[n2][1] + α * dx32 + y2 = coord[n2][2] + α * dy32 + end + push!(intersection_points, SVector{2, Tc}((x1, y1))) + push!(intersection_points, SVector{2, Tc}((x2, y2))) + end + end +end + + +function marching_triangles(triangle_iterators::Vector{LinearSimplexIterator{2}}, + levels; + Tc = Float32, + Tp = SVector{2, Tc}) + threads=map(triangle_iterators) do triangles + Threads.@spawn begin + local intersection_points = Vector{Tp}(undef, 0) + for triangle in triangles + pushintersection!(intersection_points,triangle) + end + intersection_points + end + end + fetch.(threads) +end From 18eadb065a422db05d2a957079a4dc5d489aade0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Wed, 24 Jan 2024 19:27:20 +0100 Subject: [PATCH 2/9] don't spawn if nthreads=1 --- src/marching_iterators.jl | 26 +++++++++++++++++--------- 1 file changed, 17 insertions(+), 9 deletions(-) diff --git a/src/marching_iterators.jl b/src/marching_iterators.jl index 338fc44..a7b67b6 100644 --- a/src/marching_iterators.jl +++ b/src/marching_iterators.jl @@ -1,4 +1,4 @@ -function pushintersection!(intersection_points,triangle,levels) +function pushintersection!(intersection_points,triangle::LinearSimplex{2},levels) f = values(triangle) coord=points(triangle) @@ -41,19 +41,27 @@ function pushintersection!(intersection_points,triangle,levels) end end +function intersections(triangles::T,levels, Tp) where T<:LinearSimplexIterator{2} + local intersection_points = Vector{Tp}(undef, 0) + for triangle in triangles + pushintersection!(intersection_points,triangle) + end + intersection_points +end function marching_triangles(triangle_iterators::Vector{LinearSimplexIterator{2}}, levels; Tc = Float32, Tp = SVector{2, Tc}) - threads=map(triangle_iterators) do triangles - Threads.@spawn begin - local intersection_points = Vector{Tp}(undef, 0) - for triangle in triangles - pushintersection!(intersection_points,triangle) - end - intersection_points + + if Threads.nthreads()==1 + map(triangle_iterators) do triangles + intersections(triangles,levels,Tp) + end + else + threads=map(triangle_iterators) do triangles + Threads.@spawn intersections(triangles,levels,Tp) end + fetch.(threads) end - fetch.(threads) end From e0acfcee9ca515625ff89efdd8292d550910b6ee Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Wed, 24 Jan 2024 23:19:17 +0100 Subject: [PATCH 3/9] first working plot --- src/linearsimplex.jl | 22 +++++++++++++++++++++- src/marching_iterators.jl | 20 ++++++++------------ 2 files changed, 29 insertions(+), 13 deletions(-) diff --git a/src/linearsimplex.jl b/src/linearsimplex.jl index 52ba6e0..d533972 100644 --- a/src/linearsimplex.jl +++ b/src/linearsimplex.jl @@ -31,6 +31,25 @@ function LinearSimplex(::Type{Val{3}},points::AbstractMatrix, values::AbstractVe LinearSimplex(spoints,svalues) end + +function LinearSimplex(::Type{Val{1}},points::AbstractMatrix, values::AbstractVector,coordscale) + @views spoints=SVector{2,Point{1,Float32}}(points[:,1]*coordscale,points[:,2]*coordscale) + svalues=SVector{2,Float32}(values[1],values[2]) + LinearSimplex(spoints,svalues) +end + +function LinearSimplex(::Type{Val{2}},points::AbstractMatrix, values::AbstractVector,coordscale) + @views spoints=SVector{3,Point{2,Float32}}(points[:,1]*coordscale,points[:,2]*coordscale,points[:,3]*coordscale) + svalues=SVector{3,Float32}(values[1],values[2],values[3]) + LinearSimplex(spoints,svalues) +end + +function LinearSimplex(::Type{Val{3}},points::AbstractMatrix, values::AbstractVector,coordscale) + @views spoints=SVector{4,Point{3,Float32}}(points[:,1]*coordscale,points[:,2]*coordscale,points[:,3]*coordscale,points[:,4]*coordscale) + svalues=SVector{4,Float32}(values[1],values[2],values[3],values[4]) + LinearSimplex(spoints,svalues) +end + LinearEdge(points,values)=LinearSimplex(Val{1},points,values) LinearTriangle(points,values)=LinearSimplex(Val{2},points,values) LinearTetrahedron(points,values)=LinearSimplex(Val{3},points,values) @@ -58,7 +77,8 @@ Sum up all values passed via the iterator. Designed for testing iterators. Useful for: - Consistency test - e.g. when passing constant ones, the sum can be calculated by other means and compared to the return value of testloop -- Allocation test. `testloop` itself does not allocate, so any allocations from a +- Allocation test. `testloop` itself does allocate only a few (<1000) bytes in + very few (<10) allocations. So any allocations beyond this from a call to `testloop` hint at possibilities to improve an iterator implementation. """ function testloop(iterators::AbstractVector{T}) where T<:LinearSimplexIterator diff --git a/src/marching_iterators.jl b/src/marching_iterators.jl index a7b67b6..1a85305 100644 --- a/src/marching_iterators.jl +++ b/src/marching_iterators.jl @@ -35,32 +35,28 @@ function pushintersection!(intersection_points,triangle::LinearSimplex{2},levels x2 = coord[n2][1] + α * dx32 y2 = coord[n2][2] + α * dy32 end - push!(intersection_points, SVector{2, Tc}((x1, y1))) - push!(intersection_points, SVector{2, Tc}((x2, y2))) + push!(intersection_points, Point{2,Float32}(x1, y1)) + push!(intersection_points, Point{2,Float32}(x2, y2)) end end end -function intersections(triangles::T,levels, Tp) where T<:LinearSimplexIterator{2} - local intersection_points = Vector{Tp}(undef, 0) +function intersections(triangles::T,levels) where T<:LinearSimplexIterator{2} + local intersection_points = Vector{Point{2,Float32}}(undef, 0) for triangle in triangles - pushintersection!(intersection_points,triangle) + pushintersection!(intersection_points,triangle,levels) end intersection_points end -function marching_triangles(triangle_iterators::Vector{LinearSimplexIterator{2}}, - levels; - Tc = Float32, - Tp = SVector{2, Tc}) - +function marching_triangles(triangle_iterators::Vector{T},levels) where T<:LinearSimplexIterator{2} if Threads.nthreads()==1 map(triangle_iterators) do triangles - intersections(triangles,levels,Tp) + intersections(triangles,levels) end else threads=map(triangle_iterators) do triangles - Threads.@spawn intersections(triangles,levels,Tp) + Threads.@spawn intersections(triangles,levels) end fetch.(threads) end From 22ea4c50b8728521be723a0f62cd29f7ed018d1f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Tue, 17 Jun 2025 12:45:52 +0200 Subject: [PATCH 4/9] runic in iterator-api branch --- docs/make.jl | 22 +++-- src/colors.jl | 26 +++--- src/extraction.jl | 49 +++++----- src/linearsimplex.jl | 88 +++++++++--------- src/marching.jl | 189 +++++++++++++++++++++----------------- src/marching_iterators.jl | 43 ++++----- src/markerpoints.jl | 4 +- src/planeslevels.jl | 14 +-- 8 files changed, 234 insertions(+), 201 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index 06b6d6b..ae4e542 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -2,16 +2,18 @@ using Documenter, GridVisualizeTools, ColorTypes function mkdocs() DocMeta.setdocmeta!(GridVisualizeTools, :DocTestSetup, :(using GridVisualizeTools, ColorTypes, Colors); recursive = true) - makedocs(; sitename = "GridVisualizeTools.jl", - modules = [GridVisualizeTools], - clean = false, - doctest = true, - authors = "J. Fuhrmann", - repo = "https://github.com/j-fu/GridVisualizeTools.jl", - pages = [ - "Home" => "index.md", - ]) - if !isinteractive() + makedocs(; + sitename = "GridVisualizeTools.jl", + modules = [GridVisualizeTools], + clean = false, + doctest = true, + authors = "J. Fuhrmann", + repo = "https://github.com/j-fu/GridVisualizeTools.jl", + pages = [ + "Home" => "index.md", + ] + ) + return if !isinteractive() deploydocs(; repo = "github.com/j-fu/GridVisualizeTools.jl.git", devbranch = "main") end end diff --git a/src/colors.jl b/src/colors.jl index 15c5bca..21657ab 100644 --- a/src/colors.jl +++ b/src/colors.jl @@ -15,11 +15,13 @@ RGB{Float64}(0.85,0.6,0.6) """ function region_cmap(n) - ColorSchemes.distinguishable_colors(max(5, n), - [Colors.RGB(0.85, 0.6, 0.6), Colors.RGB(0.6, 0.85, 0.6), Colors.RGB(0.6, 0.6, 0.85)]; - lchoices = range(70; stop = 80, length = 5), - cchoices = range(25; stop = 65, length = 15), - hchoices = range(20; stop = 360, length = 15)) + return ColorSchemes.distinguishable_colors( + max(5, n), + [Colors.RGB(0.85, 0.6, 0.6), Colors.RGB(0.6, 0.85, 0.6), Colors.RGB(0.6, 0.6, 0.85)]; + lchoices = range(70; stop = 80, length = 5), + cchoices = range(25; stop = 65, length = 15), + hchoices = range(20; stop = 360, length = 15) + ) end """ @@ -38,11 +40,13 @@ RGB{Float64}(1.0,0.0,0.0) """ function bregion_cmap(n) - ColorSchemes.distinguishable_colors(max(5, n), - [Colors.RGB(1.0, 0.0, 0.0), Colors.RGB(0.0, 1.0, 0.0), Colors.RGB(0.0, 0.0, 1.0)]; - lchoices = range(50; stop = 75, length = 10), - cchoices = range(75; stop = 100, length = 10), - hchoices = range(20; stop = 360, length = 30)) + return ColorSchemes.distinguishable_colors( + max(5, n), + [Colors.RGB(1.0, 0.0, 0.0), Colors.RGB(0.0, 1.0, 0.0), Colors.RGB(0.0, 0.0, 1.0)]; + lchoices = range(50; stop = 75, length = 10), + cchoices = range(75; stop = 100, length = 10), + hchoices = range(20; stop = 360, length = 30) + ) end """ @@ -56,7 +60,7 @@ RGB{Float64}(1.0,0.0,0.0) """ function Colors.RGB(c::String) c64 = Colors.color_names[c] - Colors.RGB(c64[1] / 255, c64[2] / 255, c64[3] / 255) + return Colors.RGB(c64[1] / 255, c64[2] / 255, c64[3] / 255) end """ diff --git a/src/extraction.jl b/src/extraction.jl index 4601ff1..e60c331 100644 --- a/src/extraction.jl +++ b/src/extraction.jl @@ -7,16 +7,18 @@ Extract visible tetrahedra - those intersecting with the planes Return corresponding points and facets for each region for drawing as mesh (Makie,MeshCat) or trisurf (pyplot) """ -function extract_visible_cells3D(coord, cellnodes, cellregions, nregions, xyzcut; - primepoints = zeros(0, 0), Tp = SVector{3, Float32}, Tf = SVector{3, Int32}) +function extract_visible_cells3D( + coord, cellnodes, cellregions, nregions, xyzcut; + primepoints = zeros(0, 0), Tp = SVector{3, Float32}, Tf = SVector{3, Int32} + ) all_lt = ones(Bool, 3) all_gt = ones(Bool, 3) function take(coord, simplex, xyzcut, all_lt, all_gt) - for idim = 1:3 + for idim in 1:3 all_lt[idim] = true all_gt[idim] = true - for inode = 1:4 + for inode in 1:4 c = coord[idim, simplex[inode]] - xyzcut[idim] all_lt[idim] = all_lt[idim] && (c < 0.0) all_gt[idim] = all_gt[idim] && (c > 0.0) @@ -25,22 +27,22 @@ function extract_visible_cells3D(coord, cellnodes, cellregions, nregions, xyzcut tke = false tke = tke || (!all_lt[1]) && (!all_gt[1]) && (!all_gt[2]) && (!all_gt[3]) tke = tke || (!all_lt[2]) && (!all_gt[2]) && (!all_gt[1]) && (!all_gt[3]) - tke = tke || (!all_lt[3]) && (!all_gt[3]) && (!all_gt[1]) && (!all_gt[2]) + return tke = tke || (!all_lt[3]) && (!all_gt[3]) && (!all_gt[1]) && (!all_gt[2]) end - faces = [Vector{Tf}(undef, 0) for iregion = 1:nregions] - points = [Vector{Tp}(undef, 0) for iregion = 1:nregions] + faces = [Vector{Tf}(undef, 0) for iregion in 1:nregions] + points = [Vector{Tp}(undef, 0) for iregion in 1:nregions] - for iregion = 1:nregions - for iprime = 1:size(primepoints, 2) + for iregion in 1:nregions + for iprime in 1:size(primepoints, 2) @views push!(points[iregion], Tp(primepoints[:, iprime])) end end tet = zeros(Int32, 4) - for itet = 1:size(cellnodes, 2) + for itet in 1:size(cellnodes, 2) iregion = cellregions[itet] - for i = 1:4 + for i in 1:4 tet[i] = cellnodes[i, itet] end if take(coord, tet, xyzcut, all_lt, all_gt) @@ -57,7 +59,7 @@ function extract_visible_cells3D(coord, cellnodes, cellregions, nregions, xyzcut end end end - points, faces + return points, faces end """ @@ -69,15 +71,17 @@ Extract visible boundary faces - those not cut off by the planes Return corresponding points and facets for each region for drawing as mesh (Makie,MeshCat) or trisurf (pyplot) """ -function extract_visible_bfaces3D(coord, bfacenodes, bfaceregions, nbregions, xyzcut; - primepoints = zeros(0, 0), Tp = SVector{3, Float32}, Tf = SVector{3, Int32}) +function extract_visible_bfaces3D( + coord, bfacenodes, bfaceregions, nbregions, xyzcut; + primepoints = zeros(0, 0), Tp = SVector{3, Float32}, Tf = SVector{3, Int32} + ) nbfaces = size(bfacenodes, 2) cutcoord = zeros(3) function take(coord, simplex, xyzcut) - for idim = 1:3 + for idim in 1:3 all_gt = true - for inode = 1:3 + for inode in 1:3 c = coord[idim, simplex[inode]] - xyzcut[idim] all_gt = all_gt && c > 0 end @@ -91,10 +95,10 @@ function extract_visible_bfaces3D(coord, bfacenodes, bfaceregions, nbregions, xy Tc = SVector{3, eltype(coord)} xcoord = reinterpret(Tc, reshape(coord, (length(coord),))) - faces = [Vector{Tf}(undef, 0) for iregion = 1:nbregions] - points = [Vector{Tp}(undef, 0) for iregion = 1:nbregions] - for iregion = 1:nbregions - for iprime = 1:size(primepoints, 2) + faces = [Vector{Tf}(undef, 0) for iregion in 1:nbregions] + points = [Vector{Tp}(undef, 0) for iregion in 1:nbregions] + for iregion in 1:nbregions + for iprime in 1:size(primepoints, 2) @views push!(points[iregion], Tp(primepoints[:, iprime])) end end @@ -102,7 +106,7 @@ function extract_visible_bfaces3D(coord, bfacenodes, bfaceregions, nbregions, xy # remove some type instability here function collct(points, faces) trinodes = [1, 2, 3] - for i = 1:nbfaces + for i in 1:nbfaces iregion = bfaceregions[i] trinodes[1] = bfacenodes[1, i] trinodes[2] = bfacenodes[2, i] @@ -115,7 +119,8 @@ function extract_visible_bfaces3D(coord, bfacenodes, bfaceregions, nbregions, xy @views push!(faces[iregion], (npts + 1, npts + 2, npts + 3)) end end + return end collct(points, faces) - points, faces + return points, faces end diff --git a/src/linearsimplex.jl b/src/linearsimplex.jl index d533972..11e8c23 100644 --- a/src/linearsimplex.jl +++ b/src/linearsimplex.jl @@ -1,58 +1,58 @@ -struct LinearSimplex{D,N,Tv} - points::SVector{N,Point{D,Tv}} - values::SVector{N,Tv} +struct LinearSimplex{D, N, Tv} + points::SVector{N, Point{D, Tv}} + values::SVector{N, Tv} end -values(s::LinearSimplex)=s.values -points(s::LinearSimplex)=s.points +values(s::LinearSimplex) = s.values +points(s::LinearSimplex) = s.points -function LinearSimplex(::Type{Val{D}},points::Union{Vector,Tuple}, values::Union{Vector,Tuple}) where {D} - spoints=SVector{D+1,Point{D,Float32}}(points...) - svalues=SVector{D+1,Float32}(values...) - LinearSimplex(spoints,svalues) +function LinearSimplex(::Type{Val{D}}, points::Union{Vector, Tuple}, values::Union{Vector, Tuple}) where {D} + spoints = SVector{D + 1, Point{D, Float32}}(points...) + svalues = SVector{D + 1, Float32}(values...) + return LinearSimplex(spoints, svalues) end -function LinearSimplex(::Type{Val{1}},points::AbstractMatrix, values::AbstractVector) - @views spoints=SVector{2,Point{1,Float32}}(points[:,1],points[:,2]) - svalues=SVector{2,Float32}(values[1],values[2]) - LinearSimplex(spoints,svalues) +function LinearSimplex(::Type{Val{1}}, points::AbstractMatrix, values::AbstractVector) + @views spoints = SVector{2, Point{1, Float32}}(points[:, 1], points[:, 2]) + svalues = SVector{2, Float32}(values[1], values[2]) + return LinearSimplex(spoints, svalues) end -function LinearSimplex(::Type{Val{2}},points::AbstractMatrix, values::AbstractVector) - @views spoints=SVector{3,Point{2,Float32}}(points[:,1],points[:,2],points[:,3]) - svalues=SVector{3,Float32}(values[1],values[2],values[3]) - LinearSimplex(spoints,svalues) +function LinearSimplex(::Type{Val{2}}, points::AbstractMatrix, values::AbstractVector) + @views spoints = SVector{3, Point{2, Float32}}(points[:, 1], points[:, 2], points[:, 3]) + svalues = SVector{3, Float32}(values[1], values[2], values[3]) + return LinearSimplex(spoints, svalues) end -function LinearSimplex(::Type{Val{3}},points::AbstractMatrix, values::AbstractVector) - @views spoints=SVector{4,Point{3,Float32}}(points[:,1],points[:,2],points[:,3],points[:,4]) - svalues=SVector{4,Float32}(values[1],values[2],values[3],values[4]) - LinearSimplex(spoints,svalues) +function LinearSimplex(::Type{Val{3}}, points::AbstractMatrix, values::AbstractVector) + @views spoints = SVector{4, Point{3, Float32}}(points[:, 1], points[:, 2], points[:, 3], points[:, 4]) + svalues = SVector{4, Float32}(values[1], values[2], values[3], values[4]) + return LinearSimplex(spoints, svalues) end -function LinearSimplex(::Type{Val{1}},points::AbstractMatrix, values::AbstractVector,coordscale) - @views spoints=SVector{2,Point{1,Float32}}(points[:,1]*coordscale,points[:,2]*coordscale) - svalues=SVector{2,Float32}(values[1],values[2]) - LinearSimplex(spoints,svalues) +function LinearSimplex(::Type{Val{1}}, points::AbstractMatrix, values::AbstractVector, coordscale) + @views spoints = SVector{2, Point{1, Float32}}(points[:, 1] * coordscale, points[:, 2] * coordscale) + svalues = SVector{2, Float32}(values[1], values[2]) + return LinearSimplex(spoints, svalues) end -function LinearSimplex(::Type{Val{2}},points::AbstractMatrix, values::AbstractVector,coordscale) - @views spoints=SVector{3,Point{2,Float32}}(points[:,1]*coordscale,points[:,2]*coordscale,points[:,3]*coordscale) - svalues=SVector{3,Float32}(values[1],values[2],values[3]) - LinearSimplex(spoints,svalues) +function LinearSimplex(::Type{Val{2}}, points::AbstractMatrix, values::AbstractVector, coordscale) + @views spoints = SVector{3, Point{2, Float32}}(points[:, 1] * coordscale, points[:, 2] * coordscale, points[:, 3] * coordscale) + svalues = SVector{3, Float32}(values[1], values[2], values[3]) + return LinearSimplex(spoints, svalues) end -function LinearSimplex(::Type{Val{3}},points::AbstractMatrix, values::AbstractVector,coordscale) - @views spoints=SVector{4,Point{3,Float32}}(points[:,1]*coordscale,points[:,2]*coordscale,points[:,3]*coordscale,points[:,4]*coordscale) - svalues=SVector{4,Float32}(values[1],values[2],values[3],values[4]) - LinearSimplex(spoints,svalues) +function LinearSimplex(::Type{Val{3}}, points::AbstractMatrix, values::AbstractVector, coordscale) + @views spoints = SVector{4, Point{3, Float32}}(points[:, 1] * coordscale, points[:, 2] * coordscale, points[:, 3] * coordscale, points[:, 4] * coordscale) + svalues = SVector{4, Float32}(values[1], values[2], values[3], values[4]) + return LinearSimplex(spoints, svalues) end -LinearEdge(points,values)=LinearSimplex(Val{1},points,values) -LinearTriangle(points,values)=LinearSimplex(Val{2},points,values) -LinearTetrahedron(points,values)=LinearSimplex(Val{3},points,values) +LinearEdge(points, values) = LinearSimplex(Val{1}, points, values) +LinearTriangle(points, values) = LinearSimplex(Val{2}, points, values) +LinearTetrahedron(points, values) = LinearSimplex(Val{3}, points, values) """ @@ -81,15 +81,15 @@ Useful for: very few (<10) allocations. So any allocations beyond this from a call to `testloop` hint at possibilities to improve an iterator implementation. """ -function testloop(iterators::AbstractVector{T}) where T<:LinearSimplexIterator - threads=map(iterators) do iterator - begin - local x=0.0 +function testloop(iterators::AbstractVector{T}) where {T <: LinearSimplexIterator} + threads = map(iterators) do iterator + begin + local x = 0.0 for vt in iterator - x+=sum(vt.values) + x += sum(vt.values) end - x - end + x + end end - sum(fetch.(threads)) + return sum(fetch.(threads)) end diff --git a/src/marching.jl b/src/marching.jl index b559858..6709160 100644 --- a/src/marching.jl +++ b/src/marching.jl @@ -25,28 +25,32 @@ This method can be used both for the evaluation of plane sections and for the evaluation of function isosurfaces. """ -function tet_x_plane!(ixcoord, - ixvalues, - pointlist, - node_indices, - planeq_values, - function_values; - tol = 0.0) +function tet_x_plane!( + ixcoord, + ixvalues, + pointlist, + node_indices, + planeq_values, + function_values; + tol = 0.0 + ) # If all nodes lie on one side of the plane, no intersection - @fastmath if (mapreduce(a -> a < -tol, *, planeq_values) || - mapreduce(a -> a > tol, *, planeq_values)) + @fastmath if ( + mapreduce(a -> a < -tol, *, planeq_values) || + mapreduce(a -> a > tol, *, planeq_values) + ) return 0 end # Interpolate coordinates and function_values according to # evaluation of the plane equation nxs = 0 - @inbounds @simd for n1 = 1:3 + @inbounds @simd for n1 in 1:3 N1 = node_indices[n1] - @inbounds @fastmath @simd for n2 = (n1 + 1):4 + @inbounds @fastmath @simd for n2 in (n1 + 1):4 N2 = node_indices[n2] if planeq_values[n1] != planeq_values[n2] && - planeq_values[n1] * planeq_values[n2] < tol + planeq_values[n1] * planeq_values[n2] < tol nxs += 1 t = planeq_values[n1] / (planeq_values[n1] - planeq_values[n2]) ixcoord[1, nxs] = pointlist[1, N1] + t * (pointlist[1, N2] - pointlist[1, N1]) @@ -117,41 +121,47 @@ These can be readily turned into a mesh with function values on it. Caveat: points with similar coordinates are not identified, e.g. an intersection of a plane and an edge will generate as many edge intersection points as there are tetrahedra adjacent to that edge. As a consequence, normal calculations for visualization alway will end up with facet normals, not point normals, and the visual impression of a rendered isosurface will show its piecewise linear genealogy. """ -function marching_tetrahedra(coord::Matrix{Tc}, - cellnodes::Matrix{Ti}, - func, - planes, - flevels; - tol = 1.0e-12, - primepoints = zeros(0, 0), - primevalues = zeros(0), - Tv = Float32, - Tp = SVector{3, Float32}, - Tf = SVector{3, Int32}) where {Tc, Ti} - marching_tetrahedra([coord], - [cellnodes], - [func], - planes, - flevels; - tol, - primepoints, - primevalues, - Tv, - Tp, - Tf) +function marching_tetrahedra( + coord::Matrix{Tc}, + cellnodes::Matrix{Ti}, + func, + planes, + flevels; + tol = 1.0e-12, + primepoints = zeros(0, 0), + primevalues = zeros(0), + Tv = Float32, + Tp = SVector{3, Float32}, + Tf = SVector{3, Int32} + ) where {Tc, Ti} + return marching_tetrahedra( + [coord], + [cellnodes], + [func], + planes, + flevels; + tol, + primepoints, + primevalues, + Tv, + Tp, + Tf + ) end -function marching_tetrahedra(allcoords::Vector{Matrix{Tc}}, - allcellnodes::Vector{Matrix{Ti}}, - allfuncs, - planes, - flevels; - tol = 1.0e-12, - primepoints = zeros(0, 0), - primevalues = zeros(0), - Tv = Float32, - Tp = SVector{3, Float32}, - Tf = SVector{3, Int32}) where {Tc, Ti} +function marching_tetrahedra( + allcoords::Vector{Matrix{Tc}}, + allcellnodes::Vector{Matrix{Ti}}, + allfuncs, + planes, + flevels; + tol = 1.0e-12, + primepoints = zeros(0, 0), + primevalues = zeros(0), + Tv = Float32, + Tp = SVector{3, Float32}, + Tf = SVector{3, Int32} + ) where {Tc, Ti} # We could rewrite this for Meshing.jl # CellNodes::Vector{Ttet}, Coord::Vector{Tpt} @@ -163,8 +173,8 @@ function marching_tetrahedra(allcoords::Vector{Matrix{Tc}}, all_ixcoord = Vector{Tp}(undef, 0) all_ixvalues = Vector{Tv}(undef, 0) - @assert(length(primevalues)==size(primepoints, 2)) - for iprime = 1:size(primepoints, 2) + @assert(length(primevalues) == size(primepoints, 2)) + for iprime in 1:size(primepoints, 2) @views push!(all_ixcoord, primepoints[:, iprime]) @views push!(all_ixvalues, primevalues[iprime]) end @@ -180,9 +190,9 @@ function marching_tetrahedra(allcoords::Vector{Matrix{Tc}}, function pushtris(ns, ixcoord, ixvalues) # number of intersection points can be 3 or 4 - if ns >= 3 + return if ns >= 3 last_i = length(all_ixvalues) - for is = 1:ns + for is in 1:ns @views push!(all_ixcoord, ixcoord[:, is]) push!(all_ixvalues, ixvalues[is]) # todo consider nan_replacement here end @@ -193,7 +203,7 @@ function marching_tetrahedra(allcoords::Vector{Matrix{Tc}}, end end - for igrid = 1:length(allcoords) + for igrid in 1:length(allcoords) coord = allcoords[igrid] cellnodes = allcellnodes[igrid] func = allfuncs[igrid] @@ -202,7 +212,7 @@ function marching_tetrahedra(allcoords::Vector{Matrix{Tc}}, all_planeq = Vector{Float32}(undef, nnodes) function calcxs() - @inbounds for itet = 1:ntet + return @inbounds for itet in 1:ntet node_indices[1] = cellnodes[1, itet] node_indices[2] = cellnodes[2, itet] node_indices[3] = cellnodes[3, itet] @@ -211,33 +221,39 @@ function marching_tetrahedra(allcoords::Vector{Matrix{Tc}}, planeq[2] = all_planeq[node_indices[2]] planeq[3] = all_planeq[node_indices[3]] planeq[4] = all_planeq[node_indices[4]] - nxs = tet_x_plane!(ixcoord, - ixvalues, - coord, - node_indices, - planeq, - func; - tol = tol) + nxs = tet_x_plane!( + ixcoord, + ixvalues, + coord, + node_indices, + planeq, + func; + tol = tol + ) pushtris(nxs, ixcoord, ixvalues) end end - @inbounds for iplane = 1:nplanes - @views @inbounds map!(inode -> plane_equation(planes[iplane], coord[:, inode]), - all_planeq, - 1:nnodes) + @inbounds for iplane in 1:nplanes + @views @inbounds map!( + inode -> plane_equation(planes[iplane], coord[:, inode]), + all_planeq, + 1:nnodes + ) calcxs() end # allocation free (besides push!) - @inbounds for ilevel = 1:nlevels - @views @inbounds @fastmath map!(inode -> (func[inode] - flevels[ilevel]), - all_planeq, - 1:nnodes) + @inbounds for ilevel in 1:nlevels + @views @inbounds @fastmath map!( + inode -> (func[inode] - flevels[ilevel]), + all_planeq, + 1:nnodes + ) calcxs() end end - all_ixcoord, all_ixfaces, all_ixvalues + return all_ixcoord, all_ixfaces, all_ixvalues end """ @@ -245,24 +261,28 @@ end Collect isoline snippets on triangles ready for linesegments! """ -function marching_triangles(coord::Matrix{Tv}, - cellnodes::Matrix{Ti}, - func, - levels; - Tc = Float32, - Tp = SVector{2, Tc}) where {Tv <: Number, Ti <: Number} - marching_triangles([coord], [cellnodes], [func], levels; Tc, Tp) +function marching_triangles( + coord::Matrix{Tv}, + cellnodes::Matrix{Ti}, + func, + levels; + Tc = Float32, + Tp = SVector{2, Tc} + ) where {Tv <: Number, Ti <: Number} + return marching_triangles([coord], [cellnodes], [func], levels; Tc, Tp) end -function marching_triangles(coords::Vector{Matrix{Tv}}, - cellnodes::Vector{Matrix{Ti}}, - funcs, - levels; - Tc = Float32, - Tp = SVector{2, Tc}) where {Tv <: Number, Ti <: Number} +function marching_triangles( + coords::Vector{Matrix{Tv}}, + cellnodes::Vector{Matrix{Ti}}, + funcs, + levels; + Tc = Float32, + Tp = SVector{2, Tc} + ) where {Tv <: Number, Ti <: Number} points = Vector{Tp}(undef, 0) - for igrid = 1:length(coords) + for igrid in 1:length(coords) func = funcs[igrid] coord = coords[igrid] @@ -289,7 +309,7 @@ function marching_triangles(coords::Vector{Matrix{Tv}}, df21 = f[i2] != f[i1] ? 1 / (f[i2] - f[i1]) : 0.0 df32 = f[i3] != f[i2] ? 1 / (f[i3] - f[i2]) : 0.0 - for level ∈ levels + for level in levels if (f[i1] <= level) && (level < f[i3]) α = (level - f[i1]) * df31 x1 = coord[1, n1] + α * dx31 @@ -308,11 +328,12 @@ function marching_triangles(coords::Vector{Matrix{Tv}}, push!(points, SVector{2, Tc}((x2, y2))) end end + return end - for itri = 1:size(cellnodes[igrid], 2) + for itri in 1:size(cellnodes[igrid], 2) @views isect(cellnodes[igrid][:, itri]) end end - points + return points end diff --git a/src/marching_iterators.jl b/src/marching_iterators.jl index 1a85305..b5f68e4 100644 --- a/src/marching_iterators.jl +++ b/src/marching_iterators.jl @@ -1,31 +1,31 @@ -function pushintersection!(intersection_points,triangle::LinearSimplex{2},levels) +function pushintersection!(intersection_points, triangle::LinearSimplex{2}, levels) f = values(triangle) - coord=points(triangle) - + coord = points(triangle) + (n1, n2, n3) = (1, 2, 3) - + f[1] <= f[2] ? (n1, n2) = (1, 2) : (n1, n2) = (2, 1) f[n2] <= f[3] ? n3 = 3 : (n2, n3) = (3, n2) f[n1] > f[n2] ? (n1, n2) = (n2, n1) : nothing - + dx31 = coord[n3][1] - coord[n1][1] dx21 = coord[n2][1] - coord[n1][1] dx32 = coord[n3][1] - coord[n2][1] - + dy31 = coord[n3][2] - coord[n1][2] dy21 = coord[n2][2] - coord[n1][2] dy32 = coord[n3][2] - coord[n2][2] - + df31 = f[n3] != f[n1] ? 1 / (f[n3] - f[n1]) : 0.0 df21 = f[n2] != f[n1] ? 1 / (f[n2] - f[n1]) : 0.0 df32 = f[n3] != f[n2] ? 1 / (f[n3] - f[n2]) : 0.0 - - for level ∈ levels + + for level in levels if (f[n1] <= level) && (level < f[n3]) α = (level - f[n1]) * df31 x1 = coord[n1][1] + α * dx31 y1 = coord[n1][2] + α * dy31 - + if (level < f[n2]) α = (level - f[n1]) * df21 x2 = coord[n1][1] + α * dx21 @@ -35,28 +35,29 @@ function pushintersection!(intersection_points,triangle::LinearSimplex{2},levels x2 = coord[n2][1] + α * dx32 y2 = coord[n2][2] + α * dy32 end - push!(intersection_points, Point{2,Float32}(x1, y1)) - push!(intersection_points, Point{2,Float32}(x2, y2)) + push!(intersection_points, Point{2, Float32}(x1, y1)) + push!(intersection_points, Point{2, Float32}(x2, y2)) end end + return end -function intersections(triangles::T,levels) where T<:LinearSimplexIterator{2} - local intersection_points = Vector{Point{2,Float32}}(undef, 0) +function intersections(triangles::T, levels) where {T <: LinearSimplexIterator{2}} + local intersection_points = Vector{Point{2, Float32}}(undef, 0) for triangle in triangles - pushintersection!(intersection_points,triangle,levels) + pushintersection!(intersection_points, triangle, levels) end - intersection_points + return intersection_points end -function marching_triangles(triangle_iterators::Vector{T},levels) where T<:LinearSimplexIterator{2} - if Threads.nthreads()==1 +function marching_triangles(triangle_iterators::Vector{T}, levels) where {T <: LinearSimplexIterator{2}} + return if Threads.nthreads() == 1 map(triangle_iterators) do triangles - intersections(triangles,levels) + intersections(triangles, levels) end else - threads=map(triangle_iterators) do triangles - Threads.@spawn intersections(triangles,levels) + threads = map(triangle_iterators) do triangles + Threads.@spawn intersections(triangles, levels) end fetch.(threads) end diff --git a/src/markerpoints.jl b/src/markerpoints.jl index 2ebfde7..0f94dff 100644 --- a/src/markerpoints.jl +++ b/src/markerpoints.jl @@ -11,7 +11,7 @@ function markerpoints(points, nmarkers, transform) dist(p1, p2) = norm(transform * (p1 - p2)) llen = 0.0 - for i = 2:length(points) + for i in 2:length(points) llen += dist(points[i], points[i - 1]) end @@ -36,5 +36,5 @@ function markerpoints(points, nmarkers, transform) lnext = lnext + mdist end end - push!(mpoints, points[end]) + return push!(mpoints, points[end]) end diff --git a/src/planeslevels.jl b/src/planeslevels.jl index 4f1a7fb..a4cfb7d 100644 --- a/src/planeslevels.jl +++ b/src/planeslevels.jl @@ -1,5 +1,5 @@ function makeplanes(mmin, mmax, n) - if isa(n, Number) + return if isa(n, Number) if n == 0 return [Inf] end @@ -24,21 +24,21 @@ function makeplanes(xyzmin, xyzmax, x, y, z) Y = makeplanes(xyzmin[2], xyzmax[2], y) Z = makeplanes(xyzmin[3], xyzmax[3], z) - for i = 1:length(X) + for i in 1:length(X) x = X[i] x > xyzmin[1] && x < xyzmax[1] && push!(planes, [1, 0, 0, -x]) end - for i = 1:length(Y) + for i in 1:length(Y) y = Y[i] y > xyzmin[2] && y < xyzmax[2] && push!(planes, [0, 1, 0, -y]) end - for i = 1:length(Z) + for i in 1:length(Z) z = Z[i] z > xyzmin[3] && z < xyzmax[3] && push!(planes, [0, 0, 1, -z]) end - planes + return planes end """ @@ -52,7 +52,7 @@ Update levels, limits, colorbartics based on vector given in func. otherwise, if it is a number, replace it with a linear range of corresponding length """ function makeisolevels(func::Vector{T}, levels, limits, colorbarticks) where {T <: Number} - makeisolevels([func], levels, limits, colorbarticks) + return makeisolevels([func], levels, limits, colorbarticks) end function makeisolevels(funcs::Vector{T}, levels, limits, colorbarticks) where {T <: AbstractVector} @@ -72,5 +72,5 @@ function makeisolevels(funcs::Vector{T}, levels, limits, colorbarticks) where {T end # map(t->round(t,sigdigits=4),levels),limits,map(t->round(t,sigdigits=4),colorbarticks) - levels, limits, colorbarticks + return levels, limits, colorbarticks end From e88960ed05173036f17cb2a6be0157f4a3e44625 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Tue, 17 Jun 2025 14:16:15 +0200 Subject: [PATCH 5/9] Re-introduce GeometryBasics for Point Replace this later by SArray. --- Project.toml | 2 ++ src/GridVisualizeTools.jl | 1 + 2 files changed, 3 insertions(+) diff --git a/Project.toml b/Project.toml index dc15685..bee2040 100644 --- a/Project.toml +++ b/Project.toml @@ -7,6 +7,7 @@ version = "3.0.2" ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4" Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae" +GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" @@ -14,6 +15,7 @@ StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c" ColorSchemes = "3" Colors = "0.12, 0.13" DocStringExtensions = "0.8, 0.9" +GeometryBasics = "0.5.9" StaticArrays = "1.9.13" StaticArraysCore = "1.4" julia = "1.6" diff --git a/src/GridVisualizeTools.jl b/src/GridVisualizeTools.jl index 58d9a7b..d6545a7 100644 --- a/src/GridVisualizeTools.jl +++ b/src/GridVisualizeTools.jl @@ -11,6 +11,7 @@ import ColorSchemes using DocStringExtensions: SIGNATURES, TYPEDEF, TYPEDSIGNATURES using StaticArraysCore: SVector using StaticArrays: @MArray +using GeometryBasics # replace by SArray! include("colors.jl") export region_cmap, bregion_cmap, rgbtuple, rgbcolor From 8054a5fe22f900f7ce3f9af45835cda3d0fb8786 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Tue, 17 Jun 2025 14:17:05 +0200 Subject: [PATCH 6/9] runic --- src/GridVisualizeTools.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/GridVisualizeTools.jl b/src/GridVisualizeTools.jl index d6545a7..a3fe542 100644 --- a/src/GridVisualizeTools.jl +++ b/src/GridVisualizeTools.jl @@ -11,7 +11,7 @@ import ColorSchemes using DocStringExtensions: SIGNATURES, TYPEDEF, TYPEDSIGNATURES using StaticArraysCore: SVector using StaticArrays: @MArray -using GeometryBasics # replace by SArray! +using GeometryBasics # replace by SArray! include("colors.jl") export region_cmap, bregion_cmap, rgbtuple, rgbcolor From b946834dd34f34e5ebf907efc1da34aa57c680eb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Fuhrmann?= Date: Tue, 17 Jun 2025 14:25:32 +0200 Subject: [PATCH 7/9] address copilot's remarks --- src/linearsimplex.jl | 4 ++++ src/marching.jl | 2 +- 2 files changed, 5 insertions(+), 1 deletion(-) diff --git a/src/linearsimplex.jl b/src/linearsimplex.jl index 11e8c23..783b6c7 100644 --- a/src/linearsimplex.jl +++ b/src/linearsimplex.jl @@ -1,3 +1,7 @@ +""" + LinearSimplex +Linear simplex with point coordinates and node values. +""" struct LinearSimplex{D, N, Tv} points::SVector{N, Point{D, Tv}} values::SVector{N, Tv} diff --git a/src/marching.jl b/src/marching.jl index 2111c52..7acc106 100644 --- a/src/marching.jl +++ b/src/marching.jl @@ -61,7 +61,7 @@ end This method can be used both for the evaluation of plane sections and for the evaluation of function isosurfaces. """ -function calculatCe_plane_tetrahedron_intersection!( +function calculate_plane_tetrahedron_intersection!( ixcoord, ixvalues, coordinates, From 19982efab8b7506271f83e1e4561edf69b6e1f90 Mon Sep 17 00:00:00 2001 From: Liameloy Date: Mon, 30 Jun 2025 16:15:25 +0200 Subject: [PATCH 8/9] Set-up of experimental folder Moved LinearSimplex.jl to experimental folder Created GridSimplex.jl in experimental folder --- src/experimental/gridsimplex.jl | 58 +++++++++++++++++++++++++ src/{ => experimental}/linearsimplex.jl | 16 +++---- 2 files changed, 66 insertions(+), 8 deletions(-) create mode 100644 src/experimental/gridsimplex.jl rename src/{ => experimental}/linearsimplex.jl (80%) diff --git a/src/experimental/gridsimplex.jl b/src/experimental/gridsimplex.jl new file mode 100644 index 0000000..e91d9dd --- /dev/null +++ b/src/experimental/gridsimplex.jl @@ -0,0 +1,58 @@ +""" + GridSimplex +Grid simplex with point coordinates, visible faces and regions. +""" +struct GridSimplex{D, N, Tv, Ti} + points::SVector{N, SVector{D, Tv}} + visibleFaces::SVector{N, Bool} + region::Ti +end + +points(s::GridSimplex) = s.points +visibleFaces(s::GridSimplex) = s.visibleFaces +region(s::GridSimplex) = s.region + +function GridSimplex(::Type{Val{D}}, points::Union{Vector, Tuple}, visibleFaces::Union{Vector, Tuple}, region::Ti) where {D} + simplexPoints = SVector{D + 1, SVector{D, Float32}}(points...) + simplexVisibleFaces = SVector{D + 1, Bool}(visibleFaces...) + return GridSimplex(simplexPoints, simplexVisibleFaces, region) +end + +function GridSimplex(::Type{Val{1}}, points::AbstractMatrix, visibleFaces::AbstractVector, region::Ti) + @views simplexPoints = SVector{2, SVector{1, Float32}}(points[:,1], points[:,2]) + simplexVisibleFaces = SVector{2, Bool}(visibleFaces[1], visibleFaces[2]) + return GridSimplex(simplexPoints, simplexVisibleFaces, region) +end + +function GridSimplex(::Type{Val{2}}, points::AbstractMatrix, visibleFaces::AbstractVector, region::Ti) + @views simplexPoints = SVector{3, SVector{2, Float32}}(points[:,1], points[:,2], points[:,3]) + simplexVisibleFaces = SVector{3, Bool}(visibleFaces[1], visibleFaces[2], visibleFaces[3]) + return GridSimplex(simplexPoints, simplexVisibleFaces, region) +end + +function GridSimplex(::Type{Val{3}}, points::AbstractMatrix, visibleFaces::AbstractVector, region::Ti) + @views simplexPoints = SVector{4, SVector{3, Float32}}(points[:,1], points[:,2], points[:,3], points[:,4]) + simplexVisibleFaces = SVector{4, Bool}(visibleFaces[1], visibleFaces[2], visibleFaces[3], visibleFaces[4]) + return GridSimplex(simplexPoints, simplexVisibleFaces, region) +end + +function GridSimplex(::Type{Val{1}}, points::AbstractMatrix, visibleFaces::AbstractVector, region::Ti, coordscale) + @views simplexPoints = SVector{2, SVector{1, Float32}}(points[:,1] * coordscale, points[:,2] * coordscale) + simplexVisibleFaces = SVector{2, Bool}(visibleFaces[1], visibleFaces[2]) + return GridSimplex(simplexPoints, simplexVisibleFaces, region) +end + +function GridSimplex(::Type{Val{2}}, points::AbstractMatrix, visibleFaces::AbstractVector, region::Ti, coordscale) + @views simplexPoints = SVector{3, SVector{2, Float32}}(points[:,1] * coordscale, points[:,2] * coordscale, points[:,3] * coordscale) + simplexVisibleFaces = SVector{3, Bool}(visibleFaces[1], visibleFaces[2], visibleFaces[3]) + return GridSimplex(simplexPoints, simplexVisibleFaces, region) +end + +function GridSimplex(::Type{Val{3}}, points::AbstractMatrix, visibleFaces::AbstractVector, region::Ti, coordscale) + @views simplexPoints = SVector{4, SVector{3, Float32}}(points[:,1] * coordscale, points[:,2] * coordscale, points[:,3] * coordscale, points[:,4] * coordscale) + simplexVisibleFaces = SVector{4, Bool}(visibleFaces[1], visibleFaces[2], visibleFaces[3], visibleFaces[4]) + return GridSimplex(simplexPoints, simplexVisibleFaces, region) +end + +GridEdge(points, visibleFaces, region) = GridSimplex(Val{1}, points, visibleFaces, region) +GridTriangle(points, visibleFaces, region) = GridSimplex(Val) diff --git a/src/linearsimplex.jl b/src/experimental/linearsimplex.jl similarity index 80% rename from src/linearsimplex.jl rename to src/experimental/linearsimplex.jl index 783b6c7..29f6167 100644 --- a/src/linearsimplex.jl +++ b/src/experimental/linearsimplex.jl @@ -3,7 +3,7 @@ Linear simplex with point coordinates and node values. """ struct LinearSimplex{D, N, Tv} - points::SVector{N, Point{D, Tv}} + points::SVector{N, SVector{D, Tv}} values::SVector{N, Tv} end @@ -12,44 +12,44 @@ points(s::LinearSimplex) = s.points function LinearSimplex(::Type{Val{D}}, points::Union{Vector, Tuple}, values::Union{Vector, Tuple}) where {D} - spoints = SVector{D + 1, Point{D, Float32}}(points...) + spoints = SVector{D + 1, SVector{D, Float32}}(points...) svalues = SVector{D + 1, Float32}(values...) return LinearSimplex(spoints, svalues) end function LinearSimplex(::Type{Val{1}}, points::AbstractMatrix, values::AbstractVector) - @views spoints = SVector{2, Point{1, Float32}}(points[:, 1], points[:, 2]) + @views spoints = SVector{2, SVector{1, Float32}}(points[:, 1], points[:, 2]) svalues = SVector{2, Float32}(values[1], values[2]) return LinearSimplex(spoints, svalues) end function LinearSimplex(::Type{Val{2}}, points::AbstractMatrix, values::AbstractVector) - @views spoints = SVector{3, Point{2, Float32}}(points[:, 1], points[:, 2], points[:, 3]) + @views spoints = SVector{3, SVector{2, Float32}}(points[:, 1], points[:, 2], points[:, 3]) svalues = SVector{3, Float32}(values[1], values[2], values[3]) return LinearSimplex(spoints, svalues) end function LinearSimplex(::Type{Val{3}}, points::AbstractMatrix, values::AbstractVector) - @views spoints = SVector{4, Point{3, Float32}}(points[:, 1], points[:, 2], points[:, 3], points[:, 4]) + @views spoints = SVector{4, SVector{3, Float32}}(points[:, 1], points[:, 2], points[:, 3], points[:, 4]) svalues = SVector{4, Float32}(values[1], values[2], values[3], values[4]) return LinearSimplex(spoints, svalues) end function LinearSimplex(::Type{Val{1}}, points::AbstractMatrix, values::AbstractVector, coordscale) - @views spoints = SVector{2, Point{1, Float32}}(points[:, 1] * coordscale, points[:, 2] * coordscale) + @views spoints = SVector{2, SVector{1, Float32}}(points[:, 1] * coordscale, points[:, 2] * coordscale) svalues = SVector{2, Float32}(values[1], values[2]) return LinearSimplex(spoints, svalues) end function LinearSimplex(::Type{Val{2}}, points::AbstractMatrix, values::AbstractVector, coordscale) - @views spoints = SVector{3, Point{2, Float32}}(points[:, 1] * coordscale, points[:, 2] * coordscale, points[:, 3] * coordscale) + @views spoints = SVector{3, SVector{2, Float32}}(points[:, 1] * coordscale, points[:, 2] * coordscale, points[:, 3] * coordscale) svalues = SVector{3, Float32}(values[1], values[2], values[3]) return LinearSimplex(spoints, svalues) end function LinearSimplex(::Type{Val{3}}, points::AbstractMatrix, values::AbstractVector, coordscale) - @views spoints = SVector{4, Point{3, Float32}}(points[:, 1] * coordscale, points[:, 2] * coordscale, points[:, 3] * coordscale, points[:, 4] * coordscale) + @views spoints = SVector{4, SVector{3, Float32}}(points[:, 1] * coordscale, points[:, 2] * coordscale, points[:, 3] * coordscale, points[:, 4] * coordscale) svalues = SVector{4, Float32}(values[1], values[2], values[3], values[4]) return LinearSimplex(spoints, svalues) end From 3f6973e25e31d2f8fc7591c54a88beb214c78fd1 Mon Sep 17 00:00:00 2001 From: Liameloy Date: Tue, 1 Jul 2025 15:53:18 +0200 Subject: [PATCH 9/9] Code quality fixes --- src/experimental/gridsimplex.jl | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/experimental/gridsimplex.jl b/src/experimental/gridsimplex.jl index e91d9dd..f8d086a 100644 --- a/src/experimental/gridsimplex.jl +++ b/src/experimental/gridsimplex.jl @@ -14,44 +14,44 @@ region(s::GridSimplex) = s.region function GridSimplex(::Type{Val{D}}, points::Union{Vector, Tuple}, visibleFaces::Union{Vector, Tuple}, region::Ti) where {D} simplexPoints = SVector{D + 1, SVector{D, Float32}}(points...) - simplexVisibleFaces = SVector{D + 1, Bool}(visibleFaces...) + simplexVisibleFaces = SVector{D + 1, Bool}(visibleFaces...) return GridSimplex(simplexPoints, simplexVisibleFaces, region) end function GridSimplex(::Type{Val{1}}, points::AbstractMatrix, visibleFaces::AbstractVector, region::Ti) - @views simplexPoints = SVector{2, SVector{1, Float32}}(points[:,1], points[:,2]) + @views simplexPoints = SVector{2, SVector{1, Float32}}(points[:, 1], points[:, 2]) simplexVisibleFaces = SVector{2, Bool}(visibleFaces[1], visibleFaces[2]) return GridSimplex(simplexPoints, simplexVisibleFaces, region) end function GridSimplex(::Type{Val{2}}, points::AbstractMatrix, visibleFaces::AbstractVector, region::Ti) - @views simplexPoints = SVector{3, SVector{2, Float32}}(points[:,1], points[:,2], points[:,3]) + @views simplexPoints = SVector{3, SVector{2, Float32}}(points[:, 1], points[:, 2], points[:, 3]) simplexVisibleFaces = SVector{3, Bool}(visibleFaces[1], visibleFaces[2], visibleFaces[3]) return GridSimplex(simplexPoints, simplexVisibleFaces, region) end function GridSimplex(::Type{Val{3}}, points::AbstractMatrix, visibleFaces::AbstractVector, region::Ti) - @views simplexPoints = SVector{4, SVector{3, Float32}}(points[:,1], points[:,2], points[:,3], points[:,4]) + @views simplexPoints = SVector{4, SVector{3, Float32}}(points[:, 1], points[:, 2], points[:, 3], points[:, 4]) simplexVisibleFaces = SVector{4, Bool}(visibleFaces[1], visibleFaces[2], visibleFaces[3], visibleFaces[4]) - return GridSimplex(simplexPoints, simplexVisibleFaces, region) + return GridSimplex(simplexPoints, simplexVisibleFaces, region) end function GridSimplex(::Type{Val{1}}, points::AbstractMatrix, visibleFaces::AbstractVector, region::Ti, coordscale) - @views simplexPoints = SVector{2, SVector{1, Float32}}(points[:,1] * coordscale, points[:,2] * coordscale) + @views simplexPoints = SVector{2, SVector{1, Float32}}(points[:, 1] * coordscale, points[:, 2] * coordscale) simplexVisibleFaces = SVector{2, Bool}(visibleFaces[1], visibleFaces[2]) return GridSimplex(simplexPoints, simplexVisibleFaces, region) end function GridSimplex(::Type{Val{2}}, points::AbstractMatrix, visibleFaces::AbstractVector, region::Ti, coordscale) - @views simplexPoints = SVector{3, SVector{2, Float32}}(points[:,1] * coordscale, points[:,2] * coordscale, points[:,3] * coordscale) + @views simplexPoints = SVector{3, SVector{2, Float32}}(points[:, 1] * coordscale, points[:, 2] * coordscale, points[:, 3] * coordscale) simplexVisibleFaces = SVector{3, Bool}(visibleFaces[1], visibleFaces[2], visibleFaces[3]) return GridSimplex(simplexPoints, simplexVisibleFaces, region) end function GridSimplex(::Type{Val{3}}, points::AbstractMatrix, visibleFaces::AbstractVector, region::Ti, coordscale) - @views simplexPoints = SVector{4, SVector{3, Float32}}(points[:,1] * coordscale, points[:,2] * coordscale, points[:,3] * coordscale, points[:,4] * coordscale) + @views simplexPoints = SVector{4, SVector{3, Float32}}(points[:, 1] * coordscale, points[:, 2] * coordscale, points[:, 3] * coordscale, points[:, 4] * coordscale) simplexVisibleFaces = SVector{4, Bool}(visibleFaces[1], visibleFaces[2], visibleFaces[3], visibleFaces[4]) - return GridSimplex(simplexPoints, simplexVisibleFaces, region) + return GridSimplex(simplexPoints, simplexVisibleFaces, region) end GridEdge(points, visibleFaces, region) = GridSimplex(Val{1}, points, visibleFaces, region)