Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 9255)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 9256)
@@ -2336,287 +2336,4 @@
 }
 /*}}}*/
-/*FUNCTION Tria::GetNodeIndex {{{1*/
-int Tria::GetNodeIndex(Node* node){
-
-	_assert_(nodes);
-	for(int i=0;i<NUMVERTICES;i++){
-		if(node==nodes[i])
-		 return i;
-	}
-	_error_("Node provided not found among element nodes");
-}
-/*}}}*/
-/*FUNCTION Tria::GetParameterListOnVertices(double* pvalue,int enumtype) {{{1*/
-void Tria::GetParameterListOnVertices(double* pvalue,int enumtype){
-
-	/*Intermediaries*/
-	double     value[NUMVERTICES];
-	GaussTria *gauss              = NULL;
-
-	/*Recover input*/
-	Input* input=inputs->GetInput(enumtype);
-	if (!input) _error_("Input %s not found in element",EnumToStringx(enumtype));
-
-	/*Checks in debugging mode*/
-	_assert_(pvalue);
-
-	/* Start looping on the number of vertices: */
-	gauss=new GaussTria();
-	for (int iv=0;iv<NUMVERTICES;iv++){
-		gauss->GaussVertex(iv);
-		input->GetParameterValue(&pvalue[iv],gauss);
-	}
-
-	/*clean-up*/
-	delete gauss;
-}
-/*}}}*/
-/*FUNCTION Tria::GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue) {{{1*/
-void Tria::GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue){
-
-	double     value[NUMVERTICES];
-	GaussTria *gauss = NULL;
-	Input     *input = inputs->GetInput(enumtype);
-
-	/*Checks in debugging mode*/
-	_assert_(pvalue);
-
-	/* Start looping on the number of vertices: */
-	if (input){
-		gauss=new GaussTria();
-		for (int iv=0;iv<NUMVERTICES;iv++){
-			gauss->GaussVertex(iv);
-			input->GetParameterValue(&pvalue[iv],gauss);
-		}
-	}
-	else{
-		for (int iv=0;iv<NUMVERTICES;iv++) pvalue[iv]=defaultvalue;
-	}
-
-	/*clean-up*/
-	delete gauss;
-}
-/*}}}*/
-/*FUNCTION Tria::GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue,int index) TO BE REMOVED{{{1*/
-void Tria::GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue,int index){
-
-	double     value[NUMVERTICES];
-	GaussTria *gauss = NULL;
-	Input     *input = inputs->GetInput(enumtype);
-
-	/*Checks in debugging mode*/
-	_assert_(pvalue);
-
-	/* Start looping on the number of vertices: */
-	if (input){
-		gauss=new GaussTria();
-		for (int iv=0;iv<NUMVERTICES;iv++){
-			gauss->GaussVertex(iv);
-			input->GetParameterValue(&pvalue[iv],gauss,index);
-		}
-	}
-	else{
-		for (int iv=0;iv<NUMVERTICES;iv++) pvalue[iv]=defaultvalue;
-	}
-
-	/*clean-up*/
-	delete gauss;
-}
-/*}}}*/
-/*FUNCTION Tria::GetParameterValue(double* pvalue,Node* node,int enumtype) {{{1*/
-void Tria::GetParameterValue(double* pvalue,Node* node,int enumtype){
-
-	Input* input=inputs->GetInput(enumtype);
-	if(!input) _error_("No input of type %s found in tria",EnumToStringx(enumtype));
-
-	GaussTria* gauss=new GaussTria();
-	gauss->GaussVertex(this->GetNodeIndex(node));
-
-	input->GetParameterValue(pvalue,gauss);
-	delete gauss;
-}
-/*}}}*/
-/*FUNCTION Tria::GetSidList {{{1*/
-void  Tria::GetSidList(int* sidlist){
-
-	int i;
-	for(i=0;i<NUMVERTICES;i++) sidlist[i]=nodes[i]->GetSidList();
-
-}
-/*}}}*/
-/*FUNCTION Tria::GetSolutionFromInputs{{{1*/
-void  Tria::GetSolutionFromInputs(Vec solution){
-
-	/*retrive parameters: */
-	int analysis_type;
-	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
-
-	/*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */
-	if (analysis_type==DiagnosticHorizAnalysisEnum)
-	 GetSolutionFromInputsDiagnosticHoriz(solution);
-	else if (analysis_type==DiagnosticHutterAnalysisEnum)
-	 GetSolutionFromInputsDiagnosticHutter(solution);
-	else if (analysis_type==HydrologyAnalysisEnum)
-	 GetSolutionFromInputsHydrology(solution);
-	else
-	 _error_("analysis: %s not supported yet",EnumToStringx(analysis_type));
-
-}
-/*}}}*/
-/*FUNCTION Tria::GetSolutionFromInputsDiagnosticHoriz{{{1*/
-void  Tria::GetSolutionFromInputsDiagnosticHoriz(Vec solution){
-
-	const int    numdof=NDOF2*NUMVERTICES;
-
-	int          i;
-	int*         doflist=NULL;
-	double       vx,vy;
-	double       values[numdof];
-	GaussTria*   gauss=NULL;
-
-	/*Get dof list: */
-	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-
-	/*Get inputs*/
-	Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
-	Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
-
-	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
-	/*P1 element only for now*/
-	gauss=new GaussTria();
-	for(i=0;i<NUMVERTICES;i++){
-
-		gauss->GaussVertex(i);
-
-		/*Recover vx and vy*/
-		vx_input->GetParameterValue(&vx,gauss);
-		vy_input->GetParameterValue(&vy,gauss);
-		values[i*NDOF2+0]=vx;
-		values[i*NDOF2+1]=vy;
-	}
-
-	VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES);
-
-	/*Free ressources:*/
-	delete gauss;
-	xfree((void**)&doflist);
-}
-/*}}}*/
-/*FUNCTION Tria::GetSolutionFromInputsDiagnosticHutter{{{1*/
-void  Tria::GetSolutionFromInputsDiagnosticHutter(Vec solution){
-
-	const int    numdof=NDOF2*NUMVERTICES;
-
-	int i,dummy;
-	int*         doflist=NULL;
-	double       vx,vy;
-	double       values[numdof];
-	GaussTria*   gauss=NULL;
-
-	/*Get dof list: */
-	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-
-	/*Get inputs*/
-	Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
-	Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
-
-	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
-	/*P1 element only for now*/
-	gauss=new GaussTria();
-	for(i=0;i<NUMVERTICES;i++){
-
-		gauss->GaussVertex(i);
-
-		/*Recover vx and vy*/
-		vx_input->GetParameterValue(&vx,gauss);
-		vy_input->GetParameterValue(&vy,gauss);
-		values[i*NDOF2+0]=vx;
-		values[i*NDOF2+1]=vy;
-	}
-
-	VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES);
-
-	/*Free ressources:*/
-	delete gauss;
-	xfree((void**)&doflist);
-}
-/*}}}*/
-/*FUNCTION Tria::GetSolutionFromInputsHydrology{{{1*/
-void  Tria::GetSolutionFromInputsHydrology(Vec solution){
-
-	const int    numdof=NDOF1*NUMVERTICES;
-
-	int i,dummy;
-	int*         doflist=NULL;
-	double       watercolumn;
-	double       values[numdof];
-	GaussTria*   gauss=NULL;
-
-	/*Get dof list: */
-	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-
-	/*Get inputs*/
-	Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); _assert_(watercolumn_input);
-
-	/*Ok, we have watercolumn values, fill in watercolumn array: */
-	/*P1 element only for now*/
-	gauss=new GaussTria();
-	for(i=0;i<NUMVERTICES;i++){
-
-		gauss->GaussVertex(i);
-
-		/*Recover watercolumn*/
-		watercolumn_input->GetParameterValue(&watercolumn,gauss);
-		values[i]=watercolumn;
-	}
-
-	VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES);
-
-	/*Free ressources:*/
-	delete gauss;
-	xfree((void**)&doflist);
-}
-/*}}}*/
-/*FUNCTION Tria::GetStrainRate2d(double* epsilon,double* xyz_list, GaussTria* gauss, Input* vx_input, Input* vy_input){{{1*/
-void Tria::GetStrainRate2d(double* epsilon,double* xyz_list, GaussTria* gauss, Input* vx_input, Input* vy_input){
-	/*Compute the 2d Strain Rate (3 components):
-	 * epsilon=[exx eyy exy] */
-
-	int i;
-	double epsilonvx[3];
-	double epsilonvy[3];
-
-	/*Check that both inputs have been found*/
-	if (!vx_input || !vy_input){
-		_error_("Input missing. Here are the input pointers we have for vx: %p, vy: %p\n",vx_input,vy_input);
-	}
-
-	/*Get strain rate assuming that epsilon has been allocated*/
-	vx_input->GetVxStrainRate2d(epsilonvx,xyz_list,gauss);
-	vy_input->GetVyStrainRate2d(epsilonvy,xyz_list,gauss);
-
-	/*Sum all contributions*/
-	for(i=0;i<3;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i];
-}
-/*}}}*/
-/*FUNCTION Tria::GetVectorFromInputs{{{1*/
-void  Tria::GetVectorFromInputs(Vec vector,int input_enum){
-
-	int doflist1[NUMVERTICES];
-
-	/*Get out if this is not an element input*/
-	if (!IsInput(input_enum)) return;
-
-	/*Prepare index list*/
-	this->GetDofList1(&doflist1[0]);
-
-	/*Get input (either in element or material)*/
-	Input* input=inputs->GetInput(input_enum);
-	if(!input) _error_("Input %s not found in element",EnumToStringx(input_enum));
-
-	/*We found the enum.  Use its values to fill into the vector, using the vertices ids: */
-	input->GetVectorFromInputs(vector,&doflist1[0]);
-}
-/*}}}*/
 /*FUNCTION Tria::GetHydrologyK {{{1*/
 void Tria::GetHydrologyK(double* K,double* xyz_list,GaussTria* gauss){
@@ -2689,4 +2406,287 @@
 	v[0]=pow(w,2)/(12 * mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx);
 	v[1]=pow(w,2)/(12 * mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy);
+}
+/*}}}*/
+/*FUNCTION Tria::GetNodeIndex {{{1*/
+int Tria::GetNodeIndex(Node* node){
+
+	_assert_(nodes);
+	for(int i=0;i<NUMVERTICES;i++){
+		if(node==nodes[i])
+		 return i;
+	}
+	_error_("Node provided not found among element nodes");
+}
+/*}}}*/
+/*FUNCTION Tria::GetParameterListOnVertices(double* pvalue,int enumtype) {{{1*/
+void Tria::GetParameterListOnVertices(double* pvalue,int enumtype){
+
+	/*Intermediaries*/
+	double     value[NUMVERTICES];
+	GaussTria *gauss              = NULL;
+
+	/*Recover input*/
+	Input* input=inputs->GetInput(enumtype);
+	if (!input) _error_("Input %s not found in element",EnumToStringx(enumtype));
+
+	/*Checks in debugging mode*/
+	_assert_(pvalue);
+
+	/* Start looping on the number of vertices: */
+	gauss=new GaussTria();
+	for (int iv=0;iv<NUMVERTICES;iv++){
+		gauss->GaussVertex(iv);
+		input->GetParameterValue(&pvalue[iv],gauss);
+	}
+
+	/*clean-up*/
+	delete gauss;
+}
+/*}}}*/
+/*FUNCTION Tria::GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue) {{{1*/
+void Tria::GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue){
+
+	double     value[NUMVERTICES];
+	GaussTria *gauss = NULL;
+	Input     *input = inputs->GetInput(enumtype);
+
+	/*Checks in debugging mode*/
+	_assert_(pvalue);
+
+	/* Start looping on the number of vertices: */
+	if (input){
+		gauss=new GaussTria();
+		for (int iv=0;iv<NUMVERTICES;iv++){
+			gauss->GaussVertex(iv);
+			input->GetParameterValue(&pvalue[iv],gauss);
+		}
+	}
+	else{
+		for (int iv=0;iv<NUMVERTICES;iv++) pvalue[iv]=defaultvalue;
+	}
+
+	/*clean-up*/
+	delete gauss;
+}
+/*}}}*/
+/*FUNCTION Tria::GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue,int index) TO BE REMOVED{{{1*/
+void Tria::GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue,int index){
+
+	double     value[NUMVERTICES];
+	GaussTria *gauss = NULL;
+	Input     *input = inputs->GetInput(enumtype);
+
+	/*Checks in debugging mode*/
+	_assert_(pvalue);
+
+	/* Start looping on the number of vertices: */
+	if (input){
+		gauss=new GaussTria();
+		for (int iv=0;iv<NUMVERTICES;iv++){
+			gauss->GaussVertex(iv);
+			input->GetParameterValue(&pvalue[iv],gauss,index);
+		}
+	}
+	else{
+		for (int iv=0;iv<NUMVERTICES;iv++) pvalue[iv]=defaultvalue;
+	}
+
+	/*clean-up*/
+	delete gauss;
+}
+/*}}}*/
+/*FUNCTION Tria::GetParameterValue(double* pvalue,Node* node,int enumtype) {{{1*/
+void Tria::GetParameterValue(double* pvalue,Node* node,int enumtype){
+
+	Input* input=inputs->GetInput(enumtype);
+	if(!input) _error_("No input of type %s found in tria",EnumToStringx(enumtype));
+
+	GaussTria* gauss=new GaussTria();
+	gauss->GaussVertex(this->GetNodeIndex(node));
+
+	input->GetParameterValue(pvalue,gauss);
+	delete gauss;
+}
+/*}}}*/
+/*FUNCTION Tria::GetSidList {{{1*/
+void  Tria::GetSidList(int* sidlist){
+
+	int i;
+	for(i=0;i<NUMVERTICES;i++) sidlist[i]=nodes[i]->GetSidList();
+
+}
+/*}}}*/
+/*FUNCTION Tria::GetSolutionFromInputs{{{1*/
+void  Tria::GetSolutionFromInputs(Vec solution){
+
+	/*retrive parameters: */
+	int analysis_type;
+	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+
+	/*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */
+	if (analysis_type==DiagnosticHorizAnalysisEnum)
+	 GetSolutionFromInputsDiagnosticHoriz(solution);
+	else if (analysis_type==DiagnosticHutterAnalysisEnum)
+	 GetSolutionFromInputsDiagnosticHutter(solution);
+	else if (analysis_type==HydrologyAnalysisEnum)
+	 GetSolutionFromInputsHydrology(solution);
+	else
+	 _error_("analysis: %s not supported yet",EnumToStringx(analysis_type));
+
+}
+/*}}}*/
+/*FUNCTION Tria::GetSolutionFromInputsDiagnosticHoriz{{{1*/
+void  Tria::GetSolutionFromInputsDiagnosticHoriz(Vec solution){
+
+	const int    numdof=NDOF2*NUMVERTICES;
+
+	int          i;
+	int*         doflist=NULL;
+	double       vx,vy;
+	double       values[numdof];
+	GaussTria*   gauss=NULL;
+
+	/*Get dof list: */
+	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
+
+	/*Get inputs*/
+	Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
+	Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
+
+	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
+	/*P1 element only for now*/
+	gauss=new GaussTria();
+	for(i=0;i<NUMVERTICES;i++){
+
+		gauss->GaussVertex(i);
+
+		/*Recover vx and vy*/
+		vx_input->GetParameterValue(&vx,gauss);
+		vy_input->GetParameterValue(&vy,gauss);
+		values[i*NDOF2+0]=vx;
+		values[i*NDOF2+1]=vy;
+	}
+
+	VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES);
+
+	/*Free ressources:*/
+	delete gauss;
+	xfree((void**)&doflist);
+}
+/*}}}*/
+/*FUNCTION Tria::GetSolutionFromInputsDiagnosticHutter{{{1*/
+void  Tria::GetSolutionFromInputsDiagnosticHutter(Vec solution){
+
+	const int    numdof=NDOF2*NUMVERTICES;
+
+	int i,dummy;
+	int*         doflist=NULL;
+	double       vx,vy;
+	double       values[numdof];
+	GaussTria*   gauss=NULL;
+
+	/*Get dof list: */
+	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
+
+	/*Get inputs*/
+	Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
+	Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
+
+	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
+	/*P1 element only for now*/
+	gauss=new GaussTria();
+	for(i=0;i<NUMVERTICES;i++){
+
+		gauss->GaussVertex(i);
+
+		/*Recover vx and vy*/
+		vx_input->GetParameterValue(&vx,gauss);
+		vy_input->GetParameterValue(&vy,gauss);
+		values[i*NDOF2+0]=vx;
+		values[i*NDOF2+1]=vy;
+	}
+
+	VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES);
+
+	/*Free ressources:*/
+	delete gauss;
+	xfree((void**)&doflist);
+}
+/*}}}*/
+/*FUNCTION Tria::GetSolutionFromInputsHydrology{{{1*/
+void  Tria::GetSolutionFromInputsHydrology(Vec solution){
+
+	const int    numdof=NDOF1*NUMVERTICES;
+
+	int i,dummy;
+	int*         doflist=NULL;
+	double       watercolumn;
+	double       values[numdof];
+	GaussTria*   gauss=NULL;
+
+	/*Get dof list: */
+	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
+
+	/*Get inputs*/
+	Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); _assert_(watercolumn_input);
+
+	/*Ok, we have watercolumn values, fill in watercolumn array: */
+	/*P1 element only for now*/
+	gauss=new GaussTria();
+	for(i=0;i<NUMVERTICES;i++){
+
+		gauss->GaussVertex(i);
+
+		/*Recover watercolumn*/
+		watercolumn_input->GetParameterValue(&watercolumn,gauss);
+		values[i]=watercolumn;
+	}
+
+	VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES);
+
+	/*Free ressources:*/
+	delete gauss;
+	xfree((void**)&doflist);
+}
+/*}}}*/
+/*FUNCTION Tria::GetStrainRate2d(double* epsilon,double* xyz_list, GaussTria* gauss, Input* vx_input, Input* vy_input){{{1*/
+void Tria::GetStrainRate2d(double* epsilon,double* xyz_list, GaussTria* gauss, Input* vx_input, Input* vy_input){
+	/*Compute the 2d Strain Rate (3 components):
+	 * epsilon=[exx eyy exy] */
+
+	int i;
+	double epsilonvx[3];
+	double epsilonvy[3];
+
+	/*Check that both inputs have been found*/
+	if (!vx_input || !vy_input){
+		_error_("Input missing. Here are the input pointers we have for vx: %p, vy: %p\n",vx_input,vy_input);
+	}
+
+	/*Get strain rate assuming that epsilon has been allocated*/
+	vx_input->GetVxStrainRate2d(epsilonvx,xyz_list,gauss);
+	vy_input->GetVyStrainRate2d(epsilonvy,xyz_list,gauss);
+
+	/*Sum all contributions*/
+	for(i=0;i<3;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i];
+}
+/*}}}*/
+/*FUNCTION Tria::GetVectorFromInputs{{{1*/
+void  Tria::GetVectorFromInputs(Vec vector,int input_enum){
+
+	int doflist1[NUMVERTICES];
+
+	/*Get out if this is not an element input*/
+	if (!IsInput(input_enum)) return;
+
+	/*Prepare index list*/
+	this->GetDofList1(&doflist1[0]);
+
+	/*Get input (either in element or material)*/
+	Input* input=inputs->GetInput(input_enum);
+	if(!input) _error_("Input %s not found in element",EnumToStringx(input_enum));
+
+	/*We found the enum.  Use its values to fill into the vector, using the vertices ids: */
+	input->GetVectorFromInputs(vector,&doflist1[0]);
 }
 /*}}}*/
Index: /issm/trunk/src/c/objects/Elements/Tria.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.h	(revision 9255)
+++ /issm/trunk/src/c/objects/Elements/Tria.h	(revision 9256)
@@ -181,4 +181,5 @@
 		int     GetElementType(void);
 		void	  GetDofList(int** pdoflist,int approximation_enum,int setenum);
+		void    GetHydrologyWaterVelocity(double* v, GaussTria* gauss);
 		void	  GetDofList1(int* doflist);
 		void    GetSidList(int* sidlist);