From bf8d2e958d7265738aaf10ba58fa9f75456acaed Mon Sep 17 00:00:00 2001 From: Patrick Jaap Date: Mon, 27 Jan 2025 16:53:25 +0100 Subject: [PATCH 1/4] Add line argument to marching_triangles and return additional values The method now also returns the values of the objective function at the intersection points and the adjacency information. Therefore, it is possible to construct a 1D grid from the return values. --- CHANGELOG.md | 8 +++ src/marching.jl | 130 +++++++++++++++++++++++++++++++++++------------- 2 files changed, 103 insertions(+), 35 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index a3afbfe..2f5bdd5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,14 @@ All notable changes to this project will be documented in this file. +## [unreleased] + +### Changed + +- `marching_triangles` takes now an additional `lines` argument to compute intersections with given lines and + returns also values and adjacency information to be able to construct a 1D grid + To restore the old behavior, use `lines = []` and neglect the second and third return value + ## [1.1.1] - 2024-11-02 - Update links after move to WIAS-PDELib org diff --git a/src/marching.jl b/src/marching.jl index 49fba57..8023846 100644 --- a/src/marching.jl +++ b/src/marching.jl @@ -259,38 +259,70 @@ end """ $(SIGNATURES) -Collect isoline snippets on triangles ready for linesegments! +March through the given grid and extract points and values for given iso-line levels and/or given intersection lines. +From the returned point list and value list a line plot can be created. + +Input: + coord: matrix storing the coordinates of the grid + cellnodes: connectivity matrix + func: function on the grid nodes to be evaluated + lines: vector of line definitions [a,b,c], s.t., ax + by + c = 0 defines a line + levels: vector of levels for the iso-surface + +Output: + points: vector of 2D points of the intersections of the grid with the iso-surfaces or lines + adjacencies: vector of 2D vectors storing connected points in the grid + value: interpolated values of `func` at the intersection points + +Note that passing both nonempty `lines` and `levels` will create a result with both types of points mixed. """ function marching_triangles( - coord::Matrix{Tv}, + coord::Matrix{Tc}, cellnodes::Matrix{Ti}, func, + lines, levels; - Tc = Float32, - Tp = SVector{2, Tc} - ) where {Tv <: Number, Ti <: Number} - return marching_triangles([coord], [cellnodes], [func], levels; Tc, Tp) + Tv = Float32, + Tp = SVector{2, Tv} + ) where {Tc <: Number, Ti <: Number} + return marching_triangles([coord], [cellnodes], [func], lines, levels; Tv, Tp) end + +""" + $(SIGNATURES) + + +Variant of `marching_tetrahedra` with multiple grid input +""" function marching_triangles( - coords::Vector{Matrix{Tv}}, + coords::Vector{Matrix{Tc}}, cellnodes::Vector{Matrix{Ti}}, funcs, + lines, levels; - Tc = Float32, - Tp = SVector{2, Tc} - ) where {Tv <: Number, Ti <: Number} + Tv = Float32, + Tp = SVector{2, Tv}, + ) where {Tc <: Number, Ti <: Number} points = Vector{Tp}(undef, 0) + values = Vector{Tv}(undef, 0) + adjacencies = Vector{SVector{2, Ti}}(undef, 0) for igrid in 1:length(coords) func = funcs[igrid] coord = coords[igrid] - function isect(nodes) + # pre-allcate memory + objective_values = Vector{Tv}(undef, size(coord, 2)) + + # the objective_func is used to determine the intersection (line equation or iso levels) + # the value_func is used to interpolate values at the intersections + function isect(nodes, objective_func, value_func) (i1, i2, i3) = (1, 2, 3) - f = (func[nodes[1]], func[nodes[2]], func[nodes[3]]) + f = (objective_func[nodes[1]], objective_func[nodes[2]], objective_func[nodes[3]]) + # sort f[i1] ≤ f[i2] ≤ f[i3] f[1] <= f[2] ? (i1, i2) = (1, 2) : (i1, i2) = (2, 1) f[i2] <= f[3] ? i3 = 3 : (i2, i3) = (3, i2) f[i1] > f[i2] ? (i1, i2) = (i2, i1) : nothing @@ -305,35 +337,63 @@ function marching_triangles( dy21 = coord[2, n2] - coord[2, n1] dy32 = coord[2, n3] - coord[2, n2] - df31 = f[i3] != f[i1] ? 1 / (f[i3] - f[i1]) : 0.0 - 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 in levels - if (f[i1] <= level) && (level < f[i3]) - α = (level - f[i1]) * df31 - x1 = coord[1, n1] + α * dx31 - y1 = coord[2, n1] + α * dy31 - - if (level < f[i2]) - α = (level - f[i1]) * df21 - x2 = coord[1, n1] + α * dx21 - y2 = coord[2, n1] + α * dy21 - else - α = (level - f[i2]) * df32 - x2 = coord[1, n2] + α * dx32 - y2 = coord[2, n2] + α * dy32 - end - push!(points, SVector{2, Tc}((x1, y1))) - push!(points, SVector{2, Tc}((x2, y2))) + df31 = f[i3] != f[i1] ? 1 / (f[i1] - f[i3]) : 0.0 + df21 = f[i2] != f[i1] ? 1 / (f[i1] - f[i2]) : 0.0 + df32 = f[i3] != f[i2] ? 1 / (f[i2] - f[i3]) : 0.0 + + if (f[i1] <= 0) && (0 < f[i3]) + α = f[i1] * df31 + x1 = coord[1, n1] + α * dx31 + y1 = coord[2, n1] + α * dy31 + value1 = value_func[n1] + α * (value_func[n3] - value_func[n1]) + + if (0 < f[i2]) + α = f[i1] * df21 + x2 = coord[1, n1] + α * dx21 + y2 = coord[2, n1] + α * dy21 + value2 = value_func[n1] + α * (value_func[n2] - value_func[n1]) + else + α = f[i2] * df32 + x2 = coord[1, n2] + α * dx32 + y2 = coord[2, n2] + α * dy32 + value2 = value_func[n2] + α * (value_func[n3] - value_func[n2]) end + + push!(points, SVector{2, Tc}((x1, y1))) + push!(points, SVector{2, Tc}((x2, y2))) + push!(values, value1) + push!(values, value2) + # connect last two points + push!(adjacencies, SVector{2, Ti}((length(points) - 1, length(points)))) end + return end for itri in 1:size(cellnodes[igrid], 2) - @views isect(cellnodes[igrid][:, itri]) + for level in levels + # objective func is iso-level equation + @views @fastmath map!( + inode -> (func[inode] - level), + objective_values, + 1:size(coord, 2) + ) + @views isect(cellnodes[igrid][:, itri], objective_values, func) + end + + for line in lines + @fastmath line_equation(line, coord) = coord[1] * line[1] + coord[2] * line[2] + line[3] + + # objective func is iso-level equation + @views @fastmath map!( + inode -> (line_equation(line, coord[:, inode])), + objective_values, + 1:size(coord, 2) + ) + @views isect(cellnodes[igrid][:, itri], objective_values, func) + end + end end - return points + return points, adjacencies, values end From 1b95b8e3557bbb97325678bfd6a6b06c0291d52d Mon Sep 17 00:00:00 2001 From: Patrick Jaap Date: Mon, 3 Mar 2025 09:47:33 +0100 Subject: [PATCH 2/4] marching_triangles: improve loop speed by restriction to necessary nodes --- src/marching.jl | 23 ++++++++++++++--------- 1 file changed, 14 insertions(+), 9 deletions(-) diff --git a/src/marching.jl b/src/marching.jl index 8023846..eb37582 100644 --- a/src/marching.jl +++ b/src/marching.jl @@ -312,22 +312,23 @@ function marching_triangles( func = funcs[igrid] coord = coords[igrid] - # pre-allcate memory - objective_values = Vector{Tv}(undef, size(coord, 2)) + # pre-allcate memory for triangle values (3 nodes per triangle) + objective_values = Vector{Tv}(undef, 3) # the objective_func is used to determine the intersection (line equation or iso levels) # the value_func is used to interpolate values at the intersections - function isect(nodes, objective_func, value_func) + function isect(tri_nodes, objective_func, value_func) (i1, i2, i3) = (1, 2, 3) - f = (objective_func[nodes[1]], objective_func[nodes[2]], objective_func[nodes[3]]) + # 3 values of the objective function + f = objective_func # sort f[i1] ≤ f[i2] ≤ f[i3] f[1] <= f[2] ? (i1, i2) = (1, 2) : (i1, i2) = (2, 1) f[i2] <= f[3] ? i3 = 3 : (i2, i3) = (3, i2) f[i1] > f[i2] ? (i1, i2) = (i2, i1) : nothing - (n1, n2, n3) = (nodes[i1], nodes[i2], nodes[i3]) + (n1, n2, n3) = (tri_nodes[i1], tri_nodes[i2], tri_nodes[i3]) dx31 = coord[1, n3] - coord[1, n1] dx21 = coord[1, n2] - coord[1, n1] @@ -371,14 +372,18 @@ function marching_triangles( end for itri in 1:size(cellnodes[igrid], 2) + + # nodes of the current triangle + tri_nodes = @views cellnodes[igrid][:, itri] + for level in levels # objective func is iso-level equation @views @fastmath map!( inode -> (func[inode] - level), objective_values, - 1:size(coord, 2) + tri_nodes ) - @views isect(cellnodes[igrid][:, itri], objective_values, func) + @views isect(tri_nodes, objective_values, func) end for line in lines @@ -388,9 +393,9 @@ function marching_triangles( @views @fastmath map!( inode -> (line_equation(line, coord[:, inode])), objective_values, - 1:size(coord, 2) + tri_nodes ) - @views isect(cellnodes[igrid][:, itri], objective_values, func) + @views isect(tri_nodes, objective_values, func) end end From 51f3b60c16e0a85b59f208131c6b3bf9486645a3 Mon Sep 17 00:00:00 2001 From: Patrick Jaap Date: Mon, 3 Mar 2025 12:00:16 +0100 Subject: [PATCH 3/4] More user control over the types in marching_triangles --- src/marching.jl | 25 +++++++++++++++---------- 1 file changed, 15 insertions(+), 10 deletions(-) diff --git a/src/marching.jl b/src/marching.jl index eb37582..5686e94 100644 --- a/src/marching.jl +++ b/src/marching.jl @@ -268,6 +268,9 @@ Input: func: function on the grid nodes to be evaluated lines: vector of line definitions [a,b,c], s.t., ax + by + c = 0 defines a line levels: vector of levels for the iso-surface + Tc: scalar type of coordinates + Tp: vector type of coordinates + Tv: scalar type of function values Output: points: vector of 2D points of the intersections of the grid with the iso-surfaces or lines @@ -277,15 +280,16 @@ Output: Note that passing both nonempty `lines` and `levels` will create a result with both types of points mixed. """ function marching_triangles( - coord::Matrix{Tc}, + coord::Matrix{T}, cellnodes::Matrix{Ti}, func, lines, levels; - Tv = Float32, - Tp = SVector{2, Tv} - ) where {Tc <: Number, Ti <: Number} - return marching_triangles([coord], [cellnodes], [func], lines, levels; Tv, Tp) + Tc = T, + Tp = SVector{2, Tc}, + Tv = Float64 + ) where {T <: Number, Ti <: Number} + return marching_triangles([coord], [cellnodes], [func], lines, levels; Tc, Tp, Tv) end @@ -293,17 +297,18 @@ end $(SIGNATURES) -Variant of `marching_tetrahedra` with multiple grid input +Variant of `marching_triangles` with multiple grid input """ function marching_triangles( - coords::Vector{Matrix{Tc}}, + coords::Vector{Matrix{T}}, cellnodes::Vector{Matrix{Ti}}, funcs, lines, levels; - Tv = Float32, - Tp = SVector{2, Tv}, - ) where {Tc <: Number, Ti <: Number} + Tc = T, + Tp = SVector{2, Tc}, + Tv = Float64 + ) where {T <: Number, Ti <: Number} points = Vector{Tp}(undef, 0) values = Vector{Tv}(undef, 0) adjacencies = Vector{SVector{2, Ti}}(undef, 0) From bd6a91f9adda9722eeb398ccc333fdfefbb2e771 Mon Sep 17 00:00:00 2001 From: Patrick Jaap Date: Mon, 3 Mar 2025 13:04:48 +0100 Subject: [PATCH 4/4] Add myself to authors and bump version --- CHANGELOG.md | 2 +- Project.toml | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 2f5bdd5..4fe035d 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,7 +2,7 @@ All notable changes to this project will be documented in this file. -## [unreleased] +## [3.0.0] - 2025-03-03 ### Changed diff --git a/Project.toml b/Project.toml index 507a1ba..0645e39 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "GridVisualizeTools" uuid = "5573ae12-3b76-41d9-b48c-81d0b6e61cc5" -authors = ["Jürgen Fuhrmann "] -version = "2.0.0" +authors = ["Jürgen Fuhrmann ", "Patrick Jaap