Gmsh

We provide Gmsh support, making it possible to represent more complicated meshes. (This is less relevant now that we have constrained triangulations and mesh refinement without needing Gmsh.) The discussion that follows assume you have installed Gmsh, and defined a corresponding GMSH_PATH. I have used the default,

julia> GMSH_PATH = "./gmsh-4.11.1-Windows64/gmsh.exe"

The methods for generating meshes with Gmsh are shown below.

DelaunayTriangulation.generate_meshFunction
generate_mesh(x, y, ref;
    mesh_algorithm = 6, 
    gmsh_path = "./gmsh-4.11.1-Windows64/gmsh.exe",
    verbosity = 0, 
    convert_result = true, 
    add_ghost_triangles = false, 
    check_args = true)

Using Gmsh, generates a mesh of the domain defined by (x, y).

Arguments

  • x, y: These are the coordinates defining the curves that define the boundaries of the domain. All curves are to be positively oriented, meaning the outermost boundary should be counter-clockwise while the interior boundaries should be clockwise. The accepted forms of x, y are outlined in in the boundary handling section of the docs.

  • ref: The refinement parameter – smaller ref means more elements.

Keyword Arguments

  • mesh_algorithm = 6: The meshing algorithm to use. See the Gmsh documentation for more information. The default 6 means Frontal-Delaunay.
  • gmsh_path = "./gmsh-4.11.1-Windows64/gmsh.exe": The path to the Gmsh executable.
  • verbosity = 0: The verbosity of the Gmsh output.
  • convert_result = true: Whether to convert the Gmsh output to a Triangulation.
  • add_ghost_triangles = false: Whether to add ghost triangles to the triangulation.
  • check_args = true: Whether to check the validity of the arguments.

Outputs

If convert_result, then the final result is a Triangulation type for the mesh. Otherwise, the following values are returned:

  • elements: The triangular elements of the mesh.
  • nodes: The nodes in the mesh.
  • boundary_nodes: The bonudary nodes in the mesh. All boundaries are positively oriented relative to the interior, meaning the outermost boundary is counter-clockwise while the interior boundaries are clockwise, and match the form of (x, y).

Extended help

The function proceeds in four steps:

  1. Mesh generation

Here, we write a file "meshgeometry.geo" in the working directory. This file takes the form

r = ref;
Mesh.Algorithm = mesh_algorithm; 
Mesh.Format = 1;
General.Verbosity = 0;
Point(<point index>) = {<x>, <y>, 0, r}; # For each point
Line(<line index>) = {<initial point>, <final point>}; # For each line 
Curve Loop(<boundary index>) = {<line 1>, <line 2>, ...}; # For each boundary 
Plane Surface(1) = {<curve 1>, <curve 2>, ...}; # <curve 1> = 1 and is the outermost boundary, while <curve i> = i, i > 1, are boundaries of interior holes 
Physical Curve(<last line index + i>) = {<line 1>, <line 2>, ...}; # For i ranging over the number of segments, and the lines represent that segment 
Physical Surface(1) = {1};

Most importantly, every edge input into the function generate_mesh will be included in the mesh. An older version of this function previously used cubic splines for defining boundary curves, but this has the consequence that (1) not every edge put into the function is included, and (2) the boundary is not exactly represented.

The function that handles this generation is write_gmsh.

  1. Mesh writing

The "meshgeometry.geo" file is then used to mesh the domain, running the terminal command

gmsh_path "meshgeometry.geo" -2 -format msh2

This creates a file "meshgeometry.msh" in the same working directory.

The function that handles this writing is run_gmsh.

  1. Mesh reading

Once "meshgeometry.msh" is created, we need to read it. The format used (MeshFormat) is 2.2, but note that as of writing (13/01/2013), the most modern format is 4.1.

The "meshgeometry.msh" file is split into groups:

3a. MeshFormat

This just reads off the format of the file used. This part of the file is read using read_mesh_format!.

3b. Nodes

This lists the node indices and all the coordinates of the nodes, with the first line giving the number of nodes. A single line in this section, after the first, takes the form

<node index> <x> <y> 0

and we read this using read_node_line. This entire part of the file is read using read_nodes!.

3c. Elements

The first line in this part of the file is the number of elements, though here elements refer to both the lines and the triangles. The lines (edges) are listed first, with each line taking the form

<line index> 1 2 <boundary index> <> <left node> <right node>

and will be in counter-clockwise order. After the lines are listed, all the triangles follow, with each line in this part taking the form

<triangle index> 2 2 1 1 <node 1> <node 2> <node 3>

with each triangle positively oriented. These lines are read using read_element_line. The entire part of the file is read using read_elements!.

  1. Conversion to Triangulation

Once the file "meshgeometry.geo" has been read, we have a list of triangular elements, nodes, and boundary nodes. These need to all be converted into a Triangulation type, and a constructor of Triangulation is used to accomplish this.

generate_mesh(a, b, c, d, ref; 
    mesh_algorithm=6,
    gmsh_path="./gmsh-4.11.1-Windows64/gmsh.exe",
    verbosity=0,
    single_boundary=true,
    convert_result=true,
    add_ghost_triangles=false)

Generates a mesh of a rectangle [a, b] × [c, d]. Use single_boundary=true if each side of the rectangle should be treated the same, and single_boundary=false if you want the boundary nodes to be segmented each side of the rectangle.

See the main function generate_mesh for a description of the other arguments.

Let's give some examples.

Example I: Contiguous boundary

Let us mesh a domain with a single non-segmented boundary curve.

using DelaunayTriangulation, CairoMakie
a = 4 / 5
t = LinRange(0, 2π, 100)
x = @. a * (2cos(t) + cos(2t))
y = @. a * (2sin(t) - sin(2t))
tri = generate_mesh(x, y, 0.1)
tri2 = generate_mesh(x, y, 1.0)
fig = Figure()
ax = Axis(fig[1, 1], xlabel=L"x", ylabel=L"y", width=300, height=300,
    title=L"(a):$ $ Dense mesh", titlealign=:left)
triplot!(ax, tri, show_convex_hull=true, show_constrained_edges=true)
ax = Axis(fig[1, 2], xlabel=L"x", ylabel=L"y", width=300, height=300,
    title=L"(b):$ $  Coarse mesh", titlealign=:left)
triplot!(ax, tri2, show_convex_hull=true, show_constrained_edges=true)
resize_to_layout!(fig)
Triangulation

In the figure, the red curve shows the convex hull. We note that we now have information in tri.boundary_nodes:

julia> get_boundary_nodes(tri)
178-element Vector{Int}:
  1
  2
  3
  4
  ⋮
 97
 98
 99
  1

Similarly, tri.boundary_map is now populated:

julia> get_boundary_map(tri)
OrderedDict{Int, Vector{Int}} with 1 entry:
  -1 => [1, 2, 3, 4, 5, 6, 100, 7, 101, 8  …  93, 176, 94, 177, 95, 96, 97, 98, 99, 1]

We now also have tri.boundary_edge_map:

julia> tri.boundary_edge_map
Dict{Tuple{Int, Int}, Tuple{Vector{Int}, Int}} with 177 entries:
  (116, 20)  => ([1, 2, 3, 4, 5, 6, 100, 7, 101, 8  …  93, 176, 94, 177, 95, 96, 97, 98, 99, 1], 36)
  (78, 158)  => ([1, 2, 3, 4, 5, 6, 100, 7, 101, 8  …  93, 176, 94, 177, 95, 96, 97, 98, 99, 1], 136)
  (11, 105)  => ([1, 2, 3, 4, 5, 6, 100, 7, 101, 8  …  93, 176, 94, 177, 95, 96, 97, 98, 99, 1], 16)
  (106, 13)  => ([1, 2, 3, 4, 5, 6, 100, 7, 101, 8  …  93, 176, 94, 177, 95, 96, 97, 98, 99, 1], 19)
  (103, 10)  => ([1, 2, 3, 4, 5, 6, 100, 7, 101, 8  …  93, 176, 94, 177, 95, 96, 97, 98, 99, 1], 13)
  (145, 56)  => ([1, 2, 3, 4, 5, 6, 100, 7, 101, 8  …  93, 176, 94, 177, 95, 96, 97, 98, 99, 1], 101)
  (169, 87)  => ([1, 2, 3, 4, 5, 6, 100, 7, 101, 8  …  93, 176, 94, 177, 95, 96, 97, 98, 99, 1], 156)
  (110, 111) => ([1, 2, 3, 4, 5, 6, 100, 7, 101, 8  …  93, 176, 94, 177, 95, 96, 97, 98, 99, 1], 27)
  (128, 42)  => ([1, 2, 3, 4, 5, 6, 100, 7, 101, 8  …  93, 176, 94, 177, 95, 96, 97, 98, 99, 1], 70)
  (43, 130)  => ([1, 2, 3, 4, 5, 6, 100, 7, 101, 8  …  93, 176, 94, 177, 95, 96, 97, 98, 99, 1], 73)
  (30, 31)   => ([1, 2, 3, 4, 5, 6, 100, 7, 101, 8  …  93, 176, 94, 177, 95, 96, 97, 98, 99, 1], 56)
  (156, 77)  => ([1, 2, 3, 4, 5, 6, 100, 7, 101, 8  …  93, 176, 94, 177, 95, 96, 97, 98, 99, 1], 133)
  (3, 4)     => ([1, 2, 3, 4, 5, 6, 100, 7, 101, 8  …  93, 176, 94, 177, 95, 96, 97, 98, 99, 1], 3)
  (112, 113) => ([1, 2, 3, 4, 5, 6, 100, 7, 101, 8  …  93, 176, 94, 177, 95, 96, 97, 98, 99, 1], 30)
  (41, 128)  => ([1, 2, 3, 4, 5, 6, 100, 7, 101, 8  …  93, 176, 94, 177, 95, 96, 97, 98, 99, 1], 69)
  (153, 74)  => ([1, 2, 3, 4, 5, 6, 100, 7, 101, 8  …  93, 176, 94, 177, 95, 96, 97, 98, 99, 1], 127)
  (133, 47)  => ([1, 2, 3, 4, 5, 6, 100, 7, 101, 8  …  93, 176, 94, 177, 95, 96, 97, 98, 99, 1], 80)
  ⋮          => ⋮

In this case, each output of (i, j) is the Tuple (get_boundary_nodes(tri), k). For example,

julia> u, v = 133, 47;

julia> pos = get_boundary_edge_map(tri, u, v);

julia> segment_nodes = get_boundary_nodes(tri, pos[1]);

julia> get_boundary_nodes(segment_nodes, pos[2]) == u # edges start at the left
true

julia> get_boundary_nodes(segment_nodes, pos[2]+1) == v
true

Example II: Single boundary curve with multiple segments

Let us now give an example where we still have just a single boundary curve, but we split it into multiple segments. Importantly, each segment must be counter-clockwise and join with the previous segment.

using DelaunayTriangulation, CairoMakie
# The first segment 
t = LinRange(0, 1 / 4, 25)
x1 = cos.(2π * t)
y1 = sin.(2π * t)
# The second segment 
t = LinRange(0, -3, 25)
x2 = collect(t)
y2 = repeat([1.0], length(t))
# The third segment 
t = LinRange(1, 0, 25)
x3 = -3.0 .+ (1 .- t) .* sin.(t)
y3 = collect(t)
# The fourth segment 
t = LinRange(0, 1, 25)
x4 = collect(-3.0(1 .- t))
y4 = collect(0.98t)
# The fifth segment 
x5 = [0.073914, 0.0797, 0.1522, 0.1522, 0.2, 0.28128, 0.3659, 0.4127, 0.3922, 0.4068, 0.497, 0.631, 0.728, 0.804, 0.888, 1.0]
y5 = [0.8815, 0.8056, 0.80268, 0.73258, 0.6, 0.598, 0.5777, 0.525, 0.4346, 0.3645, 0.3032, 0.2886, 0.2623, 0.1367, 0.08127, 0.0]
# Now combine the vectors 
x = [x1, x2, x3, x4, x5]
y = [y1, y2, y3, y4, y5]
# Mesh 
tri = generate_mesh(x, y, 0.05)
fig = Figure()
ax = Axis(fig[1, 1], xlabel=L"x", ylabel=L"y", width=600, height=300)
triplot!(ax, tri, show_convex_hull=true)
colors = [:red, :blue, :orange, :purple, :darkgreen]
bn_map = get_boundary_map(tri)
for (i, segment_index) in enumerate(values(bn_map))
    bn_nodes = get_boundary_nodes(tri, segment_index)
    lines!(ax, get_points(tri)[:, bn_nodes], color=colors[i], linewidth=4)
end
resize_to_layout!(fig)
Triangulation

An important feature to note for this case is that the index now used to refer to boundaries is no longer constant. In particular, the nth segment will map to -n, as we can verify using:

julia> get_boundary_map(tri)
OrderedDict{Int, Int} with 5 entries:
  -1 => 1
  -2 => 2
  -3 => 3
  -4 => 4
  -5 => 5

This map makes it simple to iterate over all parts of a boundary, as we show in the above code when plotting. The tri.boundary_edge_map in this case is:

julia> get_boundary_edge_map(tri)
Dict{Tuple{Int, Int}, Tuple{Int, Int}} with 262 entries:
  (118, 8)   => (1, 14)
  (55, 56)   => (3, 10)
  (34, 154)  => (2, 28)
  (213, 214) => (4, 35)
  (223, 224) => (4, 50)
  (143, 29)  => (2, 12)
  (24, 135)  => (1, 47)
  (178, 179) => (2, 65)
  (132, 22)  => (1, 42)
  (46, 178)  => (2, 64)
  (169, 42)  => (2, 51)
  (154, 155) => (2, 29)
  (115, 5)   => (1, 8)
  (43, 172)  => (2, 55)
  (261, 262) => (5, 38)
  (146, 147) => (2, 17)
  (49, 184)  => (3, 1)
  ⋮          => ⋮

For example,

julia> let pos = get_boundary_edge_map(tri, 115, 5)
       segment_nodes = get_boundary_nodes(tri, pos[1])
       u′ = get_boundary_nodes(segment_nodes, pos[2])
       v′ = get_boundary_nodes(segment_nodes, pos[2]+1)
       u′ == 115 && v′ == 5
       end
true

Example III: Multiple boundaries

Now let us give a more complicated example, meshing a multiply-connected domain. In this case, we provide the outer-most boundary in a counter-clockwise order, while all the inner boundaries are in a clockwise order.

using DelaunayTriangulation, CairoMakie
x1 = [collect(LinRange(0, 2, 4)),
    collect(LinRange(2, 2, 4)),
    collect(LinRange(2, 0, 4)),
    collect(LinRange(0, 0, 4))]
y1 = [collect(LinRange(0, 0, 4)),
    collect(LinRange(0, 6, 4)),
    collect(LinRange(6, 6, 4)),
    collect(LinRange(6, 0, 4))]
r = 0.5
h = k = 0.6
θ = LinRange(2π, 0, 50)
x2 = [h .+ r .* cos.(θ)]
y2 = [k .+ r .* sin.(θ)]
r = 0.2
h = 1.5
k = 0.5
x3 = [h .+ r .* cos.(θ)]
y3 = [k .+ r .* sin.(θ)]
x4 = reverse(reverse.([collect(LinRange(1, 1.5, 4)),
    collect(LinRange(1.5, 1.5, 4)),
    collect(LinRange(1.5, 1, 4)),
    collect(LinRange(1, 1, 4))]))
y4 = reverse(reverse.([collect(LinRange(2, 2, 4)),
    collect(LinRange(2, 5, 4)),
    collect(LinRange(5, 5, 4)),
    collect(LinRange(5, 2, 4))]))
x5 = [reverse([0.2, 0.5, 0.75, 0.75, 0.2, 0.2])]
y5 = [reverse([2.0, 2.0, 3.0, 4.0, 5.0, 2.0])]
x = [x1, x2, x3, x4, x5]
y = [y1, y2, y3, y4, y5]
tri = generate_mesh(x, y, 0.2)
fig, ax, sc = triplot(tri; show_convex_hull=true, show_ghost_edges=true, show_constrained_edges=true, convex_hull_linestyle=:solid, convex_hull_linewidth=4)
Triangulation

The blue edges show the interpretation of the ghost edges (you can delete via delete_ghost_triangles! if you want). For the outer boundary, these edges are pointing away from the interior, collinear with a point in the center, as we can obtain via:

julia> DelaunayTriangulation.get_representative_point_coordinates(tri, 1)

or, alternatively,

julia> DelaunayTriangulation.get_representative_point_list(tri)
Dict{Int, DelaunayTriangulation.RepresentativeCoordinates{Int, Float64}} with 5 entries:
  5 => RepresentativeCoordinates{Int, Float64}(0.475, 3.5, 0)
  4 => RepresentativeCoordinates{Int, Float64}(1.25, 3.5, 0)
  2 => RepresentativeCoordinates{Int, Float64}(0.6, 0.6, 0)
  3 => RepresentativeCoordinates{Int, Float64}(1.5, 0.5, 0)
  1 => RepresentativeCoordinates{Int, Float64}(1.5, 1.5, 0)

The keys are the indices for the boundary curve. These coordinates are visual centers, obtained via the pole of inaccessibility function; see the sidebar. For the inner boundaries, the ghost edges are no longer infinite and so they connect directly with these representative coordinates.

To access more of the boundary information, we could first consider boundary_nodes:

julia> get_boundary_nodes(tri)
5-element Vector{Vector{Vector{Int}}}:
 [[1, 128, 129, 130, 2, 131, 132, 133, 3, 134, 135, 136, 4], [4, 137, 138, 139, 140, 141, 142, 143, 144, 145  …  155, 156, 157, 158, 159, 160, 161, 162, 163, 7], [7, 164, 165, 166, 8, 167, 168, 169, 9, 170, 171, 172, 10], [10, 173, 174, 175, 176, 177, 178, 179, 180, 181  …  191, 192, 193, 194, 195, 196, 197, 198, 199, 1]]
 [[13, 14, 15, 16, 17, 18, 19, 20, 21, 22  …  53, 54, 55, 56, 57, 58, 59, 60, 61, 13]]
 [[62, 63, 64, 65, 66, 67, 68, 69, 70, 71  …  102, 103, 104, 105, 106, 107, 108, 109, 110, 62]]
 [[111, 200, 201, 202, 203, 112, 204, 205, 206, 207, 113, 208, 209, 210, 211, 114], [114, 115, 116, 117], [117, 212, 213, 214, 215, 118, 216, 217, 218, 219, 119, 220, 221, 222, 223, 120], [120, 121, 122, 111]]  
 [[123, 224, 225, 226, 227, 228, 229, 230, 231, 232  …  246, 126, 247, 248, 249, 250, 251, 127, 252, 123]]

This is simply a vector of curves, with each curve storing its segments. This vector itself does not tell us what boundary index corresponds to what segment of what curve, and this could be obtained from the boundary map:

julia> get_boundary_map(tri)
OrderedDict{Int, Tuple{Int, Int}} with 11 entries:
  -1  => (1, 1)
  -2  => (1, 2)
  -3  => (1, 3)
  -4  => (1, 4)
  -5  => (2, 1)
  -6  => (3, 1)
  -7  => (4, 1)
  -8  => (4, 2)
  -9  => (4, 3)
  -10 => (4, 4)
  -11 => (5, 1)

So, for example, the boundary index -8 comes from the second segment of the fourth curve. As before, this boundary map makes it simple to iterate over each segment as follows:

bn_map = get_boundary_map(tri)
for segment_index in values(bn_map)
    bn_nodes = get_boundary_nodes(tri, segment_index)
    nedges = num_boundary_edges(bn_nodes) # Note that nedges = length(bn_nodes) - 1
    for edge_idx in 1:nedges 
        node = get_boundary_node(bn_nodes, edge_idx)
        ...
    end
end

The form above is generic, and ignores the last part of each segment (since it is duplicated for the next segment). Of course, a version like

bn_map = get_boundary_map(tri)
for segment_index in values(bn_map)
    bn_nodes = get_boundary_nodes(tri, segment_index)
    for i in bn_nodes
        node = get_boundary_node(bn_nodes, i)
        ...
    end
end

(which includes the last part of each segment) could be used. It is up to you based on your interface how you prefer to write this. Notice also that in the previous example we used a similar style, using get_boundary_nodes(tri, segment_index) also. The function get_boundary_nodes can be used with either single integers or Tuples, making it simple to iterate with this exact pattern whether we have a contiguous boundary curve, a segmented boundary curve, or multiple boundaries.

Another feature to note is tri.boundary_index_ranges, which will tell us what other boundary indices belong to a curve given a known boundary index for that curve. This can be useful if we want to rotate around a boundary curve based on a given boundary index (see e.g. how it is used in the get_left_boundary_node and get_right_boundary_node functions). This field is a major part of making point location work in these inner boundaries, making get_adjacent work properly in this case (see e.g. the code in _safe_get_adjacent).

julia> get_boundary_index_ranges(tri)
OrderedDict{Int, UnitRange{Int}} with 11 entries:
  -1  => -4:-1
  -2  => -4:-1
  -3  => -4:-1
  -4  => -4:-1
  -5  => -5:-5
  -6  => -6:-6
  -7  => -10:-7
  -8  => -10:-7
  -9  => -10:-7
  -10 => -10:-7
  -11 => -11:-11

So, for example, we see tha the boundary index -3 belongs to a curve that also has boundary indices -1, -2, and -4. If we wanted to go from a boundary index to the index for the curve, this is what the boundary map is also for:

julia> DelaunayTriangulation.get_curve_index(tri, -3)
1

The last feature to show is the new tri.boundary_edge_map for this case, given by

julia> get_boundary_edge_map(tri)
Dict{Tuple{Int, Int}, Tuple{Tuple{Int, Int}, Int}} with 252 entries:
  (55, 56)   => ((2, 1), 43)
  (130, 2)   => ((1, 1), 4)
  (92, 93)   => ((3, 1), 31)
  (213, 214) => ((4, 3), 3)
  (14, 15)   => ((2, 1), 2)
  (172, 10)  => ((1, 3), 12)
  (203, 112) => ((4, 1), 5)
  (178, 179) => ((1, 4), 7)
  (121, 122) => ((4, 4), 2)
  (151, 152) => ((1, 2), 17)
  (26, 27)   => ((2, 1), 14)
  (171, 172) => ((1, 3), 11)
  (88, 89)   => ((3, 1), 27)
  (132, 133) => ((1, 1), 7)
  (133, 3)   => ((1, 1), 8)
  (146, 147) => ((1, 2), 12)
  (46, 47)   => ((2, 1), 34)
  ⋮          => ⋮

As before, (u, v) = (get_boundary_nodes(segment_nodes, pos[2]), get_boundary_nodes(segment_nodes, pos[2]+1), where segment_nodes = get_boundary_nodes(tri, pos[1]) and pos = get_boundary_edge_map(tri, u, v). This pattern is true for any form of boundary nodes, in fact.