Skip to content

Sjk/opt2 #20

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
wants to merge 3 commits into
base: master
Choose a base branch
from
Draft

Sjk/opt2 #20

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
efbe7e9
add inbounds in some spots
sjkelly Feb 16, 2020
df83e14
some performance improvements, reduce worst-case performance
sjkelly Feb 16, 2020
3b4f567
initial sketch of a better allocation strategy
sjkelly Feb 16, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
30 changes: 18 additions & 12 deletions src/decompositions.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -48,21 +48,27 @@ end
"""
function tet_to_edges!(pair::Vector, pair_set::Set, t)
empty!(pair_set)
for i in eachindex(t)
iszero(length(pair)) && resize!(pair, length(t)*6) # allocate memory on first iteration
ind = 1
@inbounds for i in eachindex(t)
for ep in 1:6
p1 = t[i][tetpairs[ep][1]]
p2 = t[i][tetpairs[ep][2]]
push!(pair_set, p1 > p2 ? (p2,p1) : (p1,p2))
elt = p1 > p2 ? (p2,p1) : (p1,p2)
if !in(bitpack(elt[1],elt[2]), pair_set)
push!(pair_set, bitpack(elt[1],elt[2]))
pair[ind] = elt
ind += 1
end
end
end
resize!(pair, length(pair_set))
# copy pair set to array since sets are not sortable
i = 1
for elt in pair_set
pair[i] = elt
i = i + 1
end

# sort the edge pairs for better point lookup
sort!(pair)
end
# TODO This seems to be really good for some geoemetries,
# but intoduces larger performance regressions for others
#sort!(pair)
return ind-1
end

function bitpack(xi,yi)
(unsafe_trunc(UInt64, xi) << 32) | unsafe_trunc(UInt64, yi)
end
25 changes: 14 additions & 11 deletions src/distmeshnd.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -84,7 +84,7 @@ function distmesh(fdist::Function,
stats ? DistMeshStatistics() : nothing)

# initialize arrays
pair_set = Set{Tuple{Int32,Int32}}() # set used for ensure we have a unique set of edges
pair_set = Set{UInt64}() # 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
Expand All @@ -95,6 +95,7 @@ function distmesh(fdist::Function,
# information on each iteration
lcount = 0 # iteration counter
triangulationcount = 0 # triangulation counter
num_pairs = 0 # number of edge pairs

@inbounds while true
# if large move, retriangulation
Expand All @@ -103,18 +104,20 @@ function distmesh(fdist::Function,
# compute a new delaunay triangulation
retriangulate!(fdist, result, geps, setup, triangulationcount, pt_dists)

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) # 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 triangulationcount == 0
resize!(bars, length(result.tetrahedra)*6)
resize!(L, length(result.tetrahedra)*6)
non_uniform && resize!(L0, length(result.tetrahedra)*6)
end

triangulationcount += 1
stats && push!(result.stats.retriangulations, lcount)
end

compute_displacements!(fh, dp, pair, L, L0, bars, result.points, setup, VertType)
compute_displacements!(fh, dp, pair, num_pairs, L, L0, bars, result.points, setup, VertType)

# Zero out forces on fix points
if have_fixed
Expand Down Expand Up @@ -205,7 +208,7 @@ function retriangulate!(fdist, result::DistMeshResult, geps, setup, triangulatio
# 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
@inbounds 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)
Expand All @@ -215,7 +218,7 @@ function retriangulate!(fdist, result::DistMeshResult, geps, setup, triangulatio
end


function compute_displacements!(fh, dp, pair, L, L0, bars, p, setup,
function compute_displacements!(fh, dp, pair, num_pairs, L, L0, bars, p, setup,
::Type{VertType}) where {VertType}

non_uniform = isa(typeof(L0), AbstractVector)
Expand All @@ -224,7 +227,7 @@ function compute_displacements!(fh, dp, pair, L, L0, bars, p, setup,
# Lp norm (p=3) is partially computed here
Lsum = zero(eltype(L))
L0sum = non_uniform ? zero(eltype(L0)) : length(pair)
for i in eachindex(pair)
@inbounds for i in 1:num_pairs
pb = pair[i]
b1 = p[pb[1]]
b2 = p[pb[2]]
Expand All @@ -237,7 +240,7 @@ function compute_displacements!(fh, dp, pair, L, L0, bars, p, setup,
end

# zero out force at each node
for i in eachindex(dp)
@inbounds for i in eachindex(dp)
dp[i] = zero(VertType)
end

Expand All @@ -246,7 +249,7 @@ function compute_displacements!(fh, dp, pair, L, L0, bars, p, setup,
lscbrt = (1+(0.4/2^2))*cbrt(Lsum/L0sum)

# Move mesh points based on edge lengths L and forces F
for i in eachindex(pair)
@inbounds for i in 1:num_pairs
if non_uniform && L[i] < L0[i]*lscbrt || L[i] < lscbrt
L0_f = non_uniform ? L0[i].*lscbrt : lscbrt
# compute force vectors
Expand Down