diff --git a/src/decompositions.jl b/src/decompositions.jl
index 4714ff7..e7c4aaa 100644
--- a/src/decompositions.jl
+++ b/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)
+    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]]
             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]))
+                pair[ind] = elt
+                ind += 1
+            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
\ No newline at end of file
+    # 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)
+    (unsafe_trunc(UInt64, xi) << 32) | unsafe_trunc(UInt64, yi)
+end
diff --git a/src/distmeshnd.jl b/src/distmeshnd.jl
index c3fabb2..c9424d9 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
@@ -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
@@ -205,7 +208,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)
@@ -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)
-    for i in eachindex(pair)
+    @inbounds for i in 1:num_pairs
         pb = pair[i]
         b1 = p[pb[1]]
         b2 = p[pb[2]]
@@ -237,7 +240,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 +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)
+    @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