# Amira Abdel-Rahman
# (c) Massachusetts Institute of Technology 2020
#################################################################
function runMetaVoxels!(setup,folderPath,sys="GPU")
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
for node in nodes
if node["parent"]=="11" || node["parent"]=="" #to ignore child nodes that are there for visualization
# 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]=DOF(node["restrained_degrees_of_freedom"][1],node["restrained_degrees_of_freedom"][2],node["restrained_degrees_of_freedom"][3],node["restrained_degrees_of_freedom"][4],node["restrained_degrees_of_freedom"][5],node["restrained_degrees_of_freedom"][6])
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)
nomSize=nomSize*2.0
mass=round(nomSize*nomSize*nomSize *rho;digits=10)
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
scale =0
if haskey(setup,"multiscale")
if setup["multiscale"]
scale=node["scale"]
end
nomSize=nomSize*scale
mass=round(nomSize*nomSize*nomSize *rho;digits=10)
end
N_material[i]=voxelMaterial(E,mass,nu,rho,zetaInternal,zetaGlobal,zetaCollision,muStatic,muKinetic,nomSize,linear,poisson,cTE,scale)
i=i+1
else
node["position"]["x"]= node["position"]["x"]*setup["scale"] #for export
node["position"]["y"]= node["position"]["y"]*setup["scale"] #for export
node["position"]["z"]= node["position"]["z"]*setup["scale"] #for export
end
end
i=1
for edge in edges
# element=parse(Int,edge["id"][2:end])
# find the nodes that the elements connects
fromNode = nodes[edge["source"]+1]
toNode = nodes[edge["target"]+1]
sourceRelPos=Vector3(0,0,0)
targetRelPos=Vector3(0,0,0)
if haskey(setup,"multiscale")
E_sourceNodalCoordinate[i]=edge["sourceNodalCoordinate"]
E_targetNodalCoordinate[i]=edge["targetNodalCoordinate"]
sourceRelPos=getRelativePosition(E_sourceNodalCoordinate[i],setup["voxelSize"]/2.0)
targetRelPos=getRelativePosition(E_targetNodalCoordinate[i],setup["voxelSize"]/2.0)
end
node1 = [fromNode["position"]["x"]+sourceRelPos.x fromNode["position"]["y"]+sourceRelPos.y fromNode["position"]["z"]+sourceRelPos.z]
node2 = [toNode["position"]["x"]+targetRelPos.x toNode["position"]["y"]+targetRelPos.y toNode["position"]["z"]+targetRelPos.z]
length=norm(node2-node1)
axis=normalize(collect(Iterators.flatten(node2-node1))) # pre-calculate the axis
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
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
if haskey(edge,"material")
if haskey(edge["material"],"cTE")
if edge["material"]["cTE"]>0.0
cTE=edge["material"]["cTE"]
end
end
end
E_material[i]=edgeMaterial(E,mass,nu,rho,b,h,length,loaded,strainRatio,linear,poisson,cTE)
#non linear from data
# todo cleanup
if !linear
plasticModulus=1e4
# yieldStress=1e5
# failureStress=-1
yieldStress=1e5
# failureStress=23e4
failureStress=1.1e5
# E_material[i]=setModelBilinear(E_material[i], plasticModulus, yieldStress,failureStress)
strainData=[0.0, 0.1,0.5, 1.1]
stressData=[0.0, 100000.0,500000.0, 110000.0]
strainData=[0.0, 0.1, 0.5, 0.9, -1]
stressData=[0.0, 100000.0,500000.0,450000.0, -1]
strainData=[0.0, 0.1, 0.3, 0.7, -1]
stressData=[0.0, 1e5,5e5,6e5, -1]
strainData=[0.0, 0.1, 0.5, 0.9, -1]
stressData=[0.0, 1e5,5e5,10e5, -1]
strainData=[0.0,0.1,0.25, 0.3, 0.95, -1]
stressData=[0.0,1e5,25e4, 5e5, 6e5, -1]
strainData=[0.0,0.1,0.25, 0.3, 0.5, -1]
stressData=[0.0,1e5,25e4, 25e4, 25e4, -1]
strainData=[0.0, 1.0*0.25, 1.0*0.5, 1.0*0.75, 1.0*0.9, -1]
stressData=[0.0, 1e5*0.25, 1e5*0.5, 1e5*0.75, 1e5*0.9, -1]
strainData=[0.0, 1.0*0.1, 1.0*0.2, 1.0*0.4, 1.0*0.75, -1]
stressData=[0.0, 1e5*0.1, 1e5*0.2, 5e5*0.4, 5e5*0.75, -1]
dataPointCount=8
strainData=[0.0, 1.0*0.1, 1.0*0.2, 1.0*0.4, 1.0*0.6, 1.0*0.7, 1.0*0.9, -1]
stressData=[0.0, 1e5*0.1, 1e5*0.2, 2e5 , 3.6e5 , 5e5 , 7.2e5 , -1]
# strainData=[0.0, 1.0*0.1, 1.0*0.3, 1.0*0.5, 1.0*0.7, 1.0*0.9, 1.0*0.95, -1]
# youngsData=[0.0, 0.0, 1e5, 9e5, 8e5, 1.3e6, 1.15e6,4e5,1e6]
strainData=[0.0, 1.0*0.1, 1.0*0.2, 1.0*0.4, 1.0*0.6, 1.0*0.8, 1.0*0.95, 1.0]
youngsData=[0.0, 0.0, 1e6, 1e6, 1e6, 9e5, 8e5, 7e5]
stressData=[0.0, 1e5*0.1, 1e5*0.2, 1e5*2.2, 1e5*4.1, 1e5*6.0, 1e5*7.1, 1e5*7.5]
strainData=[0.0, 1.0*0.1, 1.0*0.2, 1.0*0.4, 1.0*0.6, 1.0*0.8, 1.0*0.95, 1.0]
youngsData=[0.0, 0.0, 1e6, 1e6, 1e6, 9e5, 8e5, 7e5]
stressData=[0.0, 1e5*0.1, 1e5*0.2]
for i= 4:dataPointCount
append!(stressData,[youngsData[i]*(strainData[i]-strainData[i-1])+stressData[i-1]])
end
# append!(stressData,[-1])
# append!(stressData,[-1])
strainData=[0.0, 1.0*0.1, 1.0*0.2, 1.0*0.4, 1.0*0.6, 1.0*0.8, 1.0*0.95, -1.0]
youngsData=[0.0, 0.0, 1e4, 5e4, 1e5, 5e5, 8e5, 1e5]
stressData=[0.0, 1e5*0.1]
for i= 3:dataPointCount-1
append!(stressData,[youngsData[i]*(strainData[i]-strainData[i-1])+stressData[i-1]])
end
append!(stressData,[-1])
if i==1
println(stressData)
end
strainData=[0.0, 1.0*0.1, 1.0*0.2, 1.0*0.4, 1.0*0.6, 1.0*0.8, 1.0*0.95, -1.0]
stressData=[0.0, 0.05, 0.1, 0.4, 0.6, 0.9, 1.2, -1.0e-5].*1e5
strainData=[0.0, 1.0*0.05, 1.0*0.1, 1.0*0.2, 1.0*0.3, 1.0*0.4, 1.0*0.5, -1]
stressData=[0.0, 1e5*0.05, 1e5*0.1, 1e5*0.2, 1e5*0.3, 1e5*0.4, 1e5*0.5, -1]
strainData=[0.0, 0.1, 0.3, 0.7,0.75,0.8, 0.9, -1]
stressData=[0.0, 1e5, 5e5, 6e5,6.7e5,7e5,7.5e5, -1]
dataPointCount=6
# strainData=[0.0, 0.2, 0.4, 0.7, 1.2, 1.8]
# stressData=[0.0, 2e5 ,4e5, 7e5, 10e5, 18e5]
# strainData=[0.0, 0.2, 0.3, 0.6, 0.9, 1.8]
# stressData=[0.0, 2e5 ,3e5, 8e5, 11e5, 18e5]
# strainData=[0.0, 0.2, 0.4, 0.6, 0.8, 1.8]
# stressData=[0.0, 2e5 ,4e5, 9e5, 8e5, 10e5]
strainData=[0.0, 0.2, 0.4, 0.7, 1.2, 1.8]
stressData=[0.0, 2e5 ,4e5, 4.5e5, 5e5, 6e5]
strainData=[0.0, 0.1, 0.2, 0.5, 0.6, 1.2]
stressData=[0.0, 1e5 ,2e5, 2.5e5, 1.0e5, 0.7e5]
strainData=[0.0, 0.1, 0.2, 0.5, 0.8, 1.2]
stressData=[0.0, 1e5 ,2e5, 3.0e5, 7.0e5, 12e5]
# if i==1
# strainn=0.1
# nu=0.3
# DataCount =dataPointCount
# poisson=true
# transverseStrainSumm=0.001
# for i = 3:DataCount #(i=2; i<DataCount; i++) #go through each segment in the material model (skipping the first segment because it has already been handled.
# println(i)
# if (strainn <= strainData[i] || i==DataCount) #if in the segment ending with this point (or if this is the last point extrapolate out)
# Perc = (strainn-strainData[i-1])/(strainData[i]-strainData[i-1]);
# println("Perc: $Perc")
# basicStress = stressData[i-1] + Perc*(stressData[i]-stressData[i-1]);
# println("basicStress: $basicStress")
# if (!poisson || nu == 0.0)
# basicStress;
# else #accounting for volumetric effects
# modulus = (stressData[i]-stressData[i-1])/(strainData[i]-strainData[i-1]);
# println("modulus: $modulus")
# modulusHat = modulus/((1.0-2.0*nu)*(1.0+nu));
# effectiveStrain = basicStress/modulus; #this is the strain at which a simple linear stress strain line would hit this point at the definied modulus
# effectiveTransverseStrainSum = transverseStrainSumm*(effectiveStrain/strainn);
# m=modulusHat*((1.0-nu)*effectiveStrain + nu*effectiveTransverseStrainSum);
# println("m: $m")
# println()
# end
# end
# end
# end
E_material[i]=setModel(E_material[i],dataPointCount, strainData, stressData)
# println(E_material[i].strainData)
# println(E_material[i].stressData)
# println(E_material[i].epsilonFail)
# println(length)
end
E_source[i]=edge["source"]+1
E_target[i]=edge["target"]+1
i=i+1
end
end
function simulateParallel!(metavoxel,maxNumTimeSteps,dt,sys="GPU")
for i in 1:maxNumTimeSteps
doTimeStep!(metavoxel,dt,i,sys)
if(mod(i,saveEvery)==0)
#append!(displacements,[Array(metavoxel["N_displacementGPU"])])
updateDataAndSave!(metavoxel,setup,"$(folderPath)$(numFile).json",sys)
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=24
############# nodes
N_position=fill(Vector3(),voxCount)
N_restrained=fill(DOF(), 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)
############# edges
E_source=fill(0,linkCount)
E_target=fill(0,linkCount)
E_sourceNodalCoordinate=fill(0,linkCount)
E_targetNodalCoordinate=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)
#################################################################
if (sys=="GPU")
########################## 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_sourceNodalCoordinateGPU= CuArray(E_sourceNodalCoordinate)
E_targetNodalCoordinateGPU= CuArray(E_targetNodalCoordinate)
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_sourceNodalCoordinateGPU" =>E_sourceNodalCoordinateGPU,
"E_targetNodalCoordinateGPU" =>E_targetNodalCoordinateGPU,
"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
)
else
# CPU
metavoxel = Dict(
"N_position" => N_position,
"N_restrained" => N_restrained,
"N_displacement" => N_displacement,
"N_angle" => N_angle,
"N_currPosition" => N_currPosition,
"N_linMom" => N_linMom,
"N_angMom" => N_angMom,
"N_intForce" => N_intForce,
"N_intMoment" => N_intMoment,
"N_moment" => N_moment,
"N_force" => N_force,
"N_fixedDisplacement"=>N_fixedDisplacement,
"N_orient" => N_orient,
"N_edgeID" => N_edgeID,
"N_edgeFirst" => N_edgeFirst,
"N_material"=> N_material,
"N_poissonStrain"=> N_poissonStrain,
"E_source" =>E_source,
"E_target" =>E_target,
"E_sourceNodalCoordinate" =>E_sourceNodalCoordinate,
"E_targetNodalCoordinate" =>E_targetNodalCoordinate,
"E_stress" =>E_stress,
"E_axis" =>E_axis,
"E_currentRestLength" =>E_currentRestLength,
"E_pos2" =>E_pos2,
"E_angle1v" =>E_angle1v,
"E_angle2v" =>E_angle2v,
"E_angle1" =>E_angle1,
"E_angle2" =>E_angle2,
"E_strain" =>E_strain,
"E_maxStrain" =>E_maxStrain,
"E_strainOffset"=>E_strainOffset,
"E_currentTransverseArea" =>E_currentTransverseArea,
"E_currentTransverseStrainSum" =>E_currentTransverseStrainSum,
"E_intForce1" =>E_intForce1,
"E_intMoment1" =>E_intMoment1,
"E_intForce2" =>E_intForce2,
"E_intMoment2" =>E_intMoment2,
"E_damp" =>E_damp,
"E_smallAngle" =>E_smallAngle,
"E_material" =>E_material
)
end
#########################################
#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
# s=setup["voxelSize"]
mass=N_material[1].mass
if(!setup["hierarchical"])
s=E_material[1].L
mass=round(setup["edges"][1]["material"]["area"]*E_material[1].L *setup["edges"][1]["material"]["density"];digits=10)
end
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))
if (logging)
println("dt: $dt, s: $s, mass: $mass, momentInertiaInverse: $(N_material[1].momentInertiaInverse)")
end
numFile=0
numSaves=0
t=@timed doTimeStep!(metavoxel,dt,0,sys)
time=t[2]
if (logging)
println("first timestep took $time seconds")
end
t=@timed simulateParallel!(metavoxel,maxNumTimeSteps-1,dt,sys)
time=t[2]
setup["maxNumFiles"]=numSaves
if (logging)
println("ran $voxCount nodes and $linkCount edges for $maxNumTimeSteps time steps took $time seconds")
end
return
end
########################################
function doTimeStep!(metavoxel,dt,currentTimeStep,sys="GPU")
if (sys=="GPU")
# update forces: go through edges, get currentposition from nodes, calc pos2 and update stresses and interior forces of nodes
run_updateEdgesGPU!(dt,currentTimeStep,
metavoxel["E_sourceGPU"],
metavoxel["E_targetGPU"],
metavoxel["E_sourceNodalCoordinateGPU"],
metavoxel["E_targetNodalCoordinateGPU"],
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_updateNodesGPU!(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"],
metavoxel["E_maxStrainGPU"]
)
else
# CPU
# update forces: go through edges, get currentposition from nodes, calc pos2 and update stresses and interior forces of nodes
updateEdgesCPU!(dt,currentTimeStep,
metavoxel["E_source"],
metavoxel["E_target"],
metavoxel["E_sourceNodalCoordinate"],
metavoxel["E_targetNodalCoordinate"],
metavoxel["E_stress"],
metavoxel["E_axis"],
metavoxel["E_currentRestLength"],
metavoxel["E_pos2"],
metavoxel["E_angle1v"],
metavoxel["E_angle2v"],
metavoxel["E_angle1"],
metavoxel["E_angle2"],
metavoxel["E_intForce1"],
metavoxel["E_intMoment1"],
metavoxel["E_intForce2"],
metavoxel["E_intMoment2"],
metavoxel["E_damp"],
metavoxel["E_smallAngle"],
metavoxel["E_material"],
metavoxel["E_strain"],
metavoxel["E_maxStrain"],
metavoxel["E_strainOffset"],
metavoxel["E_currentTransverseArea"],
metavoxel["E_currentTransverseStrainSum"],
metavoxel["N_currPosition"],
metavoxel["N_orient"],
metavoxel["N_poissonStrain"]
)
# update forces: go through nodes and update interior force (according to int forces from edges), integrate and update currpos
updateNodesCPU!(dt,currentTimeStep,
metavoxel["N_position"],
metavoxel["N_restrained"],
metavoxel["N_displacement"],
metavoxel["N_angle"],
metavoxel["N_currPosition"],
metavoxel["N_linMom"],
metavoxel["N_angMom"],
metavoxel["N_intForce"],
metavoxel["N_intMoment"],
metavoxel["N_force"],
metavoxel["N_fixedDisplacement"],
metavoxel["N_moment"],
metavoxel["N_orient"],
metavoxel["N_edgeID"],
metavoxel["N_edgeFirst"],
metavoxel["N_material"],
metavoxel["N_poissonStrain"],
metavoxel["E_intForce1"],
metavoxel["E_intMoment1"],
metavoxel["E_intForce2"],
metavoxel["E_intMoment2"],
metavoxel["E_axis"],
metavoxel["E_strain"],
metavoxel["E_material"],
metavoxel["E_maxStrain"]
)
end
end