Merge pull request #9 from sjkelly/sjk/sort2

sjkelly · web-flow · commit da4a2143665d · 2019-11-30T16:22:56.000-05:00
Sjk/sort2
diff --git a/bench/hilbert_perf.jl b/bench/hilbert_perf.jl
@@ -0,0 +1,12 @@
+using DistMesh
+using BenchmarkTools
+using StaticArrays
+
+a = [rand(SVector{3,Float64}) for i = 1:50000]
+
+function test_hilbert(a)
+    b = copy(a)
+    DistMesh.hilbertsort!(b)
+end
+
+@benchmark test_hilbert(a)
diff --git a/bench/perf.jl b/bench/perf.jl
@@ -18,7 +18,10 @@ println("Benchmarking DistMesh.jl...")
 # Algorithms to benchmark
 #
 
-algos = [DistMeshSetup(deltat=0.1, distribution=:packed),DistMeshSetup(distribution=:packed)]
+algos = [DistMeshSetup(deltat=0.1, distribution=:packed),
+         DistMeshSetup(distribution=:packed),
+         DistMeshSetup(deltat=0.1, distribution=:packed,sort=true),
+         DistMeshSetup(distribution=:packed,sort=true)]
 
 fn_torus(v) = (sqrt(v[1]^2+v[2]^2)-0.5)^2 + v[3]^2 - 0.25 # torus
 fn_sphere(v) = sqrt(sum(v.^2)) -1
diff --git a/src/DistMesh.jl b/src/DistMesh.jl
@@ -13,27 +13,38 @@ abstract type AbstractDistMeshAlgorithm end
 """
     DistMeshSetup
 
-    iso (default: 0): Value of which to extract the iso surface, inside negative
+    Takes Keyword arguments as follows:
+
+    iso (default: 0): Value of which to extract the isosurface, inside surface is negative
     deltat (default: 0.1): the fraction of edge displacement to apply each iteration
+    sort (default:false): If true and no fixed points, sort points using a hilbert sort.
+    sort_interval (default:20) Use hilbert sort after the specified retriangulations
+    distribution (default: :regular) Initial point distribution, either :regular or :packed.
 """
 struct DistMeshSetup{T} <: AbstractDistMeshAlgorithm
     iso::T
     deltat::T
     ttol::T
     ptol::T
+    sort::Bool # use hilbert sort to cache-localize points
+    sort_interval::Int # retriangulations before resorting
     distribution::Symbol # intial point distribution
 end
 
 function DistMeshSetup(;iso=0,
                         ptol=.001,
                         deltat=0.05,
                         ttol=0.02,
+                        sort=false,
+                        sort_interval=20,
                         distribution=:regular)
     T = promote_type(typeof(iso),typeof(ptol),typeof(deltat), typeof(ttol))
     DistMeshSetup{T}(iso,
                      deltat,
                      ttol,
                      ptol,
+                     sort,
+                     sort_interval,
                      distribution)
 end
 
@@ -92,6 +103,7 @@ include("distmeshnd.jl")
 include("tetgen.jl")
 include("quality.jl")
 include("decompositions.jl")
+include("hilbertsort.jl")
 
 #export distmeshsurface
 export distmesh, DistMeshSetup, DistMeshStatistics, HUniform
diff --git a/src/distmeshnd.jl b/src/distmeshnd.jl
@@ -89,147 +89,20 @@ function distmesh(fdist::Function,
     # 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
 
-            delaunayn!(fdist, p, t, geps) # compute a new delaunay triangulation
-
-            tet_to_edges!(pair, pair_set, t) # 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))
-
-            stats && push!(statsdata.retriangulations, lcount)
-        end
-
-        compute_displacements!(fh, dp, pair, L, L0, bars, p, setup, VertType)
-
-        # Zero out forces on fix points
-        if have_fixed
-            for i in eachindex(fix)
-                dp[i] = zero(VertType)
-            end
-        end
-
-        # apply point forces and
-        # bring outside points back to the boundary
-        maxdp = typemin(eltype(VertType))
-        maxmove = typemin(eltype(VertType))
-        for i in eachindex(p)
-
-            p0 = p[i] # store original point location
-            p[i] = p[i].+setup.deltat.*dp[i] # apply displacements to points
-
-            # Check if we are verifiably within the bounds and use this value
-            # to avoid recomputing fdist. This increases performance greatly on
-            # 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
-            pt_dists[i] = d                         # store distance
-
-            if d < -geps
-                maxdp = max(maxdp, setup.deltat*sqrt(sum(dp[i].^2)))
+            # use hilbert sort to improve cache locality of points
+            if setup.sort && iszero(triangulationcount % setup.sort_interval)
+                hilbertsort!(p)
             end
 
-            if d <= setup.iso
-                maxmove = max(move,maxmove)
-                continue
-            end
-
-            # bring points back to boundary if outside using central difference
-            p[i] = p[i] .- centraldiff(fdist,p[i]).*(d+setup.iso)
-            maxmove = max(sqrt(sum((p[i]-p0).^2)), maxmove)
-            pt_dists[i] = setup.iso # ideally
-        end
-
-        # increment iteration counter
-        lcount = lcount + 1
-
-        # save iteration stats
-        if stats
-            push!(statsdata.maxmove,maxmove)
-            push!(statsdata.maxdp,maxdp)
-            triangle_qualities!(tris,triset,qualities,p,t)
-            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)
-        end
-
-        # Termination criterion
-        if maxdp<setup.ptol*h
-            return p, t, statsdata
-        end
-    end
-end
-
-function distmesh(fdist::Function,
-                  fh,
-                  h::Number,
-                  setup::DistMeshQuality,
-                  origin,
-                  widths,
-                  fix,
-                  ::Val{stats},
-                  ::Type{VertType}) where {VertType, stats}
-
-    geps=1e-1*h+setup.iso # parameter for filtering tets outside bounds and considering for max displacment of a node
-
-    # static parameter info
-    non_uniform = isa(fh, Function) # so we can elide fh calls
-    have_fixed = !isa(fix, Nothing)
-
-    #ptol=.001; ttol=.1; L0mult=1+.4/2^(dim-1); deltat=.2; geps=1e-1*h;
-
-    # initialize Vertex Arrays
-    # the first N points in the array correspond to'fix' points that do not move
-    if have_fixed
-        p = copy(fix)
-    else
-        p = VertType[]
-    end
-
-    pt_dists = map(fdist, p) # cache to store point locations so we can minimize fdist calls
-
-    # add points to p based on the initial distribution
-    if setup.distribution === :regular
-        simplecubic!(fdist, p, pt_dists, h, setup.iso, origin, widths, VertType)
-    elseif setup.distribution === :packed
-        # face-centered cubic point distribution
-        facecenteredcubic!(fdist, p, pt_dists, h, setup.iso, origin, widths, VertType)
-    end
-
-    # 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
-    non_uniform && (L0 = eltype(VertType)[])    # 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
-    tris = NTuple{3,Int32}[]        # triangles used for quality checks
-    triset = Set{NTuple{3,Int32}}() # set for triangles to ensure uniqueness
-    qualities = eltype(VertType)[]
-
-    # information on each iteration
-    statsdata = DistMeshStatistics()
-    lcount = 0 # iteration counter
-
-    @inbounds while true
-        # if large move, retriangulation
-        if maxmove>setup.ttol*h
-
-            delaunayn!(fdist, p, t, geps) # compute a new delaunay triangulation
+            # compute a new delaunay triangulation
+            # we use the default random insertion method
+            delaunayn!(fdist, p, t, geps, false)
 
             tet_to_edges!(pair, pair_set, t) # Describe each edge by a unique pair of nodes
 
@@ -238,6 +111,13 @@ function distmesh(fdist::Function,
             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)
         end
 
diff --git a/src/hilbertsort.jl b/src/hilbertsort.jl
@@ -0,0 +1,135 @@
+# implementing scale-free Hilbert ordering. Real all about it here:
+# http://doc.cgal.org/latest/Spatial_sorting/index.html
+
+# original implementation in:
+# https://github.com/JuliaGeometry/GeometricalPredicates.jl
+# The GeometricalPredicates.jl package is licensed under the MIT "Expat" License:
+#    Copyright (c) 2014: Ariel Keselman.
+
+# modifications for StaticArrays: sjkelly
+
+const coordinatex = 1
+const coordinatey = 2
+const coordinatez = 3
+next2d(c) = c % 2 + 1
+next3d(c) = c % 3 + 1
+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
+    #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]
+                end
+                return #v[k]
+            end
+            pivot = v[(lo+hi)>>>1]
+            i, j = lo, hi
+            while true
+                pivot_elt = pivot[coord]
+                while v[i][coord] < pivot_elt; i += 1; end
+                while pivot_elt < v[j][coord]; j -= 1; end
+                i <= j || break
+                v[i], v[j] = v[j], v[i]
+                i += 1; j -= 1
+            end
+            if k <= j
+                hi = j
+            elseif i <= k
+                lo = i
+            else
+                return #pivot
+            end
+        end
+    else
+        @inbounds while lo < hi
+            if isone(hi-lo)
+                if v[hi][coord] > v[lo][coord]
+                    v[lo], v[hi] = v[hi], v[lo]
+                end
+                return #v[k]
+            end
+            pivot = v[(lo+hi)>>>1]
+            i, j = lo, hi
+            while true
+                pivot_elt = pivot[coord]
+                while v[i][coord] > pivot_elt; i += 1; end
+                while pivot_elt > v[j][coord]; j -= 1; end
+                i <= j || break
+                v[i], v[j] = v[j], v[i]
+                i += 1; j -= 1
+            end
+            if k <= j
+                hi = j
+            elseif i <= k
+                lo = i
+            else
+                return #pivot
+            end
+        end
+    end
+    #return v[lo]
+    nothing
+end
+
+
+# 2D version
+# function hilbertsort!(directionx::AbstractDirection, directiony::AbstractDirection, coordinate::AbstractCoordinate, a::Array{T,1}, lo::Int64, hi::Int64, lim::Int64=4) where T<:AbstractPoint2D
+#     hi-lo <= lim && return a
+
+#     i2 = (lo+hi)>>>1
+#     i1 = (lo+i2)>>>1
+#     i3 = (i2+hi)>>>1
+
+#     select!(directionx, coordinate, a, i2, lo, hi)
+#     select!(directiony, next2d(coordinate), a, i1, lo, i2)
+#     select!(!directiony, next2d(coordinate), a, i3, i2, hi)
+
+#     hilbertsort!(directiony, directionx, next2d(coordinate), a, lo, i1, lim)
+#     hilbertsort!(directionx, directiony, coordinate, a, i1, i2, lim)
+#     hilbertsort!(directionx, directiony, coordinate, a, i2, i3, lim)
+#     hilbertsort!(!directiony, !directionx, next2d(coordinate), a, i3, hi, lim)
+
+#     return a
+# end
+
+function hilbertsort!(directionx, directiony, directionz, coordinate, a::Vector, lo::Integer, hi::Integer, lim::Integer=8)
+    hi-lo <= lim && return a
+
+    i4 = (lo+hi)>>>1
+    i2 = (lo+i4)>>>1
+    i1 = (lo+i2)>>>1
+    i3 = (i2+i4)>>>1
+    i6 = (i4+hi)>>>1
+    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)
+
+    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))
+hilbertsort!(a::Vector, lo::Int64, hi::Int64, lim::Int64) = hilbertsort!(backward, backward, backward, coordinatey, a, lo, hi, lim)
diff --git a/src/tetgen.jl b/src/tetgen.jl
@@ -3,8 +3,13 @@ function delaunayn(points)
     tetio
 end
 
-function delaunayn!(fdist, p, t, geps)
-    triangulation = delaunayn(p)
+function delaunayn_nosort(points)
+    tetio = tetrahedralize(TetGen.TetgenIO(points), "Qb/1") # Q- Quiet
+    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
diff --git a/test/runtests.jl b/test/runtests.jl
diff --git a/test/vals.jl b/test/vals.jl
