Fix issue with PointLocationHistory when ghost triangles are included

JuliaGeometry · Sep 27, 2023 · 40c6c34 · 40c6c34
1 parent 01eb7a2
commit 40c6c34
Show file tree

Hide file tree

Showing 2 changed files with 22 additions and 19 deletions.
diff --git a/src/point_location/initialisers.jl b/src/point_location/initialisers.jl
@@ -477,14 +477,14 @@ function check_for_intersections_with_interior_edges_adjacent_to_boundary_node(
                 # Here, q is on the edge, so we consider it as being inside the triangle.
                 if is_true(store_history)
                     add_triangle!(history, i, j, k)
-                    add_edge!(history, i, j)
+                    pᵢ ≠ q && pⱼ ≠ q && add_edge!(history, i, j) # be careful not to add a collinear edge that is just starting at q
                 end
                 return i, j, Cert.On, Cert.Inside
             elseif is_right(q_cert)
                 # Here, q is to the right of ppⱼ, which just means that it is inside the triangulation. 
                 if is_true(store_history)
                     add_triangle!(history, i, j, k)
-                    add_edge!(history, k, j)
+                    p ≠ q && pⱼ ≠ q && add_edge!(history, k, j) # be careful not to add a collinear edge that is just starting at q
                 end
                 return j, i, Cert.Right, Cert.Outside # flip i and j to get the correct orientation
             elseif is_left(q_cert)

diff --git a/test/data_structures/triangulation.jl b/test/data_structures/triangulation.jl
@@ -605,23 +605,26 @@ end
       n = 60
       for _ in 1:10
             n += rand(1:125)
-            tri = triangulate(12randn(2, n))
-            for qi in each_solid_vertex(tri)
-                  for k in each_solid_vertex(tri)
-                        q = get_point(tri, qi)
-                        history = DT.PointLocationHistory{NTuple{3,Int},NTuple{2,Int},Int}()
-                        jump_and_march(tri, q;
-                              k,
-                              store_history=true,
-                              history)
-                        visited_triangles = history.triangles
-                        collinear_segments = history.collinear_segments
-                        @test all(T -> DT.is_positively_oriented(DT.triangle_orientation(tri, T)), visited_triangles)
-                        @test all(!DT.is_none, [DT.triangle_line_segment_intersection(tri, T, (qi, k)) for T in visited_triangles])
-                        @test allunique(visited_triangles)
-                        if !isempty(collinear_segments)
-                              @test all(E -> DT.is_collinear(DT.point_position_relative_to_line(tri, qi, k, E[1])), collinear_segments)
-                              @test all(E -> DT.is_collinear(DT.point_position_relative_to_line(tri, qi, k, E[2])), collinear_segments)
+            tri1 = triangulate(12randn(2, n), delete_ghosts=false)
+            tri2 = triangulate(12randn(2, n), delete_ghosts=true)
+            for tri in (tri1, tri2)
+                  for qi in each_solid_vertex(tri)
+                        for k in each_solid_vertex(tri)
+                              q = get_point(tri, qi)
+                              history = DT.PointLocationHistory{NTuple{3,Int},NTuple{2,Int},Int}()
+                              jump_and_march(tri, q;
+                                    k,
+                                    store_history=true,
+                                    history)
+                              visited_triangles = history.triangles
+                              collinear_segments = history.collinear_segments
+                              @test all(T -> DT.is_positively_oriented(DT.triangle_orientation(tri, T)), visited_triangles)
+                              @test all(!DT.is_none, [DT.triangle_line_segment_intersection(tri, T, (qi, k)) for T in visited_triangles])
+                              @test allunique(visited_triangles)
+                              if !isempty(collinear_segments)
+                                    @test all(E -> DT.is_collinear(DT.point_position_relative_to_line(tri, qi, k, E[1])), collinear_segments)
+                                    @test all(E -> DT.is_collinear(DT.point_position_relative_to_line(tri, qi, k, E[2])), collinear_segments)
+                              end
                         end
                   end
             end