Skip to content

Commit e130de1

Browse files
sjkellysjkelly
committed
sort both points and distance cache simultaneously to avoid resamples
1 parent 3e38796 commit e130de1

File tree

2 files changed

+26
-28
lines changed

2 files changed

+26
-28
lines changed

‎src/distmeshnd.jl

Lines changed: 3 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -106,7 +106,7 @@ function distmesh(fdist::Function,
106106
if maxmove>setup.ttol*h
107107

108108
# compute a new delaunay triangulation
109-
retriangulate!(fdist, result, geps, setup, triangulationcount)
109+
retriangulate!(fdist, result, geps, setup, triangulationcount, pt_dists)
110110

111111
tet_to_edges!(pair, pair_set, result.tetrahedra) # Describe each edge by a unique pair of nodes
112112

@@ -115,12 +115,6 @@ function distmesh(fdist::Function,
115115
resize!(L, length(pair))
116116
non_uniform && resize!(L0, length(pair))
117117

118-
# if the points were sorted we need to update the distance cache
119-
if setup.sort && iszero(triangulationcount % setup.sort_interval)
120-
for i in eachindex(p)
121-
pt_dists[i] = fdist(result.points[i])
122-
end
123-
end
124118
triangulationcount += 1
125119
stats && push!(result.stats.retriangulations, lcount)
126120
end
@@ -189,10 +183,10 @@ function distmesh(fdist::Function,
189183
end
190184
end
191185

192-
function retriangulate!(fdist, result::DistMeshResult, geps, setup, triangulationcount)
186+
function retriangulate!(fdist, result::DistMeshResult, geps, setup, triangulationcount, pt_dists)
193187
# use hilbert sort to improve cache locality of points
194188
if setup.sort && iszero(triangulationcount % setup.sort_interval)
195-
hilbertsort!(result.points)
189+
hilbertsort!(result.points, pt_dists)
196190
end
197191

198192
t = result.tetrahedra

‎src/hilbertsort.jl

Lines changed: 23 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -20,13 +20,14 @@ nextnext3d(c) = (c + 1) % 3 + 1
2020
const forward = true
2121
const backward = false
2222

23-
function select!(direction, coord, v::Array{T,1}, k::Integer, lo::Integer, hi::Integer) where T<:AbstractVector
23+
function select!(direction, coord, v::Array{T,1}, k::Integer, lo::Integer, hi::Integer, carry::CT) where {T<:AbstractVector, CT}
2424
#lo <= k <= hi || error("select index $k is out of range $lo:$hi")
2525
if direction == forward
2626
@inbounds while lo < hi
2727
if isone(hi-lo)
2828
if v[hi][coord] < v[lo][coord]
2929
v[lo], v[hi] = v[hi], v[lo]
30+
if CT !== Nothing; carry[lo], carry[hi] = carry[hi], carry[lo]; end
3031
end
3132
return #v[k]
3233
end
@@ -38,6 +39,7 @@ function select!(direction, coord, v::Array{T,1}, k::Integer, lo::Integer, hi::I
3839
while pivot_elt < v[j][coord]; j -= 1; end
3940
i <= j || break
4041
v[i], v[j] = v[j], v[i]
42+
if CT !== Nothing; carry[i], carry[j] = carry[j], carry[i]; end
4143
i += 1; j -= 1
4244
end
4345
if k <= j
@@ -53,6 +55,7 @@ function select!(direction, coord, v::Array{T,1}, k::Integer, lo::Integer, hi::I
5355
if isone(hi-lo)
5456
if v[hi][coord] > v[lo][coord]
5557
v[lo], v[hi] = v[hi], v[lo]
58+
if CT !== Nothing; carry[lo], carry[hi] = carry[hi], carry[lo]; end
5659
end
5760
return #v[k]
5861
end
@@ -64,6 +67,7 @@ function select!(direction, coord, v::Array{T,1}, k::Integer, lo::Integer, hi::I
6467
while pivot_elt > v[j][coord]; j -= 1; end
6568
i <= j || break
6669
v[i], v[j] = v[j], v[i]
70+
if CT !== Nothing; carry[i], carry[j] = carry[j], carry[i]; end
6771
i += 1; j -= 1
6872
end
6973
if k <= j
@@ -100,7 +104,7 @@ end
100104
# return a
101105
# end
102106

103-
function hilbertsort!(directionx, directiony, directionz, coordinate, a::Vector, lo::Integer, hi::Integer, lim::Integer=8)
107+
function hilbertsort!(directionx, directiony, directionz, coordinate, a::Vector, lo::Integer, hi::Integer, lim::Integer, carry)
104108
hi-lo <= lim && return a
105109

106110
i4 = (lo+hi)>>>1
@@ -111,27 +115,27 @@ function hilbertsort!(directionx, directiony, directionz, coordinate, a::Vector,
111115
i5 = (i4+i6)>>>1
112116
i7 = (i6+hi)>>>1
113117

114-
select!(directionx, coordinate, a, i4, lo, hi)
115-
select!(directiony, next3d(coordinate), a, i2, lo, i4)
116-
select!(directionz, nextnext3d(coordinate), a, i1, lo, i2)
117-
select!(!directionz, nextnext3d(coordinate), a, i3, i2, i4)
118-
select!(!directiony, next3d(coordinate), a, i6, i4, hi)
119-
select!(directionz, nextnext3d(coordinate), a, i5, i4, i6)
120-
select!(!directionz, nextnext3d(coordinate), a, i7, i6, hi)
121-
122-
hilbertsort!( directionz, directionx, directiony, nextnext3d(coordinate), a, lo, i1, lim)
123-
hilbertsort!( directiony, directionz, directionx, next3d(coordinate), a, i1, i2, lim)
124-
hilbertsort!( directiony, directionz, directionx, next3d(coordinate), a, i2, i3, lim)
125-
hilbertsort!( directionx, !directiony, !directionz, coordinate, a, i3, i4, lim)
126-
hilbertsort!( directionx, !directiony, !directionz, coordinate, a, i4, i5, lim)
127-
hilbertsort!(!directiony, directionz, !directionx, next3d(coordinate), a, i5, i6, lim)
128-
hilbertsort!(!directiony, directionz, !directionx, next3d(coordinate), a, i6, i7, lim)
129-
hilbertsort!(!directionz, !directionx, directiony, nextnext3d(coordinate), a, i7, hi, lim)
118+
select!(directionx, coordinate, a, i4, lo, hi, carry)
119+
select!(directiony, next3d(coordinate), a, i2, lo, i4, carry)
120+
select!(directionz, nextnext3d(coordinate), a, i1, lo, i2, carry)
121+
select!(!directionz, nextnext3d(coordinate), a, i3, i2, i4, carry)
122+
select!(!directiony, next3d(coordinate), a, i6, i4, hi, carry)
123+
select!(directionz, nextnext3d(coordinate), a, i5, i4, i6, carry)
124+
select!(!directionz, nextnext3d(coordinate), a, i7, i6, hi, carry)
125+
126+
hilbertsort!( directionz, directionx, directiony, nextnext3d(coordinate), a, lo, i1, lim, carry)
127+
hilbertsort!( directiony, directionz, directionx, next3d(coordinate), a, i1, i2, lim, carry)
128+
hilbertsort!( directiony, directionz, directionx, next3d(coordinate), a, i2, i3, lim, carry)
129+
hilbertsort!( directionx, !directiony, !directionz, coordinate, a, i3, i4, lim, carry)
130+
hilbertsort!( directionx, !directiony, !directionz, coordinate, a, i4, i5, lim, carry)
131+
hilbertsort!(!directiony, directionz, !directionx, next3d(coordinate), a, i5, i6, lim, carry)
132+
hilbertsort!(!directiony, directionz, !directionx, next3d(coordinate), a, i6, i7, lim, carry)
133+
hilbertsort!(!directionz, !directionx, directiony, nextnext3d(coordinate), a, i7, hi, lim, carry)
130134

131135
return a
132136
end
133137

134138
#hilbertsort!(a::Array{T,1}) where {T<:AbstractPoint2D} = hilbertsort!(backward, backward, coordinatey, a, 1, length(a))
135139
#hilbertsort!(a::Array{T,1}, lo::Int64, hi::Int64, lim::Int64) where {T<:AbstractPoint2D} = hilbertsort!(backward, backward, coordinatey, a, lo, hi, lim)
136-
hilbertsort!(a::Vector) = hilbertsort!(backward, backward, backward, coordinatez, a, 1, length(a))
140+
hilbertsort!(a::Vector, carry=nothing) = hilbertsort!(backward, backward, backward, coordinatez, a, 1, length(a), 8, carry)
137141
hilbertsort!(a::Vector, lo::Int64, hi::Int64, lim::Int64) = hilbertsort!(backward, backward, backward, coordinatey, a, lo, hi, lim)

0 commit comments

Comments
 (0)