run.jl

# 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