dynamicFriction.jl

end

function RotateVec3D(a::Quaternion, f::Vector3)   
    fx= (f.x==-0) ? 0 : f.x
    fy= (f.y==-0) ? 0 : f.y
    fz= (f.z==-0) ? 0 : f.z
    # fx= f.x
    # fy= f.y
    # fz= f.z
    tw = fx*a.x + fy*a.y + fz*a.z
    tx = fx*a.w - fy*a.z + fz*a.y
    ty = fx*a.z + fy*a.w - fz*a.x
    tz = -fx*a.y + fy*a.x + fz*a.w

    return Vector3((a.w*tx+a.x*tw+a.y*tz-a.z*ty),(a.w*ty-a.x*tz+a.y*tw+a.z*tx),(a.w*tz+a.x*ty-a.y*tx+a.z*tw))
end
#!< Returns a vector representing the specified vector "f" rotated by this quaternion. @param[in] f The vector to transform.

function RotateVec3DInv(a::Quaternion, f::Vector3)  
    fx=f.x
    fy=f.y
    fz=f.z
    tw = a.x*fx + a.y*fy + a.z*fz
    tx = a.w*fx - a.y*fz + a.z*fy
    ty = a.w*fy + a.x*fz - a.z*fx
    tz = a.w*fz - a.x*fy + a.y*fx
    return Vector3((tw*a.x + tx*a.w + ty*a.z - tz*a.y),(tw*a.y - tx*a.z + ty*a.w + tz*a.x),(tw*a.z + tx*a.y - ty*a.x + tz*a.w))
end
#!< Returns a vector representing the specified vector "f" rotated by the inverse of this quaternion. This is the opposite of RotateVec3D. @param[in] f The vector to transform.

function ToRotationVector(a::Quaternion)  
    if (a.w >= 1.0 || a.w <= -1.0) 
        return Vector3(0.0,0.0,0.0)
    end
    squareLength = 1.0-a.w*a.w; # because x*x + y*y + z*z + w*w = 1.0, but more susceptible to w noise (when 
    SLTHRESH_ACOS2SQRT= 2.4e-3; # SquareLength threshhold for when we can use square root optimization for acos. From SquareLength = 1-w*w. (calculate according to 1.0-W_THRESH_ACOS2SQRT*W_THRESH_ACOS2SQRT

    if (squareLength < SLTHRESH_ACOS2SQRT) # ???????
        x=a.x*(2.0*CUDAnative.sqrt((2-2*a.w)/squareLength))
        y=a.y*(2.0*CUDAnative.sqrt((2-2*a.w)/squareLength))
        z=a.z*(2.0*CUDAnative.sqrt((2-2*a.w)/squareLength))
        x=convert(Float64,x)
        y=convert(Float64,y)
        z=convert(Float64,z)
 
        return Vector3(x,y,z) ; # acos(w) = sqrt(2*(1-x)) for w close to 1. for w=0.001, error is 1.317e-6
    else 
        x=a.x*(2.0*CUDAnative.acos(a.w)/CUDAnative.sqrt(squareLength))
        y=a.y*(2.0*CUDAnative.acos(a.w)/CUDAnative.sqrt(squareLength))
        z=a.z*(2.0*CUDAnative.acos(a.w)/CUDAnative.sqrt(squareLength))
        x=convert(Float64,x)
        y=convert(Float64,y)
        z=convert(Float64,z)

        return Vector3(x,y,z)
    end                                    
end 
# !< Returns a rotation vector representing this quaternion rotation. Adapted from http://www.euclideanspace.com/maths/geometry/rotations/conversions/quaternionToAngle/

function FromRotationVector(VecIn::Vector3)
    theta=VecIn*Vector3(0.5,0.5,0.5)
    ntheta=CUDAnative.sqrt((theta.x * theta.x) + (theta.y * theta.y) + (theta.z * theta.z))
    thetaMag2=ntheta*ntheta
    
    DBL_EPSILONx24 =5.328e-15
    if thetaMag2*thetaMag2 < DBL_EPSILONx24
        qw=1.0 - 0.5*thetaMag2
		s=1.0 - thetaMag2 / 6.0
    else
        thetaMag = CUDAnative.sqrt(thetaMag2)
		qw=CUDAnative.cos(thetaMag)
		s=CUDAnative.sin(thetaMag) / thetaMag
    end
    qx=theta.x*s
    qy=theta.y*s
    qz=theta.z*s
    
    qx=convert(Float64,qx)
    qy=convert(Float64,qy)
    qz=convert(Float64,qz)
    qw=convert(Float64,qw)
    
    return Quaternion(qx,qy,qz,qw)
end

function multiplyQuaternions(q::Quaternion,f::Quaternion)
    x=q.x
    y=q.y
    z=q.z
    w=q.w
    x1=w*f.x + x*f.w + y*f.z - z*f.y 
    y1=w*f.y - x*f.z + y*f.w + z*f.x
    z1=w*f.z + x*f.y - y*f.x + z*f.w
    w1=w*f.w - x*f.x - y*f.y - z*f.z

	return Quaternion(x1,y1,z1,w1 ); #!< overload quaternion multiplication.
end

#####################################################
function updateStrain( axialStrain,E) # ?from where strain
    strain = axialStrain # redundant?
    currentTransverseStrainSum=0.0 # ??? todo
    linear=true
    maxStrain=1000000000000000;# ?? todo later change
    if linear
        if axialStrain > maxStrain
            maxStrain = axialStrain # remember this maximum for easy reference
        end
        return stress(axialStrain,E)
    else 
        if (axialStrain > maxStrain) # if new territory on the stress/strain curve
            maxStrain = axialStrain # remember this maximum for easy reference
            returnStress = stress(axialStrain,E) # ??currentTransverseStrainSum
            if (nu != 0.0) 
                strainOffset = maxStrain-stress(axialStrain,E)/(_eHat*(1.0-nu)) # precalculate strain offset for when we back off
            else 
                strainOffset = maxStrain-returnStress/E # precalculate strain offset for when we back off
            end
        else  # backed off a non-linear material, therefore in linear region.
            relativeStrain = axialStrain-strainOffset #  treat the material as linear with a strain offset according to the maximum plastic deformation
            if (nu != 0.0) 
                returnStress = stress(relativeStrain,E)
            else 
                returnStress = E*relativeStrain
            end
        end
        return returnStress
    end
end

function stress( strain , E ) #end,transverseStrainSum, forceLinear){
    #  reference: http://www.colorado.edu/engineering/CAS/courses.d/Structures.d/IAST.Lect05.d/IAST.Lect05.pdf page 10
    #  if (isFailed(strain)) return 0.0f; //if a failure point is set and exceeded, we've broken!
    #   var E =setup.edges[0].stiffness; //todo change later to material ??
    #   var E=1000000;//todo change later to material ??
    #   var scaleFactor=1;
    return E*strain;

    #  #   if (strain <= strainData[1] || linear || forceLinear){ //for compression/first segment and linear materials (forced or otherwise), simple calculation

        #   if (nu==0.0) return E*strain;
        #   else return _eHat*((1-nu)*strain + nu*transverseStrainSum); 
        #  else return eHat()*((1-nu)*strain + nu*transverseStrainSum); 
    #  #  }

      #//the non-linear feature with non-zero poissons ratio is currently experimental
      #int DataCount = modelDataPoints();
      #for (int i=2; i<DataCount; i++){ //go through each segment in the material model (skipping the first segment because it has already been handled.
      #  if (strain <= strainData[i] || i==DataCount-1){ //if in the segment ending with this point (or if this is the last point extrapolate out) 
      #      float Perc = (strain-strainData[i-1])/(strainData[i]-strainData[i-1]);
      #      float basicStress = stressData[i-1] + Perc*(stressData[i]-stressData[i-1]);
      #      if (nu==0.0f) return basicStress;
      #      else { //accounting for volumetric effects
      #          float modulus = (stressData[i]-stressData[i-1])/(strainData[i]-strainData[i-1]);
      #          float modulusHat = modulus/((1-2*nu)*(1+nu));
      #          float effectiveStrain = basicStress/modulus; //this is the strain at which a simple linear stress strain line would hit this point at the definied modulus
      #          float effectiveTransverseStrainSum = transverseStrainSum*(effectiveStrain/strain);
      #          return modulusHat*((1-nu)*effectiveStrain + nu*effectiveTransverseStrainSum);
      #      }
      #  }
      #}

    #  assert(false); //should never reach this point
    #  return 0.0f;
end 

function axialStrain( positiveEnd,strain)
	#strainRatio = pVPos->material()->E/pVNeg->material()->E;
	strainRatio=1.0;
	return positiveEnd ? 2.0 *strain*strainRatio/(1.0+strainRatio) : 2.0*strain/(1.0+strainRatio)
end
#####################################################
function force(N_intForce,N_orient,N_force,static,currentTimeStep) 
    # forces from internal bonds
    totalForce=Vector3(0,0,0)
    # new THREE.Vector3(node.force.x,node.force.y,node.force.z);
    #  todo 


    totalForce=totalForce+N_intForce

    #  for (int i=0; i<6; i++){ 
    #  	if (links[i]) totalForce += links[i]->force(isNegative((linkDirection)i)); # total force in LCS
    #  }
    totalForce = RotateVec3D(N_orient,totalForce); # from local to global coordinates


    # assert(!(totalForce.x != totalForce.x) || !(totalForce.y != totalForce.y) || !(totalForce.z != totalForce.z)); //assert non QNAN

    # other forces
    if(static)
        totalForce=totalForce+N_force
    #  }else if(currentTimeStep<50){
    #  	totalForce.add(new THREE.Vector3(node.force.x,node.force.y,node.force.z));
    else
        #  var ex=0.1;
        #  if(node.force.y!=0){
        #  	var f=400*Math.sin(currentTimeStep*ex);
        #  	totalForce.add(new THREE.Vector3(0,f,0));

        #  }
        #x=N_position[node][3]
        #t=currentTimeStep
        #wave=getForce(x,t)
        #totalForce=totalForce+[0 wave 0]
    end


    #  if (externalExists()) totalForce += external()->force(); //external forces
    #  totalForce -= velocity()*mat->globalDampingTranslateC(); //global damping f-cv
    #  totalForce.z += mat->gravityForce(); //gravity, according to f=mg

    #  if (isCollisionsEnabled()){
    #  	for (std::vector<CVX_Collision*>::iterator it=colWatch->begin(); it!=colWatch->end(); it++){
    #  		totalForce -= (*it)->contactForce(this);
    #  	}
    #  }
    # todo make internal forces 0 again
    # N_intForce[node]=[0 0 0] # do i really need it?

    #  node.force.x=0;
    #  node.force.y=0;
    #  node.force.z=0;


    return totalForce
end

function moment(intMoment,orient,moment) 
    #moments from internal bonds
    totalMoment=Vector3(0,0,0)
    # for (int i=0; i<6; i++){ 
    # 	if (links[i]) totalMoment += links[i]->moment(isNegative((linkDirection)i)); //total force in LCS
    # }

    totalMoment=totalMoment+intMoment
    
    

    totalMoment = RotateVec3D(orient,totalMoment);
    
    

    totalMoment=totalMoment+moment


    #other moments
    # if (externalExists()) totalMoment += external()->moment(); //external moments
    # totalMoment -= angularVelocity()*mat->globalDampingRotateC(); //global damping

    return totalMoment
end
#####################################################
function updateDataAndSave!(metavoxel,setup,fileName,displacements)
    nodes      = setup["nodes"]
    edges      = setup["edges"]
    
    setup["animation"]["showDisplacement"]=true
    voxCount=size(nodes)[1]
    linkCount=size(edges)[1]
    
    N_displacement=Array(metavoxel["N_displacementGPU"])
    N_angle=Array(metavoxel["N_angleGPU"])
    E_stress=Array(metavoxel["E_stressGPU"])
    
    setup["viz"]["maxStress"]=maximum(E_stress)
    setup["viz"]["minStress"]=minimum(E_stress) 


    i=1
	for edge in edges
        edge["stress"]=E_stress[i]
        i=i+1

    end
    
 
    i=1          
	for node in nodes
        node["posTimeSteps"]=[]
        node["angTimeSteps"]=[]
        
        
        node["displacement"]["x"]=N_displacement[i].x/15
        node["displacement"]["y"]=N_displacement[i].y/15
        node["displacement"]["z"]=N_displacement[i].z/15
        
        
        node["angle"]["x"]=N_angle[i].x
        node["angle"]["y"]=N_angle[i].y
        node["angle"]["z"]=N_angle[i].z
        i=i+1

    end
   
        
    for j in 1:length(displacements)
        i=1          
        for node in nodes
            d=displacements[j][i]
            dis = Dict{String, Float64}("x" => d.x/15, "y" => d.y/15,"z" => d.z/15) 
            append!(node["posTimeSteps"],[dis])
            ang = Dict{String, Float64}("x" => N_angle[i].x, "y" => N_angle[i].y,"z" => N_angle[i].z) 
            append!(node["angTimeSteps"],[ang])
            i=i+1
        end
    end
    
    # pass data as a json string (how it shall be displayed in a file)
    stringdata = JSON.json(setup)

    # write the file with the stringdata variable information
    open(fileName, "w") do f
            write(f, stringdata)
         end
    
end
#####################################################
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
            # element=parse(Int,node["id"][2:end])
            N_position[i]=Vector3(node["position"]["x"]*15.0,node["position"]["y"]*15.0,node["position"]["z"]*15.0)
            N_restrained[i]=node["restrained_degrees_of_freedom"][1] ## todo later consider other degrees of freedom
            N_displacement[i]=Vector3(node["displacement"]["x"]*15,node["displacement"]["y"]*15,node["displacement"]["z"]*15)
            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_currPosition[i]=Vector3(node["position"]["x"]*15.0,node["position"]["y"]*15.0,node["position"]["z"]*15.0)

            # for dynamic simulations
            # append!(N_posTimeSteps,[[]])
            # append!(N_angTimeSteps,[[]])

            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"]*15.0 fromNode["position"]["y"]*15.0 fromNode["position"]["z"]*15.0]
            node2 = [toNode["position"]["x"]*15.0 toNode["position"]["y"]*15.0 toNode["position"]["z"]*15.0]

            length=norm(node2-node1)
            axis=normalize(collect(Iterators.flatten(node2-node1)))

            E_source[i]=edge["source"]+1
            E_target[i]=edge["target"]+1
            E_area[i]=edge["area"]
            E_density[i]=edge["density"]
            E_stiffness[i]=edge["stiffness"]
            E_axis[i]=Vector3(axis[1],axis[2],axis[3])
            E_currentRestLength[i]=length #?????? todo change
            # E_currentRestLength[i]=75/sqrt(2)
            

            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


            # for dynamic simulations
            # append!(E_stressTimeSteps,[[]])

            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_posTimeSteps=[]
    # N_angTimeSteps=[]
    N_force=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)

    ############# edges
    E_source=fill(0,linkCount)
    E_target=fill(0,linkCount)
    E_area=fill(0.0f0,linkCount)
    E_density=fill(0.0f0,linkCount)
    E_stiffness=fill(0.0f0,linkCount)
    E_stress=fill(0.0f0,linkCount)
    E_axis=fill(Vector3(1.0,0.0,0.0),linkCount)
    E_currentRestLength=fill(0.0f0,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_currentTransverseStrainSum=fill(0.0f0,linkCount)# TODO remove ot incorporate
    # E_stressTimeSteps=[]


    #################################################################
    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_orientGPU=      CuArray(N_orient)       
    N_edgeIDGPU=      CuArray(N_edgeID)         
    N_edgeFirstGPU=   CuArray(N_edgeFirst)         


    ############# edges
    E_sourceGPU=                    CuArray(E_source)   
    E_targetGPU=                    CuArray(E_target)
    E_areaGPU=                      CuArray(E_area)                             
    E_densityGPU=                   CuArray(E_density)
    E_stiffnessGPU=                 CuArray(E_stiffness)
    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_currentTransverseStrainSumGPU=CuArray(E_currentTransverseStrainSum)
    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_stressTimeSteps=[]


    #########################################
    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_orientGPU" => N_orientGPU,      
        "N_edgeIDGPU" => N_edgeIDGPU,      
        "N_edgeFirstGPU" => N_edgeFirstGPU,
        "E_sourceGPU" =>E_sourceGPU,                    
        "E_targetGPU" =>E_targetGPU,                    
        "E_areaGPU" =>E_areaGPU,                      
        "E_densityGPU" =>E_densityGPU,                   
        "E_stiffnessGPU" =>E_stiffnessGPU,                 
        "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_currentTransverseStrainSumGPU" =>E_currentTransverseStrainSumGPU,
        "E_intForce1GPU" =>E_intForce1GPU,                 
        "E_intMoment1GPU" =>E_intMoment1GPU,                
        "E_intForce2GPU" =>E_intForce2GPU,                 
        "E_intMoment2GPU" =>E_intMoment2GPU,                
        "E_dampGPU" =>E_dampGPU                      
    )

    #########################################
    

    dt=0.0251646
    E = 2000  # MPa
    s=2.38
    mass=10  
    
    
    
    MaxFreq2=E*s/mass
    dt= 1/(6.283185*sqrt(MaxFreq2))
    # dt=0.0001646
    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",displacements)
    end
    println("ran latticeSize $latticeSize with $voxCount voxels 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 getSetup(latticeSize)
    setup = Dict()
    name=string("../json/setupTestUni$latticeSize",".json")
    # open("../json/setupValid2.json", "r") do f
    # open("../json/setupTest.json", "r") do f
    # open("../json/trialJulia.json", "r") do f
    # open("../json/setupTestUni4.json", "r") do f
    # open("../json/setupChiral.json", "r") do f
    # open("../json/setupTestCubeUni10.json", "r") do f
    open(name, "r") do f
        # global setup
        dicttxt = String(read(f))  # file information to string
        setup=JSON.parse(dicttxt)  # parse and transform data
    end

    setup=setup["setup"]
    return setup
end


#####################################################
DDisplacements=[[],[],[],[],[]]
Es=[0.0,0,0,0,0]
latticeSize =2
setup=getSetup(latticeSize)
numTimeSteps=5000
displacements=[]
save=true
returnEvery=1
runMetavoxelGPU!(setup,numTimeSteps,latticeSize,displacements,returnEvery,true)
mmm=length(displacements)
println("num displacements $mmm")
numTimeStepsRecorded=length(displacements)
d=[]
dFEA=[]
j=length(displacements[end])
step=100
for i in 1:step:numTimeStepsRecorded
    append!(d,displacements[i][j].y)
end
DDisplacements[latticeSize]=d

# E2=getYoungsModulus(latticeSize,75,displacements[end],Load,topNodesIndices)

# Es[latticeSize]=E2



println("converged displacement= $(displacements[numTimeStepsRecorded][j].y)")
plot(1:step:numTimeStepsRecorded,d,label="Dynamic",xlabel="timestep",ylabel="displacement",title="$latticeSize Voxel Convergence Study")
# savefig("4_voxel_convergence")


#####################################################

# n=[]
# nFEA=[]
# j=length(displacements[end])
# for i in 1:j
#     append!(n,displacements[end][i].y)
# end
# scatter(1:j,n,label="Dyanmic",xlabel="Node ID",ylabel="displacement",title="Node Displacement")


#####################################################



#####################################################