Skip to content
Snippets Groups Projects
run.jl 23.8 KiB
Newer Older
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
# Amira Abdel-Rahman
# (c) Massachusetts Institute of Technology 2020

Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed


#################################################################
function runMetavoxelGPULive!(setup,folderPath)

    maxNumTimeSteps=setup["numTimeSteps"]
    maxNumFiles=setup["maxNumFiles"]

    saveEvery=round(maxNumTimeSteps/maxNumFiles)
    maxNumFiles=round(maxNumTimeSteps/saveEvery)-1
    setup["maxNumFiles"]=maxNumFiles

Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
    function initialize!(setup)
        nodes      = setup["nodes"]
        edges      = setup["edges"]

        i=1
        # pre-calculate current position
        for node in nodes
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
            # element=parse(Int,node["id"][2:end])
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
            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"])
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
            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])
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
            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"])

Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
            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)

            # 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)
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
            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)))

Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
            E_source[i]=edge["source"]+1
            E_target[i]=edge["target"]+1
            E_axis[i]=Vector3(axis[1],axis[2],axis[3])
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
            E_currentRestLength[i]=length #?????? todo change
            
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed

            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
            
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
            E=edge["material"]["stiffness"]
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
            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
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
            rho=edge["material"]["density"]
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
            rho=(N_material[E_source[i]].rho+N_material[E_target[i]].rho)/2.0
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed

            
            # E_material[i]=edgeMaterial()

            # E=edge["material"]["stiffness"]
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
            area=edge["material"]["area"]
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
            # mass=1e-6
            
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
            loaded=0.0
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
            if(haskey(edge, "loaded"))
                loaded=edge["loaded"]
            end
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed

            strainRatio=(N_material[E_source[i]].E / N_material[E_target[i]].E)
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
            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
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed

            linear=true
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
            if haskey(setup,"linear")
                linear=setup["linear"]
            end
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
            cTE=0.0 #Coefficient of thermal expansion
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
            if haskey(edge,"cTE")
                cTE=edge["cTE"]
            end
            cTE=N_material[E_source[i]].cTE+N_material[E_target[i]].cTE      
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed

Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
            E_material[i]=edgeMaterial(E,mass,nu,rho,b,h,length,loaded,strainRatio,linear,poisson,cTE)
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed

Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
            if !linear
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
                plasticModulus=1e4
                # yieldStress=1e5
                # failureStress=-1

Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
                yieldStress=1e5
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
                # 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)

Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
            i=i+1
        end 
    end
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
    function simulateParallel!(metavoxel,maxNumTimeSteps,dt)
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
        # initialize(setup)
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
        for i in 1:maxNumTimeSteps
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
            doTimeStep!(metavoxel,dt,i)
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
            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
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
            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)
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
    N_restrained=fill(DOF(), voxCount)
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
    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)
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
    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)
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
    N_fixedDisplacementGPU= CuArray(N_fixedDisplacement)           
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
    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,                    
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
        "E_angle2GPU" =>E_angle2GPU,  
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed

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

    #########################################
    
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
    #todo make recommended timestep a function
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
    dt=0.0251646
    E = setup["edges"][1]["material"]["stiffness"]  # MPa
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
    s=round(sqrt(setup["edges"][1]["material"]["area"] );digits=10)
    # s=E_material[1].L
    # s=setup["voxelSize"]
    mass=N_material[1].mass
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
    MaxFreq2=E*s/mass
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
    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, s: $s, mass: $mass, momentInertiaInverse: $(N_material[1].momentInertiaInverse)")
    
    numFile=0
    numSaves=0
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
    
    
    t=@timed doTimeStep!(metavoxel,dt,0)
    time=t[2]
    println("first timestep took $time seconds")
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
    t=@timed simulateParallel!(metavoxel,maxNumTimeSteps-1,dt)
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
    time=t[2]
    
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
    setup["maxNumFiles"]=numSaves
    

    println("ran $voxCount nodes and $linkCount edges for $maxNumTimeSteps time steps took $time seconds")
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
    return
end

Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
########################################
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed

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"],
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
        metavoxel["E_materialGPU"],
        metavoxel["E_maxStrainGPU"]