Skip to content

Commit 4036e52

Browse files
sjkellysjkelly
committed
pull out edge displacement calculation into a function
1 parent ac83daa commit 4036e52

File tree

3 files changed

+52
-81
lines changed

3 files changed

+52
-81
lines changed

‎src/DistMesh.jl

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -85,6 +85,7 @@ Uniform edge length function.
8585
"""
8686
struct HUniform end
8787

88+
include("common.jl")
8889
include("diff.jl")
8990
include("pointdistribution.jl")
9091
include("distmeshnd.jl")

‎src/common.jl

Lines changed: 48 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,48 @@
1+
2+
function compute_displacements!(fh, dp, pair, L, L0, bars, p, setup,
3+
::Type{VertType}) where {VertType}
4+
5+
non_uniform = isa(typeof(L0), AbstractVector)
6+
7+
# compute edge lengths (L) and adaptive edge lengths (L0)
8+
# Lp norm (p=3) is partially computed here
9+
Lsum = zero(eltype(L))
10+
L0sum = non_uniform ? zero(eltype(L0)) : length(pair)
11+
for i in eachindex(pair)
12+
pb = pair[i]
13+
b1 = p[pb[1]]
14+
b2 = p[pb[2]]
15+
barvec = b1 - b2 # bar vector
16+
bars[i] = barvec
17+
L[i] = sqrt(sum(barvec.^2)) # length
18+
non_uniform && (L0[i] = fh((b1+b2)./2))
19+
Lsum = Lsum + L[i].^3
20+
non_uniform && (L0sum = L0sum + L0[i].^3)
21+
end
22+
23+
# zero out force at each node
24+
for i in eachindex(dp)
25+
dp[i] = zero(VertType)
26+
end
27+
28+
# this is not hoisted correctly in the loop so we initialize here
29+
# finish computing the Lp norm (p=3)
30+
lscbrt = (1+(0.4/2^2))*cbrt(Lsum/L0sum)
31+
32+
# Move mesh points based on edge lengths L and forces F
33+
for i in eachindex(pair)
34+
if non_uniform && L[i] < L0[i]*lscbrt || L[i] < lscbrt
35+
L0_f = non_uniform ? L0[i].*lscbrt : lscbrt
36+
# compute force vectors
37+
#F = setup.nonlinear ? (L[i]+L0_f)*(L0_f-L[i])/(2*L0_f) : L0_f-L[i]
38+
F = L0_f-L[i]
39+
# edges are not allowed to pull, only repel
40+
FBar = bars[i].*F./L[i]
41+
# add the force vector to the node
42+
b1 = pair[i][1]
43+
b2 = pair[i][2]
44+
dp[b1] = dp[b1] .+ FBar
45+
dp[b2] = dp[b2] .- FBar
46+
end
47+
end
48+
end

‎src/distmeshnd.jl

Lines changed: 3 additions & 81 deletions
Original file line numberDiff line numberDiff line change
@@ -78,7 +78,7 @@ function distmesh(fdist::Function,
7878
dp = zeros(VertType, length(p)) # displacement at each node
7979
bars = VertType[] # the vector of each edge
8080
L = eltype(VertType)[] # vector length of each edge
81-
non_uniform && (L0 = eltype(VertType)[]) # desired edge length computed by dh (edge length function)
81+
L0 = non_uniform ? eltype(VertType)[] : nothing # desired edge length computed by dh (edge length function)
8282
t = GeometryBasics.SimplexFace{4,Int32}[] # tetrahedra indices from delaunay triangulation
8383
maxmove = typemax(eltype(VertType)) # stores an iteration max movement for retriangulation
8484

@@ -107,46 +107,7 @@ function distmesh(fdist::Function,
107107
stats && push!(statsdata.retriangulations, lcount)
108108
end
109109

110-
# compute edge lengths (L) and adaptive edge lengths (L0)
111-
# Lp norm (p=3) is partially computed here
112-
Lsum = zero(eltype(L))
113-
L0sum = non_uniform ? zero(eltype(L0)) : length(pair)
114-
for i in eachindex(pair)
115-
pb = pair[i]
116-
b1 = p[pb[1]]
117-
b2 = p[pb[2]]
118-
barvec = b1 - b2 # bar vector
119-
bars[i] = barvec
120-
L[i] = sqrt(sum(barvec.^2)) # length
121-
non_uniform && (L0[i] = fh((b1+b2)./2))
122-
Lsum = Lsum + L[i].^3
123-
non_uniform && (L0sum = L0sum + L0[i].^3)
124-
end
125-
126-
# zero out force at each node
127-
for i in eachindex(dp)
128-
dp[i] = zero(VertType)
129-
end
130-
131-
# this is not hoisted correctly in the loop so we initialize here
132-
# finish computing the Lp norm (p=3)
133-
lscbrt = (1+(0.4/2^2))*cbrt(Lsum/L0sum)
134-
135-
# Move mesh points based on edge lengths L and forces F
136-
for i in eachindex(pair)
137-
L0_f = non_uniform ? L0[i].*lscbrt : lscbrt
138-
# compute force vectors
139-
F = L0_f-L[i]
140-
# edges are not allowed to pull, only repel
141-
if F > zero(eltype(L))
142-
FBar = bars[i].*F./L[i]
143-
# add the force vector to the node
144-
b1 = pair[i][1]
145-
b2 = pair[i][2]
146-
dp[b1] = dp[b1] .+ FBar
147-
dp[b2] = dp[b2] .- FBar
148-
end
149-
end
110+
compute_displacements!(fh, dp, pair, L, L0, bars, p, setup, VertType)
150111

151112
# Zero out forces on fix points
152113
if have_fixed
@@ -281,46 +242,7 @@ function distmesh(fdist::Function,
281242
stats && push!(statsdata.retriangulations, lcount)
282243
end
283244

284-
# compute edge lengths (L) and adaptive edge lengths (L0)
285-
# Lp norm (p=3) is partially computed here
286-
Lsum = zero(eltype(L))
287-
L0sum = non_uniform ? zero(eltype(L0)) : length(pair)
288-
for i in eachindex(pair)
289-
pb = pair[i]
290-
b1 = p[pb[1]]
291-
b2 = p[pb[2]]
292-
barvec = b1 - b2 # bar vector
293-
bars[i] = barvec
294-
L[i] = sqrt(sum(barvec.^2)) # length
295-
non_uniform && (L0[i] = fh((b1+b2)./2))
296-
Lsum = Lsum + L[i].^3
297-
non_uniform && (L0sum = L0sum + L0[i].^3)
298-
end
299-
300-
# zero out force at each node
301-
for i in eachindex(dp)
302-
dp[i] = zero(VertType)
303-
end
304-
305-
# this is not hoisted correctly in the loop so we initialize here
306-
# finish computing the Lp norm (p=3)
307-
lscbrt = (1+(0.4/2^2))*cbrt(Lsum/L0sum)
308-
309-
# Move mesh points based on edge lengths L and forces F
310-
for i in eachindex(pair)
311-
if non_uniform && L[i] < L0[i] || L[i] < lscbrt
312-
L0_f = non_uniform ? L0[i].*lscbrt : lscbrt
313-
# compute force vectors
314-
F = L0_f-L[i]
315-
# edges are not allowed to pull, only repel
316-
FBar = bars[i].*F./L[i]
317-
# add the force vector to the node
318-
b1 = pair[i][1]
319-
b2 = pair[i][2]
320-
dp[b1] = dp[b1] .+ FBar
321-
dp[b2] = dp[b2] .- FBar
322-
end
323-
end
245+
compute_displacements!(fh, dp, pair, L, L0, bars, p, setup, VertType)
324246

325247
# Zero out forces on fix points
326248
if have_fixed

0 commit comments

Comments
 (0)