parallelFea.js

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

//BASED ON https://github.com/jonhiller/Voxelyze


var DBL_EPSILONx24 =5.328e-15; //DBL_EPSILON*24
var DISCARD_ANGLE_RAD= 1e-7; //Anything less than this angle can be considered 0
var SMALL_ANGLE_RAD= 1.732e-2; //Angles less than this get small angle approximations. To get: Root solve atan(t)/t-1+MAX_ERROR_PERCENT. From: MAX_ERROR_PERCENT = (t-atan(t))/t 
var SMALL_ANGLE_W =0.9999625; //quaternion W value corresponding to a SMALL_ANGLE_RAD. To calculate, cos(SMALL_ANGLE_RAD*0.5).
var W_THRESH_ACOS2SQRT= 0.9988; //Threshhold of w above which we can approximate acos(w) with sqrt(2-2w). To get: Root solve 1-sqrt(2-2wt)/acos(wt) - MAX_ERROR_PERCENT. From MAX_ERROR_PERCENT = (acos(wt)-sqrt(2-2wt))/acos(wt)
var X_AXIS = 0; 			//!< X Axis
var Y_AXIS = 1; 			//!< Y Axis
var Z_AXIS = 2; 		    //!< Z Axis
var currentTime=0;
var maxStrain=0;




function simulateParallel(setup,numTimeSteps,dt,static=true,saveInterval=10){
	// var instuctionsDiv=document.getElementById("footer2").innerHTML;
	initialize(setup);
    for(var i=0;i<numTimeSteps;i++){
        var t0 = performance.now();
        doTimeStep(setup,dt,static,i,saveInterval);
        var t1 = performance.now();
		console.log("TimeStep "+i+" took " + (t1 - t0) + " milliseconds.");
		
		// document.getElementById("footer2").innerHTML = "Timestep "+i +" out of "+ numTimeSteps+".";
	}
	// updateColors();
}

function initialize(setup){
	
	// pre-calculate current position
	var voxCount = setup.nodes.length;
	for(var i=0;i<voxCount;i++){
		setup.nodes[i].currPosition=new THREE.Vector3(setup.nodes[i].position.x,setup.nodes[i].position.y,setup.nodes[i].position.z);
		setup.nodes[i].orient= new THREE.Quaternion();
        setup.nodes[i].linMom=new THREE.Vector3(0,0,0);
        setup.nodes[i].angMom=new THREE.Vector3(0,0,0);
        setup.nodes[i].intForce=new THREE.Vector3(0,0,0);
        setup.nodes[i].intMoment=new THREE.Vector3(0,0,0);
		setup.nodes[i].moment={ x: 0, y: 0,z:0 };
		setup.nodes[i].displacement={ x: 0, y: 0,z:0 };
		//for dynamic simulations
		setup.nodes[i].posTimeSteps=[];
		setup.nodes[i].angTimeSteps=[];
		setup.nodes[i].nomSize=1.0;
		setup.nodes[i].massInverse=8e-6;
		setup.nodes[i].mass=1/8e-6;
		setup.nodes[i].FloorStaticFriction=false;
	} 
	// pre-calculate the axis
	var linkCount = setup.edges.length;
	for(var i=0;i<linkCount;i++){
		var node1=setup.nodes[setup.edges[i].source];
		var node2=setup.nodes[setup.edges[i].target];
		var pVNeg=new THREE.Vector3(node1.position.x,node1.position.y,node1.position.z);
		var pVPos=new THREE.Vector3(node2.position.x,node2.position.y,node2.position.z);
		var axis=pVPos.clone().sub(pVNeg).normalize();
		setup.edges[i].axis=axis.clone();


		setup.edges[i].currentRestLength=0;
        setup.edges[i].pos2= new THREE.Vector3(0,0,0);
        setup.edges[i].angle1v= new THREE.Vector3(0,0,0);
        setup.edges[i].angle2v= new THREE.Vector3(0,0,0);
        setup.edges[i].angle1=new THREE.Quaternion();
        setup.edges[i].angle2=new THREE.Quaternion();
        setup.edges[i].currentTransverseArea=0;
		setup.edges[i].currentTransverseStrainSum=0;
		//todo update stresses
		setup.edges[i].stressTimeSteps=[];
		
	} 
	
}

function doTimeStep(setup,dt,static=true,currentTimeStep,saveInterval){
	if (dt==0) 
		return true;
	else if (dt<0) 
		dt = recommendedTimeStep();

	// if (collisions) updateCollisions();
	

	var collisions=false;

	//Euler integration:
	var Diverged = false;
	var linkCount = setup.edges.length;
	for(var i=0;i<linkCount;i++){
        updateForces(setup,setup.edges[i],setup.nodes[setup.edges[i].source],setup.nodes[setup.edges[i].target],static);
        // todo: update forces and whatever
		if (axialStrain(setup.edges[i]) > 100) {
			Diverged = true; //catch divergent condition! (if any thread sets true we will fail, so don't need mutex...
		}

	}
    if (Diverged){
		console.log("Divergedd!!!!!")
		return false;

	}
	var voxCount = setup.nodes.length;

	for(var i=0;i<voxCount;i++){
		timeStep(dt,setup.nodes[i],static,currentTimeStep);
		if(!static&& currentTimeStep%saveInterval==0){
			setup.nodes[i].posTimeSteps.push(setup.nodes[i].displacement);
			setup.nodes[i].angTimeSteps.push(setup.nodes[i].angle);

		}
		// todo: update linMom,angMom, orient and whatever
	} 

	
	currentTime += dt;
	return true;

}

function updateForces(setup,edge,node1,node2,static=true){
	
	var pVNeg=new THREE.Vector3(node1.position.x,node1.position.y,node1.position.z);
	var pVPos=new THREE.Vector3(node2.position.x,node2.position.y,node2.position.z);
    var currentRestLength=pVPos.clone().sub(pVNeg).length();
	edge.currentRestLength=currentRestLength; //todo make sure updated
	

	pVNeg=node1.currPosition.clone();
	pVPos=node2.currPosition.clone();

	// Vec3D<double> three
	var oldPos2 = edge.pos2.clone();//??
	var oldAngle1v = edge.angle1v.clone();
	var oldAngle2v = edge.angle2v.clone(); //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,node1,node2); //sets pos2, angle1, angle2 /*restLength*/

	var dPos2=edge.pos2.clone().sub(oldPos2).multiplyScalar(0.5);
	var dAngle1=edge.angle1v.clone().sub(oldAngle1v).multiplyScalar(0.5);
	var dAngle2=edge.angle2v.clone().sub(oldAngle2v).multiplyScalar(0.5);

	
	
	
	//if volume effects...
    //if (!mat->isXyzIndependent() || currentTransverseStrainSum != 0) 
    //updateTransverseInfo(); //currentTransverseStrainSum != 0 catches when we disable poissons mid-simulation
	
	
	var _stress=updateStrain((edge.pos2.x/edge.currentRestLength),edge.stiffness);
	// var _stress=updateStrain(1.0);
	edge.stress = _stress;
	if(!static){
		edge.stressTimeSteps.push(_stress);
	}

	// console.log("Stress:"+edge.stress)
	if(setup.viz.minStress>edge.stress){
		setup.viz.minStress=edge.stress;
	}else if (setup.viz.maxStress<edge.stress){
		setup.viz.maxStress=edge.stress;
	}

	// 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
	
	var l   = currentRestLength;//??
	var rho = edge.density;
	var A = edge.area;
	var E = edge.stiffness;// youngs modulus
	var G=1.0;//todo shear_modulus
	var ixx = 1.0;//todo section ixx
	var I=1.0;
	var iyy = 1.0;//todo section.iyy//
	// var l0=length.dataSync();
	var J=1.0;//todo check
	// var l02 = l0 * l0;
	// var l03 = l0 * l0 * l0;
	var b1= 12*E*I/(l*l*l);
	var b2= 6*E*I/(l*l);
	var b3= 2*E*I/(l);
	var a1= E*A/l;
	var a2= G*J/l;
	var nu=0;

	// var b1= 5e6;
	// var b2= 1.25e7;
	// var b3= 2.08333e+07;
	// var a1= E*A/l;
	// var a2= 1.04167e+07;