Skip to content

Add line definitions to marching_triangles #20

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 4 commits into from
Mar 3, 2025
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 8 additions & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
@@ -2,6 +2,14 @@

All notable changes to this project will be documented in this file.

## [3.0.0] - 2025-03-03

### 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

4 changes: 2 additions & 2 deletions Project.toml
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
name = "GridVisualizeTools"
uuid = "5573ae12-3b76-41d9-b48c-81d0b6e61cc5"
authors = ["Jürgen Fuhrmann <juergen-fuhrmann@web.de>"]
version = "2.0.0"
authors = ["Jürgen Fuhrmann <juergen-fuhrmann@web.de>", "Patrick Jaap <patrick.jaap@wias-berlin.de"]
version = "2.1.0"

[deps]
ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
142 changes: 106 additions & 36 deletions src/marching.jl
Original file line number Diff line number Diff line change
@@ -259,43 +259,81 @@ 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
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
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{T},
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)
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


"""
$(SIGNATURES)


Variant of `marching_triangles` with multiple grid input
"""
function marching_triangles(
coords::Vector{Matrix{Tv}},
coords::Vector{Matrix{T}},
cellnodes::Vector{Matrix{Ti}},
funcs,
lines,
levels;
Tc = Float32,
Tp = SVector{2, Tc}
) where {Tv <: 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)

for igrid in 1:length(coords)
func = funcs[igrid]
coord = coords[igrid]

function isect(nodes)
# 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(tri_nodes, objective_func, value_func)
(i1, i2, i3) = (1, 2, 3)

f = (func[nodes[1]], func[nodes[2]], 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]
@@ -305,35 +343,67 @@ 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])

# 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,
tri_nodes
)
@views isect(tri_nodes, 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,
tri_nodes
)
@views isect(tri_nodes, objective_values, func)
end

end
end
return points
return points, adjacencies, values
end
Loading
Oops, something went wrong.