return a DistMeshResults struct

Author: sjkelly

src/distmeshnd.jl

@@ -35,6 +35,13 @@ function distmesh(fdist::Function,
     distmesh(fdist, fh, h, setup, o, w, fp, Val(stats), VT)
 end
 
+
+struct DistMeshResult{PT, TT, STATS}
+    points::Vector{PT}
+    tetrahedra::Vector{TT}
+    stats::STATS
+end
+
 function distmesh(fdist::Function,
                   fh,
                   h::Number,
@@ -71,14 +78,16 @@ function distmesh(fdist::Function,
         facecenteredcubic!(fdist, p, pt_dists, h, setup.iso, origin, widths, VertType)
     end
 
+    # Result struct for holding points, simplices, and iteration statistics
+    result = DistMeshResult(p, GeometryBasics.SimplexFace{4,Int32}[], DistMeshStatistics())
+
     # initialize arrays
     pair_set = Set{Tuple{Int32,Int32}}()        # 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
     L = eltype(VertType)[]                      # vector length of each edge
     L0 = non_uniform ? eltype(VertType)[] : nothing # desired edge length computed by dh (edge length function)
-    t = GeometryBasics.SimplexFace{4,Int32}[]   # tetrahedra indices from delaunay triangulation
     maxmove = typemax(eltype(VertType))         # stores an iteration max movement for retriangulation
 
     # arrays for tracking quality metrics
@@ -87,7 +96,6 @@ function distmesh(fdist::Function,
     qualities = eltype(VertType)[]
 
     # information on each iteration
-    statsdata = DistMeshStatistics()
     lcount = 0 # iteration counter
     triangulationcount = 0 # triangulation counter
 
@@ -102,9 +110,9 @@ function distmesh(fdist::Function,
 
             # compute a new delaunay triangulation
             # we use the default random insertion method
-            delaunayn!(fdist, p, t, geps, false)
+            delaunayn!(fdist, result, geps, false)
 
-            tet_to_edges!(pair, pair_set, t) # Describe each edge by a unique pair of nodes
+            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))
@@ -114,14 +122,14 @@ function distmesh(fdist::Function,
             # if the points were sorted we need to update the distance cache
             if setup.sort && iszero(triangulationcount % setup.sort_interval)
                 for i in eachindex(p)
-                    pt_dists[i] = fdist(p[i])
+                    pt_dists[i] = fdist(result.points[i])
                 end
             end
             triangulationcount += 1
-            stats && push!(statsdata.retriangulations, lcount)
+            stats && push!(result.stats.retriangulations, lcount)
         end
 
-        compute_displacements!(fh, dp, pair, L, L0, bars, p, setup, VertType)
+        compute_displacements!(fh, dp, pair, L, L0, bars, result.points, setup, VertType)
 
         # Zero out forces on fix points
         if have_fixed
@@ -144,7 +152,7 @@ function distmesh(fdist::Function,
             # complex distance functions and large node counts
             move = sqrt(sum((p[i]-p0).^2))          # compute movement from the displacement
             d_est = pt_dists[i] + move              # apply the movement to our cache point
-            d = d_est < -geps ? d_est : fdist(p[i]) # determine if we need correct or approximate distance
+            d = d_est < -geps ? d_est : fdist(result.points[i]) # determine if we need correct or approximate distance
             pt_dists[i] = d                         # store distance
 
             if d < -geps
@@ -167,20 +175,20 @@ function distmesh(fdist::Function,
 
         # save iteration stats
         if stats
-            push!(statsdata.maxmove,maxmove)
-            push!(statsdata.maxdp,maxdp)
-            triangle_qualities!(tris,triset,qualities,p,t)
+            push!(result.stats.maxmove,maxmove)
+            push!(result.stats.maxdp,maxdp)
+            triangle_qualities!(tris,triset,qualities,result.points,result.tetrahedra)
             sort!(qualities) # sort for median calc and robust summation
             mine, maxe = extrema(qualities)
-            push!(statsdata.average_qual, sum(qualities)/length(qualities))
-            push!(statsdata.median_qual, qualities[round(Int,length(qualities)/2)])
-            push!(statsdata.minimum_qual, mine)
-            push!(statsdata.maximum_qual, maxe)
+            push!(result.stats.average_qual, sum(qualities)/length(qualities))
+            push!(result.stats.median_qual, qualities[round(Int,length(qualities)/2)])
+            push!(result.stats.minimum_qual, mine)
+            push!(result.stats.maximum_qual, maxe)
         end
 
         # Termination criterion
         if maxdp<setup.ptol*h
-            return p, t, statsdata
+            return result
         end
     end
 end

src/tetgen.jl

@@ -15,7 +15,9 @@ function delaunayn_nosort(points)
     tetio
 end
 
-function delaunayn!(fdist, p, t, geps, sorted_pts)
+function delaunayn!(fdist, result::DistMeshResult, geps, sorted_pts)
+    t = result.tetrahedra
+    p = result.points
     triangulation = sorted_pts ? delaunayn_nosort(p) : delaunayn(p)
     t_d = triangulation.tetrahedra
     resize!(t, length(t_d))

test/runtests.jl

@@ -65,30 +65,32 @@ end
 
 @testset "distmesh 3D" begin
     d(p) = sqrt(sum(p.^2))-1
-    p,t,_ = distmesh(d,HUniform(),0.2)
-    @test length(p) == 485
-    @test length(t) == 2207
+    result = distmesh(d,HUniform(),0.2)
+    @test length(result.points) == 485
+    @test length(result.tetrahedra) == 2207
 
-    p,t,_ = distmesh(d,HUniform(),0.2, DistMeshSetup(distribution=:packed))
-    @test length(p) == 742
-    @test length(t) == 3472
+    result = distmesh(d,HUniform(),0.2, DistMeshSetup(distribution=:packed))
+    @test length(result.points) == 742
+    @test length(result.tetrahedra) == 3472
 
     # test stats is not messing
-    p,t,s = distmesh(d,HUniform(),0.2, stats=true)
-    @test length(p) == 485
-    @test length(t) == 2207
+    result = distmesh(d,HUniform(),0.2, stats=true)
+    @test length(result.points) == 485
+    @test length(result.tetrahedra) == 2207
 
-    p,t,s = distmesh(d,HUniform(),0.4, stats=true)
-    @test length(p) == 56
-    @test length(t) == 186
-    for fn in fieldnames(typeof(s))
-        @test isapprox(getproperty(s,fn), getproperty(stat_04,fn))
+    result = distmesh(d,HUniform(),0.4, stats=true)
+    @test length(result.points) == 56
+    @test length(result.tetrahedra) == 186
+    for fn in fieldnames(typeof(result.stats))
+        @test isapprox(getproperty(result.stats,fn), getproperty(stat_04,fn))
     end
 end
 
 @testset "dihedral metrics" begin
     d(p) = sqrt(sum(p.^2))-1
-    p,t,_ = distmesh(d,HUniform(),0.2)
+    result = distmesh(d,HUniform(),0.2)
+    p = result.points
+    t = result.tetrahedra
     all_angs =  DistMesh.dihedral_angles(p,t)
     min_angs =  DistMesh.min_dihedral_angles(p,t)
     ax = extrema(all_angs)