Merge pull request #15 from JuliaGeometry/sjk/refactor1

Sjk/refactor1

Author: sjkelly

src/distmeshnd.jl

@@ -35,6 +35,18 @@ function distmesh(fdist::Function,
     distmesh(fdist, fh, h, setup, o, w, fp, Val(stats), VT)
 end
 
+"""
+    DistMeshResult
+
+A struct returned from the `distmesh` function that includes point, simplex,
+and interation statistics.
+"""
+struct DistMeshResult{PT, TT, STATS}
+    points::Vector{PT}
+    tetrahedra::Vector{TT}
+    stats::STATS
+end
+
 function distmesh(fdist::Function,
                   fh,
                   h::Number,
@@ -71,14 +83,18 @@ 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}[],
+                            stats ? DistMeshStatistics() : nothing)
+
     # 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,41 +103,28 @@ function distmesh(fdist::Function,
     qualities = eltype(VertType)[]
 
     # information on each iteration
-    statsdata = DistMeshStatistics()
     lcount = 0 # iteration counter
     triangulationcount = 0 # triangulation counter
 
     @inbounds while true
         # if large move, retriangulation
         if maxmove>setup.ttol*h
 
-            # use hilbert sort to improve cache locality of points
-            if setup.sort && iszero(triangulationcount % setup.sort_interval)
-                hilbertsort!(p)
-            end
-
             # compute a new delaunay triangulation
-            # we use the default random insertion method
-            delaunayn!(fdist, p, t, geps, false)
+            retriangulate!(fdist, result, geps, setup, triangulationcount, pt_dists)
 
-            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))
             resize!(L, length(pair))
             non_uniform && resize!(L0, length(pair))
 
-            # 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])
-                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 +147,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 +170,59 @@ 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
+
+"""
+    retriangulate!
+
+
+Given a point set, generate a delaunay triangulation, and other requirements.
+This includes:
+    - Spatial sorting of points
+    - Delaunay triangulation
+    - Filtering of invalid tetrahedra outside the boundary
+"""
+function retriangulate!(fdist, result::DistMeshResult, geps, setup, triangulationcount, pt_dists)
+    # use hilbert sort to improve cache locality of points
+    if setup.sort && iszero(triangulationcount % setup.sort_interval)
+        hilbertsort!(result.points, pt_dists)
+    end
+
+    t = result.tetrahedra
+    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
+
+    # 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
+        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)
+    end
+    j <= lastindex(t) && resize!(t, j-1)
+    nothing
+end

src/hilbertsort.jl

@@ -20,13 +20,14 @@ nextnext3d(c) = (c + 1) % 3 + 1
 const forward = true
 const backward = false
 
-function select!(direction, coord, v::Array{T,1}, k::Integer, lo::Integer, hi::Integer) where T<:AbstractVector
+function select!(direction, coord, v::Array{T,1}, k::Integer, lo::Integer, hi::Integer, carry::CT) where {T<:AbstractVector, CT}
     #lo <= k <= hi || error("select index $k is out of range $lo:$hi")
     if direction == forward
         @inbounds while lo < hi
             if isone(hi-lo)
                 if v[hi][coord] < v[lo][coord]
                     v[lo], v[hi] = v[hi], v[lo]
+                    if CT !== Nothing; carry[lo], carry[hi] = carry[hi], carry[lo]; end
                 end
                 return #v[k]
             end
@@ -38,6 +39,7 @@ function select!(direction, coord, v::Array{T,1}, k::Integer, lo::Integer, hi::I
                 while pivot_elt < v[j][coord]; j -= 1; end
                 i <= j || break
                 v[i], v[j] = v[j], v[i]
+                if CT !== Nothing; carry[i], carry[j] = carry[j], carry[i]; end
                 i += 1; j -= 1
             end
             if k <= j
@@ -53,6 +55,7 @@ function select!(direction, coord, v::Array{T,1}, k::Integer, lo::Integer, hi::I
             if isone(hi-lo)
                 if v[hi][coord] > v[lo][coord]
                     v[lo], v[hi] = v[hi], v[lo]
+                    if CT !== Nothing; carry[lo], carry[hi] = carry[hi], carry[lo]; end
                 end
                 return #v[k]
             end
@@ -64,6 +67,7 @@ function select!(direction, coord, v::Array{T,1}, k::Integer, lo::Integer, hi::I
                 while pivot_elt > v[j][coord]; j -= 1; end
                 i <= j || break
                 v[i], v[j] = v[j], v[i]
+                if CT !== Nothing; carry[i], carry[j] = carry[j], carry[i]; end
                 i += 1; j -= 1
             end
             if k <= j
@@ -100,7 +104,7 @@ end
 #     return a
 # end
 
-function hilbertsort!(directionx, directiony, directionz, coordinate, a::Vector, lo::Integer, hi::Integer, lim::Integer=8)
+function hilbertsort!(directionx, directiony, directionz, coordinate, a::Vector, lo::Integer, hi::Integer, lim::Integer, carry)
     hi-lo <= lim && return a
 
     i4 = (lo+hi)>>>1
@@ -111,27 +115,31 @@ function hilbertsort!(directionx, directiony, directionz, coordinate, a::Vector,
     i5 = (i4+i6)>>>1
     i7 = (i6+hi)>>>1
 
-    select!(directionx, coordinate, a, i4, lo, hi)
-    select!(directiony, next3d(coordinate), a, i2, lo, i4)
-    select!(directionz, nextnext3d(coordinate), a, i1, lo, i2)
-    select!(!directionz, nextnext3d(coordinate), a, i3, i2, i4)
-    select!(!directiony, next3d(coordinate), a, i6, i4, hi)
-    select!(directionz, nextnext3d(coordinate), a, i5, i4, i6)
-    select!(!directionz, nextnext3d(coordinate), a, i7, i6, hi)
-
-    hilbertsort!( directionz,  directionx,  directiony, nextnext3d(coordinate), a, lo, i1, lim)
-    hilbertsort!( directiony,  directionz,  directionx, next3d(coordinate),     a, i1, i2, lim)
-    hilbertsort!( directiony,  directionz,  directionx, next3d(coordinate),     a, i2, i3, lim)
-    hilbertsort!( directionx, !directiony, !directionz, coordinate,             a, i3, i4, lim)
-    hilbertsort!( directionx, !directiony, !directionz, coordinate,             a, i4, i5, lim)
-    hilbertsort!(!directiony,  directionz, !directionx, next3d(coordinate),     a, i5, i6, lim)
-    hilbertsort!(!directiony,  directionz, !directionx, next3d(coordinate),     a, i6, i7, lim)
-    hilbertsort!(!directionz, !directionx,  directiony, nextnext3d(coordinate), a, i7, hi, lim)
+    select!(directionx, coordinate, a, i4, lo, hi, carry)
+    select!(directiony, next3d(coordinate), a, i2, lo, i4, carry)
+    select!(directionz, nextnext3d(coordinate), a, i1, lo, i2, carry)
+    select!(!directionz, nextnext3d(coordinate), a, i3, i2, i4, carry)
+    select!(!directiony, next3d(coordinate), a, i6, i4, hi, carry)
+    select!(directionz, nextnext3d(coordinate), a, i5, i4, i6, carry)
+    select!(!directionz, nextnext3d(coordinate), a, i7, i6, hi, carry)
+
+    hilbertsort!( directionz,  directionx,  directiony, nextnext3d(coordinate), a, lo, i1, lim, carry)
+    hilbertsort!( directiony,  directionz,  directionx, next3d(coordinate),     a, i1, i2, lim, carry)
+    hilbertsort!( directiony,  directionz,  directionx, next3d(coordinate),     a, i2, i3, lim, carry)
+    hilbertsort!( directionx, !directiony, !directionz, coordinate,             a, i3, i4, lim, carry)
+    hilbertsort!( directionx, !directiony, !directionz, coordinate,             a, i4, i5, lim, carry)
+    hilbertsort!(!directiony,  directionz, !directionx, next3d(coordinate),     a, i5, i6, lim, carry)
+    hilbertsort!(!directiony,  directionz, !directionx, next3d(coordinate),     a, i6, i7, lim, carry)
+    hilbertsort!(!directionz, !directionx,  directiony, nextnext3d(coordinate), a, i7, hi, lim, carry)
 
     return a
 end
 
 #hilbertsort!(a::Array{T,1}) where {T<:AbstractPoint2D} = hilbertsort!(backward, backward, coordinatey, a, 1, length(a))
 #hilbertsort!(a::Array{T,1}, lo::Int64, hi::Int64, lim::Int64) where {T<:AbstractPoint2D} = hilbertsort!(backward, backward, coordinatey, a, lo, hi, lim)
-hilbertsort!(a::Vector) = hilbertsort!(backward, backward, backward, coordinatez, a, 1, length(a))
+"""
+Hilbert Sorting. If `carry` is specified, this array will be permuted in line with the
+specified array.
+"""
+hilbertsort!(a::Vector, carry=nothing) = hilbertsort!(backward, backward, backward, coordinatez, a, 1, length(a), 8, carry)
 hilbertsort!(a::Vector, lo::Int64, hi::Int64, lim::Int64) = hilbertsort!(backward, backward, backward, coordinatey, a, lo, hi, lim)

src/tetgen.jl

@@ -15,28 +15,6 @@ function delaunayn_nosort(points)
     tetio
 end
 
-function delaunayn!(fdist, p, t, geps, sorted_pts)
-    triangulation = sorted_pts ? delaunayn_nosort(p) : 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
-
-    # 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
-        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)
-    end
-    j <= lastindex(t) && resize!(t, j-1)
-    nothing
-end
-
 # needs some tweaks, gives garbage results, might need to twek the julia wrapper?
 function reconstruct!(points, tets)
     tetio = tetrahedralize(TetGen.TetgenIO(points,tetrahedrons=tets), "Qr") # Q- Quiet, r- retriangulate

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)