diff --git a/CHANGELOG.md b/CHANGELOG.md index 25d93a4..bbc25d8 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/Project.toml b/Project.toml index 2741644..dc15685 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "GridVisualizeTools" uuid = "5573ae12-3b76-41d9-b48c-81d0b6e61cc5" -authors = ["Jürgen Fuhrmann ", "Patrick Jaap ", "Patrick Jaap ", "Liam Johnen "] +version = "3.0.2" [deps] ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4" diff --git a/docs/src/index.md b/docs/src/index.md index d7e66b4..8576b7a 100644 --- a/docs/src/index.md +++ b/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! ``` diff --git a/src/marching.jl b/src/marching.jl index 1f9b2b9..85899ab 100644 --- a/src/marching.jl +++ b/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, diff --git a/test/runtests.jl b/test/runtests.jl index 6d8477d..af736bf 100644 --- a/test/runtests.jl +++ b/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