[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.

27 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.

1 Like
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.

6 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

3 Likes