gmsh for rectangles

Author: DanielVandH

Project.toml

@@ -1,7 +1,7 @@
 name = "DelaunayTriangulation"
 uuid = "927a84f5-c5f4-47a5-9785-b46e178433df"
 authors = ["Daniel VandenHeuvel <[email protected]>"]
-version = "0.3.1"
+version = "0.3.2"
 
 [deps]
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"

README.md

@@ -123,7 +123,7 @@ poly!(ax, pts, [collect(T)[i][j] for i in 1:length(T), j in 1:3], color = (:whit
 
 ## Gmsh 
 
-Support is also added for a simple mesh generator with Gmsh (see https://gmsh.info/), tested up to v4.9.4 on Windows 64. The function for this is `generate_mesh`, and accepts inputs of boundary points that are in counter-clockwise order. This is especially useful for e.g. finite volume codes. Currently I only have code working for simply connected domains - it would be nice to have an alternative, but this is the best I can do with the time I have (the alternative would require me to think a lot more about ghost nodes, boundary edges, etc. when the domain has holes, and the impact this would have on the existing code and existing data structures).
+Support is also added for a simple mesh generator with Gmsh (see https://gmsh.info/), tested up to v4.9.4 on Windows 64. The function for this is `generate_mesh`, and accepts inputs of boundary points that are in counter-clockwise order. This is especially useful for e.g. finite volume codes. Currently I only have code working for simply connected domains - it would be nice to have an alternative, but this is the best I can do with the time I have (the alternative would require me to think a lot more about ghost nodes, boundary edges, etc. when the domain has holes, and the impact this would have on the existing code and existing data structures) -- this will eventually be supported.
 
 Let me give an example. In my directory, I have downloaded `gmsh` and saved it as `gmsh-4.9.4-Windows64`, so I define
 ```julia
@@ -185,6 +185,8 @@ lines!(ax, pts[:, BN[5]], color=:darkgreen, linewidth=4)
 
 As you can see, the triangulation successfully works on this non-convex geometry, and we have been able to clearly identify where in `pts` all the boundary nodes are. This makes it especially useful for example for applying boundary conditions in a finite volume code.
 
+You can also use `generate_mesh(a, b, c, d, ref)` to generate the mesh for a rectangle $[a, b] \times [c, d]$.
+
 # Customisation 
 
 The package is built to allow for customisation. For either `triangulate_bowyer` or `triangulate_berg`, we have the following keywords (these are not all of them):

src/data_structures.jl

@@ -421,4 +421,64 @@ function Base.iterate(tri::Triangulation, state=1)
     elseif state == 5
         return (get_points(tri), nothing)
     end
-end
+end
+
+"""
+    triangulate(triangles,
+        nodes::AbstractMatrix,
+        boundary_nodes; 
+        IntegerType::Type{I}=Int64,
+        EdgeType::Type{E}=NTuple{2,IntegerType},
+        TriangleType::Type{V}=NTuple{3,IntegerType},
+        EdgesType::Type{Es}=Set{EdgeType},
+        TrianglesType::Type{Ts}=Set{TriangleType},
+        sort_boundary=true) where {I,E,V,Es,Ts}
+
+Given an existing set of `triangles` (all in counter-clockwise order), `nodes`, and `boundary_nodes` (listed in clockwise order),
+puts them into a [`Triangulation`](@ref) struct.
+"""
+function triangulate(triangles,
+    nodes::AbstractMatrix,
+    boundary_nodes; # BOUNDARY NODES MUST BE CLOCKWISE
+    IntegerType::Type{I}=Int64,
+    EdgeType::Type{E}=NTuple{2,IntegerType},
+    TriangleType::Type{V}=NTuple{3,IntegerType},
+    EdgesType::Type{Es}=Set{EdgeType},
+    TrianglesType::Type{Ts}=Set{TriangleType},
+    sort_boundary=true) where {I,E,V,Es,Ts}
+    boundary_nodes = deepcopy(boundary_nodes)
+    boundary_nodes = reduce(vcat, boundary_nodes)
+    sort_boundary && unique!(boundary_nodes)
+    T = Ts()
+    adj = Adjacent{I,E}()
+    adj2v = Adjacent2Vertex{I,Es,E}()
+    DG = DelaunayGraph{I}()
+    itr = triangles isa AbstractMatrix ? eachcol(triangles) : triangles
+    for (i, j, k) in itr
+        add_triangle!(i, j, k, T, adj, adj2v, DG; protect_boundary=true)
+    end
+    for i in boundary_nodes
+        add_neighbour!(DG, I(BoundaryIndex), i)
+    end
+    if sort_boundary
+        cx, cy = sum(nodes; dims=2) / size(nodes, 2)
+        θ = zeros(length(boundary_nodes))
+        for (j, i) in pairs(boundary_nodes)
+            p = @view nodes[:, i]
+            x, y = p
+            θ[j] = atan(y - cy, x - cx)
+        end
+        sort_idx = sortperm(θ)
+        permute!(boundary_nodes, sort_idx)
+        reverse!(boundary_nodes) # cw 
+    end
+    n = length(boundary_nodes)
+    push!(boundary_nodes, boundary_nodes[begin])
+    for i in 1:n
+        u = boundary_nodes[i]
+        v = boundary_nodes[i+1]
+        add_edge!(adj2v, I(BoundaryIndex), u, v)
+        add_edge!(adj, u, v, I(BoundaryIndex))
+    end
+    return T, adj, adj2v, DG
+end