# Amira Abdel-Rahman
# (c) Massachusetts Institute of Technology 2020
function runMetavoxelGPU!(setup,numTimeSteps,latticeSize,displacements,returnEvery,save)
function initialize!(setup)
nodes = setup["nodes"]
edges = setup["edges"]
i=1
# pre-calculate current position
for node in nodes
node["position"]["y"]=node["position"]["y"]+0.0
node["position"]["z"]=node["position"]["z"]+0.0
# element=parse(Int,node["id"][2:end])
N_position[i]=Vector3(node["position"]["x"],node["position"]["y"],node["position"]["z"])
N_restrained[i]=node["restrained_degrees_of_freedom"][1] ## todo later consider other degrees of freedom
N_displacement[i]=Vector3(node["displacement"]["x"],node["displacement"]["y"],node["displacement"]["z"])
N_angle[i]=Vector3(node["angle"]["x"],node["angle"]["y"],node["angle"]["z"])
N_force[i]=Vector3(node["force"]["x"],node["force"]["y"],node["force"]["z"])
N_fixedDisplacement[i]=Vector3(node["fixedDisplacement"]["x"],node["fixedDisplacement"]["y"],node["fixedDisplacement"]["z"])
N_currPosition[i]=Vector3(node["position"]["x"],node["position"]["y"],node["position"]["z"])
#voxelMaterial(E,mass,nu,rho,zeta,zetaCollision,muStatic,muKinetic,nomSize)
i=i+1
end
i=1
# pre-calculate the axis
for edge in edges
# element=parse(Int,edge["id"][2:end])
# find the nodes that the lements connects
fromNode = nodes[edge["source"]+1]
toNode = nodes[edge["target"]+1]
node1 = [fromNode["position"]["x"] fromNode["position"]["y"] fromNode["position"]["z"]]
node2 = [toNode["position"]["x"] toNode["position"]["y"] toNode["position"]["z"]]
length=norm(node2-node1)
axis=normalize(collect(Iterators.flatten(node2-node1)))
E_source[i]=edge["source"]+1
E_target[i]=edge["target"]+1
E_axis[i]=Vector3(axis[1],axis[2],axis[3])
E_currentRestLength[i]=length #?????? todo change
N_edgeID[E_source[i],N_currEdge[E_source[i]]]=i
N_edgeFirst[E_source[i],N_currEdge[E_source[i]]]=true
N_currEdge[E_source[i]]+=1
N_edgeID[E_target[i],N_currEdge[E_target[i]]]=i
N_edgeFirst[E_target[i],N_currEdge[E_target[i]]]=false
N_currEdge[E_target[i]]+=1
E=(N_material[E_source[i]].E+N_material[E_target[i]].E)/2.0
mass=(N_material[E_source[i]].mass+N_material[E_target[i]].mass)/2.0
nu=(N_material[E_source[i]].nu+N_material[E_target[i]].nu)/2.0
rho=(N_material[E_source[i]].rho+N_material[E_target[i]].rho)/2.0
area=edge["material"]["area"]
E_material[i]=edgeMaterial()
loaded=0.0
strainRatio=(N_material[E_source[i]].E / N_material[E_target[i]].E)
# linear=(N_material[E_source[i]].linear && N_material[E_target[i]].linear)
# poisson=(N_material[E_source[i]].poisson || N_material[E_target[i]].poisson)
linear=true
poisson=false
cTE=0.0 #Coefficient of thermal expansion
# epsilonFail=(N_material[E_source[i]].epsilonFail+N_material[E_target[i]].epsilonFail)/2.0 #TODO CHANGE TO SMALLEST
E_material[i]=edgeMaterial(E,mass,nu,rho,sqrt(area),sqrt(area),length,loaded,strainRatio,linear,poisson,cTE)
#edgeMaterial(E,mass,nu,rho,b,h,L)
i=i+1
end
end
function simulateParallel!(metavoxel,numTimeSteps,dt,returnEvery)
# initialize(setup)
for i in 1:numTimeSteps
#println("Timestep:",i)
doTimeStep!(metavoxel,dt,i)
if(mod(i,returnEvery)==0)
append!(displacements,[Array(metavoxel["N_displacementGPU"])])
end
end
end
########
voxCount=0
linkCount=0
nodes = setup["nodes"]
edges = setup["edges"]
voxCount=size(nodes)[1]
linkCount=size(edges)[1]
strain =0 #todooo moveeee
maxNumEdges=10
########
voxCount=0
linkCount=0
nodes = setup["nodes"]
edges = setup["edges"]
voxCount=size(nodes)[1]
linkCount=size(edges)[1]
strain =0 #todooo moveeee
############# nodes
N_position=fill(Vector3(),voxCount)
N_restrained=zeros(Bool, voxCount)
N_displacement=fill(Vector3(),voxCount)
N_angle=fill(Vector3(),voxCount)
N_currPosition=fill(Vector3(),voxCount)
N_linMom=fill(Vector3(),voxCount)
N_angMom=fill(Vector3(),voxCount)
N_intForce=fill(Vector3(),voxCount)
N_intMoment=fill(Vector3(),voxCount)
N_moment=fill(Vector3(),voxCount)
N_force=fill(Vector3(),voxCount)
N_fixedDisplacement=fill(Vector3(),voxCount)
N_orient=fill(Quaternion(),voxCount)
N_edgeID=fill(-1,(voxCount,maxNumEdges))
N_edgeFirst=fill(true,(voxCount,maxNumEdges))
N_currEdge=fill(1,voxCount)
N_material=fill(voxelMaterial(),voxCount)
N_poissonStrain=fill(Vector3(),voxCount)
#voxelMaterial(E,mass,nu,rho,zeta,zetaCollision,muStatic,muKinetic,nomSize)
############# edges
E_source=fill(0,linkCount)
E_target=fill(0,linkCount)
E_stress=fill(0.0,linkCount)
E_axis=fill(Vector3(1.0,0.0,0.0),linkCount)
E_currentRestLength=fill(0.0,linkCount)
E_pos2=fill(Vector3(),linkCount)
E_angle1v=fill(Vector3(),linkCount)
E_angle2v=fill(Vector3(),linkCount)
E_angle1=fill(Quaternion(),linkCount)
E_angle2=fill(Quaternion(),linkCount)
E_intForce1=fill(Vector3(),linkCount)
E_intMoment1=fill(Vector3(),linkCount)
E_intForce2=fill(Vector3(),linkCount)
E_intMoment2=fill(Vector3(),linkCount)
E_damp=fill(false,linkCount)
E_smallAngle=fill(true,linkCount)
E_material=fill(edgeMaterial(),linkCount)
E_strain=fill(0.0,linkCount)
E_maxStrain=fill(0.0,linkCount)
E_strainOffset=fill(0.0,linkCount)
E_currentTransverseArea=fill(0.0,linkCount)
E_currentTransverseStrainSum=fill(0.0,linkCount)# TODO remove ot incorporate
#################################################################
initialize!(setup)
#################################################################
########################## turn to cuda arrays
############# nodes
N_positionGPU= CuArray(N_position)
N_restrainedGPU= CuArray(N_restrained)
N_displacementGPU= CuArray(N_displacement)
N_angleGPU= CuArray(N_angle)
N_currPositionGPU= CuArray(N_currPosition)
N_linMomGPU= CuArray(N_linMom)
N_angMomGPU= CuArray(N_angMom)
N_intForceGPU= CuArray(N_intForce)
N_intMomentGPU= CuArray(N_intMoment)
N_momentGPU= CuArray(N_moment)
N_forceGPU= CuArray(N_force)
N_fixedDisplacementGPU=CuArray(N_fixedDisplacement)
N_orientGPU= CuArray(N_orient)
N_edgeIDGPU= CuArray(N_edgeID)
N_edgeFirstGPU= CuArray(N_edgeFirst)
N_materialGPU= CuArray(N_material)
N_poissonStrainGPU= CuArray(N_poissonStrain)
############# edges
E_sourceGPU= CuArray(E_source)
E_targetGPU= CuArray(E_target)
E_stressGPU= CuArray(E_stress)
E_axisGPU= CuArray(E_axis)
E_currentRestLengthGPU= CuArray(E_currentRestLength)
E_pos2GPU= CuArray(E_pos2)
E_angle1vGPU= CuArray(E_angle1v)
E_angle2vGPU= CuArray(E_angle2v)
E_angle1GPU= CuArray(E_angle1)
E_angle2GPU= CuArray(E_angle2)
E_strainGPU= CuArray(E_strain)
E_maxStrainGPU= CuArray(E_maxStrain)
E_strainOffsetGPU= CuArray(E_strainOffset)
E_currentTransverseAreaGPU= CuArray(E_currentTransverseArea)
E_currentTransverseStrainSumGPU=CuArray(E_currentTransverseStrainSum)# TODO remove ot incorporate
E_intForce1GPU= CuArray(E_intForce1)
E_intMoment1GPU= CuArray(E_intMoment1)
E_intForce2GPU= CuArray(E_intForce2)
E_intMoment2GPU= CuArray(E_intMoment2)
E_dampGPU= CuArray(E_damp)
E_smallAngleGPU= CuArray(E_smallAngle)
E_materialGPU= CuArray(E_material)
#########################################
metavoxel = Dict(
"N_positionGPU" => N_positionGPU,
"N_restrainedGPU" => N_restrainedGPU,
"N_displacementGPU" => N_displacementGPU,
"N_angleGPU" => N_angleGPU,
"N_currPositionGPU" => N_currPositionGPU,
"N_linMomGPU" => N_linMomGPU,
"N_angMomGPU" => N_angMomGPU,
"N_intForceGPU" => N_intForceGPU,
"N_intMomentGPU" => N_intMomentGPU,
"N_momentGPU" => N_momentGPU,
"N_forceGPU" => N_forceGPU,
"N_fixedDisplacementGPU"=>N_fixedDisplacementGPU,
"N_orientGPU" => N_orientGPU,
"N_edgeIDGPU" => N_edgeIDGPU,
"N_edgeFirstGPU" => N_edgeFirstGPU,
"N_materialGPU"=> N_materialGPU,
"N_poissonStrainGPU"=> N_poissonStrainGPU,
"E_sourceGPU" =>E_sourceGPU,
"E_targetGPU" =>E_targetGPU,
"E_stressGPU" =>E_stressGPU,
"E_axisGPU" =>E_axisGPU,
"E_currentRestLengthGPU" =>E_currentRestLengthGPU,
"E_pos2GPU" =>E_pos2GPU,
"E_angle1vGPU" =>E_angle1vGPU,
"E_angle2vGPU" =>E_angle2vGPU,
"E_angle1GPU" =>E_angle1GPU,
"E_angle2GPU" =>E_angle2GPU,
"E_strainGPU" =>E_strainGPU,
"E_maxStrainGPU" =>E_maxStrainGPU,
"E_strainOffsetGPU"=>E_strainOffsetGPU,
"E_currentTransverseAreaGPU" =>E_currentTransverseAreaGPU,
"E_currentTransverseStrainSumGPU" =>E_currentTransverseStrainSumGPU,
"E_intForce1GPU" =>E_intForce1GPU,
"E_intMoment1GPU" =>E_intMoment1GPU,
"E_intForce2GPU" =>E_intForce2GPU,
"E_intMoment2GPU" =>E_intMoment2GPU,
"E_dampGPU" =>E_dampGPU,
"E_smallAngleGPU" =>E_smallAngleGPU,
"E_materialGPU" =>E_materialGPU
)
#########################################
#todo change to get min stiffness
dt=0.0251646
E = setup["edges"][1]["material"]["stiffness"] # MPa
s=sqrt(setup["edges"][1]["material"]["area"])
mass=1
# mass=1e-6
MaxFreq2=E*s/mass
dt= 1/(6.283185*sqrt(MaxFreq2))
println("dt: $dt")
append!(displacements,[Array(metavoxel["N_displacementGPU"])])
t=@timed doTimeStep!(metavoxel,dt,0)
append!(displacements,[Array(metavoxel["N_displacementGPU"])])
time=t[2]
println("first timestep took $time seconds")
t=@timed simulateParallel!(metavoxel,numTimeSteps-1,dt,returnEvery)
time=t[2]
if save
updateDataAndSave!(metavoxel,setup,"./../../json/trialJuliaParallelGPUDynamic.json")
end
println("ran latticeSize $latticeSize with $voxCount nodes and $linkCount edges for $numTimeSteps time steps took $time seconds")
return
end
#################################################################
function doTimeStep!(metavoxel,dt,currentTimeStep)
# update forces: go through edges, get currentposition from nodes, calc pos2 and update stresses and interior forces of nodes
run_updateEdges!(dt,currentTimeStep,
metavoxel["E_sourceGPU"],
metavoxel["E_targetGPU"],
metavoxel["E_stressGPU"],
metavoxel["E_axisGPU"],
metavoxel["E_currentRestLengthGPU"],
metavoxel["E_pos2GPU"],
metavoxel["E_angle1vGPU"],
metavoxel["E_angle2vGPU"],
metavoxel["E_angle1GPU"],
metavoxel["E_angle2GPU"],
metavoxel["E_intForce1GPU"],
metavoxel["E_intMoment1GPU"],
metavoxel["E_intForce2GPU"],
metavoxel["E_intMoment2GPU"],
metavoxel["E_dampGPU"],
metavoxel["E_smallAngleGPU"],
metavoxel["E_materialGPU"],
metavoxel["E_strainGPU"],
metavoxel["E_maxStrainGPU"],
metavoxel["E_strainOffsetGPU"],
metavoxel["E_currentTransverseAreaGPU"],
metavoxel["E_currentTransverseStrainSumGPU"],
metavoxel["N_currPositionGPU"],
metavoxel["N_orientGPU"],
metavoxel["N_poissonStrainGPU"])
# update forces: go through nodes and update interior force (according to int forces from edges), integrate and update currpos
run_updateNodes!(dt,currentTimeStep,
metavoxel["N_positionGPU"],
metavoxel["N_restrainedGPU"],
metavoxel["N_displacementGPU"],
metavoxel["N_angleGPU"],
metavoxel["N_currPositionGPU"],
metavoxel["N_linMomGPU"],
metavoxel["N_angMomGPU"],
metavoxel["N_intForceGPU"],
metavoxel["N_intMomentGPU"],
metavoxel["N_forceGPU"],
metavoxel["N_fixedDisplacementGPU"],
metavoxel["N_momentGPU"],
metavoxel["N_orientGPU"],
metavoxel["N_edgeIDGPU"],
metavoxel["N_edgeFirstGPU"],
metavoxel["N_materialGPU"],
metavoxel["N_poissonStrainGPU"],
metavoxel["E_intForce1GPU"],
metavoxel["E_intMoment1GPU"],
metavoxel["E_intForce2GPU"],
metavoxel["E_intMoment2GPU"],
metavoxel["E_axisGPU"],
metavoxel["E_strainGPU"],
metavoxel["E_materialGPU"]
)
end
#################################################################
function runMetavoxelGPULive!(setup,folderPath)
maxNumTimeSteps=setup["numTimeSteps"]
maxNumFiles=setup["maxNumFiles"]
saveEvery=round(maxNumTimeSteps/maxNumFiles)
maxNumFiles=round(maxNumTimeSteps/saveEvery)-1
setup["maxNumFiles"]=maxNumFiles
function initialize!(setup)
nodes = setup["nodes"]
edges = setup["edges"]
i=1
# pre-calculate current position
for node in nodes
# element=parse(Int,node["id"][2:end])
node["position"]["y"]=node["position"]["y"]+0.0
node["position"]["z"]=node["position"]["z"]+0.0
N_position[i]=Vector3(node["position"]["x"],node["position"]["y"],node["position"]["z"])
N_restrained[i]=node["restrained_degrees_of_freedom"][1] ## todo later consider other degrees of freedom
N_displacement[i]=Vector3(node["displacement"]["x"],node["displacement"]["y"],node["displacement"]["z"])
N_angle[i]=Vector3(node["angle"]["x"],node["angle"]["y"],node["angle"]["z"])
N_force[i]=Vector3(node["force"]["x"],node["force"]["y"],node["force"]["z"])
N_fixedDisplacement[i]=Vector3(node["fixedDisplacement"]["x"],node["fixedDisplacement"]["y"],node["fixedDisplacement"]["z"])
N_currPosition[i]=Vector3(node["position"]["x"],node["position"]["y"],node["position"]["z"])
E=2000
E = node["material"]["stiffness"] # MPa
nu=0.0
if haskey(node["material"], "poissonRatio") #todo change material data to nodes
nu= node["material"]["poissonRatio"]
end
# println(nu)
rho=1e3
rho = node["material"]["density"] # MPa
momentInertiaInverse=1.92e-6
inertia=1/momentInertiaInverse
zetaInternal=1.0
zetaGlobal=0.2
if haskey(setup,"globalDamping")
zetaGlobal=setup["globalDamping"]
end
zetaCollision=0.0
muStatic= 2.0
muKinetic= 0.1
nomSize=round(sqrt(node["material"]["area"] );digits=10)
mass=round(nomSize*nomSize*nomSize *rho;digits=10)
# E=1e6
# nu=0.3
# rho=1e3
# mass=1e-6 #nomSize*nomSize*nomSize *rho
# massInverse=1.0/mass
# momentInertiaInverse=1.92e-6
# inertia=1/momentInertiaInverse
# zetaInternal=1.0
# zetaGlobal=0.2
# zetaCollision=0.0
# muStatic= 2.0
# muKinetic= 0.1
# nomSize=0.001
linear=true
poisson=false
if haskey(setup,"linear")
linear=setup["linear"]
end
if haskey(setup,"poisson")
poisson= setup["poisson"]
end
cTE=0.0 #Coefficient of thermal expansion
if haskey(setup,"thermal") # later change for node matrial data
if setup["thermal"]
cTE=node["material"]["cTE"]
end
end
# print("poisson $poisson")
# epsilonFail=E*1000.0
N_material[i]=voxelMaterial(E,mass,nu,rho,zetaInternal,zetaGlobal,zetaCollision,muStatic,muKinetic,nomSize,linear,poisson,cTE)
#voxelMaterial(E,mass,nu,rho,zeta,zetaCollision,muStatic,muKinetic,nomSize)
i=i+1
end
i=1
# pre-calculate the axis
for edge in edges
# element=parse(Int,edge["id"][2:end])
# find the nodes that the lements connects
fromNode = nodes[edge["source"]+1]
toNode = nodes[edge["target"]+1]
node1 = [fromNode["position"]["x"] fromNode["position"]["y"] fromNode["position"]["z"]]
node2 = [toNode["position"]["x"] toNode["position"]["y"] toNode["position"]["z"]]
length=norm(node2-node1)
axis=normalize(collect(Iterators.flatten(node2-node1)))
E_source[i]=edge["source"]+1
E_target[i]=edge["target"]+1
E_axis[i]=Vector3(axis[1],axis[2],axis[3])
E_currentRestLength[i]=length #?????? todo change
N_edgeID[E_source[i],N_currEdge[E_source[i]]]=i
N_edgeFirst[E_source[i],N_currEdge[E_source[i]]]=true
N_currEdge[E_source[i]]+=1
N_edgeID[E_target[i],N_currEdge[E_target[i]]]=i
N_edgeFirst[E_target[i],N_currEdge[E_target[i]]]=false
N_currEdge[E_target[i]]+=1
E=edge["material"]["stiffness"]
E=(N_material[E_source[i]].E+N_material[E_target[i]].E)/2.0
mass=(N_material[E_source[i]].mass+N_material[E_target[i]].mass)/2.0
nu=(N_material[E_source[i]].nu+N_material[E_target[i]].nu)/2.0
rho=edge["material"]["density"]
rho=(N_material[E_source[i]].rho+N_material[E_target[i]].rho)/2.0
# E_material[i]=edgeMaterial()
# E=edge["material"]["stiffness"]
area=edge["material"]["area"]
# mass=1e-6
loaded=0.0
if(haskey(edge, "loaded"))
loaded=edge["loaded"]
end
strainRatio=(N_material[E_source[i]].E / N_material[E_target[i]].E)
linear=(N_material[E_source[i]].linear && N_material[E_target[i]].linear)
poisson=(N_material[E_source[i]].poisson || N_material[E_target[i]].poisson)
# epsilonFail=(N_material[E_source[i]].epsilonFail+N_material[E_target[i]].epsilonFail)/2.0 #TODO CHANGE TO SMALLEST
b=sqrt(area)
h=sqrt(area)
E_currentTransverseArea[i]=b*h
linear=true
if haskey(setup,"linear")
linear=setup["linear"]
end
cTE=0.0 #Coefficient of thermal expansion
if haskey(edge,"cTE")
cTE=edge["cTE"]
end
cTE=N_material[E_source[i]].cTE+N_material[E_target[i]].cTE
E_material[i]=edgeMaterial(E,mass,nu,rho,b,h,length,loaded,strainRatio,linear,poisson,cTE)
if !linear
plasticModulus=5e5
yieldStress=1e5
failureStress=-1.0
E_material[i]=setModelBilinear(E_material[i], plasticModulus, yieldStress)
end
i=i+1
end
end
function simulateParallel!(metavoxel,maxNumTimeSteps,dt)
# initialize(setup)
for i in 1:maxNumTimeSteps
doTimeStep!(metavoxel,dt,i)
if(mod(i,saveEvery)==0)
#append!(displacements,[Array(metavoxel["N_displacementGPU"])])
updateDataAndSave!(metavoxel,setup,"$(folderPath)$(numFile).json")
numFile+=1
if numFile>maxNumFiles
numFile=0
end
end
end
end
########
voxCount=0
linkCount=0
nodes = setup["nodes"]
edges = setup["edges"]
voxCount=size(nodes)[1]
linkCount=size(edges)[1]
strain =0 #todooo moveeee
maxNumEdges=10
########
voxCount=0
linkCount=0
nodes = setup["nodes"]
edges = setup["edges"]
voxCount=size(nodes)[1]
linkCount=size(edges)[1]
strain =0 #todooo moveeee
############# nodes
N_position=fill(Vector3(),voxCount)
N_restrained=zeros(Bool, voxCount)
N_displacement=fill(Vector3(),voxCount)
N_angle=fill(Vector3(),voxCount)
N_currPosition=fill(Vector3(),voxCount)
N_linMom=fill(Vector3(),voxCount)
N_angMom=fill(Vector3(),voxCount)
N_intForce=fill(Vector3(),voxCount)
N_intMoment=fill(Vector3(),voxCount)
N_moment=fill(Vector3(),voxCount)
N_force=fill(Vector3(),voxCount)
N_fixedDisplacement=fill(Vector3(),voxCount)
N_orient=fill(Quaternion(),voxCount)
N_edgeID=fill(-1,(voxCount,maxNumEdges))
N_edgeFirst=fill(true,(voxCount,maxNumEdges))
N_currEdge=fill(1,voxCount)
N_material=fill(voxelMaterial(),voxCount)
N_poissonStrain=fill(Vector3(),voxCount)
#voxelMaterial(E,mass,nu,rho,zeta,zetaCollision,muStatic,muKinetic,nomSize)
############# edges
E_source=fill(0,linkCount)
E_target=fill(0,linkCount)
E_stress=fill(0.0,linkCount)
E_axis=fill(Vector3(1.0,0.0,0.0),linkCount)
E_currentRestLength=fill(0.0,linkCount)
E_pos2=fill(Vector3(),linkCount)
E_angle1v=fill(Vector3(),linkCount)
E_angle2v=fill(Vector3(),linkCount)
E_angle1=fill(Quaternion(),linkCount)
E_angle2=fill(Quaternion(),linkCount)
E_intForce1=fill(Vector3(),linkCount)
E_intMoment1=fill(Vector3(),linkCount)
E_intForce2=fill(Vector3(),linkCount)
E_intMoment2=fill(Vector3(),linkCount)
E_damp=fill(false,linkCount)
E_smallAngle=fill(true,linkCount)
E_material=fill(edgeMaterial(),linkCount)
E_strain=fill(0.0,linkCount)
E_maxStrain=fill(0.0,linkCount)
E_strainOffset=fill(0.0,linkCount)
E_currentTransverseArea=fill(0.0,linkCount)
E_currentTransverseStrainSum=fill(0.0,linkCount)# TODO remove ot incorporate
#################################################################
initialize!(setup)
#################################################################
########################## turn to cuda arrays
############# nodes
N_positionGPU= CuArray(N_position)
N_restrainedGPU= CuArray(N_restrained)
N_displacementGPU= CuArray(N_displacement)
N_angleGPU= CuArray(N_angle)
N_currPositionGPU= CuArray(N_currPosition)
N_linMomGPU= CuArray(N_linMom)
N_angMomGPU= CuArray(N_angMom)
N_intForceGPU= CuArray(N_intForce)
N_intMomentGPU= CuArray(N_intMoment)
N_momentGPU= CuArray(N_moment)
N_forceGPU= CuArray(N_force)
N_fixedDisplacementGPU= CuArray(N_fixedDisplacement)
N_orientGPU= CuArray(N_orient)
N_edgeIDGPU= CuArray(N_edgeID)
N_edgeFirstGPU= CuArray(N_edgeFirst)
N_materialGPU= CuArray(N_material)
N_poissonStrainGPU= CuArray(N_poissonStrain)
############# edges
E_sourceGPU= CuArray(E_source)
E_targetGPU= CuArray(E_target)
E_stressGPU= CuArray(E_stress)
E_axisGPU= CuArray(E_axis)
E_currentRestLengthGPU= CuArray(E_currentRestLength)
E_pos2GPU= CuArray(E_pos2)
E_angle1vGPU= CuArray(E_angle1v)
E_angle2vGPU= CuArray(E_angle2v)
E_angle1GPU= CuArray(E_angle1)
E_angle2GPU= CuArray(E_angle2)
E_strainGPU= CuArray(E_strain)
E_maxStrainGPU= CuArray(E_maxStrain)
E_strainOffsetGPU= CuArray(E_strainOffset)
E_currentTransverseAreaGPU= CuArray(E_currentTransverseArea)
E_currentTransverseStrainSumGPU=CuArray(E_currentTransverseStrainSum)# TODO remove ot incorporate
E_intForce1GPU= CuArray(E_intForce1)
E_intMoment1GPU= CuArray(E_intMoment1)
E_intForce2GPU= CuArray(E_intForce2)
E_intMoment2GPU= CuArray(E_intMoment2)
E_dampGPU= CuArray(E_damp)
E_smallAngleGPU= CuArray(E_smallAngle)
E_materialGPU= CuArray(E_material)
#########################################
metavoxel = Dict(
"N_positionGPU" => N_positionGPU,
"N_restrainedGPU" => N_restrainedGPU,
"N_displacementGPU" => N_displacementGPU,
"N_angleGPU" => N_angleGPU,
"N_currPositionGPU" => N_currPositionGPU,
"N_linMomGPU" => N_linMomGPU,
"N_angMomGPU" => N_angMomGPU,
"N_intForceGPU" => N_intForceGPU,
"N_intMomentGPU" => N_intMomentGPU,
"N_momentGPU" => N_momentGPU,
"N_forceGPU" => N_forceGPU,
"N_fixedDisplacementGPU"=>N_fixedDisplacementGPU,
"N_orientGPU" => N_orientGPU,
"N_edgeIDGPU" => N_edgeIDGPU,
"N_edgeFirstGPU" => N_edgeFirstGPU,
"N_materialGPU"=> N_materialGPU,
"N_poissonStrainGPU"=> N_poissonStrainGPU,
"E_sourceGPU" =>E_sourceGPU,
"E_targetGPU" =>E_targetGPU,
"E_stressGPU" =>E_stressGPU,
"E_axisGPU" =>E_axisGPU,
"E_currentRestLengthGPU" =>E_currentRestLengthGPU,
"E_pos2GPU" =>E_pos2GPU,
"E_angle1vGPU" =>E_angle1vGPU,
"E_angle2vGPU" =>E_angle2vGPU,
"E_angle1GPU" =>E_angle1GPU,
"E_angle2GPU" =>E_angle2GPU,
"E_strainGPU" =>E_strainGPU,
"E_maxStrainGPU" =>E_maxStrainGPU,
"E_strainOffsetGPU"=>E_strainOffsetGPU,
"E_currentTransverseAreaGPU" =>E_currentTransverseAreaGPU,
"E_currentTransverseStrainSumGPU" =>E_currentTransverseStrainSumGPU,
"E_intForce1GPU" =>E_intForce1GPU,
"E_intMoment1GPU" =>E_intMoment1GPU,
"E_intForce2GPU" =>E_intForce2GPU,
"E_intMoment2GPU" =>E_intMoment2GPU,
"E_dampGPU" =>E_dampGPU,
"E_smallAngleGPU" =>E_smallAngleGPU,
"E_materialGPU" =>E_materialGPU
)
#########################################
#todo make recommended timestep a function
dt=0.0251646
E = setup["edges"][1]["material"]["stiffness"] # MPa
s=round(sqrt(setup["edges"][1]["material"]["area"] );digits=10)
s=E_material[1].L
mass=N_material[1].mass
MaxFreq2=E*s/mass
if(setup["poisson"])
# mat->_eHat*currentTransverseArea/((strain+1.0)*currentRestLength)
eHat=E_material[1].eHat
temp=eHat*E_material[1].b*E_material[1].h/((0.0+1.0)*E_material[1].L)
MaxFreq2=temp/mass
end
dt= 1.0/(6.283185*sqrt(MaxFreq2))
# println("E: $(E_material[1].E)")
# println("L: $(E_material[1].L)")
# println("b: $(E_material[1].b)")
# println("h: $(E_material[1].h)")
# println("a1: $(E_material[1].a1)")
# println("a2: $(E_material[1].a2)")
# println("b1: $(E_material[1].b1)")
# println("b2: $(E_material[1].b2)")
# println("b3: $(E_material[1].b3)")
# println("eHat: $(E_material[1].eHat)")
println("dt: $dt")
numFile=0
numSaves=0
t=@timed doTimeStep!(metavoxel,dt,0)
time=t[2]
println("first timestep took $time seconds")
t=@timed simulateParallel!(metavoxel,maxNumTimeSteps-1,dt)
time=t[2]
setup["maxNumFiles"]=numSaves
println("ran $voxCount nodes and $linkCount edges for $maxNumTimeSteps time steps took $time seconds")
return
end
########################################
function getYoungsModulus(latticeSize,voxelSize,disp,Load,topNodesIndices)
F=-Load
l0=voxelSize*latticeSize
A=l0*l0
δl1=-mean( x.y for x in disp[topNodesIndices])
stresses=F/A
strain=δl1/l0
println("Load=$Load")
println("stress=$stresses")
E=stresses/strain
return E
end
####################################
function runMetavoxel2DGPULive!(setup,x,nelx,nely,p,EE,nu,maxNumTimeSteps,saveEvery,maxNumFiles,save)
function initialize!(setup)
nodes = setup["nodes"]
edges = setup["edges"]
d=10.0
for ii in 0:(nelx-1)
for jj in 0:(nely-1)
i=(ii*nely +jj)+1
nodes[i]=deepcopy(nodes[1])
N_position[i]=Vector3(ii*d,jj*d + 2.0*d ,0.0)
nodes[i]["position"]["x"]=N_position[i].x
nodes[i]["position"]["y"]=N_position[i].y
nodes[i]["position"]["z"]=N_position[i].z
nodes[i]["id"]="n$(i-1)"
if((ii==0))
#if((ii==0&&jj==0)||(ii==(nelx-1)&&jj==0))
N_restrained[i]=true## todo later consider other degrees of freedom
nodes[i]["restrained_degrees_of_freedom"]=[true,true,true,true,true,true]
else
N_restrained[i]=false ## todo later consider other degrees of freedom
nodes[i]["restrained_degrees_of_freedom"]=[false,false,false,true,true,true]
end
#if((ii==round(nelx/2)&&jj==(nely-1)))
#if((ii==round(nelx/2)&&jj==(nely-1))||(ii==round(nelx/2)-1&&jj==(nely-1)))
if((ii==round(nelx-3)&&jj==(nely-1)))
N_force[i]=Vector3(0.0,-100.0,0.0)
else
N_force[i]=Vector3(0.0,0.0,0.0)
end
nodes[i]["force"]["x"]=N_force[i].x
nodes[i]["force"]["y"]=N_force[i].y
nodes[i]["force"]["z"]=N_force[i].z
nodes[i]["displacement"]["x"]=0.0
nodes[i]["displacement"]["y"]=0.0
nodes[i]["displacement"]["z"]=0.0
N_material[i]=voxelMaterial()
if ii==(nelx-1)&& jj==(nely-1)
E=EE*x[ii,jj]^p
nodes[i]["size"]=x[ii,jj]
elseif ii==(nelx-1)
E=EE*x[ii,jj+1]^p
nodes[i]["size"]=x[ii,jj+1]
elseif jj==(nely-1)
E=EE*x[ii+1,jj]^p
nodes[i]["size"]=x[ii+1,jj]
else
E=EE*x[ii+1,jj+1]^p
nodes[i]["size"]=x[ii+1,jj+1]
end
mass=10
rho=7.85e-9
momentInertiaInverse=1.92e-6
zeta=1.0
zetaCollision=0.0
zetaGlobal=0.2
zetaInternal=1.0
muStatic= 2.0
muKinetic= 0.1
nomSize=1.0
linear=true
poisson=false
cTE=0.0 #Coefficient of thermal expansion
N_material[i]=voxelMaterial(E,mass,nu,rho,zetaInternal,zetaGlobal,zetaCollision,muStatic,muKinetic,nomSize,linear,poisson,cTE)
N_currPosition[i]=deepcopy(N_position[i])
end
end
j=1
for ii in 0:(nelx-1)
for jj in 0:(nely-1)
i=(ii*nely +jj)+1
if ii<(nelx-1)
edges[j]=deepcopy(edges[1])
node1 = [ii*d jj*d 0.0]
node2 = [(ii+1)*d jj*d 0.0]
length=d
length=norm(node2-node1)
axis=normalize(collect(Iterators.flatten(node2-node1)))
E_source[j]=i
E_target[j]=((ii+1)*nely +jj)+1
edges[j]["id"]="e$(j-1)"
edges[j]["source"]=E_source[j]-1
edges[j]["target"]=E_target[j]-1
E_axis[j]=Vector3(axis[1],axis[2],axis[3])
E_currentRestLength[j]=length
N_edgeID[E_source[j],N_currEdge[E_source[j]]]=j
N_edgeFirst[E_source[j],N_currEdge[E_source[j]]]=true
N_currEdge[E_source[j]]+=1
N_edgeID[E_target[j],N_currEdge[E_target[j]]]=j
N_edgeFirst[E_target[j],N_currEdge[E_target[j]]]=false
N_currEdge[E_target[j]]+=1
#E=(N_material[E_source[j]].E+N_material[E_target[j]].E)/2.0
#mass=(N_material[E_source[j]].mass+N_material[E_target[j]].mass)/2.0
#nu=(N_material[E_source[j]].nu+N_material[E_target[j]].nu)/2.0
#rho=(N_material[E_source[j]].rho+N_material[E_target[j]].rho)/2.0
E=(N_material[E_source[j]].E)/1.0
mass=(N_material[E_source[j]].mass)/1.0
nu=(N_material[E_source[j]].nu)/1.0
rho=(N_material[E_source[j]].rho)/1.0
edges[j]["material"]["area"]=(E/EE)^(1/p)
E_material[j]=edgeMaterial()
loaded=0.0
linear=true
poisson=false
cTE=0.0 #Coefficient of thermal expansion
E_material[j]=edgeMaterial(E,mass,nu,rho,2.38,2.38,length,loaded,1.0,linear,poisson,cTE)
j+=1
end
if jj<(nely-1)
edges[j]=deepcopy(edges[1])
node1 = [ii*d jj*d 0.0]
node2 = [(ii)*d (jj+1)*d 0.0]
length=d
length=norm(node2-node1)
axis=normalize(collect(Iterators.flatten(node2-node1)))
E_source[j]=i
E_target[j]=((ii)*nely +(jj+1))+1
edges[j]["id"]="e$(j-1)"
edges[j]["source"]=E_source[j]-1
edges[j]["target"]=E_target[j]-1
E_axis[j]=Vector3(axis[1],axis[2],axis[3])
E_currentRestLength[j]=length
N_edgeID[E_source[j],N_currEdge[E_source[j]]]=j
N_edgeFirst[E_source[j],N_currEdge[E_source[j]]]=true
N_currEdge[E_source[j]]+=1
N_edgeID[E_target[j],N_currEdge[E_target[j]]]=j
N_edgeFirst[E_target[j],N_currEdge[E_target[j]]]=false
N_currEdge[E_target[j]]+=1
#E=(N_material[E_source[j]].E+N_material[E_target[j]].E)/2.0
#mass=(N_material[E_source[j]].mass+N_material[E_target[j]].mass)/2.0
#nu=(N_material[E_source[j]].nu+N_material[E_target[j]].nu)/2.0
#rho=(N_material[E_source[j]].rho+N_material[E_target[j]].rho)/2.0
E=(N_material[E_source[j]].E)/1.0
mass=(N_material[E_source[j]].mass)/1.0
nu=(N_material[E_source[j]].nu)/1.0
rho=(N_material[E_source[j]].rho)/1.0
edges[j]["material"]["area"]=(E/EE)^(1/p)
E_material[j]=edgeMaterial()
loaded=0.0
poisson=false
linear=true
cTE=0.0 #Coefficient of thermal expansion
E_material[j]=edgeMaterial(E,mass,nu,rho,2.38,2.38,length,loaded,1.0,linear,poisson,cTE)
j+=1
end
if jj>0&&ii<(nelx-1)
edges[j]=deepcopy(edges[1])
node1 = [ii*d jj*d 0.0]
node2 = [(ii+1)*d (jj-1)*d 0.0]
length=d
length=norm(node2-node1)
axis=normalize(collect(Iterators.flatten(node2-node1)))
E_source[j]=i
E_target[j]=((ii+1)*nely +(jj-1))+1
edges[j]["id"]="e$(j-1)"
edges[j]["source"]=E_source[j]-1
edges[j]["target"]=E_target[j]-1
E_axis[j]=Vector3(axis[1],axis[2],axis[3])
E_currentRestLength[j]=length
N_edgeID[E_source[j],N_currEdge[E_source[j]]]=j
N_edgeFirst[E_source[j],N_currEdge[E_source[j]]]=true
N_currEdge[E_source[j]]+=1
N_edgeID[E_target[j],N_currEdge[E_target[j]]]=j
N_edgeFirst[E_target[j],N_currEdge[E_target[j]]]=false
N_currEdge[E_target[j]]+=1
#E=(N_material[E_source[j]].E+N_material[E_target[j]].E)/2.0
#mass=(N_material[E_source[j]].mass+N_material[E_target[j]].mass)/2.0
#nu=(N_material[E_source[j]].nu+N_material[E_target[j]].nu)/2.0
#rho=(N_material[E_source[j]].rho+N_material[E_target[j]].rho)/2.0
E=(N_material[E_source[j]].E)/1.0
mass=(N_material[E_source[j]].mass)/1.0
nu=(N_material[E_source[j]].nu)/1.0
rho=(N_material[E_source[j]].rho)/1.0
edges[j]["material"]["area"]=(E/EE)^(1/p)
E_material[j]=edgeMaterial()
loaded=0.0
linear=true
poisson=false
cTE=0.0 #Coefficient of thermal expansion
E_material[j]=edgeMaterial(E,mass,nu,rho,2.38,2.38,length,loaded,1.0,linear,poisson,cTE)
j+=1
end
if jj<(nely-1)&&ii<(nelx-1)
edges[j]=deepcopy(edges[1])
node1 = [ii*d jj*d 0.0]
node2 = [(ii+1)*d (jj+1)*d 0.0]
length=d
length=norm(node2-node1)
axis=normalize(collect(Iterators.flatten(node2-node1)))
E_source[j]=i
E_target[j]=((ii+1)*nely +(jj+1))+1
edges[j]["id"]="e$(j-1)"
edges[j]["source"]=E_source[j]-1
edges[j]["target"]=E_target[j]-1
E_axis[j]=Vector3(axis[1],axis[2],axis[3])
E_currentRestLength[j]=length
N_edgeID[E_source[j],N_currEdge[E_source[j]]]=j
N_edgeFirst[E_source[j],N_currEdge[E_source[j]]]=true
N_currEdge[E_source[j]]+=1
N_edgeID[E_target[j],N_currEdge[E_target[j]]]=j
N_edgeFirst[E_target[j],N_currEdge[E_target[j]]]=false
N_currEdge[E_target[j]]+=1
#E=(N_material[E_source[j]].E+N_material[E_target[j]].E)/2.0
#mass=(N_material[E_source[j]].mass+N_material[E_target[j]].mass)/2.0
#nu=(N_material[E_source[j]].nu+N_material[E_target[j]].nu)/2.0
#rho=(N_material[E_source[j]].rho+N_material[E_target[j]].rho)/2.0
E=(N_material[E_source[j]].E)/1.0
mass=(N_material[E_source[j]].mass)/1.0
nu=(N_material[E_source[j]].nu)/1.0
rho=(N_material[E_source[j]].rho)/1.0
edges[j]["material"]["area"]=(E/EE)^(1/p)
E_material[j]=edgeMaterial()
loaded=0.0
linear=true
poisson=false
cTE=0.0 #Coefficient of thermal expansion
E_material[j]=edgeMaterial(E,mass,nu,rho,2.38,2.38,length,loaded,1.0,linear,poisson,cTE)
j+=1
end
end
end
end
function simulateParallel!(metavoxel,maxNumTimeSteps,dt,saveEvery)
# initialize(setup)
for i in 1:maxNumTimeSteps
doTimeStep!(metavoxel,dt,i)
if(mod(i,saveEvery)==0&&save)
updateDataAndSave!(metavoxel,setup,"./../../json/live/$(numFile).json")
numFile+=1
if numFile>maxNumFiles
numFile=0
end
end
end
end
########
maxNumEdges=10
########
#nelx=20
#nely=10
########
voxCount=(nelx)*(nely)
linkCount=(nelx-1)*(nely)+(nelx)*(nely-1)
linkCount=(nelx-1)*(nely)+(nelx)*(nely-1)+(nelx-1)*(nely-1)
linkCount=(nelx-1)*(nely)+(nelx)*(nely-1)+(nelx-1)*(nely-1)+(nelx-1)*(nely-1)
setup["nodes"] = fill(deepcopy(setup["nodes"][1]),voxCount)
setup["edges"] = fill(deepcopy(setup["edges"][1]),linkCount)
strain =0 #todooo moveeee
############# nodes
N_position=fill(Vector3(),voxCount)
N_restrained=zeros(Bool, voxCount)
N_displacement=fill(Vector3(),voxCount)
N_angle=fill(Vector3(),voxCount)
N_currPosition=fill(Vector3(),voxCount)
N_linMom=fill(Vector3(),voxCount)
N_angMom=fill(Vector3(),voxCount)
N_intForce=fill(Vector3(),voxCount)
N_intMoment=fill(Vector3(),voxCount)
N_moment=fill(Vector3(),voxCount)
N_force=fill(Vector3(),voxCount)
N_fixedDisplacement=fill(Vector3(),voxCount)
N_orient=fill(Quaternion(),voxCount)
N_edgeID=fill(-1,(voxCount,maxNumEdges))
N_edgeFirst=fill(true,(voxCount,maxNumEdges))
N_currEdge=fill(1,voxCount)
N_material=fill(voxelMaterial(),voxCount)
N_poissonStrain=fill(Vector3(),voxCount)
#voxelMaterial(E,mass,nu,rho,zeta,zetaCollision,muStatic,muKinetic,nomSize)
############# edges
E_source=fill(0,linkCount)
E_target=fill(0,linkCount)
E_stress=fill(0.0,linkCount)
E_axis=fill(Vector3(1.0,0.0,0.0),linkCount)
E_currentRestLength=fill(0.0,linkCount)
E_pos2=fill(Vector3(),linkCount)
E_angle1v=fill(Vector3(),linkCount)
E_angle2v=fill(Vector3(),linkCount)
E_angle1=fill(Quaternion(),linkCount)
E_angle2=fill(Quaternion(),linkCount)
E_intForce1=fill(Vector3(),linkCount)
E_intMoment1=fill(Vector3(),linkCount)
E_intForce2=fill(Vector3(),linkCount)
E_intMoment2=fill(Vector3(),linkCount)
E_damp=fill(false,linkCount)
E_smallAngle=fill(true,linkCount)
E_material=fill(edgeMaterial(),linkCount)
E_strain=fill(0.0,linkCount)
E_maxStrain=fill(0.0,linkCount)
E_strainOffset=fill(0.0,linkCount)
E_currentTransverseArea=fill(0.0,linkCount)
E_currentTransverseStrainSum=fill(0.0,linkCount)# TODO remove ot incorporate
#################################################################
initialize!(setup)
#################################################################
########################## turn to cuda arrays
############# nodes
N_positionGPU= CuArray(N_position)
N_restrainedGPU= CuArray(N_restrained)
N_displacementGPU= CuArray(N_displacement)
N_angleGPU= CuArray(N_angle)
N_currPositionGPU= CuArray(N_currPosition)
N_linMomGPU= CuArray(N_linMom)
N_angMomGPU= CuArray(N_angMom)
N_intForceGPU= CuArray(N_intForce)
N_intMomentGPU= CuArray(N_intMoment)
N_momentGPU= CuArray(N_moment)
N_forceGPU= CuArray(N_force)
N_fixedDisplacementGPU=CuArray(N_fixedDisplacement)
N_orientGPU= CuArray(N_orient)
N_edgeIDGPU= CuArray(N_edgeID)
N_edgeFirstGPU= CuArray(N_edgeFirst)
N_materialGPU= CuArray(N_material)
N_poissonStrainGPU= CuArray(N_poissonStrain)
############# edges
E_sourceGPU= CuArray(E_source)
E_targetGPU= CuArray(E_target)
E_stressGPU= CuArray(E_stress)
E_axisGPU= CuArray(E_axis)
E_currentRestLengthGPU= CuArray(E_currentRestLength)
E_pos2GPU= CuArray(E_pos2)
E_angle1vGPU= CuArray(E_angle1v)
E_angle2vGPU= CuArray(E_angle2v)
E_angle1GPU= CuArray(E_angle1)
E_angle2GPU= CuArray(E_angle2)
E_strainGPU= CuArray(E_strain)
E_maxStrainGPU= CuArray(E_maxStrain)
E_strainOffsetGPU= CuArray(E_strainOffset)
E_currentTransverseAreaGPU= CuArray(E_currentTransverseArea)
E_currentTransverseStrainSumGPU=CuArray(E_currentTransverseStrainSum)# TODO remove ot incorporate
E_intForce1GPU= CuArray(E_intForce1)
E_intMoment1GPU= CuArray(E_intMoment1)
E_intForce2GPU= CuArray(E_intForce2)
E_intMoment2GPU= CuArray(E_intMoment2)
E_dampGPU= CuArray(E_damp)
E_smallAngleGPU= CuArray(E_smallAngle)
E_materialGPU= CuArray(E_material)
#########################################
metavoxel = Dict(
"N_positionGPU" => N_positionGPU,
"N_restrainedGPU" => N_restrainedGPU,
"N_displacementGPU" => N_displacementGPU,
"N_angleGPU" => N_angleGPU,
"N_currPositionGPU" => N_currPositionGPU,
"N_linMomGPU" => N_linMomGPU,
"N_angMomGPU" => N_angMomGPU,
"N_intForceGPU" => N_intForceGPU,
"N_intMomentGPU" => N_intMomentGPU,
"N_momentGPU" => N_momentGPU,
"N_forceGPU" => N_forceGPU,
"N_fixedDisplacementGPU"=>N_fixedDisplacementGPU,
"N_orientGPU" => N_orientGPU,
"N_edgeIDGPU" => N_edgeIDGPU,
"N_edgeFirstGPU" => N_edgeFirstGPU,
"N_materialGPU"=> N_materialGPU,
"N_poissonStrainGPU"=> N_poissonStrainGPU,
"E_sourceGPU" =>E_sourceGPU,
"E_targetGPU" =>E_targetGPU,
"E_stressGPU" =>E_stressGPU,
"E_axisGPU" =>E_axisGPU,
"E_currentRestLengthGPU" =>E_currentRestLengthGPU,
"E_pos2GPU" =>E_pos2GPU,
"E_angle1vGPU" =>E_angle1vGPU,
"E_angle2vGPU" =>E_angle2vGPU,
"E_angle1GPU" =>E_angle1GPU,
"E_angle2GPU" =>E_angle2GPU,
"E_strainGPU" =>E_strainGPU,
"E_maxStrainGPU" =>E_maxStrainGPU,
"E_strainOffsetGPU"=>E_strainOffsetGPU,
"E_currentTransverseAreaGPU" =>E_currentTransverseAreaGPU,
"E_currentTransverseStrainSumGPU" =>E_currentTransverseStrainSumGPU,
"E_intForce1GPU" =>E_intForce1GPU,
"E_intMoment1GPU" =>E_intMoment1GPU,
"E_intForce2GPU" =>E_intForce2GPU,
"E_intMoment2GPU" =>E_intMoment2GPU,
"E_dampGPU" =>E_dampGPU,
"E_smallAngleGPU"=>E_smallAngleGPU,
"E_materialGPU" =>E_materialGPU
)
#########################################
dt=0.0251646
s=2.38
E=EE*maximum(x)
mass=1
# mass=1e-6
MaxFreq2=E*s/mass
dt= 1/(6.283185*sqrt(MaxFreq2))
numFile=0
numSaves=0
t=@timed doTimeStep!(metavoxel,dt,0)
time=t[2]
if save
println("dt: $dt")
println("first timestep took $time seconds")
end
t=@timed simulateParallel!(metavoxel,maxNumTimeSteps-1,dt,saveEvery)
time=t[2]
setup["maxNumFiles"]=numSaves
if save
println("ran $voxCount nodes and $linkCount edges for $maxNumTimeSteps time steps took $time seconds")
end
return metavoxel,Array(metavoxel["N_displacementGPU"])
end