Fixed calculation for intersections

Adjusted the logic which determines an intersection between tetrahedron and plane
Added new test case for specific data which caused issues beforehand

Author: Liameloy

src/marching.jl

@@ -25,51 +25,77 @@
   This method can be used both for the evaluation of plane sections and for
   the evaluation of function isosurfaces.
 """
-function tet_x_plane!(
+function assign_coordinates_at_point_in_plane!(
         ixcoord,
         ixvalues,
-        pointlist,
+        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
+
+function calculate_plane_tetrahedron_intersection!(
+        ixcoord,
+        ixvalues,
+        coordinates,
         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)
-        )
-        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 +257,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

@@ -4,9 +4,9 @@ DocMeta.setdocmeta!(GridVisualizeTools, :DocTestSetup, :(using GridVisualizeTool
 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 +15,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