diff --git a/CHANGELOG.md b/CHANGELOG.md index 4fe035d..d26b06a 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/Project.toml b/Project.toml index 5c8e008..cb35534 100644 --- a/Project.toml +++ b/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" diff --git a/src/GridVisualizeTools.jl b/src/GridVisualizeTools.jl index 0a80328..a2c2563 100644 --- a/src/GridVisualizeTools.jl +++ b/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 diff --git a/src/marching.jl b/src/marching.jl index 5686e94..1f9b2b9 100644 --- a/src/marching.jl +++ b/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 > 4 + @warn "computed $intersection_count intersection points of a tetrahedron and a plane. Expected at most 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) diff --git a/test/runtests.jl b/test/runtests.jl index 021eb52..6d8477d 100644 --- a/test/runtests.jl +++ b/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))