Merge pull request #21 from JuliaGeometry/sjk/alloc_fix

Sjk/alloc fix

Author: sjkelly

src/decompositions.jl

@@ -48,21 +48,27 @@ end
 """
 function tet_to_edges!(pair::Vector, pair_set::Set, t)
     empty!(pair_set)
-    for i in eachindex(t)
+    length(pair) < length(t)*6 && resize!(pair, length(t)*6)
+    num_pair = 0
+    @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(elt, pair_set)
+                push!(pair_set, elt)
+                num_pair += 1
+                pair[num_pair] = elt
+            end
         end
     end
-    resize!(pair, length(pair_set))
-    # copy pair set to array since sets are not sortable
-    i = 1
-    for elt in pair_set
-        pair[i] = elt
-        i = i + 1
-    end
 
     # sort the edge pairs for better point lookup
-    sort!(pair)
-end
+    sort!(view(pair, 1:num_pair))
+
+    return num_pair # return the number of pairs
+end
+
+function bitpack(xi,yi)
+    (unsafe_trunc(UInt64, xi) << 32) | unsafe_trunc(UInt64, yi)
+end

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
 
     @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))
+            # resize arrays for new pair count
+            if length(L) < num_pairs
+                resize!(bars, num_pairs)
+                resize!(L, num_pairs)
+                non_uniform && resize!(L0, num_pairs)
+            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
@@ -214,16 +217,17 @@ function retriangulate!(fdist, result::DistMeshResult, geps, setup, triangulatio
 end
 
 
-function compute_displacements!(fh, dp, pair, L, L0, bars, p, setup,
-    ::Type{VertType}) where {VertType}
+
+function compute_displacements!(fh, dp, pair, num_pairs, L, L0, bars, p, setup,
+                                ::Type{VertType}) where {VertType}
 
     non_uniform = isa(typeof(L0), AbstractVector)
 
     # compute edge lengths (L) and adaptive edge lengths (L0)
     # 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)
+    L0sum = non_uniform ? zero(eltype(L0)) : num_pairs
+    for i in 1:num_pairs
         pb = pair[i]
         b1 = p[pb[1]]
         b2 = p[pb[2]]
@@ -245,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
-    for i in eachindex(pair)
+    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