Merge pull request #25 from Liameloy/dev/tetrahedronSliceFixing

Fixed calculation for intersections

Author: pjaap

CHANGELOG.md

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

Project.toml

@@ -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"

docs/src/index.md

@@ -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!
 ```

src/marching.jl

@@ -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
@@ -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
 
@@ -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;
@@ -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
 
 """
@@ -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,

test/runtests.jl

@@ -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],
@@ -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