Index: /issm/trunk-jpl/src/jl/solve/inputs.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/inputs.jl	(revision 26698)
+++ /issm/trunk-jpl/src/jl/solve/inputs.jl	(revision 26699)
@@ -112,5 +112,5 @@
 
 	if input.interp==P0Enum
-		return input.element_values
+		return input.element_values[1]
 	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]
Index: /issm/trunk-jpl/src/jl/solve/modules.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/modules.jl	(revision 26698)
+++ /issm/trunk-jpl/src/jl/solve/modules.jl	(revision 26699)
@@ -158,5 +158,5 @@
 	end
 
-	return ug
+	return uf
 
 end#}}}
Index: /issm/trunk-jpl/src/jl/solve/solutionsequences.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/solutionsequences.jl	(revision 26698)
+++ /issm/trunk-jpl/src/jl/solve/solutionsequences.jl	(revision 26699)
@@ -29,6 +29,7 @@
 		Mergesolutionfromftogx(ug, uf, ys, femmodel.nodes)
 
-		print(ug)
-		error("compare with ISSM...")
+		#Check for convergence
+		converged = convergence(Kff,pf,uf,old_uf,restol,reltol,abstol)
+		InputUpdateFromSolutionx(analysis,ug,femmodel)
 
 		#Increase count
@@ -40,5 +41,68 @@
 	end
 
+	print("\n   total number of iterations: ",  count,  "\n")
+
 	error("STOP")
 
 end# }}}
+function convergence(Kff::IssmMatrix, pf::IssmVector, uf::IssmVector, old_uf::IssmVector, restol::Float64, reltol::Float64, abstol::Float64)#{{{
+
+	print("   checking convergence\n");
+
+	#If solution vector is empty, return true
+	if(IsEmpty(uf))
+		return true
+	end
+
+	#Convergence criterion #1: force equilibrium (Mandatory)
+	#compute K[n]U[n-1] - F
+	KUold  = Duplicate(uf);    MatMult!(Kff,old_uf,KUold)
+	KUoldF = Duplicate(KUold); VecCopy!(KUold, KUoldF); AXPY!(KUoldF, -1.0, pf)
+	nKUoldF = Norm(KUoldF,2)
+	nF      = Norm(pf,2)
+	res = nKUoldF/nF
+	if ~isfinite(res)
+		println("norm nf = ", nF, " and norm kuold = ",nKUoldF)
+		error("mechanical equilibrium convergence criterion is not finite!")
+	end
+	if(res<restol)
+		print("   mechanical equilibrium convergence criterion ", res*100, " < ", restol*100, " %\n")
+		converged=true
+	else
+		print("   mechanical equilibrium convergence criterion ", res*100, " > ", restol*100, " %\n")
+		converged=false;
+	end
+
+	#Convergence criterion #2: norm(du)/norm(u)
+	if ~isnan(reltol)
+		duf = Duplicate(old_uf); VecCopy!(old_uf,duf); AXPY!(duf, -1.0, uf)
+		ndu = Norm(duf, 2); nu = Norm(old_uf, 2)
+		if ~isfinite(ndu) | ~isfinite(nu) 
+			error("convergence criterion is not finite!")
+		end
+		if((ndu/nu)<reltol)
+			print("   Convergence criterion: norm(du)/norm(u)      ", ndu/nu*100, " < ", reltol*100, " %\n")
+		else
+			print("   Convergence criterion: norm(du)/norm(u)      ", ndu/nu*100, " > ", reltol*100, " %\n")
+			converged=false;
+		end
+	end
+
+	#Convergence criterion #3: max(du)
+	if ~isnan(abstol)
+		duf = Duplicate(old_uf); VecCopy!(old_uf,duf); AXPY!(duf, -1.0, uf)
+		nduinf= Norm(duf, 3)
+		if ~isfinite(nduinf) 
+			error("convergence criterion is not finite!")
+		end
+		if(nduinf<abstol)
+			print("   Convergence criterion: max(du)               ", nduinf, " < ", abstol, "\n")
+		else
+			print("   Convergence criterion: max(du)               ", nduinf, " > ", abstol, "\n")
+			converged=false;
+		end
+	end
+
+	return converged
+
+end#}}}
Index: /issm/trunk-jpl/src/jl/solve/toolkits.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve/toolkits.jl	(revision 26698)
+++ /issm/trunk-jpl/src/jl/solve/toolkits.jl	(revision 26699)
@@ -58,4 +58,42 @@
 
 end#}}}
+function IsEmpty(vector::IssmVector)#{{{
+
+	return GetSize(vector)==0
+
+end#}}}
+function Duplicate(vector::IssmVector)#{{{
+
+	#Copy data structure
+	M=GetSize(vector)
+	return IssmVector(M)
+
+end#}}}
+function VecCopy!(x::IssmVector,y::IssmVector)#{{{
+
+	y.vector = x.vector
+
+end#}}}
+function Norm(x::IssmVector,type::Int64)#{{{
+
+	norm = 0
+
+	if type==2
+		for i in 1:length(x.vector)
+			norm += x.vector[i]^2
+		end
+		norm = sqrt(norm)
+	elseif type==3
+		#Infinite norm
+		for i in 1:length(x.vector)
+			if(abs(x.vector[i])>norm) norm = abs(x.vector[i]) end
+		end
+	else
+		error("type ",type," not supported yet")
+	end
+
+	return norm
+
+end#}}}
 
 #Operations