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

Package Announcements
,
1

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

qTuczvg

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.

36 Likes
2

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

2 Likes
3

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

1 Like
4

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
5

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

5 Likes
6

Thank you :slight_smile: 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 :sweat_smile:

1 Like
7

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.

7 Likes
8

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
9

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
R₁, Rβ‚‚ = 2.0, 3.0 # perturbed interface radius, outer radius
Ξ΅ = 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(ΞΈ))
end # interface radius function
Ο΅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]
add_point!(tri, xin[1], yin[1])
for i in 2:length(xin)
    add_point!(tri, xin[i], yin[i])
    n = DelaunayTriangulation.num_points(tri)
    add_segment!(tri, n - 1, n)
end
n = DelaunayTriangulation.num_points(tri)
add_segment!(tri, n - 1, n)
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

KzZRS5N
KzZRS5N728Γ—382 89 KB

6 Likes
10

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
11

I’ve now released DelaunayTriangulation.jl 1.3.0 which is a pretty big release, bringing with it three really nice new features:

  1. Weighted Delaunay triangulations. These are triangulations of point sets where each point has an associated weight.
  2. Power diagrams. These are dual to the weighted Delaunay triangulations. Here, the Voronoi cells are defined in terms of the power distance metric instead of the Euclidean distance, i.e. the metric \pi(p, q) = d(p, q)^2 - w_p - w_q, where d(p, q) is the Euclidean distance and w_p and w_q are the weights of p and q, respectively. This could be a nice way to add non-Euclidean metric support to NaturalNeighbours.jl if anyone was interested in that (note that calls to triangle_circumcenter would be replaced with appropriate calls to triangle_orthocenter when is_weighted(tri), in addition to whatever else is needed.)
  3. Clipped Voronoi tessellations to convex polygons other than just the convex hull or rectangles (still no support for non-convex boundaries). The limitation to convex polygons is because the Sutherland-Hodgman algorithm is used. Hopefully, in the future, I can eventually use the much better VoroCrust algorithm.

I show examples of these three features below. Please see the docs for more examples Introduction Β· DelaunayTriangulation.jl.

image
image1003Γ—349 75.8 KB

using DelaunayTriangulation, CairoMakie, StableRNGs
using DelaunayTriangulation: EllipticalArc 
# Weighted triangulation 
rng = StableRNG(123)
points = rand(rng, 2, 50)
weights = randn(rng, 50)
tri1 = triangulate(points; weights, rng)

# Power diagram: Just pass a weighted triangulation to voronoi! 
vorn2 = voronoi(tri1; rng)

# Clipping to an ellipse 
points = 5randn(rng, 2, 100)
weights = rand(rng, 100)
tri3 = triangulate(points; rng, weights)
ellipse = EllipticalArc((2.0, 0.0), (2.0, 0.0), (0.0, 0.0), 2.0, 4.0, 0.0)
t = LinRange(0, 1, 200)
clip_points = ellipse.(t)
clip_vertices = [1:(length(clip_points)-1); 1]
clip_polygon = (clip_points, clip_vertices)
vorn3 = voronoi(tri3; rng, clip = true, clip_polygon)

fig = Figure()
ax1 = Axis(fig[1, 1], title = "Weighted", width = 300, height = 300)
ax2 = Axis(fig[1, 2], title = "Power", width = 300, height = 300)
ax3 = Axis(fig[1, 3], title = "Clipped", width = 300, height = 300)
triplot!(ax1, tri1)
voronoiplot!(ax2, vorn2, show_generators = false, clip = (-5, 5, -5, 5))
voronoiplot!(ax3, vorn3, show_generators = false)
resize_to_layout!(fig)
fig
11 Likes
12

Out of shear curiosity: can of above be used to generate meshes with boundary layer refinement to capture the use of wall function in Reynolds-averaged turbulent flow? Thx!

13

I would probably need an example of what you need to be sure, but the refine! function allows you to pass a custom refinement function that could in principle allow you to refine more carefully around boundaries.

I have two examples of custom refinement:

Do those help? There are also tutorials that show you how to iterate properly over the boundaries if needed.

1 Like
14

Thanks for this info!

1 Like