Skip to content
Snippets Groups Projects
parallelFEA.jl 39.5 KiB
Newer Older
Amira Abdel-Rahman's avatar
Amira Abdel-Rahman committed

	# assert(!(curForce.x != curForce.x) || !(curForce.y != curForce.y) || !(curForce.z != curForce.z)); //assert non QNAN
	linMom=linMom+curForce*dt
    

	massInverse=8e-6 # todo ?? later change
	translate=linMom*(dt*massInverse) # ??massInverse

    #  //	we need to check for friction conditions here (after calculating the translation) and stop things accordingly
	# if (isFloorEnabled() && floorPenetration() >= 0){ //we must catch a slowing voxel here since it all boils down to needing access to the dt of this timestep.
	#  	double work = fricForce.x*translate.x + fricForce.y*translate.y; //F dot disp
	#  	double hKe = 0.5*mat->_massInverse*(linMom.x*linMom.x + linMom.y*linMom.y); //horizontal kinetic energy

	#  	if(hKe + work <= 0) setFloorStaticFriction(true); //this checks for a change of direction according to the work-energy principle

	#  	if (isFloorStaticFriction()){ //if we're in a state of static friction, zero out all horizontal motion
	#  		linMom.x = linMom.y = 0;
	#  		translate.x = translate.y = 0;
	#  	}
	#  }
	#  else setFloorStaticFriction(false);
    pos=pos+translate
    
    
    N_currPosition[node]=copy(pos)
    N_displacement[node]=N_displacement[node]+translate

	# Rotation
	curMoment = moment(node)
    angMom=angMom+curMoment*dt
    

    
	momentInertiaInverse=1.0 # todo ?? later change
    # orient=orient.*FromRotationVector(copy(angMom)*(dt*momentInertiaInverse))
    
    temp=FromRotationVector(copy(angMom)*(dt*momentInertiaInverse))
    
    orient=Quat(orient.w*temp.w,orient.x*temp.x,orient.y*temp.y,orient.z*temp.z)
    
	# orient.multiply(FromRotationVector(angMom.clone().multiplyScalar((dt*momentInertiaInverse)))) # tupdate the orientation //momentInertiaInverse

    N_orient[node]=copy(orient)
    
	
	eul = ToEulerAngles(orient) # TODO I THINK WRONG
    
    N_angle[node]=copy(eul)

    N_linMom[node]=copy(linMom)
    N_angMom[node]=copy(angMom)
    
    
	
	#  if (ext){//?? todo fix 
	#  	var size = 1;//change
	#  	if (ext->isFixed(X_TRANSLATE)) {pos.x = ix*size + ext->translation().x; linMom.x=0;}
	#  	if (ext->isFixed(Y_TRANSLATE)) {pos.y = iy*size + ext->translation().y; linMom.y=0;}
	#  	if (ext->isFixed(Z_TRANSLATE)) {pos.z = iz*size + ext->translation().z; linMom.z=0;}
	#  	if (ext->isFixedAnyRotation()){ //if any rotation fixed, all are fixed
	#  		if (ext->isFixedAllRotation()){
	#  			orient = ext->rotationQuat();
	#  			angMom = Vec3D<double>();
	#  		}
	#  		else { //partial fixes: slow!
	#  			Vec3D<double> tmpRotVec = orient.ToRotationVector();
	#  			if (ext->isFixed(X_ROTATE)){ tmpRotVec.x=0; angMom.x=0;}
	#  			if (ext->isFixed(Y_ROTATE)){ tmpRotVec.y=0; angMom.y=0;}
	#  			if (ext->isFixed(Z_ROTATE)){ tmpRotVec.z=0; angMom.z=0;}
	#  			orient.FromRotationVector(tmpRotVec);
	#  		}
	#  	}
	#  }

	#  poissonsStrainInvalid = true;
end

function force(node,static,currentTimeStep) 
	# forces from internal bonds
	totalForce=[0 0 0]
	# new THREE.Vector3(node.force.x,node.force.y,node.force.z);
	#  todo 
    
    
    
	totalForce=totalForce+N_intForce[node]

	#  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[node],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[node]
	#  }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(node) 
	#moments from internal bonds
	totalMoment=[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[node]
    
	totalMoment = RotateVec3D(N_orient[node],totalMoment);
    
    totalMoment=totalMoment+N_moment[node]

	#other moments
	# if (externalExists()) totalMoment += external()->moment(); //external moments
	# totalMoment -= angularVelocity()*mat->globalDampingRotateC(); //global damping
	
    N_intMoment[node]=[0 0 0] # do i really need it?
    
	return totalMoment
end

function updateDataAndSave(setup,fileName)
    nodes      = setup["nodes"]
    edges      = setup["edges"]
    voxCount=size(nodes)[1]
    linkCount=size(edges)[1]

    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][1]
        node["displacement"]["y"]=N_displacement[i][2]
        node["displacement"]["z"]=N_displacement[i][3]
        
        node["angle"]["x"]=N_angle[i][1]
        node["angle"]["y"]=N_angle[i][2]
        node["angle"]["z"]=N_angle[i][3]
        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


###############################################################################################
latticeSize=9
numTimeSteps=200
save=false

setup = Dict()
open("../json/setupTestUni9.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"]
############# nodes
N_position=[]
N_degrees_of_freedom=[]
N_restrained_degrees_of_freedom=[]
N_displacement=[]
N_angle=[]
N_currPosition=[]
N_linMom=[]
N_angMom=[]
N_intForce=[]
N_intMoment=[]
N_moment=[]
N_posTimeSteps=[]
N_angTimeSteps=[]
N_force=[]
N_orient=[]

############# edges
E_source=[]
E_target=[]
E_area=[]
E_density=[]
E_stiffness=[]
E_stress=[]
E_axis=[]
E_currentRestLength=[]
E_pos2=[]
E_angle1v=[]
E_angle2v=[]
E_angle1=[]
E_angle2=[]
E_currentTransverseStrainSum=[]
E_stressTimeSteps=[]

########
voxCount=0
linkCount=0
nodes      = setup["nodes"]
edges      = setup["edges"]
voxCount=size(nodes)[1]
linkCount=size(edges)[1]
strain =0 #todooo moveeee


######## ######## ######## ######## 
# println(toAxisXVector3([-5,0,-5 ],[-0.7071067811865475,0,-0.7071067811865475 ]))
# println(toAxisXVector3([5,-5,0  ],[0.7071067811865475,-0.7071067811865475,0  ]))
# println(toAxisXVector3([-5,5,0  ],[-0.7071067811865475,0.7071067811865475,0  ]))
# println(toAxisXVector3([-5,-5,0 ],[-0.7071067811865475,-0.7071067811865475,0 ]))
# println(toAxisXVector3([5,0,5   ],[0.7071067811865475,0,0.7071067811865475   ]))
# println(toAxisXVector3([0,5,5   ],[0,0.7071067811865475,0.7071067811865475   ]))

# println()
# println(toAxisOriginalVector3([1 0 0],[0 0 1]))
# println(toAxisOriginalVector3([0 1 0],[0 0 1]))
# println(toAxisOriginalVector3([0 0 1],[0 0 1]))

# println()

# println(toAxisXVector3([1 0 0],[0 0 1]))
# println(toAxisXVector3([0 1 0],[0 0 1]))
# println(toAxisXVector3([0 0 1],[0 0 1]))

# println()
dt=0.0251646

t=@timed simulateParallel(setup,numTimeSteps,dt,true,10)
time=t[2]
println("ran latticeSize $latticeSize with $voxCount voxels and $linkCount edges for $numTimeSteps time steps took $time seconds")
# setup["animation"]["exaggeration"]=75444

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