Index: /issm/trunk/src/m/solutions/control_core.m
===================================================================
--- /issm/trunk/src/m/solutions/control_core.m	(revision 4133)
+++ /issm/trunk/src/m/solutions/control_core.m	(revision 4134)
@@ -25,91 +25,58 @@
 
 %initialize control parameters
-param_g=model.parameters.param_g;
+param_g=femmodel.parameters.param_g;
 
 %set optimization options.
-options=ControlOptions(model.parameters);
+options=ControlOptions(femmodel.parameters);
 
-for n=1:model.parameters.NSteps,
+for n=1:nsteps,
 
 	%set options
-	options=optimset(options,'MaxFunEvals',model.parameters.MaxIter(n));
+	options=optimset(options,'MaxFunEvals',femmodel.parameters.MaxIter(n));
 
-	disp(sprintf('\n%s%s%s%s\n',['   control method step ' num2str(n) '/' num2str(model.parameters.NSteps)]));
+	disp(sprintf('\n%s%s%s%s\n',['   control method step ' num2str(n) '/' num2str(femmodel.parameters.NSteps)]));
+	[femmodel.elements,femmodel.loads]=InputUpdateFromConstant(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,fit(n),FitEnum);
 
 	%In case we are running a steady state control method, compute new temperature field using new parameter distribution: 
-	if model.parameters.ControlSteady;
-		steadystate_results=steadystate_core(models); t_g=steadystate_results.t_g; 
-		inputs=add(inputs,'temperature',t_g,'doublevec',1,model.parameters.NumberOfNodes);
+	if control_steady;
+		femmodel=steadystate_core(femmodel);
 	end
 
-	%update inputs with new fit
-	inputs=add(inputs,'fit',model.parameters.Fit(n),'double');
-	inputs=add(inputs,model.parameters.ControlType,param_g,'doublevec',1,model.parameters.NumberOfNodes);
-
 	%Update inputs in datasets
-	[model.elements,model.nodes,model.vertices,model.loads,model.materials,model.parameters]=UpdateFromInputs(model.elements,model.nodes,model.vertices,model.loads,model.materials,model.parameters);
+	[femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=UpdateFromInputs(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
 
 	displaystring(verbose,'\n%s',['      computing gradJ...']);
-	results_grad=gradjcompute_core(models);
-	u_g=results_grad.u_g; c(n).grad_g=results_grad.grad_g;
-	if dim==3,
-		if isstokes,
-			inputs=add(inputs,'velocity',u_g,'doublevec',4,model.parameters.NumberOfNodes);
-		else
-			if model.parameters.ControlSteady;
-				inputs=add(inputs,'velocity',u_g,'doublevec',3,model.parameters.NumberOfNodes);
-			else
-				inputs=add(inputs,'velocity',u_g,'doublevec',2,model.parameters.NumberOfNodes);
-			end
-		end
-	else
-		inputs=add(inputs,'velocity',u_g,'doublevec',2,model.parameters.NumberOfNodes);
-	end
+	femmodel=gradient_core(femmodel);
 
-	if n>=2 & search_scalar==0,
-		displaystring(verbose,'\n%s',['      normalizing directions...']);
-		c(n).grad_g=Orth(c(n).grad_g,c(n-1).grad_g);
-	else
-		c(n).grad_g=Orth(c(n).grad_g,{});
-	end
-
-	%visualize direction.
-	if model.parameters.Plot
-		plot_direction;
+	%Return gradient if asked
+	if cm_gradient,
+		InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel.parameters,GradientEnum);
+		%%%%%%%%%%%%%%%%%%%%%% PROBLEM MUST GO TO THE END (SEE PARALLEL)
+		error
+		break;
 	end
 
 	displaystring(verbose,'\n%s',['      optimizing along gradient direction...']);
-	[search_scalar c(n).J]=ControlOptimization('objectivefunctionC',0,1,options,models,param_g,c(n).grad_g,n,model.parameters);
+	[search_scalar c(n).J]=ControlOptimization('objectivefunctionC',0,1,options,femmodel,param_g,c(n).grad_g,n,femmodel.parameters);
 
 	displaystring(verbose,'\n%s',['      updating parameter using optimized search scalar...']);
-	param_g=param_g+search_scalar*model.parameters.Optscal(n)*c(n).grad_g;
+	param_g=param_g+search_scalar*femmodel.parameters.Optscal(n)*c(n).grad_g;
 
 	displaystring(verbose,'\n%s',['      constraining the new distribution...']);
-	param_g=ControlConstrain(param_g,model.parameters);
+	param_g=ControlConstrain(param_g,femmodel.parameters);
 	
 	disp(['      value of misfit J after optimization #' num2str(n) ':' num2str(c(n).J)]);
 
 	%Has convergence been reached?
-	convergence=0;
-	if ~isnan(eps_cm),
-		i=n-2;
-		%go through the previous misfits(starting from n-2)
-		while (i>=1),
-			if (fit(i)==fit(n)),
-				%convergence test only if we have the same misfits
-				if ((c(i).J-c(n).J)/c(n).J <= eps_cm),
-					%convergence if convergence criteria fullfilled
-					convergence=1;
-					displaystring(verbose,'\n%s%g%s%g\n','      Convergence criterion: dJ/J = ',(c(i).J-c(n).J)/c(n).J,'<',eps_cm);
-				else
-					displaystring(verbose,'\n%s%g%s%g\n','      Convergence criterion: dJ/J = ',(c(i).J-c(n).J)/c(n).J,'>',eps_cm);
-				end
-				break;
-			end
-			i=i-1;                                                                                                                                         
-		end
+	converged=controlconvergence(J,fit,eps_cm,n);
+	if converged,
+		break;
 	end
-	%stop if convergence has been reached                                                                                                               
-	if (convergence), break; end
+
+	%Temporary saving every five iterations
+	if mod(n+1,5)==0;
+		displaystring(verbose,'\n%s',['      saving temporary results...']);
+		controlrestart(femmmodel,J);
+	end
 
 end
@@ -117,32 +84,14 @@
 %generate output
 displaystring(verbose,'\n%s',['      preparing final velocity solution...']);
-
-%compute final velocity from diagnostic_core (horiz+vertical)
-if model.parameters.ControlSteady;
-	inputs=add(inputs,model.parameters.ControlType,param_g,'doublevec',1,model.parameters.NumberOfNodes);
-	steadystate_results=steadystate_core(models); t_g=steadystate_results.t_g; 
-	u_g=steadystate_results.u_g;
-	t_g=steadystate_results.t_g;
-	m_g=steadystate_results.m_g;
+if control_steady,
+	femmodel=steadystate_core(femmodel);
 else
-	inputs=add(inputs,model.parameters.ControlType,param_g,'doublevec',1,model.parameters.NumberOfNodes);
-	results_diag=diagnostic_core(models);
-	u_g=results_diag.u_g;
+	femmodel=diagnostic_core(femmodel);
 end
 
-%Recover misfit at each iteration of the control method 
-J=zeros(length(c),1);
-for i=1:length(c),
-	J(i)=c(i).J;
-end
+%Some results not computed by diagnostic or steadystate
+InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type);
+femmodel.results.JEnum=J;
+femmodel.results.ControlTypeEnum=EnumAsString(control_type);
 
-%build results
-results.time=0;
-results.step=1;
-results.J=J;
-results.param_g=param_g;
-results.u_g=u_g;
-if model.parameters.ControlSteady,
-	results.t_g=t_g;
-	results.m_g=m_g;
-end
+end %end function
Index: /issm/trunk/src/m/solutions/controlconvergence.m
===================================================================
--- /issm/trunk/src/m/solutions/controlconvergence.m	(revision 4134)
+++ /issm/trunk/src/m/solutions/controlconvergence.m	(revision 4134)
@@ -0,0 +1,27 @@
+function converged=controlconvergence(J,fit,eps_cm,n)
+%CONTROLCONVERGENCE - determine the convergence of control_core solution
+%
+%   Usage:
+%       converged=controlconvergence(J,fit,eps_cm,n);
+
+	convergence=0;
+	if ~isnan(eps_cm),
+		i=n-2;
+		%go through the previous misfits(starting from n-2)
+		while (i>=1),
+			if (fit(i)==fit(n)),
+				%convergence test only if we have the same misfits
+				if ((c(i).J-c(n).J)/c(n).J <= eps_cm),
+					%convergence if convergence criteria fullfilled
+					convergence=1;
+					displaystring(verbose,'\n%s%g%s%g\n','      Convergence criterion: dJ/J = ',(c(i).J-c(n).J)/c(n).J,'<',eps_cm);
+				else
+					displaystring(verbose,'\n%s%g%s%g\n','      Convergence criterion: dJ/J = ',(c(i).J-c(n).J)/c(n).J,'>',eps_cm);
+				end
+				break;
+			end
+			i=i-1;                                                                                                                                         
+		end
+	end
+
+end %end function
Index: /issm/trunk/src/m/solutions/controlrestart.m
===================================================================
--- /issm/trunk/src/m/solutions/controlrestart.m	(revision 4134)
+++ /issm/trunk/src/m/solutions/controlrestart.m	(revision 4134)
@@ -0,0 +1,19 @@
+function femmodel=controlrestart(femmodel,J)
+%CONTROLRESTART - save as much as possible to be able to restart the control_core solution
+%
+%   Usage:
+%       femmodel=controlrestart(femmodel,J);
+
+	%recover parameters common to all solutions
+	nsteps=femmodel.parameters.NSteps;
+	control_type=femmodel.parameters.ControlType;
+	
+	%Save J and the parameter
+	InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel.parameters,control_type);
+	femmodel.results.JEnum=J;
+	femmodel.results.ControlTypeEnum=EnumAsString(control_type);
+
+	%Write to disk
+	OutputResults(femmodel.elements, femmodel.loads, femmodel.nodes, femmodel.vertices, femmodel.materials, femmodel.parameters);
+
+end %end function