Index: /issm/trunk/src/c/Makefile.am
===================================================================
--- /issm/trunk/src/c/Makefile.am	(revision 4047)
+++ /issm/trunk/src/c/Makefile.am	(revision 4048)
@@ -393,4 +393,6 @@
 					./modules/ScaleInputx/ScaleInputx.h\
 					./modules/ScaleInputx/ScaleInputx.cpp\
+					./modules/AXPYInputx/AXPYInputx.h\
+					./modules/AXPYInputx/AXPYInputx.cpp\
 					./modules/ControlConstrainx/ControlConstrainx.h\
 					./modules/ControlConstrainx/ControlConstrainx.cpp\
@@ -891,4 +893,6 @@
 					./modules/ScaleInputx/ScaleInputx.h\
 					./modules/ScaleInputx/ScaleInputx.cpp\
+					./modules/AXPYInputx/AXPYInputx.h\
+					./modules/AXPYInputx/AXPYInputx.cpp\
 					./modules/ControlConstrainx/ControlConstrainx.h\
 					./modules/ControlConstrainx/ControlConstrainx.cpp\
Index: /issm/trunk/src/c/modules/AXPYInputx/AXPYInputx.cpp
===================================================================
--- /issm/trunk/src/c/modules/AXPYInputx/AXPYInputx.cpp	(revision 4048)
+++ /issm/trunk/src/c/modules/AXPYInputx/AXPYInputx.cpp	(revision 4048)
@@ -0,0 +1,24 @@
+/*!\file AXPYInputx
+ * \brief: Y=Y+aX operation on inputs.
+ */
+
+#include "./AXPYInputx.h"
+#include "../../shared/shared.h"
+#include "../../include/include.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+
+void AXPYInputx(DataSet* elements,DataSet* nodes,DataSet* vertices,DataSet* loads,DataSet* materials,Parameters* parameters,int YEnum, double scalar, int XEnum){
+
+	/*intermediary:*/
+	int      i;
+
+	/*First, get elements and nodes configured: */
+	elements->Configure(elements,loads, nodes,vertices, materials,parameters);
+
+	/*Go through elemnets, and ask to carry out the AXPY operation on inputs: */
+	for(i=0;i<elements->Size();i++){
+		Element* element=(Element*)elements->GetObjectByOffset(i);
+		element->AXPYInput(YEnum, scalar, XEnum);
+	}
+}
Index: /issm/trunk/src/c/modules/AXPYInputx/AXPYInputx.h
===================================================================
--- /issm/trunk/src/c/modules/AXPYInputx/AXPYInputx.h	(revision 4048)
+++ /issm/trunk/src/c/modules/AXPYInputx/AXPYInputx.h	(revision 4048)
@@ -0,0 +1,14 @@
+/*!\file:  AXPYInputx.h
+ * \brief header file for field extrusion
+ */ 
+
+#ifndef _AXPYINPUTX_H
+#define _AXPYINPUTX_H
+
+#include "../../DataSet/DataSet.h"
+
+/* local prototypes: */
+void AXPYInputx(DataSet* elements,DataSet* nodes,DataSet* vertices,DataSet* loads,DataSet* materials,Parameters* parameters,int YEnum, double scalar, int XEnum){
+
+#endif  /* _AXPYINPUTX_H */
+
Index: /issm/trunk/src/c/modules/ControlConstrainx/ControlConstrainx.cpp
===================================================================
--- /issm/trunk/src/c/modules/ControlConstrainx/ControlConstrainx.cpp	(revision 4047)
+++ /issm/trunk/src/c/modules/ControlConstrainx/ControlConstrainx.cpp	(revision 4048)
@@ -1,8 +1,7 @@
-/*!\file ControlConstrainx
- * \brief constrain optimization parameters during inverse method
+/*!\file ControlContrainInputx
+ * \brief: Y=Y+aX operation on inputs.
  */
 
-#include "./ControlConstrainx.h"
-
+#include "./ControlContrainInputx.h"
 #include "../../shared/shared.h"
 #include "../../include/include.h"
@@ -10,24 +9,19 @@
 #include "../../EnumDefinitions/EnumDefinitions.h"
 
-void ControlConstrainx( double* p_g, int numdofnodes, double cm_min, double cm_max,int control_type){
+void ControlContrainInputx(DataSet* elements,DataSet* nodes,DataSet* vertices,DataSet* loads,DataSet* materials,Parameters* parameters,int control_type,double cm_min, double cm_max){
 
-	int i;
-	
-	if(isnan(cm_min) & isnan(cm_max)){
-		/*do nothing*/
+	/*intermediary:*/
+	int      i;
+
+    /*some early returns: */
+	if(isnan(cm_min) & isnan(cm_max))return;
+
+	/*First, get elements and nodes configured: */
+	elements->Configure(elements,loads, nodes,vertices, materials,parameters);
+
+	/*Go through elemnets, and ask to carry out the ControlContrain operation on inputs: */
+	for(i=0;i<elements->Size();i++){
+		Element* element=(Element*)elements->GetObjectByOffset(i);
+		element->ControlContrainInput(control_type,cm_min,cm_max);
 	}
-	else{
-		for(i=0;i<numdofnodes;i++){
-			if(isnan(p_g[i])){
-				ISSMERROR("NaN found in parameter p_g[%i]",i);
-			}
-			if(!isnan(cm_min)){
-				if (p_g[i]<cm_min)p_g[i]=cm_min;
-			}
-			if(!isnan(cm_max)){
-				if (p_g[i]>cm_max)p_g[i]=cm_max;
-			}
-		}
-	}
-		
 }
Index: /issm/trunk/src/c/modules/ControlConstrainx/ControlConstrainx.h
===================================================================
--- /issm/trunk/src/c/modules/ControlConstrainx/ControlConstrainx.h	(revision 4047)
+++ /issm/trunk/src/c/modules/ControlConstrainx/ControlConstrainx.h	(revision 4048)
@@ -1,14 +1,14 @@
-/*!\file:  ControlConstrainx.h
- * \brief constrain optimization parameters during inverse method
+/*!\file:  ControlContrainInputx.h
+ * \brief header file for field extrusion
  */ 
 
-#ifndef _CONTROLCONSTRAINX_H
-#define _CONTROLCONSTRAINX_H
+#ifndef _CONTROLCONTRAININPUTX_H
+#define _CONTROLCONTRAININPUTX_H
 
 #include "../../DataSet/DataSet.h"
 
 /* local prototypes: */
-void ControlConstrainx( double* p_g, int numberofnodes, double cm_min, double cm_max,int control_type);
+void ControlContrainInputx(DataSet* elements,DataSet* nodes,DataSet* vertices,DataSet* loads,DataSet* materials,Parameters* parameters,int control_type,double cm_min, double cm_max);
 
-#endif  /* _CONTROLCONSTRAINX_H */
+#endif  /* _ControlContrainINPUTX_H */
 
Index: /issm/trunk/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp
===================================================================
--- /issm/trunk/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp	(revision 4048)
+++ /issm/trunk/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp	(revision 4048)
@@ -0,0 +1,39 @@
+/*!\file GetVectorFromInputsx
+ * \brief retrieve vector from inputs in elements
+ */
+
+#include "./GetVectorFromInputsx.h"
+#include "../../shared/shared.h"
+#include "../../include/include.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+
+void GetVectorFromInputsx( Vec* pvector, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters, int NameEnum, int TypeEnum){
+
+
+	Vec vector=NULL;
+
+	/*First, get elements and loads configured: */
+	elements->Configure(elements,loads, nodes,vertices, materials,parameters);
+	loads->Configure(elements, loads, nodes,vertices, materials,parameters);
+	nodes->Configure(elements, loads, nodes,vertices, materials,parameters);
+	parameters->Configure(elements,loads, nodes,vertices, materials,parameters);
+
+
+	if(TypeEnum==VertexEnum){
+
+		vector=NewVec(vertices->Size());
+		for(i=0;i<elements->Size();i++){
+			Element* element=(Element*)elements->GetObjectByOffset(i);
+			element->GetVectorFromInputs(vector,NameEnum);
+		}
+	}
+	else ISSMERROR("%s%s%s"," vector type: ",EnumAsString(TypeEnum)," not supported yet!");
+
+	VecAssemblyBegin(vector);
+	VecAssemblyEnd(vector);
+
+	/*Assign output pointers:*/
+	*pvector=vector;
+
+}
Index: /issm/trunk/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.h
===================================================================
--- /issm/trunk/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.h	(revision 4048)
+++ /issm/trunk/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.h	(revision 4048)
@@ -0,0 +1,14 @@
+/*!\file:  GetVectorFromInputsx.h
+ */ 
+
+#ifndef _GETVECTORFROMINPUTSXX_H
+#define _GETVECTORFROMINPUTSXX_H
+
+#include "../../objects/objects.h"
+#include "../../DataSet/DataSet.h"
+
+/* local prototypes: */
+void	GetVectorFromInputsx( Vec* pvector, DataSet* elements,DataSet* nodes, DataSet* vertices,DataSet* loads, DataSet* materials,  Parameters* parameters,int NameEnum,int TypeEnum);
+
+#endif  /* _GETVECTORFROMINPUTSXX_H */
+
Index: /issm/trunk/src/c/modules/Gradjx/Gradjx.cpp
===================================================================
--- /issm/trunk/src/c/modules/Gradjx/Gradjx.cpp	(revision 4047)
+++ /issm/trunk/src/c/modules/Gradjx/Gradjx.cpp	(revision 4048)
@@ -24,5 +24,4 @@
 	parameters->FindParam(&dim,DimEnum);
 
-
 	/*Compute gradients: */
 	for (i=0;i<elements->Size();i++){
Index: /issm/trunk/src/c/modules/Orthx/Orthx.cpp
===================================================================
--- /issm/trunk/src/c/modules/Orthx/Orthx.cpp	(revision 4047)
+++ /issm/trunk/src/c/modules/Orthx/Orthx.cpp	(revision 4048)
@@ -6,4 +6,6 @@
 
 void	Orthx( Vec* pnewgradj, Vec gradj, Vec oldgradj){
+	
+	GetSolutionFromInputsx( &ug, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,  femmodel->parameters, analysis_type,sub_analysis_type);
 
 	/*output: */
Index: /issm/trunk/src/c/modules/modules.h
===================================================================
--- /issm/trunk/src/c/modules/modules.h	(revision 4047)
+++ /issm/trunk/src/c/modules/modules.h	(revision 4048)
@@ -80,4 +80,6 @@
 #include "./DuplicateInputx/DuplicateInputx.h"
 #include "./ScaleInputx/ScaleInputx.h"
+#include "./AXPYInputx/AXPYInputx.h"
+#include "./ControlContrainx/ControlContrainx.h"
 
 #endif
Index: /issm/trunk/src/c/objects/Elements/Beam.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Beam.cpp	(revision 4047)
+++ /issm/trunk/src/c/objects/Elements/Beam.cpp	(revision 4048)
@@ -1002,2 +1002,54 @@
 }
 /*}}}*/
+/*FUNCTION Beam::AXPY(int YEnum, double scalar, int XEnum);{{{1*/
+void  Beam::AXPY(int YEnum, double scalar, int XEnum){
+
+	Input* xinput=NULL;
+	Input* yinput=NULL;
+
+	/*Find x and y inputs: */
+	xinput=(Input*)this->inputs->GetInput(XEnum);
+	yinput=(Input*)this->inputs->GetInput(YEnum);
+
+	/*some checks: */
+	if(!xinput || !yinput)ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," or input ",EnumAsString(YEnum)," could not be found!");
+	if(xinput->Enum()!=yinput->Enum())ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," and input ",EnumAsString(YEnum)," are not of the same type!");
+
+	/*Scale: */
+	yinput->AXPY(xinput,scalar);
+}
+/*}}}*/
+/*FUNCTION Beam::ControlConstrain(int control_type, double cm_min, double cm_max){{{1*/
+void  Beam::ControlConstrain(int control_type, double cm_min, double cm_max){
+
+	Input* input=NULL;
+
+	/*Find input: */
+	input=(Input*)this->inputs->GetInput(control_type);
+	
+	/*Do nothing if we  don't find it: */
+	if(!input)return;
+
+	/*Constrain input using cm_min and cm_max: */
+	input->Constrain(cm_min,cm_max);
+
+}
+/*}}}*/
+/*FUNCTION Beam::GetVectorFromInputs(Vec vector,int NameEnum){{{1*/
+void  Beam::GetVectorFromInputs(Vec vector,int NameEnum){
+
+	const int numvertices=2;
+	int doflist1[numvertices];
+
+	/*Find NameEnum input in the inputs dataset, and get it to fill in the vector: */
+	for(i=0;i<this->inputs->Size();i++){
+		Input* input=(Input*)this->inputs->GetObjectByOffset(i);
+		if(input->EnumType==NameEnum){
+			/*We found the enum.  Use its values to fill into the vector, using the vertices ids: */
+			this->GetDofList1(&doflist1[0]);
+			input->GetVectorFromInputs(vector,&doflist1[0]);
+			break;
+		}
+	}
+}
+/*}}}*/
Index: /issm/trunk/src/c/objects/Elements/Beam.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Beam.h	(revision 4047)
+++ /issm/trunk/src/c/objects/Elements/Beam.h	(revision 4048)
@@ -96,4 +96,6 @@
 		void  DuplicateInput(int original_enum,int new_enum);
 		void  ScaleInput(int enum_type,double scale_factor);
+		void  AXPY(int YEnum, double scalar, int XEnum);
+		void  ControlConstrain(int control_type,double cm_min, double cm_max);
 
 		/*}}}*/
@@ -114,4 +116,6 @@
 		void  GetJacobianDeterminant(double* pJdet,double* z_list, double gauss_coord);
 		double MassFlux(double* segment);
+		void GetVectorFromInputs(Vec vector,int NameEnum);
+
 		/*}}}*/
 
Index: /issm/trunk/src/c/objects/Elements/Element.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Element.h	(revision 4047)
+++ /issm/trunk/src/c/objects/Elements/Element.h	(revision 4048)
@@ -66,6 +66,7 @@
 		virtual void   DuplicateInput(int original_enum,int new_enum)=0;
 		virtual void   ScaleInput(int enum_type,double scale_factor)=0;
-
-
+		virtual void   GetVectorFromInputs(Vec vector,int NameEnum)=0;
+		virtual void   AXPY(int YEnum, double scalar, int XEnum)=0;
+		virtual void   ControlConstrain(int control_type,double cm_min, double cm_max)=0;
 		/*Implementation: */
 
Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 4047)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 4048)
@@ -5352,2 +5352,54 @@
 }
 /*}}}*/
+/*FUNCTION Penta::AXPY(int YEnum, double scalar, int XEnum);{{{1*/
+void  Penta::AXPY(int YEnum, double scalar, int XEnum){
+
+	Input* xinput=NULL;
+	Input* yinput=NULL;
+
+	/*Find x and y inputs: */
+	xinput=(Input*)this->inputs->GetInput(XEnum);
+	yinput=(Input*)this->inputs->GetInput(YEnum);
+
+	/*some checks: */
+	if(!xinput || !yinput)ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," or input ",EnumAsString(YEnum)," could not be found!");
+	if(xinput->Enum()!=yinput->Enum())ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," and input ",EnumAsString(YEnum)," are not of the same type!");
+
+	/*Scale: */
+	yinput->AXPY(xinput,scalar);
+}
+/*}}}*/
+/*FUNCTION Penta::ControlConstrain(int control_type, double cm_min, double cm_max){{{1*/
+void  Penta::ControlConstrain(int control_type, double cm_min, double cm_max){
+
+	Input* input=NULL;
+
+	/*Find input: */
+	input=(Input*)this->inputs->GetInput(control_type);
+	
+	/*Do nothing if we  don't find it: */
+	if(!input)return;
+
+	/*Constrain input using cm_min and cm_max: */
+	input->Constrain(cm_min,cm_max);
+
+}
+/*}}}*/
+/*FUNCTION Penta::GetVectorFromInputs(Vec vector,int NameEnum){{{1*/
+void  Penta::GetVectorFromInputs(Vec vector,int NameEnum){
+
+	const int numvertices=6;
+	int doflist1[numvertices];
+
+	/*Find NameEnum input in the inputs dataset, and get it to fill in the vector: */
+	for(i=0;i<this->inputs->Size();i++){
+		Input* input=(Input*)this->inputs->GetObjectByOffset(i);
+		if(input->EnumType==NameEnum){
+			/*We found the enum.  Use its values to fill into the vector, using the vertices ids: */
+			this->GetDofList1(&doflist1[0]);
+			input->GetVectorFromInputs(vector,&doflist1[0]);
+			break;
+		}
+	}
+}
+/*}}}*/
Index: /issm/trunk/src/c/objects/Elements/Penta.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.h	(revision 4047)
+++ /issm/trunk/src/c/objects/Elements/Penta.h	(revision 4048)
@@ -163,4 +163,7 @@
 		void  DuplicateInput(int original_enum,int new_enum);
 		void  ScaleInput(int enum_type,double scale_factor);
+		void  AXPY(int YEnum, double scalar, int XEnum);
+		void  ControlConstrain(int control_type,double cm_min, double cm_max);
+		void  GetVectorFromInputs(Vec vector,int NameEnum);
 
 
Index: /issm/trunk/src/c/objects/Elements/Sing.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Sing.cpp	(revision 4047)
+++ /issm/trunk/src/c/objects/Elements/Sing.cpp	(revision 4048)
@@ -708,2 +708,54 @@
 }
 /*}}}*/
+/*FUNCTION Sing::AXPY(int YEnum, double scalar, int XEnum);{{{1*/
+void  Sing::AXPY(int YEnum, double scalar, int XEnum){
+
+	Input* xinput=NULL;
+	Input* yinput=NULL;
+
+	/*Find x and y inputs: */
+	xinput=(Input*)this->inputs->GetInput(XEnum);
+	yinput=(Input*)this->inputs->GetInput(YEnum);
+
+	/*some checks: */
+	if(!xinput || !yinput)ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," or input ",EnumAsString(YEnum)," could not be found!");
+	if(xinput->Enum()!=yinput->Enum())ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," and input ",EnumAsString(YEnum)," are not of the same type!");
+
+	/*Scale: */
+	yinput->AXPY(xinput,scalar);
+}
+/*}}}*/
+/*FUNCTION Sing::ControlConstrain(int control_type, double cm_min, double cm_max){{{1*/
+void  Sing::ControlConstrain(int control_type, double cm_min, double cm_max){
+
+	Input* input=NULL;
+
+	/*Find input: */
+	input=(Input*)this->inputs->GetInput(control_type);
+	
+	/*Do nothing if we  don't find it: */
+	if(!input)return;
+
+	/*Constrain input using cm_min and cm_max: */
+	input->Constrain(cm_min,cm_max);
+
+}
+/*}}}*/
+/*FUNCTION Sing::GetVectorFromInputs(Vec vector,int NameEnum){{{1*/
+void  Sing::GetVectorFromInputs(Vec vector,int NameEnum){
+
+	const int numvertices=1;
+	int doflist1[numvertices];
+
+	/*Find NameEnum input in the inputs dataset, and get it to fill in the vector: */
+	for(i=0;i<this->inputs->Size();i++){
+		Input* input=(Input*)this->inputs->GetObjectByOffset(i);
+		if(input->EnumType==NameEnum){
+			/*We found the enum.  Use its values to fill into the vector, using the vertices ids: */
+			this->GetDofList1(&doflist1[0]);
+			input->GetVectorFromInputs(vector,&doflist1[0]);
+			break;
+		}
+	}
+}
+/*}}}*/
Index: /issm/trunk/src/c/objects/Elements/Sing.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Sing.h	(revision 4047)
+++ /issm/trunk/src/c/objects/Elements/Sing.h	(revision 4048)
@@ -95,4 +95,6 @@
 		void  DuplicateInput(int original_enum,int new_enum);
 		void  ScaleInput(int enum_type,double scale_factor);
+		void  AXPY(int YEnum, double scalar, int XEnum);
+		void  ControlConstrain(int control_type,double cm_min, double cm_max);
 
 
@@ -111,4 +113,5 @@
 		double CostFunction(void);
 		double MassFlux(double* segment);
+		void  GetVectorFromInputs(Vec vector,int NameEnum);
 		/*}}}*/
 
Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 4047)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 4048)
@@ -5450,2 +5450,54 @@
 }
 /*}}}*/
+/*FUNCTION Tria::AXPY(int YEnum, double scalar, int XEnum);{{{1*/
+void  Tria::AXPY(int YEnum, double scalar, int XEnum){
+
+	Input* xinput=NULL;
+	Input* yinput=NULL;
+
+	/*Find x and y inputs: */
+	xinput=(Input*)this->inputs->GetInput(XEnum);
+	yinput=(Input*)this->inputs->GetInput(YEnum);
+
+	/*some checks: */
+	if(!xinput || !yinput)ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," or input ",EnumAsString(YEnum)," could not be found!");
+	if(xinput->Enum()!=yinput->Enum())ISSMERROR("%s%s%s%s%s"," input ",EnumAsString(XEnum)," and input ",EnumAsString(YEnum)," are not of the same type!");
+
+	/*Scale: */
+	yinput->AXPY(xinput,scalar);
+}
+/*}}}*/
+/*FUNCTION Tria::ControlConstrain(int control_type, double cm_min, double cm_max){{{1*/
+void  Tria::ControlConstrain(int control_type, double cm_min, double cm_max){
+
+	Input* input=NULL;
+
+	/*Find input: */
+	input=(Input*)this->inputs->GetInput(control_type);
+	
+	/*Do nothing if we  don't find it: */
+	if(!input)return;
+
+	/*Constrain input using cm_min and cm_max: */
+	input->Constrain(cm_min,cm_max);
+
+}
+/*}}}*/
+/*FUNCTION Tria::GetVectorFromInputs(Vec vector,int NameEnum){{{1*/
+void  Tria::GetVectorFromInputs(Vec vector,int NameEnum){
+
+	const int numvertices=3;
+	int doflist1[numvertices];
+
+	/*Find NameEnum input in the inputs dataset, and get it to fill in the vector: */
+	for(i=0;i<this->inputs->Size();i++){
+		Input* input=(Input*)this->inputs->GetObjectByOffset(i);
+		if(input->EnumType==NameEnum){
+			/*We found the enum.  Use its values to fill into the vector, using the vertices ids: */
+			this->GetDofList1(&doflist1[0]);
+			input->GetVectorFromInputs(vector,&doflist1[0]);
+			break;
+		}
+	}
+}
+/*}}}*/
Index: /issm/trunk/src/c/objects/Elements/Tria.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.h	(revision 4047)
+++ /issm/trunk/src/c/objects/Elements/Tria.h	(revision 4048)
@@ -141,4 +141,7 @@
 		void  DuplicateInput(int original_enum,int new_enum);
 		void  ScaleInput(int enum_type,double scale_factor);
+		void  AXPY(int YEnum, double scalar, int XEnum);
+		void  ControlConstrain(int control_type,double cm_min, double cm_max);
+		void  GetVectorFromInputs(Vec vector,int NameEnum);
 
 
Index: /issm/trunk/src/c/objects/Inputs/BeamVertexInput.cpp
===================================================================
--- /issm/trunk/src/c/objects/Inputs/BeamVertexInput.cpp	(revision 4047)
+++ /issm/trunk/src/c/objects/Inputs/BeamVertexInput.cpp	(revision 4048)
@@ -259,2 +259,37 @@
 }
 /*}}}*/
+/*FUNCTION BeamVertexInput::AXPY(Input* xinput,double scalar);{{{1*/
+void BeamVertexInput::AXPY(Input* xinput,double scalar){
+
+	int i;
+	const int numgrids=2;
+	BeamVertexInput*  xbeamvertexinput=NULL;
+
+	/*xinput is of the same type, so cast it: */
+	xbeamvertexinput=(BeamVertexInput)xinput;
+
+	/*Carry out the AXPY operation:*/
+	for(i=0;i<numgrids;i++)this->values[i]=this->values[i]+scalar*xbeamvertexinput->values[i];
+
+}
+/*}}}*/
+/*FUNCTION BeamVertexInput::Constrain(double cm_min, double cm_max){{{1*/
+void BeamVertexInput::Constrain(double cm_min, double cm_max){
+
+	int i;
+	const int numgrids=2;
+		
+	if(!isnan(cm_min)) for(i=0;i<numgrids;i++)if (this->values[i]<cm_min)this->values[i]=cm_min;
+	if(!isnan(cm_max)) for(i=0;i<numgrids;i++)if (this->values[i]>cm_max)this->values[i]=cm_max;
+
+}
+/*}}}*/
+/*FUNCTION BeamVertexInput::GetVectorFromInputs(Vec vector,int* doflist){{{1*/
+void BeamVertexInput::GetVectorFromInputs(Vec vector,int* doflist){
+
+	const int numvertices=2;
+	VecSetValues(vector,numvertices,doflist,(const double*)this->values,ADD_VALUES);
+
+
+}
+/*}}}*/
Index: /issm/trunk/src/c/objects/Inputs/BeamVertexInput.h
===================================================================
--- /issm/trunk/src/c/objects/Inputs/BeamVertexInput.h	(revision 4047)
+++ /issm/trunk/src/c/objects/Inputs/BeamVertexInput.h	(revision 4048)
@@ -75,4 +75,7 @@
 		void SquareMin(double* psquaremin, bool process_units,Parameters* parameters);
 		void Scale(double scale_factor);
+		void AXPY(Input* xinput,double scalar);
+		void Constrain(double cm_min, double cm_max);
+		void GetVectorFromInputs(Vec vector,int* doflist);
 		/*}}}*/
 
Index: /issm/trunk/src/c/objects/Inputs/BoolInput.cpp
===================================================================
--- /issm/trunk/src/c/objects/Inputs/BoolInput.cpp	(revision 4047)
+++ /issm/trunk/src/c/objects/Inputs/BoolInput.cpp	(revision 4048)
@@ -228,2 +228,30 @@
 }
 /*}}}*/
+/*FUNCTION BoolInput::AXPY(Input* xinput,double scalar);{{{1*/
+void BoolInput::AXPY(Input* xinput,double scalar){
+
+	BoolInput*  xboolinput=NULL;
+
+	/*xinput is of the same type, so cast it: */
+	xboolinput=(BoolInput)xinput;
+
+	/*Carry out the AXPY operation:*/
+	this->value=this->value+scalar*xboolinput->value;
+
+}
+/*}}}*/
+/*FUNCTION BoolInput::Constrain(double cm_min, double cm_max){{{1*/
+void BoolInput::Constrain(double cm_min, double cm_max){
+
+	if(!isnan(cm_min)) if (this->value<cm_min)this->value=cm_min;
+	if(!isnan(cm_max)) if (this->value>cm_max)this->value=cm_max;
+
+}
+/*}}}*/
+/*FUNCTION BoolInput::GetVectorFromInputs(Vec vector,int* doflist){{{1*/
+void BoolInput::GetVectorFromInputs(Vec vector,int* doflist){
+
+	ISSMERROR(" not supporte yet!");
+
+}
+/*}}}*/
Index: /issm/trunk/src/c/objects/Inputs/BoolInput.h
===================================================================
--- /issm/trunk/src/c/objects/Inputs/BoolInput.h	(revision 4047)
+++ /issm/trunk/src/c/objects/Inputs/BoolInput.h	(revision 4048)
@@ -75,4 +75,7 @@
 		void SquareMin(double* psquaremin, bool process_units,Parameters* parameters);
 		void Scale(double scale_factor);
+		void AXPY(Input* xinput,double scalar);
+		void Constrain(double cm_min, double cm_max);
+		void GetVectorFromInputs(Vec vector,int* doflist);
 		/*}}}*/
 
Index: /issm/trunk/src/c/objects/Inputs/DoubleInput.cpp
===================================================================
--- /issm/trunk/src/c/objects/Inputs/DoubleInput.cpp	(revision 4047)
+++ /issm/trunk/src/c/objects/Inputs/DoubleInput.cpp	(revision 4048)
@@ -239,2 +239,30 @@
 }
 /*}}}*/
+/*FUNCTION DoubleInput::AXPY(Input* xinput,double scalar);{{{1*/
+void DoubleInput::AXPY(Input* xinput,double scalar){
+
+	DoubleInput*  xdoubleinput=NULL;
+
+	/*xinput is of the same type, so cast it: */
+	xdoubleinput=(DoubleInput)xinput;
+
+	/*Carry out the AXPY operation:*/
+	this->value=this->value+scalar*xdoubleinput->value;
+
+}
+/*}}}*/
+/*FUNCTION DoubleInput::Constrain(double cm_min, double cm_max){{{1*/
+void DoubleInput::Constrain(double cm_min, double cm_max){
+
+	if(!isnan(cm_min)) if (this->value<cm_min)this->value=cm_min;
+	if(!isnan(cm_max)) if (this->value>cm_max)this->value=cm_max;
+
+}
+/*}}}*/
+/*FUNCTION DoubleInput::GetVectorFromInputs(Vec vector,int* doflist){{{1*/
+void DoubleInput::GetVectorFromInputs(Vec vector,int* doflist){
+
+	ISSMERROR(" not supporte yet!");
+
+}
+/*}}}*/
Index: /issm/trunk/src/c/objects/Inputs/DoubleInput.h
===================================================================
--- /issm/trunk/src/c/objects/Inputs/DoubleInput.h	(revision 4047)
+++ /issm/trunk/src/c/objects/Inputs/DoubleInput.h	(revision 4048)
@@ -76,4 +76,7 @@
 		void SquareMin(double* psquaremin, bool process_units,Parameters* parameters);
 		void Scale(double scale_factor);
+		void AXPY(Input* xinput,double scalar);
+		void Constrain(double cm_min, double cm_max);
+		void GetVectorFromInputs(Vec vector,int* doflist);
 		/*}}}*/
 
Index: /issm/trunk/src/c/objects/Inputs/Input.h
===================================================================
--- /issm/trunk/src/c/objects/Inputs/Input.h	(revision 4047)
+++ /issm/trunk/src/c/objects/Inputs/Input.h	(revision 4048)
@@ -50,6 +50,9 @@
 		virtual Input* SpawnTriaInput(int* indices)=0;
 		virtual Result* SpawnResult(int step, double time)=0;
-		virtual void   SquareMin(double* psquaremin, bool process_units,Parameters* parameters)=0;
-		virtual void   Scale(double scale_factor)=0;
+		virtual void SquareMin(double* psquaremin, bool process_units,Parameters* parameters)=0;
+		virtual void Scale(double scale_factor)=0;
+		virtual void AXPY(Input* xinput,double scalar)=0;
+		virtual void Constrain(double cm_min, double cm_max)=0;
+		virtual void GetVectorFromInputs(Vec vector,int* doflist)=0;
 		/*}}}*/
 
Index: /issm/trunk/src/c/objects/Inputs/IntInput.cpp
===================================================================
--- /issm/trunk/src/c/objects/Inputs/IntInput.cpp	(revision 4047)
+++ /issm/trunk/src/c/objects/Inputs/IntInput.cpp	(revision 4048)
@@ -226,2 +226,30 @@
 }
 /*}}}*/
+/*FUNCTION IntInput::AXPY(Input* xinput,int scalar);{{{1*/
+void IntInput::AXPY(Input* xinput,int scalar){
+
+	IntInput*  xintinput=NULL;
+
+	/*xinput is of the same type, so cast it: */
+	xintinput=(IntInput)xinput;
+
+	/*Carry out the AXPY operation:*/
+	this->value=this->value+scalar*xintinput->value;
+
+}
+/*}}}*/
+/*FUNCTION IntInput::Constrain(double cm_min, double cm_max){{{1*/
+void IntInput::Constrain(double cm_min, double cm_max){
+
+	if(!isnan(cm_min)) if (this->value<cm_min)this->value=cm_min;
+	if(!isnan(cm_max)) if (this->value>cm_max)this->value=cm_max;
+
+}
+/*}}}*/
+/*FUNCTION IntInput::GetVectorFromInputs(Vec vector,int* doflist){{{1*/
+void IntInput::GetVectorFromInputs(Vec vector,int* doflist){
+
+	ISSMERROR(" not supporte yet!");
+
+}
+/*}}}*/
Index: /issm/trunk/src/c/objects/Inputs/IntInput.h
===================================================================
--- /issm/trunk/src/c/objects/Inputs/IntInput.h	(revision 4047)
+++ /issm/trunk/src/c/objects/Inputs/IntInput.h	(revision 4048)
@@ -75,4 +75,7 @@
 		void SquareMin(double* psquaremin, bool process_units,Parameters* parameters);
 		void Scale(double scale_factor);
+		void AXPY(Input* xinput,double scalar);
+		void Constrain(double cm_min, double cm_max);
+		void GetVectorFromInputs(Vec vector,int* doflist);
 		/*}}}*/
 
Index: /issm/trunk/src/c/objects/Inputs/PentaVertexInput.cpp
===================================================================
--- /issm/trunk/src/c/objects/Inputs/PentaVertexInput.cpp	(revision 4047)
+++ /issm/trunk/src/c/objects/Inputs/PentaVertexInput.cpp	(revision 4048)
@@ -908,2 +908,37 @@
 }
 /*}}}*/
+/*FUNCTION PentaVertexInput::AXPY(Input* xinput,double scalar);{{{1*/
+void PentaVertexInput::AXPY(Input* xinput,double scalar){
+
+	int i;
+	const int numgrids=6;
+	PentaVertexInput*  xpentavertexinput=NULL;
+
+	/*xinput is of the same type, so cast it: */
+	xpentavertexinput=(PentaVertexInput)xinput;
+
+	/*Carry out the AXPY operation:*/
+	for(i=0;i<numgrids;i++)this->values[i]=this->values[i]+scalar*xpentavertexinput->values[i];
+
+}
+/*}}}*/
+/*FUNCTION PentaVertexInput::Constrain(double cm_min, double cm_max){{{1*/
+void PentaVertexInput::Constrain(double cm_min, double cm_max){
+
+	int i;
+	const int numgrids=6;
+		
+	if(!isnan(cm_min)) for(i=0;i<numgrids;i++)if (this->values[i]<cm_min)this->values[i]=cm_min;
+	if(!isnan(cm_max)) for(i=0;i<numgrids;i++)if (this->values[i]>cm_max)this->values[i]=cm_max;
+
+}
+/*}}}*/
+/*FUNCTION PentaVertexInput::GetVectorFromInputs(Vec vector,int* doflist){{{1*/
+void PentaVertexInput::GetVectorFromInputs(Vec vector,int* doflist){
+
+	const int numvertices=6;
+	VecSetValues(vector,numvertices,doflist,(const double*)this->values,ADD_VALUES);
+
+
+}
+/*}}}*/
Index: /issm/trunk/src/c/objects/Inputs/PentaVertexInput.h
===================================================================
--- /issm/trunk/src/c/objects/Inputs/PentaVertexInput.h	(revision 4047)
+++ /issm/trunk/src/c/objects/Inputs/PentaVertexInput.h	(revision 4048)
@@ -84,4 +84,7 @@
 		void SquareMin(double* psquaremin, bool process_units,Parameters* parameters);
 		void Scale(double scale_factor);
+		void AXPY(Input* xinput,double scalar);
+		void Constrain(double cm_min, double cm_max);
+		void GetVectorFromInputs(Vec vector,int* doflist);
 		/*}}}*/
 
Index: /issm/trunk/src/c/objects/Inputs/SingVertexInput.cpp
===================================================================
--- /issm/trunk/src/c/objects/Inputs/SingVertexInput.cpp	(revision 4047)
+++ /issm/trunk/src/c/objects/Inputs/SingVertexInput.cpp	(revision 4048)
@@ -229,2 +229,32 @@
 }
 /*}}}*/
+/*FUNCTION SingVertexInput::AXPY(Input* xinput,double scalar);{{{1*/
+void SingVertexInput::AXPY(Input* xinput,double scalar){
+
+	SingVertexInput*  xsingvertexinput=NULL;
+
+	/*xinput is of the same type, so cast it: */
+	xsingvertexinput=(SingVertexInput)xinput;
+
+	/*Carry out the AXPY operation:*/
+	this->value=this->value+scalar*xngvertexinput->value;
+
+}
+/*}}}*/
+/*FUNCTION SingVertexInput::Constrain(double cm_min, double cm_max){{{1*/
+void SingVertexInput::Constrain(double cm_min, double cm_max){
+
+	if(!isnan(cm_min)) if (this->value<cm_min)this->value=cm_min;
+	if(!isnan(cm_max)) if (this->value>cm_max)this->value=cm_max;
+
+}
+/*}}}*/
+/*FUNCTION SingVertexInput::GetVectorFromInputs(Vec vector,int* doflist){{{1*/
+void SingVertexInput::GetVectorFromInputs(Vec vector,int* doflist){
+
+	const int numvertices=1;
+	VecSetValues(vector,numvertices,doflist,(const double*)&this->value,ADD_VALUES);
+
+
+}
+/*}}}*/
Index: /issm/trunk/src/c/objects/Inputs/SingVertexInput.h
===================================================================
--- /issm/trunk/src/c/objects/Inputs/SingVertexInput.h	(revision 4047)
+++ /issm/trunk/src/c/objects/Inputs/SingVertexInput.h	(revision 4048)
@@ -74,4 +74,7 @@
 		void SquareMin(double* psquaremin, bool process_units,Parameters* parameters);
 		void Scale(double scale_factor);
+		void AXPY(Input* xinput,double scalar);
+		void Constrain(double cm_min, double cm_max);
+		void GetVectorFromInputs(Vec vector,int* doflist);
 		/*}}}*/
 
Index: /issm/trunk/src/c/objects/Inputs/TriaVertexInput.cpp
===================================================================
--- /issm/trunk/src/c/objects/Inputs/TriaVertexInput.cpp	(revision 4047)
+++ /issm/trunk/src/c/objects/Inputs/TriaVertexInput.cpp	(revision 4048)
@@ -482,2 +482,37 @@
 }
 /*}}}*/
+/*FUNCTION TriaVertexInput::AXPY(Input* xinput,double scalar);{{{1*/
+void TriaVertexInput::AXPY(Input* xinput,double scalar){
+
+	int i;
+	const int numgrids=3;
+	TriaVertexInput*  xtriavertexinput=NULL;
+
+	/*xinput is of the same type, so cast it: */
+	xtriavertexinput=(TriaVertexInput)xinput;
+
+	/*Carry out the AXPY operation:*/
+	for(i=0;i<numgrids;i++)this->values[i]=this->values[i]+scalar*xtriavertexinput->values[i];
+
+}
+/*}}}*/
+/*FUNCTION TriaVertexInput::Constrain(double cm_min, double cm_max){{{1*/
+void TriaVertexInput::Constrain(double cm_min, double cm_max){
+
+	int i;
+	const int numgrids=3;
+		
+	if(!isnan(cm_min)) for(i=0;i<numgrids;i++)if (this->values[i]<cm_min)this->values[i]=cm_min;
+	if(!isnan(cm_max)) for(i=0;i<numgrids;i++)if (this->values[i]>cm_max)this->values[i]=cm_max;
+
+}
+/*}}}*/
+/*FUNCTION TriaVertexInput::GetVectorFromInputs(Vec vector,int* doflist){{{1*/
+void TriaVertexInput::GetVectorFromInputs(Vec vector,int* doflist){
+
+	const int numvertices=3;
+	VecSetValues(vector,numvertices,doflist,(const double*)this->values,ADD_VALUES);
+
+
+}
+/*}}}*/
Index: /issm/trunk/src/c/objects/Inputs/TriaVertexInput.h
===================================================================
--- /issm/trunk/src/c/objects/Inputs/TriaVertexInput.h	(revision 4047)
+++ /issm/trunk/src/c/objects/Inputs/TriaVertexInput.h	(revision 4048)
@@ -81,4 +81,7 @@
 		void SquareMin(double* psquaremin, bool process_units,Parameters* parameters);
 		void Scale(double scale_factor);
+		void AXPY(Input* xinput,double scalar);
+		void Constrain(double cm_min, double cm_max);
+		void GetVectorFromInputs(Vec vector,int* doflist);
 		/*}}}*/
 
Index: /issm/trunk/src/c/objects/OptArgs.h
===================================================================
--- /issm/trunk/src/c/objects/OptArgs.h	(revision 4047)
+++ /issm/trunk/src/c/objects/OptArgs.h	(revision 4048)
@@ -25,7 +25,5 @@
 class Model;
 struct OptArgs{
-	Model* model;
-	double* param_g;
-	double* grad_g;
+	FemModel* femmodel;
 	int n;
 };
Index: /issm/trunk/src/c/solutions/ControlRestart.cpp
===================================================================
--- /issm/trunk/src/c/solutions/ControlRestart.cpp	(revision 4047)
+++ /issm/trunk/src/c/solutions/ControlRestart.cpp	(revision 4048)
@@ -1,5 +1,4 @@
-
 /*!\file: ControlRestart.cpp
- * \brief: core of the control solution 
+ * \brief: save as much as possible, to be able to restart the control_core solution
  */ 
 
@@ -8,39 +7,20 @@
 #include "../EnumDefinitions/EnumDefinitions.h"
 
-void ControlRestart(Model* model,double* param_g){
+void ControlRestart(FemModel* femmodel){
 
-	extern int my_rank;
-
-	/*output: */
-	Results* results=NULL;
 	char*    outputfilename=NULL;
-
-	/*Intermediary: */
-	int      i;
-	int      numberofnodes;
-	double* param_g_copy;
-
-	/*Recover parameters used throughout the solution:*/
-	model->FindParam(&numberofnodes,NumberOfNodesEnum);
+	
+	/*retrieve output file name: */
 	model->FindParam(&outputfilename,OutputFileNameEnum);
 
-	/*Plug COPYS of the results into output dataset: 
-	 * only the pointer is given to temporary_results and at the
-	 * end of ProcessResultsx the pointer is deleted. That would 
-	 * destroy param_g*/
+	/*we essentially want J and the parameter: */
+	InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type); //the parameter itself!
+	femmodel->otherresults->AddObject(new DoubleResult(femmodel->otherresults->Size()+1,0,1,"J",J,nsteps));
+	femmodel->otherresults->AddObject(new StringResult(results->Size()+1,ControlTypeEnum,0,1,EnumAsString(control_type)));
 
-	param_g_copy=(double*)xcalloc(numberofnodes,sizeof(double));
-	for(i=0;i<numberofnodes;i++) param_g_copy[i]=param_g[i];
-
-	results=new Results();
-	results->AddObject(new Result(results->Size()+1,0,1,"param_g",param_g_copy,numberofnodes));
-	results->AddObject(new StringResult(results->Size()+1,AnalysisTypeEnum,0,1,EnumAsString(DiagnosticAnalysisEnum)));
-
-	//Write results on disk
-	OutputResults(results,outputfilename);
+	/*write to disk: */
+	OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters, outputfilename);
 	
 	/*Free ressources:*/
-	delete results;
 	xfree((void**)&outputfilename);
-	xfree((void**)&param_g_copy);
 }
Index: /issm/trunk/src/c/solutions/control_core.cpp
===================================================================
--- /issm/trunk/src/c/solutions/control_core.cpp	(revision 4047)
+++ /issm/trunk/src/c/solutions/control_core.cpp	(revision 4048)
@@ -14,36 +14,30 @@
 void control_core(FemModel* femmodel){
 
-	Vec     u_g=NULL;
-	Vec     t_g=NULL;
-	Vec     m_g=NULL;
-	double  search_scalar;
+	int     i,n;
+	
+	/*parameters: */
+	int     verbose=0;
 	int     control_type;
+	int     control_steady;
+	int     nsteps;
+	double  eps_cm;
+	double  tolx;
+	double  cm_min;
+	double  cm_max;
+	int     cm_gradient;
+
 	double* fit=NULL;
 	double* optscal=NULL;
 	double* maxiter=NULL;
 	double* cm_jump=NULL;
-	double  eps_cm;
-	double  tolx;
-	double* param_g=NULL;
-	Vec     grad_g=NULL;
-	Vec     new_grad_g=NULL;
-	Vec     grad_g_old=NULL;
-	double* grad_g_double=NULL;
-	double  cm_min;
-	double  cm_max;
-	int     cm_gradient;
-	int     nsteps,n,i;
-	double* J=NULL;
+
+		
+	/*intermediary: */
+	double  search_scalar;
 	OptArgs optargs;
 	OptPars optpars;
-	Param*  param=NULL;
 
-	/*flags: */
-	int analysis_type;
-	int sub_analysis_type;
-	int     control_steady;
-	int verbose=0;
-	int converged=0;
-	int numberofnodes;
+	/*output: */
+	double* J=NULL;
 
 	/*Process models*/
@@ -62,22 +56,20 @@
 	femmodel->parameters->FindParam(&cm_max,CmMaxEnum);
 	femmodel->parameters->FindParam(&cm_gradient,CmGradientEnum);
-	femmodel->parameters->FindParam(&param_g,NULL,NULL,ControlParameterEnum);
-	femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
-	femmodel->parameters->FindParam(&sub_analysis_type,SubAnalysisTypeEnum);
-	femmodel->parameters->FindParam(&numberofnodes,NumberOfNodesEnum);
 	femmodel->parameters->FindParam(&control_steady,ControlSteadyEnum);
 
 	/*Initialize misfit: */
 	J=(double*)xmalloc(nsteps*sizeof(double));
-
+		
+	/*Initialize some of the BrentSearch arguments: */
+	optargs.femmodel=femmodel;
+	optpars.xmin=0; optpars.xmax=1; optpars.tolerance=tolx ;
+	
 	/*Start looping: */
 	for(n=0;n<nsteps;n++){
 
 		_printf_("\n%s%i%s%i\n","   control method step ",n+1,"/",nsteps);
-		inputs->Add(control_type,param_g,1,numberofnodes);
 		femmodel->UpdateInputsFromConstant(fit[n],FitEnum);
 		
-		/*In case we are running a steady state control method, compute new temperature field using new parameter 
-		 * distribution: */
+		/*In case we are running a steady state control method, compute new temperature field using new parameter * distribution: */
 		if (control_steady) steadystate_core(model);
 	
@@ -87,109 +79,43 @@
 		/*Return gradient if asked: */
 		if (cm_gradient){
-			/*Transfer gradient from input to results: */
 			InputToResultx(femodel->elements,femodel->nodes,femodel->vertices,femodel->loads,femodel->materials,femodel->parameters,GradientEnum);
-			return;
+			goto cleanup_and_return;
 		}
 
-		/*Normalize if last gradient not satisfying (search_scalar==0)*/
-		if (n>0 && search_scalar==0){
-			_printf_("%s","      orthogonalization...");
-			Orthx(&new_grad_g,grad_g,grad_g_old);
-			_printf_("%s\n"," done.");
-		}
-		else{ 
-			_printf_("%s","      normalizing directions...");
-			Orthx(&new_grad_g,grad_g,NULL);
-			_printf_("%s\n"," done.");
-		}
-		VecFree(&grad_g); VecFree(&grad_g_old); 
-		grad_g_old=new_grad_g;
-		VecToMPISerial(&grad_g_double,new_grad_g);
-
 		_printf_("%s\n","      optimizing along gradient direction");
-		optargs.model=model;
-		optargs.param_g=param_g; optargs.grad_g=grad_g_double; optargs.n=n;
-		optpars.xmin=0; optpars.xmax=1; optpars.tolerance=tolx; optpars.maxiter=(int)maxiter[n];optpars.cm_jump=cm_jump[n];
+		optargs.n=n; optpars.maxiter=(int)maxiter[n]; optpars.cm_jump=cm_jump[n];
 		BrentSearch(&search_scalar,J+n,&optpars,&objectivefunctionC,&optargs);
 
 		_printf_("%s","      updating parameter using optimized search scalar...");
-		for(i=0;i<numberofnodes;i++)param_g[i]=param_g[i]+search_scalar*optscal[n]*grad_g_double[i];
-		_printf_("%s\n"," done.");
+		AXPYInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type,search_scalar*optscal[n],ControlParameterEnum);
 
 		_printf_("%s","      constraining the new distribution...");    
-		ControlConstrainx(param_g,numberofnodes,cm_min,cm_max,control_type);
-		_printf_("%s\n"," done.");
+		ControlConstrainx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type,cm_min,cm_max);
+		
+		_printf_("%s","      save new parameter...");
+		DuplicateInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type,ControlParameterEnum);
 		
 		_printf_("%s%i%s%g\n","      value of misfit J after optimization #",n+1,": ",J[n]);
+		
+		if(controlconvergence)goto cleanup_and_return;
 
-		/*some freeing:*/
-		xfree((void**)&grad_g_double);
-
-		/*Has convergence been reached?*/
-		if (!isnan(eps_cm)){
-			i=n-2;
-			//go through the previous misfits(starting from n-2)
-			while(i>=0){
-				if (fit[i]==fit[n]){
-					//convergence test only if we have the same misfits
-					if ((J[i]-J[n])/J[n] <= eps_cm){
-						//convergence if convergence criteria fullfilled
-						converged=1;
-						_printf_("%s%g%s%g\n","      Convergence criterion: dJ/J = ",(J[i]-J[n])/J[n],"<",eps_cm);
-					}
-					else{
-						_printf_("%s%g%s%g\n","      Convergence criterion: dJ/J = ",(J[i]-J[n])/J[n],">",eps_cm);
-					}
-					break;
-				}
-				i=i-1;
-			}
-		}
-		//stop if convergence has been reached
-		if(converged) break;
-
-		//some temporary saving
+		/*Temporary saving every 5 control steps: */
 		if (((n+1)%5)==0){
 			_printf_("%s","      saving temporary results...");
-			ControlRestart(model,param_g);
-			_printf_("%s\n"," done.");
+			ControlRestart(femmodel);
 		}
 	}
 
-	/*Write results to disk: */
 	_printf_("%s","      preparing final velocity solution");
-	/*Launch diagnostic with the last parameter distribution*/
-	if (control_steady){
-		model->UpdateInputsFromVector(param_g,control_type,VertexEnum);
-		steadystate_results=steadystate_core(model);
+	if (control_steady) steadystate_core(femmodel);
+	else diagnostic_core(femmodel);
 
-		//extract u_g ,t_g and m_g from steadystate results, and erase diagnostic_results;
-		steadystate_results->FindResult(&u_g,"u_g");
-		steadystate_results->FindResult(&m_g,"m_g");
-		steadystate_results->FindResult(&t_g,"t_g");
-		delete steadystate_results;
-	}
-	else{
-		model->UpdateInputsFromVector(param_g,control_type,VertexEnum);
-		diagnostic_results=diagnostic_core(model);
+	/*some results not computed by steadystate_core or diagnostic_core: */
+	InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type); //the parameter itself!
+	femmodel->otherresults->AddObject(new DoubleResult(femmodel->otherresults->Size()+1,0,1,"J",J,nsteps));
+	femmodel->otherresults->AddObject(new StringResult(results->Size()+1,ControlTypeEnum,0,1,EnumAsString(control_type)));
 
-		//extract u_g from diagnostic_results, and erase diagnostic_results;
-		diagnostic_results->FindResult(&u_g,"u_g");
-		delete diagnostic_results;
-	}
-
-	/*Plug results into output dataset: */
-	results->AddObject(new Result(results->Size()+1,0,1,"u_g",u_g));
-	results->AddObject(new Result(results->Size()+1,0,1,"param_g",param_g,numberofnodes));
-	results->AddObject(new Result(results->Size()+1,0,1,"J",J,nsteps));
-	if (control_steady){
-		results->AddObject(new Result(results->Size()+1,0,1,"t_g",t_g));
-		results->AddObject(new Result(results->Size()+1,0,1,"m_g",m_g));
-	}
+	cleanup_and_return:
 	
-	/*Add analysis_type and control_type to results: */
-	results->AddObject(new StringResult(results->Size()+1,AnalysisTypeEnum,0,1,EnumAsString(analysis_type)));
-	results->AddObject(new StringResult(results->Size()+1,ControlTypeEnum,0,1,EnumAsString(control_type)));
-
 	/*Free ressources: */
 	xfree((void**)&control_type);
@@ -198,13 +124,4 @@
 	xfree((void**)&maxiter);
 	xfree((void**)&cm_jump);
-	VecFree(&new_grad_g); //do not VecFree grad_g and grad_g_old, they point to new_grad_g
-	xfree((void**)&grad_g_double);
-	xfree((void**)&param_g);
-	VecFree(&u_g);
-	VecFree(&t_g);
-	VecFree(&m_g);
 	xfree((void**)&J);
-
-	//return: 
-	return results;
 }
Index: /issm/trunk/src/c/solutions/controlconvergence.cpp
===================================================================
--- /issm/trunk/src/c/solutions/controlconvergence.cpp	(revision 4048)
+++ /issm/trunk/src/c/solutions/controlconvergence.cpp	(revision 4048)
@@ -0,0 +1,41 @@
+/*!\file: controlconvergence.cpp
+ * \brief: determine convergence of control_core solution
+ */ 
+
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+int controlconvergence(double* J, double* fit, double eps_cm, int n){
+
+	int i;
+	int converged=0;
+
+	/*Has convergence been reached?*/
+	if (!isnan(eps_cm)){
+		i=n-2;
+
+		//go through the previous misfits(starting from n-2)
+		while(i>=0){
+			if (fit[i]==fit[n]){
+				//convergence test only if we have the same misfits
+				if ((J[i]-J[n])/J[n] <= eps_cm){
+					//convergence if convergence criteria fullfilled
+					converged=1;
+					_printf_("%s%g%s%g\n","      Convergence criterion: dJ/J = ",(J[i]-J[n])/J[n],"<",eps_cm);
+				}
+				else{
+					_printf_("%s%g%s%g\n","      Convergence criterion: dJ/J = ",(J[i]-J[n])/J[n],">",eps_cm);
+				}
+				break;
+			}
+			i=i-1;
+		}
+	}
+
+
+	return converged;
+}
+
Index: /issm/trunk/src/c/solutions/gradient_core.cpp
===================================================================
--- /issm/trunk/src/c/solutions/gradient_core.cpp	(revision 4047)
+++ /issm/trunk/src/c/solutions/gradient_core.cpp	(revision 4048)
@@ -13,6 +13,5 @@
 
 
-void gradient_core(FemModel* femmodel){
-	
+void gradient_core(FemModel* femmodel,int step, double search_scalar){ //step defaults to 0, search_scalar defaults to 0
 	
 	/*parameters: */
@@ -20,4 +19,7 @@
 	int control_type;
 	int verbose;
+	Vec new_gradient=NULL;
+	Vec gradient=NULL;
+	Vec old_gradient=NULL;
 
 	/*retrieve parameters:*/
@@ -30,7 +32,31 @@
 	
 	if(verbose)_printf_("%s\n","      compute gradient:");
-	Gradjx( femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters, control_type);
+	Gradjx(&gradient, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters, control_type,AdjointEnum);
+	
+	if(verbose)_printf_("%s\n","      retrieve old gradient:");
+	GetVectorFromInputsx(&old_gradient, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters, OldGradientEnum,VertexEnum);
 
 	if(control_steady)diagnostic_core(model);
+	
+	if (step>0 && search_scalar==0){
+		_printf_("%s","      orthogonalization...");
+		Orthx(&new_gradient,gradient,old_gradient);
+	}
+	else{ 
+		_printf_("%s","      normalizing directions...");
+		Orthx(&new_gradient,gradient,NULL);
+	}
+	_printf_("%s\n"," done.");
 
+	
+	/*point gradient and old_gradient to new_gradient: */
+	VecFree(&gradient); gradient=new_gradient;
+	VecFree(&old_gradient); old_gradient=new_gradient;
+
+	/*plug back into inputs: */
+	UpdateInputsFromVectorx( femmodel-> elements,femmodel-> nodes, femmodel-> vertices,femmodel-> loads, femmodel-> materials,  femmoel->parameters,gradient,GradientEnum,VertexEnum);
+	UpdateInputsFromVectorx( femmodel-> elements,femmodel-> nodes, femmodel-> vertices,femmodel-> loads, femmodel-> materials,  femmoel->parameters,old_gradient,OldGradientEnum,VertexEnum);
+
+	/*Free ressources and return:*/
+	VecFree(&new_gradient);
 }
Index: /issm/trunk/src/c/solutions/objectivefunctionC.cpp
===================================================================
--- /issm/trunk/src/c/solutions/objectivefunctionC.cpp	(revision 4047)
+++ /issm/trunk/src/c/solutions/objectivefunctionC.cpp	(revision 4048)
@@ -21,48 +21,21 @@
 	
 	/*parameters: */
-	Model*    model=NULL;
 	FemModel* femmodel=NULL;
-	Results* diagnostic_results=NULL;
-	double* param_g=NULL;
-	double* grad_g=NULL;
-	int numberofdofspernode;
 	int n;
-
-	/*intermediary:*/
-	int gsize;
 	double* optscal=NULL;
 	double* fit=NULL;
 	double  cm_min;
 	double  cm_max;
-	int   control_type;
-	double* param_g_copy=NULL;
+	int     control_type;
 	int     analysis_type;
 	int     sub_analysis_type;
 	int     control_steady;
-	Vec     u_g=NULL;
-	Vec     u_g_full=NULL;
-	double* u_g_double=NULL;
-	double* vx=NULL;
-	double* vy=NULL;
-	double* vz=NULL;
-	int     numberofnodes;
 		
-	/*steadystate: */
-	int     dt=0;
-	int     isstokes=0;
-	Results* results_steadystate=NULL;
-	int dofs01[2]={0,1};
-	double* dofset=NULL;
-
-	/*Recover active model: */
-	model=optargs->model;
-	femmodel=model->GetActiveFormulation();
+	/*Recover finite element model: */
+	femmodel=optargs->femmodel;
 
 	/*Recover parameters: */
-	param_g=optargs->param_g;
-	grad_g=optargs->grad_g;
 	n=optargs->n;
 
-	gsize=femmodel->nodesets->GetGSize();
 	femmodel->parameters->FindParam(&optscal,NULL,NULL,OptScalEnum);
 	femmodel->parameters->FindParam(&fit,NULL,NULL,FitEnum);
@@ -73,41 +46,17 @@
 	femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
 	femmodel->parameters->FindParam(&sub_analysis_type,SubAnalysisTypeEnum);
-	femmodel->parameters->FindParam(&isstokes,IsStokesEnum);
-	femmodel->parameters->FindParam(&numberofnodes,NumberOfNodesEnum);
-	femmodel->parameters->FindParam(&numberofdofspernode,NumberOfDofsPerNodeEnum);
 
-	/*First copy param_g so we don't modify it: */
-	param_g_copy=(double*)xmalloc(numberofnodes*sizeof(double));
-	memcpy(param_g_copy,param_g,numberofnodes*sizeof(double));
-
-	/*First, update param_g using search_scalar: */
-	for(i=0;i<numberofnodes;i++)param_g_copy[i]=param_g_copy[i]+search_scalar*optscal[n]*grad_g[i];
+	/*Use ControlParameterEnum input to  reinitialize our input parameter: */
+	DuplicateInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ControlParameterEnum,control_type);
+	
+	/*Use search scalar to shoot parameter in the gradient direction: */
+	AXPYInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type,search_scalar*optscal[n],ControlParameterEnum);
 
 	/*Constrain:*/
-	ControlConstrainx(param_g_copy,numberofnodes,cm_min,cm_max,control_type);
+	ControlConstrainx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type,cm_min,cm_max);
 
-	/*Add new parameter to inputs: */
-	UpdateInputsFromVectorx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,param_g_copy,control_type,VertexEnum);
-
-	/*Run diagnostic with updated parameters.*/
-	if(!control_steady){
-		diagnostic_core_nonlinear(&u_g,NULL,NULL,NULL,femmodel,DiagnosticAnalysisEnum,sub_analysis_type); 
-		UpdateInputsFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,u_g,DiagnosticAnalysisEnum,sub_analysis_type);
-		VecFree(&u_g);
-	}
-	else{
-		//We need a 3D velocity!! (vz is required for the next thermal run)
-		diagnostic_results=	diagnostic_core(model);
-
-		//extract u_g and add it to input (3d velocity needed by thermal_core)
-		diagnostic_results->FindResult(&u_g,"u_g");
-		
-		SplitSolutionVectorx(u_g,numberofnodes,3,&vx,&vy,&vz);
-		UpdateInputsFromVectorx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,vx,VxEnum,VertexEnum);
-		UpdateInputsFromVectorx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,vy,VyEnum,VertexEnum);
-		UpdateInputsFromVectorx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,vz,VzEnum,VertexEnum);
-		
-		delete diagnostic_results;
-	}
+	/*Run diagnostic with updated inputs: */
+	if(!control_steady) solver_diagnostic_nonlinear(NULL,NULL,NULL,true,femmodel,DiagnosticAnalysisEnum,sub_analysis_type);  //true means we conserve loads at each diagnostic run
+	else                diagnostic_core(femmodel);	//We need a 3D velocity!! (vz is required for the next thermal run)
 
 	/*Compute misfit for this velocity field.*/
@@ -118,10 +67,4 @@
 	xfree((void**)&fit);
 	xfree((void**)&optscal);
-	xfree((void**)&param_g_copy);
-	xfree((void**)&dofset);
-	xfree((void**)&vx);
-	xfree((void**)&vy);
-	xfree((void**)&vz);
-
 	return J;
 }
Index: /issm/trunk/src/c/solutions/solutions.h
===================================================================
--- /issm/trunk/src/c/solutions/solutions.h	(revision 4047)
+++ /issm/trunk/src/c/solutions/solutions.h	(revision 4048)
@@ -25,5 +25,5 @@
 Results* transient_core_3d(FemModel* model);
 void adjoint_core(FemModel* model);
-void gradient_core(FemModel* model);
+void gradient_core(FemModel* model,int step=0, double search_scalar=0);
 void diagnostic_core(FemModel* model);
 void thermal_core(FemModel* model);
@@ -33,4 +33,5 @@
 //int GradJOrth(WorkspaceParams* workspaceparams);
 void convergence(int* pconverged, Mat K_ff,Vec p_f,Vec u_f,Vec u_f_old,Parameters* parameters);
+int controlconvergence(double* J, double* fit, double eps_cm, int n);
 
 int GoldenSearch(double* psearch_scalar,double* pJ,double xa, double xb, double tolerance, int maxiter, double fit,double optscal,double (*f)(double*,double,double,FemModel*),FemModel* femmodel);
@@ -47,5 +48,5 @@
 
 void ControlInitialization(FemModel* model);
-void ControlRestart(FemModel* model,double* param_g);
+void ControlRestart(FemModel* femmodel);
 
 void CreateFemModel(FemModel* femmodel,ConstDataHandle MODEL,int analysis_type,int sub_analysis_type);