# [ANN] DelaunayTriangulation v1.0: Curved domains and improved docs/code

Hi all,

I have recently released the first major version of my package DelaunayTriangulation.jl, a package for computing Delaunay triangulations and Voronoi tessellations in two dimensions. The most important new feature here is the ability to now triangulate and refine domains defined by curves rather than only piecewise linear boundaries. Thus, the packageβs major features are now:

• Unconstrained and constrained triangulations.
• Mesh refinement, with support for generic constraints.
• Voronoi tessellations.
• Voronoi tessellations clipped to a rectangle or to the convex hull.
• Centroidal Voronoi tessellations.
• Triangulations of curve-bounded domains.

More features are listed in the README and in the documentation. In addition to this new feature of curve-bounded domains, the documentation has been completely rewritten and should hopefully be a lot easier to navigate. I give some better examples of using the code, some example applications, and give better details of the mathematics involved.

The curves I provide support for are:

• `BSpline`
• `BezierCurve`
• `CircularArc`
• `CatmullRomSpline`
• `LineSegment`
• `EllipticalArc`

But you can also define your own; see the docstring for `DelaunayTriangulation.AbstractParametricCurve` or the docs.

Here is an example of triangulating and refining a domain defined by an annulus.

``````using DelaunayTriangulation
using CairoMakie
Rβ = 1.0
Rβ = 2.0
outer_circle = CircularArc((Rβ, 0.0), (Rβ, 0.0), (0.0, 0.0))
inner_circle = CircularArc((Rβ, 0.0), (Rβ, 0.0), (0.0, 0.0), positive=false)
points = NTuple{2,Float64}[]
tri = triangulate(points; boundary_nodes=[[[outer_circle]], [[inner_circle]]])
A = 2Ο * (Rβ^2 - Rβ^2)
refine!(tri; max_area=2e-3A, min_angle=33.0)
fig, ax, sc = triplot(tri)
fig
``````

If you want some more examples other than just those in the docs, I will soon be updating FiniteVolumeMethod.jl to use these new curve-bounded domains rather than e.g. manually discretising a circle and triangulating that

Many other improvements have been made. For example, you no longer need to do `check_args=false` for more complicated domains, as I have now implemented an automatic way for the package to check that your domain is valid even with nested domains and disjoint domains. The public API has now also been fully defined for the first time. See the full changelog here.

31 Likes

Cool to see this hit 1.0! Congrats! How hard would it be to extend this to 3D?

2 Likes

Thanks for this great contribution. The documentation is very well written.

1 Like

3D seems to be a much bigger challenge. I did at some point spend a few weeks trying to see if it could work easily. I got as far as making an unconstrained triangulator working, but the corner cases in 3D are so much more annoying that I gave up (2D is bad enough - just look at how long the point location code isβ¦). Not sure where that 3D code is now unfortunately. Another challenge is just in designing the data structures for working efficiently with 3D geometries, especially for constrained geometries. None of my work requires 3D triangulations unfortunately, else I probably wouldβve kept trying.

The book I used for all this (here) has a lot of great discussion on 3D, and more general cases like triangulating on a surface. If someone wanted to take up the mantle and try and get a 3D code working, this book would be where Iβd start. For now, I think TetGen.jl is probably the best to use in Julia? Never used it.

3 Likes

I just want to say I love the level of detail in the docstrings

5 Likes

Thank you Iβve tried to make most functions have detailed (enough) docstrings to make everything easier to understand, especially for internal functions incase anyone else wants to contribute to the package - rather than me being the only one with any idea whatβs going on

1 Like

I love thorough docstrings, and IMO tricky internal functions are often most direly in need of them.

[Shameless self-promotion] Iβm actually working on a package to help check that all functions/structs/global variables have decent attached docstrings: GitHub - tecosaur/CheckDoc.jl: Documentation linting, it still needs more work before Iβm ready to (more) publically announce it though.

6 Likes

Is there a GitHub shortcut to star every repo a user have ever written or ever will? That would save me some work with people like @Lilith and you

9 Likes

FiniteVolumeMethod.jl has now been updated to now use curve-bounded domains where applicable. As an example of how you could use these new features to specify a mesh, here I specify an annulus where the outer circle has been split into two parts to allow for different boundary conditions on each part. Extra segments are added to define a perturbed interface so that the diffusivity is different on each side. This example comes from this example.

``````using DelaunayTriangulation, CairoMakie
# Some diffusion_parameters
Ξ΅ = 0.05 # perturbation parameter for the perturbed interface
g = ΞΈ -> sin(3ΞΈ) + cos(5ΞΈ) # perturbation function
R1_f = let Rβ = Rβ, Ξ΅ = Ξ΅, g = g # use let for type stability
ΞΈ -> Rβ * (1.0 + Ξ΅ * g(ΞΈ))
Ο΅r = 0.25 # interior hole radius
# Define the boundary:
#   1. The boundary should be absorbing for most of it,
#   2. but have a small reflecting arc on the outer circle,
#   3. and a small absorbing hole in the interior.
dirichlet = CircularArc((Rβ * cos(Ο΅r), Rβ * sin(Ο΅r)), (Rβ * cos(2Ο - Ο΅r), Rβ * sin(2Ο - Ο΅r)), (0.0, 0.0))
neumann = CircularArc((Rβ * cos(2Ο - Ο΅r), Rβ * sin(2Ο - Ο΅r)), (Rβ * cos(Ο΅r), Rβ * sin(Ο΅r)), (0.0, 0.0))
hole = CircularArc((0.0, 1.0), (0.0, 1.0), (0.0, 0.0), positive=false)
boundary_nodes = [[[dirichlet], [neumann]], [[hole]]]
points = [(-2.0, 0.0), (0.0, 2.95)] # also add small point holes that serve as obstacles
tri = triangulate(points; boundary_nodes)
# Now add the perturbed interface. There is no capability for adding
# curves that are not the boundary, so we need to do this manually.
ΞΈ = LinRange(0, 2Ο, 250)
xin = @views (@. R1_f(ΞΈ) * cos(ΞΈ))[begin:end-1]
yin = @views (@. R1_f(ΞΈ) * sin(ΞΈ))[begin:end-1]
for i in 2:length(xin)
n = DelaunayTriangulation.num_points(tri)
end
n = DelaunayTriangulation.num_points(tri)
tri1 = deepcopy(tri) # for plotting after
refine!(tri; max_area=1e-3get_area(tri));
# Plot
fig = Figure()
ax = Axis(fig[1, 1], width=400, height=400, title="Pre-refinement")
triplot!(ax, tri1)
ax = Axis(fig[1, 2], width=400, height=400, title="Post-refinement")
triplot!(ax, tri)
resize_to_layout!(fig)
fig
``````

5 Likes

Now that AdaptivePredicates.jl has been released, Iβve made it the default kernel for computing predicates in DelaunayTriangulation as of v1.1.0. This gives a significant improvement in performance (sometimes). We can now also give robust computations of triangle areas which has fixed some bugs. Moreover, since AdaptivePredicates.jl supports `Float32`, DelaunayTriangulation.jl now supports `Float32` without having to convert to `Float64` everywhere (unless the `ExactKernel()` is used).

Additionally, instead of forcing users to use ExactPredicates.jl or AdaptivePredicates.jl, `triangulate` (and `voronoi`, `refine!`, etc.) has a `predicates` keyword argument that you can use for specifying the kernels `FastKernel()`, `AdaptiveKernel()`, or `ExactKernel()`. The difference between these choices is described here.

Here is a quick example comparing the kernels applied to the triangulation of points on a lattice (meaning there are lots of collinear and cocircular points, pretty much the worst case scenario). Here, `FastKernel()` is not even likely to give a valid triangulation but I include it anyway.

``````julia> points = DelaunayTriangulation.get_lattice_points(0, 1, 0, 1, 25, 25, LinearIndices((1:25, 1:25))); # not public, fyi

julia> @benchmark triangulate(\$points; predicates = \$AdaptiveKernel(), randomise = \$false) samples = 1000
BenchmarkTools.Trial: 431 samples with 1 evaluation.
Range (min β¦ max):   9.882 ms β¦ 158.046 ms  β GC (min β¦ max): 0.00% β¦ 92.85%
Time  (median):     10.969 ms               β GC (median):    0.00%
Time  (mean Β± Ο):   11.596 ms Β±   7.199 ms  β GC (mean Β± Ο):  3.74% Β±  6.06%

β  ββββββββββ
ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ β
9.88 ms         Histogram: frequency by time         18.3 ms <

Memory estimate: 2.92 MiB, allocs estimate: 8295.

julia> @benchmark triangulate(\$points; predicates = \$ExactKernel(), randomise = \$false) samples = 1000
BenchmarkTools.Trial: 238 samples with 1 evaluation.
Range (min β¦ max):  14.463 ms β¦ 209.383 ms  β GC (min β¦ max):  0.00% β¦ 56.42%
Time  (median):     16.461 ms               β GC (median):     0.00%
Time  (mean Β± Ο):   21.007 ms Β±  17.397 ms  β GC (mean Β± Ο):  11.87% Β± 12.37%

βββββ
ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ β
14.5 ms       Histogram: log(frequency) by time      72.1 ms <

Memory estimate: 8.18 MiB, allocs estimate: 187394.

julia> @benchmark triangulate(\$points; predicates = \$FastKernel(), randomise = \$false) samples = 1000
BenchmarkTools.Trial: 332 samples with 1 evaluation.
Range (min β¦ max):   9.369 ms β¦ 207.123 ms  β GC (min β¦ max): 0.00% β¦ 88.74%
Time  (median):     15.539 ms               β GC (median):    0.00%
Time  (mean Β± Ο):   15.080 ms Β±  11.215 ms  β GC (mean Β± Ο):  4.28% Β±  5.91%

βββββββ                βββββ    ββ
ββββββββββββββββββββββββββββββββββββββββββββββββββββββββββββ β
9.37 ms         Histogram: frequency by time         24.8 ms <

Memory estimate: 2.83 MiB, allocs estimate: 8306.
``````

There are some other benchmarks Iβve done but, typically, `AdaptiveKernel()` being comparable to `FastKernel()` is expected and `ExactKernel()` lags a bit behind with many more allocations. When there are no issues like collinear points or cocircular points, all the kernels are about the same.

(Not all predicates when using `AdaptiveKernel()` use AdaptivePredicates.jl yet. `angle_is_acute` and `parallelorder` still use ExactPredicates.jl. See Other predicates Β· Issue #18 Β· JuliaGeometry/AdaptivePredicates.jl Β· GitHub)

7 Likes