Index: /issm/trunk-jpl/src/jl/solve/analyses/stressbalanceanalysis.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/analyses/stressbalanceanalysis.jl	(revision 26658)
+++ /issm/trunk-jpl/src/jl/solve/analyses/stressbalanceanalysis.jl	(revision 26659)
@@ -68,6 +68,32 @@
 function CreateKMatrix(analysis::StressbalanceAnalysis,element::Tria)# {{{
 
-	#Initialize Element matrix
+	#Internmediaries
+	numnodes = 3
+	
+	#Initialize Element matrix and basis function derivatives
 	Ke = ElementMatrix(element.nodes)
+	dbasis = Matrix{Float64}(undef,numnodes,2)
+
+	#Retrieve all inputs and parameters
+	xyz_list = GetVerticesCoordinates(element.vertices)
+	vx_input = GetInput(element, VxEnum)
+	vy_input = GetInput(element, VyEnum)
+	H_input  = GetInput(element, ThicknessEnum)
+
+	#Start integrating
+	gauss = GaussTria(2)
+	for i in 1:gauss.numgauss
+
+		Jdet = JacobianDeterminant(xyz_list, gauss)
+
+		println(Jdet)
+
+		error("STOP")
+		vx = GetInputValue(vx_input, gauss, i)
+		vy = GetInputValue(vy_input, gauss, i)
+	end
+
+
+	error("STOP")
 
 	return Ke
Index: /issm/trunk-jpl/src/jl/solve/elements.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/elements.jl	(revision 26658)
+++ /issm/trunk-jpl/src/jl/solve/elements.jl	(revision 26659)
@@ -102,2 +102,35 @@
 
 end # }}}
+
+#Finite Element stuff
+function JacobianDeterminant(xyz_list::Matrix{Float64}, gauss::GaussTria) #{{{
+
+	#Get Jacobian Matrix
+	J = Jacobian(xyz_list)
+
+	#Get its determinant
+	Jdet = Matrix2x2Determinant(J)
+
+	#check and return
+	if(Jdet<0) error("negative Jacobian Determinant") end
+	return Jdet
+
+end#}}}
+function Jacobian(xyz_list::Matrix{Float64}) #{{{
+
+	J = Matrix{Float64}(undef,2,2)
+
+	x1 = xyz_list[1,1]
+	y1 = xyz_list[1,2]
+	x2 = xyz_list[2,1]
+	y2 = xyz_list[2,2]
+	x3 = xyz_list[3,1]
+	y3 = xyz_list[3,2]
+
+	J[1,1] = .5*(x2-x1)
+	J[1,2] = sqrt(3)/6*(2*x3 -x1 -x2)
+	J[2,1] = .5*(y2-y1)
+	J[2,2] = sqrt(3)/6*(2*y3 -y1 -y2)
+
+	return J
+end#}}}
Index: /issm/trunk-jpl/src/jl/solve/gauss.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/gauss.jl	(revision 26658)
+++ /issm/trunk-jpl/src/jl/solve/gauss.jl	(revision 26659)
@@ -23,4 +23,5 @@
 		coords3 = [0.333333333333333]
 	elseif(order==2)
+		npoints = 3
 		weights = [0.577350269189625; 0.577350269189625; 0.577350269189625]
 		coords1 = [0.666666666666667; 0.166666666666667; 0.166666666666667]
@@ -28,5 +29,5 @@
 		coords3 = [0.166666666666667; 0.166666666666667; 0.666666666666667]
 	elseif(order==3)
-		npoints=4
+		npoints = 4
 		weights = [-0.974278579257493; 0.902109795608790; 0.902109795608790; 0.902109795608790]
 		coords1 = [ 0.333333333333333; 0.600000000000000; 0.200000000000000; 0.200000000000000]
Index: /issm/trunk-jpl/src/jl/solve/modules.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/modules.jl	(revision 26658)
+++ /issm/trunk-jpl/src/jl/solve/modules.jl	(revision 26659)
@@ -69,5 +69,5 @@
 	count = 0
 	for i in 1:md.mesh.numberofvertices
-		push!(vertices,Vertex(i,x[i],y[i]))
+		push!(vertices,Vertex(i,x[i],y[i],0.))
 	end
 
Index: /issm/trunk-jpl/src/jl/solve/nodes.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/nodes.jl	(revision 26658)
+++ /issm/trunk-jpl/src/jl/solve/nodes.jl	(revision 26659)
@@ -124,5 +124,4 @@
 		count = GetDofList(nodes[i],doflist,count,setenum)
 	end
-	println(count," ",ndofs)
 	@assert count==ndofs
 
Index: /issm/trunk-jpl/src/jl/solve/solve.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/solve.jl	(revision 26658)
+++ /issm/trunk-jpl/src/jl/solve/solve.jl	(revision 26659)
@@ -14,4 +14,5 @@
 include("./modules.jl")
 include("./toolkits.jl")
+include("./utils.jl")
 include("./elementmatrix.jl")
 
Index: /issm/trunk-jpl/src/jl/solve/utils.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/utils.jl	(revision 26659)
+++ /issm/trunk-jpl/src/jl/solve/utils.jl	(revision 26659)
@@ -0,0 +1,5 @@
+function Matrix2x2Determinant(A::Matrix{Float64}) #{{{
+
+	return  A[1,1]*A[2,2]-A[2,1]*A[1,2]
+
+end#}}}
Index: /issm/trunk-jpl/src/jl/solve/vertices.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/vertices.jl	(revision 26658)
+++ /issm/trunk-jpl/src/jl/solve/vertices.jl	(revision 26659)
@@ -4,3 +4,23 @@
 	x::Float64
 	y::Float64
+	z::Float64
 end# }}}
+
+#vertices functions
+function GetVerticesCoordinates(vertices::Vector{Vertex}) #{{{
+
+	#Intermediaries
+	nbv = length(vertices)
+
+	#Allocate
+	xyz_list = Matrix{Float64}(undef,nbv,3)
+
+	#Assign value to xyz_list
+	for i in 1:nbv
+		xyz_list[i,1]=vertices[i].x
+		xyz_list[i,2]=vertices[i].y
+		xyz_list[i,3]=vertices[i].z
+	end
+
+	return xyz_list
+end #}}}
Index: /issm/trunk-jpl/src/jl/test101.jl
===================================================================
--- /issm/trunk-jpl/src/jl/test101.jl	(revision 26658)
+++ /issm/trunk-jpl/src/jl/test101.jl	(revision 26659)
@@ -6,4 +6,5 @@
 #md = ISSM.triangle(md,"../../test/Exp/Square.exp",180000.)
 md = ISSM.triangle(md,"../../test/Exp/Square.exp",250000.)
+#md = ISSM.triangle(md,"../../test/Exp/Square.exp",1000.)
 md = ISSM.setmask(md,"all","")