parallelFEA.jl

# Amira Abdel-Rahman
# (c) Massachusetts Institute of Technology 2020

using LinearAlgebra
using Plots
import JSON
using StaticArrays, Rotations
# BASED ON https://github.com/jonhiller/Voxelyze


function simulateParallel(setup,numTimeSteps,dt,static=true,saveInterval=10)
	initialize(setup)
    for i in 1:numTimeSteps
        # println("Timestep:",i)
        doTimeStep(setup,dt,static,i,saveInterval);
    end
end

function initialize(setup)
	nodes      = setup["nodes"]
    edges      = setup["edges"]
	# pre-calculate current position
	for node in nodes
        # element=parse(Int,node["id"][2:end])

        append!(N_position,[[node["position"]["x"] node["position"]["y"] node["position"]["z"]]])
        append!(N_degrees_of_freedom,[node["degrees_of_freedom"]])
        append!(N_restrained_degrees_of_freedom, [node["restrained_degrees_of_freedom"]])
        append!(N_displacement,[[node["displacement"]["x"] node["displacement"]["y"] node["displacement"]["z"]]])
        append!(N_angle,[[node["angle"]["x"] node["angle"]["y"] node["angle"]["z"]]])
        append!(N_force,[[node["force"]["x"] node["force"]["y"] node["force"]["z"]]])
        append!(N_currPosition,[[node["position"]["x"] node["position"]["y"] node["position"]["z"]]])
        append!(N_orient,[Quat(1.0,0.0,0.0,0.0)])#quat
        append!(N_linMom,[[0 0 0]])
        append!(N_angMom,[[0 0 0]])
        append!(N_intForce,[[0 0 0]])
        append!(N_intMoment,[[0 0 0]])
        append!(N_moment,[[0 0 0]])
        
        # for dynamic simulations
        append!(N_posTimeSteps,[[]])
        append!(N_angTimeSteps,[[]])
        
	end 
	# 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)))
        
        append!(E_source,[edge["source"]+1])
        append!(E_target,[edge["target"]+1])
        append!(E_area,[edge["area"]])
        append!(E_density,[edge["density"]])
        append!(E_stiffness,[edge["stiffness"]])
        append!(E_stress,[0])
        append!(E_axis,[axis])
        append!(E_currentRestLength,[length])
        append!(E_pos2,[[0 0 0]])
        append!(E_angle1v,[[0 0 0]])
        append!(E_angle2v,[[0 0 0]])
        append!(E_angle1,[Quat(1.0,0,0,0)]) #quat
        append!(E_angle2,[Quat(1.0,0,0,0)]) #quat
        append!(E_currentTransverseStrainSum,[0])
        
        # for dynamic simulations
        append!(E_stressTimeSteps,[[]])
        
	end 
	
end

function doTimeStep(setup,dt,static,currentTimeStep,saveInterval)
    
    nodes      = setup["nodes"]
    edges      = setup["edges"]
    voxCount=size(nodes)[1]
    linkCount=size(edges)[1]

	if dt==0 
		return true
	elseif (dt<0) 
		dt = recommendedTimeStep()
    end

	# if (collisions) updateCollisions();
	collisions=false

	# Euler integration:
	Diverged = false
    #  for edge in edges
    
	for i in 1:linkCount
        # 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"]]
        # updateForces(setup,edge,node1,node2,static)# element numbers??
        updateForces(setup,i,static)# element numbers??
        #  todo: update forces and whatever
		if axialStrain(true) > 100
			Diverged = true; # catch divergent condition! (if any thread sets true we will fail, so don't need mutex...
        end
    end
    
    if Diverged
		println("Divergedd!!!!!")
		return false
    end
                
	for i in 1:voxCount
		timeStep(dt,i,static,currentTimeStep)
        # timeStep(dt,node,static,currentTimeStep)
		if(!static&& currentTimeStep%saveInterval==0)
            append!(N_posTimeSteps[i],[N_displacement[i]])
            append!(N_angTimeSteps[i],[N_angle[i]])

        end
		#  todo: update linMom,angMom, orient and whatever
    end

	
	currentTimeStep = currentTimeStep+dt
	return true
end

function updateForces(setup,edge,static=true)
    
	# pVNeg=new THREE.Vector3(node1.position.x,node1.position.y,node1.position.z);
	# pVPos=new THREE.Vector3(node2.position.x,node2.position.y,node2.position.z);
    # currentRestLength=pVPos.clone().sub(pVNeg).length();
	# edge.currentRestLength=currentRestLength; # todo make sure updated
    
    
	node1=E_source[edge]
    node2=E_target[edge]
    
    currentRestLength=E_currentRestLength[edge]
    
    
	pVNeg=copy(N_currPosition[node1])# todo change to be linked to edge not node 
	pVPos=copy(N_currPosition[node2])# todo change to be linked to edge not node
    
    
    
	#  Vec3D<double> three
	oldPos2 = copy(E_pos2[edge])
    
	oldAngle1v = copy(E_angle1v[edge])
	oldAngle2v =  copy(E_angle2v[edge])# remember the positions/angles from last timestep to calculate velocity
	#  var oldAngle1v=new THREE.Vector3(node1.angle.x,node1.angle.y,node1.angle.z);//?
	#  var oldAngle2v=new THREE.Vector3(node2.angle.x,node2.angle.y,node2.angle.z); //??
    
    
	totalRot= orientLink(edge) # sets pos2, angle1, angle2 /*restLength*/
    
    

	dPos2=0.5*(copy(E_pos2[edge])-oldPos2)
	dAngle1=0.5*(copy(E_angle1v[edge])-oldAngle1v)
	dAngle2=0.5*(copy(E_angle2v[edge])-oldAngle2v)
    
    

	# if volume effects...
    # if (!mat->isXyzIndependent() || currentTransverseStrainSum != 0) 
    # updateTransverseInfo(); //currentTransverseStrainSum != 0 catches when we disable poissons mid-simulation

	
    _stress=updateStrain((E_pos2[edge][1]/E_currentRestLength[edge]),E_stiffness[edge])
    #  var _stress=updateStrain(1.0);
    

	E_stress[edge] = _stress
	if !static
        append!(E_stressTimeSteps[edge],[_stress])
    end
    
    ######### check this
	if setup["viz"]["minStress"]>_stress
		setup["viz"]["minStress"]=_stress
	elseif setup["viz"]["maxStress"]<_stress
		setup["viz"]["maxStress"]=_stress
    end
    

	#  if (isFailed()){forceNeg = forcePos = momentNeg = momentPos = Vec3D<double>(0,0,0); return;}

	#  var b1=mat->_b1, b2=mat->_b2, b3=mat->_b3, a2=mat->_a2; //local copies //todo get from where i had
	
	l   = currentRestLength # ??