From 2fefc7272b18dda0adbfb071066aace2a9c7d814 Mon Sep 17 00:00:00 2001
From: amiraa <amiraa@mit.edu>
Date: Tue, 10 Sep 2024 16:57:42 -0400
Subject: [PATCH] CHI

---
 .../assembly/assembly.js                      |  11 +-
 .../demos/indexSimulationCHI.html             | 263 ++++++++++++++++++
 .../setup/globals.js                          |   5 +-
 .../simulation/app.js                         |  78 +++++-
 .../simulation/fea/beamFea.js                 | 141 ++++++++--
 .../simulation/visualization/geometry.js      |   3 +
 .../simulation/visualization/utils.js         | 151 +++++++++-
 .../threejs/grid.js                           |  13 +-
 8 files changed, 634 insertions(+), 31 deletions(-)
 create mode 100644 01_Code/physical_computing_interface/demos/indexSimulationCHI.html

diff --git a/01_Code/physical_computing_interface/assembly/assembly.js b/01_Code/physical_computing_interface/assembly/assembly.js
index 0337e2d..70b50fd 100644
--- a/01_Code/physical_computing_interface/assembly/assembly.js
+++ b/01_Code/physical_computing_interface/assembly/assembly.js
@@ -23,6 +23,7 @@ function Assembler(three,GLOBALS,numberOfRobots,speed,robotLocations,depositLoca
 	this.color2=GLOBALS.color2;//kohly
 	this.color3=GLOBALS.color3;///*teal*/
 	this.color4= GLOBALS.color4; //red/orange
+	this.viewRobotArm=GLOBALS.viewRobotArm;
 
 	///////////////////////////voxels///////////////////////////////
 	this.voxel;
@@ -442,7 +443,11 @@ Assembler.prototype.defaultRobot=function(robotIndex) {
 	  };
 	  
 	this.THREESimulationRobot[robotIndex] = new THREE.Group();
+	
+	this.THREESimulationRobot[robotIndex].visible=this.viewRobotArm;
+	
 	this.scene.add(this.THREESimulationRobot[robotIndex]);
+	
 
 	  ////////////////
 
@@ -643,6 +648,7 @@ Assembler.prototype.getNormalAdjustment=function(robotIndex,n,vnormal,forward){/
 ////////////////////////Taget Control////////////////////////////
 Assembler.prototype.targetControl=function(robotIndex){
 	this.target[robotIndex] = new THREE.Group();
+	this.target[robotIndex].visible=this.viewRobotArm;
 	this.scene.add(this.target[robotIndex]);
 
 	this.control[robotIndex] = new THREE.TransformControls(this.camera, this.renderer.domElement);
@@ -673,6 +679,7 @@ Assembler.prototype.targetControl=function(robotIndex){
 	});
 	this.control[robotIndex].attach(this.target[robotIndex]);
 
+	this.control[robotIndex].visible=this.viewRobotArm;
 	this.scene.add(this.control[robotIndex]);
 	// control[robotIndex].visible = false;
 }
@@ -1151,7 +1158,7 @@ Assembler.prototype.buildVoxelAt=function(loc){
 	// }
 	
 
-
+	
 	this.scene.add( object1 );
 }
 //////////////////////////path planing///////////////////////////
@@ -1230,6 +1237,7 @@ Assembler.prototype.buildHelperMeshes=function(robotIndex){
 	this.targetPositionMesh[robotIndex].scale.x=0.8;
 	this.targetPositionMesh[robotIndex].scale.y=0.8;
 	this.targetPositionMesh[robotIndex].scale.z=0.8;
+	// this.targetPositionMesh[robotIndex].visible=this.viewRobotArm;
 	this.scene.add(this.targetPositionMesh[robotIndex]);
 
 	for (var count=0; count<this.numberOfRobots;count++)
@@ -1239,6 +1247,7 @@ Assembler.prototype.buildHelperMeshes=function(robotIndex){
 		mesh.position.x=this.startLocations[count].x;
 		mesh.position.y=this.startLocations[count].y;
 		mesh.position.z=this.startLocations[count].z;
+		mesh.visible=this.viewRobotArm;
 		this.scene.add(mesh);
 	}
 
diff --git a/01_Code/physical_computing_interface/demos/indexSimulationCHI.html b/01_Code/physical_computing_interface/demos/indexSimulationCHI.html
new file mode 100644
index 0000000..de4ab6d
--- /dev/null
+++ b/01_Code/physical_computing_interface/demos/indexSimulationCHI.html
@@ -0,0 +1,263 @@
+<html>
+
+<head>
+    <title>Physical Computing Interface</title>
+    <link rel="stylesheet" type="text/css" href="../css/style.css" media="screen"/>
+    <link rel="stylesheet" type="text/css" href="../lib/jsoneditor/jsoneditor.css" >
+
+    <!-- assembler control css -->
+    <link rel="stylesheet" href="../css/bootstrap.min.css">
+    <link rel="stylesheet" type="text/css" href="../css/styleAssemblerControl.css" media="screen"/>
+
+    <!-- <link href="https://unpkg.com/font-awesome@5.8.0/css/font-awesome.min.css" rel="stylesheet" type="text/css" /> -->
+    <!-- <link href="//netdna.bootstrapcdn.com/font-awesome/3.2.1/css/font-awesome.css" rel="stylesheet"> -->
+    <script src="https://kit.fontawesome.com/99c302ff33.js" crossorigin="anonymous"></script>
+    <!-- <link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/4.7.0/css/font-awesome.min.css"> -->
+
+
+</head>
+        
+<body>
+
+    <div id="threejs">
+        <div id="threejs1">
+            <div class="header1">
+                    <i> Assembly</i>
+                    <!-- Three.js -->
+            </div>
+            
+            <div id="webgl"></div>
+        </div>
+        <div class="slidecontainer">
+                <input type="range" min="0" max="0" value="1" class="slider" id="time">
+        </div>
+            
+    </div>
+
+    <div id="simulation">
+        <div id="threejs1">
+            <div class="header2">
+                    
+                    <i> Simulation     </i>
+
+                    <!-- <div  class="dropdown" style="float: right;" >
+                        <button onclick="simulationFidelity()" class="dropbtn" >Fidelity</button>
+                        <div id="myDropdown" class="dropdown-content" >
+                          <a class="High" href="#">High (slow)</a>
+                          <a class="Medium" href="#">Medium</a>
+                          <a class="Low" href="#">Low (Fast)</a>
+                        </div>
+                    </div > -->
+            </div>
+            <div class="dragbar2"></div>
+            <div id="webgl1"></div>
+            
+            
+        </div>
+            
+    </div>
+    
+
+
+    <div id="graph">
+            <div class="header2">
+                    <div class="dragbar"></div> 
+                <i> Graph</i>
+            </div>
+            <div id=jsondiveditor>
+                
+                <div id="cy"></div>
+            </div>
+            <div class="dragbar"></div> 
+            
+            
+    </div>
+
+    <div id="json">
+            
+
+            <div class="header2">
+                    <div class="dragbar"></div> 
+                <i> Node</i>
+            </div>
+            
+            
+            <div id=jsondiveditor>
+                    
+                <br></br>
+                    <!-- <p>
+                        <button class="button" id="setJSON">Get Info</button>
+                        <button class="button" id="getJSON">Set Info</button>
+                    </p> -->
+                
+                    <div id="jsoneditor"></div>
+            </div>
+            <div class="dragbar"></div> 
+            <div class="dragbar1"></div>
+            
+            
+            
+    </div>
+
+    <div class="footer1" id="instructions0">
+        <strong>left-click</strong>: place voxel/orbit, <strong>right-click</strong>: radial menu with more options
+    </div>
+
+    <div class="footer2" id="instructions">
+        <!-- update change to more instructions/feedback -->
+        <i>instructions</i>
+    </div>
+
+</body>
+
+
+
+<!-- TODO: 
+            Clean structure to modules?
+            Add another footer
+             
+-->
+
+
+<!-- libraries -->
+<script src="https://ajax.googleapis.com/ajax/libs/jquery/3.4.1/jquery.min.js"></script>
+<script src="../lib/cytoscape.min.js"></script>
+<script src="../lib/cytoscape-cxtmenu1.js"></script>
+<script src="https://unpkg.com/layout-base/layout-base.js"></script>
+<script src="https://unpkg.com/cose-base/cose-base.js"></script>
+<script src="https://unpkg.com/cytoscape-cose-bilkent/cytoscape-cose-bilkent.js"></script>
+<script src="../lib/cytoscape-expand-collapse.js"></script>
+
+
+<script src="../lib/jsoneditor/jsoneditor.js"></script>
+
+<script src="../lib/opencv.js"></script>
+
+
+<script src="../lib/three.min.js"></script>
+<script src="../lib/OrbitControls.js"></script>
+<script src="../lib/dat.gui.min.js"></script>
+<script src="../lib/TransformControls.js"></script>
+
+<script src="../assembly/InverseKinematic.js"></script><!-- TODO LOCATION -->
+<script src="../assembly/voxel.js"></script><!-- TODO CHANGE TO DICE PIECES -->
+
+<!-- simulation stuff -->
+<script src="https://ajax.googleapis.com/ajax/libs/jquery/3.4.1/jquery.min.js"></script>
+
+<script type="text/javascript" src="https://files.mcneel.com/rhino3dm/js/latest/rhino3dm.js"></script>
+<!-- <script type="text/javascript" src="../simulation/demos/rhino3dm.js"></script> -->
+
+
+<!-- <script src="https://cdn.jsdelivr.net/npm/@tensorflow/tfjs@1.0.0/dist/tf.min.js"></script> -->
+<script src="https://cdn.jsdelivr.net/npm/@tensorflow/tfjs@4.15.0/dist/tf.min.js"></script>
+
+<!-- simulation -->
+<script src="../simulation/lib/stats.js"></script>
+<script src="../simulation/lib/SVGRenderer.js"></script>
+<script src="../simulation/lib/js-colormaps.js"></script>
+<script src="../simulation/lib/LineSegments2.js"></script>
+<script src="../simulation/lib/LineSegmentsGeometry.js"></script>
+<script src="../simulation/lib/Line2.js"></script>
+<script src="../simulation/lib/LineMaterial.js"></script>
+<script src="../simulation/lib/LineGeometry.js"></script>
+<script src="../simulation/lib/GeometryUtils.js"></script>
+
+
+<!-- <script src="./metavoxels/01_Code/191115_NodeJsJulia/fea/beamFea.js"></script>  -->
+<!-- <script src="./simulation/visualization/geometry.js"></script>  -->
+
+<script src="../simulation/visualization/utils.js"></script> 
+<script src="../simulation/visualization/geometry.js"></script> 
+
+<script src="../simulation/visualization/main.js"></script> 
+
+<script src="../simulation/fea/parallelFea.js"></script> 
+<script src="../simulation/fea/beamFea.js"></script> 
+<!-- assembler control -->
+<!-- <script src="./assembly/assemblerControl/utils.js" type="text/javascript"></script>
+<script src="./assembly/assemblerControl/setup.js" type="text/javascript"></script>
+<script src="./assembly/assemblerControl/processing.js" type="text/javascript"></script>
+<script src="./assembly/assemblerControl/serialMonitor.js" type="text/javascript"></script>
+<script type="text/javascript" src="./assembly/assemblerControl/bootstrap-multiselect.js"></script> -->
+<script>
+    var i = 0;
+    
+     
+    var viewRobotArm=false;
+    
+    var startPos=0.8
+    $('#threejs').css("width",window.innerWidth*startPos);
+    $('#simulation').css("width",window.innerWidth*startPos);
+    $('#graph').css("left",window.innerWidth*startPos);
+    $('#json').css("left",window.innerWidth*startPos);
+    $('.footer1').css("width",window.innerWidth*startPos);
+    $('.footer2').css("left",window.innerWidth*startPos);
+    // onWindowResize();
+</script>
+
+<!-- code -->
+<script src="../setup/globals.js"></script> <!-- event handling and GLOBALS,UTILS -->
+<script src="../threejs/grid.js"></script><!-- threejs visualization -->
+<script src="../assembly/assembly.js"></script><!-- robot assembly -->
+<script src="../assembly/replay.js"></script><!-- assembly and timestep handling -->
+<script src="../graph/graph.js"></script><!-- graph flow visualization-->
+<script src="../json/json.js"></script><!-- json -->
+
+<script src="../simulation/app.js"></script><!-- threejs visualization -->
+
+<!-- windows control -->
+<script>
+    
+
+   $('.dragbar').mousedown(function(e){
+       
+        e.preventDefault();
+        $(document).mousemove(function(e){
+        //   $('#position').html(e.pageX +', '+ e.pageY);
+          $('#threejs').css("width",e.pageX+2);
+
+          $('#simulation').css("width",e.pageX+2);
+
+          $('#graph').css("left",e.pageX+2);
+          $('#json').css("left",e.pageX+2);
+          $('.footer1').css("width",e.pageX+2);
+          $('.footer2').css("left",e.pageX+2);
+          
+       })
+    //    onWindowResize();//todo change location
+    });
+
+    $('.dragbar1').mousedown(function(e){
+       e.preventDefault();
+       // $('#mousestatus').html("mousedown" + i++);
+       $(document).mousemove(function(e){
+        $('#graph').css("height",e.pageY+2);
+        $('#json').css("top",e.pageY+2);
+         
+      })
+   });
+
+   $('.dragbar2').mousedown(function(e){
+       e.preventDefault();
+       // $('#mousestatus').html("mousedown" + i++);
+       $(document).mousemove(function(e){
+        $('#threejs').css("height",e.pageY+2);
+        $('#simulation').css("top",e.pageY+2);
+         
+      })
+    //   onWindowResize();//todo change location
+   });
+
+   $(document).mouseup(function(e){
+       $(document).unbind('mousemove');
+    });
+
+    
+
+</script>
+
+
+
+
+</html>
\ No newline at end of file
diff --git a/01_Code/physical_computing_interface/setup/globals.js b/01_Code/physical_computing_interface/setup/globals.js
index 56150d8..b761c8c 100644
--- a/01_Code/physical_computing_interface/setup/globals.js
+++ b/01_Code/physical_computing_interface/setup/globals.js
@@ -1,6 +1,8 @@
 // Amira Abdel-Rahman
 // (c) Massachusetts Institute of Technology 2019
-
+if (viewRobotArm === undefined) {
+    var viewRobotArm=true;
+}
 function globals(utils){
     this.gridSize=20;
     this.voxelSpacing=50;
@@ -11,6 +13,7 @@ function globals(utils){
     this.color5=0x380152; //purple
     this.color6=0x696767; //grey
     this.color7=0x42f2f5; //light blue
+    this.viewRobotArm=viewRobotArm;
     this.gridPresets={
         Cubic :  {
         
diff --git a/01_Code/physical_computing_interface/simulation/app.js b/01_Code/physical_computing_interface/simulation/app.js
index 32dc56a..6f64774 100644
--- a/01_Code/physical_computing_interface/simulation/app.js
+++ b/01_Code/physical_computing_interface/simulation/app.js
@@ -72,7 +72,7 @@ var static=true;
 var latticeSize=10;
 var voxelSize=5;
 var globalSetup={
-    exaggeration:10e1,
+    exaggeration:1,
     speed:3.0,
     updateStress:false
 
@@ -84,9 +84,10 @@ setup.voxelSize=voxelSize;
 setup.gridSize=GLOBALS.gridSize
 
 document.addEventListener('runNode', function (e) { 
-    x = document.getElementById("instructions");x.innerHTML = "Simulation Started, wait a few seconds for the results!";
+    // x = document.getElementById("instructions");x.innerHTML = "Simulation Started, wait a few seconds for the results!";
 
-    solveParallel(GLOBALS.globalJson.simulationSteps);
+    // solveParallel(GLOBALS.globalJson.simulationSteps);
+    solveFEA();
 
 }, false);
 
@@ -100,7 +101,7 @@ document.addEventListener('changeNodeColor', function (e) {
              voxelSize/2+voxelSize*e.detail.z,
              voxelSize/2+voxelSize*e.detail.y]);
         console.log(load1)
-        loads.push([load1,{x:0,y:-100*5,z:0}]);
+        loads.push([load1,{x:0,y:-100,z:0}]);
     }
     resetMetavoxels(e);
 }, false);
@@ -152,6 +153,8 @@ function resetMetavoxels(e){
     
 
     sim.update(setup);
+    setup.animation.exaggeration=1.0;
+
 
 
     
@@ -169,6 +172,12 @@ rhino3dm().then(async m => {
 		density:0.028,
 		stiffness:10000000
     };
+
+    material={
+		area:1.6*6,
+		density:0.028,
+		stiffness:10000
+    };
     
     var material2={
 		area:1.0,
@@ -201,6 +210,7 @@ rhino3dm().then(async m => {
     var dof=[true,true,true,true,true,true];
     supports=[[support,dof]];
     loads=[[load,{x:0,y:-100,z:0}]];
+    loads=[[load,{x:0,y:-1,z:0}]];
     // var diffMaterialBox=[[matB,material2]];
     setup.loads=loads
     setup.supports=supports
@@ -219,10 +229,14 @@ rhino3dm().then(async m => {
     // restrainFromBox(setup,supports);
 
     setup.viz.colorMaps=[YlGnBu,coolwarm, winter ,jet];
+    // setup.viz.colorMaps=[coolwarm,YlGnBu, winter ,jet];
     setup.viz.minStress=10e6;
     setup.viz.maxStress=-10e6;
+    
+   
 
     setup.viz.exaggeration=globalSetup.exaggeration;
+    
     setup.solve=solveParallel;
     setup.numTimeSteps=100;
     // setup.supports=supports;
@@ -231,6 +245,7 @@ rhino3dm().then(async m => {
    
     sim= new threejs1(setup,"webgl1","graph",static);
     sim.init();
+    
     // three.drawConstraintBoundingBoxes([support],[load]);
    
 
@@ -272,3 +287,58 @@ function solveParallel(numTimeSteps){
     /////////////////////
 
 }
+
+function solveFEA(){
+
+    // const a = tf.tensor2d([1, 3, 3, 1, 4, 3, 1, 3, 4], [3, 3])
+    // const b = tf.tensor2d([1, 3, 3], [3, 1])
+    // console.log("matrix a")
+    // a.print()
+    // console.log("matrix d")
+    // var d=det11(a)
+    // console.log(d)
+    // var d=det(a)
+    // console.log(d)
+    // // d.print()
+
+    // const i = invertMatrix(a)
+    // console.log("inverse i")
+    // i.print()
+    // const prod = i.dot(b)
+    // console.log("b * i")
+    // prod.print()
+
+    // var displacements=lusolve(a.arraySync(), b.arraySync(),false);
+    // console.log("b * i")
+    // console.log(displacements)
+
+    // console.log(setup.viz.colorMap)
+    // console.log(sim.setup.viz.colorMap)
+    sim.setup.viz.colorMap=1; //coolwarm
+    setup.animation.exaggeration=10;
+
+    console.log("exaggeration "+setup.animation.exaggeration)
+    var ins = document.getElementById("instructions0");
+	ins.innerHTML = "Simulation running, wait a few seconds for the results!";
+    setTimeout(function() { 
+        var ins = document.getElementById("instructions0");
+	    ins.innerHTML = "Simulation running, wait a few seconds for the results!";
+    }, 1);
+
+    setTimeout(function() { 
+        solveFea(setup,2)
+        // simulateParallel( setup,numTimeSteps,dt,static,2);
+    }, 2);
+
+    setTimeout(function() { 
+        var alumnimYieldStress=120;
+		var ins = document.getElementById("instructions0");
+        var txt = "Done! Displacement min:"+setup.viz.minDisplacement.toFixed(3)+" , max:"+ setup.viz.maxDisplacement.toFixed(3)+" mm.";
+        var txt1 = " Stress min:"+setup.viz.minStress.toFixed(3)+" , max:"+ setup.viz.maxStress.toFixed(3)+" MPa.";
+        var safety=alumnimYieldStress/Math.max(Math.abs(setup.viz.minStress),Math.abs(setup.viz.maxStress))
+        var txt2 = " Safety factor is :"+safety.toFixed(2)+".";
+		ins.innerHTML = txt+txt1+txt2;
+	
+	}, 3);
+
+}
diff --git a/01_Code/physical_computing_interface/simulation/fea/beamFea.js b/01_Code/physical_computing_interface/simulation/fea/beamFea.js
index dfab386..4141afa 100644
--- a/01_Code/physical_computing_interface/simulation/fea/beamFea.js
+++ b/01_Code/physical_computing_interface/simulation/fea/beamFea.js
@@ -2,7 +2,9 @@
 // (c) Massachusetts Institute of Technology 2019
 
 ////////////3D bar-truss fea/////////////////////
-function solveFea(){
+function solveFea(setup,method){
+
+	
 	// # determine the global matrices
 	var M,K,F;
 	var t0 = performance.now();
@@ -18,24 +20,56 @@ function solveFea(){
 	// M, K, F = get_matrices(properties)
 	
 	// [Q,R]=tf.linalg.qr(K);
-	t0 = performance.now();
-	var displacements=lusolve(K.arraySync(), F.arraySync(),false);
-	t1 = performance.now();
-	console.log("Call to invert " + (t1 - t0) + " milliseconds.");
-	// console.log(displacements);
+	if (method==0){ //not working todo fix
+		t0 = performance.now();
+		console.log(K.shape)
+		K.print();
+		F.print();
+		// var d = det(K);
+		// console.log(d);
+		// d.print()
+		var inv = invertM(K);
+		console.log("hena");
+		var prod = inv.dot(F);
+		var displacements=prod.arraySync();
+		t1 = performance.now();
+		console.log("Call to invert jordan-gauss " + (t1 - t0) + " milliseconds.");
+
+
+	}else if( (method==1)){ //not working todo fix
+		t0 = performance.now();
+		var inv = invertMatrix(K)
+		var prod = inv.dot(F)
+		var displacements=prod.arraySync();
+		t1 = performance.now();
+		console.log("Call to invert matrix adjoint method " + (t1 - t0) + " milliseconds.");
+
+
+	}else if( (method==2)){
+		t0 = performance.now();
+		var displacements=lusolve(K.arraySync(), F.arraySync(),false);
+		t1 = performance.now();
+		console.log("Call to invert (math.js) " + (t1 - t0) + " milliseconds.");
+		// console.log(displacements);
+	}
+	
+	
+
+	
 
 	updateDisplacement(displacements);
 	var X= tf.tensor(displacements);//change later to more effecient solver
 	console.log("Displacements:");
-	X.print();
+	// X.print();
 
 	// # determine the stresses in each element
 	var stresses = get_stresses(setup, displacements);
-	updateStresses(stresses.arraySync());
-
 	console.log("Stresses:");
-	stresses.print();
+	// stresses.print();
+	updateStresses(stresses.arraySync());
 
+	
+	
 	// console.log(setup);
 
 	// console.log(setup);
@@ -64,11 +98,27 @@ function getMatrices(setup){
 		var rho = tf.scalar(element.density);
 		var A = element.area;
 		var E = element.stiffness;// youngs modulus
-		var G=1.0;//todo shear_modulus
-		var ixx = 1.0;//todo section ixx
-		var iyy = 1.0;//todo section.iyy//
+
+
+		var b=Math.sqrt(A);
+		var h=b;
+
+		var ixx=b*h**3/12;//rectangle filled
+		var iyy=h*b**3/12; //rectangle filled
+
+		var G = E * 1 / 3;
+		var j=b*h*(b*b+h*h)/12;//rectangle filled #polar moment of inertia
+
+
+
+		// var G=1.0;//todo shear_modulus
+		// var ixx = 1.0;//todo section ixx
+		// var iyy = 1.0;//todo section.iyy//
+		// var j=1.0;//todo check
+
+
 		var l0=length.dataSync();
-		var j=1.0;//todo check
+		
 		var l02 = l0 * l0;
         var l03 = l0 * l0 * l0;
 
@@ -94,12 +144,14 @@ function getMatrices(setup){
             [0, 0, -6*E*iyy/l02, 0, 2*E*iyy/l0, 0, 0, 0, 6*E*iyy/l02, 0, 4*E*iyy/l0, 0],
             [0, 6*E*ixx/l02, 0, 0, 0, 2*E*ixx/l0, 0, -6*E*ixx/l02, 0, 0, 0, 4*E*ixx/l0]
         ]);
+		
 
 		var x_axis=tf.tensor ([1,0,0]);
 		var y_axis=tf.tensor ([0,1,0]);
 		var z_axis=tf.tensor ([0,0,1]);
 
 		// # find rotated mass and stifness matrices
+		var tau = rotation_matrix_old(element_vector, x_axis,y_axis,z_axis);
 		var tau = rotation_matrix(element_vector, x_axis,y_axis,z_axis);
 		// var m_r=tf.matMul(tau,tf.matMul(m,tau) ,true ,false) //tau.T.dot(m).dot(tau)
 		var k_r=tf.matMul(tau,tf.matMul(k,tau) ,true ,false) //tau.T.dot(k).dot(tau)
@@ -188,6 +240,7 @@ function get_stresses(setup,X){
 		var z_axis=tf.tensor ([0,0,1]);
 
 		// # find rotated mass and stifness matrices
+		var tau = rotation_matrix_old(element_vector, x_axis,y_axis,z_axis);
 		var tau = rotation_matrix(element_vector, x_axis,y_axis,z_axis);
 
 		var global_displacements=tf.tensor([setup.nodes[element.source].displacement.x,setup.nodes[element.source].displacement.y,setup.nodes[element.source].displacement.z,
@@ -214,7 +267,7 @@ function get_stresses(setup,X){
 
 ///////utils////////
 
-function rotation_matrix(element_vector,x_axis, y_axis, z_axis){ //todo check again 
+function rotation_matrix_old(element_vector,x_axis, y_axis, z_axis){ //todo check again 
 
 	var L=element_vector.norm().arraySync();
 	element_vector=element_vector.arraySync();
@@ -268,10 +321,66 @@ function rotation_matrix(element_vector,x_axis, y_axis, z_axis){ //todo check ag
 	// transMatrix[3:6, 3:6] = dirCos
 	// transMatrix[6:9, 6:9] = dirCos
 	// transMatrix[9:12, 9:12] = dirCos
-	
+	// console.log(transMatrix)
+	// transMatrix.print()
 	return transMatrix;
 }
 
+
+
 function direction_cosine(vec1, vec2){
 	return tf.div(tf.dot(vec1,vec2) ,tf.mul( vec1.norm() , vec2.norm()));
+}
+
+function rotation_matrix(dx,x_axis, y_axis, z_axis){
+	var l0 = dx.norm().arraySync();
+	element_vector=dx.arraySync();
+    // c = dx / l0;
+
+	var cxx = element_vector[0]/ l0;
+    var cyx = element_vector[1]/ l0;
+    var czx = element_vector[2]/ l0;
+    var d = Math.sqrt(cxx * cxx + cyx * cyx + czx * czx); //todo check if correct later
+	//try to use this to fix calcAelem https://github.com/latture/threed-beam-fea/blob/master/src/threed_beam_fea.cpp
+	// var d = Math.sqrt(cxx * cxx + cyx * cyx);
+    var cxy = -cyx / d;
+    var cyy = cxx / d;
+    var czy = 0;
+    var cxz = -cxx * czx / d;
+    var cyz = -cyx * czx / d;
+    var czz = d;
+
+	// if (isNaN(cxy)) {
+	// 	var dd=direction_cosine(dx, x_axis)
+	// 	dd.print()
+	// 	var dd=direction_cosine(dx, y_axis)
+	// 	dd.print()
+	// 	var dd=direction_cosine(dx, z_axis)
+	// 	dd.print()
+	// 	console.log(d)
+	// 	console.log(element_vector)
+	// 	var gamma=tf.tensor([
+	// 		[cxx,cyx,czx],
+	// 		[cxy,cyy,czy],
+	// 		[cxz,cyz,czz]])
+	// 	gamma.print()
+	// }
+
+	var T=tf.tensor([
+		[cxx,cyx,czx,  0,  0,  0,  0,  0,  0,  0,  0,  0],
+		[cxy,cyy,czy,  0,  0,  0,  0,  0,  0,  0,  0,  0],
+		[cxz,cyz,czz,  0,  0,  0,  0,  0,  0,  0,  0,  0],
+		[  0,  0,  0,cxx,cyx,czx,  0,  0,  0,  0,  0,  0],
+		[  0,  0,  0,cxy,cyy,czy,  0,  0,  0,  0,  0,  0],
+		[  0,  0,  0,cxz,cyz,czz,  0,  0,  0,  0,  0,  0],
+		[  0,  0,  0,  0,  0,  0,cxx,cyx,czx,  0,  0,  0],
+		[  0,  0,  0,  0,  0,  0,cxy,cyy,czy,  0,  0,  0],
+		[  0,  0,  0,  0,  0,  0,cxz,cyz,czz,  0,  0,  0],
+		[  0,  0,  0,  0,  0,  0,  0,  0,  0,cxx,cyx,czx],
+		[  0,  0,  0,  0,  0,  0,  0,  0,  0,cxy,cyy,czy],
+		[  0,  0,  0,  0,  0,  0,  0,  0,  0,cxz,cyz,czz]
+		]);
+	// T.print()
+    return T
+
 }
\ No newline at end of file
diff --git a/01_Code/physical_computing_interface/simulation/visualization/geometry.js b/01_Code/physical_computing_interface/simulation/visualization/geometry.js
index 72e36cd..ac31d6a 100644
--- a/01_Code/physical_computing_interface/simulation/visualization/geometry.js
+++ b/01_Code/physical_computing_interface/simulation/visualization/geometry.js
@@ -649,6 +649,9 @@ function linkNodes(setup,nodesIndices,area, density, stiffness){ //create square
 
     }
     add_Edge(setup,nodesIndices[nodesIndices.length-1], nodesIndices[0], area, density , stiffness);
+	//chi voxel
+	add_Edge(setup,nodesIndices[0], nodesIndices[2], area, density , stiffness);
+	add_Edge(setup,nodesIndices[1], nodesIndices[3], area, density , stiffness);
 }
 
 function restrain(setup,name){
diff --git a/01_Code/physical_computing_interface/simulation/visualization/utils.js b/01_Code/physical_computing_interface/simulation/visualization/utils.js
index d86fdb2..6868964 100644
--- a/01_Code/physical_computing_interface/simulation/visualization/utils.js
+++ b/01_Code/physical_computing_interface/simulation/visualization/utils.js
@@ -93,6 +93,10 @@ function edgeNeeded(setup,source,target){
 
 //////////////
 function updateDisplacement(X){
+
+	setup.viz.minDisplacement=Math.min(...X);
+	setup.viz.maxDisplacement=Math.max(...X);
+
 	var count=0;
 	for(var i=0;i<setup.nodes.length;i++){
 		if(!setup.nodes[i].restrained_degrees_of_freedom[0]){
@@ -120,17 +124,19 @@ function updateDisplacement(X){
 
 function updateStresses(S){
 	
-	setup.viz.minStress=Math.min(...S)*1.1;
-	setup.viz.maxStress=Math.max(...S)*1.1;
+	setup.viz.minStress=Math.min(...S);
+	setup.viz.maxStress=Math.max(...S);
 
 	var count=0;
 	for(var ii=0;ii<setup.edges.length;ii++){
 		var element=setup.edges[ii];
 		element.stress=S[ii];
 	}
-	if(!node){
-		three.colorEdges(); //todo check!!
-	}
+	updateColors()
+	// three.colorEdges(); //todo check!!
+	// if(!node){
+	// 	three.colorEdges(); //todo check!!
+	// }
 	
 }
 
@@ -243,4 +249,139 @@ function lubksb(lu, b, update) {
 		b[i] = sum / A[i][i]
 	}
 	return b // solution vector x
+}
+
+/////////////////////////
+
+//from https://stackoverflow.com/questions/51805389/how-do-i-invert-a-matrix-in-tensorflow-js
+
+
+// calculate the determinant of a matrix m
+function det11(m) {
+    return tf.tidy(() => {
+        const [r, _] = m.shape
+        if (r === 2) {
+            const t = m.as1D()
+            const a = t.slice([0], [1]).dataSync()[0]
+            const b = t.slice([1], [1]).dataSync()[0]
+            const c = t.slice([2], [1]).dataSync()[0]
+            const d = t.slice([3], [1]).dataSync()[0]
+            result = a * d - b * c
+            return result
+
+        } else {
+            let s = 0;
+            rows = [...Array(r).keys()]
+            for (let i = 0; i < r; i++) {
+                sub_m = m.gather(tf.tensor1d(rows.filter(e => e !== i), 'int32'))
+                sli = sub_m.slice([0, 1], [r - 1, r - 1])
+                s += Math.pow(-1, i) * det11(sli)
+            }
+            return s
+        }
+    })
+}
+
+function det( /*tf.Tensor*/a ) {
+	const a_shape = a.shape,
+			[m,n] = a_shape.slice(-2);
+	if( m !== n )
+	  throw new Error('det(): Matrix/matrices must be square.');
+  
+	return tf.tidy(() => {
+		return tf.stack(
+		a.reshape([-1,n,n]).unstack().map( a => {
+			const [q,r] = tf.linalg.qr(a),
+				diag_r = r.flatten().stridedSlice([0],[n*n],[n+1]),
+				det_r = diag_r.prod(),
+				det_q = n%2===0 ? +1 : -1; // <- q is product of n Householder
+											//    reflections, i.e. det(q) = (-1)^n
+			return tf.mul(det_q, det_r);
+		})
+		).reshape( a_shape.slice(0,-2) ).arraySync();
+	})
+}
+
+// the inverse of the matrix : jordan-gauss method
+function invertM(m) {
+    return tf.tidy(() => {
+        if (det(m) === 0) {
+            return
+        }
+        const [r, _] = m.shape
+        let inv = m.concat(tf.eye(r), 1)
+        rows = [...Array(r).keys()]
+        for (let i = 0; i < r; i++) {
+            inv = tf.tidy(() => {
+                for (let j = i + 1; j < r; j++) {
+                    const elt = inv.slice([j, i], [1, 1]).as1D().asScalar()
+                    const pivot = inv.slice([i, i], [1, 1]).as1D().asScalar()
+                    let newrow
+                    if (elt.dataSync()[0] !== 0) {
+                        newrow = inv.gather(tf.tensor1d([i], 'int32')).sub(inv.gather(tf.tensor1d([j], 'int32')).div(elt).mul(pivot)).as1D()
+                        const sli = inv.gather(rows.filter(e => e !== j))
+                        const arr = []
+                        if (j === 0) {
+                            arr.push(newrow)
+                        }
+                        sli.unstack().forEach((t, ind) => {
+                            if (ind !== j) {
+                                arr.push(t)
+                            } else {
+                                arr.push(newrow)
+                                arr.push(t)
+                            }
+                        })
+                        if (j === r - 1) {
+                            arr.push(newrow)
+                        }
+                        inv = tf.stack(arr)
+                    }
+                }
+                return inv
+            })
+        }
+        const trian = tf.unstack(inv)
+        len = trian.length
+        trian[len - 1] = trian[len - 1].div(trian[len - 1].slice(trian[len - 1].shape[0] - 1, 1).asScalar())
+        for (let i = r - 2; i > -1; i--) {
+            for (j = r - 1; j > i; j--) {
+                trian[i] = trian[i].sub(trian[j].mul(trian[i].slice(j, 1).asScalar()))
+            }
+            trian[i] = trian[i].div(trian[i].slice(i, 1).asScalar())
+        }
+        return tf.split(tf.stack(trian), 2, 1)[1]
+    })
+}
+
+// the inverse of the matrix : matrix adjoint method
+function invertMatrix(m) {
+	return tf.tidy(() => {
+		const d = det(m)
+		if (d === 0) {
+			return
+		}
+		[r, _] = m.shape
+		rows = [...Array(r).keys()]
+		dets = [];
+		for (let i = 0; i < r; i++) {
+			for (let j = 0; j < r; j++) {
+				const sub_m = m.gather(tf.tensor1d(rows.filter(e => e !== i), 'int32'))
+				let sli
+				if (j === 0) {
+					sli = sub_m.slice([0, 1], [r - 1, r - 1])
+				} else if (j === r - 1) {
+					sli = sub_m.slice([0, 0], [r - 1, r - 1])
+				} else {
+					const [a, b, c] = tf.split(sub_m, [j, 1, r - (j + 1)], 1)
+					sli = tf.concat([a, c], 1)
+				}
+				dets.push(Math.pow(-1, (i + j)) * det(sli))
+			}
+		}
+		com = tf.tensor2d(dets, [r, r])
+		tr_com = com.transpose()
+		inv_m = tr_com.div(tf.scalar(d))
+		return inv_m
+	})
 }
\ No newline at end of file
diff --git a/01_Code/physical_computing_interface/threejs/grid.js b/01_Code/physical_computing_interface/threejs/grid.js
index f4b0679..9cfba02 100644
--- a/01_Code/physical_computing_interface/threejs/grid.js
+++ b/01_Code/physical_computing_interface/threejs/grid.js
@@ -114,6 +114,7 @@ function threejs(GLOBALS,utils,containerName,container1Name){
                         type:"force",
                         locked:true
                     }
+                    setup.animation.exaggeration=0.1;
                     GLOBALS.changeNodeColor(three.selectedNodeID.x,three.selectedNodeID.y,three.selectedNodeID.z,false,data);
                     GLOBALS.selectNode(three.selectedNodeID.x,three.selectedNodeID.y,three.selectedNodeID.z);
 
@@ -647,7 +648,7 @@ function onDocumentMouseDownThree( event ) {
             var obj=utils.getXYZfromName(intersect.object.name);
             obj=utils.getXYZfromName(three.rollOverMesh.name);
             GLOBALS.addNode (obj.x, obj.y,obj.z);
-            GLOBALS.selectNode (obj.x, obj.y,obj.z);
+            GLOBALS.selectNode (obj.x, obj.y, obj.z);
         }
         if ( intersects1.length > 0 && event.which==3){//right click
                 var intersect = intersects1[ 0 ];
@@ -779,8 +780,11 @@ document.addEventListener('changeNodeColor', function (e) {
 
 document.addEventListener('selectNode', function (e) { 
     var obj=e.detail;
-    three.selectedNodeID=obj
-    [p_x ,p_y ,p_z ,s_x ,s_y,s_z,r_y]=utils.getTransforms(three.grid,obj.x, obj.y,obj.z);
+    
+
+    three.selectedNodeID=obj;
+    
+    [p_x ,p_y ,p_z ,s_x ,s_y,s_z,r_y]=utils.getTransforms(three.grid, obj.x, obj.y,obj.z);
         
     // x = document.getElementById("instructions");x.innerHTML = "Selected node at x="+e.detail.x+", y="+e.detail.y+", z="+e.detail.y+".";
             
@@ -797,7 +801,8 @@ document.addEventListener('selectNode', function (e) {
     three.highlight.position.z=p_z;
 
     var selected = GLOBALS.globalJson.nodes.filter(selected => selected.x == obj.x&&selected.y == obj.y&&selected.z == obj.z)[0];
-    // console.log(selected)
+    
+    
     
 }, false);
 
-- 
GitLab