Skip to content
Snippets Groups Projects
parallelFEA.jl 39.5 KiB
Newer Older
  • Learn to ignore specific revisions
  • 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