diff --git a/examples/ex12/ex12.cpp b/examples/ex12/ex12.cpp
index 07655fd..c727064 100644
--- a/examples/ex12/ex12.cpp
+++ b/examples/ex12/ex12.cpp
@@ -47,7 +47,12 @@ int main(int argc, char **argv) {
   // double beta = 0;
   // double gamma = 0.5;
 
+
   ShapeFunctions();
+
+  //printf("execut until here\n");
+  //exit (EXIT_FAILURE);
+
   ReadMaterialProperties();
   /*  Step-1: Calculate the mass matrix similar to that of belytschko. */
   Assembly((char *)"mass"); // Add Direct-lumped as an option
diff --git a/femtech_config.sh b/femtech_config.sh
new file mode 100644
index 0000000..94e4fd1
--- /dev/null
+++ b/femtech_config.sh
@@ -0,0 +1,24 @@
+#!/bin/bash
+#
+#UNAMEX="ubuntu"
+#cd /home/$UNAMEX
+#sudo apt-get update
+#sudo apt-get install -y build-essential
+#sudo apt-get install -y cmake-curses-gui
+#sudo apt-get install -y libblas-dev liblapack-dev
+## sudo apt-get install -y openmpi-bin openmpi-common openssh-client openssh-server libopenmpi1.10
+#sudo apt-get install -y openmpi-bin openmpi-common libopenmpi-dev
+#git clone https://github.com/PSUCompBio/FemTech
+#cd FemTech
+mkdir build
+cd build
+cmake .. -DEXAMPLES=ON -DEXAMPLE12=ON -DEXAMPLE11=ON
+make -j8
+
+#make -j8 ex12 In the examples/ex12/ there is also a makefile. use it
+
+#mpirun -n 1 ex11 1-elt-cube.k > debug_log.txt
+#mpirun -n 1 ex12 1-elt-truss.k > debug_log.txt
+#paraview 1-elt-cube.pvd  here the *.pvd file
+
+#fileellelkfalfhf
diff --git a/include/FemTech.h b/include/FemTech.h
index e4b51e8..81adc7d 100644
--- a/include/FemTech.h
+++ b/include/FemTech.h
@@ -32,22 +32,29 @@ int LineToArray(const bool IntOrFloat, const bool CheckLastVal, \
     const char *Delim = " \t", void **Array = NULL);
 bool ReadInputFile(const char *FileName);
 bool PartitionMesh();
+
 void GaussQuadrature3D(int element, int nGaussPoint, double *Chi,double *GaussWeights);
+
 void ShapeFunctions();
 void ShapeFunction_C3D8(int e, int gp, double *Chi, double *detJ);
 void ShapeFunction_C3D4(int e, int gp, double *Chi, double *detJ);
 void ShapeFunction_T3D2(int e, int gp, double *Chi, double *detJ);
+void get_cos(int e, double *cx, double *cy, double *cz);
+
 void ReadMaterialProperties();
 void ReadBoundaryCondition(void);
+
 void AllocateArrays();
 
 void Assembly(char *operation);
 void StiffnessElementMatrix(double* Ke, int e);
 void MassElementMatrix(double* Me, int e);
+
 void WriteVTU(const char* FileName, int step, double time);
 void WritePVD(const char* FileName, int step, double time);
+
 void FreeArrays();
-void ReadMaterialProperties();
+
 void ApplySteadyBoundaryConditions(void);
 void SolveSteadyImplicit(void);
 void SolveUnsteadyNewmarkImplicit(double beta, double gamma, double dt, \
diff --git a/python_code/Main.py b/python_code/Main.py
new file mode 100644
index 0000000..72f31f6
--- /dev/null
+++ b/python_code/Main.py
@@ -0,0 +1,225 @@
+# ****************************************************************************
+# Name           : Main.py
+# Author         : Andres Valdez
+# Version        : 1.0
+# Description    : 3D truss element, pull-back matrix operations
+#                  This set of functions receives Nodal displaments, and
+#                  provides the transformation matrix, from global to local
+# Data	     : 05-05-2019
+# ****************************************************************************
+
+from __future__ import unicode_literals
+import os, sys
+import numpy as np
+import scipy as scp
+from scipy import optimize
+from scipy.linalg import solve
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+import matplotlib.mlab as mlab
+from math import pi, sin, sqrt, pow
+
+def get_cos(x,y,z,l0):
+    '''
+    Receives the variations of lengths between two consecutive mechanical
+    states. And computes the direction towards the initial state.
+    '''
+    cx = (x[1] - x[0]) / l0
+    cy = (y[1] - y[0]) / l0
+    cz = (z[1] - z[0]) / l0
+    return cx, cy, cz
+
+def map_global_to_local(cx,cy,cz):
+    '''
+    Here we define the 3d to 2d mapping. (x,y,z)--->(xo,xf)
+    '''
+    T = np.zeros((2,6))
+    T[0,0], T[0,1], T[0,2] = cx, cy, cz
+    T[1,3], T[1,4], T[1,5] = cx, cy, cz
+    
+    return T
+
+def vec_global_to_local(T,f):
+    '''
+    Receives any global vector (elemental) and the transformation matrix.
+    Returns the local vector for Isoparametric implementation
+    '''
+    return np.matmul(T,f)
+ 
+def vec_local_to_global(T,f):
+    '''
+    Receives the local vector and the global to local mapping, and returns
+    the global vector.
+    '''
+    return np.matmul(T.T,f)
+   
+def gauss_point(npts):
+    '''
+    This function returns the localization of a gauss points and weights
+    consideering the ref domain (-1,1).
+    '''
+    pg, wg = np.zeros(npts) , np.zeros(npts)
+    
+    if(npts==1):
+        pg[0] = 0.0
+        wg[0] = 2
+    if(npts==2):
+        pg[0] , pg[1] = 1/sqrt(3) , -1/sqrt(3)
+        wg[0] , wg[1] = 1 , 1
+    return pg, wg
+
+def shape_func_local(npts):
+    '''
+    Here we define the local shape functions. Isoparametric (-1,1).
+    Just for a 2-node line
+    shl(0,i,l) = x-derivative of shape function i
+    shl(1,i,l) =  shape function defined at node i
+    l: gauss point
+    '''
+    shl = np.zeros((2,2,npts))
+    
+    pg , wg = gauss_point(npts)
+    
+    for k in range(npts):
+        shl[0,0,k] = -0.5
+        shl[1,0,k] = -0.5 * (pg[k] - 1.0)
+        
+        shl[0,1,k] = 0.5
+        shl[1,1,k] = -0.5 * (pg[k] + 1.0)
+    
+    return shl
+
+def shape_func_global(npts,x):
+    '''
+    Here we define the global shape functions
+    shg(0,i,l) = x-derivative of shape function i
+    shg(1,i,l) =  shape function defined at node i
+    x[0],x[1] local-coordinates of line nodes
+    
+    '''
+    det = np.zeros(npts)
+    shg = np.zeros((2,2,npts))
+    
+    shl = shape_func_local(npts)
+    
+    for k in range(npts):
+        det[k] = sqrt((x[1]-x[0])**2)
+        
+        for i in range(2):
+            shg[0,i,k] = shl[0,i,k] / det[k]
+            shg[1,i,k] = shl[1,i,k]
+        
+    return shg,det
+
+def pos_truss(shg,x,y,z,l):
+    '''
+    Here we define the x,y,z-variables for integration @ each gauss-point
+    '''
+    px = 0.0
+    
+    for k in range(2):
+        px = px + shg[1,k,l] * x[k]
+    return px,py,pz
+
+def elemental_force(force,x,npts):
+    '''
+    Here we define the elemental force vector.
+    Receives local forces, local coordinates
+    '''
+    fe = np.zeros(2)
+    pg , wg = gauss_point(npts)
+    shg , det = shape_func_global(npts,x)
+    for l in range(npts):
+        for k in range(2):
+            fe[k] = fe[k] + force[k] * wg[l] * shg[1,k,l] * det[l]
+        
+    return fe
+
+def Green_Lagrange(npts,x,u):
+    '''
+    Here we define the elemental Green-Lagrange tensor.
+    Receives number of gauss points and local coordinates and a local
+    displacement.
+    '''
+    du = np.zeros(2)
+    Eu = np.zeros((2,2))
+    pg , wg = gauss_point(npts)
+    shg , det = shape_func_global(npts,x)
+    for l in range(npts):
+        for k in range(2):
+            du[k] = du[k] + shg[0,k,l] * wg[l] * det[l] * u[k]
+    
+    Eu[0,0] = du[0] + 0.5 * du[0]**2
+    Eu[1,1] = du[1] + 0.5 * du[1]**2
+    return Eu
+
+def elemental_int_force(Eu,E,A):
+    '''
+    Receives the Green-Lagrange strain tensor, the Young modulus, and the
+    transversal area. Returns the internal-efforts at both extremal nodes.
+    '''
+    fi = np.zeros(2)
+    fi[0] = E * A * (Eu[0,0] + Eu[0,1])
+    fi[1] = E * A * (Eu[1,0] + Eu[1,1])
+    #Note that Eu is always diagonal. Its a truss, the kinematics is penalized.
+    return fi
+    
+if __name__ == "__main__":
+
+    #
+    # Workflow
+    #
+    ###################################################################
+    ## 0) Definition of a mechanical state
+    ###################################################################
+    x , y , z = np.array([0,1]) , np.array([0,1]), np.array([0,1]) # This is the initial state
+    ###################################################################
+    ## 1) From the previous two-solutions, compute the new direction
+    ###################################################################
+    l = np.sqrt((x[0] - x[1])**2 + (y[0] - y[1])**2 + (z[0] - z[1])**2) # This is the element's length
+    cos_theta = get_cos(x,y,z,l)
+    print('the angles',type(cos_theta),cos_theta)
+    ###################################################################
+    ## 2) Evaluate the transformation matrix
+    ###################################################################
+    A = map_global_to_local(cos_theta[0],cos_theta[1],cos_theta[2])
+    print('transformation matrix')
+    print(A)
+    ###################################################################
+    ## 3) Convert the global force in 3d space to a 2d space.
+    ###################################################################
+    fg = 9.8 * np.array([0,0,-1,0,0,-1]) # This is a hand made external loading
+                                         # in our case gravity, negative z direction
+    
+    fl = vec_global_to_local(A,fg)
+    print('global force',fg)
+    print('local force',fl)
+    ###################################################################
+    ## 4) Convert the global coordinates in 3d space to a 2d space.
+    ###################################################################
+    ug = np.array([x[0],x[1],y[0],y[1],z[0],z[1]])
+    ul = vec_global_to_local(A,ug)
+    print('global coords',ug)
+    print('local coords',ul)
+    ###################################################################
+    ## 5) Here starts the isoparametric evaluations. Integrate the external-forces
+    ###################################################################
+    npts = 1 # Number of Gauss points
+    int_fl = elemental_force(fl,ul,npts)
+    
+    print('integrate force',int_fl)
+    ###################################################################
+    ## 5) Integrate the internal-forces
+    ###################################################################
+    E , A = 1 , 1
+    int_def = Green_Lagrange(npts,ul,2*ul)
+    print('The green lagrange strain tensor')
+    print(int_def)
+    
+    int_force = elemental_int_force(int_def,E,A)
+    
+    print('The internal forces')
+    print(int_force)
+    
+    sys.exit()
+    
diff --git a/python_code/truss-3d.pdf b/python_code/truss-3d.pdf
new file mode 100644
index 0000000..d65c0b8
Binary files /dev/null and b/python_code/truss-3d.pdf differ
diff --git a/src/fem/ShapeFunctions/ShapeFunction_T3D2.cpp b/src/fem/ShapeFunctions/ShapeFunction_T3D2.cpp
index 4b7b188..ccf7097 100644
--- a/src/fem/ShapeFunctions/ShapeFunction_T3D2.cpp
+++ b/src/fem/ShapeFunctions/ShapeFunction_T3D2.cpp
@@ -1,88 +1,129 @@
 #include "FemTech.h"
 
+void get_cos(int e, double *cx, double *cy, double *cz){
+      // int e :: elements' id
+      // double cx,cy,cz :: truss direction
+      // int eptr[e+1] :: number of nodes per element
+      // For each element e, you obtain the involved nodes looping until e+1.
+      double l0 = {0};
+      double x[2] = {0};
+      double y[2] = {0};
+      double z[2] = {0};
+      
+      for (int k=0; k < eptr[e+1]; k++){
+            x[k] = coordinates[ndim*connectivity[eptr[e]+k]+0];
+            y[k] = coordinates[ndim*connectivity[eptr[e]+k]+1];
+            z[k] = coordinates[ndim*connectivity[eptr[e]+k]+2];
+      }
+//      Checking the coordinates of each node, for element id:e
+
+       for (int k=0; k < eptr[e+1]; k++){
+            printf("nodes per elem: %d, ide node: %d\n",eptr[e+1],k);
+            printf("coordinates node: %d, are [%5.3e,%5.3e,%5.3e]\n",k,x[k],y[k],y[k]);
+       }
+	 printf("\n");
+      
+      l0 = sqrt(pow(x[1]-x[0],2) + pow(y[1]-y[0],2) + pow(z[1]-z[0],2));
+      
+      printf("length of the truss: %5.3e\n",l0);
+      
+      *cx = ((x[1] - x[0]) / l0);
+      *cy = ((y[1] - y[0]) / l0);
+      *cz = ((z[1] - z[0]) / l0);
+      
+      printf("the truss direction: [%5.3e,%5.3e,%5.3e]\n",*cx,*cy,*cz);
+      
+      return;
+}
+
+//https://www.sanfoundry.com/cpp-programming-examples-numerical-problems-algorithms/
+//https://www.sanfoundry.com/c-programming-examples-numerical-problems-algorithms/
+//https://www.sanfoundry.com/1000-c-algorithms-problems-programming-examples/
 
 void ShapeFunction_T3D2(int e, int gp, double *Chi, double *detJ){
+      // int e :: element ID
+      // int gp :: Gauss point ID
+      // double *Chi :: Gauss points array coordinates
+      // double *detJ :: Jacs determinant array
+
+      double chi, eta, iota;
+      chi =  Chi[ndim*gp + 0]; // local x
+      eta =  Chi[ndim*gp + 1]; // local y
+      iota = Chi[ndim*gp + 2]; // local z
+
+      // The shape functions
+      int g = GaussPoints[e];
+      const int indexShp = gptr[e] + gp * g;
+      const int indexDshp = dsptr[e] + gp * g * ndim;
+
+      shp[indexShp + 0] = 0.5*(1.0-chi);
+      shp[indexShp + 1] = 0.5*(1.0+chi);
+
+      // The first derivatives
+
+      // with respect to local x
+      dshp[indexDshp + ndim * 0 + 0] = -0.5; // derivative of N_1
+      dshp[indexDshp + ndim * 1 + 0] =  0.5; // derivative of N_2
+
+      // with respect to local y
+      dshp[indexDshp + ndim * 0 + 1] =  0.0; // derivative of N_1
+      dshp[indexDshp + ndim * 1 + 1] =  0.0; // derivative of N_2
+
+      // with respect to local z
+      dshp[indexDshp + ndim * 0 + 2] =  0.0; // derivative of N_1
+      dshp[indexDshp + ndim * 1 + 2] =  0.0; // derivative of N_2
+
+      // Here we start the transformation 6-dof to 2-dof.
+
+      double cx,cy,cz;
+      double T[12];
+
+      cx = 0.0;
+      cy = 0.0;
+      cz = 0.0;
+
+      printf("the initial truss direction: [%5.3e,%5.3e,%5.3e]\n",cx,cy,cz);
+
+      get_cos(e, &cx, &cy, &cz);
+
+      printf("the truss direction: [%5.3e,%5.3e,%5.3e]\n",cx,cy,cz);
+
+      printf("Number of spatial dimensions: %d\n",ndim);
+
+      // Compute the Jacobian's determinant
+
+      double xl[2] = {0};
+      double x[2] = {0};
+      double y[2] = {0};
+      double z[2] = {0};
+
+      for (int k=0; k < eptr[e+1]; k++){
+            x[k] = coordinates[ndim*connectivity[eptr[e]+k]+0];
+            y[k] = coordinates[ndim*connectivity[eptr[e]+k]+1];
+            z[k] = coordinates[ndim*connectivity[eptr[e]+k]+2];
+      }
+
+      for (int k=0; k < eptr[e+1]; k++){
+            xl[k] = x[k]*cx + y[k]*cy + z[k]*cz;
+      }
+
+      printf("Elem id: %d, coord-0: %5.3e, coord-1: %5.3e\n",e,xl[0],xl[1]);
+
+      detJ[gp] = xl[1] - xl[0];
+
+      // The Global first derivatives. Only for direction x. 
+
+      // with respect to local x
+      dshp[indexDshp + ndim * 0 + 0] = dshp[indexDshp + ndim * 0 + 0] / detJ[gp]; // derivative of N_1
+      dshp[indexDshp + ndim * 1 + 0] = dshp[indexDshp + ndim * 1 + 0] / detJ[gp]; // derivative of N_2
 
-	double chi, eta, iota;
-	chi =  Chi[ndim*gp + 0];
-	eta =  Chi[ndim*gp + 1];
-	iota = Chi[ndim*gp + 2];
-
-	// The shape functions
-	int g = GaussPoints[e];
-	const int indexShp = gptr[e] + gp * g;
-  const int indexDshp = dsptr[e] + gp * g * ndim;
-
-	shp[indexShp + 0] = 0.5*(1.0-chi);
-	shp[indexShp + 1] = 0.5*(1.0+chi);
-	// shp[indexShp + 2] = iota;
-	// shp[indexShp + 3] = 1.0 - eta - iota - chi;
-
-
-	// The first derivatives
-
-	// with respect to chi
-	dshp[indexDshp + ndim * 0 + 0] = -0.5; // derivative of N_1
- 	dshp[indexDshp + ndim * 1 + 0] =  0.5; // derivative of N_2
-	//dshp[indexDshp + ndim * 2 + 0] =  1.0;
-	//dshp[indexDshp + ndim * 3 + 0] =  1.0;
-
-	// with respect to eta
-	dshp[indexDshp + ndim * 0 + 1] =  0.0; // derivative of N_1
-	dshp[indexDshp + ndim * 1 + 1] =  0.0; // derivative of N_2
-	// dshp[indexDshp + ndim * 2 + 1] =  0.0;
-	// dshp[indexDshp + ndim * 3 + 1] = -1.0;
-	// with respect to iota
-	 dshp[indexDshp + ndim * 0 + 2] =  0.0; // derivative of N_1
-	 dshp[indexDshp + ndim * 1 + 2] =  0.0; // derivative of N_2
-	// dshp[indexDshp + ndim * 2 + 2] =  1.0;
-	// dshp[indexDshp + ndim * 3 + 2] = -1.0;
-
-	//  Compute jacobian transformation
- 	//  Equ 9. 16 3D stress Analysis
-
-	//  x14 y14 z14
-	//  x24 y24 z24
-	//  x34 y34 z34
-	double xs[ndim*ndim];
-
-	// // first row
-	// xs[0]=coordinates[ndim*connectivity[index+node1]+x] - coordinates[ndim*connectivity[index+node4]+x];
-	// xs[3]=coordinates[ndim*connectivity[index+node1]+y] - coordinates[ndim*connectivity[index+node4]+y];
-	// xs[6]=coordinates[ndim*connectivity[index+node1]+z] - coordinates[ndim*connectivity[index+node4]+z];
-	// // second row
-	// xs[1]=coordinates[ndim*connectivity[index+node2]+x] - coordinates[ndim*connectivity[index+node4]+x];
-	// xs[4]=coordinates[ndim*connectivity[index+node2]+y] - coordinates[ndim*connectivity[index+node4]+y];
-	// xs[7]=coordinates[ndim*connectivity[index+node2]+z] - coordinates[ndim*connectivity[index+node4]+z];
-	// // third row
-	// xs[2]=coordinates[ndim*connectivity[index+node3]+x] - coordinates[ndim*connectivity[index+node4]+x];
-	// xs[5]=coordinates[ndim*connectivity[index+node3]+y] - coordinates[ndim*connectivity[index+node4]+y];
-	// xs[8]=coordinates[ndim*connectivity[index+node3]+z] - coordinates[ndim*connectivity[index+node4]+z];
-  // double det, J_Inv[9];
-	//
-  // inverse3x3Matrix(xs, J_Inv, &det);
-  // detJ[gp] = det;
-	//
-  // // Transform derivatives to global co-ordinates
-  // double c1, c2, c3;
-  // int baseIndex;
-  // for (int i = 0; i < nShapeFunctions[e]; ++i) {
-  //   baseIndex = dsptr[e] + gp * g*ndim + ndim * i;
-  //   c1 = dshp[baseIndex]*J_Inv[0]+dshp[baseIndex+1]*J_Inv[3]+dshp[baseIndex+2]*J_Inv[6];
-  //   c2 = dshp[baseIndex]*J_Inv[1]+dshp[baseIndex+1]*J_Inv[4]+dshp[baseIndex+2]*J_Inv[7];
-  //   c3 = dshp[baseIndex]*J_Inv[2]+dshp[baseIndex+1]*J_Inv[5]+dshp[baseIndex+2]*J_Inv[8];
-	//
-  //   dshp[baseIndex] = c1;
-  //   dshp[baseIndex+1] = c2;
-  //   dshp[baseIndex+2] = c3;
-  // }
-	// //for debugging can be removed...
-	// if (debug && 1==0) {
-  //   printf("DEBUG J Inv\n");
-	// 	printf("%8.4e %8.4e %8.4e\n", J_Inv[0], J_Inv[3], J_Inv[6]);
-	// 	printf("%8.4e %8.4e %8.4e\n", J_Inv[1], J_Inv[4], J_Inv[7]);
-	// 	printf("%8.4e %8.4e %8.4e\n", J_Inv[2], J_Inv[5], J_Inv[8]);
-	// 	printf("\n");
-	// }
-   return;
+      // //for debugging can be removed...
+      // if (debug && 1==0) {
+      //   printf("DEBUG J Inv\n");
+      // 	printf("%8.4e %8.4e %8.4e\n", J_Inv[0], J_Inv[3], J_Inv[6]);
+      // 	printf("%8.4e %8.4e %8.4e\n", J_Inv[1], J_Inv[4], J_Inv[7]);
+      // 	printf("%8.4e %8.4e %8.4e\n", J_Inv[2], J_Inv[5], J_Inv[8]);
+      // 	printf("\n");
+      // }
+      return;
 }
diff --git a/src/fem/SolidMechanics/GetForce.cpp b/src/fem/SolidMechanics/GetForce.cpp
index d3a139b..8b1ce0d 100644
--- a/src/fem/SolidMechanics/GetForce.cpp
+++ b/src/fem/SolidMechanics/GetForce.cpp
@@ -12,8 +12,9 @@
 void GetForce() {
 
 	if(ndim == 1){
-		printf("GetForce function not yet implemented for 1D. Bye.\n");
-		exit(0);
+		GetForce_3D();
+		printf("Here we are outside of the gates of valhala\n");
+		exit (EXIT_FAILURE);
 	}
 	else if (ndim == 2){
 		printf("GetForce function not yet implemented for 2D. Bye.\n");
diff --git a/src/fem/SolidMechanics/GetForce_3D.cpp b/src/fem/SolidMechanics/GetForce_3D.cpp
index ab1c813..3c25724 100644
--- a/src/fem/SolidMechanics/GetForce_3D.cpp
+++ b/src/fem/SolidMechanics/GetForce_3D.cpp
@@ -21,20 +21,23 @@ void GetForce_3D() {
 		// force calculaton for hexes, tets and quads
 		for(int j=0; j<GaussPoints[i]; j++) {
 			// truss elements are unique b/c the 3D formulation is not like solid
-			// elmeents like quads or hexes, or tets.
-	    if(strcmp(ElementType[i], "T3D2") == 0){
+			// elements like quads or hexes, or tets.
+	            if(strcmp(ElementType[i], "T3D2") == 0){
+				//printf("Calculate the forces for truss\n");
+				//exit (EXIT_FAILURE);
 				TrussStressForceUpdate(i,j,fintLocal);
 				//CalculateDeformationGradient(i, j);
 				//exit(0);
-			}else{
-	      // calculate F^n
+		      }
+		      else{
+	            // calculate F^n
 				CalculateDeformationGradient(i, j);
-	      // Calculate Determinant of F
+	            // Calculate Determinant of F
 				DeterminateF(i, j);
-	      // Calculate sigma^n
+	            // Calculate sigma^n
 				StressUpdate(i, j);
-			  InternalForceUpdate(i, j, fintLocal);
-			} // else
+			      InternalForceUpdate(i, j, fintLocal);
+		      } // else
 		} //loop on gauss points
     // Move Local internal for to global force
     for (int k = 0; k < nNodes; ++k) {
diff --git a/src/fem/SolidMechanics/TrussStressForceUpdate.cpp b/src/fem/SolidMechanics/TrussStressForceUpdate.cpp
index b8412f8..cfb7668 100644
--- a/src/fem/SolidMechanics/TrussStressForceUpdate.cpp
+++ b/src/fem/SolidMechanics/TrussStressForceUpdate.cpp
@@ -2,316 +2,83 @@
 #include "blas.h"
 
 void TrussStressForceUpdate(int e, int gp, double *force){
-
-	int x = 0, y = 1, z = 2;
-	int node1 = 0, node2 = 1;
-	double x01,x02,y01,y02,z01,z02;
-	double x1,x2,y1,y2,z1,z2;
-	int index=eptr[e];
-	x01 = coordinates[ndim*connectivity[index+node1]+x];
-	x02 = coordinates[ndim*connectivity[index+node2]+x];
-	y01 = coordinates[ndim*connectivity[index+node1]+y];
-	y02 = coordinates[ndim*connectivity[index+node2]+y];
-	z01 = coordinates[ndim*connectivity[index+node1]+z];
-	z02 = coordinates[ndim*connectivity[index+node2]+z];
-
-	x1 = coordinates[ndim*connectivity[index+node1]+x]+displacements[ndim*connectivity[index+node1]+x];
-	x2 = coordinates[ndim*connectivity[index+node2]+x]+displacements[ndim*connectivity[index+node2]+x];
-	y1 = coordinates[ndim*connectivity[index+node1]+y]+displacements[ndim*connectivity[index+node1]+y];
-	y2 = coordinates[ndim*connectivity[index+node2]+y]+displacements[ndim*connectivity[index+node2]+y];
-	z1 = coordinates[ndim*connectivity[index+node1]+z]+displacements[ndim*connectivity[index+node1]+z];
-	z2 = coordinates[ndim*connectivity[index+node2]+z]+displacements[ndim*connectivity[index+node2]+z];
-  double le0= sqrt( pow((x02-x01),2.0) + pow((y02-y01),2.0) + pow((z02-z01),2.0) );
-	double le = sqrt( pow((x2-x1),  2.0) +  pow((y2-y1), 2.0) + pow((z2-z1),  2.0) );
-	// axial stretch
-	double lambda = le/le0;
-	printf("coords = %3.3e\n", coordinates[ndim*connectivity[index+node1]+x]);
-	printf("disp = %3.3e\n", displacements[ndim*connectivity[index+node1]+x]);
-	printf("x2 = %3.3f\n",x2 );
-printf("x1 = %3.3f\n",x1 );
-	printf("(x2-x1)^2 = %3.3f\n",pow((x2-x1),  2.0) );
-	printf("(y2-y1)^2 = %3.3f\n",pow((y2-y1), 2.0));
-printf("(z2-z1)^2 = %3.3f\n",pow((z2-z1), 2.0));
-printf("le0 = %3.3f\n",le0);
-printf("le = %3.3f\n",le);
-	printf("lambda = %3.3f\n",lambda);
-	// direction cosines
-	double l = (x2-x1)/le;
-	double m = (y2-y1)/le;
-	double n = (z2-z1)/le;
-
-	printf("l = %3.3f\n",l);
-	printf("m = %3.3f\n",m);
-		printf("n = %3.3f\n",n);
-	//create transformation matrix. following logan, section 3.7, page 93, equation 3.7.6
-	double T[6];
-	memset(T, 0.0, sizeof(T));
-	T[0]=-1.0*l;
-	T[1]=-1.0*m;
-	T[2]=-1.0*n;
-	T[3]=l;
-	T[4]=m;
-	T[5]=n;
-	//T[0][1]=0;T[1][1]=0;T[2][1]=0;T[3][1]=l;T[4][1]=m;T[5][1]=n;
-
-	//create local displacemtn vector
-	double d[6];
-	memset(d, 0.0, sizeof(d));
-	d[0]=displacements[ndim*connectivity[index+node1]+x];
-	d[1]=displacements[ndim*connectivity[index+node1]+y];
-	d[2]=displacements[ndim*connectivity[index+node1]+z];
-	d[3]=displacements[ndim*connectivity[index+node2]+x];
-	d[4]=displacements[ndim*connectivity[index+node2]+y];
-	d[5]=displacements[ndim*connectivity[index+node2]+z];
-
-	// Muitiple T * d
-	double sum = 0.0;
-	for (int i = 0; i < 6; i++) {
-		sum = sum + T[i] * d[i];
-	}
-
-
-  //green Lagrange strain
-  double gls=0.5*(pow(lambda,2.0)-1.0);
-	double area=1.0;
-	double youngs = 1.0;
-  double cauchystress = youngs * gls;
-	double truss_axial_force  = cauchystress*area;
-	printf("gls = %3.3f\n",gls);
-
-	printf("s = %3.3f, f = %3.3f\n",cauchystress, truss_axial_force);
-
-	// 6 values saved per gauss point for 3d
-	// in voigt notation, sigma11
-	cauchy[cptr[e] + 6 * gp + 0] = cauchystress;
-	// in voigt notation, sigma22
-	cauchy[cptr[e] + 6 * gp + 1] = 0;
-	// in voigt notation, sigma33
-	cauchy[cptr[e] + 6 * gp + 2] = 0;
-	// in voigt notation, sigma23
-	cauchy[cptr[e] + 6 * gp + 3] = 0;
-	// in voigt notation, sigma13
-	cauchy[cptr[e] + 6 * gp + 4] = 0;
-	// in voigt notation, sigma12
-	cauchy[cptr[e] + 6 * gp + 5] = 0;
-
-	force[0] = truss_axial_force;
-	force[1] = 0.0;
-	force[2] = 0.0;
-	force[3] = 0.0;
-	force[4] = 0.0;
-	force[5] = 0.0;
-
-
-	//
-  // // Calculate Tranformation
-	// double dir_truss_x[3];
-	// double dir_truss_y[3];
-	// double dir_truss_z[3];
-	// dir_truss_x[0] = x2 - x1;
-	// dir_truss_x[1] = y2 - y1;
-	// dir_truss_x[2] = z2 - z1;
-	// if (dir_truss_x[0] != 0) {
-	// 	dir_truss_y[1] = 1;
-	// 	dir_truss_y[2] = 0;
-	// 	dir_truss_y[0] = -(((dir_truss_y[1] * dir_truss_x[1]) + (dir_truss_y[2] * dir_truss_x[2])) / dir_truss_x[0]);
-	// }
-	// else {
-	// 	if ((dir_truss_x[0] == 0) && (dir_truss_x[1] != 0)) {
-	// 		dir_truss_y[0] = 1;
-	// 		dir_truss_y[2] = 1;
-	// 		dir_truss_y[1] = -(((dir_truss_y[2] * dir_truss_x[2]) + (dir_truss_y[0] * dir_truss_x[0])) / dir_truss_x[1]);
-	// 	}
-	// 	else if ((dir_truss_x[0] == 0) && (dir_truss_x[2] != 0)) {
-	// 		dir_truss_y[0] = 0;
-	// 		dir_truss_y[1] = 1;
-	// 		dir_truss_y[2] = -(((dir_truss_y[0] * dir_truss_x[0]) + (dir_truss_y[1] * dir_truss_x[1])) / dir_truss_x[2]);
-	// 	}
-	// 	else {
-	// 		printf("Truss element is not practically possible.\n");
-	// 	}
-	// }
-  // // cross product: of x cross y to get z
-	// double a1 = dir_truss_x[0];
-	// double a2 = dir_truss_x[1];
-	// double a3 = dir_truss_x[2];
-	// double b1 = dir_truss_y[0];
-	// double b2 = dir_truss_y[1];
-	// double b3 = dir_truss_y[2];
-  // dir_truss_z[0] = (a2*b3 - a3*b2);
-	// dir_truss_z[1] = (a1*b3 - a3*b1);
-	// dir_truss_z[2] = (a1*b1 - a2*b1);
-	//
-	// //normalize each vector
-	// double norm_x = sqrt(dir_truss_x[0]*dir_truss_x[0] + dir_truss_x[1]*dir_truss_x[1] + dir_truss_x[2]*dir_truss_x[2] );
-	// double norm_y = sqrt(dir_truss_y[0]*dir_truss_y[0] + dir_truss_y[1]*dir_truss_y[1] + dir_truss_y[2]*dir_truss_y[2] );
-	// double norm_z = sqrt(dir_truss_z[0]*dir_truss_z[0] + dir_truss_z[1]*dir_truss_z[1] + dir_truss_z[2]*dir_truss_z[2] );
-	// // x
-  // dir_truss_x[0] =  dir_truss_x[0]/norm_x;
-	// dir_truss_x[1] =  dir_truss_x[1]/norm_x;
-	// dir_truss_x[2] =  dir_truss_x[2]/norm_x;
- 	// // y
-	// dir_truss_y[0] =  dir_truss_y[0]/norm_y;
-	// dir_truss_y[1] =  dir_truss_y[1]/norm_y;
-	// dir_truss_y[2] =  dir_truss_y[2]/norm_y;
-	// // z
-	// dir_truss_z[0] =  dir_truss_z[0]/norm_z;
-	// dir_truss_z[1] =  dir_truss_z[1]/norm_z;
-	// dir_truss_z[2] =  dir_truss_z[2]/norm_z;
-	//
-	// // define global coordinate system
-	// double dir_global_x[3];
-	// dir_global_x[0]=1;
-	// dir_global_x[1]=0;
-	// dir_global_x[2]=0;
-	// double dir_global_y[3];
-	// dir_global_y[0]=0;
-	// dir_global_y[1]=1;
-	// dir_global_y[2]=0;
-	// double dir_global_z[3];
-	// dir_global_z[0]=0;
-	// dir_global_z[1]=0;
-	// dir_global_z[2]=1;
-	// //normalize each global vector - not really needed but left in b/c it can
-	// // accomadate different global orientations in future  - IS THIS NEEDED??
-	// double gnorm_x = sqrt(dir_global_x[0]*dir_global_x[0] + dir_global_x[1]*dir_global_x[1] + dir_global_x[2]*dir_global_x[2] );
-	// double gnorm_y = sqrt(dir_global_y[0]*dir_global_y[0] + dir_global_y[1]*dir_global_y[1] + dir_global_y[2]*dir_global_y[2] );
-	// double gnorm_z = sqrt(dir_global_z[0]*dir_global_z[0] + dir_global_z[1]*dir_global_z[1] + dir_global_z[2]*dir_global_z[2] );
-	//
-	// // define direction cosines
-	// double dot_global_x_truss_x = dir_global_x[0]*dir_truss_x[0] + dir_global_x[1]*dir_truss_x[1] + dir_global_x[2]*dir_truss_x[2];
-	// double dot_global_x_truss_y = dir_global_x[0]*dir_truss_y[0] + dir_global_x[1]*dir_truss_y[1] + dir_global_x[2]*dir_truss_y[2];
-	// double dot_global_x_truss_z = dir_global_x[0]*dir_truss_z[0] + dir_global_x[1]*dir_truss_z[1] + dir_global_x[2]*dir_truss_z[2];
-	// double dot_global_y_truss_x = dir_global_y[0]*dir_truss_x[0] + dir_global_y[1]*dir_truss_x[1] + dir_global_y[2]*dir_truss_x[2];
-	// double dot_global_y_truss_y = dir_global_y[0]*dir_truss_y[0] + dir_global_y[1]*dir_truss_y[1] + dir_global_y[2]*dir_truss_y[2];
-	// double dot_global_y_truss_z = dir_global_y[0]*dir_truss_z[0] + dir_global_y[1]*dir_truss_z[1] + dir_global_y[2]*dir_truss_z[2];
-	// double dot_global_z_truss_x = dir_global_z[0]*dir_truss_x[0] + dir_global_z[1]*dir_truss_x[1] + dir_global_z[2]*dir_truss_x[2];
-	// double dot_global_z_truss_y = dir_global_z[0]*dir_truss_y[0] + dir_global_z[1]*dir_truss_y[1] + dir_global_z[2]*dir_truss_y[2];
-	// double dot_global_z_truss_z = dir_global_z[0]*dir_truss_z[0] + dir_global_z[1]*dir_truss_z[1] + dir_global_z[2]*dir_truss_z[2];
-	// // l
-	// double l_1 = dot_global_x_truss_x / (gnorm_x * norm_x);
-	// double l_2 = dot_global_x_truss_y / (gnorm_x * norm_y);
-	// double l_3 = dot_global_x_truss_z / (gnorm_x * norm_z);
-  // // m
-	// double m_1 = dot_global_y_truss_x / (gnorm_y * norm_x);
-	// double m_2 = dot_global_y_truss_y / (gnorm_y * norm_y);
-	// double m_3 = dot_global_y_truss_z / (gnorm_y * norm_z);
-	// // n
-	// double n_1 = dot_global_z_truss_x / (gnorm_z * norm_x);
-	// double n_2 = dot_global_z_truss_y / (gnorm_z * norm_y);
-	// double n_3 = dot_global_z_truss_z / (gnorm_z * norm_z);
-	//
-	// // choice - 1: Stress Transformation Matrix
-	// // choice - 2: Strain Transformation Matrix
-	// // choice - 3: 3X3 transformation matrix as produced in my thesis in electrophysiology chapter
-	// // int choice = 3;
-	// // if (choice == 2) {
-	// // 	double T[6][6];
-	// // 	memset(T, 0.0, sizeof(T));
-	// // 	T[0][0]=l_1*l_1;T[1][0]=m_1*m_1;T[2][0]=n_1*n_1;T[3][0]=l_1*m_1;T[4][0]=m_1*n_1;T[5][0]=n_1*l_1;
-	// // 	T[0][1]=l_2*l_2;T[1][1]=m_2*m_2;T[2][1]=n_2*n_2;T[3][1]=l_2*m_2;T[4][1]=m_2*n_2;T[5][1]=n_2*l_2;
-	// // 	T[0][2]=l_3*l_3;T[1][2]=m_3*m_3;T[2][2]=n_3*n_3;T[3][2]=l_3*m_3;T[4][2]=m_3*n_3;T[5][2]=n_3*l_3;
-	// // 	T[0][3]=2.0*l_1*l_2;T[1][3]=2.0*m_1*m_2;T[2][3]=2.0*n_1*n_2;T[3][3]=(l_1*m_2)+(l_2*m_1);T[4][3]=(m_1*n_2)+(m_2*n_1);T[5][3]=(l_1*n_2)+(l_2*n_1);
-	// // 	T[0][4]=2.0*l_2*l_3;T[1][4]=2.0*m_2*m_3;T[2][4]=2.0*n_2*n_3;T[3][4]=(l_2*m_3)+(l_3*m_2);T[4][4]=(m_2*n_3)+(m_3*n_2);T[5][4]=(l_2*n_3)+(l_3*n_2);
-	// // 	T[0][5]=2.0*l_1*l_3;T[1][5]=2.0*m_1*m_3;T[2][5]=2.0*n_1*n_3;T[3][5]=(l_1*m_3)+(l_3*m_1);T[4][5]=(m_1*n_3)+(m_3*n_1);T[5][5]=(l_1*n_3)+(l_3*n_1);
-	// 	// T.row(0) << (l_1 * l_1), (m_1 * m_1), (n_1 * n_1), (l_1 * m_1), (m_1 * n_1), (n_1 * l_1);
-	// 	// T.row(1) << (l_2 * l_2), (m_2 * m_2), (n_2 * n_2), (l_2 * m_2), (m_2 * n_2), (n_2 * l_2);
-	// 	// T.row(2) << (l_3 * l_3), (m_3 * m_3), (n_3 * n_3), (l_3 * m_3), (m_3 * n_3), (n_3 * l_3);
-	// 	// T.row(3) << (2 * l_1 * l_2), (2 * m_1 * m_2), (2 * n_1 * n_2), ((l_1 * m_2) + (l_2 * m_1)), ((m_1 * n_2) + (m_2 * n_1)), ((l_1 * n_2) + (l_2 * n_1));
-	// 	// T.row(4) << (2 * l_2 * l_3), (2 * m_2 * m_3), (2 * n_2 * n_3), ((l_2 * m_3) + (l_3 * m_2)), ((m_2 * n_3) + (m_3 * n_2)), ((l_2 * n_3) + (l_3 * n_2));
-	// 	// T.row(5) << (2 * l_1 * l_3), (2 * m_1 * m_3), (2 * n_1 * n_3), ((l_1 * m_3) + (l_3 * m_1)), ((m_1 * n_3) + (m_3 * n_1)), ((l_1 * n_3) + (l_3 * n_1));
-	// // }
-	// // else if (choice == 1) {
-	// // 	double T[6][6];
-	// // 	memset(T, 0.0, sizeof(T));
-	// // 	T[0][0]=l_1*l_1;T[1][0]=m_1*m_1;T[2][0]=n_1*n_1; T[3][0]=2.0*l_1*m_1;T[4][0]=2.0*m_1*n_1;T[5][0]=2.0*n_1*l_1;
-	// // 	T[0][1]=l_2*l_2;T[1][1]=m_2*m_2;T[2][1]=n_2*n_2; T[3][1]=2.0*l_2*m_2;T[4][1]=2.0*m_2*n_2;T[5][1]=2.0*n_2*l_2;
-	// // 	T[0][2]=l_3*l_3;T[1][2]=m_3*m_3;T[2][2]=n_3*n_3; T[3][2]=2.0*l_3*m_3;T[4][2]=2.0*m_3*n_3;T[5][2]=2.0*n_3*l_3;
-	// // 	T[0][3]=l_1*l_2;T[1][3]=m_1*m_2;T[2][3]=n_1*n_2;T[3][3]=(l_1*m_2)+(l_2*m_1);T[4][3]=(m_1*n_2)+(m_2*n_1);T[5][3]=(l_1*n_2)+(l_2*n_1);
-	// // 	T[0][4]=l_2*l_3;T[1][4]=m_2*m_3;T[2][4]=n_2*n_3;T[3][4]=(l_2*m_3)+(l_3*m_2);T[4][4]=(m_2*n_3)+(m_3*n_2);T[5][4]=(l_2*n_3)+(l_3*n_2);
-	// // 	T[0][5]=l_1*l_3;T[1][5]=m_1*m_3;T[2][5]=n_1*n_3;T[3][5]=(l_1*m_3)+(l_3*m_1);T[4][5]=(m_1*n_3)+(m_3*n_1);T[5][5]=(l_1*n_3)+(l_3*n_1);
-	// // 	// T.row(0) << pow(l_1, 2), pow(m_1, 2), pow(n_1, 2), (2 * l_1 * m_1), (2 * m_1 * n_1), (2 * n_1 * l_1);
-	// // 	// T.row(1) << pow(l_2, 2), pow(m_2, 2), pow(n_2, 2), (2 * l_2 * m_2), (2 * m_2 * n_2), (2 * n_2 * l_2);
-	// // 	// T.row(2) << pow(l_3, 2), pow(m_3, 2), pow(n_3, 2), (2 * l_3 * m_3), (2 * m_3 * n_3), (2 * n_3 * l_3);
-	// // 	// T.row(3) << (1 * l_1 * l_2), (1 * m_1 * m_2), (1 * n_1 * n_2), ((l_1 * m_2) + (l_2 * m_1)), ((m_1 * n_2) + (m_2 * n_1)), ((l_1 * n_2) + (l_2 * n_1));
-	// // 	// T.row(4) << (1 * l_2 * l_3), (1 * m_2 * m_3), (1 * n_2 * n_3), ((l_2 * m_3) + (l_3 * m_2)), ((m_2 * n_3) + (m_3 * n_2)), ((l_2 * n_3) + (l_3 * n_2));
-	// // 	// T.row(5) << (1 * l_1 * l_3), (1 * m_1 * m_3), (1 * n_1 * n_3), ((l_1 * m_3) + (l_3 * m_1)), ((m_1 * n_3) + (m_3 * n_1)), ((l_1 * n_3) + (l_3 * n_1));
-	// // }
-	// // else if (choice == 3) {
-	// 	double T[3][3];
-	// 	memset(T, 0.0, sizeof(T));
-	// 	T[0][0]=l_1;T[1][0]=l_2;T[2][0]=l_3;
-	// 	T[0][1]=m_1;T[1][1]=m_2;T[2][1]=m_3;
-	// 	T[0][2]=n_1;T[1][2]=n_2;T[2][2]=n_3;
-	// 	// T.row(0) << l_1, l_2, l_3;
-	// 	// T.row(1) << m_1, m_2, m_3;
-	// 	// T.row(2) << n_1, n_2, n_3;
-	// //}
-	// // else {
-	// // 	printf("WRONG CHOICE OF TRANSFORMATION MATRIX\n");
-	// // }
-	//
-	// // Compute inverse of T
-	// double T_inv[3][3];
-	// double T_inv_transpose[3][3];
-	// memset(T, 0.0, sizeof(T_inv));
-	// // first compute determinant
-	// double tmp1 = T[0][0]*(T[1][1]*T[2][2] - T[1][2]*T[2][1]);
-	// double tmp2 = T[0][1]*(T[1][0]*T[2][2] - T[1][2]*T[2][0]);
-	// double tmp3 = T[0][2]*(T[1][0]*T[2][1] - T[1][1]*T[2][0]);
-	// double detT = tmp1 - tmp2 + tmp3;
-	// // now compute inverse
-	// T_inv[0][0] = 1/detT*(T[1][1]*T[2][2] - T[1][2]*T[2][1]);
-	// T_inv[0][1] = 1/detT*(T[0][2]*T[2][1] - T[0][1]*T[2][2]);
-	// T_inv[0][2] = 1/detT*(T[0][1]*T[1][2] - T[0][2]*T[1][1]);
-	//
-	// T_inv[1][0] = 1/detT*(T[1][2]*T[2][0] - T[1][0]*T[2][2]);
-	// T_inv[1][1] = 1/detT*(T[0][0]*T[2][2] - T[0][2]*T[2][0]);
-	// T_inv[1][2] = 1/detT*(T[0][2]*T[1][0] - T[0][0]*T[1][2]);
-	//
-	// T_inv[2][0] = 1/detT*(T[1][0]*T[2][1] - T[1][1]*T[2][0]);
-	// T_inv[2][1] = 1/detT*(T[0][1]*T[2][0] - T[0][0]*T[2][1]);
-	// T_inv[2][2] = 1/detT*(T[0][0]*T[1][1] - T[0][1]*T[1][0]);
-	// // Finding the transpose of matrix T_inv
-  // for(int i=0; i<ndim; ++i){
-  //   for(int j=0; j<ndim; ++j){
-  //    	T_inv_transpose[j][i] = T_inv[i][j];
-  //   }
-	// }
-	//
-	// // Compute deformation gradient of truss  - assumes incompressibilty
-	// // because F is typically stored at gauss points, the loop over gauss points
-	// // is included. The size of F is defined in Shapefunctions.cpp. It is already
-  // // created but it needs to referenced correctly.
-	// index = fptr[e] + ndim*ndim*gp;
-	// F[index+ndim*0+0] = lambda;
-	// F[index+ndim*1+1] = 1.0/sqrt(lambda);
-	// F[index+ndim*2+2] = 1.0/sqrt(lambda);
-	//
-	// // Muitiple T_inv_transpose * F
-	// double sum = 0;
-	// double multiply[3][3];
-	// for (int c = 0; c < ndim; c++) {
-	// 	for (int d = 0; d < ndim; d++) {
-	// 		for (int k = 0; k < ndim; k++) {
-	// 			sum = sum + T_inv_transpose[c][k]*F[index+ndim*d+k];
-	// 		}
-	// 		multiply[c][d] = sum;
-	// 		sum = 0;
-	// 	}
-	// }
-
-	// // mulitipl F * T_inv
-	// sum = 0;
-	// for (int c = 0; c < ndim; c++) {
-	// 	for (int d = 0; d < ndim; d++) {
-	// 		for (int k = 0; k < ndim; k++) {
-	// 			sum = sum + multiply[c][k]*T_inv[d][k];
-	// 		}
-	// 		F[c][d] = sum;
-	// 		sum = 0;
-	// 	}
-	// }
+      // int e :: element ID
+      // int gp :: Gauss point ID
+      // double *force :: the truss force (internal efforts)
+
+      //
+      // initial coordinates
+      //
+      double x[2] = {0};
+      double y[2] = {0};
+      double z[2] = {0};
+      //
+      // displacements
+      //
+      double dx[2] = {0};
+      double dy[2] = {0};
+      double dz[2] = {0};
+
+      for (int k=0; k < eptr[e+1]; k++){
+            x[k] = coordinates[ndim*connectivity[eptr[e]+k]+0];
+            y[k] = coordinates[ndim*connectivity[eptr[e]+k]+1];
+            z[k] = coordinates[ndim*connectivity[eptr[e]+k]+2];
+
+            dx[k] = x[k] + displacements[ndim*connectivity[eptr[e]+k]+0];
+            dy[k] = y[k] + displacements[ndim*connectivity[eptr[e]+k]+1];
+            dz[k] = z[k] + displacements[ndim*connectivity[eptr[e]+k]+2];
+      }
+
+      //
+      // Turn the 6-dof to 2-dof truss
+      //
+      double cx = 0.0;
+      double cy = 0.0;
+      double cz = 0.0;
+
+      get_cos(e, &cx, &cy, &cz);
+
+      double dl[2] = {0};
+
+      for (int k=0; k < eptr[e+1]; k++){
+            dl[k] = dx[k]*cx + dy[k]*cy + dz[k]*cz;
+      }
+
+      //
+      // Green-Lagrange strain tensor
+      //
+      int wIndex = gpPtr[e] + gp;
+      const int indexDshp = dsptr[e] + gp * GaussPoints[e] * ndim;
+      const double preFactor = gaussWeights[wIndex]*detJacobian[wIndex];
+
+      double Du[2] = {0};
+      double Eu[2] = {0}; /*Since its always diagonal store just two of them*/
+
+      for (int k=0; k < eptr[e+1]; k++){
+            Du[k] = Du[k] + dshp[indexDshp + ndim * k + 0] * preFactor * dl[k]; /*We only use the x-derivatives*/
+            Eu[k] = Du[k] + 0.5 * pow(Du[k],2);
+      }
+
+      //
+      // Define here the local forces (2-dofs)
+      //
+      double fi[2] = {0};
+      double young = 1.0;
+      double area = 1.0;
+      for (int k=0; k < eptr[e+1]; k++){
+            fi[k] = fi[k] + young * area * Eu[k];
+      }
+
+      //
+      // Reconstruct here the local forces (6-dofs) Need to review this
+      //
+      force[0] = cx * fi[0];
+      force[1] = cy * fi[0];
+      force[2] = cz * fi[0];
+
+      force[3] = cx * fi[1];
+      force[4] = cy * fi[1];
+      force[5] = cz * fi[1];
 
 	if(debug && 1==0){
 		printf("--------F for gauss point %d --------\n",gp);
@@ -324,8 +91,5 @@ printf("le = %3.3f\n",le);
 		} //loop on i
 	} // if debug
 
-
-
-
 	return ;
 }