Index: /issm/trunk/src/m/solvers/solver_stokescoupling_nonlinear.m
===================================================================
--- /issm/trunk/src/m/solvers/solver_stokescoupling_nonlinear.m	(revision 6022)
+++ /issm/trunk/src/m/solvers/solver_stokescoupling_nonlinear.m	(revision 6023)
@@ -14,7 +14,5 @@
 	femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum);
 	ug_horiz=GetSolutionFromInputs(femmodel.elements, femmodel.nodes, femmodel.vertices, femmodel.loads, femmodel.materials, femmodel.parameters);
-
-	femmodel=SetCurrentConfiguration(femmodel,DiagnosticVertAnalysisEnum);
-	ug_vert=GetSolutionFromInputs(femmodel.elements, femmodel.nodes, femmodel.vertices, femmodel.loads, femmodel.materials, femmodel.parameters);
+	uf_horiz=Reducevectorgtof( ug_horiz, femmodel.nodesets,femmodel.parameters);
 
 	while(~converged),
@@ -22,8 +20,4 @@
 		%First compute the horizontal velocity
 		femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum);
-
-		%Start non-linear iteration using input velocity: 
-		ug_horiz=GetSolutionFromInputs(femmodel.elements, femmodel.nodes, femmodel.vertices, femmodel.loads, femmodel.materials, femmodel.parameters);
-		uf_horiz=Reducevectorgtof( ug_horiz, femmodel.nodesets,femmodel.parameters);
 
 		%Update the solution to make sure that vx and vxold are similar
@@ -35,13 +29,13 @@
 
 		if kffpartitioning,
-			[K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
-			p_f = Reduceload( p_f, K_fs, femmodel.ys,femmodel.parameters);
+			[K_gg_horiz,K_ff_horiz,K_fs_horiz,p_g_horiz,p_f_horiz,kmax_horiz]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
+			p_f_horiz = Reduceload( p_f_horiz, K_fs_horiz, femmodel.ys,femmodel.parameters);
 		else
-			[K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
-			[K_ff, K_fs] = Reducematrixfromgtof( K_gg, femmodel.nodesets,femmodel.parameters); 
-			p_f = Reduceloadfromgtof( p_g, K_fs, femmodel.ys, femmodel.nodesets,femmodel.parameters);
+			[K_gg_horiz,K_ff_horiz,K_fs_horiz,p_g_horiz,p_f_horiz,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
+			[K_ff_horiz, K_fs_horiz] = Reducematrixfromgtof( K_gg_horiz, femmodel.nodesets,femmodel.parameters); 
+			p_f_horiz = Reduceloadfromgtof( p_g_horiz, K_fs_horiz, femmodel.ys, femmodel.nodesets,femmodel.parameters);
 		end
 
-		uf_horiz=Solver(K_ff,p_f,old_uf_horiz,femmodel.parameters);
+		uf_horiz=Solver(K_ff_horiz,p_f_horiz,old_uf_horiz,femmodel.parameters);
 		ug_horiz= Mergesolutionfromftog( uf_horiz, femmodel.ys, femmodel.nodesets,femmodel.parameters); 
 
@@ -49,24 +43,23 @@
 
 		%Figure out if convergence have been reached
-		converged=convergence(K_ff,p_f,uf_horiz,old_uf_horiz,femmodel.parameters);
-		converged=1;
+		converged=convergence(K_ff_horiz,p_f_horiz,uf_horiz,old_uf_horiz,femmodel.parameters);
 
-%		%Then compute vertical velocity
-%		femmodel=SetCurrentConfiguration(femmodel,DiagnosticVertAnalysisEnum);
-%	if kffpartitioning,
-%		[K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
-%		p_f = Reduceload( p_f, K_fs, femmodel.ys,femmodel.parameters);
-%	else
-%		[K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
-%		[K_ff, K_fs] = Reducematrixfromgtof( K_gg, femmodel.nodesets,femmodel.parameters); 
-%		p_f = Reduceloadfromgtof( p_g, K_fs, femmodel.ys, femmodel.nodesets,femmodel.parameters);
-%	end
-%
-%		displaystring(femmodel.parameters.Verbose>1,'%s%g','      condition number of stiffness matrix: ',condest(K_ff));
-%		uf_vert=Solver(K_ff,p_f,[],femmodel.parameters);
-%
-%		ug_vert= Mergesolutionfromftog( uf_vert, femmodel.ys, femmodel.nodesets,femmodel.parameters); 
-%
-%		[femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ug_vert);
+		%Then compute vertical velocity
+		femmodel=SetCurrentConfiguration(femmodel,DiagnosticVertAnalysisEnum);
+
+		if kffpartitioning,
+			[K_gg_vert,K_ff_vert,K_fs_vert,p_g_vert,p_f_vert,kmax_vert]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
+			p_f_vert = Reduceload( p_f_vert, K_fs_vert, femmodel.ys,femmodel.parameters);
+		else
+			[K_gg_vert,K_ff_vert,K_fs_vert,p_g_vert,p_f_vert,kmax_vert]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
+			[K_ff_vert, K_fs_vert] = Reducematrixfromgtof( K_gg_vert, femmodel.nodesets,femmodel.parameters); 
+			p_f_vert = Reduceloadfromgtof( p_g_vert, K_fs_vert, femmodel.ys, femmodel.nodesets,femmodel.parameters);
+		end
+
+		uf_vert=Solver(K_ff_vert,p_f_vert,[],femmodel.parameters);
+
+		ug_vert= Mergesolutionfromftog( uf_vert, femmodel.ys, femmodel.nodesets,femmodel.parameters); 
+
+		[femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ug_vert);
 
 	end