Index: ../trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.h
===================================================================
--- ../trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.h	(revision 18928)
+++ ../trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.h	(revision 18929)
@@ -12,26 +12,26 @@
 
 	public:
 		/*Model processing*/
+		void CreateConstraints(Constraints* constraints,IoModel* iomodel);
+		void CreateLoads(Loads* loads, IoModel* iomodel);
+		void CreateNodes(Nodes* nodes,IoModel* iomodel);
 		int  DofsPerNode(int** doflist,int domaintype,int approximation);
+		void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
 		void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum);
-		void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
-		void CreateNodes(Nodes* nodes,IoModel* iomodel);
-		void CreateConstraints(Constraints* constraints,IoModel* iomodel);
-		void CreateLoads(Loads* loads, IoModel* iomodel);
 
 		/*Finite element Analysis*/
 		void           Core(FemModel* femmodel);
 		ElementVector* CreateDVector(Element* element);
 		ElementMatrix* CreateJacobianMatrix(Element* element);
 		ElementMatrix* CreateKMatrix(Element* element);
+		ElementMatrix* CreateKMatrixBed(Element* element);
+		ElementMatrix* CreateKMatrixSurface(Element* element);
 		ElementMatrix* CreateKMatrixVolume(Element* element);
-		ElementMatrix* CreateKMatrixSurface(Element* element);
-		ElementMatrix* CreateKMatrixBed(Element* element);
 		ElementVector* CreatePVector(Element* element);
-		void GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
-		void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
-		void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
-		void InputUpdateFromSolution(IssmDouble* solution,Element* element);
-		void UpdateConstraints(FemModel* femmodel);
+		void           GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
+		void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
+		void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
+		void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
+		void           UpdateConstraints(FemModel* femmodel);
 };
 #endif
Index: ../trunk-jpl/src/c/analyses/GiaAnalysis.h
===================================================================
--- ../trunk-jpl/src/c/analyses/GiaAnalysis.h	(revision 18928)
+++ ../trunk-jpl/src/c/analyses/GiaAnalysis.h	(revision 18929)
@@ -12,12 +12,12 @@
 
 	public:
 		/*Model processing*/
+		void CreateConstraints(Constraints* constraints,IoModel* iomodel);
+		void CreateLoads(Loads* loads, IoModel* iomodel);
+		void CreateNodes(Nodes* nodes,IoModel* iomodel);
 		int  DofsPerNode(int** doflist,int domaintype,int approximation);
+		void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
 		void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum);
-		void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
-		void CreateNodes(Nodes* nodes,IoModel* iomodel);
-		void CreateConstraints(Constraints* constraints,IoModel* iomodel);
-		void CreateLoads(Loads* loads, IoModel* iomodel);
 
 		/*Finite element Analysis*/
 		void           Core(FemModel* femmodel);
@@ -25,9 +25,9 @@
 		ElementMatrix* CreateJacobianMatrix(Element* element);
 		ElementMatrix* CreateKMatrix(Element* element);
 		ElementVector* CreatePVector(Element* element);
-		void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
-		void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
-		void InputUpdateFromSolution(IssmDouble* solution,Element* element);
-		void UpdateConstraints(FemModel* femmodel);
+		void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
+		void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
+		void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
+		void           UpdateConstraints(FemModel* femmodel);
 };
 #endif
Index: ../trunk-jpl/src/c/analyses/ExtrudeFromTopAnalysis.h
===================================================================
--- ../trunk-jpl/src/c/analyses/ExtrudeFromTopAnalysis.h	(revision 18928)
+++ ../trunk-jpl/src/c/analyses/ExtrudeFromTopAnalysis.h	(revision 18929)
@@ -12,26 +12,26 @@
 
 	public:
 		/*Model processing*/
+		void CreateConstraints(Constraints* constraints,IoModel* iomodel);
+		void CreateLoads(Loads* loads, IoModel* iomodel);
+		void CreateNodes(Nodes* nodes,IoModel* iomodel);
 		int  DofsPerNode(int** doflist,int domaintype,int approximation);
+		void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
 		void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum);
-		void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
-		void CreateNodes(Nodes* nodes,IoModel* iomodel);
-		void CreateConstraints(Constraints* constraints,IoModel* iomodel);
-		void CreateLoads(Loads* loads, IoModel* iomodel);
 
 		/*Finite element Analysis*/
 		void           Core(FemModel* femmodel);
 		ElementVector* CreateDVector(Element* element);
 		ElementMatrix* CreateJacobianMatrix(Element* element);
 		ElementMatrix* CreateKMatrix(Element* element);
+		ElementMatrix* CreateKMatrixBed(Element* element);
+		ElementMatrix* CreateKMatrixSurface(Element* element);
 		ElementMatrix* CreateKMatrixVolume(Element* element);
-		ElementMatrix* CreateKMatrixSurface(Element* element);
-		ElementMatrix* CreateKMatrixBed(Element* element);
 		ElementVector* CreatePVector(Element* element);
-		void GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
-		void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
-		void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
-		void InputUpdateFromSolution(IssmDouble* solution,Element* element);
-		void UpdateConstraints(FemModel* femmodel);
+		void           GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
+		void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
+		void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
+		void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
+		void           UpdateConstraints(FemModel* femmodel);
 };
 #endif
Index: ../trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.cpp	(revision 18928)
+++ ../trunk-jpl/src/c/analyses/ExtrudeFromBaseAnalysis.cpp	(revision 18929)
@@ -5,11 +5,18 @@
 #include "../modules/modules.h"
 
 /*Model processing*/
+void ExtrudeFromBaseAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
+}/*}}}*/
+void ExtrudeFromBaseAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
+}/*}}}*/
+void ExtrudeFromBaseAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
+
+	::CreateNodes(nodes,iomodel,ExtrudeFromBaseAnalysisEnum,P1Enum);
+
+}/*}}}*/
 int  ExtrudeFromBaseAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
 	return 1;
 }/*}}}*/
-void ExtrudeFromBaseAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
-}/*}}}*/
 void ExtrudeFromBaseAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
 
 	int counter=0;
@@ -25,15 +32,8 @@
 		iomodel->FetchDataToInput(elements,MeshVertexonbaseEnum);
 	}
 }/*}}}*/
-void ExtrudeFromBaseAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
-
-	::CreateNodes(nodes,iomodel,ExtrudeFromBaseAnalysisEnum,P1Enum);
-
+void ExtrudeFromBaseAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
 }/*}}}*/
-void ExtrudeFromBaseAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
-}/*}}}*/
-void ExtrudeFromBaseAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
-}/*}}}*/
 
 /*Finite Element Analysis*/
 void           ExtrudeFromBaseAnalysis::Core(FemModel* femmodel){/*{{{*/
@@ -60,47 +60,47 @@
 	delete Ke3;
 	return Ke;
 }/*}}}*/
-ElementMatrix* ExtrudeFromBaseAnalysis::CreateKMatrixVolume(Element* element){/*{{{*/
+ElementMatrix* ExtrudeFromBaseAnalysis::CreateKMatrixBed(Element* element){/*{{{*/
 
+	if(!element->IsOnBase()) return NULL;
+
 	/*Intermediaries */
-	IssmDouble  Jdet,D;
-	IssmDouble *xyz_list = NULL;
+	int         dim;
+	IssmDouble  Jdet,D,normal[3];
+	IssmDouble *xyz_list_base = NULL;
 
 	/*Get dimension*/
-	int dim;
 	element->FindParam(&dim,DomainDimensionEnum);
 
 	/*Fetch number of nodes and dof for this finite element*/
 	int numnodes = element->GetNumberOfNodes();
 
 	/*Initialize Element vector and other vectors*/
-	ElementMatrix* Ke     = element->NewElementMatrix();
-	IssmDouble*    B      = xNew<IssmDouble>(numnodes);
-	IssmDouble*    Bprime = xNew<IssmDouble>(numnodes);
+	ElementMatrix* Ke    = element->NewElementMatrix();
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
 
 	/*Retrieve all inputs and parameters*/
-	element->GetVerticesCoordinates(&xyz_list);
+	element->GetVerticesCoordinatesBase(&xyz_list_base);
+	element->NormalBase(&normal[0],xyz_list_base);
 
 	/* Start  looping on the number of gaussian points: */
-	Gauss* gauss=element->NewGauss(2);
+	Gauss* gauss=element->NewGaussBase(2);
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
 		gauss->GaussPoint(ig);
 
-		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-		element->NodalFunctions(Bprime,gauss);
-		GetB(B,element,dim,xyz_list,gauss);
-		D=gauss->weight*Jdet;
+		element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
+		element->NodalFunctions(basis,gauss);
+		D = - gauss->weight*Jdet*normal[dim-1]; 
 
-		TripleMultiply(B,1,numnodes,1,
+		TripleMultiply(basis,1,numnodes,1,
 					&D,1,1,0,
-					Bprime,1,numnodes,0,
+					basis,1,numnodes,0,
 					&Ke->values[0],1);
-	} 
+	}
 
 	/*Clean up and return*/
-	xDelete<IssmDouble>(xyz_list);
-	xDelete<IssmDouble>(B);
-	xDelete<IssmDouble>(Bprime);
+	xDelete<IssmDouble>(xyz_list_base);
+	xDelete<IssmDouble>(basis);
 	delete gauss;
 	return Ke;
 }
@@ -150,47 +150,47 @@
 	return Ke;
 }
 /*}}}*/
-ElementMatrix* ExtrudeFromBaseAnalysis::CreateKMatrixBed(Element* element){/*{{{*/
+ElementMatrix* ExtrudeFromBaseAnalysis::CreateKMatrixVolume(Element* element){/*{{{*/
 
-	if(!element->IsOnBase()) return NULL;
-
 	/*Intermediaries */
-	int         dim;
-	IssmDouble  Jdet,D,normal[3];
-	IssmDouble *xyz_list_base = NULL;
+	IssmDouble  Jdet,D;
+	IssmDouble *xyz_list = NULL;
 
 	/*Get dimension*/
+	int dim;
 	element->FindParam(&dim,DomainDimensionEnum);
 
 	/*Fetch number of nodes and dof for this finite element*/
 	int numnodes = element->GetNumberOfNodes();
 
 	/*Initialize Element vector and other vectors*/
-	ElementMatrix* Ke    = element->NewElementMatrix();
-	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+	ElementMatrix* Ke     = element->NewElementMatrix();
+	IssmDouble*    B      = xNew<IssmDouble>(numnodes);
+	IssmDouble*    Bprime = xNew<IssmDouble>(numnodes);
 
 	/*Retrieve all inputs and parameters*/
-	element->GetVerticesCoordinatesBase(&xyz_list_base);
-	element->NormalBase(&normal[0],xyz_list_base);
+	element->GetVerticesCoordinates(&xyz_list);
 
 	/* Start  looping on the number of gaussian points: */
-	Gauss* gauss=element->NewGaussBase(2);
+	Gauss* gauss=element->NewGauss(2);
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
 		gauss->GaussPoint(ig);
 
-		element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
-		element->NodalFunctions(basis,gauss);
-		D = - gauss->weight*Jdet*normal[dim-1]; 
+		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		element->NodalFunctions(Bprime,gauss);
+		GetB(B,element,dim,xyz_list,gauss);
+		D=gauss->weight*Jdet;
 
-		TripleMultiply(basis,1,numnodes,1,
+		TripleMultiply(B,1,numnodes,1,
 					&D,1,1,0,
-					basis,1,numnodes,0,
+					Bprime,1,numnodes,0,
 					&Ke->values[0],1);
-	}
+	} 
 
 	/*Clean up and return*/
-	xDelete<IssmDouble>(xyz_list_base);
-	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(B);
+	xDelete<IssmDouble>(Bprime);
 	delete gauss;
 	return Ke;
 }
@@ -198,7 +198,7 @@
 ElementVector* ExtrudeFromBaseAnalysis::CreatePVector(Element* element){/*{{{*/
 	return NULL;
 }/*}}}*/
-void ExtrudeFromBaseAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void           ExtrudeFromBaseAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 	/*	Compute B  matrix. B=[dh1/dz dh2/dz dh3/dz dh4/dz dh5/dz dh6/dz];
 		where hi is the interpolation function for node i.*/
 
@@ -218,19 +218,19 @@
 	xDelete<IssmDouble>(dbasis);
 }
 /*}}}*/
-void ExtrudeFromBaseAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
+void           ExtrudeFromBaseAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
 	   _error_("not implemented yet");
 }/*}}}*/
-void ExtrudeFromBaseAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
+void           ExtrudeFromBaseAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
 	_error_("Not implemented yet");
 }/*}}}*/
-void ExtrudeFromBaseAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
+void           ExtrudeFromBaseAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
 
 	int inputenum;
 	element->FindParam(&inputenum,InputToExtrudeEnum);
 	element->InputUpdateFromSolutionOneDof(solution,inputenum);
 }/*}}}*/
-void ExtrudeFromBaseAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
+void           ExtrudeFromBaseAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
 	/*Default, do nothing*/
 	return;
 }/*}}}*/
Index: ../trunk-jpl/src/c/analyses/DepthAverageAnalysis.h
===================================================================
--- ../trunk-jpl/src/c/analyses/DepthAverageAnalysis.h	(revision 18928)
+++ ../trunk-jpl/src/c/analyses/DepthAverageAnalysis.h	(revision 18929)
@@ -12,12 +12,12 @@
 
 	public:
 		/*Model processing*/
+		void CreateConstraints(Constraints* constraints,IoModel* iomodel);
+		void CreateLoads(Loads* loads, IoModel* iomodel);
+		void CreateNodes(Nodes* nodes,IoModel* iomodel);
 		int  DofsPerNode(int** doflist,int domaintype,int approximation);
+		void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
 		void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum);
-		void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
-		void CreateNodes(Nodes* nodes,IoModel* iomodel);
-		void CreateConstraints(Constraints* constraints,IoModel* iomodel);
-		void CreateLoads(Loads* loads, IoModel* iomodel);
 
 		/*Finite element Analysis*/
 		void           Core(FemModel* femmodel);
@@ -25,10 +25,10 @@
 		ElementMatrix* CreateJacobianMatrix(Element* element);
 		ElementMatrix* CreateKMatrix(Element* element);
 		ElementVector* CreatePVector(Element* element);
-		void GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
-		void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
-		void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
-		void InputUpdateFromSolution(IssmDouble* solution,Element* element);
-		void UpdateConstraints(FemModel* femmodel);
+		void           GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
+		void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
+		void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
+		void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
+		void           UpdateConstraints(FemModel* femmodel);
 };
 #endif
Index: ../trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.h
===================================================================
--- ../trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.h	(revision 18928)
+++ ../trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.h	(revision 18929)
@@ -12,12 +12,12 @@
 
 	public:
 		/*Model processing*/
+		void CreateConstraints(Constraints* constraints,IoModel* iomodel);
+		void CreateLoads(Loads* loads, IoModel* iomodel);
+		void CreateNodes(Nodes* nodes,IoModel* iomodel);
 		int  DofsPerNode(int** doflist,int domaintype,int approximation);
+		void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
 		void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum);
-		void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
-		void CreateNodes(Nodes* nodes,IoModel* iomodel);
-		void CreateConstraints(Constraints* constraints,IoModel* iomodel);
-		void CreateLoads(Loads* loads, IoModel* iomodel);
 
 		/*Finite element Analysis*/
 		void           Core(FemModel* femmodel);
@@ -25,11 +25,11 @@
 		ElementMatrix* CreateJacobianMatrix(Element* element);
 		ElementMatrix* CreateKMatrix(Element* element);
 		ElementVector* CreatePVector(Element* element);
-		void GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
-		void GetBprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
-		void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
-		void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
-		void InputUpdateFromSolution(IssmDouble* solution,Element* element);
-		void UpdateConstraints(FemModel* femmodel);
+		void           GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
+		void           GetBprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
+		void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
+		void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
+		void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
+		void           UpdateConstraints(FemModel* femmodel);
 };
 #endif
Index: ../trunk-jpl/src/c/analyses/ExtrudeFromTopAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/ExtrudeFromTopAnalysis.cpp	(revision 18928)
+++ ../trunk-jpl/src/c/analyses/ExtrudeFromTopAnalysis.cpp	(revision 18929)
@@ -5,11 +5,18 @@
 #include "../modules/modules.h"
 
 /*Model processing*/
+void ExtrudeFromTopAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
+}/*}}}*/
+void ExtrudeFromTopAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
+}/*}}}*/
+void ExtrudeFromTopAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
+
+	::CreateNodes(nodes,iomodel,ExtrudeFromTopAnalysisEnum,P1Enum);
+
+}/*}}}*/
 int  ExtrudeFromTopAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
 	return 1;
 }/*}}}*/
-void ExtrudeFromTopAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
-}/*}}}*/
 void ExtrudeFromTopAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
 
 	int counter=0;
@@ -25,15 +32,8 @@
 		iomodel->FetchDataToInput(elements,MeshVertexonbaseEnum);
 	}
 }/*}}}*/
-void ExtrudeFromTopAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
-
-	::CreateNodes(nodes,iomodel,ExtrudeFromTopAnalysisEnum,P1Enum);
-
+void ExtrudeFromTopAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
 }/*}}}*/
-void ExtrudeFromTopAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
-}/*}}}*/
-void ExtrudeFromTopAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
-}/*}}}*/
 
 /*Finite Element Analysis*/
 void           ExtrudeFromTopAnalysis::Core(FemModel* femmodel){/*{{{*/
@@ -198,7 +198,7 @@
 ElementVector* ExtrudeFromTopAnalysis::CreatePVector(Element* element){/*{{{*/
 	return NULL;
 }/*}}}*/
-void ExtrudeFromTopAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void           ExtrudeFromTopAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 	/*	Compute B  matrix. B=[dh1/dz dh2/dz dh3/dz dh4/dz dh5/dz dh6/dz];
 		where hi is the interpolation function for node i.*/
 
@@ -218,19 +218,19 @@
 	xDelete<IssmDouble>(dbasis);
 }
 /*}}}*/
-void ExtrudeFromTopAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
+void           ExtrudeFromTopAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
 	   _error_("not implemented yet");
 }/*}}}*/
-void ExtrudeFromTopAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
+void           ExtrudeFromTopAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
 	_error_("Not implemented yet");
 }/*}}}*/
-void ExtrudeFromTopAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
+void           ExtrudeFromTopAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
 
 	int inputenum;
 	element->FindParam(&inputenum,InputToExtrudeEnum);
 	element->InputUpdateFromSolutionOneDof(solution,inputenum);
 }/*}}}*/
-void ExtrudeFromTopAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
+void           ExtrudeFromTopAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
 	/*Default, do nothing*/
 	return;
 }/*}}}*/
Index: ../trunk-jpl/src/c/analyses/DepthAverageAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/DepthAverageAnalysis.cpp	(revision 18928)
+++ ../trunk-jpl/src/c/analyses/DepthAverageAnalysis.cpp	(revision 18929)
@@ -5,11 +5,18 @@
 #include "../modules/modules.h"
 
 /*Model processing*/
+void DepthAverageAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
+}/*}}}*/
+void DepthAverageAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
+}/*}}}*/
+void DepthAverageAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
+
+	::CreateNodes(nodes,iomodel,DepthAverageAnalysisEnum,P1Enum);
+
+}/*}}}*/
 int  DepthAverageAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
 	return 1;
 }/*}}}*/
-void DepthAverageAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
-}/*}}}*/
 void DepthAverageAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
 
 	int counter=0;
@@ -25,15 +32,8 @@
 		iomodel->FetchDataToInput(elements,MeshVertexonbaseEnum);
 	}
 }/*}}}*/
-void DepthAverageAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
-
-	::CreateNodes(nodes,iomodel,DepthAverageAnalysisEnum,P1Enum);
-
+void DepthAverageAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
 }/*}}}*/
-void DepthAverageAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
-}/*}}}*/
-void DepthAverageAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
-}/*}}}*/
 
 /*Finite Element Analysis*/
 void           DepthAverageAnalysis::Core(FemModel* femmodel){/*{{{*/
@@ -136,7 +136,7 @@
 	delete gauss;
 	return pe;
 }/*}}}*/
-void DepthAverageAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void           DepthAverageAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 	/*	Compute B  matrix. B=[dh1/dz dh2/dz dh3/dz dh4/dz dh5/dz dh6/dz];
 		where hi is the interpolation function for node i.*/
 
@@ -156,19 +156,19 @@
 	xDelete<IssmDouble>(dbasis);
 }
 /*}}}*/
-void DepthAverageAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
+void           DepthAverageAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
 	   _error_("not implemented yet");
 }/*}}}*/
-void DepthAverageAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
+void           DepthAverageAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
 	_error_("Not implemented yet");
 }/*}}}*/
-void DepthAverageAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
+void           DepthAverageAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
 
 	int inputenum;
 	element->FindParam(&inputenum,InputToDepthaverageEnum);
 	element->InputUpdateFromSolutionOneDof(solution,inputenum);
 }/*}}}*/
-void DepthAverageAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
+void           DepthAverageAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
 	/*Default, do nothing*/
 	return;
 }/*}}}*/
Index: ../trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp	(revision 18928)
+++ ../trunk-jpl/src/c/analyses/FreeSurfaceTopAnalysis.cpp	(revision 18929)
@@ -5,11 +5,58 @@
 #include "../modules/modules.h"
 
 /*Model processing*/
+void FreeSurfaceTopAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
+}/*}}}*/
+void FreeSurfaceTopAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
+
+	/*Intermediaries*/
+	int penpair_ids[2];
+	int count=0;
+	int numvertex_pairing;
+
+	/*Create Penpair for vertex_pairing: */
+	IssmDouble *vertex_pairing=NULL;
+	IssmDouble *nodeonsurface=NULL;
+	iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,MasstransportVertexPairingEnum);
+	if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(&nodeonsurface,NULL,NULL,MeshVertexonsurfaceEnum);
+	for(int i=0;i<numvertex_pairing;i++){
+
+		if(iomodel->my_vertices[reCast<int>(vertex_pairing[2*i+0])-1]){
+
+			/*In debugging mode, check that the second node is in the same cpu*/
+			_assert_(iomodel->my_vertices[reCast<int>(vertex_pairing[2*i+1])-1]);
+
+			/*Skip if one of the two is not on the bed*/
+			if(iomodel->domaintype!=Domain2DhorizontalEnum){
+				if(!(reCast<bool>(nodeonsurface[reCast<int>(vertex_pairing[2*i+0])-1])) || !(reCast<bool>(nodeonsurface[reCast<int>(vertex_pairing[2*i+1])-1]))) continue;
+			}
+
+			/*Get node ids*/
+			penpair_ids[0]=iomodel->nodecounter+reCast<int>(vertex_pairing[2*i+0]);
+			penpair_ids[1]=iomodel->nodecounter+reCast<int>(vertex_pairing[2*i+1]);
+
+			/*Create Load*/
+			loads->AddObject(new Penpair(
+							iomodel->loadcounter+count+1,
+							&penpair_ids[0],
+							FreeSurfaceTopAnalysisEnum));
+			count++;
+		}
+	}
+
+	/*free ressources: */
+	iomodel->DeleteData(vertex_pairing,MasstransportVertexPairingEnum);
+	iomodel->DeleteData(nodeonsurface,MeshVertexonsurfaceEnum);
+}/*}}}*/
+void FreeSurfaceTopAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
+
+	if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
+	::CreateNodes(nodes,iomodel,FreeSurfaceTopAnalysisEnum,P1Enum);
+	iomodel->DeleteData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
+}/*}}}*/
 int  FreeSurfaceTopAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
 	return 1;
 }/*}}}*/
-void FreeSurfaceTopAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
-}/*}}}*/
 void FreeSurfaceTopAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
 
 	/*Now, is the model 3d? otherwise, do nothing: */
@@ -50,56 +97,9 @@
 			;
 	}
 }/*}}}*/
-void FreeSurfaceTopAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
-
-	if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
-	::CreateNodes(nodes,iomodel,FreeSurfaceTopAnalysisEnum,P1Enum);
-	iomodel->DeleteData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
+void FreeSurfaceTopAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
 }/*}}}*/
-void FreeSurfaceTopAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
-}/*}}}*/
-void FreeSurfaceTopAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
 
-	/*Intermediaries*/
-	int penpair_ids[2];
-	int count=0;
-	int numvertex_pairing;
-
-	/*Create Penpair for vertex_pairing: */
-	IssmDouble *vertex_pairing=NULL;
-	IssmDouble *nodeonsurface=NULL;
-	iomodel->FetchData(&vertex_pairing,&numvertex_pairing,NULL,MasstransportVertexPairingEnum);
-	if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(&nodeonsurface,NULL,NULL,MeshVertexonsurfaceEnum);
-	for(int i=0;i<numvertex_pairing;i++){
-
-		if(iomodel->my_vertices[reCast<int>(vertex_pairing[2*i+0])-1]){
-
-			/*In debugging mode, check that the second node is in the same cpu*/
-			_assert_(iomodel->my_vertices[reCast<int>(vertex_pairing[2*i+1])-1]);
-
-			/*Skip if one of the two is not on the bed*/
-			if(iomodel->domaintype!=Domain2DhorizontalEnum){
-				if(!(reCast<bool>(nodeonsurface[reCast<int>(vertex_pairing[2*i+0])-1])) || !(reCast<bool>(nodeonsurface[reCast<int>(vertex_pairing[2*i+1])-1]))) continue;
-			}
-
-			/*Get node ids*/
-			penpair_ids[0]=iomodel->nodecounter+reCast<int>(vertex_pairing[2*i+0]);
-			penpair_ids[1]=iomodel->nodecounter+reCast<int>(vertex_pairing[2*i+1]);
-
-			/*Create Load*/
-			loads->AddObject(new Penpair(
-							iomodel->loadcounter+count+1,
-							&penpair_ids[0],
-							FreeSurfaceTopAnalysisEnum));
-			count++;
-		}
-	}
-
-	/*free ressources: */
-	iomodel->DeleteData(vertex_pairing,MasstransportVertexPairingEnum);
-	iomodel->DeleteData(nodeonsurface,MeshVertexonsurfaceEnum);
-}/*}}}*/
-
 /*Finite Element Analysis*/
 void           FreeSurfaceTopAnalysis::Core(FemModel* femmodel){/*{{{*/
 	_error_("not implemented");
@@ -306,7 +306,7 @@
 	return pe;
 
 }/*}}}*/
-void FreeSurfaceTopAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void           FreeSurfaceTopAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
 	 * For node i, Bi can be expressed in the actual coordinate system
 	 * by: 
@@ -334,7 +334,7 @@
 	/*Clean-up*/
 	xDelete<IssmDouble>(basis);
 }/*}}}*/
-void FreeSurfaceTopAnalysis::GetBprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void           FreeSurfaceTopAnalysis::GetBprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
 	 * For node i, Bi' can be expressed in the actual coordinate system
 	 * by: 
@@ -363,17 +363,17 @@
 	xDelete<IssmDouble>(dbasis);
 
 }/*}}}*/
-void FreeSurfaceTopAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
+void           FreeSurfaceTopAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
 	   _error_("not implemented yet");
 }/*}}}*/
-void FreeSurfaceTopAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
+void           FreeSurfaceTopAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
 	_error_("Not implemented yet");
 }/*}}}*/
-void FreeSurfaceTopAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
+void           FreeSurfaceTopAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
 
 	element->InputUpdateFromSolutionOneDof(solution,SurfaceEnum);
 }/*}}}*/
-void FreeSurfaceTopAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
+void           FreeSurfaceTopAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
 	/*Default, do nothing*/
 	return;
 }/*}}}*/
Index: ../trunk-jpl/src/c/analyses/Analysis.h
===================================================================
--- ../trunk-jpl/src/c/analyses/Analysis.h	(revision 18928)
+++ ../trunk-jpl/src/c/analyses/Analysis.h	(revision 18929)
@@ -25,12 +25,12 @@
 		virtual      ~Analysis(){};
 
 		/*Model processing*/
+		virtual void CreateConstraints(Constraints* constraints,IoModel* iomodel)=0;
+		virtual void CreateLoads(Loads* loads, IoModel* iomodel)=0;
+		virtual void CreateNodes(Nodes* nodes,IoModel* iomodel)=0;
 		virtual int  DofsPerNode(int** doflist,int domaintype,int approximation)=0;
+		virtual void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type)=0;
 		virtual void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum)=0;
-		virtual void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type)=0;
-		virtual void CreateNodes(Nodes* nodes,IoModel* iomodel)=0;
-		virtual void CreateConstraints(Constraints* constraints,IoModel* iomodel)=0;
-		virtual void CreateLoads(Loads* loads, IoModel* iomodel)=0;
 
 		/*Finite element Analysis*/
 		virtual void           Core(FemModel* femmodel)=0;
@@ -38,9 +38,9 @@
 		virtual ElementMatrix* CreateJacobianMatrix(Element* element)=0;
 		virtual ElementMatrix* CreateKMatrix(Element* element)=0;
 		virtual ElementVector* CreatePVector(Element* element)=0;
-		virtual void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element)=0;
-		virtual void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index)=0;
-		virtual void InputUpdateFromSolution(IssmDouble* solution,Element* element)=0;
-		virtual void UpdateConstraints(FemModel* femmodel)=0;
+		virtual void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element)=0;
+		virtual void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index)=0;
+		virtual void           InputUpdateFromSolution(IssmDouble* solution,Element* element)=0;
+		virtual void           UpdateConstraints(FemModel* femmodel)=0;
 };
 #endif
Index: ../trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.h
===================================================================
--- ../trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.h	(revision 18928)
+++ ../trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.h	(revision 18929)
@@ -12,12 +12,12 @@
 
 	public:
 		/*Model processing*/
+		void CreateConstraints(Constraints* constraints,IoModel* iomodel);
+		void CreateLoads(Loads* loads, IoModel* iomodel);
+		void CreateNodes(Nodes* nodes,IoModel* iomodel);
 		int  DofsPerNode(int** doflist,int domaintype,int approximation);
+		void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
 		void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum);
-		void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
-		void CreateNodes(Nodes* nodes,IoModel* iomodel);
-		void CreateConstraints(Constraints* constraints,IoModel* iomodel);
-		void CreateLoads(Loads* loads, IoModel* iomodel);
 
 		/*Finite element Analysis*/
 		void           Core(FemModel* femmodel);
@@ -25,11 +25,11 @@
 		ElementMatrix* CreateJacobianMatrix(Element* element);
 		ElementMatrix* CreateKMatrix(Element* element);
 		ElementVector* CreatePVector(Element* element);
-		void GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
-		void GetBprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
-		void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
-		void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
-		void InputUpdateFromSolution(IssmDouble* solution,Element* element);
-		void UpdateConstraints(FemModel* femmodel);
+		void           GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
+		void           GetBprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
+		void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
+		void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
+		void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
+		void           UpdateConstraints(FemModel* femmodel);
 };
 #endif
Index: ../trunk-jpl/src/c/analyses/ExtrapolationAnalysis.h
===================================================================
--- ../trunk-jpl/src/c/analyses/ExtrapolationAnalysis.h	(revision 18928)
+++ ../trunk-jpl/src/c/analyses/ExtrapolationAnalysis.h	(revision 18929)
@@ -12,12 +12,12 @@
 	
  public:
 	/*Model processing*/
+	void CreateConstraints(Constraints* constraints,IoModel* iomodel);
+	void CreateLoads(Loads* loads, IoModel* iomodel);
+	void CreateNodes(Nodes* nodes,IoModel* iomodel);
 	int  DofsPerNode(int** doflist,int domaintype,int approximation);
+	void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
 	void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum);
-	void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
-	void CreateNodes(Nodes* nodes,IoModel* iomodel);
-	void CreateConstraints(Constraints* constraints,IoModel* iomodel);
-	void CreateLoads(Loads* loads, IoModel* iomodel);
 
 	/*Finite element Analysis*/
 	void           Core(FemModel* femmodel);
@@ -25,12 +25,12 @@
 	ElementMatrix* CreateJacobianMatrix(Element* element);
 	ElementMatrix* CreateKMatrix(Element* element);
 	ElementVector* CreatePVector(Element* element);
-		void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
-		void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
-	void InputUpdateFromSolution(IssmDouble* solution,Element* element);
-	void UpdateConstraints(FemModel* femmodel);
-	void GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
-	void GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss);
-	void SetConstraintsOnIce(Element* element);
+	void           GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
+	void           GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss);
+	void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
+	void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
+	void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
+	void           SetConstraintsOnIce(Element* element);
+	void           UpdateConstraints(FemModel* femmodel);
 };
 #endif
Index: ../trunk-jpl/src/c/analyses/GiaAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/GiaAnalysis.cpp	(revision 18928)
+++ ../trunk-jpl/src/c/analyses/GiaAnalysis.cpp	(revision 18929)
@@ -5,11 +5,18 @@
 #include "../modules/modules.h"
 
 /*Model processing*/
+void GiaAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
+	/*No constraints*/
+}/*}}}*/
+void GiaAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
+	/*No loads*/
+}/*}}}*/
+void GiaAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
+	::CreateNodes(nodes,iomodel,GiaAnalysisEnum,P1Enum);
+}/*}}}*/
 int  GiaAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
 	return 1;
 }/*}}}*/
-void GiaAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
-}/*}}}*/
 void GiaAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
 
 	/*Update elements: */
@@ -26,15 +33,8 @@
 	iomodel->FetchDataToInput(elements,GiaMantleViscosityEnum);
 	iomodel->FetchDataToInput(elements,GiaLithosphereThicknessEnum);
 }/*}}}*/
-void GiaAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
-	::CreateNodes(nodes,iomodel,GiaAnalysisEnum,P1Enum);
+void GiaAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
 }/*}}}*/
-void GiaAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
-	/*No constraints*/
-}/*}}}*/
-void GiaAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
-	/*No loads*/
-}/*}}}*/
 
 /*Finite Element Analysis*/
 void           GiaAnalysis::Core(FemModel* femmodel){/*{{{*/
@@ -53,16 +53,16 @@
 ElementVector* GiaAnalysis::CreatePVector(Element* element){/*{{{*/
 _error_("not implemented yet");
 }/*}}}*/
-void GiaAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
+void           GiaAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
 	   _error_("not implemented yet");
 }/*}}}*/
-void GiaAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
+void           GiaAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
 	_error_("Not implemented yet");
 }/*}}}*/
-void GiaAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
+void           GiaAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
 	_error_("not implemented yet");
 }/*}}}*/
-void GiaAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
+void           GiaAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
 	/*Default, do nothing*/
 	return;
 }/*}}}*/
Index: ../trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp	(revision 18928)
+++ ../trunk-jpl/src/c/analyses/FreeSurfaceBaseAnalysis.cpp	(revision 18929)
@@ -5,49 +5,6 @@
 #include "../modules/modules.h"
 
 /*Model processing*/
-int  FreeSurfaceBaseAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
-	return 1;
-}/*}}}*/
-void FreeSurfaceBaseAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
-}/*}}}*/
-void FreeSurfaceBaseAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
-
-	/*Now, is the model 3d? otherwise, do nothing: */
-	if (iomodel->domaintype==Domain2DhorizontalEnum)return;
-
-	/*Finite element type*/
-	int finiteelement = P1Enum;
-
-	/*Update elements: */
-	int counter=0;
-	for(int i=0;i<iomodel->numberofelements;i++){
-		if(iomodel->my_elements[i]){
-			Element* element=(Element*)elements->GetObjectByOffset(counter);
-			element->Update(i,iomodel,analysis_counter,analysis_type,finiteelement);
-			counter++;
-		}
-	}
-
-	iomodel->FetchDataToInput(elements,SurfaceEnum);
-	iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
-	iomodel->FetchDataToInput(elements,BasalforcingsGroundediceMeltingRateEnum);
-	iomodel->FetchDataToInput(elements,BasalforcingsFloatingiceMeltingRateEnum);
-	iomodel->FetchDataToInput(elements,VxEnum);
-	iomodel->FetchDataToInput(elements,VyEnum);
-	if(iomodel->domaindim==3){
-		iomodel->FetchDataToInput(elements,VzEnum);
-	}
-	if(iomodel->domaintype!=Domain2DhorizontalEnum){
-		iomodel->FetchDataToInput(elements,MeshVertexonbaseEnum);
-		iomodel->FetchDataToInput(elements,MeshVertexonsurfaceEnum);
-	}
-}/*}}}*/
-void FreeSurfaceBaseAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
-
-	if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
-	::CreateNodes(nodes,iomodel,FreeSurfaceBaseAnalysisEnum,P1Enum);
-	iomodel->DeleteData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
-}/*}}}*/
 void FreeSurfaceBaseAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
 }/*}}}*/
 void FreeSurfaceBaseAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
@@ -91,7 +48,50 @@
 	iomodel->DeleteData(vertex_pairing,MasstransportVertexPairingEnum);
 	iomodel->DeleteData(nodeonbase,MeshVertexonbaseEnum);
 }/*}}}*/
+void FreeSurfaceBaseAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
 
+	if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
+	::CreateNodes(nodes,iomodel,FreeSurfaceBaseAnalysisEnum,P1Enum);
+	iomodel->DeleteData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
+}/*}}}*/
+int  FreeSurfaceBaseAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
+	return 1;
+}/*}}}*/
+void FreeSurfaceBaseAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
+
+	/*Now, is the model 3d? otherwise, do nothing: */
+	if (iomodel->domaintype==Domain2DhorizontalEnum)return;
+
+	/*Finite element type*/
+	int finiteelement = P1Enum;
+
+	/*Update elements: */
+	int counter=0;
+	for(int i=0;i<iomodel->numberofelements;i++){
+		if(iomodel->my_elements[i]){
+			Element* element=(Element*)elements->GetObjectByOffset(counter);
+			element->Update(i,iomodel,analysis_counter,analysis_type,finiteelement);
+			counter++;
+		}
+	}
+
+	iomodel->FetchDataToInput(elements,SurfaceEnum);
+	iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
+	iomodel->FetchDataToInput(elements,BasalforcingsGroundediceMeltingRateEnum);
+	iomodel->FetchDataToInput(elements,BasalforcingsFloatingiceMeltingRateEnum);
+	iomodel->FetchDataToInput(elements,VxEnum);
+	iomodel->FetchDataToInput(elements,VyEnum);
+	if(iomodel->domaindim==3){
+		iomodel->FetchDataToInput(elements,VzEnum);
+	}
+	if(iomodel->domaintype!=Domain2DhorizontalEnum){
+		iomodel->FetchDataToInput(elements,MeshVertexonbaseEnum);
+		iomodel->FetchDataToInput(elements,MeshVertexonsurfaceEnum);
+	}
+}/*}}}*/
+void FreeSurfaceBaseAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
+}/*}}}*/
+
 /*Finite Element Analysis*/
 void           FreeSurfaceBaseAnalysis::Core(FemModel* femmodel){/*{{{*/
 	_error_("not implemented");
@@ -303,7 +303,7 @@
 	return pe;
 
 }/*}}}*/
-void FreeSurfaceBaseAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void           FreeSurfaceBaseAnalysis::GetB(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
 	 * For node i, Bi can be expressed in the actual coordinate system
 	 * by: 
@@ -331,7 +331,7 @@
 	/*Clean-up*/
 	xDelete<IssmDouble>(basis);
 }/*}}}*/
-void FreeSurfaceBaseAnalysis::GetBprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void           FreeSurfaceBaseAnalysis::GetBprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
 	 * For node i, Bi' can be expressed in the actual coordinate system
 	 * by: 
@@ -360,16 +360,16 @@
 	xDelete<IssmDouble>(dbasis);
 
 }/*}}}*/
-void FreeSurfaceBaseAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
+void           FreeSurfaceBaseAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
 	   _error_("not implemented yet");
 }/*}}}*/
-void FreeSurfaceBaseAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
+void           FreeSurfaceBaseAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
 	_error_("Not implemented yet");
 }/*}}}*/
-void FreeSurfaceBaseAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
+void           FreeSurfaceBaseAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
 	element->InputUpdateFromSolutionOneDof(solution,BaseEnum);
 }/*}}}*/
-void FreeSurfaceBaseAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
+void           FreeSurfaceBaseAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
 
 	/*Intermediary*/
 	IssmDouble phi,isonbase,base;
Index: ../trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp	(revision 18928)
+++ ../trunk-jpl/src/c/analyses/ExtrapolationAnalysis.cpp	(revision 18929)
@@ -5,15 +5,26 @@
 #include "../modules/modules.h"
 #include "../solutionsequences/solutionsequences.h"
 
-int ExtrapolationAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
-	return 1;
+void ExtrapolationAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
+	// do nothing for now
+	return;
 }
 /*}}}*/
-void ExtrapolationAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
-	//do nothing for now
+void ExtrapolationAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
+	//	do nothing for now
 	return;
+}/*}}}*/
+void ExtrapolationAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
+	int finiteelement=P1Enum;
+	if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
+	::CreateNodes(nodes,iomodel,ExtrapolationAnalysisEnum,finiteelement);
+	iomodel->DeleteData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
 }
 /*}}}*/
+int  ExtrapolationAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/
+	return 1;
+}
+/*}}}*/
 void ExtrapolationAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
 	int    finiteelement;
 
@@ -35,25 +46,14 @@
 	}
 }
 /*}}}*/
-void ExtrapolationAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
-	int finiteelement=P1Enum;
-	if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
-	::CreateNodes(nodes,iomodel,ExtrapolationAnalysisEnum,finiteelement);
-	iomodel->DeleteData(2,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum);
-}
-/*}}}*/
-void ExtrapolationAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
-	// do nothing for now
+void ExtrapolationAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
+	//do nothing for now
 	return;
 }
 /*}}}*/
-void ExtrapolationAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/
-	//	do nothing for now
-	return;
-}/*}}}*/
 
 /*Finite element Analysis*/
-void ExtrapolationAnalysis::Core(FemModel* femmodel){/*{{{*/
+void           ExtrapolationAnalysis::Core(FemModel* femmodel){/*{{{*/
 
 	/* Intermediaries */
 	bool save_results;
@@ -236,28 +236,7 @@
 	if(domaintype!=Domain2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
 	return pe;
 }/*}}}*/
-void ExtrapolationAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
-	_error_("not implemented yet");
-}/*}}}*/
-void ExtrapolationAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
-	_error_("Not implemented yet");
-}/*}}}*/
-void ExtrapolationAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
-
-	int domaintype, extrapolationvariable;
-	element->FindParam(&domaintype,DomainTypeEnum);
-	element->FindParam(&extrapolationvariable, ExtrapolationVariableEnum);
-	switch(domaintype){
-		case Domain2DhorizontalEnum:
-			element->InputUpdateFromSolutionOneDof(solution,extrapolationvariable);
-			break;
-		case Domain3DEnum:
-			element->InputUpdateFromSolutionOneDofCollapsed(solution,extrapolationvariable);
-			break;
-		default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
-	}
-}/*}}}*/
-void ExtrapolationAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void           ExtrapolationAnalysis::GetB(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
 	 * For node i, Bi can be expressed in the actual coordinate system
 	 * by: 
@@ -284,7 +263,7 @@
 	/*Clean-up*/
 	xDelete<IssmDouble>(basis);
 }/*}}}*/
-void ExtrapolationAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void           ExtrapolationAnalysis::GetBprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
 	 * For node i, Bi' can be expressed in the actual coordinate system
 	 * by: 
@@ -312,8 +291,29 @@
 	xDelete<IssmDouble>(dbasis);
 
 }/*}}}*/
-void ExtrapolationAnalysis::SetConstraintsOnIce(Element* element){/*{{{*/
+void           ExtrapolationAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
+	_error_("not implemented yet");
+}/*}}}*/
+void           ExtrapolationAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
+	_error_("Not implemented yet");
+}/*}}}*/
+void           ExtrapolationAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
 
+	int domaintype, extrapolationvariable;
+	element->FindParam(&domaintype,DomainTypeEnum);
+	element->FindParam(&extrapolationvariable, ExtrapolationVariableEnum);
+	switch(domaintype){
+		case Domain2DhorizontalEnum:
+			element->InputUpdateFromSolutionOneDof(solution,extrapolationvariable);
+			break;
+		case Domain3DEnum:
+			element->InputUpdateFromSolutionOneDofCollapsed(solution,extrapolationvariable);
+			break;
+		default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
+	}
+}/*}}}*/
+void           ExtrapolationAnalysis::SetConstraintsOnIce(Element* element){/*{{{*/
+
 	int numnodes=element->GetNumberOfNodes();	
 
 	/* Intermediaries */
@@ -346,7 +346,7 @@
 	}
 	delete gauss;
 }/*}}}*/
-void ExtrapolationAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
+void           ExtrapolationAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
 
 	for(int i=0;i<femmodel->elements->Size();i++){
 		Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
