From efbe7e995b802727d0a0d40772eb4fc4e8749f98 Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Sun, 16 Feb 2020 11:02:31 -0500 Subject: [PATCH 1/3] add inbounds in some spots --- src/decompositions.jl | 6 +++--- src/distmeshnd.jl | 8 ++++---- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/decompositions.jl b/src/decompositions.jl index 4714ff7..c93585e 100644 --- a/src/decompositions.jl +++ b/src/decompositions.jl @@ -48,7 +48,7 @@ end """ function tet_to_edges!(pair::Vector, pair_set::Set, t) empty!(pair_set) - for i in eachindex(t) + @inbounds for i in eachindex(t) for ep in 1:6 p1 = t[i][tetpairs[ep][1]] p2 = t[i][tetpairs[ep][2]] @@ -58,11 +58,11 @@ function tet_to_edges!(pair::Vector, pair_set::Set, t) resize!(pair, length(pair_set)) # copy pair set to array since sets are not sortable i = 1 - for elt in pair_set + @inbounds for elt in pair_set pair[i] = elt i = i + 1 end # sort the edge pairs for better point lookup sort!(pair) -end \ No newline at end of file +end diff --git a/src/distmeshnd.jl b/src/distmeshnd.jl index c3fabb2..e957aee 100644 --- a/src/distmeshnd.jl +++ b/src/distmeshnd.jl @@ -205,7 +205,7 @@ function retriangulate!(fdist, result::DistMeshResult, geps, setup, triangulatio # TODO: can we use the point distance array to pass boundary points to # tetgen so this call is no longer required? j = firstindex(t) - for ai in t + @inbounds for ai in t t[j] = ai pm = (p[ai[1]].+p[ai[2]].+p[ai[3]].+p[ai[4]])./4 j = ifelse(fdist(pm) <= -geps, nextind(t, j), j) @@ -224,7 +224,7 @@ function compute_displacements!(fh, dp, pair, L, L0, bars, p, setup, # Lp norm (p=3) is partially computed here Lsum = zero(eltype(L)) L0sum = non_uniform ? zero(eltype(L0)) : length(pair) - for i in eachindex(pair) + @inbounds for i in eachindex(pair) pb = pair[i] b1 = p[pb[1]] b2 = p[pb[2]] @@ -237,7 +237,7 @@ function compute_displacements!(fh, dp, pair, L, L0, bars, p, setup, end # zero out force at each node - for i in eachindex(dp) + @inbounds for i in eachindex(dp) dp[i] = zero(VertType) end @@ -246,7 +246,7 @@ function compute_displacements!(fh, dp, pair, L, L0, bars, p, setup, lscbrt = (1+(0.4/2^2))*cbrt(Lsum/L0sum) # Move mesh points based on edge lengths L and forces F - for i in eachindex(pair) + @inbounds for i in eachindex(pair) if non_uniform && L[i] < L0[i]*lscbrt || L[i] < lscbrt L0_f = non_uniform ? L0[i].*lscbrt : lscbrt # compute force vectors From df83e1486615977c51e4ecd8eecaa911976d66da Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Sun, 16 Feb 2020 11:51:33 -0500 Subject: [PATCH 2/3] some performance improvements, reduce worst-case performance --- src/decompositions.jl | 23 +++++++++++++---------- src/distmeshnd.jl | 2 +- 2 files changed, 14 insertions(+), 11 deletions(-) diff --git a/src/decompositions.jl b/src/decompositions.jl index c93585e..67f260c 100644 --- a/src/decompositions.jl +++ b/src/decompositions.jl @@ -48,21 +48,24 @@ end """ function tet_to_edges!(pair::Vector, pair_set::Set, t) empty!(pair_set) + empty!(pair) @inbounds for i in eachindex(t) for ep in 1:6 p1 = t[i][tetpairs[ep][1]] p2 = t[i][tetpairs[ep][2]] - push!(pair_set, p1 > p2 ? (p2,p1) : (p1,p2)) + elt = p1 > p2 ? (p2,p1) : (p1,p2) + if !in(bitpack(elt[1],elt[2]), pair_set) + push!(pair_set, bitpack(elt[1],elt[2])) + push!(pair, elt) + end end end - resize!(pair, length(pair_set)) - # copy pair set to array since sets are not sortable - i = 1 - @inbounds for elt in pair_set - pair[i] = elt - i = i + 1 - end - # sort the edge pairs for better point lookup - sort!(pair) + # TODO This seems to be really good for some geoemetries, + # but intoduces larger performance regressions for others + #sort!(pair) +end + +function bitpack(xi,yi) + (unsafe_trunc(UInt64, xi) << 32) | unsafe_trunc(UInt64, yi) end diff --git a/src/distmeshnd.jl b/src/distmeshnd.jl index e957aee..bfcade0 100644 --- a/src/distmeshnd.jl +++ b/src/distmeshnd.jl @@ -84,7 +84,7 @@ function distmesh(fdist::Function, stats ? DistMeshStatistics() : nothing) # initialize arrays - pair_set = Set{Tuple{Int32,Int32}}() # set used for ensure we have a unique set of edges + pair_set = Set{UInt64}() # set used for ensure we have a unique set of edges pair = Tuple{Int32,Int32}[] # edge indices (Int32 since we use Tetgen) dp = zeros(VertType, length(p)) # displacement at each node bars = VertType[] # the vector of each edge From 3b4f567a7b0a0aa34bf76b8f1b7c823e7d14613f Mon Sep 17 00:00:00 2001 From: Steve Kelly Date: Sun, 16 Feb 2020 15:28:45 -0500 Subject: [PATCH 3/3] initial sketch of a better allocation strategy --- src/decompositions.jl | 7 +++++-- src/distmeshnd.jl | 19 +++++++++++-------- 2 files changed, 16 insertions(+), 10 deletions(-) diff --git a/src/decompositions.jl b/src/decompositions.jl index 67f260c..e7c4aaa 100644 --- a/src/decompositions.jl +++ b/src/decompositions.jl @@ -48,7 +48,8 @@ end """ function tet_to_edges!(pair::Vector, pair_set::Set, t) empty!(pair_set) - empty!(pair) + iszero(length(pair)) && resize!(pair, length(t)*6) # allocate memory on first iteration + ind = 1 @inbounds for i in eachindex(t) for ep in 1:6 p1 = t[i][tetpairs[ep][1]] @@ -56,7 +57,8 @@ function tet_to_edges!(pair::Vector, pair_set::Set, t) elt = p1 > p2 ? (p2,p1) : (p1,p2) if !in(bitpack(elt[1],elt[2]), pair_set) push!(pair_set, bitpack(elt[1],elt[2])) - push!(pair, elt) + pair[ind] = elt + ind += 1 end end end @@ -64,6 +66,7 @@ function tet_to_edges!(pair::Vector, pair_set::Set, t) # TODO This seems to be really good for some geoemetries, # but intoduces larger performance regressions for others #sort!(pair) + return ind-1 end function bitpack(xi,yi) diff --git a/src/distmeshnd.jl b/src/distmeshnd.jl index bfcade0..c9424d9 100644 --- a/src/distmeshnd.jl +++ b/src/distmeshnd.jl @@ -95,6 +95,7 @@ function distmesh(fdist::Function, # information on each iteration lcount = 0 # iteration counter triangulationcount = 0 # triangulation counter + num_pairs = 0 # number of edge pairs @inbounds while true # if large move, retriangulation @@ -103,18 +104,20 @@ function distmesh(fdist::Function, # compute a new delaunay triangulation retriangulate!(fdist, result, geps, setup, triangulationcount, pt_dists) - tet_to_edges!(pair, pair_set, result.tetrahedra) # Describe each edge by a unique pair of nodes + num_pairs = tet_to_edges!(pair, pair_set, result.tetrahedra) # Describe each edge by a unique pair of nodes # resize arrays for new pair counts - resize!(bars, length(pair)) - resize!(L, length(pair)) - non_uniform && resize!(L0, length(pair)) + if triangulationcount == 0 + resize!(bars, length(result.tetrahedra)*6) + resize!(L, length(result.tetrahedra)*6) + non_uniform && resize!(L0, length(result.tetrahedra)*6) + end triangulationcount += 1 stats && push!(result.stats.retriangulations, lcount) end - compute_displacements!(fh, dp, pair, L, L0, bars, result.points, setup, VertType) + compute_displacements!(fh, dp, pair, num_pairs, L, L0, bars, result.points, setup, VertType) # Zero out forces on fix points if have_fixed @@ -215,7 +218,7 @@ function retriangulate!(fdist, result::DistMeshResult, geps, setup, triangulatio end -function compute_displacements!(fh, dp, pair, L, L0, bars, p, setup, +function compute_displacements!(fh, dp, pair, num_pairs, L, L0, bars, p, setup, ::Type{VertType}) where {VertType} non_uniform = isa(typeof(L0), AbstractVector) @@ -224,7 +227,7 @@ function compute_displacements!(fh, dp, pair, L, L0, bars, p, setup, # Lp norm (p=3) is partially computed here Lsum = zero(eltype(L)) L0sum = non_uniform ? zero(eltype(L0)) : length(pair) - @inbounds for i in eachindex(pair) + @inbounds for i in 1:num_pairs pb = pair[i] b1 = p[pb[1]] b2 = p[pb[2]] @@ -246,7 +249,7 @@ function compute_displacements!(fh, dp, pair, L, L0, bars, p, setup, lscbrt = (1+(0.4/2^2))*cbrt(Lsum/L0sum) # Move mesh points based on edge lengths L and forces F - @inbounds for i in eachindex(pair) + @inbounds for i in 1:num_pairs if non_uniform && L[i] < L0[i]*lscbrt || L[i] < lscbrt L0_f = non_uniform ? L0[i].*lscbrt : lscbrt # compute force vectors