I am trying to make an octree for randomly distributed set of particles
<http://arborjs.org/docs/barnes-hut>. I have been trying to find where
Stack Overflow Error appears whole day, but it was in vain. Could you
please help me? The problem should be somewhere in (c == [0]) && (a > 0) part,
however, I have no idea where. Thank you in advance!
type particle
R::Vector{Float64} # coordinates of the particle
m::Float64 # mass of the particle
a_tree::Float64 # acceleration of particle calculated from
tree
a_exact::Float64 # acceleration of particle calculated by
exact summation
end
type node
parent::Int64 # node's parent number
children::Vector{Int64} # set of childrens' numbers
l::Float64 # size of the node
R::Vector{Float64} # radius of geometric center
m::Float64 # total mass of the node
ptcl::Int64 # pointer to the particle in it
end
function problem3_2(N::Int64)
particles = Vector{particle}(N) # Basically, our system
tree = Vector{node}(1) # Initial tree (root)
tree[1] = node(0,[0],1.0,[0.0,0.0,0.0],0.0,0)
for i in eachindex(particles)
particles[i] = particle([rand() - 1/2, rand() - 1/2, rand() - 1/2], 1.0,
0.0, 0.0)
tree = InsertParticle(i,1,tree,particles)
end
return tree
end
function InsertParticle(i::Int64,k::Int64,tree,particles)
c = tree[k].children
a = tree[k].ptcl
if (c == [0]) && (a == 0)
tree[k].ptcl = i
elseif (c == [0]) && (a > 0)
j_1 = length(tree) + DetermOctant(particles[tree[k].ptcl],tree[k])
j_2 = length(tree) + DetermOctant(particles[i],tree[k])
tree[k].children = [x for x in length(tree)+1:length(tree)+8]
for p = 1:8
push!(tree,node(k,[0],tree[k].l/2,[0.0,0.0,0.0],0.0,0))
if (p == 1) || (p == 4) || (p == 5) || (p == 8)
tree[k+p].R[1] = tree[k].R[1] + tree[k+p].l/2
else
tree[k+p].R[1] = tree[k].R[1] - tree[k+p].l/2
end
if (p == 1) || (p == 2) || (p == 5) || (p == 6)
tree[k+p].R[2] = tree[k].R[2] + tree[k+p].l/2
else
tree[k+p].R[2] = tree[k].R[2] - tree[k+p].l/2
end
if (p == 1) || (p == 2) || (p == 3) || (p == 4)
tree[k+p].R[3] = tree[k].R[3] + tree[k+p].l/2
else
tree[k+p].R[3] = tree[k].R[3] - tree[k+p].l/2
end
end
InsertParticle(tree[k].ptcl,j_1,tree,particles)
InsertParticle(i,j_2,tree,particles)
tree[k].ptcl = 0
elseif (c != [0])
j = DetermOctant(particles[i],tree[k])
InsertParticle(i,tree[k].children[j],tree,particles)
end
return tree
end
function DetermOctant(x::particle,y::node)
c1 = y.R[1]; c2 = y.R[2]; c3 = y.R[3]
x1 = sign(x.R[1] - c1); x2 = sign(x.R[2] - c2); x3 = sign(x.R[3] - c3)
if (x1 > 0) && (x2 > 0) && (x3 > 0)
n = 1
elseif (x1 < 0) && (x2 > 0) && (x3 > 0)
n = 2
elseif (x1 < 0) && (x2 < 0) && (x3 > 0)
n = 3
elseif (x1 > 0) && (x2 < 0) && (x3 > 0)
n = 4
elseif (x1 > 0) && (x2 > 0) && (x3 < 0)
n = 5
elseif (x1 < 0) && (x2 > 0) && (x3 < 0)
n = 6
elseif (x1 < 0) && (x2 < 0) && (x3 < 0)
n = 7
else
n = 8
end
return n
end
