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

Author: Liameloy

CHANGELOG.md

@@ -2,6 +2,12 @@
 
 All notable changes to this project will be documented in this file.
 
+## Future Release 
+
+### Fixed
+
+- Reworked `tet_x_plane!()` to improve tetrahedron-plane-crosssection
+
 ## [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"

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

@@ -44,23 +44,32 @@ function tet_x_plane!(
     end
     # Interpolate coordinates and function_values according to
     # evaluation of the plane equation
-    nxs = 0
+    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] &&
-                    planeq_values[n1] * planeq_values[n2] < tol
-                nxs += 1
+
+            if planeq_values[n1] * planeq_values[n2] < tol && !visited_list[n1] && !visited_list[n2]
+
+                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, 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])
+                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])
             end
         end
     end
-    return nxs
+    if !(intersection_count in [3, 4])
+        @warn "computed $intersection_count intersection points of a tetrahedron and a plane. Expected value is 3 or 4."
+    end
+    return intersection_count
 end
 
 """
@@ -190,7 +199,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 +210,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)

test/runtests.jl

@@ -3,6 +3,21 @@ using Test, Documenter, GridVisualizeTools, ColorTypes, Colors
 DocMeta.setdocmeta!(GridVisualizeTools, :DocTestSetup, :(using GridVisualizeTools, ColorTypes, Colors); recursive = true)
 doctest(GridVisualizeTools)
 
+
+@testset "tet_x_plane" 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!(
+        [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
+end
+
+
 if isdefined(Docs, :undocumented_names) # >=1.11
     @testset "UndocumentedNames" begin
         @test isempty(Docs.undocumented_names(GridVisualizeTools))