Index: /issm/trunk/src/m/solutions/cielo/ControlInitialization.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/ControlInitialization.m	(revision 1652)
+++ /issm/trunk/src/m/solutions/cielo/ControlInitialization.m	(revision 1652)
@@ -0,0 +1,76 @@
+function [inputs model]=ControlInitialization(models,inputs);
+
+%recover models
+m_dh=models.dh;
+m_dv=models.dv;
+m_ds=models.ds;
+m_dhu=models.dhu;
+m_sl=models.sl;
+
+%recover parameters common to all solutions
+debug=m_dhu.parameters.debug;
+dim=m_dhu.parameters.dim;
+ishutter=m_dhu.parameters.ishutter;
+ismacayealpattyn=m_dh.parameters.ismacayealpattyn;
+isstokes=m_ds.parameters.isstokes;
+
+%If Pattyn or MacAyeal, just return, no initialization needed
+if ~isstokes
+	model=m_dh;
+	return
+end
+
+%1: Get slopes once for all
+
+%compute slopes
+displaystring(debug,'\n%s',['computing bed slope (x and y derivatives)...']);
+slopex=diagnostic_core_linear(m_sl,inputs,'slope_compute','bedx');
+slopey=diagnostic_core_linear(m_sl,inputs,'slope_compute','bedy');
+slopex=FieldExtrude(m_sl.elements,m_sl.nodes,m_sl.loads,m_sl.materials,slopex,'slopex',0);
+slopey=FieldExtrude(m_sl.elements,m_sl.nodes,m_sl.loads,m_sl.materials,slopey,'slopey',0);
+
+%Add to inputs
+inputs=add(inputs,'bedslopex',slopex,'doublevec',m_sl.parameters.numberofdofspernode,m_sl.parameters.numberofnodes);
+inputs=add(inputs,'bedslopey',slopey,'doublevec',m_sl.parameters.numberofdofspernode,m_sl.parameters.numberofnodes);
+
+%2: update spcs using MacAyeal or Pattyn
+
+%horizontal velocities
+displaystring(debug,'\n%s',['computing horizontal velocities...']);
+u_g=diagnostic_core_nonlinear(m_dh,inputs,'diagnostic','horiz');
+displaystring(debug,'\n%s',['extruding horizontal velocities...']);
+u_g_horiz=FieldExtrude(m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials,u_g,'velocity',1);
+
+%vertical velocities
+displaystring(debug,'\n%s',['computing vertical velocities...']);
+inputs=add(inputs,'velocity',u_g_horiz,'doublevec',m_dh.parameters.numberofdofspernode,m_dh.parameters.numberofnodes);
+u_g_vert=diagnostic_core_linear(m_dv,inputs,'diagnostic','vert');
+
+%create 3d u_g
+displaystring(debug,'\n%s',['combining horizontal and vertical velocities...']);
+u_g=zeros(m_dh.parameters.numberofnodes*3,1);
+u_g(dofsetgen([1,2],3,m_dh.parameters.numberofnodes*3))=u_g_horiz;
+u_g(dofsetgen([3],3,m_dh.parameters.numberofnodes*3))=u_g_vert;
+
+% get pressure (reconditionned) and create 4d u_g
+displaystring(debug,'\n%s',['computing pressure according to Pattyn...']);
+p_g=ComputePressure(m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials,m_dh.parameters,inputs);
+p_g=p_g/m_ds.parameters.stokesreconditioning;
+u_g_stokes=zeros(m_ds.nodesets.gsize,1);
+u_g_stokes(dofsetgen([1,2,3],4,m_ds.nodesets.gsize))=u_g;
+u_g_stokes(dofsetgen([4],4,m_ds.nodesets.gsize))=p_g;
+inputs=add(inputs,'velocity',u_g_stokes,'doublevec',4,m_ds.parameters.numberofnodes);
+
+%update BC (spcs)
+displaystring(debug,'\n%s',['update boundary conditions for stokes using velocities previously computed...']);
+m_ds.y_g=zeros(m_ds.nodesets.gsize,1);
+m_ds.y_g(dofsetgen([1,2,3],4,m_ds.nodesets.gsize))=u_g;
+[m_ds.ys m_ds.ys0]=Reducevectorgtos(m_ds.y_g,m_ds.nodesets);
+
+%Compute Stokes velocities once to have a reasonably good ug in input
+displaystring(debug,'\n%s',['computing stokes velocities and pressure ...']);
+u_g=diagnostic_core_nonlinear(m_ds,inputs,'diagnostic','stokes');
+inputs=add(inputs,'velocity',u_g,'doublevec',4,m_ds.parameters.numberofnodes);
+
+%assign output
+model=m_ds;
Index: sm/trunk/src/m/solutions/cielo/ControlPrepareStokes.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/ControlPrepareStokes.m	(revision 1651)
+++ 	(revision )
@@ -1,67 +1,0 @@
-function [inputs m_ds]=ControlPrepareStokes(models,inputs);
-
-%recover models
-m_dh=models.dh;
-m_dv=models.dv;
-m_ds=models.ds;
-m_dhu=models.dhu;
-m_sl=models.sl;
-
-%recover parameters common to all solutions
-debug=m_dhu.parameters.debug;
-dim=m_dhu.parameters.dim;
-ishutter=m_dhu.parameters.ishutter;
-ismacayealpattyn=m_dh.parameters.ismacayealpattyn;
-isstokes=m_ds.parameters.isstokes;
-
-%1: Get slopes once for all
-
-%compute slopes
-displaystring(debug,'\n%s',['computing bed slope (x and y derivatives)...']);
-slopex=diagnostic_core_linear(m_sl,inputs,'slope_compute','bedx');
-slopey=diagnostic_core_linear(m_sl,inputs,'slope_compute','bedy');
-slopex=FieldExtrude(m_sl.elements,m_sl.nodes,m_sl.loads,m_sl.materials,slopex,'slopex',0);
-slopey=FieldExtrude(m_sl.elements,m_sl.nodes,m_sl.loads,m_sl.materials,slopey,'slopey',0);
-
-%Add to inputs
-inputs=add(inputs,'bedslopex',slopex,'doublevec',m_sl.parameters.numberofdofspernode,m_sl.parameters.numberofnodes);
-inputs=add(inputs,'bedslopey',slopey,'doublevec',m_sl.parameters.numberofdofspernode,m_sl.parameters.numberofnodes);
-
-%2: update spcs using MacAyeal or Pattyn
-
-%horizontal velocities
-displaystring(debug,'\n%s',['computing horizontal velocities...']);
-u_g=diagnostic_core_nonlinear(m_dh,inputs,'diagnostic','horiz');
-displaystring(debug,'\n%s',['extruding horizontal velocities...']);
-u_g_horiz=FieldExtrude(m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials,u_g,'velocity',1);
-
-%vertical velocities
-displaystring(debug,'\n%s',['computing vertical velocities...']);
-inputs=add(inputs,'velocity',u_g_horiz,'doublevec',m_dh.parameters.numberofdofspernode,m_dh.parameters.numberofnodes);
-u_g_vert=diagnostic_core_linear(m_dv,inputs,'diagnostic','vert');
-
-%create 3d u_g
-displaystring(debug,'\n%s',['combining horizontal and vertical velocities...']);
-u_g=zeros(m_dh.parameters.numberofnodes*3,1);
-u_g(dofsetgen([1,2],3,m_dh.parameters.numberofnodes*3))=u_g_horiz;
-u_g(dofsetgen([3],3,m_dh.parameters.numberofnodes*3))=u_g_vert;
-
-% get pressure (reconditionned) and create 4d u_g
-displaystring(debug,'\n%s',['computing pressure according to Pattyn...']);
-p_g=ComputePressure(m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials,m_dh.parameters,inputs);
-p_g=p_g/m_ds.parameters.stokesreconditioning;
-u_g_stokes=zeros(m_ds.nodesets.gsize,1);
-u_g_stokes(dofsetgen([1,2,3],4,m_ds.nodesets.gsize))=u_g;
-u_g_stokes(dofsetgen([4],4,m_ds.nodesets.gsize))=p_g;
-inputs=add(inputs,'velocity',u_g_stokes,'doublevec',4,m_ds.parameters.numberofnodes);
-
-%update BC (spcs)
-displaystring(debug,'\n%s',['update boundary conditions for stokes using velocities previously computed...']);
-m_ds.y_g=zeros(m_ds.nodesets.gsize,1);
-m_ds.y_g(dofsetgen([1,2,3],4,m_ds.nodesets.gsize))=u_g;
-[m_ds.ys m_ds.ys0]=Reducevectorgtos(m_ds.y_g,m_ds.nodesets);
-
-%Compute Stokes velocities once to have a reasonably good ug in input
-displaystring(debug,'\n%s',['computing stokes velocities and pressure ...']);
-u_g=diagnostic_core_nonlinear(m_ds,inputs,'diagnostic','stokes');
-inputs=add(inputs,'velocity',u_g,'doublevec',4,m_ds.parameters.numberofnodes);
Index: /issm/trunk/src/m/solutions/cielo/control_core.m
===================================================================
--- /issm/trunk/src/m/solutions/cielo/control_core.m	(revision 1651)
+++ /issm/trunk/src/m/solutions/cielo/control_core.m	(revision 1652)
@@ -6,54 +6,33 @@
 %
 
-%recover models
-m_dh=models.dh;
-m_dv=models.dv;
-m_ds=models.ds;
-m_sl=models.sl;
+%Preprocess models
+[inputs model]=ControlInitialization(models,inputs);
 
 %recover parameters common to all solutions
-debug=m_dh.parameters.debug;
-dim=m_dh.parameters.dim;
-ishutter=m_dh.parameters.ishutter;
-ismacayealpattyn=m_dh.parameters.ismacayealpattyn;
-isstokes=m_dh.parameters.isstokes;
+debug=model.parameters.debug;
 
 %initialize control parameters
-param_g=models.dh.parameters.param_g;
+param_g=model.parameters.param_g;
 
 %set optimization options.
-options=ControlOptions(m_dh.parameters);
+options=ControlOptions(model.parameters);
 
-%Take car of Stokes : compute slope and get spc once for all
-if isstokes,
-	[inputs m_ds]=ControlPrepareStokes(models,inputs);
-end
-
-for n=1:m_dh.parameters.nsteps,
+for n=1:model.parameters.nsteps,
 
 	%set options
-	options=optimset(options,'MaxFunEvals',m_dh.parameters.maxiter(n));
+	options=optimset(options,'MaxFunEvals',model.parameters.maxiter(n));
 
-	disp(sprintf('\n%s%s%s%s\n',['   control method step ' num2str(n) '/' num2str(m_dh.parameters.nsteps)]));
+	disp(sprintf('\n%s%s%s%s\n',['   control method step ' num2str(n) '/' num2str(model.parameters.nsteps)]));
 
 	%update inputs with new fit
-	inputs=add(inputs,'fit',m_dh.parameters.fit(n),'double');
-	inputs=add(inputs,m_dh.parameters.control_type,param_g,'doublevec',1,m_dh.parameters.numberofnodes);
+	inputs=add(inputs,'fit',model.parameters.fit(n),'double');
+	inputs=add(inputs,model.parameters.control_type,param_g,'doublevec',1,model.parameters.numberofnodes);
 
 	%Update inputs in datasets
-	if isstokes,
-		[m_ds.elements,m_ds.nodes,m_ds.loads,m_ds.materials]=UpdateFromInputs(m_ds.elements,m_ds.nodes,m_ds.loads,m_ds.materials,inputs);
-	else
-		[m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials]=UpdateFromInputs(m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials,inputs);
-	end
+	[model.elements,model.nodes,model.loads,model.materials]=UpdateFromInputs(model.elements,model.nodes,model.loads,model.materials,inputs);
 
 	displaystring(debug,'\n%s',['      computing gradJ...']);
-	if isstokes,
-		[u_g c(n).grad_g]=GradJCompute(m_ds,inputs,DiagnosticAnalysisEnum(),StokesAnalysisEnum());
-		inputs=add(inputs,'velocity',u_g,'doublevec',4,m_ds.parameters.numberofnodes);
-	else
-		[u_g c(n).grad_g]=GradJCompute(m_dh,inputs,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
-		inputs=add(inputs,'velocity',u_g,'doublevec',2,m_dh.parameters.numberofnodes);
-	end
+	[u_g c(n).grad_g]=GradJCompute(model,inputs,model.parameters.analysis_type,model.parameters.sub_analysis_type);
+	inputs=add(inputs,'velocity',u_g,'doublevec',2,model.parameters.numberofnodes);
 
 	displaystring(debug,'\n%s',['      normalizing directions...']);
@@ -65,23 +44,19 @@
 
 	%visualize direction.
-	if m_dh.parameters.plot
+	if model.parameters.plot
 		plot_direction;
 	end
 
 	displaystring(debug,'\n%s',['      optimizing along gradient direction...']);
-	if isstokes,
-		[search_scalar c(n).J]=ControlOptimization('objectivefunctionC',0,1,options,m_ds,inputs,param_g,c(n).grad_g,n,DiagnosticAnalysisEnum(),StokesAnalysisEnum());
-	else
-		[search_scalar c(n).J]=ControlOptimization('objectivefunctionC',0,1,options,m_dh,inputs,param_g,c(n).grad_g,n,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
-	end
+	[search_scalar c(n).J]=ControlOptimization('objectivefunctionC',0,1,options,model,inputs,param_g,c(n).grad_g,n,model.parameters.analysis_type,model.parameters.sub_analysis_type);
 
 	displaystring(debug,'\n%s',['      updating parameter using optimized search scalar...']);
-	param_g=param_g+search_scalar*m_dh.parameters.optscal(n)*c(n).grad_g;
+	param_g=param_g+search_scalar*model.parameters.optscal(n)*c(n).grad_g;
 
 	displaystring(debug,'\n%s',['      constraining the new distribution...']);
-	param_g=ControlConstrain(param_g,m_dh.parameters);
+	param_g=ControlConstrain(param_g,model.parameters);
 
 	%visualize direction.
-	if m_dh.parameters.plot,
+	if model.parameters.plot,
 		plot_newdistribution;
 	end
@@ -94,5 +69,5 @@
 
 %compute final velocity from diagnostic_core (horiz+vertical)
-inputs=add(inputs,m_dh.parameters.control_type,param_g,'doublevec',1,m_dh.parameters.numberofnodes);
+inputs=add(inputs,model.parameters.control_type,param_g,'doublevec',1,model.parameters.numberofnodes);
 results_diag=diagnostic_core(models,inputs);