parallelFEAGPU.jl

        #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(N_intMoment::Vector3,N_orient::Quaternion,N_moment::Vector3) 
    #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+N_intMoment

    totalMoment = RotateVec3D(N_orient,totalMoment);

    totalMoment=totalMoment+N_moment

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

    return totalMoment
end

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

function updateDataAndSave(setup,fileName)
    nodes      = setup["nodes"]
    edges      = setup["edges"]
    
    setup["animation"]["showDisplacement"]=false
    voxCount=size(nodes)[1]
    linkCount=size(edges)[1]
    
    N_displacement=Array(N_displacementGPU)
    N_angle=Array(N_angleGPU)
    E_stress=Array(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["displacement"]["x"]=N_displacement[i].x
        node["displacement"]["y"]=N_displacement[i].y
        node["displacement"]["z"]=N_displacement[i].z
        
        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
    
    # 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

########################################
###############################################################################################
# times=[]
# nNodes=[]
# nEdges=[]
# latticeSizes=[]
# numTimeStepss=[]


setup = Dict()

latticeSize=4
numTimeSteps=200
save=false

open("../json/setupTestUni4.json", "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"]

########
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_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_stressTimeSteps=[]


dt=0.0251646

t=@timed doTimeStep(dt,true,0)
time=t[2]
println("first timestep took $time seconds")
t=@timed simulateParallel(numTimeSteps,dt)
time=t[2]
println("ran latticeSize $latticeSize with $voxCount voxels and $linkCount edges for $numTimeSteps time steps took $time seconds")
# append!(times, time)

if save
    updateDataAndSave(setup,"../json/trialJuliaParallelGPU.json")
end



# plot(numTimeStepss,times)
###############################################################################################
# plot1 4*4*4 voxels, 300 nodes 960 edges
# numStep=[10,100,200,400,1000]
# firstStepAvg=13.1276
# performanceGPU=[0.012447,0.05713079,0.102851299,0.1864829,0.4725757]
# performanceCPU=[3.014437901,7.700779099,12.9343908,23.8421247,56.1430382]
# plot(numStep,[firstStepAvg.+performanceGPU,performanceCPU],label=["GPU" "CPU"],title="timeSteps")
# # plot2 200 timesteps
# latticeSizes=[4,5,6,7,8,9,10]
# performanceGPU=[0.102851299,0.1021396,0.1044742,0.1043413,0.1611617,0.1674361,0.2076308]
# performanceCPU=[12.934390801,22.9971574,38.838537,60.9359617,87.5866625,128.7116549,173.5449189]
# plot(latticeSizes,[firstStepAvg.+performanceGPU,performanceCPU],label=["GPU" "CPU"],title="latticeSize")