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 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