From 0df001b2862ca0cfffbb5dba3d7abf5bcc599ec4 Mon Sep 17 00:00:00 2001 From: DanielVandH Date: Sun, 7 Jul 2024 20:12:08 +0100 Subject: [PATCH 01/14] Allow for ExactPredicates to be toggled --- .github/workflows/CI.yml | 11 +- .gitignore | 3 +- Project.toml | 6 +- docs/Manifest.toml | 88 ++--- docs/Project.toml | 1 + docs/src/literate_tutorials/convex.jl | 5 +- docs/src/literate_tutorials/point_location.jl | 9 +- src/DelaunayTriangulation.jl | 70 ++++ src/geometric_primitives/points.jl | 30 +- src/predicates/exactpredicates_definitions.jl | 210 +++++++++++- src/predicates/predicates.jl | 8 - test/Project.toml | 3 + test/interfaces/points.jl | 34 +- test/operations/delete_ghost_triangles.jl | 12 +- test/operations/delete_holes.jl | 114 ++++--- test/operations/delete_point.jl | 82 ++--- test/operations/lock_convex_hull.jl | 1 - test/predicates/general.jl | 84 ++++- test/refinement/curve_bounded.jl | 8 +- test/refinement/refine.jl | 321 +++++++++--------- test/runtests.jl | 6 +- test/triangulation/bowyer_watson.jl | 17 +- test/triangulation/constrained.jl | 204 +++++------ test/voronoi/voronoi.jl | 15 +- 24 files changed, 881 insertions(+), 461 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index 669750353..f4752c309 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -12,7 +12,7 @@ concurrency: cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }} jobs: test: - name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} + name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} - ExactPredicates ${{ matrix.preferences }} runs-on: ${{ matrix.os }} env: JULIA_NUM_THREADS: ${{ matrix.julia-threads }} @@ -21,11 +21,15 @@ jobs: matrix: version: - '1' - - '1.6' + - 'lts' os: - ubuntu-latest arch: - x64 + preferences: + - 'default' + - 'true' + - 'false' steps: - name: Checkout repository for access uses: actions/checkout@v4 @@ -35,6 +39,7 @@ jobs: with: version: ${{ matrix.version }} arch: ${{ matrix.arch }} + - run: julia -e 'using InteractiveUtils; versioninfo(verbose=true)' - name: Cache artifacts, etc. uses: julia-actions/cache@v2 @@ -44,6 +49,8 @@ jobs: - name: Run the tests uses: julia-actions/julia-runtest@v1 + env: + USE_EXACTPREDICATES: ${{ matrix.preferences }} - name: Process the coverage if: always() diff --git a/.gitignore b/.gitignore index c34786994..f6b55de47 100644 --- a/.gitignore +++ b/.gitignore @@ -6,4 +6,5 @@ sw test/Manifest.toml docs/Manifest.toml /docs/build/ -/.vscode \ No newline at end of file +/.vscode +LocalPreferences.toml \ No newline at end of file diff --git a/Project.toml b/Project.toml index 4cb4537bb..1ab26c5fc 100644 --- a/Project.toml +++ b/Project.toml @@ -1,19 +1,21 @@ name = "DelaunayTriangulation" uuid = "927a84f5-c5f4-47a5-9785-b46e178433df" authors = ["Daniel VandenHeuvel "] -version = "1.0.5" +version = "1.1.0-DEV" [deps] EnumX = "4e289a0a-7415-4d19-859d-a7e5c4648b56" ExactPredicates = "429591f6-91af-11e9-00e2-59fbe8cec110" +Preferences = "21216c6a-2e73-6563-6e65-726566657250" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" [compat] EnumX = "1.0" ExactPredicates = "2.2" -julia = "1" Random = "1" Test = "1" +julia = "1" +Preferences = "1.4" [extras] Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" diff --git a/docs/Manifest.toml b/docs/Manifest.toml index ac40fa4cd..973af2626 100644 --- a/docs/Manifest.toml +++ b/docs/Manifest.toml @@ -2,7 +2,7 @@ julia_version = "1.10.4" manifest_format = "2.0" -project_hash = "1618ed588410cd4dc664fc80ecbe5eaa4e31c80c" +project_hash = "a5f97bd51af051946996607d64905eb7a7b0c1f9" [[deps.ANSIColoredPrinters]] git-tree-sha1 = "574baf8110975760d391c710b6341da1afa48d8c" @@ -58,9 +58,9 @@ version = "1.1.1" [[deps.ArrayInterface]] deps = ["Adapt", "LinearAlgebra", "SparseArrays", "SuiteSparse"] -git-tree-sha1 = "ed2ec3c9b483842ae59cd273834e5b46206d6dda" +git-tree-sha1 = "5c9b74c973181571deb6442d41e5c902e6b9f38e" uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9" -version = "7.11.0" +version = "7.12.0" [deps.ArrayInterface.extensions] ArrayInterfaceBandedMatricesExt = "BandedMatrices" @@ -89,9 +89,9 @@ uuid = "56f22d72-fd6d-98f1-02f0-08ddc0907c33" [[deps.Automa]] deps = ["PrecompileTools", "TranscodingStreams"] -git-tree-sha1 = "588e0d680ad1d7201d4c6a804dcb1cd9cba79fbb" +git-tree-sha1 = "014bc22d6c400a7703c0f5dc1fdc302440cf88be" uuid = "67c07d97-cdcb-5c2c-af73-a7f9c32a568b" -version = "1.0.3" +version = "1.0.4" [[deps.AxisAlgorithms]] deps = ["LinearAlgebra", "Random", "SparseArrays", "WoodburyMatrices"] @@ -146,10 +146,10 @@ uuid = "159f3aea-2a34-519c-b102-8c37f9878175" version = "1.0.5" [[deps.CairoMakie]] -deps = ["CRC32c", "Cairo", "Colors", "FileIO", "FreeType", "GeometryBasics", "LinearAlgebra", "Makie", "PrecompileTools"] -git-tree-sha1 = "d69c7593fe9d7d617973adcbe4762028c6899b2c" +deps = ["CRC32c", "Cairo", "Cairo_jll", "Colors", "FileIO", "FreeType", "GeometryBasics", "LinearAlgebra", "Makie", "PrecompileTools"] +git-tree-sha1 = "f84837ccd1411ba059bb0b752dab9c7f1b0b0826" uuid = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" -version = "0.11.11" +version = "0.12.4" [[deps.Cairo_jll]] deps = ["Artifacts", "Bzip2_jll", "CompilerSupportLibraries_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"] @@ -175,9 +175,9 @@ weakdeps = ["SparseArrays"] [[deps.CodecZlib]] deps = ["TranscodingStreams", "Zlib_jll"] -git-tree-sha1 = "59939d8a997469ee05c4b4944560a820f9ba0d73" +git-tree-sha1 = "b8fe8546d52ca154ac556809e10c75e6e7430ac8" uuid = "944b1d66-785c-5afd-91f1-9de20f533193" -version = "0.7.4" +version = "0.7.5" [[deps.ColorBrewer]] deps = ["Colors", "JSON", "Test"] @@ -241,9 +241,9 @@ version = "1.1.1+0" [[deps.ConcurrentUtilities]] deps = ["Serialization", "Sockets"] -git-tree-sha1 = "6cbbd4d241d7e6579ab354737f4dd95ca43946e1" +git-tree-sha1 = "ea32b83ca4fefa1768dc84e504cc0a94fb1ab8d1" uuid = "f0e56b4a-5159-44fe-b623-3e5288b988bb" -version = "2.4.1" +version = "2.4.2" [[deps.ConstructionBase]] deps = ["LinearAlgebra"] @@ -292,10 +292,10 @@ uuid = "ab62b9b5-e342-54a8-a765-a90f495de1a6" version = "1.2.0" [[deps.DelaunayTriangulation]] -deps = ["EnumX", "ExactPredicates", "Random"] +deps = ["EnumX", "ExactPredicates", "Preferences", "Random"] path = ".." uuid = "927a84f5-c5f4-47a5-9785-b46e178433df" -version = "1.0.4" +version = "1.0.6" [[deps.DiffResults]] deps = ["StaticArraysCore"] @@ -348,9 +348,9 @@ version = "0.9.3" [[deps.Documenter]] deps = ["ANSIColoredPrinters", "AbstractTrees", "Base64", "CodecZlib", "Dates", "DocStringExtensions", "Downloads", "Git", "IOCapture", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "MarkdownAST", "Pkg", "PrecompileTools", "REPL", "RegistryInstances", "SHA", "TOML", "Test", "Unicode"] -git-tree-sha1 = "5461b2a67beb9089980e2f8f25145186b6d34f91" +git-tree-sha1 = "76deb8c15f37a3853f13ea2226b8f2577652de05" uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4" -version = "1.4.1" +version = "1.5.0" [[deps.DocumenterTools]] deps = ["AbstractTrees", "Base64", "DocStringExtensions", "Documenter", "FileWatching", "Gumbo", "LibGit2", "OpenSSH_jll", "Sass"] @@ -581,9 +581,9 @@ version = "1.3.14+0" [[deps.GridLayoutBase]] deps = ["GeometryBasics", "InteractiveUtils", "Observables"] -git-tree-sha1 = "6f93a83ca11346771a93bbde2bdad2f65b61498f" +git-tree-sha1 = "fc713f007cff99ff9e50accba6373624ddd33588" uuid = "3955a311-db13-416c-9275-1d80ed98e5e9" -version = "0.10.2" +version = "0.11.0" [[deps.Grisu]] git-tree-sha1 = "53bb909d1151e57e2484c3d1b53e19552b887fb2" @@ -692,13 +692,11 @@ deps = ["Adapt", "AxisAlgorithms", "ChainRulesCore", "LinearAlgebra", "OffsetArr git-tree-sha1 = "88a101217d7cb38a7b481ccd50d21876e1d1b0e0" uuid = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" version = "0.15.1" +weakdeps = ["Unitful"] [deps.Interpolations.extensions] InterpolationsUnitfulExt = "Unitful" - [deps.Interpolations.weakdeps] - Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" - [[deps.IntervalArithmetic]] deps = ["CRlibm_jll", "MacroTools", "RoundingEmulator"] git-tree-sha1 = "433b0bb201cd76cb087b017e49244f10394ebe9c" @@ -950,16 +948,16 @@ uuid = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09" version = "0.5.13" [[deps.Makie]] -deps = ["Animations", "Base64", "CRC32c", "ColorBrewer", "ColorSchemes", "ColorTypes", "Colors", "Contour", "DelaunayTriangulation", "Distributions", "DocStringExtensions", "Downloads", "FFMPEG_jll", "FileIO", "FilePaths", "FixedPointNumbers", "Format", "FreeType", "FreeTypeAbstraction", "GeometryBasics", "GridLayoutBase", "ImageIO", "InteractiveUtils", "IntervalSets", "Isoband", "KernelDensity", "LaTeXStrings", "LinearAlgebra", "MacroTools", "MakieCore", "Markdown", "MathTeXEngine", "Observables", "OffsetArrays", "Packing", "PlotUtils", "PolygonOps", "PrecompileTools", "Printf", "REPL", "Random", "RelocatableFolders", "Scratch", "ShaderAbstractions", "Showoff", "SignedDistanceFields", "SparseArrays", "Statistics", "StatsBase", "StatsFuns", "StructArrays", "TriplotBase", "UnicodeFun"] -git-tree-sha1 = "4d49c9ee830eec99d3e8de2425ff433ece7cc1bc" +deps = ["Animations", "Base64", "CRC32c", "ColorBrewer", "ColorSchemes", "ColorTypes", "Colors", "Contour", "Dates", "DelaunayTriangulation", "Distributions", "DocStringExtensions", "Downloads", "FFMPEG_jll", "FileIO", "FilePaths", "FixedPointNumbers", "Format", "FreeType", "FreeTypeAbstraction", "GeometryBasics", "GridLayoutBase", "ImageIO", "InteractiveUtils", "IntervalSets", "Isoband", "KernelDensity", "LaTeXStrings", "LinearAlgebra", "MacroTools", "MakieCore", "Markdown", "MathTeXEngine", "Observables", "OffsetArrays", "Packing", "PlotUtils", "PolygonOps", "PrecompileTools", "Printf", "REPL", "Random", "RelocatableFolders", "Scratch", "ShaderAbstractions", "Showoff", "SignedDistanceFields", "SparseArrays", "Statistics", "StatsBase", "StatsFuns", "StructArrays", "TriplotBase", "UnicodeFun", "Unitful"] +git-tree-sha1 = "57a1a2b3d12e04f9e9fb77d61cd12571d5541c5f" uuid = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" -version = "0.20.10" +version = "0.21.4" [[deps.MakieCore]] -deps = ["Observables", "REPL"] -git-tree-sha1 = "248b7a4be0f92b497f7a331aed02c1e9a878f46b" +deps = ["ColorTypes", "GeometryBasics", "IntervalSets", "Observables"] +git-tree-sha1 = "638bc817096742e8302f7b0b972ee5701fe00e97" uuid = "20f20a25-4f0e-4fdf-b5d1-57303727442b" -version = "0.7.3" +version = "0.8.3" [[deps.MappedArrays]] git-tree-sha1 = "2dab0221fe2b0f2cb6754eaa743cc266339f527e" @@ -978,9 +976,9 @@ version = "0.1.2" [[deps.MathTeXEngine]] deps = ["AbstractTrees", "Automa", "DataStructures", "FreeTypeAbstraction", "GeometryBasics", "LaTeXStrings", "REPL", "RelocatableFolders", "UnicodeFun"] -git-tree-sha1 = "96ca8a313eb6437db5ffe946c457a401bbb8ce1d" +git-tree-sha1 = "1865d0b8a2d91477c8b16b49152a32764c7b1f5f" uuid = "0a4f8689-d25c-4efe-a92b-7142dfc1aa53" -version = "0.5.7" +version = "0.6.0" [[deps.MbedTLS]] deps = ["Dates", "MbedTLS_jll", "MozillaCACerts_jll", "NetworkOptions", "Random", "Sockets"] @@ -1050,9 +1048,9 @@ uuid = "510215fc-4207-5dde-b226-833fc4488ee2" version = "0.5.5" [[deps.OffsetArrays]] -git-tree-sha1 = "e64b4f5ea6b7389f6f046d13d4896a8f9c1ba71e" +git-tree-sha1 = "1a27764e945a152f7ca7efa04de513d473e9542e" uuid = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" -version = "1.14.0" +version = "1.14.1" weakdeps = ["Adapt"] [deps.OffsetArrays.extensions] @@ -1500,9 +1498,9 @@ version = "0.1.1" [[deps.StaticArrays]] deps = ["LinearAlgebra", "PrecompileTools", "Random", "StaticArraysCore"] -git-tree-sha1 = "20833c5b7f7edf0e5026f23db7f268e4f23ec577" +git-tree-sha1 = "eeafab08ae20c62c44c8399ccb9354a04b80db50" uuid = "90137ffa-7385-5640-81b9-e52037218182" -version = "1.9.6" +version = "1.9.7" weakdeps = ["ChainRulesCore", "Statistics"] [deps.StaticArrays.extensions] @@ -1611,9 +1609,9 @@ uuid = "731e570b-9d59-4bfa-96dc-6df516fadf69" version = "0.10.0" [[deps.TranscodingStreams]] -git-tree-sha1 = "d73336d81cafdc277ff45558bb7eaa2b04a8e472" +git-tree-sha1 = "60df3f8126263c0d6b357b9a1017bb94f53e3582" uuid = "3bb67fe8-82b1-5028-8e26-92a6c54297fa" -version = "0.10.10" +version = "0.11.0" weakdeps = ["Random", "Test"] [deps.TranscodingStreams.extensions] @@ -1647,6 +1645,20 @@ git-tree-sha1 = "53915e50200959667e78a92a418594b428dffddf" uuid = "1cfade01-22cf-5700-b092-accc4b62d6e1" version = "0.4.1" +[[deps.Unitful]] +deps = ["Dates", "LinearAlgebra", "Random"] +git-tree-sha1 = "dd260903fdabea27d9b6021689b3cd5401a57748" +uuid = "1986cc42-f94f-5a68-af5c-568840ba703d" +version = "1.20.0" + + [deps.Unitful.extensions] + ConstructionBaseUnitfulExt = "ConstructionBase" + InverseFunctionsUnitfulExt = "InverseFunctions" + + [deps.Unitful.weakdeps] + ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9" + InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112" + [[deps.WoodburyMatrices]] deps = ["LinearAlgebra", "SparseArrays"] git-tree-sha1 = "c1a7aa6219628fcd757dede0ca95e245c5cd9511" @@ -1760,10 +1772,10 @@ uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f" version = "1.6.43+1" [[deps.libsass_jll]] -deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"] -git-tree-sha1 = "941afb93587dcec07f89e511057f5efc0bec6f0d" +deps = ["Artifacts", "JLLWrappers", "Libdl"] +git-tree-sha1 = "6044ffe7e7bf0602e2039dc747c3332a097ac74b" uuid = "47bcb7c8-5119-555a-9eeb-0afcc36cd728" -version = "3.6.4+0" +version = "3.6.6+0" [[deps.libsixel_jll]] deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Pkg", "libpng_jll"] diff --git a/docs/Project.toml b/docs/Project.toml index bd7a5bf9b..9809d0605 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -10,6 +10,7 @@ InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589" +Preferences = "21216c6a-2e73-6563-6e65-726566657250" ReferenceTests = "324d217c-45ce-50fc-942e-d289b448e8cf" SimpleGraphs = "55797a34-41de-5266-9ec1-32ac4eb504d3" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" diff --git a/docs/src/literate_tutorials/convex.jl b/docs/src/literate_tutorials/convex.jl index 9f59b36a9..edffffb84 100644 --- a/docs/src/literate_tutorials/convex.jl +++ b/docs/src/literate_tutorials/convex.jl @@ -8,6 +8,7 @@ using DelaunayTriangulation using CairoMakie using ReferenceTests #src using Test #src +using Preferences #src fig_path = joinpath(@__DIR__, "../figures") #src points = [ @@ -21,7 +22,7 @@ tri = triangulate_convex(points, 1:7) #- fig, ax, sc = triplot(tri) fig -@test_reference joinpath(fig_path, "convex_ex_1.png") fig #src +load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true) && @test_reference joinpath(fig_path, "convex_ex_1.png") fig #src # This `tri` is our triangulation of the convex polygon. # The first input is the set of points, and `S` defines @@ -44,7 +45,7 @@ tri = triangulate_convex(points, S) #- fig, ax, sc = triplot(tri) fig -@test_reference joinpath(fig_path, "convex_ex_2.png") fig #src +load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true) && @test_reference joinpath(fig_path, "convex_ex_2.png") fig #src # Here is a comparison of the time it takes to triangulate this # using `triangulate_convex` or `triangulate`. diff --git a/docs/src/literate_tutorials/point_location.jl b/docs/src/literate_tutorials/point_location.jl index 2bb0ecc51..89b0caca1 100644 --- a/docs/src/literate_tutorials/point_location.jl +++ b/docs/src/literate_tutorials/point_location.jl @@ -21,6 +21,7 @@ using StableRNGs using Test #src using ReferenceTests #src fig_path = joinpath(@__DIR__, "../figures") #src +using Preferences #src points = [ (-3.0, 6.0), (5.0, 1.0), (-5.0, 3.0), (2.0, -3.0), @@ -78,7 +79,7 @@ V = jump_and_march(tri, q) fig, ax, sc = triplot(tri, show_ghost_edges=true) scatter!(ax, q) fig -@test_reference joinpath(fig_path, "point_location_ex_1.png") fig #src +load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true) && @test_reference joinpath(fig_path, "point_location_ex_1.png") fig #src # ## Region with concave boundaries and holes # Now we give an example of point location for a reason with holes. Since the @@ -104,7 +105,7 @@ refine!(tri; max_area=0.01get_area(tri), rng); # To demonstrate this, see the following plot: fig, ax, sc = triplot(tri, show_ghost_edges=true) fig -@test_reference joinpath(fig_path, "point_location_ex_2.png") fig #src +load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true) && @test_reference joinpath(fig_path, "point_location_ex_2.png") fig #src # The ghost edges now intersect the boundary, which doesn't make sense, and creates difficulties. # Let us now demonstrate how the function still works here. We try finding the blue points shown below. @@ -117,7 +118,7 @@ qs = [ fig, ax, sc = triplot(tri, show_ghost_edges=false) scatter!(ax, qs, color=:blue, markersize=16) fig -@test_reference joinpath(fig_path, "point_location_ex_3.png") fig #src +load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true) && @test_reference joinpath(fig_path, "point_location_ex_3.png") fig #src # Now let's find the triangles. Vs = [jump_and_march(tri, q; rng) for q in qs] @@ -181,7 +182,7 @@ qs = [ fig, ax, sc = triplot(tri) scatter!(ax, qs, color=:blue, markersize=16) fig -@test_reference joinpath(fig_path, "point_location_ex_4.png") fig by=psnr_equality(10) #src +load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true) && @test_reference joinpath(fig_path, "point_location_ex_4.png") fig by=psnr_equality(10) #src # Here are the `jump_and_march` results. Vs = [jump_and_march(tri, q; rng, concavity_protection=true) for q in qs] diff --git a/src/DelaunayTriangulation.jl b/src/DelaunayTriangulation.jl index 4e2f6d752..ccc052e97 100644 --- a/src/DelaunayTriangulation.jl +++ b/src/DelaunayTriangulation.jl @@ -1,9 +1,43 @@ module DelaunayTriangulation +""" + DefaultAdjacentValue = 0 + +Default value used for representing an empty result +from an adjacency query. +""" const DefaultAdjacentValue = 0 + +""" + ∅ = DefaultAdjacentValue + +Alias for [`DefaultAdjacentValue`](@ref). +""" const ∅ = DefaultAdjacentValue + +""" + GhostVertex = -1 + +Number used for representing initial ghost vertices. +All other ghost vertices are derived from subtracting from +this number. See https://juliageometry.github.io/DelaunayTriangulation.jl/stable/manual/ghost_triangles/. +""" const GhostVertex = -1 + +""" + 𝒢 = GhostVertex + +Alias for [`𝒢`](@ref). +""" const 𝒢 = GhostVertex + +""" + ε = sqrt(eps(Float64)) + +Number used as a tolerance in certain functions, e.g. +for mesh refinement when using [`check_precision`](@ref) to +avoid degenerate circumcenters. +""" const ε = sqrt(eps(Float64)) const INF_WARN = Ref(true) @@ -15,6 +49,42 @@ By default, this warning is enabled. """ toggle_inf_warn!() = (INF_WARN[] = !INF_WARN[]) +@static if VERSION ≥ v"1.6" + using Preferences +end + +@static if VERSION ≥ v"1.6" + const USE_EXACTPREDICATES = @load_preference("USE_EXACTPREDICATES", true)::Bool +else + const USE_EXACTPREDICATES = true +end +@doc """ + USE_EXACTPREDICATES + +Whether to use ExactPredicates.jl for computing predicates. By default, +this is true, but a user can change this by defining a preference with Preferences.jl, i.e. +you could do the following + +```julia-repl +julia> using Preferences: set_preferences! + +julia> set_preferences!("DelaunayTriangulation", "USE_EXACTPREDICATES" => false) + +julia> using DelaunayTriangulation # load only after setting the preference + +julia> DelaunayTriangulation.USE_EXACTPREDICATES +false +``` + +You have set `USE_EXACTPREDICATES = $USE_EXACTPREDICATES`. + +!!! note "Precision" + + Even if you have disabled ExactPredicates.jl, the predicates + are still computed in Float64 precision. +""" +USE_EXACTPREDICATES + using ExactPredicates using EnumX using Random diff --git a/src/geometric_primitives/points.jl b/src/geometric_primitives/points.jl index dbc541035..d9fd72237 100644 --- a/src/geometric_primitives/points.jl +++ b/src/geometric_primitives/points.jl @@ -32,6 +32,13 @@ julia> gety(p) """ gety(p) = p[2] +""" + getz(p) -> Number + +Get the z-coordinate of `p`. +""" +getz(p) = p[3] + """ getxy(p) -> NTuple{2, Number} @@ -49,6 +56,13 @@ julia> getxy(p) """ getxy(p) = (getx(p), gety(p)) +""" + getxyz(p) -> NTuple{3, Number} + +Given a three-dimensional `p`, returns its coordinates as a `Tuple`. +""" +getxyz(p) = (getx(p), gety(p), getz(p)) + """ _getx(p) -> Float64 @@ -93,10 +107,17 @@ julia> DelaunayTriangulation._gety(p) """ _gety(p) = Float64(gety(p)) +""" + _getz(p) -> Float64 + +Get the z-coordinate of `p` as a `Float64.` +""" +_getz(p) = Float64(getz(p)) + """ _getxy(p) -> NTuple{2, Float64} -Get the coordinates of `p` as a `Tuple`. +Get the coordinates of `p` as a `Tuple` of `Float64`s. # Examples ```jldoctest @@ -115,6 +136,13 @@ julia> DelaunayTriangulation._getxy(p) """ _getxy(p) = (_getx(p), _gety(p)) +""" + _getxyz(p) -> NTuple{3, Float64} + +Given a three-dimemsional `p`, returns its coordinates as a `Tuple` of `Float64`s. +""" +_getxyz(p) = (_getx(p), _gety(p), _getz(p)) + @doc """ getpoint(points, vertex) -> NTuple{2, Number} diff --git a/src/predicates/exactpredicates_definitions.jl b/src/predicates/exactpredicates_definitions.jl index f699c77fa..496ea5987 100644 --- a/src/predicates/exactpredicates_definitions.jl +++ b/src/predicates/exactpredicates_definitions.jl @@ -1,3 +1,78 @@ +#= +Much of the code in the file is derived or taken from ExactPredicates.jl (https://github.com/lairez/ExactPredicates.jl). +The license for ExactPredicates.jl is printed below, taken from https://github.com/lairez/ExactPredicates.jl/blob/master/LICENSE. + +Copyright (c) 2019 Pierre Lairez + +Permission is hereby granted, free of charge, to any person obtaining a copy +of this software and associated documentation files (the "Software"), to deal +in the Software without restriction, including without limitation the rights +to use, copy, modify, merge, publish, distribute, sublicense, and/or sell +copies of the Software, and to permit persons to whom the Software is +furnished to do so, subject to the following conditions: + +The above copyright notice and this permission notice shall be included in all +copies or substantial portions of the Software. + +THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR +IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, +FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE +AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER +LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, +OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE +SOFTWARE. +=# + +""" + ext(ux, uy, vx, vy) + +Computes `ExactPredicates.ext((ux, uy), (vx, vy))`, i.e. +returns `ux * vy - uy * vx`. +""" +ext(ux, uy, vx, vy) = ux * vy - uy * vx + +""" + inp(ux, uy, vx, vy) + +Computes `ExactPredicates.inp((ux, uy), (vx, vy))`, i.e. +returns `ux * vx + uy * vy`. +""" +inp(ux, uy, vx, vy) = ux * vx + uy * vy + +""" + det(a, b, c, d) + +Computes `ExactPredicates.det(a, b, c, d)`, +i.e. returns `a*d - b*c`. +""" +det(a, b, c, d) = a * d - b * c + +""" + det(a, b, c, d, e, f, g, h, i) + +Computes `ExactPredicates.det(a, b, c, d, e, f, g, h, i)`, +i.e. returns the determinant of +```math +\\det \\begin{bmatrix} a & b & c \\\\ d & e & f \\\\ g & h & i \\end{bmatrix} +``` +""" +det(a, b, c, d, e, f, g, h, i) = a * det(e, f, h, i) - d * det(b, c, h, i) + g * det(b, c, e, f) + +""" + opposite_signs(x, y) -> Bool + +From ExactPredicates.jl, returns `true` if `x` and `y` have opposite signs, and `false` otherwise. +Assumes that `x` and `y` are in `[-1, 0, 1]`. +""" +opposite_signs(x, y) = xor(x, y) == -2 + +""" + sgn(x) -> Int + +Returns `Int(sign(x))`. +""" +sgn(x) = Int(sign(x)) + """ orient_predicate(p, q, r) -> Integer @@ -7,6 +82,10 @@ Returns `ExactPredicates.orient(p, q, r)`, in particular we return: - `0`: `(p, q, r)` is collinear / degenerate. - `-1`: `(p, q, r)` is negatively oriented. +If ExactPredicates.jl has been disabled using Preferences.jl, the +determinant defining `orient` is evaluated numerically +without exact arithmetic. + # Extended help The orient predicate is defined by the determinant @@ -14,7 +93,21 @@ The orient predicate is defined by the determinant \\text{orient}(p, q, r) = \\text{sgn} \\det \\begin{bmatrix} p_x & p_y & 1 \\\\ q_x & q_y & 1 \\\\ r_x & r_y & 1 \\end{bmatrix} = \\text{sgn} \\det \\begin{bmatrix} p_x-r_x & p_y-r_y \\\\ q_x-r_x & q_y-r_y \\end{bmatrix}. ``` """ -orient_predicate(p, q, r) = orient(_getxy(p), _getxy(q), _getxy(r)) +orient_predicate(::Any, ::Any, ::Any) + +@static if USE_EXACTPREDICATES + orient_predicate(p, q, r) = orient(_getxy(p), _getxy(q), _getxy(r)) +else + orient_predicate(p, q, r) = _orient_predicate(_getxy(p), _getxy(q), _getxy(r)) +end +function _orient_predicate(p, q, r) + px, py = _getxy(p) + qx, qy = _getxy(q) + rx, ry = _getxy(r) + ux, uy = px - rx, py - ry + vx, vy = qx - rx, qy - ry + return sgn(ext(ux, uy, vx, vy)) +end """ orient_predicate(p, q, r, s) -> Integer @@ -44,14 +137,34 @@ Here, a positively oriented tetrahedron `(p, q, r, s)` takes the form 'r '\\. +If ExactPredicates.jl has been disabled using Preferences.jl, the +determinant defining `orient` is evaluated numerically +without exact arithmetic. + # Extended help The orient predicate is defined by the determinant ```math -\\text{orient}(p, q, r, s) = \\text{sgn} \\det \\begin{bmatrix} p_x & p_y & p_y & 1 \\\\ q_x & q_y & q_r & 1 \\\\ r_x & r_y & r_z & 1 \\\\ s_x & s_y & s_z & 1 \\end{bmatrix} = \\text{sgn} \\det \\begin{bmatrix} p_x - s_x & p_y - s_y & p_z - s_y \\\\ q_x - s_x & q_y - s_y & q_z - s_z \\\\ r_x - s_x & r_y - s_y & r_z - s_z \\end{bmatrix}. +\\text{orient}(p, q, r, s) = \\text{sgn} \\det \\begin{bmatrix} p_x & p_y & p_z & 1 \\\\ q_x & q_y & q_z & 1 \\\\ r_x & r_y & r_z & 1 \\\\ s_x & s_y & s_z & 1 \\end{bmatrix} = \\text{sgn} \\det \\begin{bmatrix} p_x - s_x & p_y - s_y & p_z - s_y \\\\ q_x - s_x & q_y - s_y & q_z - s_z \\\\ r_x - s_x & r_y - s_y & r_z - s_z \\end{bmatrix}. ``` """ -orient_predicate(p, q, r, s) = orient(p, q, r, s) +orient_predicate(::Any, ::Any, ::Any, ::Any) + +@static if USE_EXACTPREDICATES + orient_predicate(p, q, r, s) = orient(_getxyz(p), _getxyz(q), _getxyz(r), _getxyz(s)) +else + orient_predicate(p, q, r, s) = _orient_predicate(_getxyz(p), _getxyz(q), _getxyz(r), _getxyz(s)) +end +function _orient_predicate(p, q, r, s) + px, py, pz = _getxyz(p) + qx, qy, qz = _getxyz(q) + rx, ry, rz = _getxyz(r) + sx, sy, sz = _getxyz(s) + a, b, c = px - sx, py - sy, pz - sz + d, e, f = qx - sx, qy - sy, qz - sz + g, h, i = rx - sx, ry - sy, rz - sz + return sgn(det(a, b, c, d, e, f, g, h, i)) +end """ incircle_predicate(a, b, c, p) -> Integer @@ -62,6 +175,10 @@ Returns `ExactPredicates.incircle(a, b, c, p)`, in particular we return: - `0`: If `p` is on the circle defined by `(a, b, c)`. - `-1`: If `p` is outside the circle defined by `(a, b, c)`. +If ExactPredicates.jl has been disabled using Preferences.jl, the +determinant defining `incircle` is evaluated numerically +without exact arithmetic. + # Extended help The incircle predicate is defined by the determinant @@ -69,7 +186,27 @@ The incircle predicate is defined by the determinant \\text{incircle}(a, b, c, d) = \\text{sgn} \\det \\begin{bmatrix} a_x & a_y & a_x^2 + a_y^2 & 1 \\\\ b_x & b_y & b_x62 + b_y^2 & 1 \\\\ c_x & c_y & c_x^2 + c_y^2 & 1 \\\\ d_x & d_y & d_x^2 + d_y^2 & 1 \\end{bmatrix} = \\text{sgn} \\det \\begin{bmatrix} a_x - d_x & a_y - d_y & (a_x - d_x)^2 + (a_y - d_y)^2 \\\\ b_x - d_x & b_y - d_y & (b_x - d_x)^2 + (b_y - d_y)^2 \\\\ c_x - d_x & c_y - d_y & (c_x - d_x)^2 + (c_y - d_y)^2 \\end{bmatrix}. ``` """ -incircle_predicate(a, b, c, p) = incircle(_getxy(a), _getxy(b), _getxy(c), _getxy(p)) +incircle_predicate + +@static if USE_EXACTPREDICATES + incircle_predicate(a, b, c, p) = incircle(_getxy(a), _getxy(b), _getxy(c), _getxy(p)) +else + incircle_predicate(a, b, c, p) = _incircle_predicate(_getxy(a), _getxy(b), _getxy(c), _getxy(p)) +end +function _incircle_predicate(p, q, r, a) + px, py = _getxy(p) + qx, qy = _getxy(q) + rx, ry = _getxy(r) + ax, ay = _getxy(a) + qpx, qpy = qx - px, qy - py + rpx, rpy = rx - px, ry - py + apx, apy = ax - px, ay - py + aqx, aqy = ax - qx, ay - qy + rqx, rqy = rx - qx, ry - qy + val1 = ext(qpx, qpy, apx, apy) * inp(rpx, rpy, rqx, rqy) + val2 = ext(qpx, qpy, rpx, rpy) * inp(apx, apy, aqx, aqy) + return sgn(val1 - val2) +end """ parallelorder_predicate(a, b, p, q) -> Integer @@ -79,8 +216,27 @@ Returns `ExactPredicates.parallelorder(a, b, p, q)`, in particular we return: - `1`: `q` is closer to the line `(a, b)` than `p`. - `0`: `p` and `q` are equidistant from the line `(a, b)`. - `-1`: `p` is closer to the line `(a, b)` than `q`. + +If ExactPredicates.jl has been disabled using Preferences.jl, the +the predicates defining `parallelorder` are all evaluated numerically +without exact arithmetic. """ -parallelorder_predicate(a, b, p, q) = parallelorder(_getxy(a), _getxy(b), _getxy(p), _getxy(q)) +parallelorder_predicate + +@static if USE_EXACTPREDICATES + parallelorder_predicate(a, b, p, q) = parallelorder(_getxy(a), _getxy(b), _getxy(p), _getxy(q)) +else + parallelorder_predicate(a, b, p, q) = _parallelorder_predicate(_getxy(a), _getxy(b), _getxy(p), _getxy(q)) +end +function _parallelorder_predicate(a, b, p, q) + ax, ay = _getxy(a) + bx, by = _getxy(b) + px, py = _getxy(p) + qx, qy = _getxy(q) + δx, δy = bx - ax, by - ay + qpx, qpy = qx - px, qy - py + return sgn(ext(δx, δy, qpx, qpy)) +end """ sameside_predicate(a, b, p) -> Integer @@ -109,13 +265,6 @@ function sameside_predicate(a, b, p) end end -""" - opposite_signs(x, y) -> Bool - -From ExactPredicates.jl, returns `true` if `x` and `y` have opposite signs, and `false` otherwise. -""" -opposite_signs(x, y) = xor(x, y) == -2 - """ meet_predicate(p, q, a, b) @@ -124,6 +273,10 @@ Returns `ExactPredicates.meet(p, q, a, b)`, in particular we return: - `1`: The open line segments `(p, q)` and `(a, b)` meet in a single point. - `0`: The closed line segments `[p, q]` and `[a, b]` meet in one or several points. - `-1`: Otherwise. + +If ExactPredicates.jl has been disabled using Preferences.jl, then all +the `orient` predicates computed for this result are evaluated numerically +without exact arithmetic. """ function meet_predicate(p, q, a, b) pqa = orient_predicate(p, q, a) @@ -146,4 +299,35 @@ function meet_predicate(p, q, a, b) else return 0 end -end \ No newline at end of file +end + +@static if USE_EXACTPREDICATES + ExactPredicates.Codegen.@genpredicate function angle_is_acute(p::2, q::2, r::2) + pr = p - r + qr = q - r + ExactPredicates.Codegen.group!(pr...) + ExactPredicates.Codegen.group!(qr...) + return pr[1] * qr[1] + pr[2] * qr[2] + end +else + function angle_is_acute(p, q, r) + px, py = _getxy(p) + qx, qy = _getxy(q) + rx, ry = _getxy(r) + ux, uy = px - rx, py - ry + vx, vy = qx - rx, qy - ry + return sgn(inp(ux, uy, vx, vy)) + end +end + +@doc """ + angle_is_acute(p, q, r) + +Tests if the angle opposite `(p, q)` in the triangle `(p, q, r)`, +meaning `∠prq`, is acute, returning: + +- `1`: `∠prq` is acute. +- `0`: `∠prq` is a right angle. +- `-1`: `∠prq` is obtuse. +""" +angle_is_acute \ No newline at end of file diff --git a/src/predicates/predicates.jl b/src/predicates/predicates.jl index d004af529..d3351db2d 100644 --- a/src/predicates/predicates.jl +++ b/src/predicates/predicates.jl @@ -540,14 +540,6 @@ function point_position_relative_to_witness_plane(tri::Triangulation, i, j, k, return convert_certificate(cert, Cert.Above, Cert.On, Cert.Below) end -ExactPredicates.Codegen.@genpredicate function angle_is_acute(p::2, q::2, r::2) - pr = p - r - qr = q - r - ExactPredicates.Codegen.group!(pr...) - ExactPredicates.Codegen.group!(qr...) - return pr[1] * qr[1] + pr[2] * qr[2] -end - """ opposite_angle(p, q, r) -> Certificate diff --git a/test/Project.toml b/test/Project.toml index a75aeaaec..496f36e71 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -5,6 +5,7 @@ CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" +DelaunayTriangulation = "927a84f5-c5f4-47a5-9785-b46e178433df" DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8" @@ -17,6 +18,8 @@ GeometryBasics = "5c1252a2-5f33-56bf-86c9-59e7332b4326" InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" OrderedCollections = "bac558e1-5e72-5ebc-8fee-abe8a469f55d" +Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f" +Preferences = "21216c6a-2e73-6563-6e65-726566657250" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" ReferenceTests = "324d217c-45ce-50fc-942e-d289b448e8cf" SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f" diff --git a/test/interfaces/points.jl b/test/interfaces/points.jl index 696f5580f..fc609e339 100644 --- a/test/interfaces/points.jl +++ b/test/interfaces/points.jl @@ -1,7 +1,7 @@ using ..DelaunayTriangulation const DT = DelaunayTriangulation using Test -using StaticArraysCore +using StaticArrays using StatsBase import GeometryBasics: Point2f @@ -235,3 +235,35 @@ end ) end +@testset "getz/_getz/getxyz/_getxyz" begin + p = (1.0, 2.0, 3.0) + q = (1.0f0, 2.0f0, 3.0f0) + r = [1.0f0, 4.0f0, 5.0f0] + s = @SVector [3, 7, 5] + t = [5.0, 13.0, -5.0] + + @test DT.getz(p) === 3.0 + @test DT._getz(p) === 3.0 + @test DT.getxyz(p) === (1.0, 2.0, 3.0) + @test DT._getxyz(p) === (1.0, 2.0, 3.0) + + @test DT.getz(q) === 3.0f0 + @test DT._getz(q) === 3.0 + @test DT.getxyz(q) === (1.0f0, 2.0f0, 3.0f0) + @test DT._getxyz(q) === (1.0, 2.0, 3.0) + + @test DT.getz(r) === 5.0f0 + @test DT._getz(r) === 5.0 + @test DT.getxyz(r) === (1.0f0, 4.0f0, 5.0f0) + @test DT._getxyz(r) === (1.0, 4.0, 5.0) + + @test DT.getz(s) === 5 + @test DT._getz(s) === 5.0 + @test DT.getxyz(s) === (3, 7, 5) + @test DT._getxyz(s) === (3.0, 7.0, 5.0) + + @test DT.getz(t) === -5.0 + @test DT._getz(t) === -5.0 + @test DT.getxyz(t) === (5.0, 13.0, -5.0) + @test DT._getxyz(t) === (5.0, 13.0, -5.0) +end \ No newline at end of file diff --git a/test/operations/delete_ghost_triangles.jl b/test/operations/delete_ghost_triangles.jl index b57d6804c..0b4d963ff 100644 --- a/test/operations/delete_ghost_triangles.jl +++ b/test/operations/delete_ghost_triangles.jl @@ -3,7 +3,7 @@ const DT = DelaunayTriangulation using StructEquality using ..DelaunayTriangulation: Triangulation using StaticArrays - +using Preferences @testset "Deleting ghost triangles" begin tri, label_map, index_map = simple_geometry() @@ -15,15 +15,19 @@ using StaticArrays end @testset "Making sure the ghost triangles are correctly adjusted on disjoint sets" begin - θ = LinRange(0, 2π, 500) |> collect + θ = LinRange(0, 2π, 100) |> collect θ[end] = 0 xy = Vector{Vector{Vector{NTuple{2,Float64}}}}() cx = 0.0 for i in 1:2 - push!(xy, [[(cx + cos(θ), sin(θ)) for θ in θ]]) - push!(xy, [[(cx + 0.5cos(θ), 0.5sin(θ)) for θ in reverse(θ)]]) + push!(xy, [[(cx + cos(θ) + 1e-3rand(), sin(θ) + 1e-3rand()) for θ in θ]]) # use perturbations to also do this check without ExactPredicates loaded. iszero makes the boundary closed + push!(xy, [[(cx + 0.5cos(θ) + 1e-3rand(), 0.5sin(θ) + 1e-3rand()) for θ in reverse(θ)]]) cx += 3.0 end + xy[1][1][end] = xy[1][1][1] + xy[2][1][end] = xy[2][1][1] + xy[3][1][end] = xy[3][1][1] + xy[4][1][end] = xy[4][1][1] boundary_nodes, points = convert_boundary_points_to_indices(xy) tri = triangulate(points; boundary_nodes=boundary_nodes) @test validate_triangulation(tri) diff --git a/test/operations/delete_holes.jl b/test/operations/delete_holes.jl index 05bb77df2..8c965db82 100644 --- a/test/operations/delete_holes.jl +++ b/test/operations/delete_holes.jl @@ -3,7 +3,7 @@ const DT = DelaunayTriangulation using CairoMakie using StableRNGs using StaticArrays - +using Preferences @testset "Simple domain" begin p1 = (0.2, 0.2) @@ -142,68 +142,70 @@ end validate_statistics(tri) end -@testset "Interior holes that were already triangles" begin - p1 = (0.0, 0.0) - p2 = (1.0, 0.0) - p3 = (1.0, 1.0) - p4 = (0.0, 1.0) - p5 = (0.2, 0.2) - p6 = (0.8, 0.2) - p7 = (0.8, 0.8) - pts = [p1, p2, p3, p4, p5, p6, p7] - tri = triangulate(pts, boundary_nodes=[[[1, 2, 3, 4, 1]], [[5, 7, 6, 5]]], delete_holes=false, delete_ghosts=false) - @test num_ghost_triangles(tri) == 4 - @test num_ghost_edges(tri) == 4 - points_to_process = DT.find_all_points_to_delete(tri) - @test points_to_process == Set((DT.𝒢,)) - triangles_to_delete = DT.find_all_triangles_to_delete(tri, points_to_process) - @test DT.compare_triangle_collections(triangles_to_delete, Set(( - (1, 4, DT.𝒢), - (4, 3, DT.𝒢), - (3, 2, DT.𝒢), - (2, 1, DT.𝒢), - (5, 6, 7) - ))) - tri = triangulate(pts; boundary_nodes=[[[1, 2, 3, 4, 1]], [[5, 7, 6, 5]]]) - @test validate_triangulation(tri) - validate_statistics(tri) +if load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true) + @testset "Interior holes that were already triangles" begin + p1 = (0.0, 0.0) + p2 = (1.0, 0.0) + p3 = (1.0, 1.0) + p4 = (0.0, 1.0) + p5 = (0.2, 0.2) + p6 = (0.8, 0.2) + p7 = (0.8, 0.8) + pts = [p1, p2, p3, p4, p5, p6, p7] + tri = triangulate(pts, boundary_nodes=[[[1, 2, 3, 4, 1]], [[5, 7, 6, 5]]], delete_holes=false, delete_ghosts=false) + @test num_ghost_triangles(tri) == 4 + @test num_ghost_edges(tri) == 4 + points_to_process = DT.find_all_points_to_delete(tri) + @test points_to_process == Set((DT.𝒢,)) + triangles_to_delete = DT.find_all_triangles_to_delete(tri, points_to_process) + @test DT.compare_triangle_collections(triangles_to_delete, Set(( + (1, 4, DT.𝒢), + (4, 3, DT.𝒢), + (3, 2, DT.𝒢), + (2, 1, DT.𝒢), + (5, 6, 7) + ))) + tri = triangulate(pts; boundary_nodes=[[[1, 2, 3, 4, 1]], [[5, 7, 6, 5]]]) + @test validate_triangulation(tri) + validate_statistics(tri) - tri = triangulate(pts, boundary_nodes=[[[1, 2, 3, 4, 1]], [[5, 7], [7, 6, 5]]], delete_holes=false, delete_ghosts=false) - points_to_process = DT.find_all_points_to_delete(tri) - @test points_to_process == Set((DT.𝒢,)) - triangles_to_delete = DT.find_all_triangles_to_delete(tri, points_to_process) - @test DT.compare_triangle_collections(triangles_to_delete, Set(( - (1, 4, DT.𝒢), - (4, 3, DT.𝒢), - (3, 2, DT.𝒢), - (2, 1, DT.𝒢), - (5, 6, 7) - ))) - tri = triangulate(pts; boundary_nodes=[[[1, 2, 3, 4, 1]], [[5, 7], [7, 6, 5]]]) - @test validate_triangulation(tri) - validate_statistics(tri) + tri = triangulate(pts, boundary_nodes=[[[1, 2, 3, 4, 1]], [[5, 7], [7, 6, 5]]], delete_holes=false, delete_ghosts=false) + points_to_process = DT.find_all_points_to_delete(tri) + @test points_to_process == Set((DT.𝒢,)) + triangles_to_delete = DT.find_all_triangles_to_delete(tri, points_to_process) + @test DT.compare_triangle_collections(triangles_to_delete, Set(( + (1, 4, DT.𝒢), + (4, 3, DT.𝒢), + (3, 2, DT.𝒢), + (2, 1, DT.𝒢), + (5, 6, 7) + ))) + tri = triangulate(pts; boundary_nodes=[[[1, 2, 3, 4, 1]], [[5, 7], [7, 6, 5]]]) + @test validate_triangulation(tri) + validate_statistics(tri) - tri = triangulate(pts, boundary_nodes=[[[1, 2, 3, 4, 1]], [[5, 7, 6, 5]]], delete_holes=false, delete_ghosts=false) - points_to_process = DT.find_all_points_to_delete(tri) - @test points_to_process == Set((DT.𝒢,)) - triangles_to_delete = DT.find_all_triangles_to_delete(tri, points_to_process) - @test DT.compare_triangle_collections(triangles_to_delete, Set(( - (1, 4, DT.𝒢), - (4, 3, DT.𝒢), - (3, 2, DT.𝒢), - (2, 1, DT.𝒢), - (5, 6, 7) - ))) - tri = triangulate(pts; boundary_nodes=[[[1, 2, 3, 4, 1]], [[5, 7, 6, 5]]]) - @test validate_triangulation(tri) - validate_statistics(tri) + tri = triangulate(pts, boundary_nodes=[[[1, 2, 3, 4, 1]], [[5, 7, 6, 5]]], delete_holes=false, delete_ghosts=false) + points_to_process = DT.find_all_points_to_delete(tri) + @test points_to_process == Set((DT.𝒢,)) + triangles_to_delete = DT.find_all_triangles_to_delete(tri, points_to_process) + @test DT.compare_triangle_collections(triangles_to_delete, Set(( + (1, 4, DT.𝒢), + (4, 3, DT.𝒢), + (3, 2, DT.𝒢), + (2, 1, DT.𝒢), + (5, 6, 7) + ))) + tri = triangulate(pts; boundary_nodes=[[[1, 2, 3, 4, 1]], [[5, 7, 6, 5]]]) + @test validate_triangulation(tri) + validate_statistics(tri) + end end @testset "A previously broken example" begin a = 4 / 5 t = LinRange(0, 2π, 6) x = @. a * (2cos(t) + cos(2t)) - y = @. a * (2sin(t) - sin(2t)) # this used to be tri = generate_mesh(x, y, 5.0); points = get_points(tri); bn_nodes = get_boundary_nodes(tri) + y = @. a * (2sin(t) - sin(2t)) points = [2.4 -0.152786404500042 -1.047213595499958 -1.047213595499958 -0.1527864045000426; 0.0 1.051462224238267 1.70130161670408 -1.70130161670408 -1.051462224238268] bn_nodes = [1, 2, 3, 4, 5, 1] tri = triangulate(points; boundary_nodes=bn_nodes, delete_ghosts=false, delete_holes=false) @@ -226,7 +228,7 @@ end a = 4 / 5 t = LinRange(0, 2π, 100) x = @. a * (2cos(t) + cos(2t)) - y = @. a * (2sin(t) - sin(2t)) # this used to be tri = generate_mesh(x, y, 15.0); points = get_points(tri); bn_nodes = get_boundary_nodes(tri) + y = @. a * (2sin(t) - sin(2t)) points = [2.4 2.390342532619651 2.361486661646358 2.313780262033189 2.247797423889998 2.16432999441036 2.064375923309772 1.949124594926212 1.819939377400058 1.678337662931399 1.525968712312314 1.364589651123909 1.196039993627106 1.022215093005981 0.8450389328830836 0.6664366846599798 0.4883074580921845 0.3124976685434042 0.1407754336467903 -0.02519360519360872 -0.1838686646186263 -0.333854062108855 -0.473917012361638 -0.6030027096667079 -0.7202464727404799 -0.824982776042241 -0.9167510419119772 -0.9952981201288541 -1.060577434849961 -1.112744832492037 -1.152151217102105 -1.179332111274063 -1.194994329879372 -1.2 -1.19534820274028 -1.182154550372295 -1.161629044931084 -1.13505259139796 -1.103752559560399 -1.069077803180918 -1.032373553014235 -0.994956601357817 -0.958091190190919 -0.9229660026456057 -0.8906726387622712 -0.8621859315182511 -0.8383464283873544 -0.8198453276894746 -0.8072121183067701 -0.8008051266355886 -0.8008051266355887 -0.8072121183067701 -0.8198453276894746 -0.8383464283873544 -0.8621859315182508 -0.890672638762271 -0.9229660026456055 -0.9580911901909194 -0.9949566013578169 -1.032373553014235 -1.069077803180918 -1.103752559560399 -1.135052591397961 -1.161629044931084 -1.182154550372295 -1.19534820274028 -1.2 -1.194994329879372 -1.179332111274063 -1.152151217102105 -1.112744832492038 -1.060577434849961 -0.9952981201288543 -0.9167510419119774 -0.8249827760422405 -0.7202464727404799 -0.6030027096667084 -0.473917012361639 -0.3338540621088567 -0.1838686646186263 -0.02519360519360934 0.1407754336467914 0.3124976685434047 0.4883074580921845 0.6664366846599791 0.8450389328830821 1.022215093005981 1.196039993627106 1.364589651123908 1.525968712312312 1.678337662931397 1.819939377400058 1.949124594926211 2.064375923309773 2.16432999441036 2.247797423889998 2.313780262033188 2.361486661646357 2.390342532619651 0.3170611563256605 0.3222593190661471 -0.6183519451946706 -0.6497375598908914 -0.03601883360826612 0.6813132651727021 -0.6442990922286358 0.677301056434819 -0.03212472989883586 1.049402030184789 -0.6932272745230642 -0.3561747556617257 1.050855112561903 -0.3586613134669976 -0.6921937990949052 -0.7738275096391242 1.420372082515759 -0.6465445728766354 -0.7401785386842585 -0.8130282799705698 1.553206818654827 -0.4939148325245979 -0.5486837455157717 -0.3610262116797264 -0.73633174592415 -0.5446205627063301 -0.5160523518705824 1.230246578448748 1.093304308222101 0.8980444010246713 1.226019949979472 -0.4997667391156656 -0.726253210863807 -0.6697924507551998 -0.4891744968539103 -0.4634488008923348 -0.301600826224346 -0.2783275986241442 -0.1177692287900377 -0.09393337228059367 -0.2546427830837441 -0.0703514373220399 0.08836743861011277 -0.2313237479039875 0.1137973151090641 0.2724811315356115 0.2996278178466897 0.4682688930880384 0.1380153936605681 0.406551246715369 0.6046552182453121 -0.426784127732563 -0.3923263038801806 -0.208465408072328 -0.3607313375439971 -0.1857892526190982 0.4960419968524275 -0.1415763437967015 0.03591030790050782 0.8685411342185969 -0.1926902902309472 -0.8592155035553608 1.716817787430798 -0.8576022838754369 -0.1921452564055306 -0.6683334786985783 0.8654953794762578 0.6952267131341282 -0.6287806122165458 0.4999252549098672 0.3218569037134379 -0.5986272435711641 0.1449039324769594 2.203335189185611 -1.101667594592808 -1.101667594592807 -0.0286998603756622 -0.3499716813694415 -0.4450559779573338 -0.04783795748406422 0.09480479796065325 2.114352958860066 -1.054192316972474 -1.054192316972475 -0.3317870997991825 1.39665732177026 -0.6206155967973439 -0.7760417249729167 0.4988340023743217 0.1461550962121355 0.3242397406248529 0.9036042833648704 -0.5300851339677813 -0.3652221764727946 -0.9373851776825108 1.875039225311675 -0.9376540476291644 -0.5066666554551625 0.2368123037070368 -0.02440461564986861 0.150646683235148 1.991415825682697 -0.9957149117466657 -0.9957009139360333 0.5336878756578392 0.7210970221965526 -0.1928809320498361 0.7689122343492278 -0.6467071096582683 1.285560087566004 -0.6388529779077368 -0.1969978949942105 0.1875732054682595 -0.4975220279896986 -0.7186754363850769 1.577166091039188 -0.8584906546541111; 0.0 0.000204308591503799 0.001629535973135443 0.005472026448394285 0.01287939060935179 0.02492716987386698 0.04259671987083449 0.06675468790926908 0.09813443380113843 0.1373197116453414 0.1847308933321656 0.2406139730886964 0.3050325470249727 0.3778629130904994 0.4587923858949879 0.5473208683088024 0.6427656684861544 0.7442694978083142 0.8508115330838226 0.9612213760111374 1.074195695220302 1.188317291935738 1.302076290158998 1.413893116908774 1.522142908049466 1.625180951076655 1.721368758301788 1.809100352482636 1.886828342268937 1.953089366954632 2.006528498919954 2.04592220767042 2.070199511290803 2.078460969082653 2.069995202699299 2.044292671697285 2.001056472471559 1.94020997634528 1.86190117239507 1.766503632611802 1.654614070392519 1.527046517275516 1.384823196404126 1.229162223576608 1.061462317070302 0.8832847449107668 0.6963327821298029 0.5024289901161495 0.3034906647750211 0.1015038293221599 -0.1015038293221605 -0.3034906647750203 -0.50242899011615 -0.6963327821298022 -0.8832847449107647 -1.061462317070301 -1.229162223576608 -1.384823196404126 -1.527046517275515 -1.654614070392518 -1.766503632611802 -1.86190117239507 -1.940209976345281 -2.001056472471559 -2.044292671697285 -2.069995202699299 -2.078460969082653 -2.070199511290803 -2.04592220767042 -2.006528498919954 -1.953089366954633 -1.886828342268938 -1.809100352482636 -1.721368758301788 -1.625180951076654 -1.522142908049467 -1.413893116908775 -1.302076290158999 -1.18831729193574 -1.074195695220302 -0.9612213760111378 -0.8508115330838222 -0.7442694978083139 -0.6427656684861544 -0.5473208683088027 -0.4587923858949887 -0.3778629130904991 -0.3050325470249729 -0.2406139730886968 -0.184730893332166 -0.137319711645342 -0.09813443380113851 -0.06675468790926926 -0.04259671987083431 -0.02492716987386698 -0.01287939060935179 -0.005472026448394374 -0.001629535973135488 -0.000204308591503799 -0.5473096449852685 0.5447620596126443 0.01428450391194042 -0.4309321583994598 0.7536074241825893 0.3536290172290472 0.4171992061645116 -0.354274248521617 -0.7697113089269275 0.1945973624956351 0.8115101356751746 -1.00610749817081 -0.1925650703607583 1.006349758355742 -0.8137846879949843 1.19333488738332 0.07348683780306886 -1.266821725186389 1.366146471157314 -1.324086653415268 -0.04174935346989145 1.135404524959176 0.9456563651429131 0.8229544359251958 -0.995445054751826 -0.9480022448317108 -0.741992356781049 -0.13995947020735 0.002345879688797254 -0.08111410706343182 0.130762025431598 -1.127145434944549 0.9963834095129508 0.6153059195085407 0.5329795234649529 0.3351323238863241 0.4528041803207281 0.2552428653861011 0.3751002341424148 0.1760741054849401 0.05622243244883154 -0.02310920543998908 0.0916097709503603 -0.1431517599572751 -0.1033611549989752 0.01112539166303821 -0.1796378655622734 -0.06423821903495724 -0.2933071885887684 0.1346998725151015 0.06389662809144629 -0.06697977586584479 -0.263147702667072 -0.3418086892373844 -0.4602230597841197 -0.5395791962871652 -0.260267079463506 0.5739282454666178 0.5007002518322551 0.2757413735876992 -0.8939205770864943 -1.486342121175829 0.0009313928164653013 1.487273513992295 0.8909822595409715 -0.6270118467477362 -0.2692763527002932 -0.1472140205276644 0.2179450297267495 0.4510351900656056 0.3329413377177828 -0.2182665505265999 -0.6631900778915942 0.0 1.908144246886934 -1.908144246886936 -0.4106426344765784 -0.6482825657476072 0.1353078965823503 -0.2192301054498285 0.3004465131589842 4.336808689942018e-17 1.825914653945082 -1.825914653945081 0.6415639095768699 -0.08968357289590935 1.254408379505216 -1.164673062563948 -0.4517076066327135 0.6541117112279523 -0.36426907951982 0.09997393244042915 0.731016147780223 -0.8325312305826258 1.623909218280243 0.0001552321360775167 -1.623753986144165 -0.5474230356714256 0.1741645772010534 -0.5950229828132197 -0.4758032489975526 8.081639737527319e-6 -1.724612653719712 1.724620735359448 0.2655807939077015 0.1867114129881211 -0.717357771233564 0.01594342553679897 -1.111060401383298 -0.004534585080420719 1.115594986463718 0.7308601982314702 0.4541306336569979 -0.3789760643386818 -1.406227077779554 0.08072235390443261 1.325504723875121] bn_nodes = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 1] tri = triangulate(points; boundary_nodes=bn_nodes) diff --git a/test/operations/delete_point.jl b/test/operations/delete_point.jl index 052a2989b..c9ac43f87 100644 --- a/test/operations/delete_point.jl +++ b/test/operations/delete_point.jl @@ -1,6 +1,7 @@ using ..DelaunayTriangulation const DT = DelaunayTriangulation using StableRNGs +using Preferences using StaticArrays @testset verbose = true "Deleting interior nodes" begin @@ -46,51 +47,54 @@ using StaticArrays end end end - @testset "Lattice with a single boundary index" begin - for j in 1:20 - rng1 = StableRNG(j) - a, b = sort(10randn(rng1, 2)) - c, d = sort(15randn(rng1, 2)) - nx = rand(rng1, 5:25) - ny = rand(rng1, 5:25) - tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=true) - points = get_points(tri) - n = nx * ny - for k in 1:(n÷10) - rng2 = StableRNG(j + k * n) - i = rand(rng2, each_solid_vertex(tri) |> collect) - while DT.is_boundary_node(tri, i)[1] + + if load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true) + @testset "Lattice with a single boundary index" begin + for j in 1:20 + rng1 = StableRNG(j) + a, b = sort(10randn(rng1, 2)) + c, d = sort(15randn(rng1, 2)) + nx = rand(rng1, 5:25) + ny = rand(rng1, 5:25) + tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=true) + points = get_points(tri) + n = nx * ny + for k in 1:(n÷10) + rng2 = StableRNG(j + k * n) i = rand(rng2, each_solid_vertex(tri) |> collect) + while DT.is_boundary_node(tri, i)[1] + i = rand(rng2, each_solid_vertex(tri) |> collect) + end + delete_point!(tri, i; rng=rng2) + @test validate_triangulation(tri) end - delete_point!(tri, i; rng=rng2) - @test validate_triangulation(tri) end + tri = triangulate_rectangle(0, 1, 0, 1, 25, 25; delete_ghosts=false, single_boundary=true) + add_segment!(tri, 7, 7 + 25) + @test_throws DT.InvalidVertexDeletionError delete_point!(tri, 2) + @test_throws DT.InvalidVertexDeletionError delete_point!(tri, 7) + @test_throws DT.InvalidVertexDeletionError delete_point!(tri, 7 + 25) + @test_throws DT.InvalidVertexDeletionError delete_point!(tri, -1) end - tri = triangulate_rectangle(0, 1, 0, 1, 25, 25; delete_ghosts=false, single_boundary=true) - add_segment!(tri, 7, 7 + 25) - @test_throws DT.InvalidVertexDeletionError delete_point!(tri, 2) - @test_throws DT.InvalidVertexDeletionError delete_point!(tri, 7) - @test_throws DT.InvalidVertexDeletionError delete_point!(tri, 7 + 25) - @test_throws DT.InvalidVertexDeletionError delete_point!(tri, -1) - end - @testset "Lattice with multiple boundary indices" begin - for j in 1:20 - rng1 = StableRNG(j) - a, b = sort(10randn(rng1, 2)) - c, d = sort(15randn(rng1, 2)) - nx = rand(rng1, 5:25) - ny = rand(rng1, 5:25) - tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=false) - points = get_points(tri) - n = nx * ny - for k in 1:(n÷10) - rng2 = StableRNG(j + k * n) - i = rand(rng2, each_solid_vertex(tri) |> collect) - while DT.is_boundary_node(tri, i)[1] + @testset "Lattice with multiple boundary indices" begin + for j in 1:20 + rng1 = StableRNG(j) + a, b = sort(10randn(rng1, 2)) + c, d = sort(15randn(rng1, 2)) + nx = rand(rng1, 5:25) + ny = rand(rng1, 5:25) + tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=false) + points = get_points(tri) + n = nx * ny + for k in 1:(n÷10) + rng2 = StableRNG(j + k * n) i = rand(rng2, each_solid_vertex(tri) |> collect) + while DT.is_boundary_node(tri, i)[1] + i = rand(rng2, each_solid_vertex(tri) |> collect) + end + delete_point!(tri, i; rng=rng2) + @test validate_triangulation(tri) end - delete_point!(tri, i; rng=rng2) - @test validate_triangulation(tri) end end end diff --git a/test/operations/lock_convex_hull.jl b/test/operations/lock_convex_hull.jl index 99aae9631..4aa922a43 100644 --- a/test/operations/lock_convex_hull.jl +++ b/test/operations/lock_convex_hull.jl @@ -4,7 +4,6 @@ using CairoMakie using DataStructures using StaticArrays - @testset "lock_convex_hull" begin for _ in 1:150 pts = rand(2, 500) diff --git a/test/predicates/general.jl b/test/predicates/general.jl index 65da51d8e..e72a94879 100644 --- a/test/predicates/general.jl +++ b/test/predicates/general.jl @@ -4,7 +4,41 @@ using StaticArrays const DT = DelaunayTriangulation using ..DelaunayTriangulation: Certificate +@testset "ext" begin + for _ in 1:50 + p = @SVector rand(2) + q = @SVector rand(2) + @test ExactPredicates.ext(p, q) == DT.ext(p..., q...) + end +end + +@testset "inp" begin + for _ in 1:50 + p = @SVector rand(2) + q = @SVector rand(2) + @test ExactPredicates.inp(p, q) == DT.inp(p..., q...) + end +end + +@testset "det" begin + for _ in 1:50 + a, b, c, d = rand(4) + @test ExactPredicates.det(a, b, c, d) == DT.det(a, b, c, d) + a, b, c, d, e, f, g, h, i = rand(9) + @test ExactPredicates.det(a, b, c, d, e, f, g, h, i) == DT.det(a, b, c, d, e, f, g, h, i) + end +end + +@testset "sgn" begin + @test DT.sgn(2.0) === 1 + @test DT.sgn(2) === 1 + @test DT.sgn(-1.0) === -1 + @test DT.sgn(-1) === -1 + @test DT.sgn(0.0) === 0 + @test DT.sgn(-0.0) === 0 + @test DT.sgn(0) === 0 +end @testset "Opposite signs" begin @test DT.opposite_signs(1, -1) @@ -22,27 +56,35 @@ end @testset "Random points" begin for _ in 1:1500 p, q, r = eachcol(rand(2, 3)) - @test DT.orient_predicate(p, q, r) == orient(p, q, r) + @test DT.orient_predicate(p, q, r) === orient(p, q, r) @inferred DT.orient_predicate(p, q, r) + @test DT._orient_predicate(p, q, r) === orient(p, q, r) + @inferred DT._orient_predicate(p, q, r) p, q, r, s = eachcol(rand(3, 4)) - @test DT.orient_predicate(p, q, r, s) == orient(p, q, r, s) + @test DT.orient_predicate(p, q, r, s) === orient(p, q, r, s) @inferred DT.orient_predicate(p, q, r, s) + @test DT._orient_predicate(p, q, r, s) === orient(p, q, r, s) + @inferred DT._orient_predicate(p, q, r, s) a, b, c, p = eachcol(rand(2, 4)) - @test DT.incircle_predicate(a, b, c, p) == incircle(a, b, c, p) + @test DT.incircle_predicate(a, b, c, p) === incircle(a, b, c, p) @inferred DT.incircle_predicate(a, b, c, p) + @test DT._incircle_predicate(a, b, c, p) === incircle(a, b, c, p) + @inferred DT._incircle_predicate(a, b, c, p) a, b, p, q = eachcol(rand(2, 4)) - @test DT.parallelorder_predicate(a, b, p, q) == parallelorder(a, b, p, q) + @test DT.parallelorder_predicate(a, b, p, q) === parallelorder(a, b, p, q) @inferred DT.parallelorder_predicate(a, b, p, q) + @test DT._parallelorder_predicate(a, b, p, q) === parallelorder(a, b, p, q) + @inferred DT._parallelorder_predicate(a, b, p, q) p, a, b = eachcol(rand(2, 3)) - @test DT.sameside_predicate(a, b, p) == sameside(p, a, b) + @test DT.sameside_predicate(a, b, p) === sameside(p, a, b) @inferred DT.sameside_predicate(a, b, p) p, q, a, b = eachcol(rand(2, 4)) - @test DT.meet_predicate(p, q, a, b) == meet(p, q, a, b) + @test DT.meet_predicate(p, q, a, b) === meet(p, q, a, b) @inferred DT.meet_predicate(p, q, a, b) end end @@ -59,15 +101,18 @@ end a = (rand(x), rand(y)) b = (rand(x), rand(y)) c = (rand(x), rand(y)) - @test DT.orient_predicate(p, q, r) == orient(p, q, r) - @test DT.incircle_predicate(a, b, c, p) == incircle(a, b, c, p) + @test DT.orient_predicate(p, q, r) === orient(p, q, r) + @test DT.incircle_predicate(a, b, c, p) === incircle(a, b, c, p) @test DT.sameside_predicate(a, b, p) == sameside(p, a, b) - @test DT.meet_predicate(p, q, a, b) == meet(p, q, a, b) + @test DT.meet_predicate(p, q, a, b) === meet(p, q, a, b) + @test DT._orient_predicate(p, q, r) === orient(p, q, r) + @test DT._incircle_predicate(a, b, c, p) === incircle(a, b, c, p) p = (rand(x), rand(y), rand(z)) q = (rand(x), rand(y), rand(z)) r = (rand(x), rand(y), rand(z)) s = (rand(x), rand(y), rand(z)) @test DT.orient_predicate(p, q, r, s) == orient(p, q, r, s) + @test DT._orient_predicate(p, q, r, s) == orient(p, q, r, s) end end @@ -75,16 +120,23 @@ end p₁ = [5.7044025422189, 1.801603986463] p₂ = [8.3797127128527, 5.8924221838871] p₃ = [2.8875415689061, 6.2038339497809] - @test DT.orient_predicate(p₁, p₂, p₃) == 1 - @test DT.orient_predicate(p₁, p₂, [10.0, 4.0]) == -1 + @test DT.orient_predicate(p₁, p₂, p₃) === 1 + @test DT.orient_predicate(p₁, p₂, [10.0, 4.0]) === -1 + @test DT._orient_predicate(p₁, p₂, p₃) === 1 + @test DT._orient_predicate(p₁, p₂, [10.0, 4.0]) === -1 p₁ = [5.0, 1.0] p₂ = [5.0, 6.0] - @test DT.orient_predicate(p₁, p₂, [5.0, 5.0]) == 0 - @test DT.orient_predicate(p₁, p₂, [5.0, 2.0]) == 0 - @test DT.orient_predicate(p₂, p₁, [5.0, 2.0]) == 0 + @test DT.orient_predicate(p₁, p₂, [5.0, 5.0]) === 0 + @test DT.orient_predicate(p₁, p₂, [5.0, 2.0]) === 0 + @test DT.orient_predicate(p₂, p₁, [5.0, 2.0]) === 0 + @test DT._orient_predicate(p₁, p₂, [5.0, 5.0]) === 0 + @test DT._orient_predicate(p₁, p₂, [5.0, 2.0]) === 0 + @test DT._orient_predicate(p₂, p₁, [5.0, 2.0]) === 0 p, q, r, s = (-4.16, -2.84, 0.0), (-4.65, -5.83, 0.0), (-1.12, -5.61, 0.0), (-7.83, 0.27, 4.6) - @test DT.orient_predicate(p, q, r, s) == -1 - @test DT.orient_predicate(p, r, q, s) == 1 + @test DT.orient_predicate(p, q, r, s) === -1 + @test DT.orient_predicate(p, r, q, s) === 1 + @test DT._orient_predicate(p, q, r, s) === -1 + @test DT._orient_predicate(p, r, q, s) === 1 end end diff --git a/test/refinement/curve_bounded.jl b/test/refinement/curve_bounded.jl index 927a1315f..9de320e17 100644 --- a/test/refinement/curve_bounded.jl +++ b/test/refinement/curve_bounded.jl @@ -10,7 +10,7 @@ using Test using StructEquality using DelimitedFiles using StableRNGs - +using Preferences const SACM = DT.SmallAngleComplexMember const SAC = DT.SmallAngleComplex @@ -1450,6 +1450,9 @@ end point_sets = deepcopy.([points_I, points_II, points_III, points_IV, points_V, points_VI, points_VII, points_VIII, points_IX, points_X, points_XI, points_XII]) curve_sets = deepcopy.([curve_I, curve_II, curve_III, curve_IV, curve_V, curve_VI, curve_VII, curve_VIII, curve_IX, curve_X, curve_XI, curve_XII]) for i in eachindex(point_sets, curve_sets) + if !load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true) && i == 4 + continue + end points, curve = deepcopy(point_sets[i]), deepcopy(curve_sets[i]) tri = triangulate(points; boundary_nodes=curve, enrich=i ≤ 3) @test validate_triangulation(tri) @@ -1517,6 +1520,9 @@ end (1e-1, 1e-2) ] for curve_idx in 4:lastindex(curve_sets) + if curve_idx ∈ (4, 6, 12) && !load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true) + continue + end for point_idx in 1:3 point_idx == 3 && curve_idx ≥ 5 && continue # no extra segments for curves ≥ 5 for (idx1, use_lens) in enumerate((false, true)) diff --git a/test/refinement/refine.jl b/test/refinement/refine.jl index 989a6d22f..41f12e931 100644 --- a/test/refinement/refine.jl +++ b/test/refinement/refine.jl @@ -8,6 +8,7 @@ using ReferenceTests using Test using StructEquality using DelimitedFiles +using Preferences using StableRNGs @struct_equal DT.InsertionEventHistory @struct_equal DT.RefinementConstraints @@ -434,7 +435,9 @@ end DT.undo_insertion!(tri, history) @test tri.boundary_edge_map == orig_tri.boundary_edge_map validate_statistics(tri) - @test validate_triangulation(tri) + if load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true) + @test validate_triangulation(tri) + end @test tri == orig_tri add_point!(tri, p, store_event_history=Val(true), event_history=history) orig_tri = deepcopy(tri) @@ -1186,7 +1189,7 @@ end tri_5 = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=false) tri_6 = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=true) @static if VERSION ≥ v"1.10" - @inferred triangulate(rand(2, 250)) + @inferred triangulate(rand(2, 250)) end for (iii, tri) in enumerate((tri_1, tri_2, tri_3, tri_4, tri_5, tri_6)) @info "Testing if encroached edges are detected. Run: $iii." @@ -1707,92 +1710,94 @@ end end end -@testset "enqueueing and splitting all encroached segments" begin - for iii in 1:100 - for use_lens in (false, true) - for pass in 1:2 - if pass == 1 - p1 = (0.0, 0.0) - p2 = (1.0, 0.0) - p3 = (0.0, 1.0) - p4 = (1.0, 1.0) - p5 = (0.5, 0.5) - pts = [p1, p2, p3, p4, p5] - C = Set{NTuple{2,Int}}() - for i in 1:15 - θ = 2π * rand() - r = 0.5sqrt(rand()) - x = 0.5 + r * cos(θ) - y = 0.5 + r * sin(θ) - push!(pts, (x, y)) - push!(C, (5, 5 + i)) - end - tri = triangulate(pts; delete_ghosts=false, boundary_nodes=[1, 2, 4, 3, 1], segments=C) - args = DT.RefinementArguments(tri, use_lens=use_lens, min_angle=25.0, seditious_angle=15.0, max_area=0.1) - DT.enqueue_all_encroached_segments!(args, tri) - @inferred DT.enqueue_all_encroached_segments!(args, tri) - manual_enqueue = PriorityQueue{NTuple{2,Int},Float64}(Base.Order.Reverse) - for e in each_edge(tri) - if DT.contains_segment(tri, e...) - flag = DT.is_encroached(tri, args, e) - if flag - manual_enqueue[e] = DT.edge_length_sqr(tri, e...) +if load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true) + @testset "enqueueing and splitting all encroached segments" begin + for iii in 1:100 + for use_lens in (false, true) + for pass in 1:2 + if pass == 1 + p1 = (0.0, 0.0) + p2 = (1.0, 0.0) + p3 = (0.0, 1.0) + p4 = (1.0, 1.0) + p5 = (0.5, 0.5) + pts = [p1, p2, p3, p4, p5] + C = Set{NTuple{2,Int}}() + for i in 1:15 + θ = 2π * rand() + r = 0.5sqrt(rand()) + x = 0.5 + r * cos(θ) + y = 0.5 + r * sin(θ) + push!(pts, (x, y)) + push!(C, (5, 5 + i)) + end + tri = triangulate(pts; delete_ghosts=false, boundary_nodes=[1, 2, 4, 3, 1], segments=C) + args = DT.RefinementArguments(tri, use_lens=use_lens, min_angle=25.0, seditious_angle=15.0, max_area=0.1) + DT.enqueue_all_encroached_segments!(args, tri) + @inferred DT.enqueue_all_encroached_segments!(args, tri) + manual_enqueue = PriorityQueue{NTuple{2,Int},Float64}(Base.Order.Reverse) + for e in each_edge(tri) + if DT.contains_segment(tri, e...) + flag = DT.is_encroached(tri, args, e) + if flag + manual_enqueue[e] = DT.edge_length_sqr(tri, e...) + end end end - end - compare_encroach_queues(args, manual_enqueue) - empty!(args.events) - DT.split_all_encroached_segments!(tri, args) - @test isempty(args.queue.segments) - @test !isempty(args.events.added_triangles) - DT.enqueue_all_encroached_segments!(args, tri) - @test isempty(args.queue.segments) - if use_lens - for e in each_segment(tri) - @test !DT.is_encroached(tri, args, e) + compare_encroach_queues(args, manual_enqueue) + empty!(args.events) + DT.split_all_encroached_segments!(tri, args) + @test isempty(args.queue.segments) + @test !isempty(args.events.added_triangles) + DT.enqueue_all_encroached_segments!(args, tri) + @test isempty(args.queue.segments) + if use_lens + for e in each_segment(tri) + @test !DT.is_encroached(tri, args, e) + end + else + @test is_conformal(tri) end + @test validate_triangulation(tri) else - @test is_conformal(tri) - end - @test validate_triangulation(tri) - else - points = [(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)] - for r in LinRange(0.1, 0.9, 15) - push!(points, (r, 0.01)) - push!(points, (r, 0.99)) - push!(points, (0.05, r)) - push!(points, (0.95, r)) - end - for _ in 1:10 - push!(points, (rand(), rand())) - end - tri = triangulate(points; boundary_nodes=[1, 2, 3, 4, 1]) - args = DT.RefinementArguments(tri, use_lens=use_lens, min_angle=29.0, seditious_angle=19.0, max_area=0.05) - DT.enqueue_all_encroached_segments!(args, tri) - manual_enqueue = PriorityQueue{NTuple{2,Int},Float64}(Base.Order.Reverse) - for e in each_edge(tri) - if DT.contains_segment(tri, e...) - flag = DT.is_encroached(tri, args, e) - if flag - manual_enqueue[e] = DT.edge_length_sqr(tri, e...) + points = [(0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0)] + for r in LinRange(0.1, 0.9, 15) + push!(points, (r, 0.01)) + push!(points, (r, 0.99)) + push!(points, (0.05, r)) + push!(points, (0.95, r)) + end + for _ in 1:10 + push!(points, (rand(), rand())) + end + tri = triangulate(points; boundary_nodes=[1, 2, 3, 4, 1]) + args = DT.RefinementArguments(tri, use_lens=use_lens, min_angle=29.0, seditious_angle=19.0, max_area=0.05) + DT.enqueue_all_encroached_segments!(args, tri) + manual_enqueue = PriorityQueue{NTuple{2,Int},Float64}(Base.Order.Reverse) + for e in each_edge(tri) + if DT.contains_segment(tri, e...) + flag = DT.is_encroached(tri, args, e) + if flag + manual_enqueue[e] = DT.edge_length_sqr(tri, e...) + end end end - end - compare_encroach_queues(args, manual_enqueue) - empty!(args.events) - DT.split_all_encroached_segments!(tri, args) - @test isempty(args.queue.segments) - @test !isempty(args.events.added_triangles) - DT.enqueue_all_encroached_segments!(args, tri) - @test isempty(args.queue.segments) - if use_lens - for e in each_segment(tri) - @test !DT.is_encroached(tri, args, e) + compare_encroach_queues(args, manual_enqueue) + empty!(args.events) + DT.split_all_encroached_segments!(tri, args) + @test isempty(args.queue.segments) + @test !isempty(args.events.added_triangles) + DT.enqueue_all_encroached_segments!(args, tri) + @test isempty(args.queue.segments) + if use_lens + for e in each_segment(tri) + @test !DT.is_encroached(tri, args, e) + end + else + @test is_conformal(tri) end - else - @test is_conformal(tri) + @test validate_triangulation(tri) end - @test validate_triangulation(tri) end end end @@ -2272,10 +2277,10 @@ end @test validate_refinement(tri, args) if _rng_num(idx1, idx2, idx3, idx4, idx5) == _rng_num(1, 3, 1, 2, 2) fig, ax, sc = triplot(tri) - @test_reference "a_simple_convex_example_circle.png" fig by=psnr_equality(15) + @test_reference "a_simple_convex_example_circle.png" fig by = psnr_equality(15) elseif _rng_num(idx1, idx2, idx3, idx4, idx5) == _rng_num(2, 3, 1, 2, 2) fig, ax, sc = triplot(tri) - @test_reference "a_simple_convex_example_lens.png" fig by=psnr_equality(15) + @test_reference "a_simple_convex_example_lens.png" fig by = psnr_equality(15) end end end @@ -2300,10 +2305,10 @@ end @test validate_refinement(tri, DT.RefinementArguments(tri; min_angle, min_area, seditious_angle, max_area, use_circumcenter=true, use_lens)) if _rng_num(idx1, idx2, idx3, idx4, idx5) == _rng_num(1, 3, 1, 3, 1) fig, ax, sc = triplot(tri) - @test_reference "triangulation_with_a_hole_circle.png" fig by=psnr_equality(15) + @test_reference "triangulation_with_a_hole_circle.png" fig by = psnr_equality(15) elseif _rng_num(idx1, idx2, idx3, idx4, idx5) == _rng_num(2, 3, 1, 3, 1) fig, ax, sc = triplot(tri) - @test_reference "triangulation_with_a_hole_lens.png" fig by=psnr_equality(15) + @test_reference "triangulation_with_a_hole_lens.png" fig by = psnr_equality(15) end end end @@ -2390,10 +2395,10 @@ end @test validate_refinement(tri; min_angle, min_area, max_area, use_circumcenter=true, seditious_angle, use_lens, check_conformal=false) if _rng_num(idx1, idx2, idx3, idx4, idx5) == _rng_num(1, 3, 1, 3, 1) fig, ax, sc = triplot(tri) - @test_reference "a_constrained_triangulation_with_multiple_holes_circle.png" fig by=psnr_equality(15) + @test_reference "a_constrained_triangulation_with_multiple_holes_circle.png" fig by = psnr_equality(15) elseif _rng_num(idx1, idx2, idx3, idx4, idx5) == _rng_num(2, 3, 1, 3, 1) fig, ax, sc = triplot(tri) - @test_reference "a_constrained_triangulation_with_multiple_holes_lens.png" fig by=psnr_equality(15) + @test_reference "a_constrained_triangulation_with_multiple_holes_lens.png" fig by = psnr_equality(15) end end end @@ -2428,10 +2433,10 @@ end @test validate_refinement(tri; min_angle, min_area, max_area, seditious_angle, use_circumcenter=true, use_lens) if _rng_num(idx1, idx2, idx3, idx4, idx5) == _rng_num(1, 4, 1, 3, 2) fig, ax, sc = triplot(tri) - @test_reference "another_square_example_circle.png" fig by=psnr_equality(15) + @test_reference "another_square_example_circle.png" fig by = psnr_equality(15) elseif _rng_num(idx1, idx2, idx3, idx4, idx5) == _rng_num(2, 4, 1, 3, 2) fig, ax, sc = triplot(tri) - @test_reference "another_square_example_lens.png" fig by=psnr_equality(15) + @test_reference "another_square_example_lens.png" fig by = psnr_equality(15) end end end @@ -2475,10 +2480,10 @@ end validate_statistics(tri) if _rng_num(idx1, idx2, idx3, idx4, idx5) == _rng_num(1, 3, 1, 2, 2) fig, ax, sc = triplot(tri) - @test_reference "avoid_bouncing_concentric_circular_shells_circle.png" fig by=psnr_equality(15) + @test_reference "avoid_bouncing_concentric_circular_shells_circle.png" fig by = psnr_equality(15) elseif _rng_num(idx1, idx2, idx3, idx4, idx5) == _rng_num(2, 3, 1, 2, 2) fig, ax, sc = triplot(tri) - @test_reference "avoid_bouncing_concentric_circular_shells_lens.png" fig by=psnr_equality(15) + @test_reference "avoid_bouncing_concentric_circular_shells_lens.png" fig by = psnr_equality(15) end end end @@ -2492,14 +2497,14 @@ end points = [(0.0, 0.0), (1.0, 0.0), (0.8, 0.2)] tri = triangulate(points) refine!(tri, use_circumcenter=true) - @test_reference "seditious_edge_testing_1.png" triplot(tri) by=psnr_equality(15) + @test_reference "seditious_edge_testing_1.png" triplot(tri) by = psnr_equality(15) @test validate_refinement(tri, use_circumcenter=true) validate_statistics(tri) points = [(0.0, 0.0), (1.0, 0.0), (0.8, 0.2)] tri = triangulate(points) refine!(tri, use_circumcenter=true, use_lens=false) - @test_reference "seditious_edge_testing_2.png" triplot(tri) by=psnr_equality(15) + @test_reference "seditious_edge_testing_2.png" triplot(tri) by = psnr_equality(15) @test validate_refinement(tri, use_circumcenter=true, use_lens=false) validate_statistics(tri) @@ -2508,7 +2513,7 @@ end args = DT.RefinementArguments(tri, use_circumcenter=true, use_lens=false, max_area=1e-3get_area(tri)) refine!(tri, args) @test validate_refinement(tri, args) - @test_reference "seditious_edge_testing_3.png" triplot(tri) by=psnr_equality(15) + @test_reference "seditious_edge_testing_3.png" triplot(tri) by = psnr_equality(15) validate_statistics(tri) end end @@ -2544,7 +2549,7 @@ end @test validate_triangulation(tri) @test validate_refinement(tri; max_area=0.001, max_points=5000, use_circumcenter=true, use_lens=true) fig, ax, sc = triplot(tri) - @test_reference "triangulating_with_an_interior_hole.png" fig by=psnr_equality(15) + @test_reference "triangulating_with_an_interior_hole.png" fig by = psnr_equality(15) end @testset "Refining disjoint sets" begin @@ -2568,50 +2573,52 @@ end validate_statistics(tri) @test validate_refinement(tri; min_area, max_area, use_circumcenter=true) fig, ax, sc = triplot(tri) - @test_reference "refining_disjoint_sets.png" fig by=psnr_equality(15) + @test_reference "refining_disjoint_sets.png" fig by = psnr_equality(15) end - @testset "Small angles" begin - ps = 0 - fig = Figure(fontsize=52) - for i in 1:12 - if i > 6 - use_lens = true - i -= 6 - ax = Axis(fig[2, i], title="Lens; $i", width=600, height=600) - else - use_lens = false - ax = Axis(fig[1, i], title="Circle; $i", width=600, height=600) - end - hidedecorations!(ax) - hidespines!(ax) - rng = StableRNG(i) - p1 = (0.0, 0.0) - p2 = (1.0, 0.0) - p3 = (0.0, 1.0) - p4 = (1.0, 1.0) - p5 = (0.5, 0.5) - pts = [p1, p2, p3, p4, p5] - C = Set{NTuple{2,Int}}() - for i in 1:20 - θ = 2π * rand(rng) - r = 0.5sqrt(rand(rng)) - x = 0.5 + r * cos(θ) - y = 0.5 + r * sin(θ) - push!(pts, (x, y)) - push!(C, (5, 5 + i)) + if load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true) + @testset "Small angles" begin + ps = 0 + fig = Figure(fontsize=52) + for i in 1:12 + if i > 6 + use_lens = true + i -= 6 + ax = Axis(fig[2, i], title="Lens; $i", width=600, height=600) + else + use_lens = false + ax = Axis(fig[1, i], title="Circle; $i", width=600, height=600) + end + hidedecorations!(ax) + hidespines!(ax) + rng = StableRNG(i) + p1 = (0.0, 0.0) + p2 = (1.0, 0.0) + p3 = (0.0, 1.0) + p4 = (1.0, 1.0) + p5 = (0.5, 0.5) + pts = [p1, p2, p3, p4, p5] + C = Set{NTuple{2,Int}}() + for i in 1:20 + θ = 2π * rand(rng) + r = 0.5sqrt(rand(rng)) + x = 0.5 + r * cos(θ) + y = 0.5 + r * sin(θ) + push!(pts, (x, y)) + push!(C, (5, 5 + i)) + end + tri = triangulate(pts; rng, boundary_nodes=[1, 2, 4, 3, 1], segments=C) + refine!(tri; min_angle=27.3, min_area=0.0, use_circumcenter=true, rng, use_lens) + stats = statistics(tri) + ps += DT.get_largest_angle(stats) ≤ max(π - 2 * deg2rad(17.0), 2asin((sqrt(3) - 1) / 2)) # Corollary 8 of "When and Why Ruppert's Algorithm Works. In Twelfth International Meshing Roundtable, pp. 91–102, Santa Fe, NM, Sept 2003." + validate_statistics(tri) + @test validate_refinement(tri; min_angle=27.3, min_area=0.0, use_circumcenter=true, warn=false, use_lens, rng) + triplot!(ax, tri) end - tri = triangulate(pts; rng, boundary_nodes=[1, 2, 4, 3, 1], segments=C) - refine!(tri; min_angle=27.3, min_area=0.0, use_circumcenter=true, rng, use_lens) - stats = statistics(tri) - ps += DT.get_largest_angle(stats) ≤ max(π - 2 * deg2rad(17.0), 2asin((sqrt(3) - 1) / 2)) # Corollary 8 of "When and Why Ruppert's Algorithm Works. In Twelfth International Meshing Roundtable, pp. 91–102, Santa Fe, NM, Sept 2003." - validate_statistics(tri) - @test validate_refinement(tri; min_angle=27.3, min_area=0.0, use_circumcenter=true, warn=false, use_lens, rng) - triplot!(ax, tri) + resize_to_layout!(fig) + @test_reference "refinement_with_small_angles.png" fig + @test ps == 12 end - resize_to_layout!(fig) - @test_reference "refinement_with_small_angles.png" fig - @test ps == 12 end @testset "Another disjoint domain" begin @@ -2643,7 +2650,7 @@ end validate_statistics(tri) @test validate_refinement(tri, max_area=1e-3A, min_area=0.0, use_circumcenter=true, use_lens=true, warn=false) fig, ax, sc = triplot(tri) - @test_reference "refining_disjoint_sets_2.png" fig by=psnr_equality(15) + @test_reference "refining_disjoint_sets_2.png" fig by = psnr_equality(15) end @testset "Tight example with a single triangle boundary interior" begin @@ -2695,7 +2702,7 @@ end validate_statistics(tri) @test validate_refinement(tri; max_area=1e-3A, rng, use_circumcenter=true, warn=true) fig, ax, sc = triplot(tri) - @test_reference "tight_example_with_a_single_triangle_boundary_interior.png" fig by=psnr_equality(15) + @test_reference "tight_example_with_a_single_triangle_boundary_interior.png" fig by = psnr_equality(15) end @testset "A complicated example with tight walls and small angles" begin @@ -2767,27 +2774,29 @@ end refine!(tri; max_area=1.0, rng, use_circumcenter=true) validate_statistics(tri) fig, ax, sc = triplot(tri) - @test_reference "complicated_example_with_tight_walls_and_small_angles.png" fig by=psnr_equality(15) + @test_reference "complicated_example_with_tight_walls_and_small_angles.png" fig by = psnr_equality(15) end - if !(get(ENV, "CI", "false") == "true") - @testset "Tasmania" begin - rng = StableRNG(123) - tassy_path = joinpath(dirname(dirname(pathof(DelaunayTriangulation))), "test", "tassy.txt") - tassy = readdlm(tassy_path) - ymax = @views maximum(tassy[:, 2]) - tassy = [(x, ymax - y) for (x, y) in eachrow(tassy)] - reverse!(tassy) - unique!(tassy) - push!(tassy, tassy[begin]) - boundary_nodes, points = convert_boundary_points_to_indices(tassy) - tri = triangulate(points; rng, boundary_nodes=boundary_nodes) - A = get_area(tri) - refine!(tri; max_area=1e-2A, use_circumcenter=true, rng) - validate_statistics(tri) - @test validate_refinement(tri; max_area=1e-2A, rng, use_circumcenter=true, warn=false) - fig, ax, sc = triplot(tri) - @test_reference "tasmania.png" fig by=psnr_equality(15) + if load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true) + if !(get(ENV, "CI", "false") == "true") + @testset "Tasmania" begin + rng = StableRNG(123) + tassy_path = joinpath(dirname(dirname(pathof(DelaunayTriangulation))), "test", "tassy.txt") + tassy = readdlm(tassy_path) + ymax = @views maximum(tassy[:, 2]) + tassy = [(x, ymax - y) for (x, y) in eachrow(tassy)] + reverse!(tassy) + unique!(tassy) + push!(tassy, tassy[begin]) + boundary_nodes, points = convert_boundary_points_to_indices(tassy) + tri = triangulate(points; rng, boundary_nodes=boundary_nodes) + A = get_area(tri) + refine!(tri; max_area=1e-2A, use_circumcenter=true, rng) + validate_statistics(tri) + @test validate_refinement(tri; max_area=1e-2A, rng, use_circumcenter=true, warn=false) + fig, ax, sc = triplot(tri) + @test_reference "tasmania.png" fig by = psnr_equality(15) + end end end @@ -3051,6 +3060,6 @@ end validate_statistics(tri) @test validate_refinement(tri; min_angle=26.45, max_area=0.005A / 9, rng, use_circumcenter=true) fig, ax, sc = triplot(tri) - @test_reference "julia_logo.png" fig by=psnr_equality(15) + @test_reference "julia_logo.png" fig by = psnr_equality(15) end end \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index bc52eb882..384b1941b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,3 +1,6 @@ +@show ENV["USE_EXACTPREDICATES"] +throw("...") + # get all the compilation out of the way using BenchmarkTools using CairoMakie @@ -30,6 +33,7 @@ using StatsBase using StructEquality using Aqua using Test +using Preferences const ALL_TEST_SCRIPTS = Set{String}() const NON_TEST_SCRIPTS = Set{String}(["helper_functions.jl", "triangulation_validation.jl", "runtests.jl", "triangulation\\weighted.jl", "triangulation/weighted.jl"]) @@ -50,7 +54,7 @@ end @testset verbose = true "DelaunayTriangulation.jl" begin @testset verbose = true "Aqua" begin - Aqua.test_all(DelaunayTriangulation; ambiguities=false, project_extras=false) # don't care about julia < 1.2 + Aqua.test_all(DelaunayTriangulation; ambiguities=false, project_extras=false, stale_deps = load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true)) # don't care about julia < 1.2 Aqua.test_ambiguities(DelaunayTriangulation) # don't pick up Base and Core... end diff --git a/test/triangulation/bowyer_watson.jl b/test/triangulation/bowyer_watson.jl index 888109f08..871dee261 100644 --- a/test/triangulation/bowyer_watson.jl +++ b/test/triangulation/bowyer_watson.jl @@ -4,8 +4,7 @@ using Random using StableRNGs using CairoMakie using DataStructures - - +using Preferences @testset "Getting the correct order" begin points = rand(2, 50) @@ -98,12 +97,14 @@ end @test validate_triangulation(tri) end -@testset "Lots of collinearity" begin - _tri = triangulate_rectangle(-3.0, 2.0, 5.0, 17.3, 23, 57; single_boundary=true) - @test validate_triangulation(_tri) - for _ in 1:100 - tri = DT.triangulate(_tri.points) - @test validate_triangulation(tri) +if load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true) + @testset "Lots of collinearity" begin + _tri = triangulate_rectangle(-3.0, 2.0, 5.0, 17.3, 23, 57; single_boundary=true) + @test validate_triangulation(_tri) + for _ in 1:100 + tri = DT.triangulate(_tri.points) + @test validate_triangulation(tri) + end end end diff --git a/test/triangulation/constrained.jl b/test/triangulation/constrained.jl index 6daee652e..53ed60f12 100644 --- a/test/triangulation/constrained.jl +++ b/test/triangulation/constrained.jl @@ -3,7 +3,7 @@ const DT = DelaunayTriangulation using CairoMakie using Random using StableRNGs - +using Preferences @testset "Random constrained Delaunay triangulations" begin for i in 1:25 @@ -91,81 +91,83 @@ end end end -@testset "Lattice" begin - for m in 1:2 - @info "Testing dense lattice. Run: $m" - rng = StableRNG(m) - a = 0.0 - b = 5.0 - c = -3.0 - d = 7.0 - nx = 13 - ny = 20 - tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=true) - add_segment!(tri, 56, 162; rng) - for e in [(1, 249), (1, 250), (1, 251), (1, 26), (1, 39), (1, 52)] - add_segment!(tri, e; rng) - end - add_segment!(tri, 190, 99; rng) - for e in [(99, 113), (113, 101), (101, 115)] - add_segment!(tri, e; rng) - end - @test validate_triangulation(tri) +if load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true) + @testset "Lattice" begin + for m in 1:2 + @info "Testing dense lattice. Run: $m" + rng = StableRNG(m) + a = 0.0 + b = 5.0 + c = -3.0 + d = 7.0 + nx = 13 + ny = 20 + tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=true) + add_segment!(tri, 56, 162; rng) + for e in [(1, 249), (1, 250), (1, 251), (1, 26), (1, 39), (1, 52)] + add_segment!(tri, e; rng) + end + add_segment!(tri, 190, 99; rng) + for e in [(99, 113), (113, 101), (101, 115)] + add_segment!(tri, e; rng) + end + @test validate_triangulation(tri) - a = -0.1 - b = 0.1 - c = -0.01 - d = 0.01 - nx = 25 - ny = 25 - tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=true) - tri = triangulate(get_points(tri)) - for i in 2:24 - add_segment!(tri, i, 600 + i; rng) - end - @test validate_triangulation(tri) - tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=true) - tri = triangulate(get_points(tri); rng) - for e in [(1, 28), (28, 103), (103, 180), (180, 625), (625, 523)] - add_segment!(tri, e; rng) - end - for e in [(437, 614), (527, 602), (528, 603), (555, 605)] - add_segment!(tri, e; rng) + a = -0.1 + b = 0.1 + c = -0.01 + d = 0.01 + nx = 25 + ny = 25 + tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=true) + tri = triangulate(get_points(tri)) + for i in 2:24 + add_segment!(tri, i, 600 + i; rng) + end + @test validate_triangulation(tri) + tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=true) + tri = triangulate(get_points(tri); rng) + for e in [(1, 28), (28, 103), (103, 180), (180, 625), (625, 523)] + add_segment!(tri, e; rng) + end + for e in [(437, 614), (527, 602), (528, 603), (555, 605)] + add_segment!(tri, e; rng) + end + @test validate_triangulation(tri) end - @test validate_triangulation(tri) - end - for m in 1:20 - rng = StableRNG(m) - @info "Testing coarse lattice. Run: $m" - a = 0.0 - b = 1.0 - c = 0.0 - d = 1.0 - nx = 2 - ny = 2 - tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=true) - add_segment!(tri, 1, 4; rng) - @test validate_triangulation(tri) + for m in 1:20 + rng = StableRNG(m) + @info "Testing coarse lattice. Run: $m" + a = 0.0 + b = 1.0 + c = 0.0 + d = 1.0 + nx = 2 + ny = 2 + tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=true) + add_segment!(tri, 1, 4; rng) + @test validate_triangulation(tri) - a = 0 - b = 1 - c = 0 - d = 5 - nx = 25 - ny = 3 - tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=true) - tri = triangulate(get_points(tri); rng) - for i in 1:(nx-1) - u = i - v = 2nx - add_segment!(tri, u, v) - end - for i in 51:75 - u = i - v = 26 - add_segment!(tri, u, v; rng) + a = 0 + b = 1 + c = 0 + d = 5 + nx = 25 + ny = 3 + tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=true) + tri = triangulate(get_points(tri); rng) + for i in 1:(nx-1) + u = i + v = 2nx + add_segment!(tri, u, v) + end + for i in 51:75 + u = i + v = 26 + add_segment!(tri, u, v; rng) + end + @test validate_triangulation(tri) end - @test validate_triangulation(tri) end end @@ -185,34 +187,36 @@ end end end -@testset "Triangulation with two curves" begin - for i in 1:10 - @info "Testing triangulation of a domain with two curves: Run $i" - rng = StableRNG(i) - pts = [(rand(rng), rand(rng)) for _ in 1:50] - x1 = LinRange(0, 1, 100) - y1 = LinRange(0, 0, 100) - x2 = LinRange(1, 1, 100) - y2 = LinRange(0, 1, 100) - x3 = LinRange(1, 0, 100) - y3 = LinRange(1, 1, 100) - x4 = LinRange(0, 0, 100) - y4 = LinRange(1, 0, 100) - x = [x1..., x2..., x3..., x4...] - y = [y1..., y2..., y3..., y4...] - outer_square = [(x, y) for (x, y) in zip(x, y)] |> unique - push!(outer_square, outer_square[begin]) - outer_square_x = [first.(outer_square)] - outer_square_y = [last.(outer_square)] - circ_pts = [(0.3cos(θ), 0.3sin(θ)) .+ 0.5 for θ in LinRange(2π, 0, 50)] - circ_pts[end] = circ_pts[1] - inner_circle_x = [first.(circ_pts)] - inner_circle_y = [last.(circ_pts)] - x = [outer_square_x, inner_circle_x] - y = [outer_square_y, inner_circle_y] - nodes, pts = convert_boundary_points_to_indices(x, y; existing_points=pts) - tri = triangulate(pts; boundary_nodes=nodes, rng) - @test validate_triangulation(tri) +if load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true) + @testset "Triangulation with two curves" begin + for i in 1:10 + @info "Testing triangulation of a domain with two curves: Run $i" + rng = StableRNG(i) + pts = [(rand(rng), rand(rng)) for _ in 1:50] + x1 = LinRange(0, 1, 100) + y1 = LinRange(0, 0, 100) + x2 = LinRange(1, 1, 100) + y2 = LinRange(0, 1, 100) + x3 = LinRange(1, 0, 100) + y3 = LinRange(1, 1, 100) + x4 = LinRange(0, 0, 100) + y4 = LinRange(1, 0, 100) + x = [x1..., x2..., x3..., x4...] + y = [y1..., y2..., y3..., y4...] + outer_square = [(x, y) for (x, y) in zip(x, y)] |> unique + push!(outer_square, outer_square[begin]) + outer_square_x = [first.(outer_square)] + outer_square_y = [last.(outer_square)] + circ_pts = [(0.3cos(θ), 0.3sin(θ)) .+ 0.5 for θ in LinRange(2π, 0, 50)] + circ_pts[end] = circ_pts[1] + inner_circle_x = [first.(circ_pts)] + inner_circle_y = [last.(circ_pts)] + x = [outer_square_x, inner_circle_x] + y = [outer_square_y, inner_circle_y] + nodes, pts = convert_boundary_points_to_indices(x, y; existing_points=pts) + tri = triangulate(pts; boundary_nodes=nodes, rng) + @test validate_triangulation(tri) + end end end diff --git a/test/voronoi/voronoi.jl b/test/voronoi/voronoi.jl index a01956611..d957bae9f 100644 --- a/test/voronoi/voronoi.jl +++ b/test/voronoi/voronoi.jl @@ -8,6 +8,7 @@ import GeometryBasics: Point2f using StaticArrays using LinearAlgebra using StructEquality +using Preferences @struct_equal DT.Queue @testset "Unconstrained test" begin @@ -194,13 +195,13 @@ end end @testset "Voronoi point location" begin - A = (-1.0, 7.0) - B = (4.0, 4.0) - C = (-2.0, -1.0) - D = (-1.0, 3.0) - E = (3.0, -1.0) - F = (1.0, 4.0) - G = (-3.0, 5.0) + A = (-1.0, 7.0) .+ (1e-6rand(), 1e-6rand()) # perturb to allow the tests to work even without ExactPredicates + B = (4.0, 4.0) .+ (1e-6rand(), 1e-6rand()) + C = (-2.0, -1.0) .+ (1e-6rand(), 1e-6rand()) + D = (-1.0, 3.0) .+ (1e-6rand(), 1e-6rand()) + E = (3.0, -1.0) .+ (1e-6rand(), 1e-6rand()) + F = (1.0, 4.0) .+ (1e-6rand(), 1e-6rand()) + G = (-3.0, 5.0) .+ (1e-6rand(), 1e-6rand()) pts = [A, B, C, D, E, F, G] tri = triangulate(pts; delete_ghosts=false, randomise=false) vor = voronoi(tri) From 4b281841c7368c9fb4982d9e55aff24fb6a3b474 Mon Sep 17 00:00:00 2001 From: DanielVandH Date: Sun, 7 Jul 2024 20:22:29 +0100 Subject: [PATCH 02/14] Use the ENV --- test/runtests.jl | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 384b1941b..6308e53f5 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,13 @@ -@show ENV["USE_EXACTPREDICATES"] +# setup LocalPreferences.toml +using Preferences +USE_EXACTPREDICATES = ENV["USE_EXACTPREDICATES"] +if USE_EXACTPREDICATES == "true" + @set_preferences!("DelaunayTriangulation", "USE_EXACTPREDICATES" => true) +elseif USE_EXACTPREDICATES == "false" + @set_preferences!("DelaunayTriangulation", "USE_EXACTPREDICATES" => false) +end # if USE_EXACTPREDICATES == "default", do nothing +@show DelaunayTriangulation.USE_EXACTPREDICATES + throw("...") # get all the compilation out of the way From 740d8bccc498c92ca6efa62231bd30191916d01a Mon Sep 17 00:00:00 2001 From: DanielVandH Date: Sun, 7 Jul 2024 20:28:40 +0100 Subject: [PATCH 03/14] fix --- test/runtests.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 6308e53f5..3b280461b 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -2,9 +2,9 @@ using Preferences USE_EXACTPREDICATES = ENV["USE_EXACTPREDICATES"] if USE_EXACTPREDICATES == "true" - @set_preferences!("DelaunayTriangulation", "USE_EXACTPREDICATES" => true) + set_preferences!("DelaunayTriangulation", "USE_EXACTPREDICATES" => true) elseif USE_EXACTPREDICATES == "false" - @set_preferences!("DelaunayTriangulation", "USE_EXACTPREDICATES" => false) + set_preferences!("DelaunayTriangulation", "USE_EXACTPREDICATES" => false) end # if USE_EXACTPREDICATES == "default", do nothing @show DelaunayTriangulation.USE_EXACTPREDICATES From 4106f386c79bbd78895b26e2c7775fd478fffab6 Mon Sep 17 00:00:00 2001 From: DanielVandH Date: Sun, 7 Jul 2024 20:37:25 +0100 Subject: [PATCH 04/14] t --- test/runtests.jl | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 3b280461b..32f2e85c1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -6,9 +6,6 @@ if USE_EXACTPREDICATES == "true" elseif USE_EXACTPREDICATES == "false" set_preferences!("DelaunayTriangulation", "USE_EXACTPREDICATES" => false) end # if USE_EXACTPREDICATES == "default", do nothing -@show DelaunayTriangulation.USE_EXACTPREDICATES - -throw("...") # get all the compilation out of the way using BenchmarkTools @@ -42,7 +39,11 @@ using StatsBase using StructEquality using Aqua using Test -using Preferences + +@show DelaunayTriangulation.USE_EXACTPREDICATES + +throw("...") + const ALL_TEST_SCRIPTS = Set{String}() const NON_TEST_SCRIPTS = Set{String}(["helper_functions.jl", "triangulation_validation.jl", "runtests.jl", "triangulation\\weighted.jl", "triangulation/weighted.jl"]) From 3b0e17f2f5f6f61c38d658fffb9eb8d1a5c8abc9 Mon Sep 17 00:00:00 2001 From: DanielVandH Date: Sun, 7 Jul 2024 20:54:50 +0100 Subject: [PATCH 05/14] Go CI --- test/runtests.jl | 5 ----- 1 file changed, 5 deletions(-) diff --git a/test/runtests.jl b/test/runtests.jl index 32f2e85c1..ef316594a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -40,11 +40,6 @@ using StructEquality using Aqua using Test -@show DelaunayTriangulation.USE_EXACTPREDICATES - -throw("...") - - const ALL_TEST_SCRIPTS = Set{String}() const NON_TEST_SCRIPTS = Set{String}(["helper_functions.jl", "triangulation_validation.jl", "runtests.jl", "triangulation\\weighted.jl", "triangulation/weighted.jl"]) include("helper_functions.jl") From 37f22d9d5aa7b46fec475dba0bc82f547e09647c Mon Sep 17 00:00:00 2001 From: DanielVandH Date: Sun, 7 Jul 2024 21:38:26 +0100 Subject: [PATCH 06/14] Document the disabling of ExactPredicates; add the triangulation validator into the package --- docs/make.jl | 1 + docs/src/api/triangulation.md | 1 + docs/src/extended/utils.md | 6 + docs/src/index.md | 2 +- docs/src/manual/disabling_ea.md | 90 +++++ docs/src/manual/predicates.md | 4 +- docs/src/math/clipped.md | 2 +- src/DelaunayTriangulation.jl | 2 + src/README.md | 1 - src/predicates/exactpredicates_definitions.jl | 2 +- .../validation.jl | 318 ++++++++++-------- test/helper_functions.jl | 2 +- 12 files changed, 275 insertions(+), 156 deletions(-) create mode 100644 docs/src/manual/disabling_ea.md delete mode 100644 src/README.md rename test/triangulation_validation.jl => src/validation.jl (78%) diff --git a/docs/make.jl b/docs/make.jl index e18dd1787..8018f5bc5 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -146,6 +146,7 @@ const _PAGES = [ "Geometrical Predicates" => "manual/predicates.md", "Triangulation Output" => "manual/triangulation_output.md", "Voronoi Tessellation Output" => "manual/voronoi_output.md", + "Disabling Exact Predicates" => "manual/disabling_ea.md", ], "API Reference" => [ "Section Overview" => "api/overview.md", diff --git a/docs/src/api/triangulation.md b/docs/src/api/triangulation.md index a49e35873..2896d9f61 100644 --- a/docs/src/api/triangulation.md +++ b/docs/src/api/triangulation.md @@ -59,4 +59,5 @@ get_insertion_order get_vertices(::Triangulation) get_edges(::Triangulation) num_neighbours +validate_triangulation ``` \ No newline at end of file diff --git a/docs/src/extended/utils.md b/docs/src/extended/utils.md index 3e00efb44..99505482b 100644 --- a/docs/src/extended/utils.md +++ b/docs/src/extended/utils.md @@ -28,4 +28,10 @@ Pages = ["src/predicates/boundaries_and_ghosts.jl"] ```@docs convert_certificate +DefaultAdjacentValue +𝒢 +GhostVertex +ε +∅ +USE_EXACTPREDICATE ``` \ No newline at end of file diff --git a/docs/src/index.md b/docs/src/index.md index 63db4fd18..6ca4a2732 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -18,7 +18,7 @@ This is a package for computing Delaunay triangulations and Voronoi tessellation - Computation of the pole of inaccessibility. - The interface for defining geometric primitives is fully customisable. -To ensure that the algorithms are robust, we use [ExactPredicates.jl](https://github.com/lairez/ExactPredicates.jl) to define all geometric predicates in this package. Much of the work in this package is derived from the book [*Delaunay Mesh Generation* by Cheng, Dey, and Shewchuk (2013)](https://people.eecs.berkeley.edu/~jrs/meshbook.html). +To ensure that the algorithms are robust, we use [ExactPredicates.jl](https://github.com/lairez/ExactPredicates.jl) to define all geometric predicates in this package. (ExactPredicates.jl can be disabled, as described [here](manual/disabling_ea.md).) Much of the work in this package is derived from the book [*Delaunay Mesh Generation* by Cheng, Dey, and Shewchuk (2013)](https://people.eecs.berkeley.edu/~jrs/meshbook.html). # Documentation Structure diff --git a/docs/src/manual/disabling_ea.md b/docs/src/manual/disabling_ea.md new file mode 100644 index 000000000..946ef54ef --- /dev/null +++ b/docs/src/manual/disabling_ea.md @@ -0,0 +1,90 @@ +```@meta +CurrentModule = DelaunayTriangulation +``` + +# Disabling Exact Predicates + +For performance reasons, you may find it useful to want to disable exact predicates using [ExactPredicates.jl](https://github.com/lairez/ExactPredicates.jl). This can be easily done using a setup with [Preferences.jl](https://github.com/JuliaPackaging/Preferences.jl), but before you consider disabling exact predicates.jl, there are a few things to be aware of. If you just want to disable them without reading a lot of information warning you about the consequences, pleas skip to the end. + +## Why use exact predicates? + +Three great resources for understanding why we need exact predicates are + +1. Jonathan Shewchuk's paper on adaptive precision floating-point arithmetic [here](https://doi.org/10.1007/PL00009321). +2. Jonathan Shewchuk's lecture notes on geometric robustness [here](https://people.eecs.berkeley.edu/~jrs/meshpapers/robnotes.pdf). +3. [This paper](https://doi.org/10.1016/j.comgeo.2007.06.003) by Kettner et al. (2008) on some examples of issues with inexact arithmetic. + +We give a simple summary here. A big component of the algorithms used in this package are what are known as _geometric predicates_, some of these being: + +- `orient(p, q, r)`: Is `r` left, right, or on the line through `pq`? +- `incircle(p, q, r, s)`: Is `s` inside, outside, or on the circle through `p`, `q`, and `r`? + +These predicates can be computed using determinants: + +```math +\begin{align*} +O_{pqr} &:= \textrm{orient}(p, q, r) = \begin{vmatrix} p_x - r_x & p_y - r_y \\ q_x - r_x & q_y - r_y \end{vmatrix}, \\ +C_{pqrs} &:= \textrm{incircle}(p, q, r, s) = \begin{vmatrix} p_x - s_x & p_y - s_y & (p_x - s_x)^2 + (p_y - s_y)^2 \\ q_x - s_x & q_y - s_y & (q_x - s_x)^2 + (q_y - s_y)^2 \\ r_x - s_x & r_y - s_y & (r_x - s_x)^2 + (r_y - s_y)^2 \end{vmatrix}. +\end{align*} +``` + +The signs of these determinants $O_{pqr}$ and $C_{pqrs}$ are used to determinant the answers to the above questions. In inexact arithmetic, it is common that the sign picked is wrong when the determinants are close to zero. The consequences of this can be catastrophic: + +1. The algorithms may hang or crash. +2. The final triangulation may be completely invalid. For example, if a point is being added into a triangulation right onto an existing edge, then in exact arithmetic we would know to split the edge to the left and to the right. In inexact arithmetic, the point may be to the left of the edge but detected as being to the right of it, thus adding a triangle that crosses an edge. +3. You may encounter `BoundsError`s from bad `Adjacent` queries where a triangle is expected to exist but doesn't. + +Another issue is due to the fact that floating point arithmetic is not associative. In exact arithmetic, we would expect for example that + +```math +O_{pqr} = O_{qrp} = O_{rpq}, +``` + +but this is not true in floating point arithmetic. This causes issues with consistency - a point may be found to be both left and right of an edge depending on the order of the points given to the `orient` predicate, inevitably leading to an invalid triangulation. With the use of exact predicates, this property is guaranteed to hold, ensuring that all the predicate results are consistent with each other. This has the following consequence: **Even if you think exact predicates are not necessary for you because none of your inputs are exact (for example), you still want them to guarantee consistency with predicates regardless of the input order**. + +## Will disabling exact predicates give me better performance? + +It is also not even the case that using inexact predicates will give you better performance than if you were to us exact predicates. ExactPredicates.jl uses clever filters to ensure that the fast path, meaning an inexact computation, is performed first. This means that, in most cases, computing the `orient` predicate would have the same speed as if you were to compute it using the determinant definition. + +One case where you could see improvements would be if there were many collinear points, where the slow computations from exact predicates would actually be needed - but this is exactly the case where you expect inexact predicates to cause you problems! + +You should always benchmark your problems to see if disabling exact predicates, if you choose to do, will actually give you better performance. + +## How do I disable exact predicates? + +If you still want to disable exact predicates, here is how you can do so. Before doing `using DelaunayTriangulation`, you can use + +```julia +using Preferences: set_preferences! +set_preferences!("DelaunayTriangulation", "USE_EXACTPREDICATES" => false) +using DelaunayTriangulation # only load after setting the preference +``` + +The `set_preferences!` call will make a `LocalPreferences.toml` file in your directory that sets this preference. Once this file exists you are free to delete `Preferences.jl`. If you want to skip using Preferences.jl entirely, you can also just create `LocalPreferences.toml` in your working directory manually and put + +``` +[DelaunayTriangulation] +USE_EXACTPREDICATES = false +``` + +into it. If you later want to re-enable exact predicates, either delete the file or write `USE_EXACTPREDICATES = true` instead. This setup ensures that there is no slowdown in the package from checking if ExactPredicates.jl is being used at runtime, as it is all done at compile time instead. + +## Can I check if my computed triangulation is valid? + +When you are not using exact predicates, you may want to check if your computed triangulation is actually a valid Delaunay triangulation. We provide the function `DelaunayTriangulation.validate_triangulation` for this purpose. This functionality is quite slow to use and is not currently optimised or well-documented (contributions towards addressing these issues are welcome), but it will work. One important note is that this check does actually use predicates in certain areas, so this check is still not guaranteed to be 100% accurate. Here is an example of its use. + +```@example validate +using DelaunayTriangulation +tri = triangulate(rand(2, 50)) +DelaunayTriangulation.validate_triangulation(tri) +``` + +```@example validate +T = first(each_solid_triangle(tri)) +DelaunayTriangulation.delete_triangle!(tri, T) # break the triangulation for this example +DelaunayTriangulation.validate_triangulation(tri) +``` + +```@example validate +DelaunayTriangulation.validate_triangulation(tri; print_result = false) +``` diff --git a/docs/src/manual/predicates.md b/docs/src/manual/predicates.md index 4d77271e8..5d4e862f6 100644 --- a/docs/src/manual/predicates.md +++ b/docs/src/manual/predicates.md @@ -4,7 +4,7 @@ CurrentModule = DelaunayTriangulation # Geometrical Predicates -This section discusses how geometrical predicates are defined in this package. The predicates in this package are primarily derived from those implemented in [ExactPredicates.jl](https://github.com/lairez/ExactPredicates.jl). The choice of exact predicates is important for the robustness of the algorithms in this package. Without using exact predicates, you may quickly find issues such as infinite loops or errors in the algorithms, as discussed for example at the start of p.3 of [these notes](https://perso.uclouvain.be/jean-francois.remacle/LMECA2170/robnotes.pdf) by Shewchuk. +This section discusses how geometrical predicates are defined in this package. The predicates in this package are primarily derived from those implemented in [ExactPredicates.jl](https://github.com/lairez/ExactPredicates.jl). The choice of exact predicates is important for the robustness of the algorithms in this package. Without using exact predicates, you may quickly find issues such as infinite loops or errors in the algorithms, as discussed for example at the start of p.3 of [these notes](https://perso.uclouvain.be/jean-francois.remacle/LMECA2170/robnotes.pdf) by Shewchuk. If you do want to disable exact predicates, see [here](disabling_ea.md). ## Certificates @@ -24,4 +24,4 @@ DelaunayTriangulation.is_positively_oriented(flag) ## Predicates -The predicates we define, as mentioned, are primarily derived from those in [ExactPredicates.jl](https://github.com/lairez/ExactPredicates.jl), where we simply extend the [`orient`](@ref orient_predicate), [`incircle`](@ref incircle_predicate), [`parallelorder`](@ref parallelorder_predicate), [`sameside`](@ref sameside_predicate), and [`meet`](@ref meet_predicate) predicates, allowing us to predicates for, for example, [`triangle_orientation`](@ref) and [`point_position_relative_to_line`](@ref). Predicates for working with the boundary and ghost vertices are also implemented, for example, [`is_ghost_edge`](@ref is_ghost_edge) and [`is_boundary_node`](@ref). \ No newline at end of file +The predicates we define, as mentioned, are primarily derived from those in [ExactPredicates.jl](https://github.com/lairez/ExactPredicates.jl), where we simply extend the [`orient`](@ref orient_predicate), [`incircle`](@ref incircle_predicate), [`parallelorder`](@ref parallelorder_predicate), [`sameside`](@ref sameside_predicate), and [`meet`](@ref meet_predicate) predicates, allowing us to define predicates for, for example, [`triangle_orientation`](@ref) and [`point_position_relative_to_line`](@ref). Predicates for working with the boundary and ghost vertices are also implemented, for example, [`is_ghost_edge`](@ref is_ghost_edge) and [`is_boundary_node`](@ref). \ No newline at end of file diff --git a/docs/src/math/clipped.md b/docs/src/math/clipped.md index 3108fe181..3af108345 100644 --- a/docs/src/math/clipped.md +++ b/docs/src/math/clipped.md @@ -318,7 +318,7 @@ new_points[u] = p new_points[v] = q new_vertices[u] = u new_vertices[v] = v -poly!(ax, [new_points; new_points[1]], color = :darkgreen, linewidth = 3) +poly!(ax, [new_points; new_points[1]], color = :darkgreen) lines!(ax, [0.0, 5.0, 5.0, 0.0, 0.0], [-15.0, -15.0, 15.0, 15.0, -15.0], color = :red) resize_to_layout!(fig) ``` diff --git a/src/DelaunayTriangulation.jl b/src/DelaunayTriangulation.jl index ccc052e97..491ed8f06 100644 --- a/src/DelaunayTriangulation.jl +++ b/src/DelaunayTriangulation.jl @@ -185,6 +185,8 @@ include("algorithms/voronoi/clipped.jl") include("algorithms/voronoi/main.jl") include("algorithms/voronoi/unbounded.jl") +include("validation.jl") + export each_triangle, each_solid_triangle, diff --git a/src/README.md b/src/README.md deleted file mode 100644 index 8b1378917..000000000 --- a/src/README.md +++ /dev/null @@ -1 +0,0 @@ - diff --git a/src/predicates/exactpredicates_definitions.jl b/src/predicates/exactpredicates_definitions.jl index 496ea5987..ea2824dc7 100644 --- a/src/predicates/exactpredicates_definitions.jl +++ b/src/predicates/exactpredicates_definitions.jl @@ -183,7 +183,7 @@ without exact arithmetic. The incircle predicate is defined by the determinant ```math -\\text{incircle}(a, b, c, d) = \\text{sgn} \\det \\begin{bmatrix} a_x & a_y & a_x^2 + a_y^2 & 1 \\\\ b_x & b_y & b_x62 + b_y^2 & 1 \\\\ c_x & c_y & c_x^2 + c_y^2 & 1 \\\\ d_x & d_y & d_x^2 + d_y^2 & 1 \\end{bmatrix} = \\text{sgn} \\det \\begin{bmatrix} a_x - d_x & a_y - d_y & (a_x - d_x)^2 + (a_y - d_y)^2 \\\\ b_x - d_x & b_y - d_y & (b_x - d_x)^2 + (b_y - d_y)^2 \\\\ c_x - d_x & c_y - d_y & (c_x - d_x)^2 + (c_y - d_y)^2 \\end{bmatrix}. +\\text{incircle}(a, b, c, d) = \\text{sgn} \\det \\begin{bmatrix} a_x & a_y & a_x^2 + a_y^2 & 1 \\\\ b_x & b_y & b_x^2 + b_y^2 & 1 \\\\ c_x & c_y & c_x^2 + c_y^2 & 1 \\\\ d_x & d_y & d_x^2 + d_y^2 & 1 \\end{bmatrix} = \\text{sgn} \\det \\begin{bmatrix} a_x - d_x & a_y - d_y & (a_x - d_x)^2 + (a_y - d_y)^2 \\\\ b_x - d_x & b_y - d_y & (b_x - d_x)^2 + (b_y - d_y)^2 \\\\ c_x - d_x & c_y - d_y & (c_x - d_x)^2 + (c_y - d_y)^2 \\end{bmatrix}. ``` """ incircle_predicate diff --git a/test/triangulation_validation.jl b/src/validation.jl similarity index 78% rename from test/triangulation_validation.jl rename to src/validation.jl index ba1cea877..7fa78b6f0 100644 --- a/test/triangulation_validation.jl +++ b/src/validation.jl @@ -1,5 +1,3 @@ -# maybe this should belong inside src for users at some point? - #= """ abstract type AbstractTriangulationState @@ -19,6 +17,22 @@ abstract type AbstractTriangulationState end test_state(state::AbstractTriangulationState) = state.flag Base.show(io::IO, state::AbstractTriangulationState) = print(io, summary(state)) +function compare_edge_vectors(E1, E2) + E1s = sort_edge_vector(collect(E1)) + E2s = sort_edge_vector(collect(E2)) + return E1s == E2s +end + +function sort_edge_vector(E) + sorted_E = similar(E) + for i in eachindex(E) + u, v = E[i] + e = (min(u, v), max(u, v)) + sorted_E[i] = e + end + return sort(sorted_E) +end + struct TriangleOrientationState <: AbstractTriangulationState flag::Bool bad_triangle::NTuple{3,Int} @@ -34,12 +48,12 @@ function Base.summary(state::TriangleOrientationState) end function test_triangle_orientation(tri) for T in each_solid_triangle(tri) - cert = DT.triangle_orientation(tri, T) - flag = DT.is_positively_oriented(cert) - orientation = DT.is_positively_oriented(cert) ? :positively : DT.is_negatively_oriented(cert) ? :negatively : :degenerately + cert = triangle_orientation(tri, T) + flag = is_positively_oriented(cert) + orientation = is_positively_oriented(cert) ? :positively : is_negatively_oriented(cert) ? :negatively : :degenerately !flag && return TriangleOrientationState(flag, Int.(triangle_vertices(T)), orientation) end - return TriangleOrientationState(true, (DT.∅, DT.∅, DT.∅), :positive) + return TriangleOrientationState(true, (∅, ∅, ∅), :positive) end struct DelaunayCriterionState <: AbstractTriangulationState @@ -52,7 +66,7 @@ function Base.summary(state::DelaunayCriterionState) if test_state(state) return "All the triangles are Delaunay." else - if state.bad_vertex !== DT.∅ + if state.bad_vertex !== ∅ return "The Delaunay criterion does not hold for the triangle-vertex pair ($(state.bad_triangle), $(state.bad_vertex))." else return "The test of the Delaunay criterion failed as there was a BoundsError when testing the visibility." @@ -62,23 +76,23 @@ end function test_delaunay_criterion(tri) try points = get_points(tri) - triangle_tree = DT.BoundaryRTree(points) - segment_tree = DT.BoundaryRTree(points) - failures = Tuple{DT.triangle_type(tri),DT.integer_type(tri)}[] + triangle_tree = BoundaryRTree(points) + segment_tree = BoundaryRTree(points) + failures = Tuple{triangle_type(tri),integer_type(tri)}[] for T in each_solid_triangle(tri) - i, j, k = DT.triangle_vertices(T) + i, j, k = triangle_vertices(T) p, q, r = get_point(tri, i, j, k) - c = DT.triangle_circumcenter(p, q, r) - cr = DT.triangle_circumradius(p, q, r) + c = triangle_circumcenter(p, q, r) + cr = triangle_circumradius(p, q, r) xmin, xmax = getx(c) - cr, getx(c) + cr ymin, ymax = gety(c) - cr, gety(c) + cr any(!isfinite, (xmin, xmax, ymin, ymax)) && continue - bbox = DT.BoundingBox(xmin, xmax, ymin, ymax) - dbx = DT.DiametralBoundingBox(bbox, (i, j)) # just store (i, j) and get k via get_adjacent + bbox = BoundingBox(xmin, xmax, ymin, ymax) + dbx = DiametralBoundingBox(bbox, (i, j)) # just store (i, j) and get k via get_adjacent insert!(triangle_tree.tree, dbx) end for e in each_segment(tri) - i, j = DT.edge_vertices(e) + i, j = edge_vertices(e) p, q = get_point(tri, i, j) px, py = getxy(p) qx, qy = getxy(q) @@ -86,32 +100,32 @@ function test_delaunay_criterion(tri) xmax = max(px, qx) ymin = min(py, qy) ymax = max(py, qy) - bbox = DT.BoundingBox(xmin, xmax, ymin, ymax) - dbx = DT.DiametralBoundingBox(bbox, (i, j)) + bbox = BoundingBox(xmin, xmax, ymin, ymax) + dbx = DiametralBoundingBox(bbox, (i, j)) insert!(segment_tree.tree, dbx) end for r in shuffle(collect(each_solid_vertex(tri))) !isempty(failures) && break - intersects = DT.get_intersections(triangle_tree, r, cache_id=1) # can't use multithreading here + intersects = get_intersections(triangle_tree, r, cache_id=1) # can't use multithreading here for box in intersects - i, j = DT.get_edge(box) + i, j = get_edge(box) k = get_adjacent(tri, i, j) any(==(r), (i, j, k)) && continue - T = DT.construct_triangle(DT.triangle_type(tri), i, j, k) - cert = DT.point_position_relative_to_circumcircle(tri, T, r) - c = DT.triangle_centroid(get_point(tri, i, j, k)...) - if DT.is_inside(cert) - DT.is_boundary_edge(tri, i, j) && DT.is_right(DT.point_position_relative_to_line(tri, i, j, r)) && continue # if it's outside of the domain relative to this edge, just continue - DT.is_boundary_edge(tri, j, k) && DT.is_right(DT.point_position_relative_to_line(tri, j, k, r)) && continue - DT.is_boundary_edge(tri, k, i) && DT.is_right(DT.point_position_relative_to_line(tri, k, i, r)) && continue - for (i, j) in DT.triangle_edges(i, j, k) - DT.contains_segment(tri, i, j) && continue # visibility is defined according to the relative interior of the simplex, which means that it's fine if a segment can see the vertex + T = construct_triangle(triangle_type(tri), i, j, k) + cert = point_position_relative_to_circumcircle(tri, T, r) + c = triangle_centroid(get_point(tri, i, j, k)...) + if is_inside(cert) + is_boundary_edge(tri, i, j) && is_right(point_position_relative_to_line(tri, i, j, r)) && continue # if it's outside of the domain relative to this edge, just continue + is_boundary_edge(tri, j, k) && is_right(point_position_relative_to_line(tri, j, k, r)) && continue + is_boundary_edge(tri, k, i) && is_right(point_position_relative_to_line(tri, k, i, r)) && continue + for (i, j) in triangle_edges(i, j, k) + contains_segment(tri, i, j) && continue # visibility is defined according to the relative interior of the simplex, which means that it's fine if a segment can see the vertex if rand() < 1 / 2 # just testing both - cert = DT.test_visibility(tri, i, j, r, shift=0.01, attractor=c) + cert = test_visibility(tri, i, j, r, shift=0.01, attractor=c) else cert = test_visibility(tri, segment_tree, i, j, r, c) end - flag = DT.is_visible(cert) + flag = is_visible(cert) flag && push!(failures, (T, r)) flag && break end @@ -119,35 +133,35 @@ function test_delaunay_criterion(tri) end end if isempty(failures) - return DelaunayCriterionState(true, (DT.∅, DT.∅, DT.∅), DT.∅) + return DelaunayCriterionState(true, (∅, ∅, ∅), ∅) else return DelaunayCriterionState(false, Int.(failures[1][1]), Int(failures[1][2])) end catch e - e isa BoundsError && return DelaunayCriterionState(false, (DT.∅, DT.∅, DT.∅), DT.∅) + e isa BoundsError && return DelaunayCriterionState(false, (∅, ∅, ∅), ∅) rethrow(e) end end function test_visibility(tri::Triangulation, segment_tree, i, j, k, centroid=nothing) - if DT.is_boundary_edge(tri, j, i) + if is_boundary_edge(tri, j, i) i, j = j, i end p, q, a = get_point(tri, i, j, k) - if DT.is_boundary_edge(tri, i, j) - side_e = DT.point_position_relative_to_line(p, q, a) - DT.is_right(side_e) && return DT.Cert.Invisible + if is_boundary_edge(tri, i, j) + side_e = point_position_relative_to_line(p, q, a) + is_right(side_e) && return Cert.Invisible end # Need to see if i or j is a boundary node without the other being a boundary node. # If this is the case, it's possible that one of them mistakenly sees k. for u in (i, j) - flag, g = DT.is_boundary_node(tri, u) + flag, g = is_boundary_node(tri, u) if flag - ℓ = DT.get_left_boundary_node(tri, u, g) - cert = DT.point_position_relative_to_line(tri, ℓ, u, k) - DT.is_right(cert) && return DT.Cert.Invisible - ℓ = DT.get_right_boundary_node(tri, u, g) - cert = DT.point_position_relative_to_line(tri, u, ℓ, k) - DT.is_right(cert) && return DT.Cert.Invisible + ℓ = get_left_boundary_node(tri, u, g) + cert = point_position_relative_to_line(tri, ℓ, u, k) + is_right(cert) && return Cert.Invisible + ℓ = get_right_boundary_node(tri, u, g) + cert = point_position_relative_to_line(tri, u, ℓ, k) + is_right(cert) && return Cert.Invisible end end flags = falses(10) @@ -163,18 +177,17 @@ function test_visibility(tri::Triangulation, segment_tree, i, j, k, centroid=not mx, my = getxy(m) ax, ay = getxy(a) xmin, ymin, xmax, ymax = min(mx, ax), min(my, ay), max(mx, ax), max(my, ay) - bbox = DT.BoundingBox(xmin, xmax, ymin, ymax) - intersections = DT.get_intersections(segment_tree, bbox; cache_id=2) + bbox = BoundingBox(xmin, xmax, ymin, ymax) + intersections = get_intersections(segment_tree, bbox; cache_id=2) for box in intersections - u, v = DT.get_edge(box) - # !DT.edges_are_disjoint((i, j), (u, v)) && continue + u, v = get_edge(box) p′, q′ = get_point(tri, u, v) - cert = DT.line_segment_intersection_type(m, a, p′, q′) - flags[idx] = !DT.has_no_intersections(cert) && !DT.is_touching(cert) + cert = line_segment_intersection_type(m, a, p′, q′) + flags[idx] = !has_no_intersections(cert) && !is_touching(cert) flags[idx] && break end end - return all(flags) ? DT.Cert.Invisible : DT.Cert.Visible + return all(flags) ? Cert.Invisible : Cert.Visible end struct EdgesHaveTwoIncidentTrianglesState <: AbstractTriangulationState @@ -191,22 +204,22 @@ function Base.summary(state::EdgesHaveTwoIncidentTrianglesState) end function test_each_edge_has_two_incident_triangles(tri) for e in each_edge(tri) - i, j = DT.edge_vertices(e) + i, j = edge_vertices(e) vᵢⱼ = get_adjacent(tri, i, j) vⱼᵢ = get_adjacent(tri, j, i) - if DT.is_boundary_edge(tri, j, i) - flag = DT.is_ghost_vertex(vᵢⱼ) && DT.edge_exists(vⱼᵢ) + if is_boundary_edge(tri, j, i) + flag = is_ghost_vertex(vᵢⱼ) && edge_exists(vⱼᵢ) !flag && return EdgesHaveTwoIncidentTrianglesState(flag, Int.(edge_vertices(e))) - elseif DT.is_boundary_edge(tri, i, j) - flag = DT.is_ghost_vertex(vⱼᵢ) && DT.edge_exists(vᵢⱼ) + elseif is_boundary_edge(tri, i, j) + flag = is_ghost_vertex(vⱼᵢ) && edge_exists(vᵢⱼ) !flag && return EdgesHaveTwoIncidentTrianglesState(flag, Int.(edge_vertices(e))) else - flag = DT.edge_exists(vᵢⱼ) && DT.edge_exists(vⱼᵢ) + flag = edge_exists(vᵢⱼ) && edge_exists(vⱼᵢ) !flag && return EdgesHaveTwoIncidentTrianglesState(flag, Int.(edge_vertices(e))) end end - DT.clear_empty_features!(tri) - return EdgesHaveTwoIncidentTrianglesState(true, (DT.∅, DT.∅)) + clear_empty_features!(tri) + return EdgesHaveTwoIncidentTrianglesState(true, (∅, ∅)) end struct AdjacentMapState <: AbstractTriangulationState @@ -225,7 +238,7 @@ function Base.summary(state::AdjacentMapState) end function test_adjacent_map_matches_triangles(tri) for T in each_triangle(tri) - u, v, w = DT.triangle_vertices(T) + u, v, w = triangle_vertices(T) flag1 = get_adjacent(tri, u, v) == w !flag1 && return AdjacentMapState(flag1, Int.(triangle_vertices(T)), Int.((u, v)), get_adjacent(tri, u, v)) flag2 = get_adjacent(tri, v, w) == u @@ -233,8 +246,8 @@ function test_adjacent_map_matches_triangles(tri) flag3 = get_adjacent(tri, w, u) == v !flag3 && return AdjacentMapState(flag3, Int.(triangle_vertices(T)), Int.((w, u)), get_adjacent(tri, w, u)) end - DT.clear_empty_features!(tri) - return AdjacentMapState(true, (DT.∅, DT.∅, DT.∅), (DT.∅, DT.∅), DT.∅) + clear_empty_features!(tri) + return AdjacentMapState(true, (∅, ∅, ∅), (∅, ∅), ∅) end struct Adjacent2VertexMapState <: AbstractTriangulationState @@ -253,20 +266,20 @@ function Base.summary(state::Adjacent2VertexMapState) end end function test_adjacent2vertex_map_matches_triangles(tri) - E = DT.edge_type(tri) + E = edge_type(tri) for T in each_triangle(tri) u, v, w = triangle_vertices(T) for (u, v, w) in ((u, v, w), (v, w, u), (w, u, v)) - vw = DT.construct_edge(E, v, w) + vw = construct_edge(E, v, w) Su = get_adjacent2vertex(tri, u) - flag = DT.contains_edge(vw, Su) + flag = contains_edge(vw, Su) if !flag - _Su = Set{NTuple{2,Int}}(Int.(DT.edge_vertices(e)) for e in each_edge(Su)) - return Adjacent2VertexMapState(flag, Int.(DT.triangle_vertices(T)), Int.((v, w)), Int(u), _Su) + _Su = Set{NTuple{2,Int}}(Int.(edge_vertices(e)) for e in each_edge(Su)) + return Adjacent2VertexMapState(flag, Int.(triangle_vertices(T)), Int.((v, w)), Int(u), _Su) end end end - return Adjacent2VertexMapState(true, (DT.∅, DT.∅, DT.∅), (DT.∅, DT.∅), DT.∅, Set{NTuple{2,Int}}()) + return Adjacent2VertexMapState(true, (∅, ∅, ∅), (∅, ∅), ∅, Set{NTuple{2,Int}}()) end struct AdjacentMapAdjacent2VertexMapState <: AbstractTriangulationState @@ -286,14 +299,14 @@ function Base.summary(state::AdjacentMapAdjacent2VertexMapState) end function test_adjacent_map_matches_adjacent2vertex_map(tri) for (k, S) in get_adjacent2vertex(get_adjacent2vertex(tri)) - _S = Set{NTuple{2,Int}}(Int.(DT.edge_vertices(e)) for e in each_edge(S)) + _S = Set{NTuple{2,Int}}(Int.(edge_vertices(e)) for e in each_edge(S)) for e in each_edge(S) flag = get_adjacent(tri, e) == k !flag && return AdjacentMapAdjacent2VertexMapState(flag, Int(k), get_adjacent(tri, e), _S, Int.(edge_vertices(e))) end end - DT.clear_empty_features!(tri) - return AdjacentMapAdjacent2VertexMapState(true, DT.∅, DT.∅, Set{NTuple{2,Int}}(), (DT.∅, DT.∅)) + clear_empty_features!(tri) + return AdjacentMapAdjacent2VertexMapState(true, ∅, ∅, Set{NTuple{2,Int}}(), (∅, ∅)) end struct Adjacent2VertexMapAdjacentMapState <: AbstractTriangulationState @@ -320,14 +333,14 @@ function test_adjacent2vertex_map_matches_adjacent_map(tri) k = get_adjacent(tri, e) if haskey(adj2v, k) S = get_adjacent2vertex(tri, k) - flag = DT.contains_edge(e, S) + flag = contains_edge(e, S) else flag = false end !flag && return Adjacent2VertexMapAdjacentMapState(flag, Int(k), Int.(edge_vertices(e)), haskey(adj2v, k)) end - DT.clear_empty_features!(tri) - return Adjacent2VertexMapAdjacentMapState(true, DT.∅, (DT.∅, DT.∅), true) + clear_empty_features!(tri) + return Adjacent2VertexMapAdjacentMapState(true, ∅, (∅, ∅), true) end struct GraphState <: AbstractTriangulationState @@ -345,14 +358,14 @@ end function test_graph_contains_all_vertices(tri) all_vertices = Set{Int}() for T in each_triangle(tri) # need a method that doesn't use the graph - i, j, k = DT.triangle_vertices(T) + i, j, k = triangle_vertices(T) push!(all_vertices, i, j, k) end for i in all_vertices - flag = DT.has_vertex(tri, i) + flag = has_vertex(tri, i) !flag && return GraphState(flag, i) end - return GraphState(true, DT.∅) + return GraphState(true, ∅) end struct GraphAdjacentMapState <: AbstractTriangulationState @@ -368,23 +381,23 @@ function Base.summary(state::GraphAdjacentMapState) end end function test_graph_matches_adjacent_map(tri) - E = DT.edge_type(tri) + E = edge_type(tri) adj_dict = get_adjacent(get_adjacent(tri)) for e in each_edge(tri) - i, j = DT.edge_vertices(e) - if DT.has_ghost_triangles(tri) || !DT.is_ghost_edge(i, j) - eᵢⱼ = DT.construct_edge(E, i, j) - eⱼᵢ = DT.construct_edge(E, j, i) - flag = if !DT.has_multiple_sections(tri) + i, j = edge_vertices(e) + if has_ghost_triangles(tri) || !is_ghost_edge(i, j) + eᵢⱼ = construct_edge(E, i, j) + eⱼᵢ = construct_edge(E, j, i) + flag = if !has_multiple_sections(tri) eᵢⱼ ∈ keys(adj_dict) && eⱼᵢ ∈ keys(adj_dict) else - DT.edge_exists(tri, i, j) && DT.edge_exists(tri, j, i) + edge_exists(tri, i, j) && edge_exists(tri, j, i) end !flag && return GraphAdjacentMapState(flag, Int.(edge_vertices(e))) end end - DT.clear_empty_features!(tri) - return GraphAdjacentMapState(true, (DT.∅, DT.∅)) + clear_empty_features!(tri) + return GraphAdjacentMapState(true, (∅, ∅)) end struct AdjacentMapGraphState <: AbstractTriangulationState @@ -409,7 +422,7 @@ function Base.summary(state::AdjacentMapGraphState) end function test_adjacent_map_matches_graph(tri) for (e, k) in get_adjacent(get_adjacent(tri)) - i, j = DT.edge_vertices(e) + i, j = edge_vertices(e) flag1 = all(∈(get_neighbours(tri, i)), (j, k)) !flag1 && return AdjacentMapGraphState(flag1, Int(i), Int.((i, j)), Int.((j, k)), false, Int(k)) flag2 = all(∈(get_neighbours(tri, j)), (k, i)) @@ -417,7 +430,7 @@ function test_adjacent_map_matches_graph(tri) flag3 = all(∈(get_neighbours(tri, k)), (i, j)) !flag3 && return AdjacentMapGraphState(flag3, Int(k), Int.((i, j)), Int.((i, j)), true, Int(k)) end - return AdjacentMapGraphState(true, DT.∅, (DT.∅, DT.∅), (DT.∅, DT.∅), true, DT.∅) + return AdjacentMapGraphState(true, ∅, (∅, ∅), (∅, ∅), true, ∅) end struct GraphTrianglesState <: AbstractTriangulationState @@ -431,7 +444,7 @@ function Base.summary(state::GraphTrianglesState) if test_state(state) return "The graph correctly matches the triangle set." else - if state.bad_vertex == DT.∅ + if state.bad_vertex == ∅ return "The graph is inconsistent with the triangle set, as one of the vertices of $(state.bad_triangle) is not a vertex in the graph." else return "The graph is inconsistent with the triangle set. The triangle $(state.bad_triangle) is in the triangle set but either $(state.bad_edge[1]) or $(state.bad_edge[2]) are not in $(state.bad_vertex)'s neighbourhood." @@ -441,7 +454,7 @@ end function test_graph_matches_triangles(tri) for T in each_triangle(tri) try - i, j, k = DT.triangle_vertices(T) + i, j, k = triangle_vertices(T) flag1 = all(∈(get_neighbours(tri, i)), (j, k)) !flag1 && return GraphTrianglesState(flag1, Int.(triangle_vertices(T)), Int.((j, k)), Int(i)) flag2 = all(∈(get_neighbours(tri, j)), (k, i)) @@ -449,11 +462,11 @@ function test_graph_matches_triangles(tri) flag3 = all(∈(get_neighbours(tri, k)), (i, j)) !flag3 && return GraphTrianglesState(flag3, Int.(triangle_vertices(T)), Int.((i, j)), Int(k)) catch e - e isa KeyError && return GraphTrianglesState(false, Int.(triangle_vertices(T)), (DT.∅, DT.∅), DT.∅) + e isa KeyError && return GraphTrianglesState(false, Int.(triangle_vertices(T)), (∅, ∅), ∅) rethrow(e) end end - return GraphTrianglesState(true, (DT.∅, DT.∅, DT.∅), (DT.∅, DT.∅), DT.∅) + return GraphTrianglesState(true, (∅, ∅, ∅), (∅, ∅), ∅) end struct SegmentState <: AbstractTriangulationState @@ -477,31 +490,31 @@ function Base.summary(state::SegmentState) end function test_segments(tri) for e in each_segment(tri) - flag = DT.edge_exists(tri, e) || DT.edge_exists(tri, DT.reverse_edge(e)) + flag = edge_exists(tri, e) || edge_exists(tri, reverse_edge(e)) !flag && return SegmentState(flag, Int.(edge_vertices(e)), 0) - flag = DT.contains_edge(e, get_interior_segments(tri)) || DT.contains_edge(DT.reverse_edge(e), get_interior_segments(tri)) + flag = contains_edge(e, get_interior_segments(tri)) || contains_edge(reverse_edge(e), get_interior_segments(tri)) if !flag - flag = DT.contains_boundary_edge(tri, e) || DT.contains_boundary_edge(tri, DT.reverse_edge(e)) + flag = contains_boundary_edge(tri, e) || contains_boundary_edge(tri, reverse_edge(e)) !flag && return SegmentState(flag, Int.(edge_vertices(e)), 1) end end for e in get_interior_segments(tri) - flag = DT.edge_exists(tri, e) || DT.edge_exists(tri, DT.reverse_edge(e)) + flag = edge_exists(tri, e) || edge_exists(tri, reverse_edge(e)) !flag && return SegmentState(flag, Int.(edge_vertices(e)), 0) end int_segs = Set{NTuple{2,Int}}() all_segs = Set{NTuple{2,Int}}() for e in get_interior_segments(tri) - u, v = DT.edge_vertices(e) + u, v = edge_vertices(e) push!(int_segs, minmax(u, v)) end for e in each_segment(tri) - u, v = DT.edge_vertices(e) + u, v = edge_vertices(e) push!(all_segs, minmax(u, v)) end flag = int_segs ⊆ all_segs - !flag && return SegmentState(flag, (DT.∅, DT.∅), 2) - return SegmentState(true, (DT.∅, DT.∅), 0) + !flag && return SegmentState(flag, (∅, ∅), 2) + return SegmentState(true, (∅, ∅), 0) end struct BoundaryEdgeMapBoundaryNodesState <: AbstractTriangulationState @@ -518,15 +531,15 @@ function Base.summary(state::BoundaryEdgeMapBoundaryNodesState) return "The boundary edge map is inconsistent with the boundary nodes. The edge $(state.bad_edge) maps to $(state.bad_pos) in the boundary edge map but $(state.bad_pos) corresponds to the edge $(state.bad_edge_2)." end end -function test_boundary_edge_map_matches_boundary_nodes(tri::DT.Triangulation) - boundary_edge_map = DT.get_boundary_edge_map(tri) +function test_boundary_edge_map_matches_boundary_nodes(tri::Triangulation) + boundary_edge_map = get_boundary_edge_map(tri) for (edge, pos) in boundary_edge_map - u = DT.get_boundary_nodes(DT.get_boundary_nodes(tri, pos[1]), pos[2]) - v = DT.get_boundary_nodes(DT.get_boundary_nodes(tri, pos[1]), pos[2] + 1) - flag = u == DT.initial(edge) + u = get_boundary_nodes(get_boundary_nodes(tri, pos[1]), pos[2]) + v = get_boundary_nodes(get_boundary_nodes(tri, pos[1]), pos[2] + 1) + flag = u == initial(edge) !flag && return BoundaryEdgeMapBoundaryNodesState(flag, Int.(edge_vertices(edge)), Int.((u, v)), pos) end - return BoundaryEdgeMapBoundaryNodesState(true, (DT.∅, DT.∅), (DT.∅, DT.∅), ()) + return BoundaryEdgeMapBoundaryNodesState(true, (∅, ∅), (∅, ∅), ()) end struct BoundaryNodesBoundaryEdgeMapState <: AbstractTriangulationState @@ -543,53 +556,53 @@ function Base.summary(state::BoundaryNodesBoundaryEdgeMapState) return "The boundary nodes are inconsistent with the boundary edge map. The edge $(state.bad_edge) maps to $(state.bad_pos_2) but is at $(state.bad_pos)." end end -function test_boundary_nodes_matches_boundary_edge_map(tri::DT.Triangulation) - boundary_nodes = DT.get_boundary_nodes(tri) - E = DT.edge_type(tri) - if DT.has_multiple_curves(tri) - nc = DT.num_curves(tri) +function test_boundary_nodes_matches_boundary_edge_map(tri::Triangulation) + boundary_nodes = get_boundary_nodes(tri) + E = edge_type(tri) + if has_multiple_curves(tri) + nc = num_curves(tri) for i in 1:nc curve_nodes = get_boundary_nodes(tri, i) - ns = DT.num_sections(curve_nodes) + ns = num_sections(curve_nodes) for j in 1:ns segment_nodes = get_boundary_nodes(curve_nodes, j) - ne = DT.num_boundary_edges(segment_nodes) + ne = num_boundary_edges(segment_nodes) for k in 1:ne left_node = get_boundary_nodes(segment_nodes, k) right_node = get_boundary_nodes(segment_nodes, k + 1) - edge = DT.construct_edge(E, left_node, right_node) - pos = DT.get_boundary_edge_map(tri, edge) + edge = construct_edge(E, left_node, right_node) + pos = get_boundary_edge_map(tri, edge) flag = pos == ((i, j), k) !flag && return BoundaryNodesBoundaryEdgeMapState(flag, Int.(edge_vertices(edge)), pos, ((i, j), k)) end end end - elseif DT.has_multiple_sections(tri) - ns = DT.num_sections(tri) + elseif has_multiple_sections(tri) + ns = num_sections(tri) for i in 1:ns segment_nodes = get_boundary_nodes(tri, i) - ne = DT.num_boundary_edges(segment_nodes) + ne = num_boundary_edges(segment_nodes) for j in 1:ne left_node = get_boundary_nodes(segment_nodes, j) right_node = get_boundary_nodes(segment_nodes, j + 1) - edge = DT.construct_edge(E, left_node, right_node) - pos = DT.get_boundary_edge_map(tri, edge) + edge = construct_edge(E, left_node, right_node) + pos = get_boundary_edge_map(tri, edge) flag = pos == (i, j) !flag && return BoundaryNodesBoundaryEdgeMapState(flag, Int.(edge_vertices(edge)), pos, (i, j)) end end else - ne = DT.num_boundary_edges(boundary_nodes) + ne = num_boundary_edges(boundary_nodes) for i in 1:ne left_node = get_boundary_nodes(boundary_nodes, i) right_node = get_boundary_nodes(boundary_nodes, i + 1) - edge = DT.construct_edge(E, left_node, right_node) - pos = DT.get_boundary_edge_map(tri, edge) + edge = construct_edge(E, left_node, right_node) + pos = get_boundary_edge_map(tri, edge) flag = pos == (get_boundary_nodes(tri), i) !flag && return BoundaryNodesBoundaryEdgeMapState(flag, Int.(edge_vertices(edge)), pos, (get_boundary_nodes(tri), i)) end end - return BoundaryNodesBoundaryEdgeMapState(true, (DT.∅, DT.∅), (), ()) + return BoundaryNodesBoundaryEdgeMapState(true, (∅, ∅), (), ()) end abstract type AbstractTriangulationIteratorState <: AbstractTriangulationState end @@ -687,8 +700,8 @@ end function GhostEdgeIteratorState(unique_flag, length_flag, output_flag, iterator_length, correct_length) return IteratorState(unique_flag, length_flag, output_flag, iterator_length, correct_length, :ghost_edge, :edge, :each_ghost_edge) end -function test_iterators(tri::DT.Triangulation) - I = DT.integer_type(tri) +function test_iterators(tri::Triangulation) + I = integer_type(tri) T = NTuple{3,I} E = NTuple{2,I} solid_triangles = T[] @@ -701,16 +714,16 @@ function test_iterators(tri::DT.Triangulation) ghost_edges = E[] all_edges = E[] for T in each_triangle(tri) - i, j, k = DT.triangle_vertices(T) + i, j, k = triangle_vertices(T) push!(all_triangles, (i, j, k)) - if DT.is_ghost_triangle(i, j, k) + if is_ghost_triangle(i, j, k) push!(ghost_triangles, (i, j, k)) else push!(solid_triangles, (i, j, k)) end - for (u, v) in DT.triangle_edges(i, j, k) + for (u, v) in triangle_edges(i, j, k) push!(all_edges, (u, v)) - if DT.is_ghost_edge(u, v) + if is_ghost_edge(u, v) push!(ghost_edges, (u, v)) else push!(solid_edges, (u, v)) @@ -718,7 +731,7 @@ function test_iterators(tri::DT.Triangulation) end for k in (i, j, k) push!(all_vertices, k) - if DT.is_ghost_vertex(k) + if is_ghost_vertex(k) push!(ghost_vertices, k) else push!(solid_vertices, k) @@ -732,15 +745,15 @@ function test_iterators(tri::DT.Triangulation) unique_ghost_triangle_flag = allunique(ghost_triangles) unique_ghost_edge_flag = allunique(ghost_edges) for (i, e) in enumerate(all_edges) - u, v = DT.edge_vertices(e) + u, v = edge_vertices(e) all_edges[i] = (min(u, v), max(u, v)) end for (i, e) in enumerate(solid_edges) - u, v = DT.edge_vertices(e) + u, v = edge_vertices(e) solid_edges[i] = (min(u, v), max(u, v)) end for (i, e) in enumerate(ghost_edges) - u, v = DT.edge_vertices(e) + u, v = edge_vertices(e) ghost_edges[i] = (min(u, v), max(u, v)) end unique!(all_edges) @@ -770,9 +783,9 @@ function test_iterators(tri::DT.Triangulation) sort!(all_vertices) sort!(solid_vertices) sort!(ghost_vertices) - triangle_output_flag = DT.compare_triangle_collections(all_triangles, collect(each_triangle(tri))) - solid_triangle_output_flag = DT.compare_triangle_collections(solid_triangles, collect(each_solid_triangle(tri))) - ghost_triangle_output_flag = DT.compare_triangle_collections(ghost_triangles, collect(each_ghost_triangle(tri))) + triangle_output_flag = compare_triangle_collections(all_triangles, each_triangle(tri)) + solid_triangle_output_flag = compare_triangle_collections(solid_triangles, each_solid_triangle(tri)) + ghost_triangle_output_flag = compare_triangle_collections(ghost_triangles, each_ghost_triangle(tri)) edge_output_flag = compare_edge_vectors(all_edges, each_edge(tri)) solid_edge_output_flag = compare_edge_vectors(solid_edges, each_solid_edge(tri)) ghost_edge_output_flag = compare_edge_vectors(ghost_edges, each_ghost_edge(tri)) @@ -813,14 +826,14 @@ function test_no_duplicate_segments(tri) interior_segments = get_interior_segments(tri) all_segments = get_all_segments(tri) for e in each_edge(interior_segments) - flag = !DT.contains_edge(DT.reverse_edge(e), interior_segments) - !flag && return DuplicateSegmentsState(flag, true, Int.(DT.edge_vertices(e))) + flag = !contains_edge(reverse_edge(e), interior_segments) + !flag && return DuplicateSegmentsState(flag, true, Int.(edge_vertices(e))) end for e in each_edge(all_segments) - flag = !DT.contains_edge(DT.reverse_edge(e), all_segments) - !flag && return DuplicateSegmentsState(flag, false, Int.(DT.edge_vertices(e))) + flag = !contains_edge(reverse_edge(e), all_segments) + !flag && return DuplicateSegmentsState(flag, false, Int.(edge_vertices(e))) end - return DuplicateSegmentsState(true, true, (DT.∅, DT.∅)) + return DuplicateSegmentsState(true, true, (∅, ∅)) end struct TriangulationState <: AbstractTriangulationState @@ -850,11 +863,11 @@ struct TriangulationState <: AbstractTriangulationState ghost_edge_iterator_state::IteratorState end -function TriangulationState(tri::DT.Triangulation) - has_ghosts = DT.has_ghost_triangles(tri) - DT.delete_ghost_triangles!(tri) - DT.add_ghost_triangles!(tri) - DT.clear_empty_features!(tri) +function TriangulationState(tri::Triangulation) + has_ghosts = has_ghost_triangles(tri) + delete_ghost_triangles!(tri) + add_ghost_triangles!(tri) + clear_empty_features!(tri) state = TriangulationState( Adjacent2VertexMapAdjacentMapState(tri), Adjacent2VertexMapState(tri), @@ -873,7 +886,7 @@ function TriangulationState(tri::DT.Triangulation) TriangleOrientationState(tri), test_iterators(tri)... ) - !has_ghosts && DT.delete_ghost_triangles!(tri) + !has_ghosts && delete_ghost_triangles!(tri) return state end @@ -898,7 +911,14 @@ function Base.show(io::IO, triangulation_state::TriangulationState) return end -function validate_triangulation(tri; print_result=true) +""" + validate_triangulation(tri::Triangulation; print_result=true) -> Bool + +Tests if `tri` is a valid `Triangulation`. Returns `true` if so, +and `false` otherwise. If `print_result=true` and `tri` is not a +valid triangulation, all the issues with `tri` will be printed. +""" +function validate_triangulation(tri::Triangulation; print_result=true) state = TriangulationState(tri) print_result && !test_state(state) && println(state) return test_state(state) diff --git a/test/helper_functions.jl b/test/helper_functions.jl index cb8d59565..c510ee51a 100644 --- a/test/helper_functions.jl +++ b/test/helper_functions.jl @@ -14,7 +14,7 @@ getxy((0.0, 0.0)) # avoid shadow const DT = DelaunayTriangulation -include("triangulation_validation.jl") +using DelaunayTriangulation: validate_triangulation function complicated_geometry() x1 = [collect(LinRange(0, 2, 4)), From c6dd8448fd22166fc9cf51f01c459ca6590c1e99 Mon Sep 17 00:00:00 2001 From: DanielVandH Date: Sun, 7 Jul 2024 21:41:24 +0100 Subject: [PATCH 07/14] Add info --- test/runtests.jl | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index ef316594a..3de0e4d75 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,12 +1,14 @@ # setup LocalPreferences.toml using Preferences -USE_EXACTPREDICATES = ENV["USE_EXACTPREDICATES"] +USE_EXACTPREDICATES = get(ENV, "USE_EXACTPREDICATES", "default") if USE_EXACTPREDICATES == "true" set_preferences!("DelaunayTriangulation", "USE_EXACTPREDICATES" => true) elseif USE_EXACTPREDICATES == "false" set_preferences!("DelaunayTriangulation", "USE_EXACTPREDICATES" => false) end # if USE_EXACTPREDICATES == "default", do nothing +@info "Testing with USE_EXACTPREDICATES = $USE_EXACTPREDICATES" + # get all the compilation out of the way using BenchmarkTools using CairoMakie From 3f493262f4cd880546bfb6b663b9ed7e0f3d38f8 Mon Sep 17 00:00:00 2001 From: DanielVandH Date: Sun, 7 Jul 2024 21:46:50 +0100 Subject: [PATCH 08/14] fix exports --- test/Project.toml | 1 - test/helper_functions.jl | 12 ++++++------ 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/test/Project.toml b/test/Project.toml index 496f36e71..676712c9d 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -5,7 +5,6 @@ CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" ColorSchemes = "35d6a980-a343-548e-a6ea-1d62b119f2f4" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" -DelaunayTriangulation = "927a84f5-c5f4-47a5-9785-b46e178433df" DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8" diff --git a/test/helper_functions.jl b/test/helper_functions.jl index c510ee51a..dd1611cfb 100644 --- a/test/helper_functions.jl +++ b/test/helper_functions.jl @@ -2235,12 +2235,12 @@ export get_child_from_tree export traverse_tree export poor_triangulation_example export example_triangulation -export test_adjacent2vertex_map_matches_adjacent_map -export test_state -export test_adjacent_map_matches_adjacent2vertex_map -export test_each_edge_has_two_incident_triangles -export test_triangle_orientation -export test_iterators +export DelaunayTriangulation.test_adjacent2vertex_map_matches_adjacent_map +export DelaunayTriangulation.test_state +export DelaunayTriangulation.test_adjacent_map_matches_adjacent2vertex_map +export DelaunayTriangulation.test_each_edge_has_two_incident_triangles +export DelaunayTriangulation.test_triangle_orientation +export DelaunayTriangulation.test_iterators export example_empty_triangulation export shewchuk_example_constrained export test_segment_triangle_intersections From f1de00922bbe2867c9ed445f1d6b3ec3ed402c43 Mon Sep 17 00:00:00 2001 From: DanielVandH Date: Sun, 7 Jul 2024 21:50:26 +0100 Subject: [PATCH 09/14] news --- NEWS.md | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/NEWS.md b/NEWS.md index 71957459b..8e26ee2ff 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,16 @@ # Changelog +## v.1.1.0 + +- Added the option to disable ExactPredicates.jl using Preferences.jl. See [#131](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/131). +- Added `DelauanyTriangulation.validate_triangulation` for validating triangulations. See [#131](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/131). +- Fixed a bug with the currently unused `orient(p, q, r, s)` predicate. See [#131](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/131). +- Added private functions `getz`, `_getz`, `getxyz`, and `_getxyz`. See [#131](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/pull/131). + +## v1.0.5 + +- Disabled `deepcopy` on `PolygonTree`s and made it a no-op. See [#129](https://github.com/JuliaGeometry/DelaunayTriangulation.jl/issues/129). + ## v1.0.4 Nothing breaking. Main changes: From 83f7e528943ef8600984f12af2cd66b033e501bf Mon Sep 17 00:00:00 2001 From: DanielVandH Date: Sun, 7 Jul 2024 21:57:24 +0100 Subject: [PATCH 10/14] fix syntax --- test/helper_functions.jl | 26 +++++++++++++++++--------- 1 file changed, 17 insertions(+), 9 deletions(-) diff --git a/test/helper_functions.jl b/test/helper_functions.jl index dd1611cfb..2862d37f5 100644 --- a/test/helper_functions.jl +++ b/test/helper_functions.jl @@ -2197,6 +2197,14 @@ function plot_small_angle_complexes(enricher) return fig end +using DelaunayTriangulation: + test_adjacent2vertex_map_matches_adjacent_map, + test_state, + test_adjacent_map_matches_adjacent2vertex_map, + test_each_edge_has_two_incident_triangles, + test_triangle_orientation, + test_iterators + export validate_triangulation export @_adj export _make_graph_from_adjacency @@ -2210,11 +2218,11 @@ export complicated_geometry export validate_refinement export validate_statistics export validate_tessellation -export compare_edge_vectors -export simple_geometry +export compare_edge_vectors +export simple_geometry export fixed_shewchuk_example_constrained export ⊢ -export is_sink +export is_sink export is_conformal export _validate_offcenter export example_with_special_corners @@ -2235,12 +2243,6 @@ export get_child_from_tree export traverse_tree export poor_triangulation_example export example_triangulation -export DelaunayTriangulation.test_adjacent2vertex_map_matches_adjacent_map -export DelaunayTriangulation.test_state -export DelaunayTriangulation.test_adjacent_map_matches_adjacent2vertex_map -export DelaunayTriangulation.test_each_edge_has_two_incident_triangles -export DelaunayTriangulation.test_triangle_orientation -export DelaunayTriangulation.test_iterators export example_empty_triangulation export shewchuk_example_constrained export test_segment_triangle_intersections @@ -2263,4 +2265,10 @@ export all_points_are_inside export all_diametral_circles_are_empty export compute_diametral_lens export get_points_in_diametral_lens +export test_adjacent2vertex_map_matches_adjacent_map +export test_state +export test_adjacent_map_matches_adjacent2vertex_map +export test_each_edge_has_two_incident_triangles +export test_triangle_orientation +export test_iterators end \ No newline at end of file From 8989cc92a056d2726d3be193d90fb8ac411a5b38 Mon Sep 17 00:00:00 2001 From: DanielVandH Date: Sun, 7 Jul 2024 22:23:49 +0100 Subject: [PATCH 11/14] Add NEWS, fix contributing, and fixes --- README.md | 3 ++- docs/src/extended/utils.md | 2 +- docs/src/index.md | 1 + docs/src/manual/disabling_ea.md | 20 +++++++++++++++----- 4 files changed, 19 insertions(+), 7 deletions(-) diff --git a/README.md b/README.md index a51553662..553338b4d 100644 --- a/README.md +++ b/README.md @@ -113,8 +113,9 @@ When contributing in the form of a pull request, there are a few important featu 2. **Document any new functions**. If your contribution involves any new functions, make sure they are well documented even if they are only for internal use. 3. **You need to include appropriate tests for your contribution**. For small changes, simple tests are fine, but for larger changes that implement or change an algorithm, understanding what tests you need to write is crucial when you are implementing a new feature. For example, if you wanted to implement a new algorithm for inserting a curve into a triangulation, you need to make sure that you test things like (1) point sets in general position, (2) point sets with many cocircular points, (3) collinear edges, and so on. In particular, degeneracies are some of the most important things to test for. This is not simple to do if you are inexperienced, so do feel free to ask for guidance. 4. **Your contribution should only do one thing**. To make sure that it is easy to track your changes, and to make it easier to track regressions across versions, please try to only make your contribution do one thing. For example, do not move around a bunch of files while at the same time implementing a function somewhere else - make two pull requests. +5. **Add to NEWS.md**. Describe your change in the NEWS.md file. -If you are just contributing something minimal, for example a typo, even a blank PR is fine so long as it is obvious what you have done. +If you are just contributing something minimal, for example a typo, even a blank PR description is fine so long as it is obvious what you have done. ## Similar Packages diff --git a/docs/src/extended/utils.md b/docs/src/extended/utils.md index 99505482b..45ee8e486 100644 --- a/docs/src/extended/utils.md +++ b/docs/src/extended/utils.md @@ -33,5 +33,5 @@ DefaultAdjacentValue GhostVertex ε ∅ -USE_EXACTPREDICATE +USE_EXACTPREDICATES ``` \ No newline at end of file diff --git a/docs/src/index.md b/docs/src/index.md index 6ca4a2732..4f621609d 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -71,6 +71,7 @@ When contributing in the form of a pull request, there are a few important featu 2. **Document any new functions**. If your contribution involves any new functions, make sure they are well documented even if they are only for internal use. 3. **You need to include appropriate tests for your contribution**. For small changes, simple tests are fine, but for larger changes that implement or change an algorithm, understanding what tests you need to write is crucial when you are implementing a new feature. For example, if you wanted to implement a new algorithm for inserting a curve into a triangulation, you need to make sure that you test things like (1) point sets in general position, (2) point sets with many cocircular points, (3) collinear edges, and so on. In particular, degeneracies are some of the most important things to test for. This is not simple to do if you are inexperienced, so do feel free to ask for guidance. 4. **Your contribution should only do one thing**. To make sure that it is easy to track your changes, and to make it easier to track regressions across versions, please try to only make your contribution do one thing. For example, do not move around a bunch of files while at the same time implementing a function somewhere else - make two pull requests. +5. **Add to NEWS.md**. Describe your change in the NEWS.md file. If you are just contributing something minimal, for example a typo, even a blank PR is fine so long as it is obvious what you have done. diff --git a/docs/src/manual/disabling_ea.md b/docs/src/manual/disabling_ea.md index 946ef54ef..3098b4d7c 100644 --- a/docs/src/manual/disabling_ea.md +++ b/docs/src/manual/disabling_ea.md @@ -4,7 +4,7 @@ CurrentModule = DelaunayTriangulation # Disabling Exact Predicates -For performance reasons, you may find it useful to want to disable exact predicates using [ExactPredicates.jl](https://github.com/lairez/ExactPredicates.jl). This can be easily done using a setup with [Preferences.jl](https://github.com/JuliaPackaging/Preferences.jl), but before you consider disabling exact predicates.jl, there are a few things to be aware of. If you just want to disable them without reading a lot of information warning you about the consequences, pleas skip to the end. +For performance reasons, you may find it useful to want to disable exact predicates using [ExactPredicates.jl](https://github.com/lairez/ExactPredicates.jl). This can be easily done using a setup with [Preferences.jl](https://github.com/JuliaPackaging/Preferences.jl), but before you consider disabling exact predicates.jl, there are a few things to be aware of. If you just want to disable them without reading a lot of information warning you about the consequences, please skip to the end. ## Why use exact predicates? @@ -73,18 +73,28 @@ into it. If you later want to re-enable exact predicates, either delete the file When you are not using exact predicates, you may want to check if your computed triangulation is actually a valid Delaunay triangulation. We provide the function `DelaunayTriangulation.validate_triangulation` for this purpose. This functionality is quite slow to use and is not currently optimised or well-documented (contributions towards addressing these issues are welcome), but it will work. One important note is that this check does actually use predicates in certain areas, so this check is still not guaranteed to be 100% accurate. Here is an example of its use. -```@example validate +```julia using DelaunayTriangulation tri = triangulate(rand(2, 50)) DelaunayTriangulation.validate_triangulation(tri) ``` - -```@example validate +```julia +true +``` +```julia T = first(each_solid_triangle(tri)) DelaunayTriangulation.delete_triangle!(tri, T) # break the triangulation for this example DelaunayTriangulation.validate_triangulation(tri) ``` +```julia +The edge (12, 40) does not have two incident triangles. +The edge (12, 40) appears as an edge in the graph but it and its reverse are not both a key of the adjacent map. -```@example validate +false +``` +```julia DelaunayTriangulation.validate_triangulation(tri; print_result = false) ``` +```julia +false +``` From 4fdb0f56731399235a890d81957c485deb823f4a Mon Sep 17 00:00:00 2001 From: DanielVandH Date: Mon, 8 Jul 2024 00:10:28 +0100 Subject: [PATCH 12/14] Fix some tests --- test/data_structures/curves.jl | 59 +++++++++++----------- test/point_location/jump_and_march.jl | 70 ++++++++++++++------------- test/voronoi/voronoi.jl | 70 ++++++++++++++------------- 3 files changed, 104 insertions(+), 95 deletions(-) diff --git a/test/data_structures/curves.jl b/test/data_structures/curves.jl index f2ee2ba2e..2f04225a3 100644 --- a/test/data_structures/curves.jl +++ b/test/data_structures/curves.jl @@ -5,6 +5,7 @@ using LinearAlgebra EllipticalArc = DT.EllipticalArc # shadow using CairoMakie using StableRNGs +using Preferences using StructEquality using ForwardDiff using ReferenceTests @@ -1978,36 +1979,38 @@ end end ## Closest point - for spl in (spl, pspl) - for _ in 1:500 - p = randn(2) * 30 |> Tuple - t′, q = DT.get_closest_point(spl, p) - @test q ⪧ spl(t′) - _t, _q = closest_point_on_curve(spl, p) - if !(spl == pspl) - @test _t ≈ t′ rtol = 1e-1 atol = 1e-1 - elseif t ≠ 0 && t ≠ 1 && t′ ≠ 0 && t′ ≠ 1 - @test _t ≈ t′ rtol = 1e-1 atol = 1e-1 + if load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true) + for spl in (spl, pspl) + for _ in 1:500 + p = randn(2) * 30 |> Tuple + t′, q = DT.get_closest_point(spl, p) + @test q ⪧ spl(t′) + _t, _q = closest_point_on_curve(spl, p) + if !(spl == pspl) + @test _t ≈ t′ rtol = 1e-1 atol = 1e-1 + elseif t ≠ 0 && t ≠ 1 && t′ ≠ 0 && t′ ≠ 1 + @test _t ≈ t′ rtol = 1e-1 atol = 1e-1 + end + @test q ⪧ _q rtol = 1e-1 atol = 1e-1 + @test DT.dist(p, _q) ≈ DT.dist(p, q) rtol = 1e-1 atol = 1e-1 end - @test q ⪧ _q rtol = 1e-1 atol = 1e-1 - @test DT.dist(p, _q) ≈ DT.dist(p, q) rtol = 1e-1 atol = 1e-1 - end - for t in LinRange(0, 1, 150) - p = spl(t) - t′, q = DT.get_closest_point(spl, p) - @test q ⪧ spl(t′) - if !(spl == pspl) - @test t ≈ t′ rtol = 1e-1 atol = 1e-1 - elseif t ≠ 0 && t ≠ 1 - @test t ≈ t′ rtol = 1e-1 atol = 1e-1 + for t in LinRange(0, 1, 150) + p = spl(t) + t′, q = DT.get_closest_point(spl, p) + @test q ⪧ spl(t′) + if !(spl == pspl) + @test t ≈ t′ rtol = 1e-1 atol = 1e-1 + elseif t ≠ 0 && t ≠ 1 + @test t ≈ t′ rtol = 1e-1 atol = 1e-1 + end + @test p ⪧ q rtol = 1e-1 atol = 1e-1 + end + @test DT.get_closest_point(spl, spl.control_points[1]) == (0.0, spl.control_points[1]) + if spl ≠ pspl + @test DT.get_closest_point(spl, spl.control_points[end]) == (1.0, spl.control_points[end]) + else + @test DT.get_closest_point(spl, spl.control_points[end]) == (0.0, spl.control_points[end]) end - @test p ⪧ q rtol = 1e-1 atol = 1e-1 - end - @test DT.get_closest_point(spl, spl.control_points[1]) == (0.0, spl.control_points[1]) - if spl ≠ pspl - @test DT.get_closest_point(spl, spl.control_points[end]) == (1.0, spl.control_points[end]) - else - @test DT.get_closest_point(spl, spl.control_points[end]) == (0.0, spl.control_points[end]) end end diff --git a/test/point_location/jump_and_march.jl b/test/point_location/jump_and_march.jl index f9b7a90cd..dae6fbb77 100644 --- a/test/point_location/jump_and_march.jl +++ b/test/point_location/jump_and_march.jl @@ -4,7 +4,7 @@ using LinearAlgebra using Random using StableRNGs using StatsBase - +using Preferences tri, label_map, index_map = simple_geometry() @@ -121,40 +121,42 @@ rep[3].y = mean([12.0, 6.0, 2.0, 4.0, 6.0, 10.0]) rep[3].y = mean([12.0, 6.0, 2.0, 4.0, 6.0, 10.0]) end - @testset "Test that we can find a point in every triangle" begin - for run in 1:6 - for V in each_triangle(tri.triangles) - rand() < 0.5 && continue # skip 50% - if !DT.is_exterior_ghost_triangle(tri, triangle_vertices(V)...) - i, j, k = triangle_vertices(V) - p, q, r = get_point(tri, i, j, k) - local c - c = (p .+ q .+ r) ./ 3 - for k in DT.each_solid_vertex(tri) - _V = jump_and_march(tri, c; k, concavity_protection=true) - @test DT.is_positively_oriented(DT.triangle_orientation(tri, _V)) - if !DT.is_ghost_triangle(_V...) - @test DT.compare_triangles(_V, V) && - DT.is_inside(DT.point_position_relative_to_triangle(tri, - _V, - c)) - else - local V1, V2 - V1 = DT.sort_triangle(V) - V2 = DT.sort_triangle(_V) - i1 = DT.geti(V1) - i2 = DT.geti(V2) - if i1 ≠ i2 - i1 = i1 - 1 - end - if i1 ≠ i2 - i1 = i1 + 1 + if load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true) + @testset "Test that we can find a point in every triangle" begin + for run in 1:6 + for V in each_triangle(tri.triangles) + rand() < 0.5 && continue # skip 50% + if !DT.is_exterior_ghost_triangle(tri, triangle_vertices(V)...) + i, j, k = triangle_vertices(V) + p, q, r = get_point(tri, i, j, k) + local c + c = (p .+ q .+ r) ./ 3 + for k in DT.each_solid_vertex(tri) + _V = jump_and_march(tri, c; k, concavity_protection=true) + @test DT.is_positively_oriented(DT.triangle_orientation(tri, _V)) + if !DT.is_ghost_triangle(_V...) + @test DT.compare_triangles(_V, V) && + DT.is_inside(DT.point_position_relative_to_triangle(tri, + _V, + c)) + else + local V1, V2 + V1 = DT.sort_triangle(V) + V2 = DT.sort_triangle(_V) + i1 = DT.geti(V1) + i2 = DT.geti(V2) + if i1 ≠ i2 + i1 = i1 - 1 + end + if i1 ≠ i2 + i1 = i1 + 1 + end + _V = DT.construct_triangle(typeof(V), i1, DT.getj(V1), DT.getk(V1)) + @test DT.compare_triangles(_V, V) && + DT.is_inside(DT.point_position_relative_to_triangle(tri, + _V, + c)) end - _V = DT.construct_triangle(typeof(V), i1, DT.getj(V1), DT.getk(V1)) - @test DT.compare_triangles(_V, V) && - DT.is_inside(DT.point_position_relative_to_triangle(tri, - _V, - c)) end end end diff --git a/test/voronoi/voronoi.jl b/test/voronoi/voronoi.jl index d957bae9f..1639af91c 100644 --- a/test/voronoi/voronoi.jl +++ b/test/voronoi/voronoi.jl @@ -221,31 +221,33 @@ end @test findmin(all_dists)[2] == u end - points = [ - (0.0, 0.0), (-1.0, 1.0), (-0.5, 1.0), (0.0, 1.0), (0.5, 1.0), (1.0, 1.0), - (1.0, 0.8), (1.0, 0.0), (1.0, -0.5), (1.0, -1.0), - (0.1, -1.0), (-0.8, -1.0), (-1.0, -1.0), - (-1.0, -0.7), (-1.0, -0.1), (-1.0, 0.6), - (-0.1, -0.8), (0.2, -0.8), - (-0.6, -0.4), (0.9, 0.0), (-0.5, 0.5), (-0.4, 0.6), (-0.1, 0.8) - ] - tri = triangulate(points, delete_ghosts=false) - vorn = voronoi(tri) - @test validate_tessellation(vorn) - xg = LinRange(-1, 1, 250) - yg = LinRange(-1, 1, 250) - x = vec([x for x in xg, _ in yg]) - y = vec([y for _ in xg, y in yg]) - for (ξ, η) in zip(x, y) - p = (ξ, η) - u = get_nearest_neighbour(vorn, p) - @inferred get_nearest_neighbour(vorn, p) - all_dists = [norm(p .- get_generator(vorn, i)) for i in sort(collect(each_generator(vorn)))] - k = findmin(all_dists)[2] - @test k == u - for m in DT.each_point_index(tri) - u = get_nearest_neighbour(vorn, p, try_points=m) - @test u == k + if load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true) + points = [ + (0.0, 0.0), (-1.0, 1.0), (-0.5, 1.0), (0.0, 1.0), (0.5, 1.0), (1.0, 1.0), + (1.0, 0.8), (1.0, 0.0), (1.0, -0.5), (1.0, -1.0), + (0.1, -1.0), (-0.8, -1.0), (-1.0, -1.0), + (-1.0, -0.7), (-1.0, -0.1), (-1.0, 0.6), + (-0.1, -0.8), (0.2, -0.8), + (-0.6, -0.4), (0.9, 0.0), (-0.5, 0.5), (-0.4, 0.6), (-0.1, 0.8) + ] + tri = triangulate(points, delete_ghosts=false) + vorn = voronoi(tri) + @test validate_tessellation(vorn) + xg = LinRange(-1, 1, 250) + yg = LinRange(-1, 1, 250) + x = vec([x for x in xg, _ in yg]) + y = vec([y for _ in xg, y in yg]) + for (ξ, η) in zip(x, y) + p = (ξ, η) + u = get_nearest_neighbour(vorn, p) + @inferred get_nearest_neighbour(vorn, p) + all_dists = [norm(p .- get_generator(vorn, i)) for i in sort(collect(each_generator(vorn)))] + k = findmin(all_dists)[2] + @test k == u + for m in DT.each_point_index(tri) + u = get_nearest_neighbour(vorn, p, try_points=m) + @test u == k + end end end end @@ -894,14 +896,16 @@ end @test flag / tot > 0.9 end -@testset "Lattice" begin - for _ in 1:100 - tri = triangulate_rectangle(0, 1, 0, 1, 11, 11, delete_ghosts=false) - tri = triangulate(tri.points) - vorn = voronoi(tri) - @test validate_tessellation(vorn; check_convex=false) - vorn = voronoi(tri, clip=true) - @test validate_tessellation(vorn; check_convex=false, check_adjacent=false) +if load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true) + @testset "Lattice" begin + for _ in 1:100 + tri = triangulate_rectangle(0, 1, 0, 1, 11, 11, delete_ghosts=false) + tri = triangulate(tri.points) + vorn = voronoi(tri) + @test validate_tessellation(vorn; check_convex=false) + vorn = voronoi(tri, clip=true) + @test validate_tessellation(vorn; check_convex=false, check_adjacent=false) + end end end From 1bc0fd36cf4739d5597a8a98ce07bf635fe65805 Mon Sep 17 00:00:00 2001 From: DanielVandH Date: Mon, 8 Jul 2024 01:25:12 +0100 Subject: [PATCH 13/14] Fix --- .../segment_location.jl | 47 ++++++++++--------- test/point_location/jump_and_march.jl | 18 +++---- 2 files changed, 35 insertions(+), 30 deletions(-) diff --git a/test/constrained_triangulation/segment_location.jl b/test/constrained_triangulation/segment_location.jl index a6f2aa5fa..c46658514 100644 --- a/test/constrained_triangulation/segment_location.jl +++ b/test/constrained_triangulation/segment_location.jl @@ -1,6 +1,7 @@ using ..DelaunayTriangulation const DT = DelaunayTriangulation using CairoMakie +using Preferences @testset "Shewchuk Example: A small example with some collinearities" begin tri = shewchuk_example_constrained() @@ -278,29 +279,31 @@ end end end -@testset "A previously broken example" begin - for m in 1:100 - a = -0.1 - b = 0.1 - c = -0.01 - d = 0.01 - nx = 25 - ny = 25 - tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=true) - tri = triangulate(get_points(tri)) - for i in 2:24 - add_segment!(tri, i, 600 + i) +if load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true) + @testset "A previously broken example" begin + for m in 1:100 + a = -0.1 + b = 0.1 + c = -0.01 + d = 0.01 + nx = 25 + ny = 25 + tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=true) + tri = triangulate(get_points(tri)) + for i in 2:24 + add_segment!(tri, i, 600 + i) + end + tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=true) + tri = triangulate(get_points(tri)) + e = (23, 71) + history = DT.PointLocationHistory{NTuple{3,Int},NTuple{2,Int},Int}() + jump_and_march(tri, get_point(tri, 71); + m=nothing, k=23, store_history=true, history=history) + collinear_segments = history.collinear_segments + DT.connect_segments!(collinear_segments) + DT.extend_segments!(collinear_segments, e) + @test collinear_segments == [(23, 47), (47, 71)] end - tri = triangulate_rectangle(a, b, c, d, nx, ny; delete_ghosts=false, single_boundary=true) - tri = triangulate(get_points(tri)) - e = (23, 71) - history = DT.PointLocationHistory{NTuple{3,Int},NTuple{2,Int},Int}() - jump_and_march(tri, get_point(tri, 71); - m=nothing, k=23, store_history=true, history=history) - collinear_segments = history.collinear_segments - DT.connect_segments!(collinear_segments) - DT.extend_segments!(collinear_segments, e) - @test collinear_segments == [(23, 47), (47, 71)] end end diff --git a/test/point_location/jump_and_march.jl b/test/point_location/jump_and_march.jl index dae6fbb77..d055e5f0a 100644 --- a/test/point_location/jump_and_march.jl +++ b/test/point_location/jump_and_march.jl @@ -164,14 +164,16 @@ rep[3].y = mean([12.0, 6.0, 2.0, 4.0, 6.0, 10.0]) end end - @testset "Test that we don't break for points already in the triangulation" begin - for _ in 1:6 - for k in DT.each_solid_vertex(tri) - rand() < 1 / 2 && continue - for j in DT.each_solid_vertex(tri) - _V = jump_and_march(tri, get_point(tri, k); k=j) - @test k ∈ triangle_vertices(_V) - @test DT.is_positively_oriented(DT.triangle_orientation(tri, _V)) + if !load_preference(DelaunayTriangulation, "USE_EXACTPREDICATES", true) && tri_idx ≠ 3 + @testset "Test that we don't break for points already in the triangulation" begin + for _ in 1:6 + for k in DT.each_solid_vertex(tri) + rand() < 1 / 2 && continue + for j in DT.each_solid_vertex(tri) + _V = jump_and_march(tri, get_point(tri, k); k=j) + @test k ∈ triangle_vertices(_V) + @test DT.is_positively_oriented(DT.triangle_orientation(tri, _V)) + end end end end From cea95caa22b107ec77eb7fdd8ae4f89b5cc43288 Mon Sep 17 00:00:00 2001 From: DanielVandH Date: Mon, 8 Jul 2024 08:24:46 +0100 Subject: [PATCH 14/14] Fix occasionally broken test --- .github/workflows/CI.yml | 2 +- test/runtests.jl | 2 +- test/utils.jl | 24 +++++++++++++----------- 3 files changed, 15 insertions(+), 13 deletions(-) diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml index f4752c309..e8550798d 100644 --- a/.github/workflows/CI.yml +++ b/.github/workflows/CI.yml @@ -12,7 +12,7 @@ concurrency: cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }} jobs: test: - name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} - ExactPredicates ${{ matrix.preferences }} + name: Julia ${{ matrix.version }} - ${{ matrix.arch }} - ${{ github.event_name }} - ExactPredicates ${{ matrix.preferences }} - ${{ matrix.os }} runs-on: ${{ matrix.os }} env: JULIA_NUM_THREADS: ${{ matrix.julia-threads }} diff --git a/test/runtests.jl b/test/runtests.jl index 3de0e4d75..5e1c98102 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -52,7 +52,7 @@ function safe_include(filename; name=filename, push=true, verbose = true) # Work push && push!(ALL_TEST_SCRIPTS, normpath(filename)) mod = @eval module $(gensym()) end @info "[$(ct())] Testing $name" - @testset verbose = verbose "$name" begin + @testset verbose = verbose "$(basename(name))" begin @eval mod using ..HelperFunctions @eval mod using ..Test Base.include(mod, filename) diff --git a/test/utils.jl b/test/utils.jl index b57e1cbd9..f837c929b 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -1224,17 +1224,19 @@ end end @testset "replace_(boundary/ghost)_triangle_with_(ghost/boundary)_triangle" begin - tri = triangulate(rand(2, 5000), delete_ghosts=false) - _tri = triangulate_rectangle(0.0, 1.0, 0.0, 1.0, 25, 25, delete_ghosts=false) - for tri in (tri, _tri) - for T in each_ghost_triangle(tri) - T = DT.sort_triangle(T) - i, j, k = triangle_vertices(T) - w = get_adjacent(tri, j, i) - V = (j, i, w) - @test DT.replace_boundary_triangle_with_ghost_triangle(tri, V) == T - @test DT.replace_ghost_triangle_with_boundary_triangle(tri, T) == V - if (i, j) ∉ ((624, 625), (1, 26), (2, 1), (625, 600)) # corners of the lattice have two associated ghosts + for _ in 1:10 + tri = triangulate(rand(2, 5000), delete_ghosts=false) + _tri = triangulate_rectangle(0.0, 1.0, 0.0, 1.0, 25, 25, delete_ghosts=false) + for tri in (tri, _tri) + for T in each_ghost_triangle(tri) + T = DT.sort_triangle(T) + i, j, k = triangle_vertices(T) + w = get_adjacent(tri, j, i) + V = (j, i, w) + kr = count(e -> DT.is_boundary_edge(tri, e...), DT.triangle_edges(V)) + kr == 2 && continue # two associated ghosts, test is not clear + @test DT.replace_boundary_triangle_with_ghost_triangle(tri, V) == T + @test DT.replace_ghost_triangle_with_boundary_triangle(tri, T) == V V = (i, w, j) @test DT.replace_boundary_triangle_with_ghost_triangle(tri, V) == T @test DT.replace_ghost_triangle_with_boundary_triangle(tri, T) == (j, i, w)