Skip to content

Fixed calculation for intersections #25

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 2 commits into from
May 27, 2025
Merged
Show file tree
Hide file tree
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
Expand Up @@ -2,6 +2,14 @@

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

## [3.0.2] - 2025-05-22

### Fixed

- Renamed `tet_x_plane!()` to `calculate_plane_tetrahedron_intersection!()`
- Adjusted logic in `calculate_plane_tetrahedron_intersection!()` to better determine intersection
- Separated functionality of intersection determination and coordinate assignment in `calculate_plane_tetrahedron_intersection!()` into subfunctions

## [3.0.1] - 2025-04-16

### Fixed
Expand Down
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 <[email protected]>", "Patrick Jaap <[email protected]"]
version = "3.0.1"
authors = ["Jürgen Fuhrmann <[email protected]>", "Patrick Jaap <[email protected]>", "Liam Johnen <[email protected]>"]
version = "3.0.2"

[deps]
ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
Expand Down
2 changes: 1 addition & 1 deletion docs/src/index.md
Original file line number Diff line number Diff line change
Expand Up @@ -122,5 +122,5 @@ makeisolevels(collect(0:0.1:10), 1, (1,-1),nothing)
## Private API

```@docs
GridVisualizeTools.tet_x_plane!
GridVisualizeTools.calculate_plane_tetrahedron_intersection!
```
90 changes: 62 additions & 28 deletions src/marching.jl
Original file line number Diff line number Diff line change
@@ -1,3 +1,39 @@
function assign_coordinates_at_point_in_plane!(
ixcoord,
ixvalues,
coordinates,
function_values,
node_index,
intersection_count
)

@views ixcoord[:, intersection_count] .= coordinates[:, node_index]
ixvalues[intersection_count] = function_values[node_index]

return nothing

end

function assign_coordinates_at_intersection!(
ixcoord,
ixvalues,
coordinates,
function_values,
planeq_values_start,
planeq_values_end,
node_index_start,
node_index_end,
intersection_count
)

t = planeq_values_start / (planeq_values_start - planeq_values_end)
@views @. ixcoord[:, intersection_count] = coordinates[:, node_index_start] + t * (coordinates[:, node_index_end] - coordinates[:, node_index_start])
ixvalues[intersection_count] = function_values[node_index_start] + t * (function_values[node_index_end] - function_values[node_index_start])

return nothing
end


"""
$(SIGNATURES)
Calculate intersections between tetrahedron with given piecewise linear
Expand All @@ -10,8 +46,8 @@
and the plane.

Input:
- pointlist: 3xN array of grid point coordinates
- node_indices: 4 element array of node indices (pointing into pointlist and function_values)
- coordinates: 3xN array of grid point coordinates
- node_indices: 4 element array of node indices (pointing into coordinates and function_values)
- planeq_values: 4 element array of plane equation evaluated at the node coordinates
- function_values: N element array of function values

Expand All @@ -20,15 +56,15 @@
- ixvalues: 4 element array of function values at plane - tetdedge intersections

Returns:
- nxs,ixcoord,ixvalues
- amount_intersections,ixcoord,ixvalues

This method can be used both for the evaluation of plane sections and for
the evaluation of function isosurfaces.
"""
function tet_x_plane!(
function calculate_plane_tetrahedron_intersection!(
ixcoord,
ixvalues,
pointlist,
coordinates,
node_indices,
planeq_values,
function_values;
Expand All @@ -42,34 +78,32 @@ function tet_x_plane!(
)
return 0
end
# Interpolate coordinates and function_values according to
# evaluation of the plane equation
intersection_count = 0
# List to check whether a node has already been visited, when checking for possible intersections
visited_list = @MArray zeros(Bool, 4)
@inbounds @simd for n1 in 1:3
N1 = node_indices[n1]
@inbounds @fastmath @simd for n2 in (n1 + 1):4
N2 = node_indices[n2]

if planeq_values[n1] * planeq_values[n2] < tol && !visited_list[n1] && !visited_list[n2]
amount_intersections = 0

abs(planeq_values[n1]) < tol && (visited_list[n1] = true)
abs(planeq_values[n2]) < tol && (visited_list[n2] = true)

intersection_count += 1
t = planeq_values[n1] / (planeq_values[n1] - planeq_values[n2])
ixcoord[1, intersection_count] = pointlist[1, N1] + t * (pointlist[1, N2] - pointlist[1, N1])
ixcoord[2, intersection_count] = pointlist[2, N1] + t * (pointlist[2, N2] - pointlist[2, N1])
ixcoord[3, intersection_count] = pointlist[3, N1] + t * (pointlist[3, N2] - pointlist[3, N1])
ixvalues[intersection_count] = function_values[N1] + t * (function_values[N2] - function_values[N1])
@inbounds for n1 in 1:4
N1 = node_indices[n1]
if abs(planeq_values[n1]) < tol
amount_intersections += 1
assign_coordinates_at_point_in_plane!(ixcoord, ixvalues, coordinates, function_values, N1, amount_intersections)
else
for n2 in (n1 + 1):4
N2 = node_indices[n2]
if (abs(planeq_values[n2]) < tol) # We do not allow the 2nd node to be in the plane
continue
end
if planeq_values[n1] * planeq_values[n2] < tol^2
amount_intersections += 1
assign_coordinates_at_intersection!(ixcoord, ixvalues, coordinates, function_values, planeq_values[n1], planeq_values[n2], N1, N2, amount_intersections)
end
end
end
end
if intersection_count > 4
@warn "computed $intersection_count intersection points of a tetrahedron and a plane. Expected at most 4."

if amount_intersections > 4
@warn "computed $(amount_intersections) intersection points of a tetrahedron and a plane. Expected at most 4."
end
return intersection_count
return amount_intersections
end

"""
Expand Down Expand Up @@ -231,7 +265,7 @@ function marching_tetrahedra(
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!(
nxs = calculate_plane_tetrahedron_intersection!(
ixcoord,
ixvalues,
coord,
Expand Down
15 changes: 12 additions & 3 deletions test/runtests.jl
Original file line number Diff line number Diff line change
Expand Up @@ -3,10 +3,9 @@ using Test, Documenter, GridVisualizeTools, ColorTypes, Colors
DocMeta.setdocmeta!(GridVisualizeTools, :DocTestSetup, :(using GridVisualizeTools, ColorTypes, Colors); recursive = true)
doctest(GridVisualizeTools)


@testset "tet_x_plane" begin
@testset "calculate_plane_tetrahedron_intersection" begin
#Testing amount of intersected edges of (x,y,z)-tetrahedron (0,0,0),(1,0,0),(1,1,0),(1,1,1) with plane x+y-1=0
@test GridVisualizeTools.tet_x_plane!(
@test GridVisualizeTools.calculate_plane_tetrahedron_intersection!(
[0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0; 0.0 0.0 0.0 0.0 0.0 0.0],
[0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
[0.0 1.0 0.0 1.0 0.0 1.0 0.0 1.0; 0.0 0.0 1.0 1.0 0.0 0.0 1.0 1.0; 0.0 0.0 0.0 0.0 1.0 1.0 1.0 1.0],
Expand All @@ -15,6 +14,16 @@ doctest(GridVisualizeTools)
[0, 0, 0, 0, 0, 0, 0, 1],
tol = 1.0e-12
) == 3
#Testing amount of intersected edges for specific tetrahedron with plane x = 0.5
@test GridVisualizeTools.calculate_plane_tetrahedron_intersection!(
[0.5 0.5 0.500000001260275 0.5000000012855993 0.0 0.0; 0.7 1.0 0.13480492277555536 0.8372051924096586 0.0 0.0; 0.0 0.1 0.19515864670162444 0.18583544507073008 0.0 0.0],
[1.460741148435986e-32, 1.2130506551371849e-32, 0.19103133044823084, 0.22244570933305352, 0.0, 0.0],
[0.5 0.5 0.48305 0.55059; 0.7165 0.8 0.76944 0.76519; 0.12501 0.0 0.16147 0.16147],
Int32[1, 2, 3, 4],
[0.0, 0.0, -0.016950000077486038, 0.05059000104665756],
[0.18823234602961633, 1.519472645376028e-32, 0.2384620788839228, 0.2550147563453821],
tol = 1.0e-12
) == 3
end


Expand Down
Loading