Fixed issue of missing shapes in sliceplots of given grids and planes

Fixed error in tet_x_plane function in marching.jl which incorrectly assessed intersections between tetrahedron edges and plane

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

Adjusted index.md for documentation

Adjusting files to make documentation buildable
Will be squashed later

Changed position of docstring for Documentation

Further adjustments for the documentation

Changed the change log

Author: Liameloy

CHANGELOG.md

@@ -2,6 +2,14 @@
 
 All notable changes to this project will be documented in this file.
 
+## Future Release 
+
+### 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.0] - 2025-03-03
 
 ### Changed

Project.toml

@@ -7,11 +7,13 @@ version = "3.0.0"
 ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4"
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
 
 [compat]
 ColorSchemes = "3"
 Colors = "0.12, 0.13"
 DocStringExtensions = "0.8, 0.9"
+StaticArrays = "1.9.13"
 StaticArraysCore = "1.4"
 julia = "1.6"

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/GridVisualizeTools.jl

@@ -10,6 +10,7 @@ import ColorSchemes
 
 using DocStringExtensions: SIGNATURES, TYPEDEF, TYPEDSIGNATURES
 using StaticArraysCore: SVector
+using StaticArrays: @MArray
 
 include("colors.jl")
 export region_cmap, bregion_cmap, rgbtuple, rgbcolor

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,47 +56,47 @@
   - 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;
         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
-    nxs = 0
-    @inbounds @simd for n1 in 1:3
+
+    amount_intersections = 0
+
+    @inbounds for n1 in 1:4
         N1 = node_indices[n1]
-        @inbounds @fastmath @simd for n2 in (n1 + 1):4
-            N2 = node_indices[n2]
-            if planeq_values[n1] != planeq_values[n2] &&
-                    planeq_values[n1] * planeq_values[n2] < tol
-                nxs += 1
-                t = planeq_values[n1] / (planeq_values[n1] - planeq_values[n2])
-                ixcoord[1, nxs] = pointlist[1, N1] + t * (pointlist[1, N2] - pointlist[1, N1])
-                ixcoord[2, nxs] = pointlist[2, N1] + t * (pointlist[2, N2] - pointlist[2, N1])
-                ixcoord[3, nxs] = pointlist[3, N1] + t * (pointlist[3, N2] - pointlist[3, N1])
-                ixvalues[nxs] = function_values[N1] + t * (function_values[N2] - function_values[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
-    return nxs
+
+    if amount_intersections > 4
+        @warn "computed $(amount_intersections) intersection points of a tetrahedron and a plane. Expected at most 4."
+    end
+    return amount_intersections
 end
 
 """
@@ -190,7 +226,7 @@ function marching_tetrahedra(
 
     function pushtris(ns, ixcoord, ixvalues)
         # number of intersection points can be 3 or 4
-        return if ns >= 3
+        if ns >= 3
             last_i = length(all_ixvalues)
             for is in 1:ns
                 @views push!(all_ixcoord, ixcoord[:, is])
@@ -201,6 +237,7 @@ function marching_tetrahedra(
                 push!(all_ixfaces, (last_i + 3, last_i + 2, last_i + 4))
             end
         end
+        return nothing
     end
 
     for igrid in 1:length(allcoords)
@@ -221,7 +258,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,6 +3,31 @@ using Test, Documenter, GridVisualizeTools, ColorTypes, Colors
 DocMeta.setdocmeta!(GridVisualizeTools, :DocTestSetup, :(using GridVisualizeTools, ColorTypes, Colors); recursive = true)
 doctest(GridVisualizeTools)
 
+
+@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.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],
+        Int32[1, 2, 4, 8],
+        [-1.0, 0.0, 1.0, 1.0],
+        [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
+
+
 if isdefined(Docs, :undocumented_names) # >=1.11
     @testset "UndocumentedNames" begin
         @test isempty(Docs.undocumented_names(GridVisualizeTools))