Index: /issm/trunk-jpl/src/jl/solve/analyses/stressbalanceanalysis.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/analyses/stressbalanceanalysis.jl	(revision 26674)
+++ /issm/trunk-jpl/src/jl/solve/analyses/stressbalanceanalysis.jl	(revision 26675)
@@ -81,4 +81,7 @@
 	H_input  = GetInput(element, ThicknessEnum)
 
+	#Prepare material object
+	material = Matice()
+	
 	#Start integrating
 	gauss = GaussTria(2)
@@ -86,6 +89,7 @@
 
 		Jdet = JacobianDeterminant(xyz_list, gauss)
+		NodalFunctionsDerivatives(element,dbasis,xyz_list,gauss)
 
-		println(Jdet)
+		H = GetInputValue(H_input, gauss, i)
 
 		error("STOP")
Index: /issm/trunk-jpl/src/jl/solve/elements.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/elements.jl	(revision 26674)
+++ /issm/trunk-jpl/src/jl/solve/elements.jl	(revision 26675)
@@ -129,8 +129,57 @@
 
 	J[1,1] = .5*(x2-x1)
-	J[1,2] = sqrt(3)/6*(2*x3 -x1 -x2)
-	J[2,1] = .5*(y2-y1)
+	J[1,2] = .5*(y2-y1)
+	J[2,1] = sqrt(3)/6*(2*x3 -x1 -x2)
 	J[2,2] = sqrt(3)/6*(2*y3 -y1 -y2)
 
 	return J
 end#}}}
+function JacobianInvert(xyz_list::Matrix{Float64}, gauss::GaussTria) #{{{
+
+	#Get Jacobian matrix
+	J = Jacobian(xyz_list)
+
+	#Get its determinant
+	Jinv = Matrix2x2Invert(J)
+
+	return Jinv
+end#}}}
+function NodalFunctionsDerivatives(element::Tria,dbasis::Matrix{Float64},xyz_list::Matrix{Float64}, gauss::GaussTria) #{{{
+
+	#Get nodal function derivatives in reference element
+	dbasis_ref = Matrix{Float64}(undef,3,2)
+	NodalFunctionsDerivativesReference(dbasis_ref,gauss,P1Enum)
+
+	#Get invert of the Jacobian
+	Jinv = JacobianInvert(xyz_list,gauss)
+
+	#Build dbasis:
+	#[ dNi/dx ] = Jinv * [dNhat_i/dr]
+	#[ dNi/dy ] =        [dNhat_i/ds]
+	for i in 1:3
+		dbasis[i,1] = Jinv[1,1]*dbasis_ref[i,1]+Jinv[1,2]*dbasis_ref[i,2]
+		dbasis[i,2] = Jinv[2,1]*dbasis_ref[i,1]+Jinv[2,2]*dbasis_ref[i,2]
+	end
+
+end#}}}
+function NodalFunctionsDerivativesReference(dbasis::Matrix{Float64}, gauss::GaussTria, finiteelement::IssmEnum) #{{{
+
+	if(finiteelement==P0Enum)
+		#Nodal function 1
+		dbasis[1,1]= 0.
+		dbasis[1,2]= 0.
+
+	elseif(finiteelement==P1Enum)
+		#Nodal function 1
+		dbasis[1,1]= -.5
+		dbasis[1,2]= -sqrt(3)/6
+		#Nodal function 2
+		dbasis[2,1]= .5
+		dbasis[2,2]= -sqrt(3)/6
+		#Nodal function 3
+		dbasis[3,1]= 0
+		dbasis[3,2]= sqrt(3)/3
+	else
+		error("Element type ",finiteelement," not supported yet")
+	end
+end#}}}
Index: /issm/trunk-jpl/src/jl/solve/matice.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/matice.jl	(revision 26675)
+++ /issm/trunk-jpl/src/jl/solve/matice.jl	(revision 26675)
@@ -0,0 +1,10 @@
+#Matice class definition
+mutable struct Matice#{{{
+end# }}}
+
+#vertices functions
+function GetViscositySSA() #{{{
+
+	error("not implemented")
+	return mu
+end #}}}
Index: /issm/trunk-jpl/src/jl/solve/solve.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/solve.jl	(revision 26674)
+++ /issm/trunk-jpl/src/jl/solve/solve.jl	(revision 26675)
@@ -8,4 +8,5 @@
 include("./elements.jl")
 include("./constraints.jl")
+include("./matice.jl")
 include("./femmodel.jl")
 include("./analyses/analysis.jl")
Index: /issm/trunk-jpl/src/jl/solve/utils.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/utils.jl	(revision 26674)
+++ /issm/trunk-jpl/src/jl/solve/utils.jl	(revision 26675)
@@ -4,2 +4,23 @@
 
 end#}}}
+function Matrix2x2Invert(A::Matrix{Float64}) #{{{
+
+	#Initialize output
+	Ainv = Matrix{Float64}(undef,2,2)
+
+	#Compute determinant
+	det = Matrix2x2Determinant(A)
+	if(abs(det)<eps(Float64)) error("Determinant smaller than machine epsilon") end
+
+	#Multiplication is faster than divsion, so we multiply by the reciprocal
+	det_reciprocal = 1/det
+
+	#compute invert matrix
+	Ainv[1,1]=   A[2,2]*det_reciprocal # =  d/det
+   Ainv[1,2]= - A[1,2]*det_reciprocal # = -b/det
+   Ainv[2,1]= - A[2,1]*det_reciprocal # = -c/det
+   Ainv[2,2]=   A[1,1]*det_reciprocal # =  a/det
+
+	return  Ainv
+
+end#}}}