Index: /issm/trunk-jpl/src/jl/solve/analyses/stressbalanceanalysis.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/analyses/stressbalanceanalysis.jl	(revision 26709)
+++ /issm/trunk-jpl/src/jl/solve/analyses/stressbalanceanalysis.jl	(revision 26710)
@@ -177,5 +177,38 @@
 
 	if(IsIcefront(element))
-		error("not supported")
+		A = zeros(6)
+
+		#Get additional parameters and inputs
+		b_input   = GetInput(element, BaseEnum)
+		rho_water = FindParam(element, MaterialsRhoSeawaterEnum)
+
+		#Get normal and ice front coordinates
+		xyz_list_front = Matrix{Float64}(undef,2,3)
+		GetIcefrontCoordinates!(element, xyz_list_front, xyz_list, MaskIceLevelsetEnum)
+		nx, ny = NormalSection(element, xyz_list_front)
+		#println("nx ",nx," ny:",ny)
+
+		gauss = GaussTria(element, xyz_list, xyz_list_front, 3)
+		for ig in 1:gauss.numgauss
+
+			Jdet = JacobianDeterminantSurface(xyz_list_front, gauss)
+			NodalFunctions(element, basis, gauss, ig, P1Enum)
+
+			H  = GetInputValue(H_input, gauss, ig)
+			b  = GetInputValue(b_input, gauss, ig)
+			sl = 0
+
+			term = 0.5*g*rho_ice*H^2 + 0.5*g*rho_water*(min(0, H+b-sl)^2 - min(0, b-sl)^2)
+
+			for i in 1:numnodes
+				#println("1: ",gauss.weights[ig]*Jdet*term*nx*basis[i]," 2: ",gauss.weights[ig]*Jdet*term*ny*basis[i])
+				pe.values[2*i-1] += gauss.weights[ig]*Jdet*term*nx*basis[i]
+				pe.values[2*i  ] += gauss.weights[ig]*Jdet*term*ny*basis[i]
+				A[2*i-1] += gauss.weights[ig]*Jdet*term*nx*basis[i]
+				A[2*i  ] += gauss.weights[ig]*Jdet*term*ny*basis[i]
+			end
+		end
+		#error("S")
+		display(A)
 	end
 
Index: /issm/trunk-jpl/src/jl/solve/elements.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/elements.jl	(revision 26709)
+++ /issm/trunk-jpl/src/jl/solve/elements.jl	(revision 26710)
@@ -176,4 +176,60 @@
 	end
 end#}}}
+function GetIcefrontCoordinates!(element::Tria, xyz_front::Matrix{Float64}, xyz_list::Matrix{Float64}, levelsetenum::IssmEnum) #{{{
+
+	#Intermediaries
+	level        = Vector{Float64}(undef,3)
+	indicesfront = Vector{Int64}(undef,3)
+
+	#Recover value of levelset for all vertices
+	GetInputListOnVertices!(element, level, levelsetenum)
+
+	#Get nodes where there is no ice
+	num_frontnodes = 0
+	for i in 1:3
+		if(level[i]>=0.)
+			num_frontnodes += 1
+			indicesfront[num_frontnodes] = i
+		end
+	end
+	@assert num_frontnodes==2
+
+	#Arrange order of frontnodes such that they are oriented counterclockwise
+	NUMVERTICES = 3
+	if((NUMVERTICES+indicesfront[1]-indicesfront[2])%NUMVERTICES != NUMVERTICES-1)
+		index=indicesfront[1]
+		indicesfront[1]=indicesfront[2]
+		indicesfront[2]=index
+	end
+
+	#Return nodes
+	xyz_front[1,:]=xyz_list[indicesfront[1],:]
+	xyz_front[2,:]=xyz_list[indicesfront[2],:]
+
+end#}}}
+function GetArea(element::Tria)#{{{
+
+	#Get xyz list
+	xyz_list = GetVerticesCoordinates(element.vertices)
+	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]
+
+	@assert x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1>0
+	return (x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1)/2
+end#}}}
+function NormalSection(element::Tria, xyz_front::Matrix{Float64}) #{{{
+
+	#Build output pointing vector
+	nx =  xyz_front[2,2] - xyz_front[1,2]
+	ny = -xyz_front[2,1] + xyz_front[1,1]
+
+	#normalize
+	norm = sqrt(nx^2 + ny^2)
+	nx = nx/norm
+	ny = ny/norm
+
+	return nx, ny
+end#}}}
 
 #Finite Element stuff
@@ -185,4 +241,15 @@
 	#Get its determinant
 	Jdet = Matrix2x2Determinant(J)
+
+	#check and return
+	if(Jdet<0) error("negative Jacobian Determinant") end
+	return Jdet
+
+end#}}}
+function JacobianDeterminantSurface(xyz_list::Matrix{Float64}, gauss::GaussTria) #{{{
+
+	x1 = xyz_list[1,1]; y1 = xyz_list[1,2]
+	x2 = xyz_list[2,1]; y2 = xyz_list[2,2]
+	Jdet = .5*sqrt((x2-x1)^2 + (y2-y1)^2)
 
 	#check and return
@@ -281,2 +348,31 @@
 	end
 end#}}}
+function GaussTria(element::Tria, xyz_list::Matrix{Float64}, xyz_list_front::Matrix{Float64}, order::Int64) #{{{
+
+	area_coordinates = Matrix{Float64}(undef,2,3)
+	GetAreaCoordinates!(element, area_coordinates, xyz_list_front, xyz_list)
+
+	return GaussTria(area_coordinates, order)
+end# }}}
+function GetAreaCoordinates!(element::Tria, area_coordinates::Matrix{Float64}, xyz_zero::Matrix{Float64}, xyz_list::Matrix{Float64})#{{{
+
+	numpoints = size(area_coordinates,1)
+	area = GetArea(element)
+
+	#Copy original xyz_list
+	xyz_bis=copy(xyz_list)
+	for i in 1:numpoints
+		for j in 1:3
+
+			#Change appropriate line
+			xyz_bis[j,:] = xyz_zero[i,:]
+
+			#Compute area fraction
+			area_portion=abs(xyz_bis[2,1]*xyz_bis[3,2] - xyz_bis[2,2]*xyz_bis[3,1] + xyz_bis[1,1]*xyz_bis[2,2] - xyz_bis[1,2]*xyz_bis[2,1] + xyz_bis[3,1]*xyz_bis[1,2] - xyz_bis[3,2]*xyz_bis[1,1])/2
+			area_coordinates[i,j] = area_portion/area
+
+			#reinitialize xyz_list
+			xyz_bis[j,:] = xyz_list[j,:]
+		end
+	end
+end #}}}
Index: /issm/trunk-jpl/src/jl/solve/gauss.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/gauss.jl	(revision 26709)
+++ /issm/trunk-jpl/src/jl/solve/gauss.jl	(revision 26710)
@@ -7,4 +7,13 @@
 	coords3::Vector{Float64}
 end #}}}
+function Base.show(io::IO, this::GaussTria)# {{{
+
+	println(io,"GaussTria:")
+	println(io,"   numgauss: ",this.numgauss)
+	println(io,"   weights:  ",this.weights)
+	println(io,"   coords1:  ",this.coords1)
+	println(io,"   coords2:  ",this.coords2)
+	println(io,"   coords3:  ",this.coords3)
+end# }}}
 
 #Gauss constructor
@@ -60,2 +69,45 @@
 	return GaussTria(npoints,weights,coords1,coords2,coords3)
 end# }}}
+function GaussTria(area_coordinates::Matrix{Float64}, order::Int64) #{{{
+	#=Gauss-Legendre quadrature points.
+
+	The recurrence coefficients for Legendre polynomials on (-1,1)
+	are defined (from the ORTHPOL subroutine RECUR with ipoly=1) as:
+
+	alpha(i)=0.
+	beta (i)=1./(4.-1./(i-1)^2))
+
+	For degree p, the required number of Gauss-Legendre points is
+	n>=(p+1)/2.=#
+
+	if(order==1)
+		npoint  = 1
+		weights = [2.000000000000000]
+		coords  = [0.000000000000000]
+	elseif(order==2)
+		npoints = 2
+		weights = [1.000000000000000, 1.000000000000000]
+		coords  = [-0.577350269189626, 0.577350269189626]
+	elseif(order==3)
+		npoints = 3
+		weights = [0.555555555555556, 0.888888888888889, 0.555555555555556]
+		coords  = [-0.774596669241483, 0.000000000000000, 0.774596669241483]
+	elseif(order==4)
+		npoints = 4
+		weights = [0.347854845137454, 0.652145154862546, 0.652145154862546, 0.347854845137454]
+		coords  = [-0.861136311594053,-0.339981043584856, 0.339981043584856, 0.861136311594053]
+	else
+      error("order ",order," not supported yet");
+	end
+
+   coords1  = Vector{Float64}(undef,npoints)
+   coords2  = Vector{Float64}(undef,npoints)
+   coords3  = Vector{Float64}(undef,npoints)
+   for i in 1:npoints
+      coords1[i]  = 0.5*(area_coordinates[1,1]+area_coordinates[2,1]) + 0.5*coords[i]*(area_coordinates[2,1]-area_coordinates[1,1]);
+      coords2[i]  = 0.5*(area_coordinates[1,2]+area_coordinates[2,2]) + 0.5*coords[i]*(area_coordinates[2,2]-area_coordinates[1,2]);
+      coords3[i]  = 0.5*(area_coordinates[1,3]+area_coordinates[2,3]) + 0.5*coords[i]*(area_coordinates[2,3]-area_coordinates[1,3]);
+   end
+
+	return GaussTria(npoints, weights, coords1, coords2, coords3)
+end# }}}
Index: /issm/trunk-jpl/src/jl/solve/solutionsequences.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/solutionsequences.jl	(revision 26709)
+++ /issm/trunk-jpl/src/jl/solve/solutionsequences.jl	(revision 26710)
@@ -35,5 +35,5 @@
 		#Increase count
 		count += 1
-		if(count>maxiter)
+		if(count>=maxiter)
 			println("   maximum number of nonlinear iterations (",maxiter,") exceeded")
 			converged = true