Index: issm/trunk/src/c/DataSet/DataSet.cpp
===================================================================
--- issm/trunk/src/c/DataSet/DataSet.cpp	(revision 4054)
+++ issm/trunk/src/c/DataSet/DataSet.cpp	(revision 4055)
@@ -684,29 +684,4 @@
 	xfree((void**)&truedofs);
 	xfree((void**)&alltruedofs);
-
-}
-/*}}}*/
-/*FUNCTION DataSet::FieldExtrude{{{1*/
-void  DataSet::FieldExtrude(Vec field,double* field_serial,char* field_name, int collapse){
-
-	vector<Object*>::iterator object;
-	Penta* penta=NULL;
-	Node*  node=NULL;
-
-	for ( object=objects.begin() ; object < objects.end(); object++ ){
-
-		if((*object)->Enum()==PentaEnum){
-
-			penta=(Penta*)(*object);
-			penta->FieldExtrude(field,field_serial,field_name,collapse);
-
-		}
-		if((*object)->Enum()==NodeEnum){
-
-			node=(Node*)(*object);
-			node->FieldExtrude(field,field_serial,field_name);
-
-		}
-	}
 
 }
Index: issm/trunk/src/c/DataSet/DataSet.h
===================================================================
--- issm/trunk/src/c/DataSet/DataSet.h	(revision 4054)
+++ issm/trunk/src/c/DataSet/DataSet.h	(revision 4055)
@@ -80,5 +80,4 @@
 		void  InputExtrude(int enum_type);
 		int   DeleteObject(Object* object);
-		void  FieldExtrude(Vec field,double* field_serial,char* field_name, int collapse);
 		void  OutputRifts(Vec riftproperties);
 		/*}}}*/
Index: issm/trunk/src/c/EnumDefinitions/EnumAsString.cpp
===================================================================
--- issm/trunk/src/c/EnumDefinitions/EnumAsString.cpp	(revision 4054)
+++ issm/trunk/src/c/EnumDefinitions/EnumAsString.cpp	(revision 4055)
@@ -26,7 +26,4 @@
 		case VerticesEnum : return "Vertices";
 		case SolutionTypeEnum : return "SolutionType";
-		case DiagnosticSolutionEnum : return "DiagnosticSolution";
-		case SurfaceSlopeSolutionEnum : return "SurfaceSlopeSolution";
-		case BedSlopeSolutionEnum : return "BedSlopeSolution";
 		case AnalysisTypeEnum : return "AnalysisType";
 		case SubAnalysisTypeEnum : return "SubAnalysisType";
@@ -46,7 +43,8 @@
 		case InverseAnalysisEnum : return "InverseAnalysis";
 		case ThermalAnalysisEnum : return "ThermalAnalysis";
-		case TransientAnalysisEnum : return "TransientAnalysis";
+		case Transient2DAnalysisEnum : return "Transient2DAnalysis";
+		case Transient3DAnalysisEnum : return "Transient3DAnalysis";
 		case SteadyAnalysisEnum : return "SteadyAnalysis";
-		case SlopeComputeAnalysisEnum : return "SlopeComputeAnalysis";
+		case SlopeAnalysisEnum : return "SlopeAnalysis";
 		case BedSlopeXAnalysisEnum : return "BedSlopeXAnalysis";
 		case BedSlopeYAnalysisEnum : return "BedSlopeYAnalysis";
@@ -220,4 +218,7 @@
 		case StringExternalResultEnum : return "StringExternalResult";
 		case JEnum : return "J";
+		case RelativeEnum : return "Relative";
+		case ResidualEnum : return "Residual";
+		case AbsoluteEnum : return "Absolute";
 		case BetaEnum : return "Beta";
 		case CmGradientEnum : return "CmGradient";
@@ -225,5 +226,7 @@
 		case CmMaxEnum : return "CmMax";
 		case CmMinEnum : return "CmMin";
+		case AdjointEnum : return "Adjoint";
 		case GradientEnum : return "Gradient";
+		case OldGradientEnum : return "OldGradient";
 		case ConnectivityEnum : return "Connectivity";
 		case ControlParameterEnum : return "ControlParameter";
Index: issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h
===================================================================
--- issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 4054)
+++ issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h	(revision 4055)
@@ -22,7 +22,4 @@
 	/*Solution types {{{1 */
 	SolutionTypeEnum,
-	DiagnosticSolutionEnum,
-	SurfaceSlopeSolutionEnum,
-	BedSlopeSolutionEnum,
 	/*}}}*/
 	/*Analysis types {{{1 */
@@ -48,8 +45,9 @@
 	ThermalAnalysisEnum,
 	//transient
-	TransientAnalysisEnum,
+	Transient2DAnalysisEnum,
+	Transient3DAnalysisEnum,
 	SteadyAnalysisEnum,
 	//slope
-	SlopeComputeAnalysisEnum,
+	SlopeAnalysisEnum,
 	BedSlopeXAnalysisEnum,
 	BedSlopeYAnalysisEnum,
@@ -254,4 +252,9 @@
 	JEnum,
 	/*}}}*/
+	/*Convergence{{{1*/
+	RelativeEnum,
+	ResidualEnum,
+	AbsoluteEnum,
+	/*}}}*/
 	/*Parameters{{{1*/
 	BetaEnum,
@@ -260,5 +263,7 @@
 	CmMaxEnum,
 	CmMinEnum,
+	AdjointEnum,
 	GradientEnum,
+	OldGradientEnum,
 	ConnectivityEnum,
 	ControlParameterEnum,
Index: issm/trunk/src/c/EnumDefinitions/StringAsEnum.cpp
===================================================================
--- issm/trunk/src/c/EnumDefinitions/StringAsEnum.cpp	(revision 4054)
+++ issm/trunk/src/c/EnumDefinitions/StringAsEnum.cpp	(revision 4055)
@@ -24,7 +24,4 @@
 	else if (strcmp(name,"Vertices")==0) return VerticesEnum;
 	else if (strcmp(name,"SolutionType")==0) return SolutionTypeEnum;
-	else if (strcmp(name,"DiagnosticSolution")==0) return DiagnosticSolutionEnum;
-	else if (strcmp(name,"SurfaceSlopeSolution")==0) return SurfaceSlopeSolutionEnum;
-	else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum;
 	else if (strcmp(name,"AnalysisType")==0) return AnalysisTypeEnum;
 	else if (strcmp(name,"SubAnalysisType")==0) return SubAnalysisTypeEnum;
@@ -44,7 +41,8 @@
 	else if (strcmp(name,"InverseAnalysis")==0) return InverseAnalysisEnum;
 	else if (strcmp(name,"ThermalAnalysis")==0) return ThermalAnalysisEnum;
-	else if (strcmp(name,"TransientAnalysis")==0) return TransientAnalysisEnum;
+	else if (strcmp(name,"Transient2DAnalysis")==0) return Transient2DAnalysisEnum;
+	else if (strcmp(name,"Transient3DAnalysis")==0) return Transient3DAnalysisEnum;
 	else if (strcmp(name,"SteadyAnalysis")==0) return SteadyAnalysisEnum;
-	else if (strcmp(name,"SlopeComputeAnalysis")==0) return SlopeComputeAnalysisEnum;
+	else if (strcmp(name,"SlopeAnalysis")==0) return SlopeAnalysisEnum;
 	else if (strcmp(name,"BedSlopeXAnalysis")==0) return BedSlopeXAnalysisEnum;
 	else if (strcmp(name,"BedSlopeYAnalysis")==0) return BedSlopeYAnalysisEnum;
@@ -218,4 +216,7 @@
 	else if (strcmp(name,"StringExternalResult")==0) return StringExternalResultEnum;
 	else if (strcmp(name,"J")==0) return JEnum;
+	else if (strcmp(name,"Relative")==0) return RelativeEnum;
+	else if (strcmp(name,"Residual")==0) return ResidualEnum;
+	else if (strcmp(name,"Absolute")==0) return AbsoluteEnum;
 	else if (strcmp(name,"Beta")==0) return BetaEnum;
 	else if (strcmp(name,"CmGradient")==0) return CmGradientEnum;
@@ -223,5 +224,7 @@
 	else if (strcmp(name,"CmMax")==0) return CmMaxEnum;
 	else if (strcmp(name,"CmMin")==0) return CmMinEnum;
+	else if (strcmp(name,"Adjoint")==0) return AdjointEnum;
 	else if (strcmp(name,"Gradient")==0) return GradientEnum;
+	else if (strcmp(name,"OldGradient")==0) return OldGradientEnum;
 	else if (strcmp(name,"Connectivity")==0) return ConnectivityEnum;
 	else if (strcmp(name,"ControlParameter")==0) return ControlParameterEnum;
Index: issm/trunk/src/c/Makefile.am
===================================================================
--- issm/trunk/src/c/Makefile.am	(revision 4054)
+++ issm/trunk/src/c/Makefile.am	(revision 4055)
@@ -510,8 +510,10 @@
 					./modules/DepthAverageInputx/DepthAverageInputx.cpp\
 					./modules/DepthAverageInputx/DepthAverageInputx.h\
-					./modules/FieldExtrudex/FieldExtrudex.cpp\
-					./modules/FieldExtrudex/FieldExtrudex.h\
+					./modules/VecExtrudex/VecExtrudex.cpp\
+					./modules/VecExtrudex/VecExtrudex.h\
 					./modules/InputToResultx/InputToResultx.cpp\
 					./modules/InputToResultx/InputToResultx.h\
+					./modules/InputConvergencex/InputConvergencex.cpp\
+					./modules/InputConvergencex/InputConvergencex.h\
 					./modules/ExtrudeInputx/ExtrudeInputx.cpp\
 					./modules/ExtrudeInputx/ExtrudeInputx.h\
@@ -1014,8 +1016,10 @@
 					./modules/DepthAverageInputx/DepthAverageInputx.cpp\
 					./modules/DepthAverageInputx/DepthAverageInputx.h\
-					./modules/FieldExtrudex/FieldExtrudex.cpp\
-					./modules/FieldExtrudex/FieldExtrudex.h\
+					./modules/VecExtrudex/VecExtrudex.cpp\
+					./modules/VecExtrudex/VecExtrudex.h\
 					./modules/InputToResultx/InputToResultx.cpp\
 					./modules/InputToResultx/InputToResultx.h\
+					./modules/InputConvergencex/InputConvergencex.cpp\
+					./modules/InputConvergencex/InputConvergencex.h\
 					./modules/ExtrudeInputx/ExtrudeInputx.cpp\
 					./modules/ExtrudeInputx/ExtrudeInputx.h\
@@ -1030,5 +1034,6 @@
 					./solutions/controlrestart.cpp\
 					./solutions/objectivefunctionC.cpp\
-					./solutions/gradjcompute_core.cpp\
+					./solutions/gradient_core.cpp\
+					./solutions/adjoint_core.cpp\
 					./solutions/prognostic_core.cpp\
 					./solutions/prognostic2_core.cpp\
@@ -1040,7 +1045,6 @@
 					./solutions/bedslope.cpp\
 					./solutions/bedslope_core.cpp\
-					./solutions/transient_core.cpp\
-					./solutions/transient_core_2d.cpp\
-					./solutions/transient_core_3d.cpp\
+					./solutions/transient2d_core.cpp\
+					./solutions/transient3d_core.cpp\
 					./solutions/steadystate_core.cpp\
 					./solutions/ResetBoundaryConditions.cpp\
@@ -1058,5 +1062,5 @@
 bin_PROGRAMS = 
 else 
-bin_PROGRAMS = DiagnosticSolution.exe ThermalSolution.exe PrognosticSolution.exe Prognostic2Solution.exe BalancedthicknessSolution.exe Balancedthickness2Solution.exe BalancedvelocitiesSolution.exe TransientSolution.exe SteadystateSolution.exe SurfaceSlopeSolution.exe BedSlopeSolution.exe 
+bin_PROGRAMS = DiagnosticSolution.exe ThermalSolution.exe PrognosticSolution.exe Prognostic2Solution.exe BalancedthicknessSolution.exe Balancedthickness2Solution.exe BalancedvelocitiesSolution.exe Transient2DSolution.exe Transient23Solution.exe SteadystateSolution.exe SurfaceSlopeSolution.exe BedSlopeSolution.exe 
 endif
 
@@ -1093,4 +1097,7 @@
 BedSlopeSolution_exe_CXXFLAGS= -fPIC -D_PARALLEL_ 
 
-TransientSolution_exe_SOURCES = solutions/transient.cpp
-TransientSolution_exe_CXXFLAGS= -fPIC -D_PARALLEL_ 
+Transient2DSolution_exe_SOURCES = solutions/transient2d.cpp
+Transient2DSolution_exe_CXXFLAGS= -fPIC -D_PARALLEL_ 
+
+Transient3DSolution_exe_SOURCES = solutions/transient3d.cpp
+Transient3DSolution_exe_CXXFLAGS= -fPIC -D_PARALLEL_ 
Index: issm/trunk/src/c/modules/CostFunctionx/CostFunctionx.h
===================================================================
--- issm/trunk/src/c/modules/CostFunctionx/CostFunctionx.h	(revision 4054)
+++ issm/trunk/src/c/modules/CostFunctionx/CostFunctionx.h	(revision 4055)
@@ -10,6 +10,5 @@
 
 /* local prototypes: */
-void CostFunctionx( double* pJ, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters, 
-			int analysis_type,int sub_analysis_type);
+void CostFunctionx( double* pJ, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters);
 
 #endif  /* _MISFITX_H */
Index: issm/trunk/src/c/modules/DepthAverageInputx/DepthAverageInputx.cpp
===================================================================
--- issm/trunk/src/c/modules/DepthAverageInputx/DepthAverageInputx.cpp	(revision 4054)
+++ issm/trunk/src/c/modules/DepthAverageInputx/DepthAverageInputx.cpp	(revision 4055)
@@ -10,5 +10,5 @@
 #include "../../EnumDefinitions/EnumDefinitions.h"
 
-void DepthAverageInputx( DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters,int enum_type){
+void DepthAverageInputx( DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters,int enum_type,int average_enum_type){
 	
 	/*Intermediary*/
@@ -22,9 +22,8 @@
 	for (i=0;i<elements->Size();i++){
 		element=(Element*)elements->GetObjectByOffset(i);
-		element->DepthAverageInputAtBase(enum_type);
+		element->DepthAverageInputAtBase(enum_type,average_enum_type);
 	}
 
 	/*Then extrude vertically the new inputs*/
-	elements->InputExtrude(enum_type);
-
+	elements->InputExtrude(average_enum_type);
 }
Index: issm/trunk/src/c/modules/DepthAverageInputx/DepthAverageInputx.h
===================================================================
--- issm/trunk/src/c/modules/DepthAverageInputx/DepthAverageInputx.h	(revision 4054)
+++ issm/trunk/src/c/modules/DepthAverageInputx/DepthAverageInputx.h	(revision 4055)
@@ -9,5 +9,5 @@
 
 /* local prototypes: */
-void DepthAverageInputx( DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters,int enum_type);
+void DepthAverageInputx( DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters,int enum_type,int average_enum_type);
 
 #endif  /* _DEPTHAVERAGEINPUTX_H */
Index: issm/trunk/src/c/modules/Dux/Dux.h
===================================================================
--- issm/trunk/src/c/modules/Dux/Dux.h	(revision 4054)
+++ issm/trunk/src/c/modules/Dux/Dux.h	(revision 4055)
@@ -10,7 +10,5 @@
 
 /* local prototypes: */
-void Dux( Vec* pdu_g, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters, 
-			int analysis_type,int sub_analysis_type);
-
+void Dux( Vec* pdu_g, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters);
 #endif  /* _DUX_H */
 
Index: issm/trunk/src/c/modules/FieldDepthAveragex/FieldDepthAveragex.cpp
===================================================================
--- issm/trunk/src/c/modules/FieldDepthAveragex/FieldDepthAveragex.cpp	(revision 4054)
+++ issm/trunk/src/c/modules/FieldDepthAveragex/FieldDepthAveragex.cpp	(revision 4055)
@@ -43,5 +43,8 @@
 
 	/*Extrude field to all nodes: */
-	nodes->FieldExtrude(field,field_serial,fieldname,0);
+	for(i=0;i<nodes->Size();i++){
+		node=(Node*)nodes->GetObjectByOffset(i);
+		node->VecExtrude(field,field_serial);
+	}
 
 	/*Assemble vector: */
Index: issm/trunk/src/c/modules/GetSolutionFromInputsx/GetSolutionFromInputsx.h
===================================================================
--- issm/trunk/src/c/modules/GetSolutionFromInputsx/GetSolutionFromInputsx.h	(revision 4054)
+++ issm/trunk/src/c/modules/GetSolutionFromInputsx/GetSolutionFromInputsx.h	(revision 4055)
@@ -10,5 +10,5 @@
 
 /* local prototypes: */
-void GetSolutionFromInputsx( Vec* psolution, DataSet* elements,DataSet* nodes, DataSet* vertices,DataSet* loads, DataSet* materials,  Parameters* parameters, int analysis_type, int sub_analysis_type);
+void GetSolutionFromInputsx( Vec* psolution, DataSet* elements,DataSet* nodes, DataSet* vertices,DataSet* loads, DataSet* materials,  Parameters* parameters);
 
 #endif  /* _GETSOLUTIONFROMINPUTSXX_H */
Index: issm/trunk/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp
===================================================================
--- issm/trunk/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp	(revision 4054)
+++ issm/trunk/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.cpp	(revision 4055)
@@ -38,2 +38,20 @@
 
 }
+
+void GetVectorFromInputsx( double* pvector, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters, int NameEnum, int TypeEnum){
+	
+	/*output: */
+	double* vector=NULL;
+	
+	/*intermediary: */
+	Vec vec_vector=NULL;
+
+	GetVectorFromInputsx( &vec_vector, elements,nodes, vertices, loads, materials, parameters, NameEnum, TypeEnum);
+	VecToMPISerial(&vector,vec_vector);
+
+	/*Free ressources:*/
+	VecFree(&vec_vector);
+
+	/*Assign output pointers:*/
+	*pvector=vector;
+}
Index: issm/trunk/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.h
===================================================================
--- issm/trunk/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.h	(revision 4054)
+++ issm/trunk/src/c/modules/GetVectorFromInputsx/GetVectorFromInputsx.h	(revision 4055)
@@ -10,4 +10,5 @@
 /* local prototypes: */
 void	GetVectorFromInputsx( Vec* pvector, DataSet* elements,DataSet* nodes, DataSet* vertices,DataSet* loads, DataSet* materials,  Parameters* parameters,int NameEnum,int TypeEnum);
+void	GetVectorFromInputsx( double* 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 4054)
+++ issm/trunk/src/c/modules/Gradjx/Gradjx.cpp	(revision 4055)
@@ -10,10 +10,11 @@
 #include "../../EnumDefinitions/EnumDefinitions.h"
 
-void Gradjx( DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters,int control_type){
+void Gradjx( Vec* pgradient, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters,int control_type){
 
 	int i;
 	int dim;
-
-
+	int numberofvertices;
+	Vec gradient=NULL;
+	
 	/*First, get elements and loads configured: */
 	elements->  Configure(elements,loads, nodes,vertices, materials,parameters);
@@ -23,13 +24,20 @@
 	/*retrieve some parameters: */
 	parameters->FindParam(&dim,DimEnum);
+	parameters->FindParam(&numberofvertices,NumberOfVerticesEnum);
 
 	/*Compute gradients: */
 	for (i=0;i<elements->Size();i++){
 		Element* element=(Element*)elements->GetObjectByOffset(i);
-		element->Gradj(control_type);
+		element->Gradj(gradient,control_type);
 	}
 
+	/*Assemble vector: */
+	VecAssemblyBegin(gradient);
+	VecAssemblyEnd(gradient);
+
 	/*Extrude if needed: */
-	if(dim==3) ExtrudeInputx( elements,nodes, vertices, loads, materials, parameters,GradientEnum);
+	if(dim==3) VecExtrudex(gradient, elements,nodes, vertices, loads, materials, parameters,0);
 
+	/*Assign output pointers: */
+	*pgradient=gradient;
 }
Index: issm/trunk/src/c/modules/Gradjx/Gradjx.h
===================================================================
--- issm/trunk/src/c/modules/Gradjx/Gradjx.h	(revision 4054)
+++ issm/trunk/src/c/modules/Gradjx/Gradjx.h	(revision 4055)
@@ -10,5 +10,5 @@
 
 /* local prototypes: */
-void Gradjx( DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials,  Parameters* parameters, int control_type);
+void Gradjx(Vec* pgrad_g, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials,  Parameters* parameters, int control_type);
 
 #endif  /* _GRADJX_H */
Index: issm/trunk/src/c/modules/InputConvergencex/InputConvergencex.cpp
===================================================================
--- issm/trunk/src/c/modules/InputConvergencex/InputConvergencex.cpp	(revision 4055)
+++ issm/trunk/src/c/modules/InputConvergencex/InputConvergencex.cpp	(revision 4055)
@@ -0,0 +1,58 @@
+/*!\file InputConvergencex
+ */
+
+#include "./InputConvergencex.h"
+#include "../../shared/shared.h"
+#include "../../include/include.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+
+int InputConvergencex(DataSet* elements,DataSet* nodes,DataSet* vertices,DataSet* loads,DataSet* materials,Parameters* parameters,int enums, int num_enums, int criterionenums, double* criterionvalues,int num_criterionenums){
+
+	/*intermediary:*/
+	int         i;
+	int         converged = 1;
+	int         node_converged=1;
+	double     *eps       = NULL;
+	extern int  my_rank;
+	int         rank_allowed;
+
+	/*First, configure our datasets: */
+	elements->Configure(elements,loads, nodes,vertices, materials,parameters);
+	nodes->     Configure(elements,loads, nodes,vertices, materials,parameters);
+	loads->     Configure(elements, loads, nodes,vertices, materials,parameters);
+	parameters->Configure(elements,loads, nodes,vertices, materials,parameters);
+
+	/*allocate dynamic memory: */
+	eps=(double*)xmalloc(num_criterionenums*sizeof(double));
+
+	/*Go through elements, and ask them to do the job: */
+	for(i=0;i<elements->Size();i++){
+		Element* element=(Element*)elements->GetObjectByOffset(i);
+		element->InputConvergence(&converged,eps,enums,num_enums,criterionenums,criterionvalues,num_criterionenums);
+		if(!converged)break;
+	}
+
+	if(!converged){ 
+		/*get convergence information to user: */
+		if(verbose>2)rank_allowed=my_rank;
+		else rank_allowed=0;
+		if(my_rank==rank_allowed){
+			printf("%s%i%s"," non convergence on element ",element->Id()," for inputs ("); for(i=0;i<num_enums-1;i++)printf("%s,",EnumAsString(enums[i])); printf("%s) ",EnumAsString(enums[num_enums-1]));
+			printf("criterions: "); for(i=0;i<num_criterionenums;i++)printf("%s: %g/%g ",EnumAsString(criterionenums[i]),eps[i],criterionvalues[i]);
+		}
+	}
+	
+	/*In parallel, we need to gather the converged status: */
+	#ifdef _PARALLEL_
+	MPI_Reduce (&converged,&node_converged,1,MPI_INT,MPI_MIN,0,MPI_COMM_WORLD );
+	MPI_Bcast(&node_converged,1,MPI_INT,0,MPI_COMM_WORLD);   
+	converged=node_converged;
+	#endif
+
+	/*Free ressources:*/
+	xfree((void**)&eps);
+
+	/*return: */
+	return converged;
+}
Index: issm/trunk/src/c/modules/InputConvergencex/InputConvergencex.h
===================================================================
--- issm/trunk/src/c/modules/InputConvergencex/InputConvergencex.h	(revision 4055)
+++ issm/trunk/src/c/modules/InputConvergencex/InputConvergencex.h	(revision 4055)
@@ -0,0 +1,13 @@
+/*!\file:  InputConvergencex.h
+ */ 
+
+#ifndef _INPUTCONVERGENCEX_H
+#define _INPUTCONVERGENCEX_H
+
+#include "../../DataSet/DataSet.h"
+
+/* local prototypes: */
+int InputConvergencex(DataSet* elements,DataSet* nodes,DataSet* vertices,DataSet* loads,DataSet* materials,Parameters* parameters,int enums, int num_enums, int criterionenums, double* criterionvalues,int num_criterionenums);
+
+#endif  /* _INPUTCONVERGENCEX_H */
+
Index: issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp
===================================================================
--- issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp	(revision 4054)
+++ issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp	(revision 4055)
@@ -57,5 +57,5 @@
 			break;
 
-		case SlopeComputeAnalysisEnum:
+		case SlopeAnalysisEnum:
 			CreateNodesSlopeCompute(pnodes, iomodel,iomodel_handle);
 			CreateConstraintsSlopeCompute(pconstraints,iomodel,iomodel_handle);
Index: issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp
===================================================================
--- issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp	(revision 4054)
+++ issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp	(revision 4055)
@@ -37,5 +37,4 @@
 	if (iomodel->dim==2) parameters->AddObject(new IntParam(DimEnum,2));
 	else parameters->AddObject(new IntParam(DimEnum,3));
-	parameters->AddObject(new StringParam(OutputFileNameEnum,iomodel->outputfilename));
 	parameters->AddObject(new BoolParam(IsHutterEnum,iomodel->ishutter));
 	parameters->AddObject(new BoolParam(IsMacAyealPattynEnum,iomodel->ismacayealpattyn));
Index: issm/trunk/src/c/modules/ModelProcessorx/SlopeCompute/CreateNodesSlopeCompute.cpp
===================================================================
--- issm/trunk/src/c/modules/ModelProcessorx/SlopeCompute/CreateNodesSlopeCompute.cpp	(revision 4054)
+++ issm/trunk/src/c/modules/ModelProcessorx/SlopeCompute/CreateNodesSlopeCompute.cpp	(revision 4055)
@@ -46,5 +46,5 @@
 			
 			/*Add node to nodes dataset: */
-			nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,iomodel,SlopeComputeAnalysisEnum));
+			nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,iomodel,SlopeAnalysisEnum));
 
 		}
Index: issm/trunk/src/c/modules/OutputResultsx/OutputResultsx.cpp
===================================================================
--- issm/trunk/src/c/modules/OutputResultsx/OutputResultsx.cpp	(revision 4054)
+++ issm/trunk/src/c/modules/OutputResultsx/OutputResultsx.cpp	(revision 4055)
@@ -15,9 +15,10 @@
 #include "../../objects/objects.h"
 		
-void OutputResults(DataSet* elements, DataSet* loads, DataSet* nodes, DataSet* vertices, DataSet* materials, Parameters* parameters, char* filename){
+void OutputResults(DataSet* elements, DataSet* loads, DataSet* nodes, DataSet* vertices, DataSet* materials, Parameters* parameters){
 
 	int i;
 
 	Patch* patch=NULL;
+	char*  filename=NULL;
 
 	int solutiontype;
@@ -37,4 +38,5 @@
 	//Recover solutiontype which will be output to disk: 
 	parameters->FindParam(&solutiontype,SolutionTypeEnum);
+	parameters->FindParam(&filename,OutputFileNameEnum);
 
 	//Process results to be output in the correct units
@@ -104,4 +106,5 @@
 	/*Free ressources:*/
 	delete patch;
+	xfree((void**)&filename);
 
 }
Index: issm/trunk/src/c/modules/OutputResultsx/OutputResultsx.h
===================================================================
--- issm/trunk/src/c/modules/OutputResultsx/OutputResultsx.h	(revision 4054)
+++ issm/trunk/src/c/modules/OutputResultsx/OutputResultsx.h	(revision 4055)
@@ -9,5 +9,5 @@
 
 /* local prototypes: */
-void OutputResults(DataSet* elements, DataSet* loads, DataSet* nodes, DataSet* vertices, DataSet* materials, Parameters* parameters, char* filename);
+void OutputResults(DataSet* elements, DataSet* loads, DataSet* nodes, DataSet* vertices, DataSet* materials, Parameters* parameters);
 
 #endif  /* _OUTPUTRESULTS_H */
Index: issm/trunk/src/c/modules/UpdateGeometryx/UpdateGeometryx.cpp
===================================================================
--- issm/trunk/src/c/modules/UpdateGeometryx/UpdateGeometryx.cpp	(revision 4054)
+++ issm/trunk/src/c/modules/UpdateGeometryx/UpdateGeometryx.cpp	(revision 4055)
@@ -11,12 +11,13 @@
 #include "../../DataSet/DataSet.h"
 
-void UpdateGeometryx(Vec* poutthickness,Vec* poutbed,Vec* poutsurface, 
-		DataSet* elements, DataSet* nodes, DataSet* vertices,DataSet* loads, DataSet* materials, Parameters* parameters,
-		Vec vec_newthickness,Vec vec_bed,Vec vec_surface){
+void UpdateGeometryx( DataSet* elements, DataSet* nodes, DataSet* vertices,DataSet* loads, DataSet* materials, Parameters* parameters){
 
 	int i;
 	int dof;
 
-	/*serialized input: */
+	/*vectorized inputs: */
+	Vec     vec_newthickness=NULL;
+	Vec     vec_bed=NULL;
+	Vec     vec_surface=NULL:
 	double* newthickness=NULL;
 	double* bed=NULL;
@@ -24,5 +25,5 @@
 
 	/*objects: */
-	Node* node=NULL;
+	Vertex* vertex=NULL;
 	Object* object=NULL;
 	Matpar* matpar=NULL;
@@ -31,14 +32,4 @@
 	double h,b,s;
 	double rho_ice,rho_water;
-
-	/*output: */
-	Vec outthickness=NULL;
-	Vec outbed=NULL;
-	Vec outsurface=NULL;
-
-	/*Duplicate inputs to outputs, do not copy values!: */
-	VecDuplicate(vec_newthickness,&outthickness);
-	VecDuplicate(vec_bed,&outbed);
-	VecDuplicate(vec_surface,&outsurface);
 
 	/*First, get elements and loads configured: */
@@ -49,5 +40,9 @@
 	parameters->Configure(elements,loads, nodes,vertices, materials,parameters);
 
-	/*Serialize inputs :*/
+	/*retrieve inputs: */
+	GetVectorFromInputsx( &vec_newthickness, elements, nodes,  vertices,  loads,  materials,  parameters, ThicknessEnum,VertexEnum);
+	GetVectorFromInputsx( &vec_bed, elements, nodes,  vertices,  loads,  materials,  parameters, ThicknessEnum,VertexEnum);
+	GetVectorFromInputsx( &vec_surface, elements, nodes,  vertices,  loads,  materials,  parameters, ThicknessEnum,VertexEnum);
+
 	VecToMPISerial(&newthickness,vec_newthickness);
 	VecToMPISerial(&bed,vec_bed);
@@ -67,14 +62,14 @@
 	rho_water=matpar->GetRhoWater();
 
-	/*Go through nodes and for each node, update the thickness, bed and surface, using the 
+	/*Go through vertices and for each vertex, update the thickness, bed and surface, using the 
 	 * new thickness computed in the prognostic_core part of the transient solution sequence: */
 
-	for(i=0;i<nodes->Size();i++){
+	for(i=0;i<vertices->Size();i++){
 
 		/*get i'th node: */
-		node=(Node*)nodes->GetObjectByOffset(i);
+		vertex=(Vertex*)vertices->GetObjectByOffset(i);
 
-		/*first, recover thickness, surface and bed for this node: */
-		dof=node->GetDofList1();
+		/*first, recover thickness, surface and bed for this vertex: */
+		dof=vertex->GetDofList1();
 		h=newthickness[dof];
 		s=surface[dof];
@@ -85,10 +80,10 @@
 
 		//For grids on ice sheet, the new bed remains the same, the new surface becomes bed+new_thickness. 
-		if (node->IsOnSheet()){
+		if (vertex->IsOnSheet()){
 			s=b+h;
 		}
 
 		//For grids on ice shelt, we have hydrostatic equilibrium (for now)
-		if (node->IsOnShelf()){
+		if (vertex->IsOnShelf()){
 			b=-rho_ice/rho_water*h;
 			s=(1-rho_ice/rho_water)*h;
@@ -96,27 +91,31 @@
 
 		/*Ok, plug values of thickness, surafce and bed into our outputs: */
-		VecSetValues(outthickness,1,&dof,&h,INSERT_VALUES);
-		VecSetValues(outbed,1,&dof,&b,INSERT_VALUES);
-		VecSetValues(outsurface,1,&dof,&s,INSERT_VALUES);
+		VecSetValues(vec_newthickness,1,&dof,&h,INSERT_VALUES);
+		VecSetValues(vec_bed,1,&dof,&b,INSERT_VALUES);
+		VecSetValues(vec_surface,1,&dof,&s,INSERT_VALUES);
 	}
 
 	/*Assemble vectors: */
-	VecAssemblyBegin(outthickness);
-	VecAssemblyEnd(outthickness);
+	VecAssemblyBegin(vec_newthickness);
+	VecAssemblyEnd(vec_newthickness);
 
-	VecAssemblyBegin(outbed);
-	VecAssemblyEnd(outbed);
+	VecAssemblyBegin(vec_bed);
+	VecAssemblyEnd(vec_bed);
 
-	VecAssemblyBegin(outsurface);
-	VecAssemblyEnd(outsurface);
+	VecAssemblyBegin(vec_surface);
+	VecAssemblyEnd(vec_surface);
+
+	/*Put back into inputs:*/
+	UpdateInputsFromVectorx( elements,nodes, vertices,loads, materials,  parameters,vec_newthickness,ThicknessEnum,VertexEnum); //clobber existing thickness input with new one
+	UpdateInputsFromVectorx( elements,nodes, vertices,loads, materials,  parameters,vec_bed,BedEnum,VertexEnum); 
+	UpdateInputsFromVectorx( elements,nodes, vertices,loads, materials,  parameters,vec_surface,SurfaceEnum,VertexEnum); 
 
 	/*Free ressources:*/
+	VecFree(&vec_newthickness);
+	VecFree(&vec_bed);
+	VecFree(&vec_surface);
 	xfree((void**)&newthickness);
 	xfree((void**)&bed);
 	xfree((void**)&surface);
 	
-	/*Assign output pointers: */
-	*poutthickness=outthickness;
-	*poutbed=outbed;
-	*poutsurface=outsurface;
 }
Index: issm/trunk/src/c/modules/UpdateGeometryx/UpdateGeometryx.h
===================================================================
--- issm/trunk/src/c/modules/UpdateGeometryx/UpdateGeometryx.h	(revision 4054)
+++ issm/trunk/src/c/modules/UpdateGeometryx/UpdateGeometryx.h	(revision 4055)
@@ -10,7 +10,5 @@
 
 /* local prototypes: */
-void UpdateGeometryx(Vec* poutthickness,Vec* poutbed,Vec* poutsurface, 
-		DataSet* elements, DataSet* nodes,DataSet* vertices,DataSet* loads, DataSet* materials, Parameters* parameters,
-		Vec newthickness,Vec bed,Vec surface);
+void UpdateGeometryx(DataSet* elements, DataSet* nodes,DataSet* vertices,DataSet* loads, DataSet* materials, Parameters* parameters);
 
 #endif  /* _UPDATEGEOMETRYX_H */
Index: issm/trunk/src/c/modules/VecExtrudex/VecExtrudex.cpp
===================================================================
--- issm/trunk/src/c/modules/VecExtrudex/VecExtrudex.cpp	(revision 4055)
+++ issm/trunk/src/c/modules/VecExtrudex/VecExtrudex.cpp	(revision 4055)
@@ -0,0 +1,56 @@
+/*!\file VecExtrudex
+ * \brief: vertical velocity extrusion
+ */
+
+#include "./VecExtrudex.h"
+
+#include "../../shared/shared.h"
+#include "../../include/include.h"
+#include "../../toolkits/toolkits.h"
+#include "../../EnumDefinitions/EnumDefinitions.h"
+
+void VecExtrudex( Vec vector, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters,int collapse){
+
+	double* vector_serial=NULL;
+	int i;
+	Element* element=NULL;
+	Penta*   penta=NULL;
+	Node*    node=NULL;
+
+	/*First, get elements and nodes configured: */
+	elements->Configure(elements,loads, nodes,vertices, materials,parameters);
+	nodes->Configure(elements,loads, nodes,vertices, materials,parameters);
+	parameters->Configure(elements,loads, nodes,vertices, materials,parameters);
+
+	/*Serialize vector: */
+	VecToMPISerial(&vector_serial,vector);
+
+	/*ok, is the extrusion dependent upon elements or just nodes? If collapse is 1, it means we only want 
+	 * to extrude the vector on elemenets that are collapsed (element operation). If collapse is 0, it means 
+	 * we want to extrude vector on every node (node operation): */
+
+	if(collapse){
+		/*Extrude vector vertically: */
+		for(i=0;i<elements->Size();i++){
+			element=(Element*)elements->GetObjectByOffset(i);
+			if(element->Enum()==PentaEnum){
+				penta=(Penta*)element;
+				penta->VecExtrude(vector,vector_serial,collapse);
+			}
+		}
+	}
+	else{
+		for(i=0;i<nodes->Size();i++){
+			node=(Node*)nodes->GetObjectByOffset(i);
+			node->VecExtrude(vector,vector_serial);
+		}
+	}
+
+	/*Assemble vector: */
+	VecAssemblyBegin(vector);
+	VecAssemblyEnd(vector);
+
+	/*Free ressources:*/
+	xfree((void**)&vector_serial);
+
+}
Index: issm/trunk/src/c/modules/VecExtrudex/VecExtrudex.h
===================================================================
--- issm/trunk/src/c/modules/VecExtrudex/VecExtrudex.h	(revision 4055)
+++ issm/trunk/src/c/modules/VecExtrudex/VecExtrudex.h	(revision 4055)
@@ -0,0 +1,14 @@
+/*!\file:  VecExtrudex.h
+ * \brief header file for field extrusion
+ */ 
+
+#ifndef _VECEXTRUDEX_H
+#define _VECEXTRUDEX_H
+
+#include "../../DataSet/DataSet.h"
+
+/* local prototypes: */
+void VecExtrudex( Vec vector, DataSet* elements,DataSet* nodes, DataSet* vertices, DataSet* loads, DataSet* materials, Parameters* parameters,int collapse);
+
+#endif  /* _VECEXTRUDEX_H */
+
Index: issm/trunk/src/c/modules/modules.h
===================================================================
--- issm/trunk/src/c/modules/modules.h	(revision 4054)
+++ issm/trunk/src/c/modules/modules.h	(revision 4055)
@@ -53,5 +53,5 @@
 #include "./ComputePressurex/ComputePressurex.h"
 #include "./ComputeStrainRatex/ComputeStrainRatex.h"
-#include "./FieldExtrudex/FieldExtrudex.h"
+#include "./VecExtrudex/VecExtrudex.h"
 #include "./Qmux/Qmux.h"
 #include "./NodeConnectivityx/NodeConnectivityx.h"
@@ -81,4 +81,6 @@
 #include "./ScaleInputx/ScaleInputx.h"
 #include "./AXPYInputx/AXPYInputx.h"
+#include "./GetVectorFromInputsx/GetVectorFromInputsx.h"
+#include "./InputConvergencex/InputConvergencex.h"
 
 #endif
Index: issm/trunk/src/c/objects/Elements/Beam.cpp
===================================================================
--- issm/trunk/src/c/objects/Elements/Beam.cpp	(revision 4054)
+++ issm/trunk/src/c/objects/Elements/Beam.cpp	(revision 4055)
@@ -562,15 +562,15 @@
 /*}}}*/
 /*FUNCTION Beam::Gradj{{{1*/
-void  Beam::Gradj(int control_type){
+void  Beam::Gradj(Vec gradient,int control_type){
 	ISSMERROR(" not supported yet!");
 }
 /*}}}*/
 /*FUNCTION Beam::GradjB{{{1*/
-void  Beam::GradjB(void){
+void  Beam::GradjB(Vec gradient){
 	ISSMERROR(" not supported yet!");
 }
 /*}}}*/
 /*FUNCTION Beam::GradjDrag{{{1*/
-void  Beam::GradjDrag(void){
+void  Beam::GradjDrag(Vec gradient){
 	ISSMERROR(" not supported yet!");
 }
@@ -1055,2 +1055,30 @@
 }
 /*}}}*/
+/*FUNCTION Beam::InputConvergence(int* pconverged, double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){{{1*/
+void  Beam::InputConvergence(int* pconverged,double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){a
+
+	Input** new_inputs=NULL;
+	Input** old_inputs=NULL;
+	int     converged=1;
+
+	new_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the new inputs
+	old_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the old inputs
+	
+	for(i=0;i<num_enums/2;i++){
+		new_inputs[i]=(Input*)this->inputs->GetInput(new_enums[2*i+0]);
+		old_inputs[i]=(Input*)this->inputs->GetInput(old_enums[2*i+1]);
+		if(!new_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(new_enums[2*i+0]));
+		if(!old_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(old_enums[2*i+0]));
+	}
+
+	/*ok, we've got the inputs (new and old), now loop throught the number of criterions and fill the eps array:*/
+	for(i=0;i<num_criterionenums;i++){
+		IsInputConverged(eps+i,new_inputs,old_inputs,num_enums/2,,criterionenums[i]);
+		if(eps[i]>criterionvalues[i]) converged=0; 
+	}
+
+	/*Assign output pointers:*/
+	*pconverged=converged;
+
+}
+/*}}}*/
Index: issm/trunk/src/c/objects/Elements/Beam.h
===================================================================
--- issm/trunk/src/c/objects/Elements/Beam.h	(revision 4054)
+++ issm/trunk/src/c/objects/Elements/Beam.h	(revision 4055)
@@ -62,5 +62,5 @@
 		void  UpdateInputsFromConstant(bool constant, int name){ISSMERROR("Not implemented yet!");}
 		void  UpdateInputsFromSolution(double* solution);
-		void  DepthAverageInputAtBase(int enum_type){ISSMERROR("not implemented yet");};
+		void  DepthAverageInputAtBase(int enum_type,int average_enum_type){ISSMERROR("not implemented yet");};
 		void  InputToResult(int enum_type,int step,double time);
 		void   ProcessResultsUnits(void);
@@ -98,4 +98,5 @@
 		void  AXPYInput(int YEnum, double scalar, int XEnum);
 		void  ControlConstrainInput(int control_type,double cm_min, double cm_max);
+		void  InputConvergence(int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums);
 
 		/*}}}*/
@@ -106,7 +107,7 @@
 		void  GetThicknessList(double* thickness_list);
 		void  Du(Vec);
-		void  Gradj(int control_type);
-		void  GradjDrag(void);
-		void  GradjB(void);
+		void  Gradj(Vec gradient,int control_type);
+		void  GradjDrag(Vec gradient);
+		void  GradjB(Vec gradient);
 		double Misfit(void);
 		double SurfaceArea(void);
Index: issm/trunk/src/c/objects/Elements/Element.h
===================================================================
--- issm/trunk/src/c/objects/Elements/Element.h	(revision 4054)
+++ issm/trunk/src/c/objects/Elements/Element.h	(revision 4055)
@@ -37,11 +37,11 @@
 		virtual void   GetBedList(double* bed_list)=0;
 		virtual void   Du(Vec du_g)=0;
-		virtual void   Gradj(int control_type)=0;
-		virtual void   GradjDrag(void)=0;
-		virtual void   GradjB(void)=0;
+		virtual void   Gradj(Vec gradient,int control_type)=0;
+		virtual void   GradjDrag(Vec gradient)=0;
+		virtual void   GradjB(Vec gradient)=0;
 		virtual double Misfit(void)=0;
 		virtual double CostFunction(void)=0;
 		virtual double SurfaceArea(void)=0;
-		virtual void   DepthAverageInputAtBase(int enum_type)=0;
+		virtual void   DepthAverageInputAtBase(int enum_type,int average_enum_type)=0;
 		virtual void   ComputeBasalStress(Vec sigma_b)=0;
 		virtual void   ComputePressure(Vec p_g)=0;
@@ -69,4 +69,6 @@
 		virtual void   AXPYInput(int YEnum, double scalar, int XEnum)=0;
 		virtual void   ControlConstrainInput(int control_type,double cm_min, double cm_max)=0;
+		virtual void   InputConvergence(int* pconverged, double* eps, int* enums,int num_enums,int* criterionenums,double criterionvalues,double* criterionvalues,int num_criterionenums)=0;
+
 		/*Implementation: */
 
Index: issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 4054)
+++ issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 4055)
@@ -658,5 +658,5 @@
 
 	}
-	else if (analysis_type==SlopeComputeAnalysisEnum){
+	else if (analysis_type==SlopeAnalysisEnum){
 
 		UpdateInputsFromSolutionSlopeCompute( solution);
@@ -1270,5 +1270,5 @@
 		else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet");
 	}
-	else if (analysis_type==SlopeComputeAnalysisEnum){
+	else if (analysis_type==SlopeAnalysisEnum){
 		CreateKMatrixSlopeCompute( Kgg);
 	}
@@ -2357,5 +2357,5 @@
 		else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet");
 	}
-	else if (analysis_type==SlopeComputeAnalysisEnum){
+	else if (analysis_type==SlopeAnalysisEnum){
 		CreatePVectorSlopeCompute( pg);
 	}
@@ -3263,12 +3263,10 @@
 }
 /*}}}*/
-/*FUNCTION Penta::FieldExtrude {{{1*/
-void  Penta::FieldExtrude(Vec field,double* field_serial,char* field_name, int iscollapsed){
+/*FUNCTION Penta::VecExtrude {{{1*/
+void  Penta::VecExtrude(Vec vector,double* vector_serial,int iscollapsed){
 
 	/* node data: */
-	const int    numgrids=6;
-	int          numberofdofspernode;
+	int   i;
 	Node* node=NULL;
-	int   i;
 	int   extrude=0;
 
@@ -3296,116 +3294,35 @@
 	if(onbed==0)extrude=0;
 
-	/*Go on and extrude field: */
+	/*Go on and extrude vector: */
 	if (extrude){
 
-		if (strcmp(field_name,"velocity")==0){
-
-			/* node data: */
-			const int    numdof=2*numgrids;
-			int          doflist[numdof];
-			int          nodedofs[2];
-			double       fieldel[2];
-
-
-			GetDofList(&doflist[0],&numberofdofspernode);
-
-			/*this penta is a collapsed macayeal. For each node on the base of this penta, 
-			 * we grab the field. Once we know the field, we follow the upper nodes, 
-			 * inserting the same field value into field, until we reach the surface: */
-			for(i=0;i<3;i++){
-
-				node=nodes[i]; //base nodes
-
-				/*get field for this base node: */
-				fieldel[0]=field_serial[doflist[numberofdofspernode*i+0]];
-				fieldel[1]=field_serial[doflist[numberofdofspernode*i+1]];
-
-				//go throfieldn all nodes which sit on top of this node, until we reach the surface, 
-				//and plfield  field in field
-				for(;;){
-
-					node->GetDofList(&nodedofs[0],&numberofdofspernode);
-					VecSetValues(field,1,&nodedofs[0],&fieldel[0],INSERT_VALUES);
-					VecSetValues(field,1,&nodedofs[1],&fieldel[1],INSERT_VALUES);
-
-					if (node->IsOnSurface())break;
-					/*get next node: */
-					node=node->GetUpperNode();
-				}
-			}
-		} //if (strcmp(field_name,"velocity")==0)
-		else if (strcmp(field_name,"gradj")==0){
-
-			/* node data: */
-			int          dof1;
-			double       fieldel;
-
-			/*this penta is a collapsed macayeal. For each node on the base of this penta, 
-			 * we grab the field. Once we know the field, we follow the upper nodes, 
-			 * inserting the same field value into field, until we reach the surface: */
-			for(i=0;i<3;i++){
-
-				node=nodes[i]; //base nodes
+		/* node data: */
+		int          dof1;
+		double       vectorel;
+
+		/*this penta is a collapsed macayeal. For each node on the base of this penta, 
+		 * we grab the vector. Once we know the vector, we follow the upper nodes, 
+		 * inserting the same vector value into vector, until we reach the surface: */
+		for(i=0;i<3;i++){
+
+			node=nodes[i]; //base nodes
+			dof1=node->GetDofList1();
+
+			/*get vector for this base node: */
+			vectorel=vector_serial[dof1];
+
+			//go throvectorn all nodes which sit on top of this node, until we reach the surface, 
+			//and plvector  vector in vector
+			for(;;){
+
 				dof1=node->GetDofList1();
-
-				/*get field for this base node: */
-				fieldel=field_serial[dof1];
-
-				//go throfieldn all nodes which sit on top of this node, until we reach the surface, 
-				//and plfield  field in field
-				for(;;){
-
-					dof1=node->GetDofList1();
-					VecSetValues(field,1,&dof1,&fieldel,INSERT_VALUES);
-
-					if (node->IsOnSurface())break;
-					/*get next node: */
-					node=node->GetUpperNode();
-				}
+				VecSetValues(vector,1,&dof1,&vectorel,INSERT_VALUES);
+
+				if (node->IsOnSurface())break;
+				/*get next node: */
+				node=node->GetUpperNode();
 			}
 		}
-		else if ( 
-					(strcmp(field_name,"thickness")==0) ||
-					(strcmp(field_name,"surface")==0)  ||
-					(strcmp(field_name,"bed")==0)  ||
-					(strcmp(field_name,"slopex")==0)  ||
-					(strcmp(field_name,"slopey")==0)
-				  ){
-
-			/* node data: */
-			const int    numdof=1*numgrids;
-			int          doflist[numdof];
-			int          nodedofs;
-			double       fieldel;
-
-			GetDofList(&doflist[0],&numberofdofspernode);
-
-			/*this penta is on the bed. For each node on the base of this penta, 
-			 * we grab the thickness. Once we know the thickness, we follow the upper nodes, 
-			 * inserting the same thickness value into tg, until we reach the surface: */
-			for(i=0;i<3;i++){
-
-				node=nodes[i]; //base nodes
-
-				/*get velocity for this base node: */
-				fieldel=field_serial[doflist[numberofdofspernode*i+0]];
-
-				//go through all nodes which sit on top of this node, until we reach the surface, 
-				//and pltg  fieldel in field:
-				for(;;){
-
-					node->GetDofList(&nodedofs,&numberofdofspernode);
-					VecSetValues(field,1,&nodedofs,&fieldel,INSERT_VALUES);
-
-					if (node->IsOnSurface())break;
-					/*get next node: */
-					node=node->GetUpperNode();
-				}
-			}
-
-		}
-		else ISSMERROR("%s%s%s"," field ",field_name," not supported yet!");
-
-	} //if (extrude)
+	} 
 }
 /*}}}*/
@@ -3451,5 +3368,5 @@
 /*}}}*/
 /*FUNCTION Penta::DepthAverageInputAtBase{{{1*/
-void  Penta::DepthAverageInputAtBase(int enum_type){
+void  Penta::DepthAverageInputAtBase(int enum_type,int average_enum_type){
 	ISSMERROR("Not implemented yet (see Penta::InputExtrude and Node::FieldDepthAverageAtBase)");
 }
@@ -4555,5 +4472,5 @@
 /*}}}*/
 /*FUNCTION Penta::Gradj {{{1*/
-void  Penta::Gradj(int control_type){
+void  Penta::Gradj(Vec gradient,int control_type){
 
 	/*inputs: */
@@ -4567,8 +4484,8 @@
 
 	if (control_type==DragCoefficientEnum){
-		GradjDrag();
+		GradjDrag(gradient);
 	}
 	else if (control_type=RheologyBEnum){
-		GradjB();
+		GradjB(gradient);
 	}
 	else ISSMERROR("%s%i","control type not supported yet: ",control_type);
@@ -4576,10 +4493,10 @@
 /*}}}*/
 /*FUNCTION Penta::GradjDrag {{{1*/
-void  Penta::GradjDrag(void){
+void  Penta::GradjDrag(Vec gradient){
 
 	int i;
 	Tria* tria=NULL;
 	TriaVertexInput* triavertexinput=NULL;
-	double gradient[6]={0,0,0,0,0,0};
+	double temp_gradient[6]={0,0,0,0,0,0};
 
 	/*inputs: */
@@ -4611,5 +4528,6 @@
 		/*MacAyeal or Pattyn*/
 		tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
-		tria->GradjDrag();
+		tria->GradjDrag(gradient);
+		delete tria;
 	}
 	else if (sub_analysis_type==StokesAnalysisEnum){
@@ -4617,24 +4535,18 @@
 		/*Stokes*/
 		tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
-		tria->GradjDragStokes();
+		tria->GradjDragStokes(gradient);
+		delete tria;
 	}
 	else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet");
 
-	/*Now, the tria  has a GradientEnum input. Take it, and make a PentaVertexInput out of it, then delete the tria: */
-	triavertexinput=(TriaVertexInput*)tria->inputs->GetInput(GradientEnum);
-	for(i=0;i<3;i++)gradient[i]=triavertexinput->values[i];
-	this->inputs->AddInput(new PentaVertexInput(GradientEnum,&gradient[0]));
-	
-	delete tria;
 
 }
 /*}}}*/
 /*FUNCTION Penta::GradjB {{{1*/
-void  Penta::GradjB(void){
+void  Penta::GradjB(Vec gradient){
 
 	int i;
 	Tria* tria=NULL;
 	TriaVertexInput* triavertexinput=NULL;
-	double gradient[6]={0,0,0,0,0,0};
 
 	/*inputs: */
@@ -4658,18 +4570,14 @@
 		 * and compute gardj*/
 		tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
-		tria->GradjB();
+		tria->GradjB(gradient);
+		delete tria;
 	}
 	else{
 		/*B is a 2d field, use MacAyeal(2d) gradient even if it is Stokes or Pattyn*/
 		tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria (lower face).
-		tria->GradjB();
+		tria->GradjB(gradient);
+		delete tria;
 	}
 	
-	/*Now, the tria  has a GradientEnum input. Take it, and make a PentaVertexInput out of it, then delete the tria: */
-	triavertexinput=(TriaVertexInput*)tria->inputs->GetInput(GradientEnum);
-	for(i=0;i<3;i++)gradient[i]=triavertexinput->values[i];
-	this->inputs->AddInput(new PentaVertexInput(GradientEnum,&gradient[0]));
-	
-	delete tria;
 
 }
@@ -5407,2 +5315,30 @@
 }
 /*}}}*/
+/*FUNCTION Penta::InputConvergence(int* pconverged, double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){{{1*/
+void  Penta::InputConvergence(int* pconverged,double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){a
+
+	Input** new_inputs=NULL;
+	Input** old_inputs=NULL;
+	int     converged=1;
+
+	new_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the new inputs
+	old_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the old inputs
+	
+	for(i=0;i<num_enums/2;i++){
+		new_inputs[i]=(Input*)this->inputs->GetInput(new_enums[2*i+0]);
+		old_inputs[i]=(Input*)this->inputs->GetInput(old_enums[2*i+1]);
+		if(!new_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(new_enums[2*i+0]));
+		if(!old_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(old_enums[2*i+0]));
+	}
+
+	/*ok, we've got the inputs (new and old), now loop throught the number of criterions and fill the eps array:*/
+	for(i=0;i<num_criterionenums;i++){
+		IsInputConverged(eps+i,new_inputs,old_inputs,num_enums/2,,criterionenums[i]);
+		if(eps[i]>criterionvalues[i]) converged=0; 
+	}
+
+	/*Assign output pointers:*/
+	*pconverged=converged;
+
+}
+/*}}}*/
Index: issm/trunk/src/c/objects/Elements/Penta.h
===================================================================
--- issm/trunk/src/c/objects/Elements/Penta.h	(revision 4054)
+++ issm/trunk/src/c/objects/Elements/Penta.h	(revision 4055)
@@ -86,7 +86,7 @@
 		bool   GetOnBed();
 		void  Du(Vec du_gg);
-		void  Gradj(int control_type);
-		void  GradjDrag(void);
-		void  GradjB(void);
+		void  Gradj(Vec gradient,int control_type);
+		void  GradjDrag(Vec gradient);
+		void  GradjB(Vec gradient);
 		double Misfit(void);
 		double SurfaceArea(void);
@@ -110,7 +110,7 @@
 		void  GetParameterDerivativeValue(double* p, double* p_list,double* xyz_list, double* gauss_coord);
 		void  GetNodalFunctions(double* l1l6, double* gauss_coord);
-		void  FieldExtrude(Vec field,double* field_serial,char* field_name, int iscollapsed);
+		void  VecExtrude(Vec vector,double* vector_serial,int iscollapsed);
 		void  InputExtrude(int enum_type);
-		void  DepthAverageInputAtBase(int enum_type);
+		void  DepthAverageInputAtBase(int enum_type,int average_enum_type);
 		void  ComputeBasalStress(Vec sigma_bg);
 		void  ComputePressure(Vec p_gg);
@@ -165,4 +165,5 @@
 		void  AXPYInput(int YEnum, double scalar, int XEnum);
 		void  ControlConstrainInput(int control_type,double cm_min, double cm_max);
+		void  InputConvergence(int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums);
 		void  GetVectorFromInputs(Vec vector,int NameEnum);
 
Index: issm/trunk/src/c/objects/Elements/Sing.cpp
===================================================================
--- issm/trunk/src/c/objects/Elements/Sing.cpp	(revision 4054)
+++ issm/trunk/src/c/objects/Elements/Sing.cpp	(revision 4055)
@@ -383,15 +383,15 @@
 /*}}}*/
 /*FUNCTION Sing::Gradj {{{1*/
-void  Sing::Gradj(int control_type){
+void  Sing::Gradj(Vec gradient,int control_type){
 	ISSMERROR(" not supported yet!");
 }
 /*}}}*/
 /*FUNCTION Sing::GradB {{{1*/
-void  Sing::GradjB(void){
+void  Sing::GradjB(Vec gradient){
 	ISSMERROR(" not supported yet!");
 }
 /*}}}*/
 /*FUNCTION Sing::GradjDrag {{{1*/
-void  Sing::GradjDrag(void){
+void  Sing::GradjDrag(Vec gradient){
 	ISSMERROR(" not supported yet!");
 }
@@ -761,2 +761,30 @@
 }
 /*}}}*/
+/*FUNCTION Sing::InputConvergence(int* pconverged, double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){{{1*/
+void  Sing::InputConvergence(int* pconverged,double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){a
+
+	Input** new_inputs=NULL;
+	Input** old_inputs=NULL;
+	int     converged=1;
+
+	new_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the new inputs
+	old_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the old inputs
+	
+	for(i=0;i<num_enums/2;i++){
+		new_inputs[i]=(Input*)this->inputs->GetInput(new_enums[2*i+0]);
+		old_inputs[i]=(Input*)this->inputs->GetInput(old_enums[2*i+1]);
+		if(!new_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(new_enums[2*i+0]));
+		if(!old_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(old_enums[2*i+0]));
+	}
+
+	/*ok, we've got the inputs (new and old), now loop throught the number of criterions and fill the eps array:*/
+	for(i=0;i<num_criterionenums;i++){
+		IsInputConverged(eps+i,new_inputs,old_inputs,num_enums/2,,criterionenums[i]);
+		if(eps[i]>criterionvalues[i]) converged=0; 
+	}
+
+	/*Assign output pointers:*/
+	*pconverged=converged;
+
+}
+/*}}}*/
Index: issm/trunk/src/c/objects/Elements/Sing.h
===================================================================
--- issm/trunk/src/c/objects/Elements/Sing.h	(revision 4054)
+++ issm/trunk/src/c/objects/Elements/Sing.h	(revision 4055)
@@ -62,5 +62,5 @@
 		void  UpdateInputsFromConstant(bool constant, int name){ISSMERROR("Not implemented yet!");}
 		void  UpdateInputsFromSolution(double* solutiong);
-		void  DepthAverageInputAtBase(int enum_type){ISSMERROR("not implemented yet");};
+		void  DepthAverageInputAtBase(int enum_type,int average_enum_type){ISSMERROR("not implemented yet");};
 		void  InputToResult(int enum_type,int step,double time);
 		void   ProcessResultsUnits(void);
@@ -96,4 +96,5 @@
 		void  ScaleInput(int enum_type,double scale_factor);
 		void  AXPYInput(int YEnum, double scalar, int XEnum);
+		void  InputConvergence(int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums);
 		void  ControlConstrainInput(int control_type,double cm_min, double cm_max);
 
@@ -106,7 +107,7 @@
 		void  GetThicknessList(double* thickness_list);
 		void  Du(Vec);
-		void  Gradj(int control_type);
-		void  GradjDrag(void);
-		void  GradjB(void);
+		void  Gradj(Vec gradient,int control_type);
+		void  GradjDrag(Vec gradient);
+		void  GradjB(Vec gradient);
 		double Misfit(void);
 		double SurfaceArea(void);
Index: issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 4054)
+++ issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 4055)
@@ -569,5 +569,5 @@
 		else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet");
 	}
-	else if (analysis_type==SlopeComputeAnalysisEnum){
+	else if (analysis_type==SlopeAnalysisEnum){
 		UpdateInputsFromSolutionSlopeCompute( solution);
 	}
@@ -979,5 +979,5 @@
 		else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet");
 	}
-	else if (analysis_type==SlopeComputeAnalysisEnum){
+	else if (analysis_type==SlopeAnalysisEnum){
 		CreateKMatrixSlopeCompute( Kgg);
 	}
@@ -2394,5 +2394,5 @@
 		else ISSMERROR("%s%i%s\n","sub_analysis: ",sub_analysis_type," not supported yet");
 	}
-	else if (analysis_type==SlopeComputeAnalysisEnum){
+	else if (analysis_type==SlopeAnalysisEnum){
 		CreatePVectorSlopeCompute( pg);
 	}
@@ -4125,5 +4125,5 @@
 /*}}}*/
 /*FUNCTION Tria::Gradj {{{1*/
-void  Tria::Gradj(int control_type){
+void  Tria::Gradj(Vec gradient,int control_type){
 
 	/*inputs: */
@@ -4137,8 +4137,8 @@
 
 	if (control_type==DragCoefficientEnum){
-		GradjDrag();
+		GradjDrag(gradient);
 	}
 	else if (control_type==RheologyBEnum){
-		GradjB();
+		GradjB(gradient);
 	}
 	else ISSMERROR("%s%i","control type not supported yet: ",control_type);
@@ -4146,5 +4146,5 @@
 /*}}}*/
 /*FUNCTION Tria::GradjDrag {{{1*/
-void  Tria::GradjDrag(void){
+void  Tria::GradjDrag(Vec gradient){
 
 
@@ -4295,6 +4295,6 @@
 
 
-	/*Add grade_g to the inputs of this element: */
-	this->inputs->AddInput(new TriaVertexInput(GradientEnum,&grade_g[0]));
+	/*Add grade_g to global vector gradient: */
+	VecSetValues(gradient,numgrids,doflist1,(const double*)grade_g,ADD_VALUES);
 
 	cleanup_and_return: 
@@ -4308,5 +4308,5 @@
 /*}}}*/
 /*FUNCTION Tria::GradjDragStokes {{{1*/
-void  Tria::GradjDragStokes(void){
+void  Tria::GradjDragStokes(Vec gradient){
 
 	int i;
@@ -4465,4 +4465,7 @@
 	}
 
+	/*Add grade_g to global vector gradient: */
+	VecSetValues(gradient,numgrids,doflist1,(const double*)grade_g,ADD_VALUES);
+
 	/*Add grade_g to the inputs of this element: */
 	this->inputs->AddInput(new TriaVertexInput(GradientEnum,&grade_g[0]));
@@ -4478,5 +4481,5 @@
 /*}}}*/
 /*FUNCTION Tria::GradjB{{{1*/
-void  Tria::GradjB(void){
+void  Tria::GradjB(Vec gradient){
 
 	int i;
@@ -4613,6 +4616,6 @@
 	}
 
-	/*Add grade_g to the inputs of this element: */
-	this->inputs->AddInput(new TriaVertexInput(GradientEnum,&grade_g[0]));
+	/*Add grade_g to global vector gradient: */
+	VecSetValues(gradient,numgrids,doflist1,(const double*)grade_g,ADD_VALUES);
 
 	cleanup_and_return: 
@@ -5014,5 +5017,5 @@
 /*}}}*/
 /*FUNCTION Tria::DepthAverageInputAtBase {{{1*/
-void  Tria::DepthAverageInputAtBase(int enum_type){
+void  Tria::DepthAverageInputAtBase(int enum_type,int average_enum_type){
 
 	/*New input*/
@@ -5026,17 +5029,5 @@
 
 	/*Assign new name (average)*/
-	switch(enum_type){
-
-		case VxEnum:
-			newinput->ChangeEnum(VxAverageEnum);
-			break;
-
-		case VyEnum:
-			newinput->ChangeEnum(VyAverageEnum);
-			break;
-
-		default:
-			ISSMERROR("enum_type: %i (%s) not suppoerted yet",enum_type,EnumAsString(enum_type));
-	}
+	newinput->ChangeEnum(average_enum_type);
 
 	/*Add new input to current element*/
@@ -5503,2 +5494,30 @@
 }
 /*}}}*/
+/*FUNCTION Tria::InputConvergence(int* pconverged, double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){{{1*/
+void  Tria::InputConvergence(int* pconverged,double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){a
+
+	Input** new_inputs=NULL;
+	Input** old_inputs=NULL;
+	int     converged=1;
+
+	new_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the new inputs
+	old_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the old inputs
+	
+	for(i=0;i<num_enums/2;i++){
+		new_inputs[i]=(Input*)this->inputs->GetInput(new_enums[2*i+0]);
+		old_inputs[i]=(Input*)this->inputs->GetInput(old_enums[2*i+1]);
+		if(!new_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(new_enums[2*i+0]));
+		if(!old_inputs[i])ISSMERROR("%s%s"," could not find input with enum ",EnumAsString(old_enums[2*i+0]));
+	}
+
+	/*ok, we've got the inputs (new and old), now loop throught the number of criterions and fill the eps array:*/
+	for(i=0;i<num_criterionenums;i++){
+		IsInputConverged(eps+i,new_inputs,old_inputs,num_enums/2,,criterionenums[i]);
+		if(eps[i]>criterionvalues[i]) converged=0; 
+	}
+
+	/*Assign output pointers:*/
+	*pconverged=converged;
+
+}
+/*}}}*/
Index: issm/trunk/src/c/objects/Elements/Tria.h
===================================================================
--- issm/trunk/src/c/objects/Elements/Tria.h	(revision 4054)
+++ issm/trunk/src/c/objects/Elements/Tria.h	(revision 4055)
@@ -55,5 +55,5 @@
 		int   MyRank(void);
 		void  SetClone(int* minranks);
-		void  DepthAverageInputAtBase(int enum_type);
+		void  DepthAverageInputAtBase(int enum_type,int average_enum_type);
 		void*  SpawnSing(int g0);
 		void*  SpawnBeam(int g0, int g1);
@@ -88,8 +88,8 @@
 		void  GetJacobian(double* J, double* xyz_list,double* gauss_l1l2l3);
 		void  Du(Vec du_g);
-		void  Gradj(int control_type);
-		void  GradjDrag(void);
-		void  GradjDragStokes(void);
-		void  GradjB(void);
+		void  Gradj(Vec gradient,int control_type);
+		void  GradjDrag(Vec gradient);
+		void  GradjDragStokes(Vec gradient);
+		void  GradjB(Vec gradient);
 		void  SurfaceNormal(double* surface_normal, double xyz_list[3][3]);
 		double Misfit(void);
@@ -142,4 +142,5 @@
 		void  ScaleInput(int enum_type,double scale_factor);
 		void  AXPYInput(int YEnum, double scalar, int XEnum);
+		void  InputConvergence(int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums);
 		void  ControlConstrainInput(int control_type,double cm_min, double cm_max);
 		void  GetVectorFromInputs(Vec vector,int NameEnum);
Index: issm/trunk/src/c/objects/FemModel.cpp
===================================================================
--- issm/trunk/src/c/objects/FemModel.cpp	(revision 4054)
+++ issm/trunk/src/c/objects/FemModel.cpp	(revision 4055)
@@ -35,17 +35,17 @@
 	/*Dynamically allocate whatever is a list of length nummodels: */
 	analysis_type_list=(int*)xmalloc(nummodels*sizeof(int));
-	Rmg=(Mat*)xmalloc(nummodels*sizeof(Mat));
-	Gmn=(Mat*)xmalloc(nummodels*sizeof(Mat));
-	nodesets=(NodeSets**)xmalloc(nummodels*sizeof(NodeSets*));
-	yg=(Vec*)xmalloc(nummodels*sizeof(Vec));
-	ys=(Vec*)xmalloc(nummodels*sizeof(Vec));
+	m_Rmg=(Mat*)xmalloc(nummodels*sizeof(Mat));
+	m_Gmn=(Mat*)xmalloc(nummodels*sizeof(Mat));
+	m_nodesets=(NodeSets**)xmalloc(nummodels*sizeof(NodeSets*));
+	m_yg=(Vec*)xmalloc(nummodels*sizeof(Vec));
+	m_ys=(Vec*)xmalloc(nummodels*sizeof(Vec));
 
 	/*Initialize: */
 	for(i=0;i<nummodels;i++)analysis_type_list[i]=analyses[i];
-	for(i=0;i<nummodels;i++)Rmg[i]=NULL;
-	for(i=0;i<nummodels;i++)Gmn[i]=NULL;
-	for(i=0;i<nummodels;i++)nodesets[i]=NULL;
-	for(i=0;i<nummodels;i++)yg[i]=NULL;
-	for(i=0;i<nummodels;i++)ys[i]=NULL;
+	for(i=0;i<nummodels;i++)m_Rmg[i]=NULL;
+	for(i=0;i<nummodels;i++)m_Gmn[i]=NULL;
+	for(i=0;i<nummodels;i++)m_nodesets[i]=NULL;
+	for(i=0;i<nummodels;i++)m_yg[i]=NULL;
+	for(i=0;i<nummodels;i++)m_ys[i]=NULL;
 
 	_printf_("   fill model with matlab workspace data\n");
@@ -64,17 +64,17 @@
 
 		_printf_("   create single point constraints: \n");
-		SpcNodesx( &yg[i], nodes,constraints,analysis_type); 
+		SpcNodesx( &m_yg[i], nodes,constraints,analysis_type); 
 
 		_printf_("   create rigid body constraints:\n");
-		MpcNodesx( &Rmg[i], nodes,constraints,analysis_type); 
+		MpcNodesx( &m_Rmg[i], nodes,constraints,analysis_type); 
 
 		_printf_("   create node sets:\n");
-		BuildNodeSetsx(&nodesets[i], nodes,analysis_type);
+		BuildNodeSetsx(&m_nodesets[i], nodes,analysis_type);
 
 		_printf_("   reducing single point constraints vector:\n");
-		Reducevectorgtosx(&ys[i], yg[i],nodesets[i]);
+		Reducevectorgtosx(&m_ys[i], m_yg[i],m_nodesets[i]);
 
 		_printf_("   normalizing rigid body constraints matrix:\n");
-		NormalizeConstraintsx(&Gmn[i], Rmg[i],nodesets[i]);
+		NormalizeConstraintsx(&m_Gmn[i], m_Rmg[i],m_nodesets[i]);
 
 		_printf_("   configuring element and loads:\n");
@@ -106,22 +106,22 @@
 
 	for(i=0;i<nummodels;i++){
-		Mat temp_Rmg=Rmg[i];
+		Mat temp_Rmg=m_Rmg[i];
 		MatFree(&temp_Rmg);
-		Mat temp_Gmn=Gmn[i];
+		Mat temp_Gmn=m_Gmn[i];
 		MatFree(&temp_Gmn);
-		NodeSets* temp_nodesets=nodesets[i];
-		delete nodesets;
-		Vec temp_yg=yg[i];
+		NodeSets* temp_nodesets=m_nodesets[i];
+		delete temp_nodesets;
+		Vec temp_yg=m_yg[i];
 		VecFree(&temp_yg);
-		Vec temp_ys=ys[i];
+		Vec temp_ys=m_ys[i];
 		VecFree(&temp_ys);
 	}
 
 	/*Delete dynamically allocated arrays: */
-	delete Rmg;
-	delete Gmn;
-	delete nodesets;
-	delete yg;
-	delete ys;
+	delete m_Rmg;
+	delete m_Gmn;
+	delete m_nodesets;
+	delete m_yg;
+	delete m_ys;
 
 }
@@ -151,4 +151,5 @@
 /*FUNCTION FemModel::SetCurrentAnalysis {{{1*/
 void FemModel::SetCurrentAnalysis(int analysis_type){
+
 	int found=-1;
 	for(int i=0;i<nummodels;i++){
@@ -160,4 +161,25 @@
 	if(found!=-1) analysis_counter=found;
 	else ISSMERROR("Could not find analysis_type %s in list of FemModel analyses",EnumAsString(analysis_type));
+
+	/*activate matrices/vectors: */
+	Rmg=m_Rmg[analysis_counter];
+	Gmn=m_Gmn[analysis_counter];
+	nodesets=m_nodesets[analysis_counter];
+	yg=m_yg[analysis_counter];
+	ys=m_ys[analysis_counter];
+
+}
+/*}}}1*/
+/*FUNCTION FemModel::SetCurrentAnalysisAlias {{{1*/
+void FemModel::SetCurrentAnalysisAlias(int base_analysis_type,int real_analysis_type){
+	int found=-1;
+	for(int i=0;i<nummodels;i++){
+		if (analysis_type_list[i]==base_analysis_type){
+			found=i;
+			break;
+		}
+	}
+	if(found!=-1) analysis_counter=found;
+	else ISSMERROR("Could not find alias for analysis_type %s in list of FemModel analyses",EnumAsString(base_analysis_type));
 }
 /*}}}1*/
Index: issm/trunk/src/c/objects/FemModel.h
===================================================================
--- issm/trunk/src/c/objects/FemModel.h	(revision 4054)
+++ issm/trunk/src/c/objects/FemModel.h	(revision 4055)
@@ -40,10 +40,17 @@
 		DofVec*             tpartition;
 		
-		//multiple  sets of matrices/vectors for each analysis_type
-		Mat*                 Rmg; //rigid body motions matrices
-		Mat*                 Gmn;
-		NodeSets**           nodesets; //boundary conditions dof sets
-		Vec*                 yg; //boundary conditions in global g-set
-		Vec*                 ys; //boundary conditions, in reduced s-set
+		//multiple  sets of matrices/vectors for each analysis_type. m stands for multiple
+		Mat*                 m_Rmg; //rigid body motions matrices
+		Mat*                 m_Gmn;
+		NodeSets**           m_nodesets; //boundary conditions dof sets
+		Vec*                 m_yg; //boundary conditions in global g-set
+		Vec*                 m_ys; //boundary conditions, in reduced s-set
+
+		//pointers to point to sets of matrices/vectors, for a certain analysis type. activated in SetCurrentAnalysis
+		Mat                  Rmg;
+		Mat                  Gmn;
+		NodeSets*            nodesets;
+		Vec                  yg;
+		Vec                  ys;
 
 		/*constructors, destructors: */
@@ -56,4 +63,5 @@
 		/*Fem: */
 		void  SetCurrentAnalysis(int analysis_type);
+		void  SetCurrentAnalysisAlias(int base_analysis_type,int real_analysis_type);
 		int   GetCurrentAnalysis(void);
 
Index: issm/trunk/src/c/objects/Inputs/BeamVertexInput.cpp
===================================================================
--- issm/trunk/src/c/objects/Inputs/BeamVertexInput.cpp	(revision 4054)
+++ issm/trunk/src/c/objects/Inputs/BeamVertexInput.cpp	(revision 4055)
@@ -294,2 +294,10 @@
 }
 /*}}}*/
+/*FUNCTION BeamVertexInput::GetValuesPtr(double* pvalues,int* pnum_values);{{{1*/
+void BeamVertexInput::GetValuesPtr(double* pvalues,int* pnum_values){
+
+	*pvalues=this->values;
+	*pnum_values=2;
+
+}
+/*}}}*/
Index: issm/trunk/src/c/objects/Inputs/BeamVertexInput.h
===================================================================
--- issm/trunk/src/c/objects/Inputs/BeamVertexInput.h	(revision 4054)
+++ issm/trunk/src/c/objects/Inputs/BeamVertexInput.h	(revision 4055)
@@ -78,4 +78,5 @@
 		void Constrain(double cm_min, double cm_max);
 		void GetVectorFromInputs(Vec vector,int* doflist);
+		void GetValuesPtr(double* pvalues,int* pnum_values);
 		/*}}}*/
 
Index: issm/trunk/src/c/objects/Inputs/BoolInput.cpp
===================================================================
--- issm/trunk/src/c/objects/Inputs/BoolInput.cpp	(revision 4054)
+++ issm/trunk/src/c/objects/Inputs/BoolInput.cpp	(revision 4055)
@@ -256,2 +256,9 @@
 }
 /*}}}*/
+/*FUNCTION BoolInput::GetValuesPtr(double* pvalues,int* pnum_values);{{{1*/
+void BoolInput::GetValuesPtr(double* pvalues,int* pnum_values){
+
+	ISSMERROR(" not supported yet!");
+
+}
+/*}}}*/
Index: issm/trunk/src/c/objects/Inputs/BoolInput.h
===================================================================
--- issm/trunk/src/c/objects/Inputs/BoolInput.h	(revision 4054)
+++ issm/trunk/src/c/objects/Inputs/BoolInput.h	(revision 4055)
@@ -78,4 +78,5 @@
 		void Constrain(double cm_min, double cm_max);
 		void GetVectorFromInputs(Vec vector,int* doflist);
+		void GetValuesPtr(double* pvalues,int* pnum_values);
 		/*}}}*/
 
Index: issm/trunk/src/c/objects/Inputs/DoubleInput.cpp
===================================================================
--- issm/trunk/src/c/objects/Inputs/DoubleInput.cpp	(revision 4054)
+++ issm/trunk/src/c/objects/Inputs/DoubleInput.cpp	(revision 4055)
@@ -267,2 +267,9 @@
 }
 /*}}}*/
+/*FUNCTION DoubleInput::GetValuesPtr(double* pvalues,int* pnum_values);{{{1*/
+void DoubleInput::GetValuesPtr(double* pvalues,int* pnum_values){
+
+	ISSMERROR(" not supported yet!");
+
+}
+/*}}}*/
Index: issm/trunk/src/c/objects/Inputs/DoubleInput.h
===================================================================
--- issm/trunk/src/c/objects/Inputs/DoubleInput.h	(revision 4054)
+++ issm/trunk/src/c/objects/Inputs/DoubleInput.h	(revision 4055)
@@ -79,4 +79,5 @@
 		void Constrain(double cm_min, double cm_max);
 		void GetVectorFromInputs(Vec vector,int* doflist);
+		void GetValuesPtr(double* pvalues,int* pnum_values);
 		/*}}}*/
 
Index: issm/trunk/src/c/objects/Inputs/Input.h
===================================================================
--- issm/trunk/src/c/objects/Inputs/Input.h	(revision 4054)
+++ issm/trunk/src/c/objects/Inputs/Input.h	(revision 4055)
@@ -55,4 +55,5 @@
 		virtual void Constrain(double cm_min, double cm_max)=0;
 		virtual void GetVectorFromInputs(Vec vector,int* doflist)=0;
+		virtual void GetValuesPtr(double* pvalues,int* pnum_values)=0;
 		/*}}}*/
 
Index: issm/trunk/src/c/objects/Inputs/IntInput.cpp
===================================================================
--- issm/trunk/src/c/objects/Inputs/IntInput.cpp	(revision 4054)
+++ issm/trunk/src/c/objects/Inputs/IntInput.cpp	(revision 4055)
@@ -257,2 +257,9 @@
 }
 /*}}}*/
+/*FUNCTION IntInput::GetValuesPtr(double* pvalues,int* pnum_values);{{{1*/
+void IntInput::GetValuesPtr(double* pvalues,int* pnum_values){
+
+	ISSMERROR(" not supported yet!");
+
+}
+/*}}}*/
Index: issm/trunk/src/c/objects/Inputs/IntInput.h
===================================================================
--- issm/trunk/src/c/objects/Inputs/IntInput.h	(revision 4054)
+++ issm/trunk/src/c/objects/Inputs/IntInput.h	(revision 4055)
@@ -78,4 +78,5 @@
 		void Constrain(double cm_min, double cm_max);
 		void GetVectorFromInputs(Vec vector,int* doflist);
+		void GetValuesPtr(double* pvalues,int* pnum_values);
 		/*}}}*/
 
Index: issm/trunk/src/c/objects/Inputs/PentaVertexInput.cpp
===================================================================
--- issm/trunk/src/c/objects/Inputs/PentaVertexInput.cpp	(revision 4054)
+++ issm/trunk/src/c/objects/Inputs/PentaVertexInput.cpp	(revision 4055)
@@ -943,2 +943,10 @@
 }
 /*}}}*/
+/*FUNCTION PentaVertexInput::GetValuesPtr(double* pvalues,int* pnum_values);{{{1*/
+void PentaVertexInput::GetValuesPtr(double* pvalues,int* pnum_values){
+
+	*pvalues=this->values;
+	*pnum_values=6;
+
+}
+/*}}}*/
Index: issm/trunk/src/c/objects/Inputs/PentaVertexInput.h
===================================================================
--- issm/trunk/src/c/objects/Inputs/PentaVertexInput.h	(revision 4054)
+++ issm/trunk/src/c/objects/Inputs/PentaVertexInput.h	(revision 4055)
@@ -87,4 +87,5 @@
 		void Constrain(double cm_min, double cm_max);
 		void GetVectorFromInputs(Vec vector,int* doflist);
+		void GetValuesPtr(double* pvalues,int* pnum_values);
 		/*}}}*/
 
Index: issm/trunk/src/c/objects/Inputs/SingVertexInput.cpp
===================================================================
--- issm/trunk/src/c/objects/Inputs/SingVertexInput.cpp	(revision 4054)
+++ issm/trunk/src/c/objects/Inputs/SingVertexInput.cpp	(revision 4055)
@@ -258,3 +258,11 @@
 
 }
-/*}}}*/
+/2*}}}*/
+/*FUNCTION SingVertexInput::GetValuesPtr(double* pvalues,int* pnum_values);{{{1*/
+void SingVertexInput::GetValuesPtr(double* pvalues,int* pnum_values){
+
+	*pvalues=&this->value;
+	*pnum_values=1;
+
+}
+/*}}}*/
Index: issm/trunk/src/c/objects/Inputs/SingVertexInput.h
===================================================================
--- issm/trunk/src/c/objects/Inputs/SingVertexInput.h	(revision 4054)
+++ issm/trunk/src/c/objects/Inputs/SingVertexInput.h	(revision 4055)
@@ -77,4 +77,5 @@
 		void Constrain(double cm_min, double cm_max);
 		void GetVectorFromInputs(Vec vector,int* doflist);
+		void GetValuesPtr(double* pvalues,int* pnum_values);
 		/*}}}*/
 
Index: issm/trunk/src/c/objects/Inputs/TriaVertexInput.cpp
===================================================================
--- issm/trunk/src/c/objects/Inputs/TriaVertexInput.cpp	(revision 4054)
+++ issm/trunk/src/c/objects/Inputs/TriaVertexInput.cpp	(revision 4055)
@@ -517,2 +517,10 @@
 }
 /*}}}*/
+/*FUNCTION TriaVertexInput::GetValuesPtr(double* pvalues,int* pnum_values);{{{1*/
+void TriaVertexInput::GetValuesPtr(double* pvalues,int* pnum_values){
+
+	*pvalues=this->values;
+	*pnum_values=3;
+
+}
+/*}}}*/
Index: issm/trunk/src/c/objects/Inputs/TriaVertexInput.h
===================================================================
--- issm/trunk/src/c/objects/Inputs/TriaVertexInput.h	(revision 4054)
+++ issm/trunk/src/c/objects/Inputs/TriaVertexInput.h	(revision 4055)
@@ -84,4 +84,5 @@
 		void Constrain(double cm_min, double cm_max);
 		void GetVectorFromInputs(Vec vector,int* doflist);
+		void GetValuesPtr(double* pvalues,int* pnum_values);
 		/*}}}*/
 
Index: issm/trunk/src/c/objects/Node.cpp
===================================================================
--- issm/trunk/src/c/objects/Node.cpp	(revision 4054)
+++ issm/trunk/src/c/objects/Node.cpp	(revision 4055)
@@ -144,5 +144,5 @@
 				analysis_type==PrognosticAnalysisEnum || 
 				analysis_type==MeltingAnalysisEnum || 
-				analysis_type==SlopeComputeAnalysisEnum || 
+				analysis_type==SlopeAnalysisEnum || 
 				analysis_type==BalancedvelocitiesAnalysisEnum || 
 				analysis_type==BalancedthicknessAnalysisEnum
@@ -668,6 +668,6 @@
 }
 /*}}}*/
-/*FUNCTION Node::FieldExtrude {{{2*/
-void  Node::FieldExtrude(Vec field,double* field_serial,char* field_name){
+/*FUNCTION Node::VecExtrude {{{2*/
+void  Node::VecExtrude(Vec vector,double* vector_serial){
 		
 	/* node data: */
@@ -684,103 +684,27 @@
 	if (onbed){
 
-		if (strcmp(field_name,"velocity")==0){
-
-			/* node data: */
-			int          vertexdof;
-			int          doflist[2];
-			double       fieldnode[2];
-
-			vertexdof=GetVertexDof();
-
-			/*get dofs for this base node velocity: we know there are two dofs in field_serial */
-			doflist[0]=2*vertexdof;
-			doflist[1]=2*vertexdof+1;
-
-			//initilaize node
-			node=this;
-		
-			/*get field for this base node: */
-			fieldnode[0]=field_serial[doflist[0]];
-			fieldnode[1]=field_serial[doflist[1]];
-
-			//go through all nodes which sit on top of this node, until we reach the surface, 
-			//and plfield  field in field
-			for(;;){
-
-				vertexdof=node->GetVertexDof();
-				doflist[0]=2*vertexdof;
-				doflist[1]=2*vertexdof+1;
-				VecSetValues(field,1,&doflist[0],&fieldnode[0],INSERT_VALUES);
-				VecSetValues(field,1,&doflist[1],&fieldnode[1],INSERT_VALUES);
-
-				if (node->IsOnSurface())break;
-				/*get next node: */
-				node=node->GetUpperNode();
-			}
-		} //if (strcmp(field_name,"velocity")==0)
-		else if (strcmp(field_name,"gradj")==0){
-
-			/* node data: */
-			int          dof1;
-			double       fieldel;
-
-			//initilaize node and get dof1
-			node=this;
+		/* node data: */
+		int          dof1;
+		double       vectorel;
+
+		//initilaize node and get dof1
+		node=this;
+		dof1=node->GetVertexDof();
+
+		/*get vector for this base node: */
+		vectorel=vector_serial[dof1];
+
+		//go throvectorn all nodes which sit on top of this node, until we reach the surface, 
+		//and plvector  vector in vector
+		for(;;){
+
 			dof1=node->GetVertexDof();
-
-			/*get field for this base node: */
-			fieldel=field_serial[dof1];
-
-			//go throfieldn all nodes which sit on top of this node, until we reach the surface, 
-			//and plfield  field in field
-			for(;;){
-
-				dof1=node->GetVertexDof();
-				VecSetValues(field,1,&dof1,&fieldel,INSERT_VALUES);
-
-				if (node->IsOnSurface())break;
-				/*get next node: */
-				node=node->GetUpperNode();
-			}
-		}
-		else if ( 
-				(strcmp(field_name,"vx")==0) ||
-				(strcmp(field_name,"vy")==0) ||
-				(strcmp(field_name,"thickness")==0) ||
-				(strcmp(field_name,"temperature")==0) ||
-				(strcmp(field_name,"surface")==0)  ||
-				(strcmp(field_name,"bed")==0)  || 
-				(strcmp(field_name,"slopex")==0)  ||
-				(strcmp(field_name,"slopey")==0) 
-				){
-
-			/* node data: */
-			int          doflist;
-			int          nodedofs;
-			double       fieldnode;
-				
-			GetDofList(&doflist,&numberofdofspernode);
-
-			//initialize node
-			node=this;
-
-			/*get field for this bed node: */
-			fieldnode=field_serial[doflist];
-
-			//go through all nodes which sit on top of this node, until we reach the surface, 
-			//and plug  fieldnode in field:
-			for(;;){
-
-				node->GetDofList(&nodedofs,&numberofdofspernode);
-				VecSetValues(field,1,&nodedofs,&fieldnode,INSERT_VALUES);
-
-				if (node->IsOnSurface())break;
-				/*get next node: */
-				node=node->GetUpperNode();
-			}
-		}
-		else ISSMERROR("%s%s%s"," field ",field_name," not supported yet!");
-
-	} //if (extrude)
+			VecSetValues(vector,1,&dof1,&vectorel,INSERT_VALUES);
+
+			if (node->IsOnSurface())break;
+			/*get next node: */
+			node=node->GetUpperNode();
+		}
+	} 
 }
 /*}}}*/
Index: issm/trunk/src/c/objects/Node.h
===================================================================
--- issm/trunk/src/c/objects/Node.h	(revision 4054)
+++ issm/trunk/src/c/objects/Node.h	(revision 4055)
@@ -89,5 +89,5 @@
 		int   IsOnShelf();
 		int   IsOnSheet();
-		void  FieldExtrude(Vec field,double* field_serial,char* field_name);
+		void  VecExtrude(Vec vector,double* vector_serial);
 		/*}}}*/
 		/*FUNCTION DofObject routines {{{1*/
Index: issm/trunk/src/c/shared/Dofs/DistributeNumDofs.cpp
===================================================================
--- issm/trunk/src/c/shared/Dofs/DistributeNumDofs.cpp	(revision 4054)
+++ issm/trunk/src/c/shared/Dofs/DistributeNumDofs.cpp	(revision 4055)
@@ -25,5 +25,5 @@
 		numdofs=2;
 	}
-	else if (analysis_type==SlopeComputeAnalysisEnum){
+	else if (analysis_type==SlopeAnalysisEnum){
 		numdofs=1;
 	}
Index: issm/trunk/src/c/shared/Numerics/IsInputConverged.cpp
===================================================================
--- issm/trunk/src/c/shared/Numerics/IsInputConverged.cpp	(revision 4055)
+++ issm/trunk/src/c/shared/Numerics/IsInputConverged.cpp	(revision 4055)
@@ -0,0 +1,54 @@
+/*!\file:  IsInputConverged.cpp
+ * \brief: figure out convergence at the input level. 
+ * We are given a list of inputs, new and old, and a criterion, and using the inputs, we 
+ * return the value of the criterion test, which will be used at the solution level 
+ * to determine convergence.
+ */ 
+
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "../../objects/objects.h"
+
+void IsInputConverged(double* peps, Input** new_inputs,Input** old_inputs,int num_inputs,int criterion_enum){
+
+	int i;
+
+	/*output: */
+	double eps;
+	
+	/*intermediary: */
+	double *newvalues     = NULL;
+	double *oldvalues     = NULL;
+	int     num_values;
+	double  ndu        = 0;
+	double  nu         = 0;
+
+	if(criterion_enum==RelativeEnum){
+
+		/*conpute ndu/du (where u could be velocity, pressure, temperature, etc ...): */
+		for(i=0;i<num_inputs;i++){
+			new_inputs[i]->GetValuesPtr(&newvalues,&num_values);
+			old_inputs[i]->GetValuesPtr(&oldvalues,&num_values);
+			for(j=0;j<num_values;j++){
+				ndu+=pow(newvalues[j]-oldvalues[j],2);
+				nu+=pow(oldvalues[j],2);
+			}
+		}
+
+		/*take square root: */
+		ndu=sqrt(ndu);
+		nu=sqrt(nu);
+
+		/*now, compute eps: */
+		if(nu)eps=ndu/nu;
+		else eps=0;
+	}
+	else ISSMERROR("%s%s%s"," convergence criterion ",EnumAsString(criterion_enum)," not supported yet!");
+
+	/*Assign output pointers:*/
+	*peps=eps;
+}
Index: issm/trunk/src/c/shared/Numerics/numerics.h
===================================================================
--- issm/trunk/src/c/shared/Numerics/numerics.h	(revision 4054)
+++ issm/trunk/src/c/shared/Numerics/numerics.h	(revision 4055)
@@ -8,4 +8,5 @@
 #include "./GaussPoints.h"
 #include "./isnan.h"
+class Input;
 
 struct OptArgs;
@@ -18,4 +19,5 @@
 void   cross(double* result,double* vector1,double* vector2);
 double norm(double* vector);
+void IsInputConverged(double* peps, Input** new_inputs,Input** old_inputs,int num_inputs,int criterion_enum);
 
 #endif //ifndef _NUMERICS_H_
Index: issm/trunk/src/c/solutions/adjoint_core.cpp
===================================================================
--- issm/trunk/src/c/solutions/adjoint_core.cpp	(revision 4054)
+++ issm/trunk/src/c/solutions/adjoint_core.cpp	(revision 4055)
@@ -20,4 +20,5 @@
 	char* solverstring=NULL;
 	bool  isstokes=false;
+	bool conserve_loads=true;
 	
 	/*intermediary: */
@@ -38,12 +39,12 @@
 
 	/*set analysis type to compute velocity: */
-	if(isstokes)femmodel->SetAnalysis(DiagnosticAnalysisEnum,DiagnosticStokesAnalysisEnum);
-	else femmodel->SetAnalysis(DiagnosticAnalysisEnum,DiagnosticHorizAnalysisEnum);
+	if(isstokes)femmodel->SetCurrentAnalysis(DiagnosticStokesAnalysisEnum);
+	else femmodel->SetCurrentAnalysis(DiagnosticHorizAnalysisEnum);
 	
 	_printf_("%s\n","      recover solution for this stiffness and right hand side:");
-	solver_diagnostic_nonlinear(NULL,&K_ff0,&K_fs0,NULL, femmodel,DiagnosticAnalysisEnum,sub_analysis_type);
+	solver_diagnostic_nonlinear(NULL,&K_ff0,&K_fs0, femmodel,conserve_loads); 
 
 	_printf_("%s\n","      buid Du, difference between observed velocity and model velocity:");
-	Dux(du_g, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
+	Dux(&du_g, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
 
 	_printf_("%s\n","      reduce adjoint load from g-set to f-set:");
@@ -66,14 +67,14 @@
 
 	/*Update inputs using adjoint solution, and same type of setup as diagnostic solution: */
-	if(isstokes)femmodel->SetAnalysisAlias(DiagnosticAnalysisEnum,DiagnosticStokesAnalysisEnum,AdjointAnalysisEnum,NoneAnalysisEnum);
-	else femmodel->SetAnalysisAlias(DiagnosticAnalysisEnum,DiagnosticHorizAnalysisEnum,AdjointAnalysisEnum,NoneAnalysisEnum);
+	if(isstokes)femmodel->SetCurrentAnalysisAlias(DiagnosticStokesAnalysisEnum,AdjointAnalysisEnum);
+	else femmodel->SetCurrentAnalysisAlias(DiagnosticHorizAnalysisEnum,AdjointAnalysisEnum);
 	
-	UpdateInputsFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,adjoint_g,AdjointAnalysisEnum, NoneAnalysisEnum);
+	UpdateInputsFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,adjoint_g);
 
 	/*Free ressources:*/
-	char* solverstring=NULL;
-	VecFree(&ug);
+	xfree((void**)&solverstring);
+	VecFree(&u_g);
 	MatFree(&K_ff0);
-	MatFree(&Mat K_fs0);
+	MatFree(&K_fs0);
 	VecFree(&du_g);
 	VecFree(&du_f);
Index: issm/trunk/src/c/solutions/balancedthickness.cpp
===================================================================
--- issm/trunk/src/c/solutions/balancedthickness.cpp	(revision 4054)
+++ issm/trunk/src/c/solutions/balancedthickness.cpp	(revision 4055)
@@ -22,17 +22,10 @@
 	char* inputfilename=NULL;
 	char* outputfilename=NULL;
+	bool  qmu_analysis=false;
 	char* lockname=NULL;
-	int   numberofnodes;
 	bool  waitonlock=false;
 
-	Model* model=NULL;
-
-	double  dt;
-	double  yts;
-	bool    qmu_analysis;
-
-	/*Results: */
-	Results* results=NULL;
-	Param*   param=NULL;
+	/*FemModel: */
+	FemModel* femmodel=NULL;
 
 	/*time*/
@@ -40,4 +33,7 @@
 	double   start_core, finish_core;
 	double   start_init, finish_init;
+
+	int analyses[1]={BalancedthicknessAnalysisEnum};
+	int solution_type=BalancedthicknessAnalysisEnum;
 
 	MODULEBOOT();
@@ -66,13 +62,14 @@
 	fid=pfopen(inputfilename,"rb");
 
-	/*Initialize model structure: */
-	model=new Model();
+	_printf_("create finite element model:\n");
+	femmodel=new FemModel(fid,solution_type,analyses,1);
 
-	_printf_("read and create finite element model:\n");
-	model->AddFormulation(fid,BalancedthicknessAnalysisEnum);
-
-	/*recover parameters: */
-	model->FindParam(&waitonlock,WaitOnLockEnum);
-	model->FindParam(&qmu_analysis,QmuAnalysisEnum);
+	/*add outputfilename in parameters: */
+	femmodel->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename);
+	
+	/*get parameters: */
+	femmodel->parameters->FindParam(&qmu_analysis,QmuAnalysisEnum);
+	femmodel->parameters->FindParam(&control_analysis,ControlAnalysisEnum);
+	femmodel->parameters->FindParam(&waitonlock,WaitOnLockEnum);
 
 	MPI_Barrier(MPI_COMM_WORLD); finish_init=MPI_Wtime();
@@ -84,6 +81,10 @@
 		_printf_("call computational core:\n");
 		MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
-		results=balancedthickness_core(model);
+		balancedthickness_core(model);
 		MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
+
+		_printf_("write results to disk:\n");
+		OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters, outputfilename);
+	}
 
 	}
@@ -95,5 +96,5 @@
 		#ifdef _HAVE_DAKOTA_ 
 		MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
-		Qmux(model,BalancedthicknessAnalysisEnum,NoneAnalysisEnum);
+		Qmux(femmodel,BalancedthicknessAnalysisEnum,NoneAnalysisEnum);
 		MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
 	 	#else
@@ -101,7 +102,4 @@
 		#endif
 	}
-
-	_printf_("write results to disk:\n");
-	OutputResults(results,outputfilename);
 
 	if (waitonlock>0){
@@ -111,6 +109,5 @@
 
 	/*Free ressources:*/
-	delete results;
-	delete model;
+	delete femmodel;
 
 	/*Get finish time and close*/
Index: issm/trunk/src/c/solutions/balancedthickness2.cpp
===================================================================
--- issm/trunk/src/c/solutions/balancedthickness2.cpp	(revision 4054)
+++ issm/trunk/src/c/solutions/balancedthickness2.cpp	(revision 4055)
@@ -23,18 +23,10 @@
 	char* inputfilename=NULL;
 	char* outputfilename=NULL;
+	bool  qmu_analysis=false;
 	char* lockname=NULL;
 	bool  waitonlock=false;
 
-	Model* model=NULL;
-
-	double  dt;
-	double  yts;
-	bool    qmu_analysis;
-
-	/*Results: */
-	Results* results=NULL;
-	Result*  result=NULL;
-
-	Param*   param=NULL;
+	/*FemModel: */
+	FemModel* femmodel=NULL;
 
 	/*time*/
@@ -42,4 +34,8 @@
 	double   start_core, finish_core;
 	double   start_init, finish_init;
+
+	int analyses[1]={Balancedthickness2AnalysisEnum};
+	int solution_type=Balancedthickness2AnalysisEnum;
+
 
 	MODULEBOOT();
@@ -68,13 +64,9 @@
 	fid=pfopen(inputfilename,"rb");
 
-	/*Initialize model structure: */
-	model=new Model();
+	_printf_("create finite element model:\n");
+	femmodel=new FemModel(fid,solution_type,analyses,1);
 
-	_printf_("read and create finite element model:\n");
-	model->AddFormulation(fid,Balancedthickness2AnalysisEnum);
-
-	/*recover parameters: */
-	model->FindParam(&waitonlock,WaitOnLockEnum);
-	model->FindParam(&qmu_analysis,QmuAnalysisEnum);
+	/*add outputfilename in parameters: */
+	femmodel->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename);
 
 	MPI_Barrier(MPI_COMM_WORLD); finish_init=MPI_Wtime();
@@ -86,6 +78,9 @@
 		_printf_("call computational core:\n");
 		MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
-		results=balancedthickness2_core(model);
+		balancedthickness2_core(femmodel);
 		MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
+
+		_printf_("write results to disk:\n");
+		OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters, outputfilename);
 
 	}
@@ -97,5 +92,5 @@
 		#ifdef _HAVE_DAKOTA_ 
 		MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
-		Qmux(model,Balancedthickness2AnalysisEnum,NoneAnalysisEnum);
+		Qmux(femmodel,Balancedthickness2AnalysisEnum,NoneAnalysisEnum);
 		MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
 	 	#else
@@ -103,7 +98,4 @@
 		#endif
 	}
-
-	_printf_("write results to disk:\n");
-	OutputResults(results,outputfilename);
 
 	if (waitonlock>0){
@@ -113,6 +105,5 @@
 
 	/*Free ressources:*/
-	delete results;
-	delete model;
+	delete femmodel;
 
 	/*Get finish time and close*/
Index: issm/trunk/src/c/solutions/balancedthickness2_core.cpp
===================================================================
--- issm/trunk/src/c/solutions/balancedthickness2_core.cpp	(revision 4054)
+++ issm/trunk/src/c/solutions/balancedthickness2_core.cpp	(revision 4055)
@@ -12,71 +12,35 @@
 #include "../solvers/solvers.h"
 
-Results* balancedthickness2_core(Model* model){
-
-	extern int my_rank;
-
-	/*output: */
-	Results* results=NULL;
-	Result* result=NULL;
-
-	/*intermediary: */
-	Vec vx_g=NULL;
-	Vec vy_g=NULL;
-
-	/*solutions: */
-	Vec h_g=NULL;
+void balancedthickness2_core(FemModel* femmodel){
 
 	/*flags: */
 	int verbose=0;
-	int numberofdofspernode;
-	int numberofnodes;
-	int numberofvertices;
-	int dofs[1]={1};
+	int dim;
 
-	/*fem balancedthickness2 model: */
-	FemModel* fem_p=NULL;
-
-	//initialize results:
-	results=new Results();
-
-	fem_p=model->GetFormulation(Balancedthickness2AnalysisEnum);
-
-	//first recover parameters common to all solutions
-	model->FindParam(&verbose,VerboseEnum);
-	model->FindParam(&numberofnodes,NumberOfNodesEnum);
-	model->FindParam(&numberofvertices,NumberOfVerticesEnum);
-	model->FindParam(&numberofdofspernode,NumberOfDofsPerNodeEnum);
+	/*activate formulation: */
+	femmodel->SetCurrentAnalysis(BalancedthicknessAnalysisEnum);
+	
+	/*recover parameters: */
+	femmodel->parameters->FindParam(&verbose,VerboseEnum);
+	femmodel->parameters->FindParam(&dim,DimEnum);
 
 	_printf_("depth averaging velocity...\n");
-	//vx_g=inputs->Get("vx",&dofs[0],1);
-	//vy_g=inputs->Get("vy",&dofs[0],1);
-	/* NOT WORKING YET....
-	FieldDepthAveragex(vx_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"vx");
-	FieldDepthAveragex(vy_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"vy");
-	*/
-	//inputs->Add("vx_average",vx_g,1,numberofvertices);
-	//inputs->Add("vy_average",vy_g,1,numberofvertices);
-	ISSMERROR("not supported yet!");
-	
+	DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VxEnum);
+	DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VyEnum);
+	if(dim==3) DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VzEnum);
+
+
 	_printf_("call computational core:\n");
-	diagnostic_core_linear(&h_g,fem_p,Balancedthickness2AnalysisEnum,NoneAnalysisEnum);
+	solver_linear(NULL,femmodel);
 
 	_printf_("Averaging over vertices:\n");
-	FieldAverageOntoVerticesx(&h_g,fem_p->elements,fem_p->nodes,fem_p->vertices,fem_p->loads,fem_p->materials,fem_p->parameters);
+	ISSMERROR(" not supported yet!");
+//FieldAverageOntoVerticesx(&h_g,fem_p->elements,fem_p->nodes,fem_p->vertices,fem_p->loads,fem_p->materials,fem_p->parameters);
 
 //	_printf_("extrude computed thickness on all layers:\n");
 //	FieldExtrudex( h_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"thickness",0);
 
-	/*Plug results into output dataset: */
-	InputToResultx(&result,fem_p->elements,fem_p->nodes,fem_p->vertices, fem_p->loads, fem_p->materials,fem_p->parameters,ThicknessEnum,results->Size()+1,0,1); results->AddObject(result);
-	results->AddObject(new StringResult(results->Size()+1,AnalysisTypeEnum,0,1,EnumAsString(Balancedthickness2AnalysisEnum)));
-
-	/*Free ressources:*/
-	VecFree(&vx_g);
-	VecFree(&vy_g);
-	VecFree(&h_g);
-	
-	/*return: */
-	return results;
+	if(verbose)_printf_("saving results:\n");
+	InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum);
 
 }
Index: issm/trunk/src/c/solutions/balancedthickness_core.cpp
===================================================================
--- issm/trunk/src/c/solutions/balancedthickness_core.cpp	(revision 4054)
+++ issm/trunk/src/c/solutions/balancedthickness_core.cpp	(revision 4055)
@@ -12,60 +12,30 @@
 #include "../solvers/solvers.h"
 
-Results* balancedthickness_core(Model* model){
+void balancedthickness_core(FemModel* femmodel){
 
-	extern int my_rank;
+	/*parameters: */
+	int verbose=0;
+	int dim;
 
-	/*output: */
-	Results* results=NULL;
-	Result* result=NULL;
-
-	/*intermediary: */
-	Vec u_g=NULL;
-
-	/*solutions: */
-	Vec h_g=NULL;
-
-	/*flags: */
-	int verbose=0;
-	int numberofdofspernode;
-	int numberofnodes;
-	int dofs[2]={1,1};
-
-	/*fem balancedthickness model: */
-	FemModel* fem_p=NULL;
-
-	//initialize results:
-	results=new Results();
-
-	/*recover fem model: */
-	fem_p=model->GetFormulation(BalancedthicknessAnalysisEnum);
-
-	//first recover parameters common to all solutions
-	model->FindParam(&verbose,VerboseEnum);
-	model->FindParam(&numberofnodes,NumberOfNodesEnum);
-	model->FindParam(&numberofdofspernode,NumberOfDofsPerNodeEnum);
+	/*activate formulation: */
+	femmodel->SetCurrentAnalysis(BalancedthicknessAnalysisEnum);
+	
+	/*recover parameters: */
+	femmodel->parameters->FindParam(&verbose,VerboseEnum);
+	femmodel->parameters->FindParam(&dim,DimEnum);
 
 	_printf_("depth averaging velocity...\n");
-	//u_g=inputs->Get("velocity",&dofs[0],2); //take (vx,vy) from inputs velocity
-	//FieldDepthAveragex( u_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"velocity");
-	//inputs->Add("velocity_average",u_g,2,numberofnodes);
-	ISSMERROR("not supported yet!");
-	
+	DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VxEnum);
+	DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VyEnum);
+	if(dim==3) DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VzEnum);
+
 	_printf_("call computational core:\n");
-	diagnostic_core_linear(&h_g,fem_p,BalancedthicknessAnalysisEnum,NoneAnalysisEnum);
+	solver_linear(NULL,femmodel);
 
 	_printf_("extrude computed thickness on all layers:\n");
-	FieldExtrudex( h_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"thickness",0);
+	ExtrudeInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,ThicknessEnum);
 
-	/*Plug results into output dataset: */
-	InputToResultx(&result,fem_p->elements,fem_p->nodes,fem_p->vertices, fem_p->loads, fem_p->materials,fem_p->parameters,ThicknessEnum,results->Size()+1,0,1); results->AddObject(result);
-	results->AddObject(new StringResult(results->Size()+1,AnalysisTypeEnum,0,1,EnumAsString(BalancedthicknessAnalysisEnum)));
-
-	/*Free ressources:*/
-	VecFree(&u_g);
-	VecFree(&h_g);
-	
-	/*return: */
-	return results;
+	if(verbose)_printf_("saving results:\n");
+	InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum);
 
 }
Index: issm/trunk/src/c/solutions/balancedvelocities.cpp
===================================================================
--- issm/trunk/src/c/solutions/balancedvelocities.cpp	(revision 4054)
+++ issm/trunk/src/c/solutions/balancedvelocities.cpp	(revision 4055)
@@ -23,16 +23,10 @@
 	char* inputfilename=NULL;
 	char* outputfilename=NULL;
+	bool   qmu_analysis;
 	char* lockname=NULL;
-	int   numberofnodes;
 	bool  waitonlock=false;
 
-	Model* model=NULL;
-
-	bool   qmu_analysis;
-
-	/*Results: */
-	Results* results=NULL;
-	Result*  result=NULL;
-	Param*   param=NULL;
+	/*FemModel: */
+	FemModel* femmodel=NULL;
 
 	/*time*/
@@ -40,4 +34,7 @@
 	double   start_core, finish_core;
 	double   start_init, finish_init;
+
+	int analyses[1]={BalancedvelocitiesAnalysisEnum};
+	int solution_type=BalancedvelocitiesAnalysisEnum;
 
 	MODULEBOOT();
@@ -66,13 +63,14 @@
 	fid=pfopen(inputfilename,"rb");
 
-	/*Initialize model structure: */
-	model=new Model();
+	_printf_("create finite element model:\n");
+	femmodel=new FemModel(fid,solution_type,analyses,1);
 
-	_printf_("read and create finite element model:\n");
-	model->AddFormulation(fid,BalancedvelocitiesAnalysisEnum);
-
-	/*recover parameters: */
-	model->FindParam(&waitonlock,WaitOnLockEnum);
-	model->FindParam(&qmu_analysis,QmuAnalysisEnum);
+	/*add outputfilename in parameters: */
+	femmodel->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename);
+	
+	/*get parameters: */
+	femmodel->parameters->FindParam(&qmu_analysis,QmuAnalysisEnum);
+	femmodel->parameters->FindParam(&control_analysis,ControlAnalysisEnum);
+	femmodel->parameters->FindParam(&waitonlock,WaitOnLockEnum);
 
 	MPI_Barrier(MPI_COMM_WORLD); finish_init=MPI_Wtime();
@@ -84,6 +82,9 @@
 		_printf_("call computational core:\n");
 		MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
-		results=balancedvelocities_core(model);
+		balancedvelocities_core(femmodel);
 		MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
+
+		_printf_("write results to disk:\n");
+		OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters, outputfilename);
 
 	}
@@ -95,5 +96,5 @@
 		#ifdef _HAVE_DAKOTA_ 
 		MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
-		Qmux(model,BalancedvelocitiesAnalysisEnum,NoneAnalysisEnum);
+		Qmux(femmodel,BalancedvelocitiesAnalysisEnum,NoneAnalysisEnum);
 		MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
 	 	#else
@@ -101,7 +102,4 @@
 		#endif
 	}
-
-	_printf_("write results to disk:\n");
-	OutputResults(results,outputfilename);
 
 	if (waitonlock>0){
@@ -111,6 +109,5 @@
 
 	/*Free ressources:*/
-	delete results;
-	delete model;
+	delete femmodel;
 
 	/*Get finish time and close*/
Index: issm/trunk/src/c/solutions/balancedvelocities_core.cpp
===================================================================
--- issm/trunk/src/c/solutions/balancedvelocities_core.cpp	(revision 4054)
+++ issm/trunk/src/c/solutions/balancedvelocities_core.cpp	(revision 4055)
@@ -11,60 +11,32 @@
 #include "../solvers/solvers.h"
 
-Results* balancedvelocities_core(Model* model){
-
-	extern int my_rank;
-
-	/*output: */
-	Results* results=NULL;
-	Result* result=NULL;
-
-	/*intermediary: */
-	Vec u_g=NULL;
-
-	/*solutions: */
-	Vec v_g=NULL;
+void balancedvelocities_core(FemModel* femmodel){
 
 	/*flags: */
 	int verbose=0;
-	int numberofdofspernode;
-	int numberofnodes;
-	int dofs[2]={1,1};
+	int dim;
 
-	/*fem balancedvelocities model: */
-	FemModel* fem_p=NULL;
-
-	//initialize results:
-	results=new Results();
-
-	/*recover fem model: */
-	fem_p=model->GetFormulation(BalancedvelocitiesAnalysisEnum);
-
-	//first recover parameters common to all solutions
-	model->FindParam(&verbose,VerboseEnum);
-	model->FindParam(&numberofnodes,NumberOfNodesEnum);
-	model->FindParam(&numberofdofspernode,NumberOfDofsPerNodeEnum);
+	/*activate formulation: */
+	femmodel->SetCurrentAnalysis(BalancedvelocitiesAnalysisEnum);
+	
+	/*recover parameters: */
+	femmodel->parameters->FindParam(&verbose,VerboseEnum);
+	femmodel->parameters->FindParam(&dim,DimEnum);
 
 	_printf_("depth averaging velocity...\n");
-	//u_g=inputs->Get("velocity",&dofs[0],2); //take (vx,vy) from inputs velocity
-	//FieldDepthAveragex( u_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"velocity");
-	//inputs->Add("velocity_average",u_g,2,numberofnodes);
-	ISSMERROR(" not supported yet!");
-	
+	DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VxEnum);
+	DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VyEnum);
+	if(dim==3) DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VzEnum);
+
 	_printf_("call computational core:\n");
-	diagnostic_core_linear(&v_g,fem_p,BalancedvelocitiesAnalysisEnum,NoneAnalysisEnum);
+	solver_linear(NULL,femmodel);
 
-	_printf_("extrude computed thickness on all layers:\n");
-	FieldExtrudex( v_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"thickness",0);
+	_printf_("extrude computed velocity on all layers:\n");
+	ExtrudeInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VxEnum);
+	ExtrudeInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VyEnum);
 
-	/*Plug results into output dataset: */
-	InputToResultx(&result,fem_p->elements,fem_p->nodes,fem_p->vertices, fem_p->loads, fem_p->materials,fem_p->parameters,VxEnum,results->Size()+1,0,1); results->AddObject(result);
-	InputToResultx(&result,fem_p->elements,fem_p->nodes,fem_p->vertices, fem_p->loads, fem_p->materials,fem_p->parameters,VyEnum,results->Size()+1,0,1); results->AddObject(result);
-	results->AddObject(new StringResult(results->Size()+1,AnalysisTypeEnum,0,1,EnumAsString(BalancedvelocitiesAnalysisEnum)));
+	if(verbose)_printf_("saving results:\n");
+	InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum);
+	InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum);
 
-	/*Free ressources:*/
-	VecFree(&u_g);
-	VecFree(&v_g);
-
-	/*return: */
-	return results;
 }
Index: issm/trunk/src/c/solutions/bedslope.cpp
===================================================================
--- issm/trunk/src/c/solutions/bedslope.cpp	(revision 4054)
+++ issm/trunk/src/c/solutions/bedslope.cpp	(revision 4055)
@@ -37,6 +37,6 @@
 	double   start_init, finish_init;
 
-	int analyses[1]={BedSlopeComputeAnalysisEnum};
-	int solution_type=BedSlopeComputeSolutionEnum;
+	int analyses[1]={SlopeAnalysisEnum};
+	int solution_type=SlopeAnalysisEnum;
 
 	MODULEBOOT();
@@ -67,4 +67,7 @@
 	_printf_("create finite element model, using analyses types statically defined above:\n");
 	femmodel=new FemModel(fid,solution_type,analyses,1);
+
+	/*add outputfilename in parameters: */
+	femmodel->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename));
 	
 	/*get parameters: */
@@ -83,5 +86,5 @@
 
 		_printf_("write results to disk:\n");
-		OutputResults(femmodel,outputfilename,DiagnosticAnalysisEnum);
+		OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters);
 	}
 	else{
@@ -92,5 +95,5 @@
 		#ifdef _HAVE_DAKOTA_ 
 		MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
-		Qmux(femmodel,SlopeComputeAnalysisEnum,NoneAnalysisEnum);
+		Qmux(femmodel,BedSlopeAnalysisEnum,NoneAnalysisEnum);
 		MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
 	 	#else
Index: issm/trunk/src/c/solutions/bedslope_core.cpp
===================================================================
--- issm/trunk/src/c/solutions/bedslope_core.cpp	(revision 4054)
+++ issm/trunk/src/c/solutions/bedslope_core.cpp	(revision 4055)
@@ -25,12 +25,19 @@
 
 	/*Call on core computations: */
-	solver_linear(NULL,femmodel,SlopeComputeAnalysisEnum,BedSlopeXAnalysisEnum);
-	solver_linear(NULL,femmodel,SlopeComputeAnalysisEnum,BedSlopeYAnalysisEnum);
+	femmodel->SetCurrentAnalysisAlias(SlopeAnalysisEnum,BedSlopeXAnalysisEnum);
+	solver_linear(NULL,femmodel);
+	femmodel->SetCurrentAnalysisAlias(SlopeAnalysisEnum,BedSlopeYAnalysisEnum);
+	solver_linear(NULL,femmodel);
 	
 	/*extrude inputs if we are in 3D: */
 	if (dim==3){
-		if(verbose)_printf_("%s\n","extruding slope in 3d...");
-		InputsExtrudex( femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedSlopeXEnum,0);
-		InputsExtrudex( femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedSlopeYEnum,0);
+		if(verbose)_printf_("%s\n","extruding bed in 3d...");
+		ExtrudeInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,BedSlopeXEnum);
+		ExtrudeInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,BedSlopeYEnum);
 	}
+	
+	if(verbose)_printf_("saving results:\n");
+	InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedSlopeXEnum);
+	InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedSlopeYEnum);
+
 }
Index: issm/trunk/src/c/solutions/control_core.cpp
===================================================================
--- issm/trunk/src/c/solutions/control_core.cpp	(revision 4054)
+++ issm/trunk/src/c/solutions/control_core.cpp	(revision 4055)
@@ -104,5 +104,5 @@
 		if (((n+1)%5)==0){
 			_printf_("%s","      saving temporary results...");
-			controlrestart(femmodel);
+			controlrestart(femmodel,J);
 		}
 	}
Index: issm/trunk/src/c/solutions/controlrestart.cpp
===================================================================
--- issm/trunk/src/c/solutions/controlrestart.cpp	(revision 4054)
+++ issm/trunk/src/c/solutions/controlrestart.cpp	(revision 4055)
@@ -7,20 +7,20 @@
 #include "../EnumDefinitions/EnumDefinitions.h"
 
-void controlrestart(FemModel* femmodel){
+void controlrestart(FemModel* femmodel,double* J){
 
-	char*    outputfilename=NULL;
-	
+	int      control_type;
+	int      nsteps;
+
 	/*retrieve output file name: */
-	model->FindParam(&outputfilename,OutputFileNameEnum);
+	femmodel->parameters->FindParam(&control_type,ControlTypeEnum);
+	femmodel->parameters->FindParam(&nsteps,NStepsEnum);
 
 	/*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)));
+	femmodel->results->AddObject(new DoubleVecExternalResult(femmodel->results->Size()+1,JEnum,J,nsteps,1,0));
+	femmodel->results->AddObject(new StringExternalResult(femmodel->results->Size()+1,ControlTypeEnum,EnumAsString(control_type),1,0));
 
 	/*write to disk: */
-	OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters, outputfilename);
-	
-	/*Free ressources:*/
-	xfree((void**)&outputfilename);
+	OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters);
+
 }
Index: issm/trunk/src/c/solutions/diagnostic.cpp
===================================================================
--- issm/trunk/src/c/solutions/diagnostic.cpp	(revision 4054)
+++ issm/trunk/src/c/solutions/diagnostic.cpp	(revision 4055)
@@ -36,5 +36,5 @@
 	double   start_init, finish_init;
 
-	int analyses[5]={DiagnosticHorizAnalysisEnum,DiagnosticVertAnalysisEnum,DiagnosticStokesAnalysisEnum,DiagnosticHutterAnalysisEnum,SlopeComputeAnalysisEnum};
+	int analyses[5]={DiagnosticHorizAnalysisEnum,DiagnosticVertAnalysisEnum,DiagnosticStokesAnalysisEnum,DiagnosticHutterAnalysisEnum,SlopeAnalysisEnum};
 	int solution_type=DiagnosticAnalysisEnum;
 
@@ -58,5 +58,4 @@
 	lockname=argv[4];
 
-	/*Initialize femmodel structure: */
 	MPI_Barrier(MPI_COMM_WORLD); start_init=MPI_Wtime();
 
@@ -66,4 +65,7 @@
 	_printf_("create finite element model:\n");
 	femmodel=new FemModel(fid,solution_type,analyses,5);
+
+	/*add outputfilename in parameters: */
+	femmodel->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename));
 	
 	/*get parameters: */
@@ -94,5 +96,5 @@
 
 		_printf_("write results to disk:\n");
-		OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters, outputfilename);
+		OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters);
 	}
 	else{
Index: issm/trunk/src/c/solutions/diagnostic_core.cpp
===================================================================
--- issm/trunk/src/c/solutions/diagnostic_core.cpp	(revision 4054)
+++ issm/trunk/src/c/solutions/diagnostic_core.cpp	(revision 4055)
@@ -22,6 +22,8 @@
 	bool isstokes=false;
 	double stokesreconditioning;
+	bool conserve_loads=true;
+	bool modify_loads=true;
 
-	/* recover parameters: {{{1*/
+	/* recover parameters:*/
 	femmodel->parameters->FindParam(&verbose,VerboseEnum);
 	femmodel->parameters->FindParam(&dim,DimEnum);
@@ -31,5 +33,4 @@
 	femmodel->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum);
 	femmodel->parameters->FindParam(&qmu_analysis,QmuAnalysisEnum);
-	/*}}}*/
 
 	/*for qmu analysis, reinitialize velocity so that fake sensitivities do not show up as a result of a different restart of the convergence at each trial.*/
@@ -48,5 +49,6 @@
 			
 		if(verbose)_printf_("%s\n"," computing hutter velocities...");
-		solver_linear(NULL,femmodel,DiagnosticAnalysisEnum,HutterAnalysisEnum); 
+		femmodel->SetCurrentAnalysis(DiagnosticHutterAnalysisEnum);
+		solver_linear(NULL,femmodel);
 		
 		if (ismacayealpattyn) ResetBoundaryConditions(femmodel,DiagnosticAnalysisEnum,HorizAnalysisEnum);
@@ -56,5 +58,6 @@
 		
 		if(verbose)_printf_("%s\n"," computing horizontal velocities...");
-		solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,false,DiagnosticAnalysisEnum,HorizAnalysisEnum); //false means we expose loads to changes inside the solution
+		femmodel->SetCurrentAnalysis(DiagnosticHorizAnalysisEnum);
+		solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,modify_loads); 
 	}
 	
@@ -63,5 +66,6 @@
 
 		if(verbose)_printf_("%s\n"," computing vertical velocities...");
-		solver_linear(NULL,femmodel,DiagnosticAnalysisEnum,VertAnalysisEnum);
+		femmodel->SetCurrentAnalysis(DiagnosticVertAnalysisEnum);
+		solver_linear(NULL,femmodel);
 
 		if (isstokes){
@@ -75,5 +79,6 @@
 
 			if(verbose)_printf_("%s\n"," computing stokes velocities and pressure ...");
-			solver_diagnostic_nonlinear(NULL,NULL,NULL,NULL,femmodel,DiagnosticAnalysisEnum,StokesAnalysisEnum);
+			femmodel->SetCurrentAnalysis(DiagnosticStokesAnalysisEnum);
+			solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,conserve_loads);
 		}
 	}
Index: issm/trunk/src/c/solutions/gradient_core.cpp
===================================================================
--- issm/trunk/src/c/solutions/gradient_core.cpp	(revision 4054)
+++ issm/trunk/src/c/solutions/gradient_core.cpp	(revision 4055)
@@ -32,10 +32,10 @@
 	
 	if(verbose)_printf_("%s\n","      compute gradient:");
-	Gradjx(&gradient, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters, control_type,AdjointEnum);
+	Gradjx(&gradient, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters, control_type);
 	
 	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(control_steady)diagnostic_core(femmodel);
 	
 	if (step>0 && search_scalar==0){
@@ -55,6 +55,6 @@
 
 	/*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);
+	UpdateInputsFromVectorx( femmodel-> elements,femmodel-> nodes, femmodel-> vertices,femmodel-> loads, femmodel-> materials,  femmodel->parameters,gradient,GradientEnum,VertexEnum);
+	UpdateInputsFromVectorx( femmodel-> elements,femmodel-> nodes, femmodel-> vertices,femmodel-> loads, femmodel-> materials,  femmodel->parameters,old_gradient,OldGradientEnum,VertexEnum);
 
 	/*Free ressources and return:*/
Index: issm/trunk/src/c/solutions/objectivefunctionC.cpp
===================================================================
--- issm/trunk/src/c/solutions/objectivefunctionC.cpp	(revision 4054)
+++ issm/trunk/src/c/solutions/objectivefunctionC.cpp	(revision 4055)
@@ -3,8 +3,5 @@
  */ 
 
-#include "../modules/modules.h"
-#include "./solutions.h"
-#include "../EnumDefinitions/EnumDefinitions.h"
-
+/*include files: {{{1*/
 #ifdef HAVE_CONFIG_H
 	#include "config.h"
@@ -12,4 +9,14 @@
 #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
 #endif
+
+#include "../toolkits/toolkits.h"
+#include "../objects/objects.h"
+#include "../shared/shared.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+#include "../solvers/solvers.h"
+#include "./solutions.h"
+#include "../modules/modules.h"
+#include "../include/include.h"
+/*}}}*/
 
 double objectivefunctionC(double search_scalar,OptArgs* optargs){
@@ -22,5 +29,5 @@
 	/*parameters: */
 	FemModel* femmodel=NULL;
-	int n;
+	int     n;
 	double* optscal=NULL;
 	double* fit=NULL;
@@ -29,7 +36,7 @@
 	int     control_type;
 	int     analysis_type;
-	int     sub_analysis_type;
 	int     control_steady;
-		
+	bool    conserve_loads=true;
+	
 	/*Recover finite element model: */
 	femmodel=optargs->femmodel;
@@ -45,5 +52,4 @@
 	femmodel->parameters->FindParam(&control_steady,ControlSteadyEnum);
 	femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
-	femmodel->parameters->FindParam(&sub_analysis_type,SubAnalysisTypeEnum);
 
 	/*Use ControlParameterEnum input to  reinitialize our input parameter: */
@@ -57,10 +63,10 @@
 
 	/*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
+	if(!control_steady) solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,conserve_loads);  //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.*/
 	UpdateInputsFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,fit[n],FitEnum);
-	CostFunctionx( &J, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters,analysis_type,sub_analysis_type);
+	CostFunctionx( &J, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters);
 
 	/*Free ressources:*/
Index: issm/trunk/src/c/solutions/prognostic.cpp
===================================================================
--- issm/trunk/src/c/solutions/prognostic.cpp	(revision 4054)
+++ issm/trunk/src/c/solutions/prognostic.cpp	(revision 4055)
@@ -24,15 +24,9 @@
 	char* outputfilename=NULL;
 	char* lockname=NULL;
-	int   numberofnodes;
+	bool  qmu_analysis=false;
 	bool  waitonlock=false;
 
-	Model* model=NULL;
-
-	bool qmu_analysis;
-
-	/*Results: */
-	Results* results=NULL;
-
-	Param*   param=NULL;
+	/*FemModel: */
+	FemModel* femmodel=NULL;
 
 	/*time*/
@@ -40,4 +34,7 @@
 	double   start_core, finish_core;
 	double   start_init, finish_init;
+
+	int analyses[1]={PrognosticAnalysisEnum};
+	int solution_type=PrognosticAnalysisEnum;
 
 	MODULEBOOT();
@@ -62,15 +59,18 @@
 	/*Initialize model structure: */
 	MPI_Barrier(MPI_COMM_WORLD); start_init=MPI_Wtime();
-	model=new Model();
 
 	/*Open handle to data on disk: */
 	fid=pfopen(inputfilename,"rb");
 
-	_printf_("read and create finite element model:\n");
-	model->AddFormulation(fid,PrognosticAnalysisEnum);
+	_printf_("create finite element model:\n");
+	femmodel=new FemModel(fid,solution_type,analyses,1);
 
-	/*recover parameters: */
-	model->FindParam(&waitonlock,WaitOnLockEnum);
-	model->FindParam(&qmu_analysis,QmuAnalysisEnum);
+	/*add outputfilename in parameters: */
+	femmodel->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename));
+
+	/*get parameters: */
+	femmodel->parameters->FindParam(&qmu_analysis,QmuAnalysisEnum);
+	femmodel->parameters->FindParam(&waitonlock,WaitOnLockEnum);
+
 
 	MPI_Barrier(MPI_COMM_WORLD); finish_init=MPI_Wtime();
@@ -82,6 +82,9 @@
 		_printf_("call computational core:\n");
 		MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
-		results=prognostic_core(model);
+		prognostic_core(femmodel);
 		MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
+		
+		_printf_("write results to disk:\n");
+		OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters);
 
 	}
@@ -93,5 +96,5 @@
 		#ifdef _HAVE_DAKOTA_ 
 		MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
-		Qmux(model,PrognosticAnalysisEnum,NoneAnalysisEnum);
+		Qmux(femmodel,PrognosticAnalysisEnum,NoneAnalysisEnum);
 		MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
 	 	#else
@@ -101,7 +104,4 @@
 
 		
-	_printf_("write results to disk:\n");
-	OutputResults(results,outputfilename);
-
 	if (waitonlock>0){
 		_printf_("write lock file:\n");
@@ -110,6 +110,5 @@
 
 	/*Free ressources:*/
-	delete results;
-	delete model;
+	delete femmodel;
 
 	/*Get finish time and close*/
Index: issm/trunk/src/c/solutions/prognostic2.cpp
===================================================================
--- issm/trunk/src/c/solutions/prognostic2.cpp	(revision 4054)
+++ issm/trunk/src/c/solutions/prognostic2.cpp	(revision 4055)
@@ -24,15 +24,9 @@
 	char* outputfilename=NULL;
 	char* lockname=NULL;
-	int   numberofnodes;
+	bool  qmu_analysis=false;
 	bool  waitonlock=false;
 
-	Model* model=NULL;
-
-	bool   qmu_analysis;
-
-	/*Results: */
-	Results* results=NULL;
-
-	Param*   param=NULL;
+	/*FemModel: */
+	FemModel* femmodel=NULL;
 
 	/*time*/
@@ -40,4 +34,7 @@
 	double   start_core, finish_core;
 	double   start_init, finish_init;
+
+	int analyses[1]={Prognostic2AnalysisEnum};
+	int solution_type=Prognosti2cAnalysisEnum;
 
 	MODULEBOOT();
@@ -62,16 +59,18 @@
 	/*Initialize model structure: */
 	MPI_Barrier(MPI_COMM_WORLD); start_init=MPI_Wtime();
-	model=new Model();
 
 	/*Open handle to data on disk: */
 	fid=pfopen(inputfilename,"rb");
 
-	_printf_("read and create finite element model:\n");
-	model->AddFormulation(fid,Prognostic2AnalysisEnum);
+	_printf_("create finite element model:\n");
+	femmodel=new FemModel(fid,solution_type,analyses,1);
+
+	/*add outputfilename in parameters: */
+	femmodel->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename));
 
 	/*recover parameters: */
 	model->FindParam(&waitonlock,WaitOnLockEnum);
 	model->FindParam(&qmu_analysis,QmuAnalysisEnum);
-	model->FindParam(&numberofnodes,NumberOfNodesEnum);
+
 
 	MPI_Barrier(MPI_COMM_WORLD); finish_init=MPI_Wtime();
@@ -83,6 +82,9 @@
 		_printf_("call computational core:\n");
 		MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
-		results=prognostic2_core(model);
+		prognostic2_core(femmodel);
 		MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
+
+		_printf_("write results to disk:\n");
+		OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters);
 
 	}
@@ -94,5 +96,5 @@
 		#ifdef _HAVE_DAKOTA_ 
 		MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
-		Qmux(model,Prognostic2AnalysisEnum,NoneAnalysisEnum);
+		Qmux(femmodel,Prognostic2AnalysisEnum,NoneAnalysisEnum);
 		MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
 	 	#else
@@ -100,7 +102,4 @@
 		#endif
 	}
-
-	_printf_("write results to disk:\n");
-	OutputResults(results,outputfilename);
 
 	if (waitonlock>0){
@@ -110,6 +109,5 @@
 
 	/*Free ressources:*/
-	delete results;
-	delete model;
+	delete femmodel;
 
 	/*Get finish time and close*/
Index: issm/trunk/src/c/solutions/prognostic2_core.cpp
===================================================================
--- issm/trunk/src/c/solutions/prognostic2_core.cpp	(revision 4054)
+++ issm/trunk/src/c/solutions/prognostic2_core.cpp	(revision 4055)
@@ -11,64 +11,29 @@
 #include "../solvers/solvers.h"
 
-Results* prognostic2_core(Model* model){
-
-	extern int my_rank;
-
-	/*output: */
-	Results* results=NULL;
-	Result* result=NULL;
-
-	/*intermediary: */
-	Vec vx_g=NULL;
-	Vec vy_g=NULL;
-
-	/*solutions: */
-	Vec h_g=NULL;
+void prognostic2_core(FemModel* femmodel){
 
 	/*flags: */
 	int verbose=0;
-	int numberofdofspernode;
-	int numberofnodes;
-	int numberofvertices;
-	int dofs[1]={1};
 
-	/*fem prognostic model: */
-	FemModel* fem_p=NULL;
+	/*activate formulation: */
+	femmodel->SetCurrentAnalysis(Prognostic2AnalysisEnum);
 
-	//initialize results:
-	results=new Results();
-
-	/*recover fem model: */
-	fem_p=model->GetFormulation(Prognostic2AnalysisEnum);
-
-	//first recover parameters common to all solutions
-	model->FindParam(&verbose,VerboseEnum);
-	model->FindParam(&numberofnodes,NumberOfNodesEnum);
-	model->FindParam(&numberofvertices,NumberOfVerticesEnum);
-	model->FindParam(&numberofdofspernode,NumberOfNodesEnum);
+	/*recover parameters: */
+	femmodel->parameters->FindParam(&verbose,VerboseEnum);
 
 	_printf_("depth averaging velocity...\n");
 	/*Where is it done?*/
+
+	_printf_("call computational core:\n");
+	solver_linear(NULL,femmodel);
 	
-	_printf_("call computational core:\n");
-	diagnostic_core_linear(&h_g,fem_p,Prognostic2AnalysisEnum,NoneAnalysisEnum);
-
 	_printf_("Averaging over vertices:\n");
-	FieldAverageOntoVerticesx(&h_g,fem_p->elements,fem_p->nodes,fem_p->vertices,fem_p->loads,fem_p->materials,fem_p->parameters);
+	//FieldAverageOntoVerticesx(&h_g,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
 
 	//_printf_("extrude computed thickness on all layers:\n");
 	//FieldExtrudex(h_g, fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,"thickness",0);
 
-	/*Plug results into output dataset: */
-	InputToResultx(&result,fem_p->elements,fem_p->nodes,fem_p->vertices, fem_p->loads, fem_p->materials,fem_p->parameters,ThicknessEnum,results->Size()+1,0,1); results->AddObject(result);
-	results->AddObject(new StringResult(results->Size()+1,AnalysisTypeEnum,0,1,EnumAsString(Prognostic2AnalysisEnum)));
-
-	/*Free ressources:*/
-	VecFree(&vx_g);
-	VecFree(&vy_g);
-	VecFree(&h_g);
-	
-	/*return: */
-	return results;
+	if(verbose)_printf_("saving results:\n");
+	InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum);
 
 }
Index: issm/trunk/src/c/solutions/prognostic_core.cpp
===================================================================
--- issm/trunk/src/c/solutions/prognostic_core.cpp	(revision 4054)
+++ issm/trunk/src/c/solutions/prognostic_core.cpp	(revision 4055)
@@ -11,52 +11,28 @@
 #include "../solvers/solvers.h"
 
-Results* prognostic_core(Model* model){
+void prognostic_core(FemModel* femmodel){
 
-	extern int my_rank;
-
-	/*output: */
-	Results* results=NULL;
-	Result*  result=NULL;
-
-	/*solutions: */
-	Vec h_g=NULL;
-
-	/*flags: */
+	/*parameters: */
 	int verbose=0;
 
-	/*fem prognostic model: */
-	FemModel* fem_p=NULL;
-
-	//initialize results
-	results=new Results();
-
-	/*recover fem model: */
-	fem_p=model->GetFormulation(PrognosticAnalysisEnum); ISSMASSERT(fem_p);
-
-	//first recover parameters common to all solutions
-	model->FindParam(&verbose,VerboseEnum);
+	/*activate formulation: */
+	femmodel->SetCurrentAnalysis(PrognosticAnalysisEnum);
+	
+	/*recover parameters: */
+	femmodel->parameters->FindParam(&verbose,VerboseEnum);
 
 	_printf_("depth averaging velocity...\n");
-	DepthAverageInputx(fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,VxEnum);
-	DepthAverageInputx(fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,VyEnum);
+	DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VxEnum,VxAverageEnum);
+	DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,VyEnum,VyAverageEnum);
 	
 	_printf_("call computational core:\n");
-	diagnostic_core_linear(&h_g,fem_p,PrognosticAnalysisEnum,NoneAnalysisEnum);
+	solver_linear(NULL,femmodel);
 		
-	_printf_("update inputs:\n");
-	UpdateInputsFromSolutionx( fem_p->elements,fem_p->nodes, fem_p->vertices, fem_p->loads, fem_p->materials, fem_p->parameters,h_g,PrognosticAnalysisEnum,NoneAnalysisEnum);
+	_printf_("extrude computed thickness on all layers:\n");
+	ExtrudeInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,ThicknessEnum);
 
-	_printf_("extrude computed thickness on all layers:\n");
-	ExtrudeInputx(fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,ThicknessEnum);
+	if(verbose)_printf_("saving results:\n");
+	InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum);
 	
-	_printf_("extract result from extruded inputs: \n");
-	InputToResultx(&result,fem_p->elements,fem_p->nodes,fem_p->vertices, fem_p->loads, fem_p->materials,fem_p->parameters,ThicknessEnum,results->Size()+1,0,1); results->AddObject(result);
-	results->AddObject(new StringResult(results->Size()+1,AnalysisTypeEnum,0,1,EnumAsString(PrognosticAnalysisEnum)));
+}
 
-	/*Free ressources:*/
-	VecFree(&h_g);
-	
-	//return: 
-	return results;
-
-}
Index: issm/trunk/src/c/solutions/solutions.h
===================================================================
--- issm/trunk/src/c/solutions/solutions.h	(revision 4054)
+++ issm/trunk/src/c/solutions/solutions.h	(revision 4055)
@@ -13,44 +13,37 @@
 
 /*cores: */
-Results* prognostic_core(FemModel* femmodel);
-Results* prognostic2_core(FemModel* femmodel);
-Results* balancedthickness_core(FemModel* femmodel);
-Results* balancedthickness2_core(FemModel* femmodel);
-Results* balancedvelocities_core(FemModel* femmodel);
-Results* slopecompute_core(FemModel* femmodel);
-Results* steadystate_core(FemModel* femmodel);
-Results* transient_core(FemModel* femmodel);
-Results* transient_core_2d(FemModel* femmodel);
-Results* transient_core_3d(FemModel* femmodel);
 void adjoint_core(FemModel* femmodel);
 void gradient_core(FemModel* femmodel,int step=0, double search_scalar=0);
 void diagnostic_core(FemModel* femmodel);
 void thermal_core(FemModel* femmodel);
-void surfaceslope_core(FemModel* femfemmodel);
-void bedslope_core(FemModel* femfemmodel);
+void thermal_core_step(FemModel* femmodel,int step, double time);
+void surfaceslope_core(FemModel* femmodel);
+void bedslope_core(FemModel* femmodel);
 void control_core(FemModel* femmodel);
+void prognostic_core(FemModel* femmodel);
+void prognostic2_core(FemModel* femmodel);
+void balancedthickness_core(FemModel* femmodel);
+void balancedthickness2_core(FemModel* femmodel);
+void balancedvelocities_core(FemModel* femmodel);
+void slopecompute_core(FemModel* femmodel);
+void steadystate_core(FemModel* femmodel);
+void transient2d_core(FemModel* femmodel);
+void transient3d_core(FemModel* femmodel);
+double objectivefunctionC(double search_scalar,OptArgs* optargs);
 
-//int GradJOrth(WorkspaceParams* workspaceparams);
+//convergence:
 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* femfemmodel);
+//optimization
+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);
+int BrentSearch(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);
+int GradJSearch(double* search_vector,FemModel* femmodel,int step);
 
-int BrentSearch(double* psearch_scalar,double* pJ,double xa, double xb, double tolerance, int maxiter, double fit,double optscal,double (*f)(double*,double,double,FemModel*),FemModel* femfemmodel);
-	
-double objectivefunctionC(double search_scalar,OptArgs* optargs);
-
-int GradJSearch(double* search_vector,FemModel* femfemmodel,int step);
-//int GradJCheck(WorkspaceParams* workspaceparams,int step,int status);
-
-//int ParameterUpdate(double* search_vector,int step, WorkspaceParams* workspaceparams,BatchParams* batchparams);
+//diverse
 void WriteLockFile(char* filename);
-
-void stokescontrolinit(FemModel* femfemmodel);
-void controlrestart(FemModel* femfemmodel);
-
-void CreateFemModel(FemModel* femfemmodel,ConstDataHandle MODEL,int analysis_type,int sub_analysis_type);
-//int BatchDebug(Mat* Kgg,Vec* pg,FemModel* femfemmodel,char* filename);
-void ResetBoundaryConditions(FemModel* femfemmodel, int analysis_type, int sub_analysis_type);
+void stokescontrolinit(FemModel* femmodel);
+void controlrestart(FemModel* femmodel,double* J);
+void ResetBoundaryConditions(FemModel* femmodel, int analysis_type, int sub_analysis_type);
 
 #endif
Index: issm/trunk/src/c/solutions/steadystate.cpp
===================================================================
--- issm/trunk/src/c/solutions/steadystate.cpp	(revision 4054)
+++ issm/trunk/src/c/solutions/steadystate.cpp	(revision 4055)
@@ -24,25 +24,10 @@
 	char* outputfilename=NULL;
 	char* lockname=NULL;
-	int   numberofnodes;
 	bool  qmu_analysis=false;
 	bool  control_analysis=false;
-	char* control_type=NULL;
+	bool waitonlock=false;
 
-	/*Model: */
-	Model* model=NULL;
-	FemModel* fem_dh=NULL;
-	FemModel* fem_ds=NULL;
-
-	/*Results: */
-	Results* results=NULL;
-	
-	bool waitonlock=false;
-	
-	double* u_g_initial=NULL;
-	double* p_g_initial=NULL;
-	double* u_g_obs=NULL;
-	double* weights=NULL;
-	double  dt;
-	BoolParam*  param=NULL;
+	/*FemModel: */
+	FemModel* femmodel=NULL;
 
 	/*time*/
@@ -50,4 +35,10 @@
 	double   start_core, finish_core;
 	double   start_init, finish_init;
+
+	/*intermediary: */
+	BoolParam* param=NULL;
+
+	int analyses[7]={DiagnosticHorizAnalysisEnum,DiagnosticVertAnalysisEnum,DiagnosticStokesAnalysisEnum,DiagnosticHutterAnalysisEnum,SlopeAnalysisEnum,ThermalAnalysisEnum,MeltingAnalysisEnum};
+	int solution_type=SteadystateAnalysisEnum;
 
 	MODULEBOOT();
@@ -70,32 +61,15 @@
 	lockname=argv[4];
 
+	/*Initialize model structure: */
+	MPI_Barrier(MPI_COMM_WORLD); start_init=MPI_Wtime();
+
 	/*Open handle to data on disk: */
 	fid=pfopen(inputfilename,"rb");
 
-	/*Initialize model structure: */
-	MPI_Barrier(MPI_COMM_WORLD); start_init=MPI_Wtime();
-	model=new Model();
+	_printf_("create finite element model:\n");
+	femmodel=new FemModel(fid,solution_type,analyses,5);
 
-	_printf_("read and create finite element model:\n");
-	
-	_printf_("\n   reading diagnostic horiz model data:\n");
-	model->AddFormulation(fid,DiagnosticAnalysisEnum,HorizAnalysisEnum);
-
-	_printf_("\n   reading diagnostic vert model data:\n");
-	model->AddFormulation(fid,DiagnosticAnalysisEnum,VertAnalysisEnum);
-	
-	_printf_("\n   reading diagnostic stokes model data:\n");
-	model->AddFormulation(fid,DiagnosticAnalysisEnum,StokesAnalysisEnum);
-	
-	_printf_("\n   reading diagnostic hutter model data:\n");
-	model->AddFormulation(fid,DiagnosticAnalysisEnum,HutterAnalysisEnum);
-	
-	_printf_("\n   reading surface and bed slope computation model data:\n");
-	model->AddFormulation(fid,SlopecomputeAnalysisEnum);
-
-	_printf_("\n   read and create thermal finite element model:\n");
-	model->AddFormulation(fid,ThermalAnalysisEnum);
-	_printf_("\n   read and create melting finite element model:\n");
-	model->AddFormulation(fid,MeltingAnalysisEnum);
+	/*add outputfilename in parameters: */
+	femmodel->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename));
 
 	/*recover parameters: */
@@ -113,5 +87,5 @@
 			_printf_("call computational core:\n");
 			MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
-			results=steadystate_core(model);
+			steadystate_core(femmodel);
 			MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
 
@@ -119,20 +93,16 @@
 		else{
 			/*change control_steady to 1 to know we are doing steadystate*/
-			fem_dh=model->GetFormulation(DiagnosticAnalysisEnum,HorizAnalysisEnum);
-			fem_ds=model->GetFormulation(DiagnosticAnalysisEnum,StokesAnalysisEnum);
-
-			param=(BoolParam*)fem_dh->parameters->FindParamObject(ControlSteadyEnum); param->value=true;
-			param=(BoolParam*)fem_ds->parameters->FindParamObject(ControlSteadyEnum); param->value=true;
+			param=(BoolParam*)femmodel->parameters->FindParamObject(ControlSteadyEnum); param->value=true;
 
 			/*run control analysis: */
 			_printf_("call computational core:\n");
 			MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
-			results=control_core(model);
+			control_core(femmodel);
 			MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
-
 		}
 
 		_printf_("write results to disk:\n");
-		OutputResults(results,outputfilename);
+		OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters);
+
 	}
 	else{
@@ -142,5 +112,5 @@
 		#ifdef _HAVE_DAKOTA_ 
 		MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
-		Qmux(model,SteadystateAnalysisEnum,NoneAnalysisEnum);
+		Qmux(femmodel,SteadystateAnalysisEnum,NoneAnalysisEnum);
 		MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
 	 	#else
@@ -155,9 +125,5 @@
 	
 	/*Free ressources */
-	xfree((void**)&u_g_initial);
-	xfree((void**)&u_g_obs);
-	xfree((void**)&p_g_initial);
-	delete model;
-	delete results;
+	delete femmodel;
 
 	/*Get finish time and close*/
Index: issm/trunk/src/c/solutions/steadystate_core.cpp
===================================================================
--- issm/trunk/src/c/solutions/steadystate_core.cpp	(revision 4054)
+++ issm/trunk/src/c/solutions/steadystate_core.cpp	(revision 4055)
@@ -9,136 +9,54 @@
 #include "./solutions.h"
 #include "../modules/modules.h"
+#include "../include/include.h"
 #include "../solvers/solvers.h"
 
-Results* steadystate_core(Model* model){
+void steadystate_core(FemModel* femmodel){
 
-	extern int my_rank;
+	/*intermediary: */
+	int step; 
 
-	/*fem models: */
-	FemModel* fem_dh=NULL;
-	FemModel* fem_dv=NULL;
-	FemModel* fem_dhu=NULL;
-	FemModel* fem_ds=NULL;
-	FemModel* fem_sl=NULL;
-	FemModel* fem_t=NULL;
-	FemModel* fem_m=NULL;
+	/*parameters: */
+	int verbose;
+	int dim;
+	
+	/* recover parameters:*/
+	femmodel->parameters->FindParam(&verbose,VerboseEnum); 
+	femmodel->parameters->FindParam(&dim,DimEnum);
 
-	/*output: */
-	Results* results=NULL;
-	Results* results_thermal=NULL;
-	Results* results_diagnostic=NULL;
-
-	/*solutions: */
-	Vec u_g=NULL;
-	Vec old_u_g=NULL;
-	Vec t_g=NULL;
-	Vec t_g_average=NULL;
-	Vec old_t_g=NULL;
-	Vec p_g=NULL;
-	Vec m_g=NULL;
-	Vec du_g=NULL;
-	Vec dt_g=NULL;
-	double ndu,nu;
-	double normdt,normt;
-	double eps_rel;
-
-	/*flags: */
-	int verbose=0;
-	int isstokes=0;
-	int numberofnodes;
-	int ndof;
-	int converged;
-	int step;
-
-	/*recover fem models: */
-	fem_dh=model->GetFormulation(DiagnosticAnalysisEnum,HorizAnalysisEnum);
-	fem_dv=model->GetFormulation(DiagnosticAnalysisEnum,VertAnalysisEnum);
-	fem_ds=model->GetFormulation(DiagnosticAnalysisEnum,StokesAnalysisEnum);
-	fem_dhu=model->GetFormulation(DiagnosticAnalysisEnum,HutterAnalysisEnum);
-	fem_sl=model->GetFormulation(SlopecomputeAnalysisEnum);
-	fem_t=model->GetFormulation(ThermalAnalysisEnum);
-	fem_m=model->GetFormulation(MeltingAnalysisEnum);
-
-
-	//first recover parameters common to all solutions
-	model->FindParam(&verbose,VerboseEnum);
-	model->FindParam(&numberofnodes,NumberOfNodesEnum);
-	model->FindParam(&eps_rel,EpsResEnum);
-	model->FindParam(&isstokes,IsStokesEnum);
-
-	//initialize: 
-	converged=0;
+	/*intialize counters: */
 	step=1;
-
-	if (isstokes)ndof=4;
-	else ndof=3;
 
 	for(;;){
 	
 		if(verbose)_printf_("%s%i\n","   computing temperature and velocity for step: ",step);
+		thermal_core(femmodel);
 
-		//first compute temperature at steady state.
-		results_thermal=thermal_core(model);
+		if(verbose)_printf_("%s\n","computing depth average temperature");
+		DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,TemperatureEnum,TemperatureAverageEnum);
 	
-		//get t_g and m_g;
-		VecFree(&t_g);results_thermal->FindResult(&t_g,"t_g");
-		VecFree(&m_g);results_thermal->FindResult(&m_g,"m_g");
-		delete results_thermal;
+		if(verbose)_printf_("%s\n","computing new velocity");
+		diagnostic_core(femmodel);
 
-		//Add temperature to inputs.
-		//compute depth averaged temperature and add to inputs
-		VecDuplicatePatch(&t_g_average,t_g); 
-		FieldDepthAveragex( t_g_average, fem_t->elements,fem_t->nodes, fem_t->vertices,fem_t->loads, fem_t->materials,fem_t->parameters,"temperature");
-		model->UpdateInputsFromVector(t_g_average,TemperatureAverageEnum,VertexEnum);
-		model->UpdateInputsFromVector(t_g,TemperatureEnum,VertexEnum);
-		VecFree(&t_g_average); //not needed anymore
-
-		//now compute diagnostic velocity using the steady state temperature.
-		results_diagnostic=diagnostic_core(model);
-
-		//get p_g and u_g
-		VecFree(&u_g);results_diagnostic->FindResult(&u_g,"u_g");
-		VecFree(&p_g);results_diagnostic->FindResult(&p_g,"p_g");
-		delete results_diagnostic;
-
-		//convergence? 
-		if(step>1){
-			VecDuplicatePatch(&du_g,old_u_g);VecAYPX(du_g,-1.0,u_g);
-			VecNorm(du_g,NORM_2,&ndu); VecNorm(old_u_g,NORM_2,&nu); VecFree(&du_g);
-
-			VecDuplicatePatch(&dt_g,old_t_g); VecAYPX(dt_g,-1.0,t_g);
-			VecNorm(dt_g,NORM_2,&normdt); VecNorm(old_t_g,NORM_2,&normt);VecFree(&dt_g);
-					
-			if (verbose) _printf_("%-60s%g\n                                     %s%g\n                                     %s%g%s\n",
-					  "      relative convergence criterion: velocity -> norm(du)/norm(u)=   ",ndu/nu*100," temperature -> norm(dt)/norm(t)=",normdt/normt*100," eps_rel:                        ",eps_rel*100," %");
+		if(verbose)_printf_("%s\n","checking velocity, temperature and pressure convergence");
+		if (step>1) if(steadystateconvergence(femmodel))goto cleanup_and_return;
 		
-			if ((ndu/nu<=eps_rel)  && (normdt/normt<=eps_rel)) converged=1;
-			else converged=0;
-		}
-		else{
-			converged=0;
-		}
-
-		VecFree(&old_u_g);VecDuplicatePatch(&old_u_g,u_g);
-		VecFree(&old_t_g);VecDuplicatePatch(&old_t_g,t_g);
-
+		if(verbose)_printf_("%s\n","saving velocity, temperature and pressure to check for convergence at next step");
+		DuplicateInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum,OldVxEnum);
+		DuplicateInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum,OldVyEnum);
+		if(dim==3)DuplicateInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum,OldVzEnum);
+		DuplicateInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,OldPressureEnum);
+		DuplicateInputx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,OldTemperatureEnum);
+		
+		//increase counter
 		step++;
-		if (converged)break;
 	}
 
-	/*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,"p_g",p_g));
-	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));
-	results->AddObject(new StringResult(results->Size()+1,AnalysisTypeEnum,0,1,EnumAsString(SteadyAnalysisEnum)));
-
-
-	/*Free ressource*/
-	VecFree(&old_u_g);
-	VecFree(&old_t_g);
-	VecFree(&u_g);
-	VecFree(&p_g);
-	VecFree(&t_g);
-	VecFree(&m_g);
+	cleanup_and_return:
+	
+	if(verbose)_printf_("saving results:\n");
+	InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum);
+	InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum);
+	InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum);
+	if(dim==3) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum);
 }
Index: issm/trunk/src/c/solutions/steadystateconvergence.cpp
===================================================================
--- issm/trunk/src/c/solutions/steadystateconvergence.cpp	(revision 4055)
+++ issm/trunk/src/c/solutions/steadystateconvergence.cpp	(revision 4055)
@@ -0,0 +1,36 @@
+/*!\file: steadystateconvergence.cpp
+ * \brief: determine convergence of steady state solution
+ */ 
+
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+int steadystateconvergence(FemModel* femmodel){
+	
+	/*output: */
+	int converged=0;
+	int velocity_converged=0;
+	int temperature_converged=0;
+
+	/*intermediary: */
+	int velocityenums[8]={VxEnum,OldVxEnum,VyEnum,OldVyEnum,VzEnum,OldVzEnum,PressureEnum,OldPressureEnum}; //pairs of enums (new and old) on which to carry out the converence tests
+	int temperatureenums[2]={TemperatureEnum,OldTemperatureEnum};
+	int convergencecriterion[1]={RelativeEnum}; //criterions for convergence, RelativeEnum or AbsoluteEnum 
+	int convergencecriterionvalue[1]; //value of criterion to be respected
+
+	/*retrieve parameters: */
+	femmodel->parameters->FindParam(&convergencecriterionvalue[0],EpsRelEnum);
+	
+	/*figure out convergence at the input level, because we don't have the solution vectors!: */
+	velocity_converged=InputConvergencex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,velocityenums,8,convergencecriterion,convergencecriterionvalue,1);
+	temperature_converged=InputConvergencex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,temperatureenums,2,convergencecriterion,convergencecriterionvalue,1);
+
+	if(velocity_converged && temperature_converged)converged=1;
+
+	/*return: */
+	return converged;
+
+}
Index: issm/trunk/src/c/solutions/stokescontrolinit.cpp
===================================================================
--- issm/trunk/src/c/solutions/stokescontrolinit.cpp	(revision 4054)
+++ issm/trunk/src/c/solutions/stokescontrolinit.cpp	(revision 4055)
@@ -18,15 +18,16 @@
 	int isstokes=0;
 	double stokesreconditioning;
-
+	bool conserve_loads=true;
+	
 	/*first recover parameters common to all solutions:*/
-	model->FindParam(&verbose,VerboseEnum);
-	model->FindParam(&isstokes,IsStokesEnum);
-	model->FindParam(&stokesreconditioning,StokesReconditioningEnum);
+	femmodel->parameters->FindParam(&verbose,VerboseEnum);
+	femmodel->parameters->FindParam(&isstokes,IsStokesEnum);
+	femmodel->parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum);
 
 	/*if no Stokes analysis carried out, assign output and return*/
-	if (!isstokes)femmodel->SetAnalysisType(DiagnosticAnalysisEnum,DiagnosticHorizAnalysisEnum);
+	if (!isstokes)femmodel->SetCurrentAnalysis(DiagnosticHorizAnalysisEnum);
 
 	/* For Stokes inverse control method, we are going to carry out the inversion only on the Stokes part. So we need 
-	 * to solve the Hutter or MacAyeal/Pattyn model here, and constrain the Stokes model using the Hutter 
+	 * to solve the Hutter or MacAyeal/Pattyn femmodel here, and constrain the Stokes femmodel using the Hutter 
 	 * or MacAyeal/Pattyn at the boundary. We don't want to have to do that at every inversion step, as 
 	 * it needs be done only once: */
@@ -36,8 +37,9 @@
 	
 	/*Run a complete diagnostic to update the Stokes spcs: */
-	solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,false,DiagnosticAnalysisEnum,HorizAnalysisEnum); //false means we expose loads to changes inside the solution
+	femmodel->SetCurrentAnalysis(DiagnosticHorizAnalysisEnum);
+	solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,conserve_loads);
 
 	//vertical velocity
-	solver_linear(NULL,femmodel,DiagnosticAnalysisEnum,VertAnalysisEnum);
+	solver_linear(NULL,femmodel);
 
 	//recondition" pressure computed previously:
@@ -49,7 +51,5 @@
 
 	if(verbose)_printf_("%s\n"," computing stokes velocities and pressure ...");
-	solver_diagnostic_nonlinear(NULL,NULL,NULL,NULL,femmodel,DiagnosticAnalysisEnum,StokesAnalysisEnum);
-
-	/*Assign output*/
-	model->SetActiveFormulation(fem_ds);
+	femmodel->SetCurrentAnalysis(DiagnosticStokesAnalysisEnum);
+	solver_diagnostic_nonlinear(NULL,NULL,NULL,femmodel,conserve_loads);
 }
Index: issm/trunk/src/c/solutions/surfaceslope.cpp
===================================================================
--- issm/trunk/src/c/solutions/surfaceslope.cpp	(revision 4054)
+++ issm/trunk/src/c/solutions/surfaceslope.cpp	(revision 4055)
@@ -37,6 +37,6 @@
 	double   start_init, finish_init;
 
-	int analyses[1]={SurfaceSlopeComputeAnalysisEnum};
-	int solution_type=SurfaceSlopeComputeSolutionEnum;
+	int analyses[1]={SlopeAnalysisEnum};
+	int solution_type=SlopeAnalysisEnum;
 
 	MODULEBOOT();
@@ -67,4 +67,7 @@
 	_printf_("create finite element model, using analyses types statically defined above:\n");
 	femmodel=new FemModel(fid,solution_type,analyses,1);
+
+	/*add outputfilename in parameters: */
+	femmodel->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename));
 	
 	/*get parameters: */
@@ -83,5 +86,5 @@
 
 		_printf_("write results to disk:\n");
-		OutputResults(femmodel,outputfilename,DiagnosticAnalysisEnum);
+		OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters);
 	}
 	else{
@@ -92,5 +95,5 @@
 		#ifdef _HAVE_DAKOTA_ 
 		MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
-		Qmux(femmodel,SlopeComputeAnalysisEnum,NoneAnalysisEnum);
+		Qmux(femmodel,SurfaceSlopeAnalysisEnum,NoneAnalysisEnum);
 		MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
 	 	#else
Index: issm/trunk/src/c/solutions/surfaceslope_core.cpp
===================================================================
--- issm/trunk/src/c/solutions/surfaceslope_core.cpp	(revision 4054)
+++ issm/trunk/src/c/solutions/surfaceslope_core.cpp	(revision 4055)
@@ -25,12 +25,19 @@
 
 	/*Call on core computations: */
-	solver_linear(NULL,femmodel,SlopeComputeAnalysisEnum,SurfaceSlopeXAnalysisEnum);
-	solver_linear(NULL,femmodel,SlopeComputeAnalysisEnum,SurfaceSlopeYAnalysisEnum);
+	femmodel->SetCurrentAnalysisAlias(SlopeAnalysisEnum,SurfaceSlopeXAnalysisEnum);
+	solver_linear(NULL,femmodel);
+	femmodel->SetCurrentAnalysisAlias(SlopeAnalysisEnum,SurfaceSlopeYAnalysisEnum);
+	solver_linear(NULL,femmodel);
 	
 	/*extrude inputs if we are in 3D: */
 	if (dim==3){
 		if(verbose)_printf_("%s\n","extruding slope in 3d...");
-		InputsExtrudex( femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceSlopeXEnum,0);
-		InputsExtrudex( femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceSlopeYEnum,0);
+		ExtrudeInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,SurfaceSlopeXEnum);
+		ExtrudeInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,SurfaceSlopeYEnum);
 	}
+	
+	if(verbose)_printf_("saving results:\n");
+	InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceSlopeXEnum);
+	InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceSlopeYEnum);
+
 }
Index: issm/trunk/src/c/solutions/thermal.cpp
===================================================================
--- issm/trunk/src/c/solutions/thermal.cpp	(revision 4054)
+++ issm/trunk/src/c/solutions/thermal.cpp	(revision 4055)
@@ -68,4 +68,7 @@
 	femmodel=new FemModel(fid,solution_type,analyses,2);
 
+	/*add outputfilename in parameters: */
+	femmodel->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename));
+
 	/*get parameters: */
 	femmodel->parameters->FindParam(&qmu_analysis,QmuAnalysisEnum);
@@ -85,5 +88,5 @@
 		
 		_printf_("write results to disk:\n");
-		OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters, outputfilename);
+		OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters);
 
 	}
Index: issm/trunk/src/c/solutions/thermal_core.cpp
===================================================================
--- issm/trunk/src/c/solutions/thermal_core.cpp	(revision 4054)
+++ issm/trunk/src/c/solutions/thermal_core.cpp	(revision 4055)
@@ -42,9 +42,6 @@
 		time=(i+1)*dt;
 
-		if(verbose)_printf_("computing temperatures:\n");
-		solver_thermal_nonlinear(NULL,NULL,femmodel,ThermalAnalysisEnum,NoneAnalysisEnum);
-
-		if(verbose)_printf_("computing melting:\n");
-		solver_linear(NULL,femmodel,MeltingAnalysisEnum,NoneAnalysisEnum);
+		/*call thermal_core_step: */
+		thermal_core_step(femmodel,i,time);
 
 		if(verbose)_printf_("saving results:\n");
Index: issm/trunk/src/c/solutions/thermal_core_step.cpp
===================================================================
--- issm/trunk/src/c/solutions/thermal_core_step.cpp	(revision 4055)
+++ issm/trunk/src/c/solutions/thermal_core_step.cpp	(revision 4055)
@@ -0,0 +1,26 @@
+/*!\file: thermal_core_step.cpp
+ * \brief: core of the thermal solution 
+ */ 
+
+#include "../toolkits/toolkits.h"
+#include "../objects/objects.h"
+#include "../shared/shared.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+#include "./solutions.h"
+#include "../modules/modules.h"
+#include "../include/include.h"
+#include "../solvers/solvers.h"
+
+void thermal_core_step(FemModel* femmodel,int step, double time){
+
+	//first recover parameters common to all solutions
+	femmodel->parameters->FindParam(&verbose,VerboseEnum);
+
+	if(verbose)_printf_("computing temperatures:\n");
+	femmodel->SetCurrentAnalysis(ThermalAnalysisEnum);
+	solver_thermal_nonlinear(NULL,NULL,femmodel);
+
+	if(verbose)_printf_("computing melting:\n");
+	femmodel->SetCurrentAnalysis(MeltingAnalysisEnum);
+	solver_linear(NULL,femmodel);
+}
Index: issm/trunk/src/c/solutions/transient.cpp
===================================================================
--- issm/trunk/src/c/solutions/transient.cpp	(revision 4054)
+++ 	(revision )
@@ -1,150 +1,0 @@
-/*!\file:  transient.cpp
- * \brief: transient solution
- */ 
-
-#ifdef HAVE_CONFIG_H
-	#include "config.h"
-#else
-#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
-#endif
-
-#include "../objects/objects.h"
-#include "../shared/shared.h"
-#include "../DataSet/DataSet.h"
-#include "../EnumDefinitions/EnumDefinitions.h"
-#include "../include/include.h"
-#include "../modules/modules.h"
-#include "./solutions.h"
-
-int main(int argc,char* *argv){
-	
-	/*I/O: */
-	FILE* fid=NULL;
-	char* inputfilename=NULL;
-	char* outputfilename=NULL;
-	char* lockname=NULL;
-	int   numberofnodes;
-	bool qmu_analysis=false;
-	bool waitonlock=false;
-
-	/*Model: */
-	Model* model=NULL;
-	int dim=-1;
-
-	/*Results: */
-	Results* results=NULL;
-	
-	Param*  param=NULL;
-
-	/*time*/
-	double   start, finish;
-	double   start_core, finish_core;
-	double   start_init, finish_init;
-
-	MODULEBOOT();
-
-	#if !defined(_PARALLEL_) || (defined(_PARALLEL_) && !defined(_HAVE_PETSC_))
-	ISSMERROR(" parallel executable was compiled without support of parallel libraries!");
-	#endif
-
-	/*Initialize Petsc and get start time*/
-	PetscInitialize(&argc,&argv,(char *)0,"");  
-	MPI_Barrier(MPI_COMM_WORLD); start=MPI_Wtime();
-
-	/*Size and rank: */
-	MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);  
-	MPI_Comm_size(MPI_COMM_WORLD,&num_procs); 
-
-	_printf_("recover , input file name and output file name:\n");
-	inputfilename=argv[2];
-	outputfilename=argv[3];
-	lockname=argv[4];
-
-	/*Open handle to data on disk: */
-	fid=pfopen(inputfilename,"rb");
-
-	/*Initialize model structure: */
-	MPI_Barrier(MPI_COMM_WORLD); start_init=MPI_Wtime();
-	model=new Model();
-
-	_printf_("read and create finite element model:\n");
-	_printf_("\n   reading diagnostic horiz model data:\n");
-	model->AddFormulation(fid,DiagnosticAnalysisEnum,HorizAnalysisEnum);
-
-	_printf_("\n   reading diagnostic vert model data:\n");
-	model->AddFormulation(fid,DiagnosticAnalysisEnum,VertAnalysisEnum);
-	
-	_printf_("\n   reading diagnostic stokes model data:\n");
-	model->AddFormulation(fid,DiagnosticAnalysisEnum,StokesAnalysisEnum);
-	
-	_printf_("\n   reading diagnostic hutter model data:\n");
-	model->AddFormulation(fid,DiagnosticAnalysisEnum,HutterAnalysisEnum);
-	
-	_printf_("\n   reading surface and bed slope computation model data:\n");
-	model->AddFormulation(fid,SlopecomputeAnalysisEnum);
-
-	_printf_("\n   reading prognositc model data:\n");
-	model->AddFormulation(fid,PrognosticAnalysisEnum);
-	
-	/*Do we run in 3d?, in which case we need thermal and melting also:*/
-	model->FindParam(&dim,DimEnum);
-	if(dim==3){
-		_printf_("read and create thermal finite element model:\n");
-		model->AddFormulation(fid,ThermalAnalysisEnum,TransientAnalysisEnum);
-		_printf_("read and create melting finite element model:\n");
-		model->AddFormulation(fid,MeltingAnalysisEnum,TransientAnalysisEnum);
-	}
-
-	/*recover parameters: */
-	model->FindParam(&waitonlock,WaitOnLockEnum);
-	model->FindParam(&qmu_analysis,QmuAnalysisEnum);
-	
-	MPI_Barrier(MPI_COMM_WORLD); finish_init=MPI_Wtime();
-
-	/*are we running the solution sequence, or a qmu wrapper around it? : */
-	if(!qmu_analysis){
-
-		/*run diagnostic analysis: */
-		_printf_("call computational core:\n");
-		MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
-		results=transient_core(model);
-		MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
-	
-		_printf_("write results to disk:\n");
-		OutputResults(results,outputfilename);
-	}
-	else{
-		/*run qmu analysis: */
-		_printf_("calling qmu analysis on transient core:\n");
-
-		#ifdef _HAVE_DAKOTA_ 
-		MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
-		Qmux(model,TransientAnalysisEnum,NoneAnalysisEnum);
-		MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
-	 	#else
-		ISSMERROR(" Dakota not present, cannot do qmu!");
-		#endif
-	}
-
-	if (waitonlock>0){
-		_printf_("write lock file:\n");
-		WriteLockFile(lockname);
-	}
-
-	/*Free ressources:*/
-	delete results;
-	delete model;
-
-	/*Get finish time and close*/
-	MPI_Barrier(MPI_COMM_WORLD); finish = MPI_Wtime( );
-	_printf_("\n   %-34s %f seconds  \n","Model initialization elapsed time:",finish_init-start_init);
-	_printf_("   %-34s %f seconds  \n","Core solution elapsed time:",finish_core-start_core);
-	_printf_("\n   %s %i hrs %i min %i sec\n\n","Total elapsed time:",int((finish-start)/3600),int(int(finish-start)%3600/60),int(finish-start)%60);
-	_printf_("closing MPI and Petsc\n");
-	PetscFinalize(); 
-	
-	/*end module: */
-	MODULEEND();
-	
-	return 0; //unix success return;
-}
Index: issm/trunk/src/c/solutions/transient2d.cpp
===================================================================
--- issm/trunk/src/c/solutions/transient2d.cpp	(revision 4055)
+++ issm/trunk/src/c/solutions/transient2d.cpp	(revision 4055)
@@ -0,0 +1,120 @@
+/*!\file:  transient2d.cpp
+ * \brief: transient2d solution
+ */ 
+
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "../objects/objects.h"
+#include "../shared/shared.h"
+#include "../DataSet/DataSet.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+#include "../include/include.h"
+#include "../modules/modules.h"
+#include "./solutions.h"
+
+int main(int argc,char* *argv){
+	
+	/*I/O: */
+	FILE* fid=NULL;
+	char* inputfilename=NULL;
+	char* outputfilename=NULL;
+	char* lockname=NULL;
+	bool qmu_analysis=false;
+	bool waitonlock=false;
+
+	/*FemModel: */
+	FemModel* femmodel=NULL;
+
+	/*time*/
+	double   start, finish;
+	double   start_core, finish_core;
+	double   start_init, finish_init;
+
+	int analyses[2]={DiagnosticHorizAnalysisEnum,PrognosticAnalysisEnum};
+	int solution_type=TransientAnalysisEnum;
+
+	MODULEBOOT();
+
+	#if !defined(_PARALLEL_) || (defined(_PARALLEL_) && !defined(_HAVE_PETSC_))
+	ISSMERROR(" parallel executable was compiled without support of parallel libraries!");
+	#endif
+
+	/*Initialize Petsc and get start time*/
+	PetscInitialize(&argc,&argv,(char *)0,"");  
+	MPI_Barrier(MPI_COMM_WORLD); start=MPI_Wtime();
+
+	/*Size and rank: */
+	MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);  
+	MPI_Comm_size(MPI_COMM_WORLD,&num_procs); 
+
+	_printf_("recover , input file name and output file name:\n");
+	inputfilename=argv[2];
+	outputfilename=argv[3];
+	lockname=argv[4];
+
+	/*Open handle to data on disk: */
+	fid=pfopen(inputfilename,"rb");
+
+	_printf_("create finite element model:\n");
+	femmodel=new FemModel(fid,solution_type,analyses,7);
+
+	/*add outputfilename in parameters: */
+	femmodel->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename));
+
+	/*get parameters: */
+	femmodel->parameters->FindParam(&qmu_analysis,QmuAnalysisEnum);
+	femmodel->parameters->FindParam(&control_analysis,ControlAnalysisEnum);
+	femmodel->parameters->FindParam(&waitonlock,WaitOnLockEnum);
+
+	MPI_Barrier(MPI_COMM_WORLD); finish_init=MPI_Wtime();
+
+	/*are we running the solution sequence, or a qmu wrapper around it? : */
+	if(!qmu_analysis){
+
+		/*run diagnostic analysis: */
+		_printf_("call computational core:\n");
+		MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
+		transient2d_core(model);
+		MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
+	
+		_printf_("write results to disk:\n");
+		OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters);
+	}
+	else{
+		/*run qmu analysis: */
+		_printf_("calling qmu analysis on transient2d core:\n");
+
+		#ifdef _HAVE_DAKOTA_ 
+		MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
+		Qmux(femmodel,Transient2DAnalysisEnum,NoneAnalysisEnum);
+		MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
+	 	#else
+		ISSMERROR(" Dakota not present, cannot do qmu!");
+		#endif
+	}
+
+	if (waitonlock>0){
+		_printf_("write lock file:\n");
+		WriteLockFile(lockname);
+	}
+
+	/*Free ressources:*/
+	delete femmodel;
+
+	/*Get finish time and close*/
+	MPI_Barrier(MPI_COMM_WORLD); finish = MPI_Wtime( );
+	_printf_("\n   %-34s %f seconds  \n","Model initialization elapsed time:",finish_init-start_init);
+	_printf_("   %-34s %f seconds  \n","Core solution elapsed time:",finish_core-start_core);
+	_printf_("\n   %s %i hrs %i min %i sec\n\n","Total elapsed time:",int((finish-start)/3600),int(int(finish-start)%3600/60),int(finish-start)%60);
+	_printf_("closing MPI and Petsc\n");
+	PetscFinalize(); 
+	
+	/*end module: */
+	MODULEEND();
+	
+	return 0; //unix success return;
+}
Index: issm/trunk/src/c/solutions/transient2d_core.cpp
===================================================================
--- issm/trunk/src/c/solutions/transient2d_core.cpp	(revision 4055)
+++ issm/trunk/src/c/solutions/transient2d_core.cpp	(revision 4055)
@@ -0,0 +1,59 @@
+/*!\file: transient_2d_core.cpp
+ * \brief: core of the transient_2d solution 
+ */ 
+
+#include "../toolkits/toolkits.h"
+#include "../objects/objects.h"
+#include "../shared/shared.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+#include "./solutions.h"
+#include "../modules/modules.h"
+
+void transient2d_core(FemModel* femmodel){
+
+	int i;
+
+	/*parameters: */
+	int verbose=0;
+	double finaltime;
+	double dt;
+	
+	/*intermediary: */
+	int    step;
+	double time;
+
+	/* recover parameters: */
+	femmodel->parameters->FindParam(&verbose,VerboseEnum);
+	femmodel->parameters->FindParam(&finaltime,NdtEnum);
+	femmodel->parameters->FindParam(&dt,DtEnum);
+
+	/*initialize: */
+	step=0;
+	time=0;
+
+	while(time<finaltime){ //make sure we run up to finaltime.
+	
+		_printf_("%s%g%s%i%s%g\n","time [yr]: ",time,"    iteration number: ",step,"/",floor(finaltime/dt));
+
+		step+=1;
+		time+=dt;
+
+		if(verbose)_printf_("%s\n","computing new velocity");
+		diagnostic_core(femmodel);
+
+		if(verbose)_printf_("%s\n","computing new thickness");
+		prognostic_core(femmodel);
+	
+		if(verbose)_printf_("%s\n","updating geometry");
+		UpdateGeometryx(femmodel->elements, femmodel->nodes,femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters); 
+
+		if(verbose)_printf_("%s\n","saving results:\n");
+		InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum,step,time); 
+		InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum,step,time);
+		InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,step,time);
+		InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum,step,time);
+		InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceEnum,step,time);
+		InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedEnum,step,time);
+	}
+
+}
Index: issm/trunk/src/c/solutions/transient3d.cpp
===================================================================
--- issm/trunk/src/c/solutions/transient3d.cpp	(revision 4055)
+++ issm/trunk/src/c/solutions/transient3d.cpp	(revision 4055)
@@ -0,0 +1,120 @@
+/*!\file:  transient.cpp
+ * \brief: transient solution
+ */ 
+
+#ifdef HAVE_CONFIG_H
+	#include "config.h"
+#else
+#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
+#endif
+
+#include "../objects/objects.h"
+#include "../shared/shared.h"
+#include "../DataSet/DataSet.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+#include "../include/include.h"
+#include "../modules/modules.h"
+#include "./solutions.h"
+
+int main(int argc,char* *argv){
+	
+	/*I/O: */
+	FILE* fid=NULL;
+	char* inputfilename=NULL;
+	char* outputfilename=NULL;
+	char* lockname=NULL;
+	bool qmu_analysis=false;
+	bool waitonlock=false;
+
+	/*FemModel: */
+	FemModel* femmodel=NULL;
+
+	/*time*/
+	double   start, finish;
+	double   start_core, finish_core;
+	double   start_init, finish_init;
+
+	int analyses[8]={DiagnosticHorizAnalysisEnum,DiagnosticVertAnalysisEnum,DiagnosticStokesAnalysisEnum,DiagnosticHutterAnalysisEnum,SlopeAnalysisEnum,PrognosticAnalysisEnum,ThermalAnalysisEnum,MeltingAnalysisEnum};
+	int solution_type=TransientAnalysisEnum;
+
+	MODULEBOOT();
+
+	#if !defined(_PARALLEL_) || (defined(_PARALLEL_) && !defined(_HAVE_PETSC_))
+	ISSMERROR(" parallel executable was compiled without support of parallel libraries!");
+	#endif
+
+	/*Initialize Petsc and get start time*/
+	PetscInitialize(&argc,&argv,(char *)0,"");  
+	MPI_Barrier(MPI_COMM_WORLD); start=MPI_Wtime();
+
+	/*Size and rank: */
+	MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);  
+	MPI_Comm_size(MPI_COMM_WORLD,&num_procs); 
+
+	_printf_("recover , input file name and output file name:\n");
+	inputfilename=argv[2];
+	outputfilename=argv[3];
+	lockname=argv[4];
+
+	/*Open handle to data on disk: */
+	fid=pfopen(inputfilename,"rb");
+
+	_printf_("create finite element model:\n");
+	femmodel=new FemModel(fid,solution_type,analyses,8);
+
+	/*add outputfilename in parameters: */
+	femmodel->parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename));
+
+	/*get parameters: */
+	femmodel->parameters->FindParam(&qmu_analysis,QmuAnalysisEnum);
+	femmodel->parameters->FindParam(&control_analysis,ControlAnalysisEnum);
+	femmodel->parameters->FindParam(&waitonlock,WaitOnLockEnum);
+
+	MPI_Barrier(MPI_COMM_WORLD); finish_init=MPI_Wtime();
+
+	/*are we running the solution sequence, or a qmu wrapper around it? : */
+	if(!qmu_analysis){
+
+		/*run diagnostic analysis: */
+		_printf_("call computational core:\n");
+		MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
+		transient3d_core(model);
+		MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
+	
+		_printf_("write results to disk:\n");
+		OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters);
+	}
+	else{
+		/*run qmu analysis: */
+		_printf_("calling qmu analysis on transient core:\n");
+
+		#ifdef _HAVE_DAKOTA_ 
+		MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
+		Qmux(femmodel,Transient3DAnalysisEnum,NoneAnalysisEnum);
+		MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
+	 	#else
+		ISSMERROR(" Dakota not present, cannot do qmu!");
+		#endif
+	}
+
+	if (waitonlock>0){
+		_printf_("write lock file:\n");
+		WriteLockFile(lockname);
+	}
+
+	/*Free ressources:*/
+	delete femmodel;
+
+	/*Get finish time and close*/
+	MPI_Barrier(MPI_COMM_WORLD); finish = MPI_Wtime( );
+	_printf_("\n   %-34s %f seconds  \n","Model initialization elapsed time:",finish_init-start_init);
+	_printf_("   %-34s %f seconds  \n","Core solution elapsed time:",finish_core-start_core);
+	_printf_("\n   %s %i hrs %i min %i sec\n\n","Total elapsed time:",int((finish-start)/3600),int(int(finish-start)%3600/60),int(finish-start)%60);
+	_printf_("closing MPI and Petsc\n");
+	PetscFinalize(); 
+	
+	/*end module: */
+	MODULEEND();
+	
+	return 0; //unix success return;
+}
Index: issm/trunk/src/c/solutions/transient3d_core.cpp
===================================================================
--- issm/trunk/src/c/solutions/transient3d_core.cpp	(revision 4055)
+++ issm/trunk/src/c/solutions/transient3d_core.cpp	(revision 4055)
@@ -0,0 +1,75 @@
+/*!\file: transient_3d_core.cpp
+ * \brief: core of the transient_3d solution 
+ */ 
+
+#include "../toolkits/toolkits.h"
+#include "../objects/objects.h"
+#include "../shared/shared.h"
+#include "../EnumDefinitions/EnumDefinitions.h"
+#include "./solutions.h"
+#include "../modules/modules.h"
+#include "../include/include.h"
+#include "../solvers/solvers.h"
+
+void transient3d_core(FemModel* femmodel){
+
+	/*parameters: */
+	int verbose=0;
+	double finaltime;
+	double dt;
+	
+	/*intermediary: */
+	int    step;
+	double time;
+
+	//first recover parameters common to all solutions
+	femmodel->parameters->FindParam(&verbose,VerboseEnum);
+	femmodel->parameters->FindParam(&finaltime,NdtEnum);
+	femmodel->parameters->FindParam(&dt,DtEnum);
+
+	/*initialize: */
+	step=0;
+	time=0;
+
+	while(time<finaltime){ //make sure we run up to finaltime.
+	
+		if(verbose)_printf_("%s%g%s%i%s%g\n","time [yr]: ",time,"    iteration number: ",step,"/",floor(finaltime/dt));
+
+		step+=1;
+		time+=dt;
+
+		if(verbose)_printf_("computing temperatures:\n");
+		thermal_core_step(femmodel,step,time);
+
+		if(verbose)_printf_("%s\n","computing depth average temperature");
+		DepthAverageInputx(femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters,TemperatureEnum,TemperatureAverageEnum);
+
+		if(verbose)_printf_("%s\n","computing new velocity");
+		diagnostic_core(femmodel);
+
+		if(verbose)_printf_("%s\n","computing new thickness");
+		prognostic_core(femmodel);
+	
+		if(verbose)_printf_("   updating geometry\n");
+		UpdateGeometryx(femmodel->elements, femmodel->nodes,femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters); 
+		
+		if(verbose)_printf_("%s\n","updating vertices positions");
+		UpdateVertexPositionsx(femmodel->elements, femmodel->nodes,femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters); 
+
+		if(verbose)_printf_("%s\n","saving results:\n");
+		InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum,step,time);
+		InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum,step,time);
+		InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,step,time);
+		InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum,step,time);
+		InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceEnum,step,time);
+		InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedEnum,step,time);
+		InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,step,time);
+		InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,MeltingRateEnum,step,time);
+
+		if (step%5==0){
+			if(verbose)_printf_("%s","      saving temporary results...");
+			OutputResults(femmodel->elements, femmodel->loads, femmodel->nodes, femmodel->vertices, femmodel->materials, femmodel->parameters);
+		}
+	}
+
+}
Index: issm/trunk/src/c/solutions/transient_core.cpp
===================================================================
--- issm/trunk/src/c/solutions/transient_core.cpp	(revision 4054)
+++ 	(revision )
@@ -1,30 +1,0 @@
-/*!\file: transient_core.cpp
- * \brief: core of the transient solution 
- */ 
-
-#include "../toolkits/toolkits.h"
-#include "../objects/objects.h"
-#include "../shared/shared.h"
-#include "../EnumDefinitions/EnumDefinitions.h"
-#include "./solutions.h"
-#include "../modules/modules.h"
-#include "../solvers/solvers.h"
-
-Results* transient_core(Model* model){
-
-	int dim=-1;
-	Results* results=NULL;
-
-	//first recover parameters common to all solutions
-	model->FindParam(&dim,DimEnum);
-
-	//branch out 
-	if(dim==2){
-		results=transient_core_2d(model);
-	}
-	else if(dim==3){
-		results=transient_core_3d(model);
-	}
-	else ISSMERROR("%s%i%s"," dimension ",dim," not supported yet!");
-
-}
Index: issm/trunk/src/c/solutions/transient_core_2d.cpp
===================================================================
--- issm/trunk/src/c/solutions/transient_core_2d.cpp	(revision 4054)
+++ 	(revision )
@@ -1,155 +1,0 @@
-/*!\file: transient_2d_core.cpp
- * \brief: core of the transient_2d solution 
- */ 
-
-#include "../toolkits/toolkits.h"
-#include "../objects/objects.h"
-#include "../shared/shared.h"
-#include "../EnumDefinitions/EnumDefinitions.h"
-#include "./solutions.h"
-#include "../modules/modules.h"
-
-Results* transient_core_2d(Model* model){
-
-	extern int my_rank;
-
-	/*fem models: */
-	FemModel* fem_dh=NULL;
-	FemModel* fem_dv=NULL;
-	FemModel* fem_dhu=NULL;
-	FemModel* fem_ds=NULL;
-	FemModel* fem_sl=NULL;
-	FemModel* fem_p=NULL;
-
-	/*output: */
-	Results* results=NULL;
-	Result * result =NULL;
-
-	/*solutions: */
-	int step,i;
-	double time;
-
-	Vec u_g=NULL;
-	Vec p_g=NULL;
-	Vec h_g=NULL;
-	Vec new_h_g=NULL;
-	Vec h_g_intermediary=NULL;
-	Vec s_g=NULL;
-	Vec new_s_g=NULL;
-	Vec b_g=NULL;
-	Vec new_b_g=NULL;
-	double* vx=NULL;
-	double* vy=NULL;
-
-	/*intermediary: */
-	double finaltime;
-	double dt;
-	Results* diagnostic_results=NULL;
-	Results* prognostic_results=NULL;
-
-
-	/*flags: */
-	int verbose=0;
-	int numberofnodes;
-
-	/*dof recovery: */
-	int dof01[2]={0,1};
-	int dof2[1]={2};
-	int dof012[3]={0,1,2};
-	int dof3[1]={3};
-
-	//initialize results
-	results=new Results();
-
-	/*recover fem models: */
-	fem_dh=model->GetFormulation(DiagnosticAnalysisEnum,HorizAnalysisEnum);
-	fem_dv=model->GetFormulation(DiagnosticAnalysisEnum,VertAnalysisEnum);
-	fem_ds=model->GetFormulation(DiagnosticAnalysisEnum,StokesAnalysisEnum);
-	fem_dhu=model->GetFormulation(DiagnosticAnalysisEnum,HutterAnalysisEnum);
-	fem_sl=model->GetFormulation(SlopecomputeAnalysisEnum);
-	fem_p=model->GetFormulation(PrognosticAnalysisEnum);
-
-
-	//first recover parameters common to all solutions
-	model->FindParam(&verbose,VerboseEnum);
-	model->FindParam(&finaltime,NdtEnum);
-	model->FindParam(&dt,DtEnum);
-	model->FindParam(&numberofnodes,NumberOfNodesEnum);
-
-	/*initialize: */
-	step=1;
-	time=0;
-
-	while(time<finaltime){ //make sure we run up to finaltime.
-	
-		_printf_("%s%g%s%i%s%g\n","time [yr]: ",time,"    iteration number: ",step,"/",floor(finaltime/dt));
-
-		step+=1;
-		time+=dt;
-
-		//Deal with velocities.
-		_printf_("%s\n","computing new velocity");
-		diagnostic_results= diagnostic_core(model);
-
-		//extract u_g and p_g from diagnostic_results, and erase diagnostic_results;
-		VecFree(&u_g); diagnostic_results->FindResult(&u_g,"u_g");
-		VecFree(&p_g); diagnostic_results->FindResult(&p_g,"p_g");
-		delete diagnostic_results;
-
-		//compute new thickness
-		_printf_("%s\n","computing new thickness");
-		
-		SplitSolutionVectorx(u_g,numberofnodes,2,&vx,&vy);
-		model->UpdateInputsFromVector(vx,VxEnum,VertexEnum);
-		model->UpdateInputsFromVector(vy,VyEnum,VertexEnum);
-		xfree((void**)&vx); xfree((void**)&vy);
-		
-		prognostic_results= prognostic_core(model);
-	
-		//extract h_g prognostic_results, and erase prognostic_results;
-		prognostic_results->FindResult(&h_g_intermediary,"h_g");
-		delete prognostic_results;
-
-		//update surface and bed using the new thickness
-		_printf_("   updating geometry\n");
-		UpdateGeometryx(&new_h_g,&new_b_g,&new_s_g, fem_p->elements, fem_p->nodes,fem_p->vertices,fem_p->loads, fem_p->materials, fem_p->parameters,h_g_intermediary,b_g,s_g); 
-
-		VecFree(&h_g);h_g=new_h_g; 
-		VecFree(&s_g);s_g=new_s_g; 
-		VecFree(&b_g);b_g=new_b_g; 
-		VecFree(&h_g_intermediary);
-
-		//plug into results.
-		results->AddObject(new Result(results->Size()+1,time,step,"u_g",u_g)); 
-		results->AddObject(new Result(results->Size()+1,time,step,"p_g",p_g)); 
-		results->AddObject(new Result(results->Size()+1,time,step,"h_g",h_g)); 
-		results->AddObject(new Result(results->Size()+1,time,step,"s_g",s_g)); 
-		results->AddObject(new Result(results->Size()+1,time,step,"b_g",b_g)); 
-
-		//update inputs
-		model->UpdateInputsFromVector(h_g,ThicknessEnum,VertexEnum);
-		model->UpdateInputsFromVector(s_g,SurfaceEnum,VertexEnum);
-		model->UpdateInputsFromVector(b_g,BedEnum,VertexEnum);
-	
-		SplitSolutionVectorx(u_g,numberofnodes,2,&vx,&vy);
-		model->UpdateInputsFromVector(vx,VxEnum,VertexEnum);
-		model->UpdateInputsFromVector(vy,VyEnum,VertexEnum);
-
-	}
-
-	/*Add analysis_type to results: */
-	results->AddObject(new StringResult(results->Size()+1,AnalysisTypeEnum,0,1,EnumAsString(TransientAnalysisEnum)));
-
-	/*Free ressources:*/
-	VecFree(&u_g);
-	VecFree(&p_g);
-	VecFree(&h_g);
-	VecFree(&s_g);
-	VecFree(&b_g);
-	xfree((void**)&vx);
-	xfree((void**)&vy);
-	
-	//return: 
-	return results;
-
-}
Index: issm/trunk/src/c/solutions/transient_core_3d.cpp
===================================================================
--- issm/trunk/src/c/solutions/transient_core_3d.cpp	(revision 4054)
+++ 	(revision )
@@ -1,206 +1,0 @@
-/*!\file: transient_3d_core.cpp
- * \brief: core of the transient_3d solution 
- */ 
-
-#include "../toolkits/toolkits.h"
-#include "../objects/objects.h"
-#include "../shared/shared.h"
-#include "../EnumDefinitions/EnumDefinitions.h"
-#include "./solutions.h"
-#include "../modules/modules.h"
-
-Results* transient_core_3d(Model* model){
-
-	extern int my_rank;
-
-	/*fem models: */
-	FemModel* fem_dh=NULL;
-	FemModel* fem_dv=NULL;
-	FemModel* fem_dhu=NULL;
-	FemModel* fem_ds=NULL;
-	FemModel* fem_sl=NULL;
-	FemModel* fem_p=NULL;
-	FemModel* fem_t=NULL;
-	FemModel* fem_m=NULL;
-
-	/*output: */
-	Results* results=NULL;
-
-	/*solutions: */
-	int step;
-	double time;
-	
-	double* vx=NULL;
-	double* vy=NULL;
-	double* vz=NULL;
-	
-	Vec u_g=NULL;
-	Vec p_g=NULL;
-	Vec h_g=NULL;
-	Vec h_g_intermediary=NULL;
-	Vec s_g=NULL;
-	Vec b_g=NULL;
-	Vec t_g=NULL;
-	Vec t_g_average=NULL;
-	double melting_offset;
-	Vec m_g=NULL;
-
-	/*intermediary: */
-	double finaltime;
-	double dt;
-	Results* diagnostic_results=NULL;
-	Results* prognostic_results=NULL;
-
-
-	/*flags: */
-	int verbose=0;
-	int numberofnodes;
-
-	/*dof recovery: */
-	int dof01[2]={0,1};
-	int dof2[1]={2};
-	int dof012[3]={0,1,2};
-	int dof3[1]={3};
-
-	//initialize results
-	results=new Results();
-
-	/*recover fem models: */
-	fem_dh=model->GetFormulation(DiagnosticAnalysisEnum,HorizAnalysisEnum);
-	fem_dv=model->GetFormulation(DiagnosticAnalysisEnum,VertAnalysisEnum);
-	fem_ds=model->GetFormulation(DiagnosticAnalysisEnum,StokesAnalysisEnum);
-	fem_dhu=model->GetFormulation(DiagnosticAnalysisEnum,HutterAnalysisEnum);
-	fem_sl=model->GetFormulation(SlopecomputeAnalysisEnum);
-	fem_p=model->GetFormulation(PrognosticAnalysisEnum);
-	fem_t=model->GetFormulation(ThermalAnalysisEnum);
-	fem_m=model->GetFormulation(MeltingAnalysisEnum);
-
-
-	//first recover parameters common to all solutions
-	model->FindParam(&verbose,VerboseEnum);
-	model->FindParam(&finaltime,NdtEnum);
-	model->FindParam(&dt,DtEnum);
-	model->FindParam(&numberofnodes,NumberOfDofsPerNodeEnum);
-
-	/*initialize: */
-	step=1;
-	time=0;
-
-	while(time<finaltime){ //make sure we run up to finaltime.
-	
-		if(verbose)_printf_("%s%g%s%i%s%g\n","time [yr]: ",time,"    iteration number: ",step,"/",floor(finaltime/dt));
-
-		step+=1;
-		time+=dt;
-
-		//Deal with temperature first 
-		if(verbose)_printf_("%s\n","computing temperature");
-		thermal_core_nonlinear(&t_g,&melting_offset,fem_t,ThermalAnalysisEnum,TransientAnalysisEnum);
-		if(verbose)_printf_("%s\n","computing melting");
-		model->UpdateInputsFromVector(t_g,TemperatureEnum,VertexEnum);
-		model->UpdateInputsFromConstant(melting_offset,MeltingOffsetEnum);
-		diagnostic_core_linear(&m_g,fem_m,MeltingAnalysisEnum,TransientAnalysisEnum);
-
-		//Compute depth averaged temperature and add to inputs
-		if(verbose)_printf_("%s\n","computing depth average temperature");
-		VecDuplicatePatch(&t_g_average,t_g); 
-		FieldDepthAveragex( t_g_average, fem_t->elements,fem_t->nodes, fem_t->vertices,fem_t->loads, fem_t->materials,fem_t->parameters,"temperature");
-		model->UpdateInputsFromVector(t_g_average,TemperatureAverageEnum,VertexEnum);
-		VecFree(&t_g_average); //not needed anymore
-
-		//Deal with velocities.
-		diagnostic_results=diagnostic_core(model);
-
-		//extract u_g and p_g from diagnostic_results, and erase diagnostic_results;
-		diagnostic_results->FindResult(&u_g,"u_g");
-		diagnostic_results->FindResult(&p_g,"p_g");
-		delete diagnostic_results;
-
-		//compute new thickness
-		if(verbose)_printf_("%s\n","computing new thickness");
-	
-		SplitSolutionVectorx(u_g,numberofnodes,3,&vx,&vy,&vz);
-		model->UpdateInputsFromVector(vx,VxEnum,VertexEnum);
-		model->UpdateInputsFromVector(vy,VyEnum,VertexEnum);
-		model->UpdateInputsFromVector(vz,VzEnum,VertexEnum);
-		xfree((void**)&vx); xfree((void**)&vy); xfree((void**)&vz);
-	
-		prognostic_results=prognostic_core(model);
-	
-		//extract h_g prognostic_results, and erase prognostic_results;
-		prognostic_results->FindResult(&h_g_intermediary,"h_g");
-		delete prognostic_results;
-
-		//update surface and bed using the new thickness
-		if(verbose)_printf_("   updating geometry\n");
-		UpdateGeometryx(&h_g,&b_g,&s_g, fem_p->elements, fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials, fem_p->parameters,h_g_intermediary,b_g,s_g); 
-		VecFree(&h_g_intermediary);
-		
-		if(verbose)_printf_("%s\n","updating node positions");
-		UpdateVertexPositionsx( fem_dh ->vertices,h_g,b_g);
-		UpdateVertexPositionsx( fem_dv ->vertices,h_g,b_g);
-		UpdateVertexPositionsx( fem_dhu->vertices,h_g,b_g);
-		UpdateVertexPositionsx( fem_ds ->vertices,h_g,b_g);
-		UpdateVertexPositionsx( fem_sl ->vertices,h_g,b_g);
-		UpdateVertexPositionsx( fem_p  ->vertices,h_g,b_g);
-		UpdateVertexPositionsx( fem_t  ->vertices,h_g,b_g);
-		UpdateVertexPositionsx( fem_m  ->vertices,h_g,b_g);
-
-		//plug into results.
-		results->AddObject(new Result(results->Size()+1,time,step,"u_g",u_g)); 
-		results->AddObject(new Result(results->Size()+1,time,step,"p_g",p_g)); 
-		results->AddObject(new Result(results->Size()+1,time,step,"h_g",h_g)); 
-		results->AddObject(new Result(results->Size()+1,time,step,"s_g",s_g)); 
-		results->AddObject(new Result(results->Size()+1,time,step,"b_g",b_g)); 
-		results->AddObject(new Result(results->Size()+1,time,step,"t_g",t_g)); 
-		results->AddObject(new Result(results->Size()+1,time,step,"m_g",m_g)); 
-
-		//update inputs
-		model->UpdateInputsFromVector(h_g,ThicknessEnum,VertexEnum);
-		model->UpdateInputsFromVector(s_g,SurfaceEnum,VertexEnum);
-		model->UpdateInputsFromVector(b_g,BedEnum,VertexEnum);
-	
-		SplitSolutionVectorx(u_g,numberofnodes,3,&vx,&vy,&vz);
-		model->UpdateInputsFromVector(vx,VxEnum,VertexEnum);
-		model->UpdateInputsFromVector(vy,VyEnum,VertexEnum);
-		model->UpdateInputsFromVector(vz,VzEnum,VertexEnum);
-		xfree((void**)&vx); xfree((void**)&vy); xfree((void**)&vz);
-		
-		model->UpdateInputsFromVector(p_g,PressureEnum,VertexEnum);
-		model->UpdateInputsFromVector(t_g,TemperatureEnum,VertexEnum);
-
-		//Save temporary results every 5 steps
-		if (step%5==0){
-			_printf_("%s","      saving temporary results...");
-
-			char*    outputfilename=NULL;
-
-			model->FindParam(&outputfilename,OutputFileNameEnum);
-			OutputResults(results,outputfilename);
-
-			xfree((void**)&outputfilename);
-
-			_printf_("%s\n"," done.");
-		}
-	}
-
-	/*Add analysis_type to results: */
-	results->AddObject(new StringResult(results->Size()+1,AnalysisTypeEnum,0,1,EnumAsString(TransientAnalysisEnum)));
-
-
-	/*Free ressources:*/
-	VecFree(&u_g);
-	VecFree(&p_g);
-	VecFree(&h_g);
-	VecFree(&s_g);
-	VecFree(&b_g);
-	VecFree(&t_g);
-	VecFree(&m_g);
-	xfree((void**)&vx); 
-	xfree((void**)&vy); 
-	xfree((void**)&vz);
-	
-	//return: 
-	return results;
-
-}
Index: issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp
===================================================================
--- issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp	(revision 4054)
+++ issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp	(revision 4055)
@@ -9,5 +9,5 @@
 #include "./solutions.h"
 
-void solver_diagnostic_nonlinear(Vec* pug,Mat* pKff0,Mat* pKfs0, FemModel* fem,bool conserve_loads, int analysis_type,int sub_analysis_type){
+void solver_diagnostic_nonlinear(Vec* pug,Mat* pKff0,Mat* pKfs0, FemModel* fem,bool conserve_loads){
 
 
@@ -61,5 +61,5 @@
 
 	/*Start non-linear iteration using input velocity: */
-	GetSolutionFromInputsx(&ug, fem->elements, fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters, analysis_type, sub_analysis_type);
+	GetSolutionFromInputsx(&ug, fem->elements, fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters);
 	Reducevectorgtofx(&uf, ug, fem->nodesets);
 
@@ -72,9 +72,9 @@
 		if (verbose) _printf_("   Generating matrices\n");
 		//*Generate system matrices
-		SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); 
+		SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters,kflag,pflag);
 
 		if (verbose) _printf_("   Generating penalty matrices\n");
 		//*Generate penalty system matrices
-		PenaltySystemMatricesx(Kgg, pg,NULL,fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); 
+		PenaltySystemMatricesx(Kgg, pg,NULL,fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters,kflag,pflag);
 
 		if (verbose) _printf_("   reducing matrix from g to f set\n");
@@ -102,9 +102,9 @@
 
 		//Update inputs using new solution:
-		UpdateInputsFromSolutionx( fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,ug,analysis_type, sub_analysis_type);
+		UpdateInputsFromSolutionx( fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,ug);
 
 		//Deal with penalty loads
 		if (verbose) _printf_("   penalty constraints\n");
-		PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters,analysis_type,sub_analysis_type); 
+		PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters);
 
 		if(verbose)_printf_("   number of unstable constraints: %i\n",num_unstable_constraints);
@@ -146,5 +146,5 @@
 		kflag=1; pflag=0; //stiffness generation only
 	
-		SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); 
+		SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,loads,fem->materials,fem->parameters,kflag,pflag);
 		Reducematrixfromgtofx(&Kff,&Kfs,Kgg,fem->Gmn,fem->nodesets);
 		MatFree(&Kgg);VecFree(&pg);
Index: issm/trunk/src/c/solvers/solver_linear.cpp
===================================================================
--- issm/trunk/src/c/solvers/solver_linear.cpp	(revision 4054)
+++ issm/trunk/src/c/solvers/solver_linear.cpp	(revision 4055)
@@ -8,5 +8,5 @@
 #include "../modules/modules.h"
 
-void solver_linear(Vec* pug, FemModel* fem,int analysis_type,int sub_analysis_type){
+void solver_linear(Vec* pug, FemModel* fem){
 
 	/*parameters:*/
@@ -33,9 +33,9 @@
 	//*Generate system matrices
 	if (verbose) _printf_("   Generating matrices\n");
-	SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); 
+	SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag);
 
 	if (verbose) _printf_("   Generating penalty matrices\n");
 	//*Generate penalty system matrices
-	PenaltySystemMatricesx(Kgg, pg,NULL,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); 
+	PenaltySystemMatricesx(Kgg, pg,NULL,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag);
 
 	/*!Reduce matrix from g to f size:*/
@@ -56,5 +56,5 @@
 
 	//Update inputs using new solution:
-	UpdateInputsFromSolutionx( fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,ug,analysis_type, sub_analysis_type);
+	UpdateInputsFromSolutionx( fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,ug);
 
 	/*free ressources: */
Index: issm/trunk/src/c/solvers/solver_thermal_nonlinear.cpp
===================================================================
--- issm/trunk/src/c/solvers/solver_thermal_nonlinear.cpp	(revision 4054)
+++ issm/trunk/src/c/solvers/solver_thermal_nonlinear.cpp	(revision 4055)
@@ -8,5 +8,5 @@
 #include "../modules/modules.h"
 
-void solver_thermal_nonlinear(Vec* ptg,double* pmelting_offset,FemModel* fem,int analysis_type,int sub_analysis_type){
+void solver_thermal_nonlinear(Vec* ptg,double* pmelting_offset,FemModel* fem){
 
 	/*solution : */
@@ -63,5 +63,5 @@
 			/*Compute Kgg_nopenalty and pg_nopenalty once for all: */
 			if (count==1){
-				SystemMatricesx(&Kgg_nopenalty, &pg_nopenalty,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); 
+				SystemMatricesx(&Kgg_nopenalty, &pg_nopenalty,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag);
 			}
 
@@ -71,10 +71,10 @@
 
 			//apply penalties each time
-			PenaltySystemMatricesx(Kgg, pg,&melting_offset,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); 
+			PenaltySystemMatricesx(Kgg, pg,&melting_offset,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag);
 		}
 		else{
-			SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); 
+			SystemMatricesx(&Kgg, &pg,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag);
 			//apply penalties
-			PenaltySystemMatricesx(Kgg, pg,&melting_offset,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag,analysis_type,sub_analysis_type); 
+			PenaltySystemMatricesx(Kgg, pg,&melting_offset,fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,kflag,pflag);
 		}
 
@@ -108,9 +108,9 @@
 		//Deal with penalty loads
 		if (verbose) _printf_("   penalty constraints\n");
-		PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters,analysis_type,sub_analysis_type); 
+		PenaltyConstraintsx(&constraints_converged, &num_unstable_constraints, fem->elements,fem->nodes,fem->vertices,fem->loads,fem->materials,fem->parameters);
 		
 		//Update inputs using new solution:
 		UpdateInputsFromVectorx( fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,tg,TemperatureEnum,VertexEnum);
-		UpdateInputsFromSolutionx(fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,tg,analysis_type, sub_analysis_type);
+		UpdateInputsFromSolutionx(fem->elements,fem->nodes, fem->vertices, fem->loads, fem->materials, fem->parameters,tg);
 
 		if (!converged){
Index: issm/trunk/src/c/solvers/solvers.h
===================================================================
--- issm/trunk/src/c/solvers/solvers.h	(revision 4054)
+++ issm/trunk/src/c/solvers/solvers.h	(revision 4055)
@@ -12,8 +12,8 @@
 class FemModel;
 
-void solver_thermal_nonlinear(Vec* ptg,double* pmelting_offset,FemModel* fem,int analysis_type,int sub_analysis_type);
-void solver_diagnostic_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, FemModel* fem,bool conserve_loads,int analysis_type,int sub_analysis_type);
-void solver_linear(Vec* pug, FemModel* femmodel,int  analysis_type,int sub_analysis_type);
-void solver_adjoint(Vec* pug, FemModel* femmodel,int  analysis_type,int sub_analysis_type);
+void solver_thermal_nonlinear(Vec* ptg,double* pmelting_offset,FemModel* femmodel);
+void solver_diagnostic_nonlinear(Vec* pug,Mat* pK_ff0,Mat* pK_fs0, FemModel* femmodel,bool conserve_loads);
+void solver_linear(Vec* pug, FemModel* femmodel);
+void solver_adjoint(Vec* pug, FemModel* femmodel);
 
 
