Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18127)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18128)
@@ -2851,4 +2851,5 @@
 
 	/*Get input (either in element or material)*/
+	if(control_enum==MaterialsRheologyBbarEnum) control_enum=MaterialsRheologyBEnum;
 	Input* input=inputs->GetInput(control_enum);
 	if(!input) _error_("Input " << EnumToStringx(control_enum) << " not found in element");
@@ -2865,7 +2866,18 @@
 
 	IssmDouble  values[NUMVERTICES];
-	int     vertexpidlist[NUMVERTICES];
-	Input  *input     = NULL;
-	Input  *new_input = NULL;
+	int         vertexpidlist[NUMVERTICES],control_init;
+	Input      *input     = NULL;
+	Input      *new_input = NULL;
+
+	/*Specific case for depth averaged quantities*/
+	control_init=control_enum;
+	if(control_enum==MaterialsRheologyBbarEnum){
+		control_enum=MaterialsRheologyBEnum;
+		if(!IsOnBase()) return;
+	}
+	if(control_enum==DamageDbarEnum){
+		control_enum=DamageDEnum;
+		if(!IsOnBase()) return;
+	}
 
 	/*Get out if this is not an element input*/
@@ -2876,5 +2888,5 @@
 
 	/*Get values on vertices*/
-	for (int i=0;i<NUMVERTICES;i++){
+	for(int i=0;i<NUMVERTICES;i++){
 		values[i]=vector[vertexpidlist[i]];
 	}
@@ -2887,4 +2899,11 @@
 
 	((ControlInput*)input)->SetInput(new_input);
+
+	if(control_init==MaterialsRheologyBbarEnum){
+		this->InputExtrude(control_enum);
+	}
+	if(control_init==DamageDbarEnum){
+		this->InputExtrude(control_enum);
+	}
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18127)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18128)
@@ -983,5 +983,4 @@
 		this->inputs->AddInput(new IntInput(ApproximationEnum,reCast<int>(iomodel->Data(FlowequationElementEquationEnum)[index])));
 	}
-
 
 	/*Control Inputs*/
Index: /issm/trunk-jpl/src/c/cores/control_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/control_core.cpp	(revision 18127)
+++ /issm/trunk-jpl/src/c/cores/control_core.cpp	(revision 18128)
@@ -11,6 +11,10 @@
 
 /*Local prototypes*/
-bool controlconvergence(IssmDouble J, IssmDouble tol_cm);
-IssmDouble objectivefunction(IssmDouble search_scalar,void* optargs);
+IssmDouble FormFunction(IssmDouble* X,void* usr);
+IssmDouble FormFunctionGradient(IssmDouble** pG,IssmDouble* X,void* usr);
+typedef struct {
+	FemModel* femmodel;
+	int       nsize;
+} AppCtx;
 
 void control_core(FemModel* femmodel){
@@ -19,5 +23,5 @@
 
 	/*parameters: */
-	int        num_controls;
+	int        num_controls,nsize;
 	int        nsteps;
 	IssmDouble tol_cm;
@@ -27,9 +31,8 @@
 
 	int        *control_type   = NULL;
-	IssmDouble *maxiter        = NULL;
+	int*        maxiter        = NULL;
 	IssmDouble *cm_jump        = NULL;
 
 	/*intermediary: */
-	IssmDouble search_scalar = 1.;
 	OptPars    optpars;
 
@@ -61,35 +64,39 @@
 	if(isFS) solutioncore(femmodel);
 
-	/*Initialize cost function: */
-	J=xNew<IssmDouble>(nsteps);
+	/*Get initial guess*/
+	Vector<IssmDouble> *Xpetsc = NULL;
+	GetVectorFromControlInputsx(&Xpetsc,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
+	IssmDouble* X0 = Xpetsc->ToMPISerial();
+	Xpetsc->GetSize(&nsize);
+	delete Xpetsc;
 
 	/*Initialize some of the BrentSearch arguments: */
-	optpars.xmin=0; optpars.xmax=1;
-
-	/*Start looping: */
-	for(int n=0;n<nsteps;n++){
-
-		/*Display info*/
-		if(VerboseControl()) _printf0_("\n" << "   control method step " << n+1 << "/" << nsteps << "\n");
-
-
-		/*In steady state inversion, compute new temperature field now*/
-		if(solution_type==SteadystateSolutionEnum) solutioncore(femmodel);
-
-		if(VerboseControl()) _printf0_("   compute adjoint state:\n");
-		adjointcore(femmodel);
-		gradient_core(femmodel,n,search_scalar==0.);
-
-		if(VerboseControl()) _printf0_("   optimizing along gradient direction\n");
-		optpars.maxiter=reCast<int,IssmDouble>(maxiter[n]); optpars.cm_jump=cm_jump[n];
-		BrentSearch(&search_scalar,J+n,&optpars,&objectivefunction,(void*)femmodel);
-
-		if(VerboseControl()) _printf0_("   updating parameter using optimized search scalar\n"); //true means update save controls
-		InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar,true);
-
-		if(controlconvergence(J[n],tol_cm)) break;
-	}
+	optpars.xmin    = 0; 
+	optpars.xmax    = 1;
+	optpars.nsteps  = nsteps;
+	optpars.nsize   = nsize;
+	optpars.maxiter = maxiter;
+	optpars.cm_jump = cm_jump;
+
+	/*Initialize function argument*/
+	AppCtx usr;
+	usr.femmodel = femmodel;
+	usr.nsize    = nsize;
+
+	/*Call Brent optimization*/
+	BrentSearch(&J,optpars,X0,&FormFunction,&FormFunctionGradient,(void*)&usr);
 
 	if(VerboseControl()) _printf0_("   preparing final solution\n");
+	IssmDouble  *XL = NULL;
+	IssmDouble  *XU = NULL;
+	GetVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
+	GetVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
+	for(long i=0;i<nsize;i++){
+		if(X0[i]>XU[i]) X0[i]=XU[i];
+		if(X0[i]<XL[i]) X0[i]=XL[i];
+	}
+	xDelete<IssmDouble>(XU);
+	xDelete<IssmDouble>(XL);
+	SetControlInputsFromVectorx(femmodel,X0);
 	femmodel->parameters->SetParam(true,SaveResultsEnum);
 	solutioncore(femmodel);
@@ -111,22 +118,11 @@
 	/*Free ressources: */
 	xDelete<int>(control_type);
-	xDelete<IssmDouble>(maxiter);
+	xDelete<int>(maxiter);
 	xDelete<IssmDouble>(cm_jump);
 	xDelete<IssmDouble>(J);
+	xDelete<IssmDouble>(X0);
 }
-bool controlconvergence(IssmDouble J, IssmDouble tol_cm){
-
-	bool converged=false;
-
-	/*Has convergence been reached?*/
-	if (!xIsNan<IssmDouble>(tol_cm) && J<tol_cm){
-		converged=true;
-		if(VerboseConvergence()) _printf0_("      Convergence criterion reached: J = " << J << " < " << tol_cm);
-	}
-
-	return converged;
-}
-
-IssmDouble objectivefunction(IssmDouble search_scalar,void* optargs){
+
+IssmDouble FormFunction(IssmDouble* X,void* usrvoid){
 
 	/*output: */
@@ -134,8 +130,10 @@
 
 	/*parameters: */
-	int        solution_type,analysis_type;
-	bool       isFS       = false;
+	int        solution_type,analysis_type,num_responses;
+	bool       isFS           = false;
 	bool       conserve_loads = true;
-	FemModel  *femmodel       = (FemModel*)optargs;
+	AppCtx*    usr = (AppCtx*)usrvoid;
+	FemModel  *femmodel  = usr->femmodel;
+	int        nsize     = usr->nsize;
 
 	/*Recover parameters: */
@@ -143,19 +141,39 @@
 	femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
 	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
-
-	/*set analysis type to compute velocity: */
+	femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
+
+	/*Constrain input vector*/
+	IssmDouble  *XL = NULL;
+	IssmDouble  *XU = NULL;
+	GetVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
+	GetVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
+	for(long i=0;i<nsize;i++){
+		if(X[i]>XU[i]) X[i]=XU[i];
+		if(X[i]<XL[i]) X[i]=XL[i];
+	}
+
+	/*Update control input*/
+	SetControlInputsFromVectorx(femmodel,X);
+
+	/*solve forward: */
 	switch(solution_type){
 		case SteadystateSolutionEnum:
+			femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
+			stressbalance_core(femmodel);	//We need a 3D velocity!! (vz is required for the next thermal run)
+			break;
 		case StressbalanceSolutionEnum:
 			femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
+			solutionsequence_nonlinear(femmodel,conserve_loads); 
 			break;
 		case BalancethicknessSolutionEnum:
 			femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum);
+			solutionsequence_linear(femmodel); 
 			break;
 		case BalancethicknessSoftSolutionEnum:
-			femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum);
+			/*NOTHING*/
 			break;
 		case Balancethickness2SolutionEnum:
 			femmodel->SetCurrentConfiguration(Balancethickness2AnalysisEnum);
+			solutionsequence_linear(femmodel); 
 			break;
 		default:
@@ -163,31 +181,123 @@
 	}
 
-	/*update parameter according to scalar: */ //false means: do not save control
-	InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar,false);
-
-	/*Run stressbalance with updated inputs: */
-	if (solution_type==SteadystateSolutionEnum){
-		stressbalance_core(femmodel);	//We need a 3D velocity!! (vz is required for the next thermal run)
-	}
-	else if (solution_type==StressbalanceSolutionEnum){
-		solutionsequence_nonlinear(femmodel,conserve_loads); 
-	}
-	else if (solution_type==BalancethicknessSolutionEnum){
-		solutionsequence_linear(femmodel); 
-	}
-	else if (solution_type==Balancethickness2SolutionEnum){
-		solutionsequence_linear(femmodel); 
-	}
-	else if (solution_type==BalancethicknessSoftSolutionEnum){
-		/*Don't do anything*/
-	}
-	else{
-		_error_("Solution " << EnumToStringx(solution_type) << " not implemented yet");
-	}
-
 	/*Compute misfit for this velocity field.*/
-	femmodel->CostFunctionx(&J,NULL,NULL);
+	IssmDouble* Jlist = NULL;
+	femmodel->CostFunctionx(&J,&Jlist,NULL);
+	//_printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<"  |  ");
+
+	/*Retrieve objective functions independently*/
+	//for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[i]);
+	//_printf0_("\n");
 
 	/*Free ressources:*/
+	xDelete<IssmDouble>(XU);
+	xDelete<IssmDouble>(XL);
+	xDelete<IssmDouble>(Jlist);
 	return J;
 }
+IssmDouble FormFunctionGradient(IssmDouble** pG,IssmDouble* X,void* usrvoid){
+
+	/*output: */
+	IssmDouble J;
+
+	/*parameters: */
+	void (*adjointcore)(FemModel*)=NULL;
+	int         solution_type,analysis_type,num_responses,num_controls,numvertices;
+	bool        isFS           = false;
+	bool        conserve_loads = true;
+	IssmDouble *scalar_list    = NULL;
+	IssmDouble *Jlist          = NULL;
+	IssmDouble *G              = NULL;
+	IssmDouble *norm_list      = NULL;
+	AppCtx     *usr            = (AppCtx*)usrvoid;
+	FemModel   *femmodel       = usr->femmodel;
+	int         nsize          = usr->nsize;
+
+	/*Recover parameters: */
+	femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
+	femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+	femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
+	femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
+	femmodel->parameters->FindParam(&scalar_list,NULL,NULL,InversionGradientScalingEnum);
+	femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);   _assert_(num_controls);
+	numvertices=femmodel->vertices->NumberOfVertices();
+
+	/*Constrain input vector*/
+	IssmDouble  *XL = NULL;
+	IssmDouble  *XU = NULL;
+	GetVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
+	GetVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
+	for(long i=0;i<nsize;i++){
+		if(X[i]>XU[i]) X[i]=XU[i];
+		if(X[i]<XL[i]) X[i]=XL[i];
+	}
+
+	/*Update control input*/
+	SetControlInputsFromVectorx(femmodel,X);
+
+	/*Compute Adjoint*/
+	AdjointCorePointerFromSolutionEnum(&adjointcore,solution_type);
+	adjointcore(femmodel);
+
+	/*Compute gradient*/
+	Gradjx(&G,&norm_list,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+
+	/*Compute scaling factor*/
+	IssmDouble scalar = scalar_list[0]/norm_list[0];
+	for(int i=1;i<num_controls;i++) scalar=min(scalar,scalar_list[i]/norm_list[i]);
+
+	/*Constrain Gradient*/
+	for(int i=0;i<num_controls;i++){
+		for(int j=0;j<numvertices;j++){
+			G[i*numvertices+j] = scalar*G[i*numvertices+j];
+		}
+	}
+
+	/*Needed for output results (FIXME: should be placed 6 lines below)*/
+	ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G);
+
+	for(long i=0;i<nsize;i++){
+		if(X[i]>=XU[i]) G[i]=0.;
+		if(X[i]<=XL[i]) G[i]=0.;
+	}
+
+	/*solve forward: (FIXME: not needed actually...)*/
+	switch(solution_type){
+		case SteadystateSolutionEnum:
+			femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
+			stressbalance_core(femmodel);	//We need a 3D velocity!! (vz is required for the next thermal run)
+			break;
+		case StressbalanceSolutionEnum:
+			femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
+			solutionsequence_nonlinear(femmodel,conserve_loads); 
+			break;
+		case BalancethicknessSolutionEnum:
+			femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum);
+			solutionsequence_linear(femmodel); 
+			break;
+		case BalancethicknessSoftSolutionEnum:
+			/*NOTHING*/
+			break;
+		case Balancethickness2SolutionEnum:
+			femmodel->SetCurrentConfiguration(Balancethickness2AnalysisEnum);
+			solutionsequence_linear(femmodel); 
+			break;
+		default:
+			_error_("Solution " << EnumToStringx(solution_type) << " not implemented yet");
+	}
+
+	/*Compute misfit for this velocity field.*/
+	femmodel->CostFunctionx(&J,&Jlist,NULL);
+	//_printf0_("f(x) = "<<setw(12)<<setprecision(7)<<*pf<<"  |  ");
+	//for(int i=0;i<num_responses;i++) _printf0_(" "<<setw(12)<<setprecision(7)<<Jlist[i]);
+	//_printf0_("\n");
+
+	/*Clean-up and return*/
+	xDelete<IssmDouble>(XU);
+	xDelete<IssmDouble>(XL);
+	xDelete<IssmDouble>(norm_list);
+	xDelete<IssmDouble>(scalar_list);
+	xDelete<IssmDouble>(Jlist);
+	*pG = G;
+	return J;
+}
Index: /issm/trunk-jpl/src/c/cores/controlm1qn3_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/controlm1qn3_core.cpp	(revision 18127)
+++ /issm/trunk-jpl/src/c/cores/controlm1qn3_core.cpp	(revision 18128)
@@ -106,5 +106,5 @@
 
 	/*Get solution*/
-	SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
+	SetControlInputsFromVectorx(femmodel,X);
 	ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G);
 	femmodel->OutputControlsx(&femmodel->results);
@@ -131,8 +131,7 @@
 	FemModel   *femmodel  = (FemModel*)dzs;
 
-	/*Recover responses*/
-	int         num_responses;
-	int        *responses = NULL;
-	femmodel->parameters->FindParam(&responses,&num_responses,InversionCostFunctionsEnum);
+	/*Recover number of cost functions responses*/
+	int num_responses;
+	femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
 
 	/*Constrain input vector*/
@@ -147,5 +146,5 @@
 
 	/*Update control input*/
-	SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
+	SetControlInputsFromVectorx(femmodel,X);
 
 	/*Recover some parameters*/
@@ -170,5 +169,4 @@
 	if(indic==0){
 		/*dry run, no gradient required*/
-		xDelete<int>(responses);
 		xDelete<IssmDouble>(XU);
 		xDelete<IssmDouble>(XL);
@@ -193,5 +191,4 @@
 
 	/*Clean-up and return*/
-	xDelete<int>(responses);
 	xDelete<IssmDouble>(XU);
 	xDelete<IssmDouble>(XL);
Index: /issm/trunk-jpl/src/c/cores/controltao_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/controltao_core.cpp	(revision 18127)
+++ /issm/trunk-jpl/src/c/cores/controltao_core.cpp	(revision 18128)
@@ -93,5 +93,5 @@
 	G=new Vector<IssmDouble>(0); VecFree(&G->pvector->vector);
 	TaoGetGradientVector(tao,&G->pvector->vector);
-	SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
+	SetControlInputsFromVectorx(femmodel,X);
 	ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G);
 	femmodel->OutputControlsx(&femmodel->results);
@@ -113,9 +113,9 @@
 	TaoFinalize();
 }
-int FormFunctionGradient(TaoSolver tao, Vec Xpetsc, IssmDouble *fcn,Vec G,void *userCtx){
+int FormFunctionGradient(TaoSolver tao, Vec Xpetsc, IssmDouble *fcn,Vec G,void *uservoid){
 
 	/*Retreive arguments*/
 	int                  solution_type;
-	AppCtx              *user            = (AppCtx *)userCtx;
+	AppCtx              *user            = (AppCtx *)uservoid;
 	FemModel            *femmodel        = user->femmodel;
 	Vector<IssmDouble>  *gradient        = NULL;
@@ -127,5 +127,5 @@
 	/*Set new variable*/
 	//VecView(X,PETSC_VIEWER_STDOUT_WORLD);
-	SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
+	SetControlInputsFromVectorx(femmodel,X);
 	delete X;
 
Index: /issm/trunk-jpl/src/c/cores/controlvalidation_core.cpp
===================================================================
--- /issm/trunk-jpl/src/c/cores/controlvalidation_core.cpp	(revision 18127)
+++ /issm/trunk-jpl/src/c/cores/controlvalidation_core.cpp	(revision 18128)
@@ -73,5 +73,5 @@
 
 		/*Calculate j(k+alpha delta k) */
-		SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X);
+		SetControlInputsFromVectorx(femmodel,X);
 		solutioncore(femmodel);
 		femmodel->CostFunctionx(&j,NULL,NULL);
Index: /issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp	(revision 18127)
+++ /issm/trunk-jpl/src/c/modules/Gradjx/Gradjx.cpp	(revision 18128)
@@ -50,4 +50,6 @@
 		delete gradient_list[i];
 	}
+	//gradient->Echo();
+	//_error_("S");
 
 	/*Check that gradient is clean*/
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp	(revision 18127)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Control/CreateParametersControl.cpp	(revision 18128)
@@ -16,8 +16,8 @@
 	int         num_cm_responses;
 	int        *control_type     = NULL;
+	int        *maxiter          = NULL;
 	int        *cm_responses     = NULL;
 	IssmDouble *cm_jump          = NULL;
 	IssmDouble *optscal          = NULL;
-	IssmDouble *maxiter          = NULL;
 
 	/*retrieve some parameters: */
@@ -56,5 +56,5 @@
 				parameters->AddObject(new DoubleMatParam(InversionGradientScalingEnum,optscal,nsteps,num_control_type));
 				parameters->AddObject(new DoubleVecParam(InversionStepThresholdEnum,cm_jump,nsteps));
-				parameters->AddObject(new DoubleVecParam(InversionMaxiterPerStepEnum,maxiter,nsteps));
+				parameters->AddObject(new IntVecParam(InversionMaxiterPerStepEnum,maxiter,nsteps));
 				break;
 			case 1:/*TAO*/
@@ -83,7 +83,7 @@
 		xDelete<int>(control_type);
 		xDelete<int>(cm_responses);
+		xDelete<int>(maxiter);
 		iomodel->DeleteData(cm_jump,InversionStepThresholdEnum);
 		iomodel->DeleteData(optscal,InversionGradientScalingEnum);
-		iomodel->DeleteData(maxiter,InversionMaxiterPerStepEnum);
 	}
 }
Index: /issm/trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.cpp	(revision 18127)
+++ /issm/trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.cpp	(revision 18128)
@@ -7,5 +7,5 @@
 #include "../../toolkits/toolkits.h"
 
-void SetControlInputsFromVectorx(Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,IssmDouble* vector){
+void SetControlInputsFromVectorx(FemModel* femmodel,IssmDouble* vector){
 
 	int  num_controls;
@@ -13,10 +13,10 @@
 
 	/*Retrieve some parameters*/
-	parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
-	parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
+	femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
+	femmodel->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
 
 	for(int i=0;i<num_controls;i++){
-		for(int j=0;j<elements->Size();j++){
-			Element* element=(Element*)elements->GetObjectByOffset(j);
+		for(int j=0;j<femmodel->elements->Size();j++){
+			Element* element=(Element*)femmodel->elements->GetObjectByOffset(j);
 			element->SetControlInputsFromVector(vector,control_type[i],i);
 		}
@@ -26,13 +26,8 @@
 }
 
-void SetControlInputsFromVectorx(Elements* elements,Nodes* nodes, Vertices* vertices, Loads* loads, Materials* materials, Parameters* parameters,Vector<IssmDouble>* vector){
+void SetControlInputsFromVectorx(FemModel* femmodel,Vector<IssmDouble>* vector){
 
-	IssmDouble* serial_vector=NULL;
-
-	serial_vector=vector->ToMPISerial();
-
-	SetControlInputsFromVectorx(elements,nodes, vertices, loads, materials, parameters,serial_vector);
-
-	/*Free ressources:*/
+	IssmDouble* serial_vector=vector->ToMPISerial();
+	SetControlInputsFromVectorx(femmodel,serial_vector);
 	xDelete<IssmDouble>(serial_vector);
 }
Index: /issm/trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.h
===================================================================
--- /issm/trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.h	(revision 18127)
+++ /issm/trunk-jpl/src/c/modules/SetControlInputsFromVectorx/SetControlInputsFromVectorx.h	(revision 18128)
@@ -8,6 +8,6 @@
 
 /* local prototypes: */
-void SetControlInputsFromVectorx(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials,  Parameters* parameters,Vector<IssmDouble>* vector);
-void SetControlInputsFromVectorx(Elements* elements,Nodes* nodes, Vertices* vertices,Loads* loads, Materials* materials,  Parameters* parameters,IssmDouble* vector);
+void SetControlInputsFromVectorx(FemModel* femmodel,Vector<IssmDouble>* vector);
+void SetControlInputsFromVectorx(FemModel* femmodel,IssmDouble* vector);
 
 #endif 
Index: /issm/trunk-jpl/src/c/shared/Numerics/BrentSearch.cpp
===================================================================
--- /issm/trunk-jpl/src/c/shared/Numerics/BrentSearch.cpp	(revision 18127)
+++ /issm/trunk-jpl/src/c/shared/Numerics/BrentSearch.cpp	(revision 18128)
@@ -20,5 +20,5 @@
 #include "./isnan.h"
 
-void BrentSearch(IssmDouble* psearch_scalar,IssmDouble* pJ,OptPars* optpars,IssmDouble (*f)(IssmDouble,void*),void* optargs){
+void BrentSearch(IssmDouble** pJ,OptPars optpars,IssmDouble* X0,IssmDouble (*f)(IssmDouble*,void*),IssmDouble (*g)(IssmDouble**,IssmDouble*,void*),void* usr){
 
 	/* This routine is optimizing a given function using Brent's method
@@ -26,198 +26,219 @@
 
 	/*Intermediary*/
+	int        iter;
 	IssmDouble si,gold,intervalgold,oldintervalgold;
-	IssmDouble parab_num,parab_den;
-	IssmDouble distance,cm_jump;
+	IssmDouble parab_num,parab_den,distance;
 	IssmDouble fxmax,fxmin,fxbest;
 	IssmDouble fx,fx1,fx2;
-	IssmDouble xmax,xmin,xbest;
-	IssmDouble x,x1,x2,xm;
+	IssmDouble x,x1,x2,xm,xbest;
 	IssmDouble tol1,tol2,seps;
 	IssmDouble tolerance = 1.e-4;
-	int        maxiter ,iter;
-	bool       loop= true,goldenflag;
 
 	/*Recover parameters:*/
-	xmin=optpars->xmin;
-	xmax=optpars->xmax;
-	maxiter=optpars->maxiter;
-	cm_jump=optpars->cm_jump;
-
-	/*initialize counter and get response at the boundaries*/
-	iter=0;
-	fxmin = (*f)(xmin,optargs);
-	if (xIsNan<IssmDouble>(fxmin)) _error_("Function evaluation returned NaN");
-	cout<<setprecision(5);
-	if(VerboseControl()) _printf0_("\n");
-	if(VerboseControl()) _printf0_("       Iteration         x           f(x)       Tolerance         Procedure\n");
-	if(VerboseControl()) _printf0_("\n");
-	if(VerboseControl()) _printf0_("           N/A    "<<setw(12)<<xmin<<"  "<<setw(12)<<fxmin<<"           N/A         boundary\n");
-	fxmax = (*f)(xmax,optargs);
-	if (xIsNan<IssmDouble>(fxmax)) _error_("Function evaluation returned NaN");
-	if(VerboseControl()) _printf0_("           N/A    "<<setw(12)<<xmax<<"  "<<setw(12)<<fxmax<<"           N/A         boundary\n");
-
-	/*test if jump option activated and xmin==0*/
-	if (!xIsNan<IssmDouble>(cm_jump) && (xmin==0) && (fxmax/fxmin)<cm_jump){
-		*psearch_scalar=xmax;
-		*pJ=fxmax;
-		return;
+	int         nsteps  = optpars.nsteps;
+	int         nsize   = optpars.nsize;
+	IssmDouble  xmin    = optpars.xmin;
+	IssmDouble  xmax    = optpars.xmax;
+	int        *maxiter = optpars.maxiter;
+	IssmDouble *cm_jump = optpars.cm_jump;
+
+	/*Initialize gradient and controls*/
+	IssmDouble* G = NULL;
+	IssmDouble* J = xNew<IssmDouble>(nsteps);
+	IssmDouble* X = xNew<IssmDouble>(nsize);
+
+	/*start iterations*/
+	for(int n=0;n<nsteps;n++){
+
+		/*Reset some variables*/
+		iter = 0;
+		xmin = 0.;
+		xmax = 1.;
+		bool loop = true;
+		cout<<setprecision(5);
+
+		/*Get current Gradient at xmin=0*/
+		if(VerboseControl()) _printf0_("\n" << "   step " << n+1 << "/" << nsteps << "\n");
+		fxmin = (*g)(&G,X0,usr); if(xIsNan<IssmDouble>(fxmin)) _error_("Function evaluation returned NaN");
+		if(VerboseControl()) _printf0_("\n");
+		if(VerboseControl()) _printf0_("       Iteration         x           f(x)       Tolerance         Procedure\n");
+		if(VerboseControl()) _printf0_("\n");
+		if(VerboseControl()) _printf0_("           N/A    "<<setw(12)<<xmin<<"  "<<setw(12)<<fxmin<<"           N/A         boundary\n");
+		
+		/*Get f(xmax)*/
+		for(int i=0;i<nsize;i++) X[i]=X0[i]+xmax*G[i];
+		fxmax = (*f)(X,usr); if (xIsNan<IssmDouble>(fxmax)) _error_("Function evaluation returned NaN");
+		if(VerboseControl()) _printf0_("           N/A    "<<setw(12)<<xmax<<"  "<<setw(12)<<fxmax<<"           N/A         boundary\n");
+
+		/*test if jump option activated and xmin==0*/
+		if(!xIsNan<IssmDouble>(cm_jump[n]) && (xmin==0) && (fxmax/fxmin)<cm_jump[n]){
+			for(int i=0;i<nsize;i++) X0[i]=X0[i]+xmax*G[i];
+			xDelete<IssmDouble>(G);
+			J[n]=fxmax;
+			continue;
+		}
+
+		/*initialize optimization variables*/
+		seps=sqrt(DBL_EPSILON);    //precision of a IssmDouble
+		distance=0.0;              //new_x=old_x + distance
+		gold=0.5*(3.0-sqrt(5.0));  //gold = 1 - golden ratio
+		intervalgold=0.0;          //distance used by Golden procedure
+
+		/*1: initialize the values of the 4 x needed (x1,x2,x,xbest)*/
+		x1=xmin+gold*(xmax-xmin);
+		x2=x1;
+		xbest=x1;
+		x=xbest;
+
+		/*2: call the function to be evaluated*/
+		for(int i=0;i<nsize;i++) X[i]=X0[i]+x*G[i];
+		fxbest = (*f)(X,usr); if(xIsNan<IssmDouble>(fxbest)) _error_("Function evaluation returned NaN");
+		iter++;
+
+		/*3: update the other variables*/
+		fx1=fxbest;
+		fx2=fxbest;
+		/*xm is always in the middle of a and b*/
+		xm=0.5*(xmin+xmax);                           
+		/*update tolerances*/
+		tol1=seps*sqrt(pow(xbest,2))+tolerance/3.0;
+		tol2=2.0*tol1;
+
+		/*4: print result*/
+		if(VerboseControl())
+		 _printf0_("         "<<setw(5)<<iter<<"    "<<setw(12)<<xbest<<"  "<<setw(12)<<fxbest<<"  "<<setw(12)<<pow(pow(xbest-xm,2),0.5)<<"         initial\n");
+		if (!xIsNan<IssmDouble>(cm_jump[n]) && (xmin==0) && ((fxbest/fxmin)<cm_jump[n])){
+			if(VerboseControl()) _printf0_("      optimization terminated: current x satisfies criteria 'cm_jump'=" << cm_jump[n] << "\n");
+			loop=false;
+		}
+
+		while(loop){
+
+			bool goldenflag=true;
+
+			// Is a parabolic fit possible ?
+			if (sqrt(pow(intervalgold,2))>tol1){
+
+				// Yes, so fit parabola
+				goldenflag=false;
+				parab_num=(xbest-x1)*(xbest-x1)*(fxbest-fx2)-(xbest-x2)*(xbest-x2)*(fxbest-fx1);;
+				parab_den=2.0*(xbest-x1)*(fxbest-fx2)-2.0*(xbest-x2)*(fxbest-fx1);
+
+				//reverse p if necessary
+				if(parab_den>0.0){ 
+					parab_num=-parab_num;
+				}
+				parab_den=sqrt(pow(parab_den,2));
+				oldintervalgold=intervalgold;
+				intervalgold=distance;
+
+				// Is the parabola acceptable (we use seps here because in some configuration parab_num==parab_den*(xmax-xbest)
+				// and the result is not repeatable anymore
+				if (( sqrt(pow(parab_num,2)) < sqrt(pow(0.5*parab_den*oldintervalgold,2))) &&
+							(parab_num>parab_den*(xmin-xbest)+seps) && 
+							(parab_num<parab_den*(xmax-xbest)-seps)){
+
+					// Yes, parabolic interpolation step
+					distance=parab_num/parab_den;
+					x=xbest+distance;
+
+					// f must not be evaluated too close to min_x or max_x
+					if (((x-xmin)<tol2) || ((xmax-x)<tol2)){
+						if ((xm-xbest)<0.0) si=-1;
+						else                si=1;
+						//compute new distance
+						distance=tol1*si;
+					}
+				}
+				else{
+					// Not acceptable, must do a golden section step
+					goldenflag=true;
+				}
+			}
+
+			//Golden procedure
+			if(goldenflag){
+				// compute the new distance d
+				if(xbest>=xm){
+					intervalgold=xmin-xbest;    
+				}
+				else{ 
+					intervalgold=xmax-xbest;  
+				}
+				distance=gold*intervalgold;
+			}
+
+			// The function must not be evaluated too close to xbest
+			if(distance<0) si=-1;
+			else           si=1;
+			if(sqrt(pow(distance,2))>tol1) x=xbest+si*sqrt(pow(distance,2));
+			else                           x=xbest+si*tol1;
+
+			//evaluate function on x
+			for(int i=0;i<nsize;i++) X[i]=X0[i]+x*G[i];
+			fx = (*f)(X,usr); if(xIsNan<IssmDouble>(fx)) _error_("Function evaluation returned NaN");
+			iter++;
+
+			// Update a, b, xm, x1, x2, tol1, tol2
+			if (fx<=fxbest){
+				if (x>=xbest) xmin=xbest;
+				else          xmax=xbest;
+				x1=x2;    fx1=fx2;
+				x2=xbest; fx2=fxbest;
+				xbest=x;  fxbest=fx;
+			}
+			else{ // fx > fxbest
+				if (x<xbest) xmin=x;
+				else         xmax=x;
+				if ((fx<=fx2) || (x2==xbest)){
+					x1=x2; fx1=fx2;
+					x2=x;  fx2=fx;
+				}
+				else if ( (fx <= fx1) || (x1 == xbest) || (x1 == x2) ){
+					x1=x;  fx1=fx;
+				}
+			}
+			xm = 0.5*(xmin+xmax);
+			tol1=seps*pow(pow(xbest,2),0.5)+tolerance/3.0;
+			tol2=2.0*tol1;
+			if(VerboseControl())
+			 _printf0_("         "<<setw(5)<<iter<<"    "<<setw(12)<<x<<"  "<<setw(12)<<fx<<"  "<<setw(12)<<pow(pow(xbest-xm,2),0.5)<<
+						 "         "<<(goldenflag?"golden":"parabolic")<<"\n");
+
+			/*Stop the optimization?*/
+			if (sqrt(pow(xbest-xm,2)) < (tol2-0.5*(xmax-xmin))){
+				if(VerboseControl()) _printf0_("      optimization terminated: current x satisfies criteria 'tolx'=" << tolerance << "\n");
+				loop=false;
+			}
+			else if (iter>=maxiter[n]){
+				if(VerboseControl()) _printf0_("      exiting: Maximum number of iterations has been exceeded  ('maxiter'=" << maxiter[n] << ")\n");
+				loop=false;
+			}
+			else if (!xIsNan<IssmDouble>(cm_jump[n]) && (xmin==0) && ((fxbest/fxmin)<cm_jump[n])){
+				if(VerboseControl()) _printf0_("      optimization terminated: current x satisfies criteria 'cm_jump'=" << cm_jump[n] << "\n");
+				loop=false;
+			}
+			else{
+				//continue
+				loop=true;
+			}
+		}//end while
+
+		//Now, check that the value on the boundaries are not better than current fxbest
+		if (fxbest>fxmin){
+			xbest=optpars.xmin; fxbest=fxmin;
+		}
+		if (fxbest>fxmax){
+			xbest=optpars.xmax; fxbest=fxmax;
+		}
+
+		/*Assign output pointers: */
+		for(int i=0;i<nsize;i++) X0[i]=X0[i]+xbest*G[i];
+		xDelete<IssmDouble>(G);
+		J[n]=fxbest;
 	}
-
-	/*initialize optimization variables*/
-	seps=sqrt(DBL_EPSILON);    //precision of a IssmDouble
-	distance=0.0;              //new_x=old_x + distance
-	gold=0.5*(3.0-sqrt(5.0));  //gold = 1 - golden ratio
-	intervalgold=0.0;          //distance used by Golden procedure
-
-	/*1: initialize the values of the 4 x needed (x1,x2,x,xbest)*/
-	x1=xmin+gold*(xmax-xmin);
-	x2=x1;
-	xbest=x1;
-	x=xbest;
-
-	/*2: call the function to be evaluated*/
-	fxbest = (*f)(x,optargs);
-	if(xIsNan<IssmDouble>(fxbest)) _error_("Function evaluation returned NaN");
-	iter=iter+1;
-
-	/*3: update the other variables*/
-	fx1=fxbest;
-	fx2=fxbest;
-	/*xm is always in the middle of a and b*/
-	xm=0.5*(xmin+xmax);                           
-	/*update tolerances*/
-	tol1=seps*sqrt(pow(xbest,2))+tolerance/3.0;
-	tol2=2.0*tol1;
-
-	/*4: print result*/
-	if(VerboseControl())
-	 _printf0_("         "<<setw(5)<<iter<<"    "<<setw(12)<<xbest<<"  "<<setw(12)<<fxbest<<"  "<<setw(12)<<pow(pow(xbest-xm,2),0.5)<<"         initial\n");
-	if (!xIsNan<IssmDouble>(cm_jump) && (xmin==0) && ((fxbest/fxmin)<cm_jump)){
-		if(VerboseControl()) _printf0_("      optimization terminated: current x satisfies criteria 'cm_jump'=" << cm_jump << "\n");
-		loop=false;
-	}
-
-	while(loop){
-
-		goldenflag=true;
-
-		// Is a parabolic fit possible ?
-		if (sqrt(pow(intervalgold,2))>tol1){
-
-			// Yes, so fit parabola
-			goldenflag=false;
-			parab_num=(xbest-x1)*(xbest-x1)*(fxbest-fx2)-(xbest-x2)*(xbest-x2)*(fxbest-fx1);;
-			parab_den=2.0*(xbest-x1)*(fxbest-fx2)-2.0*(xbest-x2)*(fxbest-fx1);
-
-			//reverse p if necessary
-			if(parab_den>0.0){ 
-				parab_num=-parab_num;
-			}
-			parab_den=sqrt(pow(parab_den,2));
-			oldintervalgold=intervalgold;
-			intervalgold=distance;
-
-			// Is the parabola acceptable (we use seps here because in some configuration parab_num==parab_den*(xmax-xbest)
-			// and the result is not repeatable anymore
-			if (( sqrt(pow(parab_num,2)) < sqrt(pow(0.5*parab_den*oldintervalgold,2))) &&
-						(parab_num>parab_den*(xmin-xbest)+seps) && 
-						(parab_num<parab_den*(xmax-xbest)-seps)){
-
-				// Yes, parabolic interpolation step
-				distance=parab_num/parab_den;
-				x=xbest+distance;
-
-				// f must not be evaluated too close to min_x or max_x
-				if (((x-xmin)<tol2) || ((xmax-x)<tol2)){
-					if ((xm-xbest)<0.0) si=-1;
-					else                si=1;
-					//compute new distance
-					distance=tol1*si;
-				}
-			}
-			else{
-				// Not acceptable, must do a golden section step
-				goldenflag=true;
-			}
-		}
-
-		//Golden procedure
-		if(goldenflag){
-			// compute the new distance d
-			if(xbest>=xm){
-				intervalgold=xmin-xbest;    
-			}
-			else{ 
-				intervalgold=xmax-xbest;  
-			}
-			distance=gold*intervalgold;
-		}
-
-		// The function must not be evaluated too close to xbest
-		if(distance<0) si=-1;
-		else           si=1;
-		if(sqrt(pow(distance,2))>tol1) x=xbest+si*sqrt(pow(distance,2));
-		else                           x=xbest+si*tol1;
-
-		//evaluate function on x
-		fx = (*f)(x,optargs);
-		if(xIsNan<IssmDouble>(fx)) _error_("Function evaluation returned NaN");
-		iter=iter+1;
-
-		// Update a, b, xm, x1, x2, tol1, tol2
-		if (fx<=fxbest){
-			if (x>=xbest) xmin=xbest;
-			else          xmax=xbest;
-			x1=x2;    fx1=fx2;
-			x2=xbest; fx2=fxbest;
-			xbest=x;  fxbest=fx;
-		}
-		else{ // fx > fxbest
-			if (x<xbest) xmin=x;
-			else         xmax=x;
-			if ((fx<=fx2) || (x2==xbest)){
-				x1=x2; fx1=fx2;
-				x2=x;  fx2=fx;
-			}
-			else if ( (fx <= fx1) || (x1 == xbest) || (x1 == x2) ){
-				x1=x;  fx1=fx;
-			}
-		}
-		xm = 0.5*(xmin+xmax);
-		tol1=seps*pow(pow(xbest,2),0.5)+tolerance/3.0;
-		tol2=2.0*tol1;
-		if(VerboseControl())
-		 _printf0_("         "<<setw(5)<<iter<<"    "<<setw(12)<<x<<"  "<<setw(12)<<fx<<"  "<<setw(12)<<pow(pow(xbest-xm,2),0.5)<<
-					 "         "<<(goldenflag?"golden":"parabolic")<<"\n");
-
-		/*Stop the optimization?*/
-		if (sqrt(pow(xbest-xm,2)) < (tol2-0.5*(xmax-xmin))){
-			if(VerboseControl()) _printf0_("      optimization terminated: current x satisfies criteria 'tolx'=" << tolerance << "\n");
-			loop=false;
-		}
-		else if (iter>=maxiter){
-			if(VerboseControl()) _printf0_("      exiting: Maximum number of iterations has been exceeded  ('maxiter'=" << maxiter << ")\n");
-			loop=false;
-		}
-		else if (!xIsNan<IssmDouble>(cm_jump) && (xmin==0) && ((fxbest/fxmin)<cm_jump)){
-			if(VerboseControl()) _printf0_("      optimization terminated: current x satisfies criteria 'cm_jump'=" << cm_jump << "\n");
-			loop=false;
-		}
-		else{
-			//continue
-			loop=true;
-		}
-	}//end while
-
-	//Now, check that the value on the boundaries are not better than current fxbest
-	if (fxbest>fxmin){
-		xbest=optpars->xmin; fxbest=fxmin;
-	}
-	if (fxbest>fxmax){
-		xbest=optpars->xmax; fxbest=fxmax;
-	}
-
-	/*Assign output pointers: */
-	*psearch_scalar=xbest;
-	*pJ=fxbest;
+	
+	/*return*/
+	xDelete<IssmDouble>(X);
+	*pJ=J;
 }
Index: /issm/trunk-jpl/src/c/shared/Numerics/OptPars.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Numerics/OptPars.h	(revision 18127)
+++ /issm/trunk-jpl/src/c/shared/Numerics/OptPars.h	(revision 18128)
@@ -10,8 +10,10 @@
 struct OptPars{
 
-	IssmDouble xmin;
-	IssmDouble xmax;
-	IssmDouble cm_jump;
-	int maxiter;
+	IssmDouble  xmin;
+	IssmDouble  xmax;
+	IssmDouble *cm_jump;
+	int* maxiter;
+	int  nsteps;
+	int  nsize;
 
 };
Index: /issm/trunk-jpl/src/c/shared/Numerics/numerics.h
===================================================================
--- /issm/trunk-jpl/src/c/shared/Numerics/numerics.h	(revision 18127)
+++ /issm/trunk-jpl/src/c/shared/Numerics/numerics.h	(revision 18128)
@@ -30,5 +30,5 @@
 int         min(int a,int b);
 int         max(int a,int b);
-void        BrentSearch(IssmDouble *psearch_scalar,IssmDouble*pJ,OptPars*optpars,IssmDouble (*f)(IssmDouble,void*),void* optargs);
+void        BrentSearch(IssmDouble** pJ,OptPars optpars,IssmDouble* X0,IssmDouble (*f)(IssmDouble*,void*),IssmDouble (*g)(IssmDouble**,IssmDouble*,void*),void* usr);
 void        cross(IssmDouble *result,IssmDouble*vector1,IssmDouble*vector2);
 void        XZvectorsToCoordinateSystem(IssmDouble *T,IssmDouble*xzvectors);
Index: /issm/trunk-jpl/src/m/classes/inversion.m
===================================================================
--- /issm/trunk-jpl/src/m/classes/inversion.m	(revision 18127)
+++ /issm/trunk-jpl/src/m/classes/inversion.m	(revision 18128)
@@ -25,53 +25,53 @@
 	end
 	methods
-         function createxml(obj,fid) % {{{
-            fprintf(fid, '<!-- inversion -->\n');            
-                    
-            % inversion parameters
-            fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="inversion parameters">','<section name="inversion" />');                    
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="iscontrol" type="',class(obj.iscontrol),'" default="',convert2str(obj.iscontrol),'">','     <section name="inversion" />','     <help> is inversion activated? </help>','  </parameter>');
-            
-            % incompleteadjoing drop-down (0 or 1)
-            fprintf(fid,'%s\n%s\n%s\n%s\n','  <parameter key ="incomplete_adjoint" type="alternative" optional="false">','     <section name="inversion" />','     <help> 1: linear viscosity, 0: non-linear viscosity </help>');
-            fprintf(fid,'%s\n','       <option value="0" type="string" default="true"> </option>');
-            fprintf(fid,'%s\n%s\n','       <option value="1" type="string" default="false"> </option>','</parameter>');
-            
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="control_parameters" type="',class(obj.control_parameters),'" default="',convert2str(obj.control_parameters),'">','     <section name="inversion" />','     <help> ex: {''FrictionCoefficient''}, or {''MaterialsRheologyBbar''} </help>','  </parameter>');
-                
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="nsteps" type="',class(obj.nsteps),'" default="',convert2str(obj.nsteps),'">','     <section name="inversion" />','     <help> number of optimization searches </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="cost_functions" type="',class(obj.cost_functions),'" default="',convert2str(obj.cost_functions),'">','     <section name="inversion" />','     <help> indicate the type of response for each optimization step  </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="cost_functions_coefficients" type="',class(obj.cost_functions_coefficients),'" default="',convert2str(obj.cost_functions_coefficients),'">','     <section name="inversion" />','     <help> cost_functions_coefficients applied to the misfit of each vertex and for each control_parameter </help>','  </parameter>');
-                
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="cost_function_threshold" type="',class(obj.cost_function_threshold),'" default="',convert2str(obj.cost_function_threshold),'">','     <section name="inversion" />','     <help> misfit convergence criterion. Default is 1%, NaN if not applied </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="maxiter_per_step" type="',class(obj.maxiter_per_step),'" default="',convert2str(obj.maxiter_per_step),'">','     <section name="inversion" />','     <help> maximum iterations during each optimization step  </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="gradient_scaling" type="',class(obj.gradient_scaling),'" default="',convert2str(obj.gradient_scaling),'">','     <section name="inversion" />','     <help> scaling factor on gradient direction during optimization, for each optimization step </help>','  </parameter>');
-               
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="step_threshold" type="',class(obj.step_threshold),'" default="',convert2str(obj.step_threshold),'">','     <section name="inversion" />','     <help> decrease threshold for misfit, default is 30% </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="min_parameters" type="',class(obj.min_parameters),'" default="',convert2str(obj.min_parameters),'">','     <section name="inversion" />','     <help> absolute minimum acceptable value of the inversed parameter on each vertex </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="max_parameters" type="',class(obj.max_parameters),'" default="',convert2str(obj.max_parameters),'">','     <section name="inversion" />','     <help> absolute maximum acceptable value of the inversed parameter on each vertex </help>','  </parameter>');
-               
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vx_obs" type="',class(obj.vx_obs),'" default="',convert2str(obj.vx_obs),'">','     <section name="inversion" />','     <help> observed velocity x component [m/yr] </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vy_obs" type="',class(obj.vy_obs),'" default="',convert2str(obj.vy_obs),'">','     <section name="inversion" />','     <help> observed velocity y component [m/yr]  </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vel_obs" type="',class(obj.vel_obs),'" default="',convert2str(obj.vel_obs),'">','     <section name="inversion" />','     <help> observed velocity magnitude [m/yr] </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="thickness_obs" type="',class(obj.thickness_obs),'" default="',convert2str(obj.thickness_obs),'">','     <section name="inversion" />','     <help> observed thickness [m]) </help>','  </parameter>');
-               
-            fprintf(fid,'%s\n%s\n','</frame>');    
-            
-            fprintf(fid,'%s\n%s\n%s\n','<frame key="2" label="Available cost functions">','<section name="inversion" />');                    
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceAbsVelMisfit" type="','string','" default="','101','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceRelVelMisfit" type="','string','" default="','102','">','     <section name="inversion" />','     <help>   </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceLogVelMisfit" type="','string','" default="','103','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
-                
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceLogVxVyMisfit" type="','string','" default="','104','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceAverageVelMisfit" type="','string','" default="','105','">','     <section name="inversion" />','     <help>   </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="ThicknessAbsMisfit" type="','string','" default="','106','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
-               
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="DragCoefficientAbsGradient" type="','string','" default="','107','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="RheologyBbarAbsGradient" type="','string','" default="','108','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
-            fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="ThicknessAbsGradient" type="','string','" default="','109','">','     <section name="inversion" />','     <help> </help>','  </parameter>');
-               
-            fprintf(fid,'%s\n%s\n','</frame>');    
-        
-        end % }}}       
+		function createxml(obj,fid) % {{{
+			fprintf(fid, '<!-- inversion -->\n');            
+
+			% inversion parameters
+			fprintf(fid,'%s\n%s\n%s\n','<frame key="1" label="inversion parameters">','<section name="inversion" />');                    
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="iscontrol" type="',class(obj.iscontrol),'" default="',convert2str(obj.iscontrol),'">','     <section name="inversion" />','     <help> is inversion activated? </help>','  </parameter>');
+
+			% incompleteadjoing drop-down (0 or 1)
+			fprintf(fid,'%s\n%s\n%s\n%s\n','  <parameter key ="incomplete_adjoint" type="alternative" optional="false">','     <section name="inversion" />','     <help> 1: linear viscosity, 0: non-linear viscosity </help>');
+			fprintf(fid,'%s\n','       <option value="0" type="string" default="true"> </option>');
+			fprintf(fid,'%s\n%s\n','       <option value="1" type="string" default="false"> </option>','</parameter>');
+
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="control_parameters" type="',class(obj.control_parameters),'" default="',convert2str(obj.control_parameters),'">','     <section name="inversion" />','     <help> ex: {''FrictionCoefficient''}, or {''MaterialsRheologyBbar''} </help>','  </parameter>');
+
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="nsteps" type="',class(obj.nsteps),'" default="',convert2str(obj.nsteps),'">','     <section name="inversion" />','     <help> number of optimization searches </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="cost_functions" type="',class(obj.cost_functions),'" default="',convert2str(obj.cost_functions),'">','     <section name="inversion" />','     <help> indicate the type of response for each optimization step  </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="cost_functions_coefficients" type="',class(obj.cost_functions_coefficients),'" default="',convert2str(obj.cost_functions_coefficients),'">','     <section name="inversion" />','     <help> cost_functions_coefficients applied to the misfit of each vertex and for each control_parameter </help>','  </parameter>');
+
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="cost_function_threshold" type="',class(obj.cost_function_threshold),'" default="',convert2str(obj.cost_function_threshold),'">','     <section name="inversion" />','     <help> misfit convergence criterion. Default is 1%, NaN if not applied </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="maxiter_per_step" type="',class(obj.maxiter_per_step),'" default="',convert2str(obj.maxiter_per_step),'">','     <section name="inversion" />','     <help> maximum iterations during each optimization step  </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="gradient_scaling" type="',class(obj.gradient_scaling),'" default="',convert2str(obj.gradient_scaling),'">','     <section name="inversion" />','     <help> scaling factor on gradient direction during optimization, for each optimization step </help>','  </parameter>');
+
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="step_threshold" type="',class(obj.step_threshold),'" default="',convert2str(obj.step_threshold),'">','     <section name="inversion" />','     <help> decrease threshold for misfit, default is 30% </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="min_parameters" type="',class(obj.min_parameters),'" default="',convert2str(obj.min_parameters),'">','     <section name="inversion" />','     <help> absolute minimum acceptable value of the inversed parameter on each vertex </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="max_parameters" type="',class(obj.max_parameters),'" default="',convert2str(obj.max_parameters),'">','     <section name="inversion" />','     <help> absolute maximum acceptable value of the inversed parameter on each vertex </help>','  </parameter>');
+
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vx_obs" type="',class(obj.vx_obs),'" default="',convert2str(obj.vx_obs),'">','     <section name="inversion" />','     <help> observed velocity x component [m/yr] </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vy_obs" type="',class(obj.vy_obs),'" default="',convert2str(obj.vy_obs),'">','     <section name="inversion" />','     <help> observed velocity y component [m/yr]  </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="vel_obs" type="',class(obj.vel_obs),'" default="',convert2str(obj.vel_obs),'">','     <section name="inversion" />','     <help> observed velocity magnitude [m/yr] </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="thickness_obs" type="',class(obj.thickness_obs),'" default="',convert2str(obj.thickness_obs),'">','     <section name="inversion" />','     <help> observed thickness [m]) </help>','  </parameter>');
+
+			fprintf(fid,'%s\n%s\n','</frame>');    
+
+			fprintf(fid,'%s\n%s\n%s\n','<frame key="2" label="Available cost functions">','<section name="inversion" />');                    
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceAbsVelMisfit" type="','string','" default="','101','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceRelVelMisfit" type="','string','" default="','102','">','     <section name="inversion" />','     <help>   </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceLogVelMisfit" type="','string','" default="','103','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
+
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceLogVxVyMisfit" type="','string','" default="','104','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="SurfaceAverageVelMisfit" type="','string','" default="','105','">','     <section name="inversion" />','     <help>   </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="ThicknessAbsMisfit" type="','string','" default="','106','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
+
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="DragCoefficientAbsGradient" type="','string','" default="','107','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="RheologyBbarAbsGradient" type="','string','" default="','108','">','     <section name="inversion" />','     <help>  </help>','  </parameter>');
+			fprintf(fid,'%s%s%s%s%s\n%s\n%s\n%s\n','  <parameter key ="ThicknessAbsGradient" type="','string','" default="','109','">','     <section name="inversion" />','     <help> </help>','  </parameter>');
+
+			fprintf(fid,'%s\n%s\n','</frame>');    
+
+		end % }}}       
 		function obj = inversion(varargin) % {{{
 			switch nargin
@@ -196,5 +196,5 @@
 			if ~obj.iscontrol, return; end
 			WriteData(fid,'object',obj,'fieldname','nsteps','format','Integer');
-			WriteData(fid,'object',obj,'fieldname','maxiter_per_step','format','DoubleMat','mattype',3);
+			WriteData(fid,'object',obj,'fieldname','maxiter_per_step','format','IntMat','mattype',3);
 			WriteData(fid,'object',obj,'fieldname','cost_functions_coefficients','format','DoubleMat','mattype',1);
 			WriteData(fid,'object',obj,'fieldname','gradient_scaling','format','DoubleMat','mattype',3);