Skip to content
Snippets Groups Projects
run.jl 48.8 KiB
Newer Older
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed
                    linear=true 
                    poisson=false
                    cTE=0.0 #Coefficient of thermal expansion
                    E_material[j]=edgeMaterial(E,mass,nu,rho,2.38,2.38,length,loaded,1.0,linear,poisson,cTE) 

                    
                    j+=1
                    
                end
                
                if jj<(nely-1)&&ii<(nelx-1)
                    edges[j]=deepcopy(edges[1])
                    node1 = [ii*d jj*d 0.0]
                    node2 = [(ii+1)*d (jj+1)*d 0.0]
                    
                    length=d
                    length=norm(node2-node1)
                    axis=normalize(collect(Iterators.flatten(node2-node1)))
                    
                    E_source[j]=i
                    E_target[j]=((ii+1)*nely +(jj+1))+1
                    
                    edges[j]["id"]="e$(j-1)"
                    
                    edges[j]["source"]=E_source[j]-1
                    edges[j]["target"]=E_target[j]-1
                    E_axis[j]=Vector3(axis[1],axis[2],axis[3])
                    E_currentRestLength[j]=length
                    
                    N_edgeID[E_source[j],N_currEdge[E_source[j]]]=j
                    N_edgeFirst[E_source[j],N_currEdge[E_source[j]]]=true
                    N_currEdge[E_source[j]]+=1

                    N_edgeID[E_target[j],N_currEdge[E_target[j]]]=j
                    N_edgeFirst[E_target[j],N_currEdge[E_target[j]]]=false
                    N_currEdge[E_target[j]]+=1

                    #E=(N_material[E_source[j]].E+N_material[E_target[j]].E)/2.0
                    #mass=(N_material[E_source[j]].mass+N_material[E_target[j]].mass)/2.0
                    #nu=(N_material[E_source[j]].nu+N_material[E_target[j]].nu)/2.0
                    #rho=(N_material[E_source[j]].rho+N_material[E_target[j]].rho)/2.0
                    E=(N_material[E_source[j]].E)/1.0
                    mass=(N_material[E_source[j]].mass)/1.0
                    nu=(N_material[E_source[j]].nu)/1.0
                    rho=(N_material[E_source[j]].rho)/1.0
                    
                    edges[j]["material"]["area"]=(E/EE)^(1/p)

                    E_material[j]=edgeMaterial()
                    loaded=0.0
                    linear=true
                    poisson=false
                    cTE=0.0 #Coefficient of thermal expansion
                    E_material[j]=edgeMaterial(E,mass,nu,rho,2.38,2.38,length,loaded,1.0,linear,poisson,cTE) 

                    
                    j+=1
                    
                end
                
            end
        end
    end
    function simulateParallel!(metavoxel,maxNumTimeSteps,dt,saveEvery)
        # initialize(setup)

        for i in 1:maxNumTimeSteps
            doTimeStep!(metavoxel,dt,i)
            if(mod(i,saveEvery)==0&&save)
                updateDataAndSave!(metavoxel,setup,"./../../json/live/$(numFile).json")
                numFile+=1
                if numFile>maxNumFiles
                    numFile=0
                end
            end
        end
    end
    
    ########

    maxNumEdges=10

    ########
    #nelx=20
    #nely=10
    ########
    voxCount=(nelx)*(nely)
    linkCount=(nelx-1)*(nely)+(nelx)*(nely-1)
    linkCount=(nelx-1)*(nely)+(nelx)*(nely-1)+(nelx-1)*(nely-1)
    linkCount=(nelx-1)*(nely)+(nelx)*(nely-1)+(nelx-1)*(nely-1)+(nelx-1)*(nely-1) 
    setup["nodes"]      =  fill(deepcopy(setup["nodes"][1]),voxCount)
    setup["edges"]      =  fill(deepcopy(setup["edges"][1]),linkCount)
    strain =0 #todooo moveeee

    ############# nodes
    N_position=fill(Vector3(),voxCount)
    N_restrained=zeros(Bool, voxCount)
    N_displacement=fill(Vector3(),voxCount)
    N_angle=fill(Vector3(),voxCount)
    N_currPosition=fill(Vector3(),voxCount)
    N_linMom=fill(Vector3(),voxCount)
    N_angMom=fill(Vector3(),voxCount)
    N_intForce=fill(Vector3(),voxCount)
    N_intMoment=fill(Vector3(),voxCount)
    N_moment=fill(Vector3(),voxCount)
    N_force=fill(Vector3(),voxCount)
    N_fixedDisplacement=fill(Vector3(),voxCount)
    N_orient=fill(Quaternion(),voxCount)
    N_edgeID=fill(-1,(voxCount,maxNumEdges))
    N_edgeFirst=fill(true,(voxCount,maxNumEdges))
    N_currEdge=fill(1,voxCount)
    N_material=fill(voxelMaterial(),voxCount)
    N_poissonStrain=fill(Vector3(),voxCount)
    #voxelMaterial(E,mass,nu,rho,zeta,zetaCollision,muStatic,muKinetic,nomSize)

    ############# edges
    E_source=fill(0,linkCount)
    E_target=fill(0,linkCount)
    E_stress=fill(0.0,linkCount)
    E_axis=fill(Vector3(1.0,0.0,0.0),linkCount)
    E_currentRestLength=fill(0.0,linkCount)
    E_pos2=fill(Vector3(),linkCount)
    E_angle1v=fill(Vector3(),linkCount)
    E_angle2v=fill(Vector3(),linkCount)
    E_angle1=fill(Quaternion(),linkCount)
    E_angle2=fill(Quaternion(),linkCount)

    E_intForce1=fill(Vector3(),linkCount)
    E_intMoment1=fill(Vector3(),linkCount) 

    E_intForce2=fill(Vector3(),linkCount)
    E_intMoment2=fill(Vector3(),linkCount)
    E_damp=fill(false,linkCount)
    E_smallAngle=fill(true,linkCount)
    E_material=fill(edgeMaterial(),linkCount)

    E_strain=fill(0.0,linkCount)
    E_maxStrain=fill(0.0,linkCount)
    E_strainOffset=fill(0.0,linkCount)
    E_currentTransverseArea=fill(0.0,linkCount)
    E_currentTransverseStrainSum=fill(0.0,linkCount)# TODO remove ot incorporate


    #################################################################
    initialize!(setup)
    #################################################################

    ########################## turn to cuda arrays
    ############# nodes
    N_positionGPU=      CuArray(N_position)      
    N_restrainedGPU=    CuArray(N_restrained)  
    N_displacementGPU=  CuArray(N_displacement)   
    N_angleGPU=         CuArray(N_angle)       
    N_currPositionGPU=  CuArray(N_currPosition)    
    N_linMomGPU=        CuArray(N_linMom)        
    N_angMomGPU=        CuArray(N_angMom)        
    N_intForceGPU=      CuArray(N_intForce)     
    N_intMomentGPU=     CuArray(N_intMoment)        
    N_momentGPU=        CuArray(N_moment)         
    N_forceGPU=         CuArray(N_force)
    N_fixedDisplacementGPU=CuArray(N_fixedDisplacement)           
    N_orientGPU=        CuArray(N_orient)       
    N_edgeIDGPU=        CuArray(N_edgeID)         
    N_edgeFirstGPU=     CuArray(N_edgeFirst)
    N_materialGPU=      CuArray(N_material)
    N_poissonStrainGPU= CuArray(N_poissonStrain)


    ############# edges
    E_sourceGPU=                    CuArray(E_source)   
    E_targetGPU=                    CuArray(E_target)
    E_stressGPU=                    CuArray(E_stress)
    E_axisGPU=                      CuArray(E_axis)          
    E_currentRestLengthGPU=         CuArray(E_currentRestLength)
    E_pos2GPU=                      CuArray(E_pos2)
    E_angle1vGPU=                   CuArray(E_angle1v)
    E_angle2vGPU=                   CuArray(E_angle2v)
    E_angle1GPU=                    CuArray(E_angle1)
    E_angle2GPU=                    CuArray(E_angle2)

    E_strainGPU=                    CuArray(E_strain)
    E_maxStrainGPU=                 CuArray(E_maxStrain)
    E_strainOffsetGPU=              CuArray(E_strainOffset)
    E_currentTransverseAreaGPU=     CuArray(E_currentTransverseArea)
    E_currentTransverseStrainSumGPU=CuArray(E_currentTransverseStrainSum)# TODO remove ot incorporate

    E_intForce1GPU=                 CuArray(E_intForce1) 
    E_intMoment1GPU=                CuArray(E_intMoment1)  
    E_intForce2GPU=                 CuArray(E_intForce2) 
    E_intMoment2GPU=                CuArray(E_intMoment2)
    E_dampGPU=                      CuArray(E_damp)
    E_smallAngleGPU=                CuArray(E_smallAngle)
    E_materialGPU=                  CuArray(E_material)


    #########################################
    metavoxel = Dict(
        "N_positionGPU" => N_positionGPU,    
        "N_restrainedGPU" => N_restrainedGPU,  
        "N_displacementGPU" => N_displacementGPU,
        "N_angleGPU" => N_angleGPU,       
        "N_currPositionGPU" => N_currPositionGPU,
        "N_linMomGPU" => N_linMomGPU,      
        "N_angMomGPU" => N_angMomGPU,      
        "N_intForceGPU" => N_intForceGPU,    
        "N_intMomentGPU" => N_intMomentGPU,   
        "N_momentGPU" => N_momentGPU,      
        "N_forceGPU" => N_forceGPU,
        "N_fixedDisplacementGPU"=>N_fixedDisplacementGPU,       
        "N_orientGPU" => N_orientGPU,      
        "N_edgeIDGPU" => N_edgeIDGPU,      
        "N_edgeFirstGPU" => N_edgeFirstGPU,
        "N_materialGPU"=>    N_materialGPU,
        "N_poissonStrainGPU"=> N_poissonStrainGPU,

        "E_sourceGPU" =>E_sourceGPU,                    
        "E_targetGPU" =>E_targetGPU,                    
        "E_stressGPU" =>E_stressGPU,                    
        "E_axisGPU" =>E_axisGPU,                      
        "E_currentRestLengthGPU" =>E_currentRestLengthGPU,         
        "E_pos2GPU" =>E_pos2GPU,                      
        "E_angle1vGPU" =>E_angle1vGPU,                   
        "E_angle2vGPU" =>E_angle2vGPU,                   
        "E_angle1GPU" =>E_angle1GPU,                    
        "E_angle2GPU" =>E_angle2GPU,

        "E_strainGPU" =>E_strainGPU,
        "E_maxStrainGPU" =>E_maxStrainGPU,
        "E_strainOffsetGPU"=>E_strainOffsetGPU,
        "E_currentTransverseAreaGPU" =>E_currentTransverseAreaGPU,
        "E_currentTransverseStrainSumGPU" =>E_currentTransverseStrainSumGPU,

        "E_intForce1GPU" =>E_intForce1GPU,                 
        "E_intMoment1GPU" =>E_intMoment1GPU,                
        "E_intForce2GPU" =>E_intForce2GPU,                 
        "E_intMoment2GPU" =>E_intMoment2GPU,                
        "E_dampGPU" =>E_dampGPU,
        "E_smallAngleGPU"=>E_smallAngleGPU,
        "E_materialGPU" =>E_materialGPU
    )

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

    dt=0.0251646
    s=2.38
    E=EE*maximum(x)
    mass=1
    # mass=1e-6  
    MaxFreq2=E*s/mass
    dt= 1/(6.283185*sqrt(MaxFreq2))
    
    numFile=0
    numSaves=0
    
    
    t=@timed doTimeStep!(metavoxel,dt,0)
    time=t[2]
    if save
        println("dt: $dt")
        println("first timestep took $time seconds")
    end
    t=@timed simulateParallel!(metavoxel,maxNumTimeSteps-1,dt,saveEvery)
    time=t[2]
    
    setup["maxNumFiles"]=numSaves
    
    if save
        println("ran $voxCount nodes and $linkCount edges for $maxNumTimeSteps time steps took $time seconds")
    end
    return metavoxel,Array(metavoxel["N_displacementGPU"])
end