Index: /issm/trunk-jpl/src/jl/solve/analyses/stressbalanceanalysis.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/analyses/stressbalanceanalysis.jl	(revision 26680)
+++ /issm/trunk-jpl/src/jl/solve/analyses/stressbalanceanalysis.jl	(revision 26681)
@@ -77,29 +77,67 @@
 	#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)
 
 	#Prepare material object
-	material = Matice()
+	material = Matice(element)
 	
 	#Start integrating
 	gauss = GaussTria(2)
-	for i in 1:gauss.numgauss
+	for ig in 1:gauss.numgauss
 
 		Jdet = JacobianDeterminant(xyz_list, gauss)
 		NodalFunctionsDerivatives(element,dbasis,xyz_list,gauss)
 
-		H = GetInputValue(H_input, gauss, i)
+		H  = GetInputValue(H_input, gauss, ig)
+		mu = ViscositySSA(material, xyz_list, gauss, ig)
 
-		error("STOP")
-		vx = GetInputValue(vx_input, gauss, i)
-		vy = GetInputValue(vy_input, gauss, i)
+		for i in 1:numnodes
+			for j in 1:numnodes
+				Ke.values[2*i-1,2*j-1] += gauss.weights[ig]*Jdet*mu*H*(4*dbasis[j,1]*dbasis[i,1] + dbasis[j,2]*dbasis[i,2])
+				Ke.values[2*i-1,2*j  ] += gauss.weights[ig]*Jdet*mu*H*(2*dbasis[j,2]*dbasis[i,1] + dbasis[j,1]*dbasis[i,2])
+				Ke.values[2*i  ,2*j-1] += gauss.weights[ig]*Jdet*mu*H*(2*dbasis[j,1]*dbasis[i,2] + dbasis[j,2]*dbasis[i,1])
+				Ke.values[2*i  ,2*j  ] += gauss.weights[ig]*Jdet*mu*H*(4*dbasis[j,2]*dbasis[i,2] + dbasis[j,1]*dbasis[i,1])
+			end
+		end
 	end
 
+	return Ke
+end #}}}
+function CreatePVector(analysis::StressbalanceAnalysis,element::Tria)# {{{
 
+	#Internmediaries
+	numnodes = 3
+
+	#Initialize Element vectro and basis functions
+	pe = ElementVector(element.nodes)
+	basis = Vector{Float64}(undef,numnodes)
+
+	#Retrieve all inputs and parameters
+	xyz_list = GetVerticesCoordinates(element.vertices)
+	H_input  = GetInput(element, ThicknessEnum)
+	s_input  = GetInput(element, SurfaceEnum)
+	rho_ice  = FindParam(element, MaterialsRhoIceEnum)
+	g        = FindParam(element, ConstantsGEnum)
+
+	#Start integrating
+	gauss = GaussTria(2)
+	for ig in 1:gauss.numgauss
+
+		Jdet = JacobianDeterminant(xyz_list, gauss)
+		NodalFunctions(element, basis, gauss, ig, P1Enum)
+
+		H  = GetInputValue(H_input, gauss, ig)
+		ds = GetInputDerivativeValue(s_input, xyz_list, gauss, ig)
+
+		for i in 1:numnodes
+			pe.values[2*i-1] += -gauss.weights[ig]*Jdet*rho_ice*g*H*ds[1]*basis[i]
+			pe.values[2*i  ] += -gauss.weights[ig]*Jdet*rho_ice*g*H*ds[2]*basis[i]
+		end
+	end
+
+	print(pe)
 	error("STOP")
 
-	return Ke
+	return Pe
 end #}}}
 function GetSolutionFromInputs(analysis::StressbalanceAnalysis,ug::Vector{Float64},element::Tria) #{{{
Index: /issm/trunk-jpl/src/jl/solve/elementmatrix.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/elementmatrix.jl	(revision 26680)
+++ /issm/trunk-jpl/src/jl/solve/elementmatrix.jl	(revision 26681)
@@ -43,5 +43,5 @@
 
 	#Initialize element matrix with zeros
-	values = zeros(nrows,nrows)
+	values = zeros(nrows)
 
 	#Get dof list
@@ -54,7 +54,5 @@
 	println(io,"ElementVector:")
 	println(io,"   nrows: ",this.nrows)
-	println(io,"   gglobaldoflist: ",this.gglobaldoflist)
 	println(io,"   fglobaldoflist: ",this.fglobaldoflist)
-	println(io,"   sglobaldoflist: ",this.sglobaldoflist)
 	print(io,"   values: ")
 	display(this.values)
Index: /issm/trunk-jpl/src/jl/solve/elements.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/elements.jl	(revision 26680)
+++ /issm/trunk-jpl/src/jl/solve/elements.jl	(revision 26681)
@@ -89,4 +89,9 @@
 
 end # }}}
+function FindParam(element::Tria,enum::IssmEnum) # {{{
+
+	return FindParam(element.parameters, enum)
+
+end # }}}
 function InputServe!(element::Tria,input::ElementInput) # {{{
 
@@ -145,9 +150,24 @@
 	return Jinv
 end#}}}
+function NodalFunctions(element::Tria,basis::Vector{Float64}, gauss::GaussTria, ig::Int64, finiteelement::IssmEnum) #{{{
+
+	if(finiteelement==P0Enum)
+		#Nodal function 1
+		basis[1]= 1
+	elseif(finiteelement==P1Enum)
+		basis[1] = gauss.coords1[ig]
+		basis[2] = gauss.coords2[ig]
+		basis[3] = gauss.coords3[ig]
+	else
+		error("Element type ",finiteelement," not supported yet")
+	end
+
+
+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)
+	NodalFunctionsDerivativesReferenceTria(dbasis_ref,gauss,P1Enum)
 
 	#Get invert of the Jacobian
@@ -163,10 +183,10 @@
 
 end#}}}
-function NodalFunctionsDerivativesReference(dbasis::Matrix{Float64}, gauss::GaussTria, finiteelement::IssmEnum) #{{{
+function NodalFunctionsDerivativesReferenceTria(dbasis::Matrix{Float64}, gauss::GaussTria, finiteelement::IssmEnum) #{{{
 
 	if(finiteelement==P0Enum)
 		#Nodal function 1
-		dbasis[1,1]= 0.
-		dbasis[1,2]= 0.
+		dbasis[1,1]= 0
+		dbasis[1,2]= 0
 
 	elseif(finiteelement==P1Enum)
@@ -184,2 +204,10 @@
 	end
 end#}}}
+function NumberofNodesTria(finiteelement) #{{{
+
+	if    (finiteelement==P0Enum) return 0
+	elseif(finiteelement==P1Enum) return 3
+	else
+		error("Element type ",finiteelement," not supported yet")
+	end
+end#}}}
Index: /issm/trunk-jpl/src/jl/solve/inputs.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/inputs.jl	(revision 26680)
+++ /issm/trunk-jpl/src/jl/solve/inputs.jl	(revision 26681)
@@ -97,8 +97,20 @@
 	input.values[indices] = values
 end#}}}
-function GetInputValue(input::ElementInput,gauss::GaussTria,i::Int64)
+function GetInputAverageValue(input::ElementInput,gauss::GaussTria,i::Int64) #{{{
+
+	numnodes = NumberofNodesTria(input.interp)
+	value = 0.0
+
+	for i in 1:numnodes
+		value+=input.element_values[i]
+	end
+
+	return value/numnodes
+
+end#}}}
+function GetInputValue(input::ElementInput,gauss::GaussTria,i::Int64) #{{{
 
 	if input.interp==P0Enum
-		return input.element_value
+		return input.element_values
 	elseif input.interp==P1Enum
 		value = input.element_values[1]*gauss.coords1[i] +  input.element_values[2]*gauss.coords2[i] +  input.element_values[3]*gauss.coords3[i]
@@ -110,2 +122,31 @@
 
 end#}}}
+function GetInputDerivativeValue(input::ElementInput,xyz_list::Matrix{Float64},gauss::GaussTria,i::Int64) #{{{
+
+	#Get nodal function derivatives in reference element
+	numnodes = NumberofNodesTria(input.interp)
+	dbasis_ref = Matrix{Float64}(undef,numnodes,2)
+	NodalFunctionsDerivativesReferenceTria(dbasis_ref,gauss,input.interp)
+
+	#Get invert of the Jacobian
+	Jinv = JacobianInvert(xyz_list,gauss)
+
+	#Build dbasis:
+	#[ dNi/dx ] = Jinv * [dNhat_i/dr]
+	#[ dNi/dy ] =        [dNhat_i/ds]
+	dbasis = Matrix{Float64}(undef,numnodes,2)
+	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
+
+	#Get derivatives: dp/dx dp/dy
+	dp = [0.0;0.0]
+	for i in 1:3
+		dp[1] += dbasis[i,1]*input.element_values[i]
+		dp[2] += dbasis[i,2]*input.element_values[i]
+	end
+
+	return dp
+
+end#}}}
Index: /issm/trunk-jpl/src/jl/solve/matice.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/matice.jl	(revision 26680)
+++ /issm/trunk-jpl/src/jl/solve/matice.jl	(revision 26681)
@@ -1,10 +1,44 @@
 #Matice class definition
-mutable struct Matice#{{{
+struct Matice#{{{
+	vx_input::ElementInput
+	vy_input::ElementInput
+	B_input::ElementInput
+	n_input::ElementInput
 end# }}}
 
+function Matice(element::Tria) #{{{
+
+	vx_input  = GetInput(element, VxEnum)
+	vy_input  = GetInput(element, VyEnum)
+	B_input   = GetInput(element, MaterialsRheologyBEnum)
+	n_input   = GetInput(element, MaterialsRheologyNEnum)
+
+	return Matice(vx_input, vy_input, B_input, n_input)
+end#}}}
+
 #vertices functions
-function GetViscositySSA() #{{{
+function ViscositySSA(matice::Matice, xyz_list::Matrix{Float64}, gauss::GaussTria, i::Int64) #{{{
 
-	error("not implemented")
+	#Get strain rate
+	dvx = GetInputDerivativeValue(matice.vx_input,xyz_list,gauss,i)
+	dvy = GetInputDerivativeValue(matice.vy_input,xyz_list,gauss,i)
+	eps_xx = dvx[1]
+	eps_yy = dvy[2]
+	eps_xy = 0.5*(dvx[2] + dvy[1])
+
+	#In SSA, eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy
+	eps_eff = sqrt(eps_xx*eps_xx + eps_yy*eps_yy + eps_xy*eps_xy + eps_xx*eps_yy)
+
+	#Get B and n
+	n = GetInputValue(matice.n_input, gauss, i)
+	B = GetInputValue(matice.B_input, gauss, i)
+
+	#Compute viscosity
+	if eps_eff==0.
+		mu = 1.e+14/2
+	else
+		mu = B/(2*eps_eff^((n-1)/n))
+	end
+
 	return mu
 end #}}}
Index: /issm/trunk-jpl/src/jl/solve/modules.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/modules.jl	(revision 26680)
+++ /issm/trunk-jpl/src/jl/solve/modules.jl	(revision 26681)
@@ -79,4 +79,5 @@
 	AddParam(parameters,md.materials.rho_water,MaterialsRhoSeawaterEnum)
 	AddParam(parameters,md.materials.rho_freshwater,MaterialsRhoFreshwaterEnum)
+	AddParam(parameters,md.constants.g,ConstantsGEnum)
 end# }}}
 function CreateInputs(inputs::Inputs,elements::Vector{Tria},md::model) #{{{