From 48d23023cf72d143e20f6aa7a163c1037731762d Mon Sep 17 00:00:00 2001 From: Daniel Sharp Date: Wed, 15 Feb 2023 16:05:35 -0500 Subject: [PATCH 1/4] Add to_mesh --- src/DistMesh2D.jl | 3 +- src/to_mesh.jl | 87 +++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 89 insertions(+), 1 deletion(-) create mode 100644 src/to_mesh.jl diff --git a/src/DistMesh2D.jl b/src/DistMesh2D.jl index 2c8333a..55c05c9 100644 --- a/src/DistMesh2D.jl +++ b/src/DistMesh2D.jl @@ -11,8 +11,9 @@ include("pointstoforces.jl") include("finalpoints.jl") include("triangulationexception.jl") include("distmesh.jl") +include("to_mesh.jl") export distmesh2d, drectangle, dcircle, TriangulationException, ddiff, dunion, -dintersect, huniform, protate +dintersect, huniform, protate, to_mesh end diff --git a/src/to_mesh.jl b/src/to_mesh.jl new file mode 100644 index 0000000..12d8f3d --- /dev/null +++ b/src/to_mesh.jl @@ -0,0 +1,87 @@ +""" + to_graph(x,y) +Assumes that `x,y` are of the form [a,b,NaN,c,d,NaN,....], representing line segments +""" +function to_graph(x,y) + Z = [x y] + Z_pts = Vector{Tuple{Float64,Float64}}(undef, 2size(Z,1)÷3) + for j in 1:size(Z,1)÷3 + Z_pts[2j-1] = (Z[3j-2,1],Z[3j-2,2]) + Z_pts[2j] = (Z[3j-1,1],Z[3j-1,1]) + end + Z_pts = [(z[1],z[2]) for z in eachrow(Z[Not(3:3:end),:])] + verts = unique(Z_pts) + g=[Int[] for _ in eachindex(verts)] + N_edges = size(Z_pts,1)÷2 + edges = Vector{Tuple{Int,Int}}(undef, N_edges) + for j in eachindex(edges) + z1,z2 = Z_pts[[2j-1,2j]] + z1_v = findfirst(==(z1), verts) + z2_v = findfirst(==(z2), verts) + push!(g[z1_v], z2_v) + push!(g[z2_v], z1_v) + end + verts, g +end + +function sort3(a,b,c) + t = a + if a>b + a = b + b = t + end + if b > c + t = b + b = c + c = t + end + if a > b + t = a + a = b + b = t + end + a,b,c +end + +function counterclockwise(v1::Int,v2::Int,v3::Int, verts) + v1,v2,v3 = sort3(v1,v2,v3) + ax,ay = verts[v1] + bx,by = verts[v2] + cx,cy = verts[v3] + det_tri = (bx-ax)*(cy-ay)-(cx-ax)*(by-ay) + is_cc = det_tri > 0 + is_cc ? (v1,v2,v3) : (v2,v1,v3) +end + +function tris_from_graph(g::Vector{Vector{Int}}, verts) + tris = NTuple{3,Int}[] + for v in eachindex(g) + for v1 in g[v] + for v2 in g[v1] + for v3 in g[v2] + if v3 == v + tri = counterclockwise(v1,v2,v3, verts) + push!(tris, tri) + end + end + end + end + end + tris +end + +""" + to_mesh(x,y) +Given line segments with coordinates in x,y of form [x1,x2,NaN,x3,x4,NaN,...] and similar for y, +construct a mesh using vertices and a conductivity matrix. +Returns: +`verts`: list of Tuple{Float64,Float64}s representing points +`tris`: list of unique triangles of the form Tuple(v1,v2,v3), where v1,v2,v3 are indices in counter clockwise order +""" +function to_mesh(x,y) + verts, g = to_graph(x,y) + tris = tris_from_graph(g,verts) + p = reduce(vcat, [v[1] v[2]] for v in verts) + t = reduce(vcat, [t[1] t[2] t[3]] for t in unique(tris)) + p,t +end From 7944b771bd99ba0e162ae97764ed42bd87be5c21 Mon Sep 17 00:00:00 2001 From: Daniel Sharp Date: Wed, 15 Feb 2023 17:02:40 -0500 Subject: [PATCH 2/4] Add docs, speed up mesh --- src/to_mesh.jl | 45 +++++++++++++++++++++++++++++++-------------- 1 file changed, 31 insertions(+), 14 deletions(-) diff --git a/src/to_mesh.jl b/src/to_mesh.jl index 12d8f3d..b0d2731 100644 --- a/src/to_mesh.jl +++ b/src/to_mesh.jl @@ -1,20 +1,23 @@ """ to_graph(x,y) -Assumes that `x,y` are of the form [a,b,NaN,c,d,NaN,....], representing line segments +Assumes that `x,y` are of the form [a,b,NaN,c,d,NaN,....], representing line segments. +Create an adjacency matrix representing the graph of the mesh from the line segments, +and all the vertices """ function to_graph(x,y) - Z = [x y] - Z_pts = Vector{Tuple{Float64,Float64}}(undef, 2size(Z,1)÷3) - for j in 1:size(Z,1)÷3 - Z_pts[2j-1] = (Z[3j-2,1],Z[3j-2,2]) - Z_pts[2j] = (Z[3j-1,1],Z[3j-1,1]) + Z_pts = Vector{Tuple{Float64,Float64}}(undef, 2length(x)÷3) + idx = 1 + @inbounds for j in eachindex(x) + if j % 3 > 0 + Z_pts[idx] = (x[j],y[j]) + idx += 1 + end end - Z_pts = [(z[1],z[2]) for z in eachrow(Z[Not(3:3:end),:])] verts = unique(Z_pts) g=[Int[] for _ in eachindex(verts)] N_edges = size(Z_pts,1)÷2 edges = Vector{Tuple{Int,Int}}(undef, N_edges) - for j in eachindex(edges) + @inbounds for j in eachindex(edges) z1,z2 = Z_pts[[2j-1,2j]] z1_v = findfirst(==(z1), verts) z2_v = findfirst(==(z2), verts) @@ -24,7 +27,11 @@ function to_graph(x,y) verts, g end -function sort3(a,b,c) +""" + sort3(a::T,b::T,c::T) +returns a NTuple{3,T} that gives (a,b,c) from smallest to largest +""" +function sort3(a::T,b::T,c::T)::NTuple{3,T} where {T} t = a if a>b a = b @@ -43,7 +50,11 @@ function sort3(a,b,c) a,b,c end -function counterclockwise(v1::Int,v2::Int,v3::Int, verts) +""" + counterclockwise(v1::Int,v2::Int,v3::Int, verts::Vector{Tuple{T,T}}) +Sorts the indices in ascending order then swaps first two to ensure counterclockwise +""" +function counterclockwise(v1::Int,v2::Int,v3::Int, verts::Vector{Tuple{T,T}}) where {T} v1,v2,v3 = sort3(v1,v2,v3) ax,ay = verts[v1] bx,by = verts[v2] @@ -53,9 +64,15 @@ function counterclockwise(v1::Int,v2::Int,v3::Int, verts) is_cc ? (v1,v2,v3) : (v2,v1,v3) end -function tris_from_graph(g::Vector{Vector{Int}}, verts) - tris = NTuple{3,Int}[] - for v in eachindex(g) +""" + tris_from_graph(g::Vector{Vector{Int}}, verts::Vector{Tuple{T,T}}) +Takes a graph representation of a mesh (represented by a sparse adjacency matrix) and finds all the triangles. +Scales poorly with the number of edges. +""" +function tris_from_graph(g::Vector{Vector{Int}}, verts::Vector{Tuple{T,T}}) where {T} + tris = Set{NTuple{3,Int}}() + sizehint!(tris, 2length(verts)) + @inbounds for v in eachindex(g) for v1 in g[v] for v2 in g[v1] for v3 in g[v2] @@ -82,6 +99,6 @@ function to_mesh(x,y) verts, g = to_graph(x,y) tris = tris_from_graph(g,verts) p = reduce(vcat, [v[1] v[2]] for v in verts) - t = reduce(vcat, [t[1] t[2] t[3]] for t in unique(tris)) + t = reduce(vcat, [t[1] t[2] t[3]] for t in tris) p,t end From 399561d6eb8b9786ff0b475c93786b4b351ffe52 Mon Sep 17 00:00:00 2001 From: Daniel Sharp Date: Wed, 15 Feb 2023 17:06:44 -0500 Subject: [PATCH 3/4] Fix to_mesh docs --- .vscode/settings.json | 3 +++ src/to_mesh.jl | 4 ++-- 2 files changed, 5 insertions(+), 2 deletions(-) create mode 100644 .vscode/settings.json diff --git a/.vscode/settings.json b/.vscode/settings.json new file mode 100644 index 0000000..2b10073 --- /dev/null +++ b/.vscode/settings.json @@ -0,0 +1,3 @@ +{ + "julia.environmentPath": "/home/dannys4/.julia/dev/DistMesh2D" +} \ No newline at end of file diff --git a/src/to_mesh.jl b/src/to_mesh.jl index b0d2731..950f695 100644 --- a/src/to_mesh.jl +++ b/src/to_mesh.jl @@ -92,8 +92,8 @@ end Given line segments with coordinates in x,y of form [x1,x2,NaN,x3,x4,NaN,...] and similar for y, construct a mesh using vertices and a conductivity matrix. Returns: -`verts`: list of Tuple{Float64,Float64}s representing points -`tris`: list of unique triangles of the form Tuple(v1,v2,v3), where v1,v2,v3 are indices in counter clockwise order +`p`: matrix where each row is a point on the mesh +`t`: matrix where each row gives a triangle via the indices of rows of p representing the 3 points """ function to_mesh(x,y) verts, g = to_graph(x,y) From d554522b2332f241b428f45baeec4e2569ea06a9 Mon Sep 17 00:00:00 2001 From: Daniel Sharp Date: Wed, 15 Feb 2023 17:08:03 -0500 Subject: [PATCH 4/4] Finish to_mesh docs --- src/to_mesh.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/to_mesh.jl b/src/to_mesh.jl index 950f695..a6777c7 100644 --- a/src/to_mesh.jl +++ b/src/to_mesh.jl @@ -91,9 +91,10 @@ end to_mesh(x,y) Given line segments with coordinates in x,y of form [x1,x2,NaN,x3,x4,NaN,...] and similar for y, construct a mesh using vertices and a conductivity matrix. + Returns: -`p`: matrix where each row is a point on the mesh -`t`: matrix where each row gives a triangle via the indices of rows of p representing the 3 points +- `p`: matrix where each row is a point on the mesh +- `t`: matrix where each row gives a triangle via the indices of rows of p representing the 3 points """ function to_mesh(x,y) verts, g = to_graph(x,y)