Index: /issm/trunk/src/m/solvers/solver_diagnostic_nonlinear.m
===================================================================
--- /issm/trunk/src/m/solvers/solver_diagnostic_nonlinear.m	(revision 4114)
+++ /issm/trunk/src/m/solvers/solver_diagnostic_nonlinear.m	(revision 4115)
@@ -1,7 +1,7 @@
-function [femmodel ug varargout]=diagnostic_core_nonlinear(femmodel,conserve_loads)
-%DIAGNOSTIC_CORE_NONLINEAR - core solution of diagnostic non-linear
+function [femmodel ug varargout]=solver_diagnostic_nonlinear(femmodel,conserve_loads)
+%SOLVER_DIAGNOSTIC_NONLINEAR - core solver of diagnostic run
 %
 %   Usage:
-%      [femmodel ug varargout]=diagnostic_core_nonlinear(femmodel,conserve_loads)
+%      [femmodel ug varargout]=solver_diagnostic_nonlinear(femmodel,conserve_loads)
 	
 
@@ -44,5 +44,5 @@
 
 		%Update inputs using new solution
-		femmodel.elements=UpdateInputsFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters,ug);
+		femmodel.elements=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters,ug);
 		
 		displaystring(femmodel.parameters.Verbose,'%s','   penalty constraints');
Index: /issm/trunk/src/m/solvers/solver_linear.m
===================================================================
--- /issm/trunk/src/m/solvers/solver_linear.m	(revision 4114)
+++ /issm/trunk/src/m/solvers/solver_linear.m	(revision 4115)
@@ -1,27 +1,25 @@
-function u_g=diagnostic_core_linear(m,analysis_type,sub_analysis_type)
-%DIAGNOSTIC_CORE_LINEAR - linear solution sequence
+function [femmodel u_g]=solver_linear(femmodel)
+%SOLVER_LINEAR - core solver of any linear solution sequence
 %
 %   Usage:
-%      u_g=diagnostic_core_linear(m,analysis_type,sub_analysis_type)
+%      [femmodel, u_g]=solver_linear(femmodel)
 
 	%stiffness and load generation only:
-	m.parameters.Kflag=1; m.parameters.Pflag=1;
+	femmodel.parameters.Kflag=1; femmodel.parameters.Pflag=1;
 
 	%system matrices
-	[K_gg, p_g]=SystemMatrices(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,analysis_type,sub_analysis_type);
-	[K_gg, p_g,kmax]=PenaltySystemMatrices(K_gg,p_g,m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,analysis_type,sub_analysis_type);
+	[K_gg, p_g]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
+	[K_gg, p_g,kmax]=PenaltySystemMatrices(K_gg,p_g,femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
 	
 	%Reduce tangent matrix from g size to f size
-	[K_ff, K_fs] = Reducematrixfromgtof( K_gg, m.Gmn, m.nodesets); 
-	displaystring(m.parameters.Verbose>1,'%s%g','      condition number of stiffness matrix: ',condest(K_ff));
+	[K_ff, K_fs] = Reducematrixfromgtof( K_gg, femmodel.Gmn, femmodel.nodesets); 
+	displaystring(femmodel.parameters.Verbose>1,'%s%g','      condition number of stiffness matrix: ',condest(K_ff));
 	
 	%Reduce load from g size to f size
-	[p_f] = Reduceloadfromgtof( p_g, m.Gmn, K_fs, m.ys, m.nodesets);
+	p_f = Reduceloadfromgtof( p_g, femmodel.Gmn, K_fs, femmodel.ys, femmodel.nodesets);
 
 	%Solve	
-	[u_f]=Solver(K_ff,p_f,[],m.parameters);
+	u_f=Solver(K_ff,p_f,[],femmodel.parameters);
 	
 	%Merge back to g set
-	[u_g]= Mergesolutionfromftog( u_f, m.Gmn, m.ys, m.nodesets ); 
-
-end %end function
+	u_g= Mergesolutionfromftog( u_f, femmodel.Gmn, femmodel.ys, femmodel.nodesets ); 
Index: /issm/trunk/src/m/solvers/solver_thermal_nonlinear.m
===================================================================
--- /issm/trunk/src/m/solvers/solver_thermal_nonlinear.m	(revision 4114)
+++ /issm/trunk/src/m/solvers/solver_thermal_nonlinear.m	(revision 4115)
@@ -1,68 +1,53 @@
-function [t_g ,loads, melting_offset]=thermal_core_nonlinear(m,analysis_type,sub_analysis_type)
-%THERMAL_CORE_NONLINEAR - core of thermal solution sequence.
-%   model is return together with temperature
+function [femmodel t_g ,melting_offset]=solver_thermal_nonlinear(femmodel)
+%SOLVER_THERMAL_NONLINEAR - core of thermal solution sequence.
+%   femmodel is returned together with temperature and melting_offset, in case loads have been modified
 %
 %   Usage:
-%      [t_g ,loads, melting_offset]=thermal_core_nonlinear(m,analysis_type,sub_analysis_type);
+%      [femmodel t_g melting_offset]=solver_thermal_nonlinear(femmodel)
 
 	count=1;
 	converged=0;
 
-%   we are going to return the loads, make them a variable of this routine
-	loads=m.loads;
+	%stiffness and load generation only:
+	femmodel.parameters.Kflag=1; femmodel.parameters.Pflag=1;
 
-	%stiffness and load generation only:
-	m.parameters.Kflag=1; m.parameters.Pflag=1;
-
-	displaystring(m.parameters.Verbose,'\n%s',['   starting direct shooting method']);
-	
+	displaystring(femmodel.parameters.Verbose,'\n%s',['   starting direct shooting method']);
 	while(~converged),
 
 		%Update inputs in datasets
 		if count==1 reset_penalties=0; else reset_penalties=1; end
-		[m.elements m.loads]=UpdateInputsFromConstant(m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters,reset_penalties,ResetPenaltiesEnum);
+		InputUpdateFromConstant(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,reset_penalties,ResetPenaltiesEnum);
 
 		%system matrices 
-		if ~m.parameters.Lowmem 
-			if count==1
-				displaystring(m.parameters.Verbose,'%s',['   system matrices']);
-				[K_gg_nopenalty, p_g_nopenalty]=SystemMatrices(m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters,analysis_type,sub_analysis_type);
-			end
-			displaystring(m.parameters.Verbose,'%s',['   penalty system matrices']);
-			[K_gg , p_g, melting_offset]=PenaltySystemMatrices(K_gg_nopenalty,p_g_nopenalty,m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters,analysis_type,sub_analysis_type);
-		else
-				m.parameters.Verbose
-			displaystring(m.parameters.Verbose,'%s',['   system matrices']);
-			[K_gg , p_g]=SystemMatrices(m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters,analysis_type,sub_analysis_type);
-			displaystring(m.parameters.Verbose,'%s',['   penalty system matrices']);
-			[K_gg , p_g, melting_offset]=PenaltySystemMatrices(K_gg,p_g,m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters,analysis_type,sub_analysis_type);
+		if count==1
+			displaystring(femmodel.parameters.Verbose,'%s',['   system matrices']);
+			[K_gg_nopenalty, p_g_nopenalty]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,analysis_type,sub_analysis_type);
 		end
+		displaystring(femmodel.parameters.Verbose,'%s',['   penalty system matrices']);
+		[K_gg , p_g, melting_offset]=PenaltySystemMatrices(K_gg_nopenalty,p_g_nopenalty,femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,analysis_type,sub_analysis_type);
 
 		%Reduce tangent matrix from g size to f size
-		[K_ff, K_fs] = Reducematrixfromgtof( K_gg, m.Gmn, m.nodesets); 
+		[K_ff, K_fs] = Reducematrixfromgtof( K_gg, femmodel.Gmn, femmodel.nodesets); 
 
 		%Reduce load from g size to f size
-		[p_f] = Reduceloadfromgtof( p_g, m.Gmn, K_fs, m.ys, m.nodesets);
+		p_f = Reduceloadfromgtof( p_g, femmodel.Gmn, K_fs, femmodel.ys, femmodel.nodesets);
 
-		%Solve	
-		displaystring(m.parameters.Verbose,'%s%g','   condition number of stiffness matrix: ',condest(K_ff));
-		[t_f]=Solver(K_ff,p_f,[],m.parameters);
+		displaystring(femmodel.parameters.Verbose,'%s%g','   condition number of stiffness matrix: ',condest(K_ff));
+		t_f=Solver(K_ff,p_f,[],femmodel.parameters);
 
-		%Merge back to g set
-		displaystring(m.parameters.Verbose,'%s',['   merging solution back to g set']);
-		[t_g]= Mergesolutionfromftog( t_f, m.Gmn, m.ys, m.nodesets ); 
+		displaystring(femmodel.parameters.Verbose,'%s',['   merging solution back to g set']);
+		t_g= Mergesolutionfromftog( t_f, femmodel.Gmn, femmodel.ys, femmodel.nodesets ); 
 
-		%Update inputs in datasets
-		displaystring(m.parameters.Verbose,'%s',['   update inputs']);
-		[m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters] = UpdateInputsFromVector(m.elements,m.nodes,m.vertices,loads,m.materials,m.parameters,t_g,TemperatureEnum, VertexEnum);
+		displaystring(femmodel.parameters.Verbose,'%s',['   update inputs']);
+		femmodel.elements=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,t_g);
 	
 		%penalty constraints
-		displaystring(m.parameters.Verbose,'%s',['   penalty constraints']);
-		[loads,constraints_converged,num_unstable_constraints] =PenaltyConstraints(m.elements,m.nodes,m.vertices,loads, m.materials,m.parameters,analysis_type,sub_analysis_type);
+		displaystring(femmodel.parameters.Verbose,'%s',['   penalty constraints']);
+		[femmodel.loads,constraints_converged,num_unstable_constraints] =PenaltyConstraints(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters,analysis_type,sub_analysis_type);
 	
 		if ~converged,
-			displaystring(m.parameters.Verbose,'%s%i','   #unstable constraints ',num_unstable_constraints);
+			displaystring(femmodel.parameters.Verbose,'%s%i','   #unstable constraints ',num_unstable_constraints);
 			
-			if num_unstable_constraints<=m.parameters.MinThermalConstraints,
+			if num_unstable_constraints<=femmodel.parameters.MinThermalConstraints,
 				converged=1;
 			end