Mesh Refinement
In this tutorial, we show how we can use mesh refinement to improve the quality of a mesh by inserting more points. In this package, we allow for constant area and angle constraints, and also for custom constraints based on a user-provided function. You can also limit the maximum number of points. Moreover, any type of triangulation can be provided, regardless of whether the triangulation is unconstrained or triangulations, and regardless of the number of holes and domains. Curve-bounded domains, not included in this tutorial, can also be refined as discussed in this tutorial.
Let us start by loading in the packages we will need.
using DelaunayTriangulation
using CairoMakie
using StableRNGs
Unconstrained triangulation
Let us start with a simple example, refining an unconstrained triangulation. We will constrain the triangulation such that the minimum angle is 30 degrees, and the maximum area of a triangulation is 1% of the triangulation's total area. Note that below we need to make sure points
is mutable, else it is not possible to push points into the triangulation. Here we use a vector, but you could also use e.g. an ElasticMatrix
from ElasticArrays.jl.
rng = StableRNG(123)
x = rand(rng, 50)
y = rand(rng, 50)
points = tuple.(x, y)
tri = triangulate(points; rng)
orig_tri = deepcopy(tri)
A = get_area(tri)
refine!(tri; min_angle = 30.0, max_area = 0.01A, rng)
Delaunay Triangulation.
Number of vertices: 293
Number of triangles: 522
Number of edges: 814
Has boundary nodes: false
Has ghost triangles: true
Curve-bounded: false
Weighted: false
Constrained: false
The refine!
function operates on tri
in-place. If we wanted to review the statistics of the refined mesh, we can use statistics
:
statistics(tri)
Delaunay Triangulation Statistics.
Triangulation area: 0.6676219360027273
Number of vertices: 294
Number of solid vertices: 293
Number of ghost vertices: 1
Number of edges: 876
Number of solid edges: 814
Number of ghost edges: 62
Number of triangles: 584
Number of solid triangles: 522
Number of ghost triangles: 62
Number of boundary segments: 0
Number of interior segments: 0
Number of segments: 0
Number of convex hull vertices: 62
Smallest angle: 30.093123918976033°
Largest angle: 115.16815010364209°
Smallest area: 2.1013155442142383e-5
Largest area: 0.006504759798840958
Smallest radius-edge ratio: 0.5783859342412219
Largest radius-edge ratio: 0.9971940823138769
As we can see, the maximum area of a triangle is about 0.0064, which is indeed less than 1% of the triangulation's area, which is about 0.0067. Moreover, the smallest angle is indeed greater than 30.
Let us compare the triangulation pre- and post-refinement.
fig, ax, sc = triplot(orig_tri, axis = (title = "Pre-refinement",))
ax = Axis(fig[1, 2], title = "Post-refinement")
triplot!(ax, tri)
fig
The triangulation is now much finer. There are still some parts with many more triangles than other regions, but these are mostly near a boundary or where was a cluster of random points. If we wanted, we could refine again to try and improve this.
refine!(tri; min_angle = 30.0, max_area = 0.001A, rng) # 0.1% instead of 1%
fig, ax, sc = triplot(tri)
fig
The quality has now been improved. We could also try improving the minimum angle further, but even 30 is a bit closer to the limit of convergence (which is about 33.9 degrees). For example, if we try a minimum angle of 35 degrees, the algorithm just doesn't even converge, instead it reaches the maximum number of points.
test_tri = deepcopy(tri)
refine!(test_tri; min_angle = 35.0, max_area = 0.001A, max_points = 5_000, rng) # 20_000 so that it doesn't just keep going
statistics(test_tri)
Delaunay Triangulation Statistics.
Triangulation area: 0.6676219360027267
Number of vertices: 5001
Number of solid vertices: 5000
Number of ghost vertices: 1
Number of edges: 14997
Number of solid edges: 14829
Number of ghost edges: 168
Number of triangles: 9998
Number of solid triangles: 9830
Number of ghost triangles: 168
Number of boundary segments: 0
Number of interior segments: 0
Number of segments: 0
Number of convex hull vertices: 168
Smallest angle: 32.569525342355476°
Largest angle: 113.43199723096218°
Smallest area: 6.547564216329479e-6
Largest area: 0.0006565996109229845
Smallest radius-edge ratio: 0.5773550586596291
Largest radius-edge ratio: 0.9288110493512934
As we can see, the smallest angle is about 29 degrees instead of 35 degrees, and there are now 5000 points in the triangulation. The resulting triangulation is given below:
fig, ax, sc = triplot(test_tri)
fig
This is certainly not a suitable triangulation.
One useful figure to look at for these plots are histograms that look at the areas and angles. Looking to tri
, we can plot these as follows:
stats = statistics(tri)
fig = Figure(fontsize = 33)
areas = get_all_stat(stats, :area) ./ A
angles = first.(get_all_stat(stats, :angles)) # the first is the smallest
ax = Axis(fig[1, 1], xlabel = "A/A(Ω)", ylabel = "Count", title = "Area histogram", width = 400, height = 400, titlealign = :left)
hist!(ax, areas, bins = 0:0.0001:0.0005)
ax = Axis(fig[1, 2], xlabel = "θₘᵢₙ", ylabel = "Count", title = "Angle histogram", width = 400, height = 400, titlealign = :left)
hist!(ax, rad2deg.(angles), bins = 20:2:60)
vlines!(ax, [30.0], color = :red)
resize_to_layout!(fig)
fig
We see that indeed many of the triangle areas are very small, and the angles are all greater than 30 degrees.
Constrained triangulation and custom constraints
We now give an example of a constrained triangulation being refined. This is the most common case where the mesh refinement is needed. For this example, we consider an example with holes, but note that any triangulation can be refined, regardless of the type. Here is the triangulation we consider.
n = 100
θ = LinRange(0, 2π, n + 1)
θ = [θ[1:n]; 0]
rev_θ = reverse(θ) # to go from ccw to cw
r₁ = 10.0
r₂ = 5.0
r₃ = 2.5
outer_x, outer_y = r₁ * cos.(θ), r₁ * sin.(θ)
inner_x, inner_y = r₂ * cos.(rev_θ), r₂ * sin.(rev_θ)
innermost_x, innermost_y = r₃ * cos.(θ), r₃ * sin.(θ)
x = [[outer_x], [inner_x], [innermost_x]]
y = [[outer_y], [inner_y], [innermost_y]]
boundary_nodes, points = convert_boundary_points_to_indices(x, y)
rng = StableRNG(456)
tri = triangulate(points; boundary_nodes, rng)
fig, ax, sc = triplot(tri)
fig
Let us now refine this triangulation.
A = get_area(tri)
refine!(tri; min_angle = 27.3, max_area = 0.01A, rng)
fig, ax, sc = triplot(tri)
fig
We inspect the plot, and we might think that it's perhaps not fine enough. Let's use finer constraints and see what happens. Since refine!
operates on tri
in-place, refining it again with the constraints below is going to take roughly the same amount of time as if we had refined it with these constraints in the first place.
refine!(tri; min_angle = 33.9, max_area = 0.001A, rng)
fig, ax, sc = triplot(tri)
fig
This is indeed much better, but notice that the inner hole is much more fine than the outer. This is because we are applying the same area constraint inside and outside, when really we should try and take note of the total contribution to the area that each part of the domain gives. To refine with this in mind, we need to use custom constraints. The function that we use for constraining the area takes the form f(tri, T)
, where T
is the triangle's vertices and tri
is the triangulation. It should return true
if the triangle should be refined, and false
otherwise. Let us define a function such that, instead of applying constraints so that the triangles are limited to 1% of the total triangulation area, we do 0.5% or 0.1% of the area of the inner or outer domain, respectively.
outer_area = π * (r₁^2 - r₂^2)
inner_area = π * r₃^2
function in_inner(p, q, r)
px, py = getxy(p)
qx, qy = getxy(q)
rx, ry = getxy(r)
cx, cy = (px + qx + rx) / 3, (py + qy + ry) / 3
rad2 = cx^2 + cy^2
return rad2 ≤ r₃^2
end
function area_constraint(_tri, T)
i, j, k = triangle_vertices(T)
p, q, r = get_point(_tri, i, j, k)
A = DelaunayTriangulation.triangle_area(p, q, r)
return in_inner(p, q, r) ? (A ≥ 0.005inner_area) : (A ≥ 0.001outer_area)
end
area_constraint (generic function with 1 method)
Let's now refine. We recompute the triangulation so that we can see the new results.
boundary_nodes, points = convert_boundary_points_to_indices(x, y)
rng = StableRNG(456)
tri = triangulate(points; boundary_nodes, rng)
refine!(tri; min_angle = 30.0, custom_constraint = area_constraint, rng)
fig, ax, sc = triplot(tri)
fig
This is now much better, and the two parts of the domain are appropriately refined. Let us extend our custom constraint function to also require that any triangle has minimum angle less than 33 degrees inside the innermost domain, and less than 20 degrees outside the innermost domain.
function angle_constraint(_tri, T)
i, j, k = triangle_vertices(T)
p, q, r = get_point(_tri, i, j, k)
θ = rad2deg(minimum(DelaunayTriangulation.triangle_angles(p, q, r)))
return in_inner(p, q, r) ? (θ ≤ 33.9) : (θ ≤ 20.0)
end
function custom_constraint(_tri, T)
return area_constraint(_tri, T) || angle_constraint(_tri, T)
end
boundary_nodes, points = convert_boundary_points_to_indices(x, y)
rng = StableRNG(456)
tri = triangulate(points; boundary_nodes, rng)
refine!(tri; custom_constraint, rng)
fig, ax, sc = triplot(tri)
fig
Indeed, the inner domain is much finer. These examples could be extended to more complicated cases, for example using adaptive mesh refinement for a numerical PDE solution so that triangles are refined based on some a posteriori error estimate, implemented using a custom area constraint like above, or even some refinement based on the triangle's location in space in case of some geospatial application.
Domains with small angles
In the examples considered, none of the boundaries had small angles. When domains have small angles, it is not always possible to satisfy the minimum angle constraints, but the algorithm will still try its best to refine in these locations. Let's consider a complicated example with many small angles. We consider the boundary of Switzerland, as obtained in this NaturalNeighbours.jl example.
using Downloads
using DelimitedFiles
boundary_url = "https://gist.githubusercontent.com/DanielVandH/13687b0918e45a416a5c93cd52c91449/raw/a8da6cdc94859fd66bcff85a2307f0f9cd57a18c/boundary.txt"
boundary_dir = Downloads.download(boundary_url)
boundary = readdlm(boundary_dir, skipstart = 6)
boundary_points = [(boundary[i, 1], boundary[i, 2]) for i in axes(boundary, 1)]
reverse!(boundary_points)
5126-element Vector{Tuple{Float64, Float64}}:
(7.133862520764, 45.871870678164)
(7.138563027059, 45.871870678164)
(7.138563027059, 45.876571184459)
(7.143263533355, 45.876571184459)
(7.14796403965, 45.876571184459)
(7.152664545945, 45.876571184459)
(7.15736505224, 45.876571184459)
(7.162065558535, 45.876571184459)
(7.162065558535, 45.871870678164)
(7.166766064831, 45.871870678164)
⋮
(7.110359989288, 45.857769159278)
(7.115060495583, 45.857769159278)
(7.119761001878, 45.857769159278)
(7.119761001878, 45.862469665573)
(7.124461508174, 45.862469665573)
(7.129162014469, 45.862469665573)
(7.129162014469, 45.867170171868)
(7.133862520764, 45.867170171868)
(7.133862520764, 45.871870678164)
Here is the boundary.
boundary_nodes, points = convert_boundary_points_to_indices(boundary_points)
rng = StableRNG(789)
tri = triangulate(points; boundary_nodes, rng)
fig, ax, sc = triplot(tri)
fig
Now let's refine.
A = get_area(tri)
refine!(tri; min_angle = 30.0, max_area = 0.001A, rng)
Delaunay Triangulation.
Number of vertices: 14165
Number of triangles: 23178
Number of edges: 37342
Has boundary nodes: true
Has ghost triangles: true
Curve-bounded: false
Weighted: false
Constrained: true
fig, ax, sc = triplot(tri)
fig
We see that the triangulation is now adequately refined. There are still triangles near the boundaries whose minimum angle is less than 30 degrees, though, because of the angles that boundary edges meet at in some places. Most of the triangles will satisfy the constraint, though, as we show below.
stats = statistics(tri)
angles = first.(get_all_stat(stats, :angles)) # the first is the smallest
fig, ax, sc = scatter(rad2deg.(angles))
hlines!(ax, [30.0], color = :red, linewidth = 4)
fig
As we can see, the vast majority of the triangles satisfy the constraint, but there are still some that do not. Here is another set of results with a lower minimum angle constraint.
boundary_nodes, points = convert_boundary_points_to_indices(boundary_points)
rng = StableRNG(789)
tri = triangulate(points; boundary_nodes, rng)
refine!(tri; min_angle = 18.73, max_area = 0.001A, rng)
fig = Figure(fontsize = 43)
ax = Axis(fig[1, 1], width = 600, height = 400)
triplot!(tri)
ax = Axis(fig[1, 2], width = 600, height = 400)
stats = statistics(tri)
angles = first.(get_all_stat(stats, :angles)) # the first is the smallest
scatter!(ax, rad2deg.(angles))
hlines!(ax, [18.73], color = :red, linewidth = 4)
resize_to_layout!(fig)
fig
In this case, all the triangles satisfy the constraint, of course at the expense of some other triangles having lesser quality.
Just the code
An uncommented version of this example is given below. You can view the source code for this file here.
using DelaunayTriangulation
using CairoMakie
using StableRNGs
rng = StableRNG(123)
x = rand(rng, 50)
y = rand(rng, 50)
points = tuple.(x, y)
tri = triangulate(points; rng)
orig_tri = deepcopy(tri)
A = get_area(tri)
refine!(tri; min_angle = 30.0, max_area = 0.01A, rng)
statistics(tri)
fig, ax, sc = triplot(orig_tri, axis = (title = "Pre-refinement",))
ax = Axis(fig[1, 2], title = "Post-refinement")
triplot!(ax, tri)
fig
refine!(tri; min_angle = 30.0, max_area = 0.001A, rng) # 0.1% instead of 1%
fig, ax, sc = triplot(tri)
fig
test_tri = deepcopy(tri)
refine!(test_tri; min_angle = 35.0, max_area = 0.001A, max_points = 5_000, rng) # 20_000 so that it doesn't just keep going
statistics(test_tri)
fig, ax, sc = triplot(test_tri)
fig
stats = statistics(tri)
fig = Figure(fontsize = 33)
areas = get_all_stat(stats, :area) ./ A
angles = first.(get_all_stat(stats, :angles)) # the first is the smallest
ax = Axis(fig[1, 1], xlabel = "A/A(Ω)", ylabel = "Count", title = "Area histogram", width = 400, height = 400, titlealign = :left)
hist!(ax, areas, bins = 0:0.0001:0.0005)
ax = Axis(fig[1, 2], xlabel = "θₘᵢₙ", ylabel = "Count", title = "Angle histogram", width = 400, height = 400, titlealign = :left)
hist!(ax, rad2deg.(angles), bins = 20:2:60)
vlines!(ax, [30.0], color = :red)
resize_to_layout!(fig)
fig
n = 100
θ = LinRange(0, 2π, n + 1)
θ = [θ[1:n]; 0]
rev_θ = reverse(θ) # to go from ccw to cw
r₁ = 10.0
r₂ = 5.0
r₃ = 2.5
outer_x, outer_y = r₁ * cos.(θ), r₁ * sin.(θ)
inner_x, inner_y = r₂ * cos.(rev_θ), r₂ * sin.(rev_θ)
innermost_x, innermost_y = r₃ * cos.(θ), r₃ * sin.(θ)
x = [[outer_x], [inner_x], [innermost_x]]
y = [[outer_y], [inner_y], [innermost_y]]
boundary_nodes, points = convert_boundary_points_to_indices(x, y)
rng = StableRNG(456)
tri = triangulate(points; boundary_nodes, rng)
fig, ax, sc = triplot(tri)
fig
A = get_area(tri)
refine!(tri; min_angle = 27.3, max_area = 0.01A, rng)
fig, ax, sc = triplot(tri)
fig
refine!(tri; min_angle = 33.9, max_area = 0.001A, rng)
fig, ax, sc = triplot(tri)
fig
outer_area = π * (r₁^2 - r₂^2)
inner_area = π * r₃^2
function in_inner(p, q, r)
px, py = getxy(p)
qx, qy = getxy(q)
rx, ry = getxy(r)
cx, cy = (px + qx + rx) / 3, (py + qy + ry) / 3
rad2 = cx^2 + cy^2
return rad2 ≤ r₃^2
end
function area_constraint(_tri, T)
i, j, k = triangle_vertices(T)
p, q, r = get_point(_tri, i, j, k)
A = DelaunayTriangulation.triangle_area(p, q, r)
return in_inner(p, q, r) ? (A ≥ 0.005inner_area) : (A ≥ 0.001outer_area)
end
boundary_nodes, points = convert_boundary_points_to_indices(x, y)
rng = StableRNG(456)
tri = triangulate(points; boundary_nodes, rng)
refine!(tri; min_angle = 30.0, custom_constraint = area_constraint, rng)
fig, ax, sc = triplot(tri)
fig
function angle_constraint(_tri, T)
i, j, k = triangle_vertices(T)
p, q, r = get_point(_tri, i, j, k)
θ = rad2deg(minimum(DelaunayTriangulation.triangle_angles(p, q, r)))
return in_inner(p, q, r) ? (θ ≤ 33.9) : (θ ≤ 20.0)
end
function custom_constraint(_tri, T)
return area_constraint(_tri, T) || angle_constraint(_tri, T)
end
boundary_nodes, points = convert_boundary_points_to_indices(x, y)
rng = StableRNG(456)
tri = triangulate(points; boundary_nodes, rng)
refine!(tri; custom_constraint, rng)
fig, ax, sc = triplot(tri)
fig
using Downloads
using DelimitedFiles
boundary_url = "https://gist.githubusercontent.com/DanielVandH/13687b0918e45a416a5c93cd52c91449/raw/a8da6cdc94859fd66bcff85a2307f0f9cd57a18c/boundary.txt"
boundary_dir = Downloads.download(boundary_url)
boundary = readdlm(boundary_dir, skipstart = 6)
boundary_points = [(boundary[i, 1], boundary[i, 2]) for i in axes(boundary, 1)]
reverse!(boundary_points)
boundary_nodes, points = convert_boundary_points_to_indices(boundary_points)
rng = StableRNG(789)
tri = triangulate(points; boundary_nodes, rng)
fig, ax, sc = triplot(tri)
fig
A = get_area(tri)
refine!(tri; min_angle = 30.0, max_area = 0.001A, rng)
fig, ax, sc = triplot(tri)
fig
stats = statistics(tri)
angles = first.(get_all_stat(stats, :angles)) # the first is the smallest
fig, ax, sc = scatter(rad2deg.(angles))
hlines!(ax, [30.0], color = :red, linewidth = 4)
fig
boundary_nodes, points = convert_boundary_points_to_indices(boundary_points)
rng = StableRNG(789)
tri = triangulate(points; boundary_nodes, rng)
refine!(tri; min_angle = 18.73, max_area = 0.001A, rng)
fig = Figure(fontsize = 43)
ax = Axis(fig[1, 1], width = 600, height = 400)
triplot!(tri)
ax = Axis(fig[1, 2], width = 600, height = 400)
stats = statistics(tri)
angles = first.(get_all_stat(stats, :angles)) # the first is the smallest
scatter!(ax, rad2deg.(angles))
hlines!(ax, [18.73], color = :red, linewidth = 4)
resize_to_layout!(fig)
fig
This page was generated using Literate.jl.