wip

Author: sjkelly

src/decompositions.jl

@@ -46,11 +46,11 @@ end
     Decompose tets to edges, using a pre-allocated array and set.
     Set ensures uniqueness, and result will be sorted.
 """
-function tet_to_edges!(pair::Vector, pair_set::Set, t)
+function tet_to_edges!(pair::Vector, pair_set::Set, t, num_tets)
     empty!(pair_set)
-    length(pair) < length(t)*6 && resize!(pair, length(t)*6)
+    length(pair) < num_tets*6 && resize!(pair, num_tets*6)
     num_pair = 1
-    @inbounds for i in eachindex(t)
+    @inbounds for i in 1:num_tets
         for ep in 1:6
             p1 = t[i][tetpairs[ep][1]]
             p2 = t[i][tetpairs[ep][2]]

src/distmeshnd.jl

@@ -96,15 +96,16 @@ function distmesh(fdist::Function,
     lcount = 0 # iteration counter
     triangulationcount = 0 # triangulation counter
     num_pairs = 0
+    num_tets = 0
 
     @inbounds while true
         # if large move, retriangulation
         if maxmove>setup.ttol*h
 
             # compute a new delaunay triangulation
-            retriangulate!(fdist, result, geps, setup, triangulationcount, pt_dists)
+            num_tets = retriangulate!(fdist, result, geps, setup, triangulationcount, pt_dists)
 
-            num_pairs = 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, num_tets) # Describe each edge by a unique pair of nodes
 
             # resize arrays for new pair count
             length(bars) < num_pairs && resize!(bars, num_pairs)
@@ -170,6 +171,7 @@ function distmesh(fdist::Function,
 
         # Termination criterion
         if maxdp<setup.ptol*h
+            resize!(result.tetrahedra, num_tets)
             return result
         end
     end
@@ -195,23 +197,23 @@ function retriangulate!(fdist, result::DistMeshResult, geps, setup, triangulatio
     p = result.points
     triangulation = delaunayn(p)
     t_d = triangulation.tetrahedra
-    resize!(t, length(t_d))
-    copyto!(t, t_d) # we need to copy since we have a shared reference with tetgen
+    length(t) < length(t_d) && resize!(t, length(t_d))
 
     # average points to get mid point of each tetrahedra
     # if the mid point of the tetrahedra is outside of
     # the boundary we remove it.
     # TODO: this is an inlined filter call. Would be good to revert
     # 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
+    j = 1
+    for ai in t_d
         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)
+        j = ifelse(fdist(pm) <= -geps, j+1, j)
     end
-    j <= lastindex(t) && resize!(t, j-1)
-    nothing
+    #j <= lastindex(t) && resize!(t, j-1)
+    j <= length(t_d) && return j-1
+    return length(t_d)
 end