Index: /issm/trunk/src/c/objects/Elements/Element.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Element.h	(revision 5632)
+++ /issm/trunk/src/c/objects/Elements/Element.h	(revision 5633)
@@ -35,8 +35,8 @@
 		virtual bool   GetShelf()=0; 
 		virtual bool   GetOnBed()=0;
-		virtual void   GetThicknessList(double* thickness_list)=0;
+		virtual void   GetParameterListOnVertices(double* pvalue,int enumtype)=0;
+		virtual void   GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue)=0;
 		virtual void   GetParameterValue(double* pvalue,Node* node,int enumtype)=0;
 		virtual void   GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,Input* input_in)=0;
-		virtual void   GetBedList(double* bed_list)=0;
 		virtual void   Gradj(Vec gradient,int control_type)=0;
 		virtual void   GradjDrag(Vec gradient)=0;
Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5632)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5633)
@@ -943,14 +943,4 @@
 }
 /*}}}*/
-/*FUNCTION Penta::GetBedList {{{1*/
-void Penta::GetBedList(double* bedlist){
-
-	const int numgrids=6;
-	double  gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
-	
-	inputs->GetParameterValues(bedlist,&gaussgrids[0][0],6,BedEnum);
-
-}
-/*}}}*/
 /*FUNCTION Penta::GetMatPar {{{1*/
 void* Penta::GetMatPar(){
@@ -1024,13 +1014,4 @@
 		ISSMERROR("analysis: %i (%s) not supported yet",analysis_type,EnumToString(analysis_type));
 	}
-}
-/*}}}*/
-/*FUNCTION Penta::GetThicknessList {{{1*/
-void Penta::GetThicknessList(double* thickness_list){
-
-	const int numgrids=6;
-	double  gaussgrids[numgrids][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
-	inputs->GetParameterValues(thickness_list,&gaussgrids[0][0],6,ThicknessEnum);
-
 }
 /*}}}*/
@@ -4835,4 +4816,51 @@
 }
 /*}}}*/
+/*FUNCTION Penta::GetParameterListOnVertices(double* pvalue,int enumtype) {{{1*/
+void Penta::GetParameterListOnVertices(double* pvalue,int enumtype){
+
+	/*Intermediaries*/
+	const int  numvertices = 6;
+	double     value[numvertices];
+	double     gauss[numvertices][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
+
+	/*Recover input*/
+	Input* input=inputs->GetInput(enumtype);
+	if (!input) ISSMERROR("Input %s not found in element",EnumToString(enumtype));
+
+	/*Checks in debugging mode*/
+	ISSMASSERT(pvalue);
+
+	/* Start looping on the number of vertices: */
+	for (int iv=0;iv<numvertices;iv++){
+		input->GetParameterValue(&pvalue[iv],&gauss[iv][0]);
+	}
+
+	/*clean-up*/
+}
+/*}}}*/
+/*FUNCTION Penta::GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue) {{{1*/
+void Penta::GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue){
+
+	/*Intermediaries*/
+	const int  numvertices = 6;
+	double     value[numvertices];
+	double     gauss[numvertices][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
+
+	/*Recover input*/
+	Input* input=inputs->GetInput(enumtype);
+
+	/*Checks in debugging mode*/
+	ISSMASSERT(pvalue);
+
+	/* Start looping on the number of vertices: */
+	if (input)
+	 for (int iv=0;iv<numvertices;iv++) input->GetParameterValue(&pvalue[iv],&gauss[iv][0]);
+	else
+	 for (int iv=0;iv<numvertices;iv++) pvalue[iv]=defaultvalue;
+
+	/*clean-up*/
+	delete gauss;
+}
+/*}}}*/
 /*FUNCTION Penta::GetParameterValue(double* pvalue,Node* node,int enumtype) {{{1*/
 void Penta::GetParameterValue(double* pvalue,Node* node,int enumtype){
Index: /issm/trunk/src/c/objects/Elements/Penta.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.h	(revision 5632)
+++ /issm/trunk/src/c/objects/Elements/Penta.h	(revision 5633)
@@ -79,5 +79,4 @@
 		void   CreatePVector(Vec pg);
 		void   DeleteResults(void);
-		void   GetBedList(double* bed_list);
 		void*  GetMatPar();
 		void   GetNodes(void** nodes);
@@ -85,5 +84,4 @@
 		bool   GetShelf(); 
 		void   GetSolutionFromInputs(Vec solution);
-		void   GetThicknessList(double* thickness_list);
 		void   GetVectorFromInputs(Vec vector,int NameEnum);
 		void   Gradj(Vec gradient,int control_type);
@@ -160,4 +158,6 @@
 		void	  GetNodalFunctionsDerivativesStokes(double* dh1dh7,double* xyz_list, double* gauss_coord);
 		void	  GetNodalFunctionsStokes(double* l1l7, double* gauss_coord);
+		void    GetParameterListOnVertices(double* pvalue,int enumtype);
+		void    GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue);
 		void      GetParameterValue(double* pvalue,Node* node,int enumtype);
 		void      GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,Input* input_in);
Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 5632)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 5633)
@@ -450,5 +450,4 @@
 	
 	const int numvertices=3;
-	double  gaussgrids[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
 	int i,j;
 
@@ -477,6 +476,6 @@
 					
 					/*retrieve inputs: */
-					inputs->GetParameterValues(&bed[0],&gaussgrids[0][0],3,BedEnum);
-					inputs->GetParameterValues(&surface[0],&gaussgrids[0][0],3,SurfaceEnum);
+					GetParameterListOnVertices(&bed[0],BedEnum);
+					GetParameterListOnVertices(&surface[0],SurfaceEnum);
 
 					/*build new thickness: */
@@ -594,5 +593,4 @@
 	double thickness[numvertices];
 	double rho_ice,g;
-	double gauss[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
 	
 	/*Get dof list on which we will plug the pressure values: */
@@ -602,11 +600,7 @@
 	rho_ice=matpar->GetRhoIce();
 	g=matpar->GetG();
-
-	/*recover value of thickness at gauss points (0,0,1), (0,1,0),(1,0,0): */
-	inputs->GetParameterValues(&thickness[0],&gauss[0][0],3,ThicknessEnum);
-
-	for(i=0;i<numvertices;i++){
-		pressure[i]=rho_ice*g*thickness[i];
-	}
+	GetParameterListOnVertices(&thickness[0],ThicknessEnum);
+
+	for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*thickness[i];
 
 	/*plug local pressure values into global pressure vector: */
@@ -849,29 +843,4 @@
 }
 /*}}}*/
-/*FUNCTION Tria::GetBedList {{{1*/
-void  Tria::GetBedList(double* bedlist){
-
-	Input     *bed_input = NULL;
-	GaussTria *gauss     = NULL;
-
-	/*Recover input*/
-	bed_input=inputs->GetInput(BedEnum);
-
-	/*Checks in debugging mode*/
-	ISSMASSERT(bedlist);
-	ISSMASSERT(bed_input);
-
-	/* Start looping on the number of vertices: */
-	gauss=new GaussTria();
-	for (int iv=0;iv<3;iv++){
-		gauss->GaussVertex(iv);
-		bed_input->GetParameterValue(&bedlist[iv],gauss);
-	}
-
-	/*clean-up*/
-	delete gauss;
-
-}
-/*}}}*/
 /*FUNCTION Tria::GetMatPar {{{1*/
 void* Tria::GetMatPar(){
@@ -923,35 +892,8 @@
 	if (analysis_type==DiagnosticHorizAnalysisEnum)
 	 GetSolutionFromInputsDiagnosticHoriz(solution);
-	else if (analysis_type==AdjointHorizAnalysisEnum)
-	 GetSolutionFromInputsAdjointHoriz(solution);
 	else if (analysis_type==DiagnosticHutterAnalysisEnum)
 	 GetSolutionFromInputsDiagnosticHutter(solution);
 	else
 	 ISSMERROR("analysis: %s not supported yet",EnumToString(analysis_type));
-
-}
-/*}}}*/
-/*FUNCTION Tria::GetThicknessList {{{1*/
-void Tria::GetThicknessList(double* thicknesslist){
-
-	Input     *thickness_input = NULL;
-	GaussTria *gauss     = NULL;
-
-	/*Recover input*/
-	thickness_input=inputs->GetInput(ThicknessEnum);
-
-	/*Checks in debugging mode*/
-	ISSMASSERT(thicknesslist);
-	ISSMASSERT(thickness_input);
-
-	/* Start looping on the number of vertices: */
-	gauss=new GaussTria();
-	for (int iv=0;iv<3;iv++){
-		gauss->GaussVertex(iv);
-		thickness_input->GetParameterValue(&thicknesslist[iv],gauss);
-	}
-
-	/*clean-up*/
-	delete gauss;
 
 }
@@ -6338,4 +6280,60 @@
 }
 /*}}}*/
+/*FUNCTION Tria::GetParameterListOnVertices(double* pvalue,int enumtype) {{{1*/
+void Tria::GetParameterListOnVertices(double* pvalue,int enumtype){
+
+	/*Intermediaries*/
+	const int  numvertices        = 3;
+	double     value[numvertices];
+	GaussTria *gauss              = NULL;
+
+	/*Recover input*/
+	Input* input=inputs->GetInput(enumtype);
+	if (!input) ISSMERROR("Input %s not found in element",EnumToString(enumtype));
+
+	/*Checks in debugging mode*/
+	ISSMASSERT(pvalue);
+
+	/* Start looping on the number of vertices: */
+	gauss=new GaussTria();
+	for (int iv=0;iv<3;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){
+
+	/*Intermediaries*/
+	const int  numvertices        = 3;
+	double     value[numvertices];
+	GaussTria *gauss              = NULL;
+
+	/*Recover input*/
+	Input* input=inputs->GetInput(enumtype);
+
+	/*Checks in debugging mode*/
+	ISSMASSERT(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::GetParameterValue(double* pvalue,Node* node,int enumtype) {{{1*/
 void Tria::GetParameterValue(double* pvalue,Node* node,int enumtype){
@@ -6495,24 +6493,30 @@
 
 	int i;
-
 	const int    numvertices=3;
 	const int    numdofpervertex=2;
 	const int    numdof=numdofpervertex*numvertices;
-	double       gauss[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
 	int*         doflist=NULL;
 	double       values[numdof];
 	double       vx;
 	double       vy;
+	GaussTria*   gauss=NULL;
 
 	/*Get dof list: */
 	GetDofList(&doflist);
 
+	/*Get inputs*/
+	Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input);
+	Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(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*/
-		inputs->GetParameterValue(&vx,&gauss[i][0],VxEnum);
-		inputs->GetParameterValue(&vy,&gauss[i][0],VyEnum);
+		vx_input->GetParameterValue(&vx,gauss);
+		vy_input->GetParameterValue(&vy,gauss);
 		values[i*numdofpervertex+0]=vx;
 		values[i*numdofpervertex+1]=vy;
@@ -6523,53 +6527,16 @@
 
 	/*Free ressources:*/
+	delete gauss;
 	xfree((void**)&doflist);
 
 }
 /*}}}*/
-/*FUNCTION Tria::GetSolutionFromInputsAdjointHoriz{{{1*/
-void  Tria::GetSolutionFromInputsAdjointHoriz(Vec solution){
+/*FUNCTION Tria::GetSolutionFromInputsDiagnosticHutter{{{1*/
+void  Tria::GetSolutionFromInputsDiagnosticHutter(Vec solution){
 
 	int i;
-
 	const int    numvertices=3;
 	const int    numdofpervertex=2;
 	const int    numdof=numdofpervertex*numvertices;
-	double       gauss[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
-	int*         doflist=NULL;
-	double       values[numdof];
-	double       lambdax;
-	double       lambday;
-
-	/*Get dof list: */
-	GetDofList(&doflist);
-
-	/*Ok, we have lambdax and lambday in values, fill in lambdax and lambday arrays: */
-	/*P1 element only for now*/
-	for(i=0;i<numvertices;i++){
-
-		/*Recover lambdax and lambday*/
-		inputs->GetParameterValue(&lambdax,&gauss[i][0],VxEnum);
-		inputs->GetParameterValue(&lambday,&gauss[i][0],VyEnum);
-		values[i*numdofpervertex+0]=lambdax;
-		values[i*numdofpervertex+1]=lambday;
-	}
-
-	/*Add value to global vector*/
-	VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES);
-
-	/*Free ressources:*/
-	xfree((void**)&doflist);
-
-}
-/*}}}*/
-/*FUNCTION Tria::GetSolutionFromInputsDiagnosticHutter{{{1*/
-void  Tria::GetSolutionFromInputsDiagnosticHutter(Vec solution){
-
-	int i;
-
-	const int    numvertices=3;
-	const int    numdofpervertex=2;
-	const int    numdof=numdofpervertex*numvertices;
-	double       gauss[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
 	int*         doflist=NULL;
 	double       values[numdof];
@@ -6577,15 +6544,23 @@
 	double       vy;
 	int          dummy;
+	GaussTria*   gauss=NULL;
 
 	/*Get dof list: */
 	GetDofList(&doflist);
 
+	/*Get inputs*/
+	Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input);
+	Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(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*/
-		inputs->GetParameterValue(&vx,&gauss[i][0],VxEnum);
-		inputs->GetParameterValue(&vy,&gauss[i][0],VyEnum);
+		vx_input->GetParameterValue(&vx,gauss);
+		vy_input->GetParameterValue(&vy,gauss);
 		values[i*numdofpervertex+0]=vx;
 		values[i*numdofpervertex+1]=vy;
@@ -6596,4 +6571,5 @@
 
 	/*Free ressources:*/
+	delete gauss;
 	xfree((void**)&doflist);
 
@@ -6892,7 +6868,4 @@
 	double       thickness[numvertices];
 	double       rho_ice,g;
-	double       gauss[numvertices][numvertices]={{1,0,0},{0,1,0},{0,0,1}};
-	Input*       vz_input=NULL;
-	double*      vz_ptr=NULL;
 	int          dummy;
 	
@@ -6912,15 +6885,5 @@
 
 	/*Get Vz*/
-	vz_input=inputs->GetInput(VzEnum);
-	if (vz_input){
-		if (vz_input->Enum()!=TriaVertexInputEnum){
-			ISSMERROR("Cannot compute Vel as Vz is of type %s",EnumToString(vz_input->Enum()));
-		}
-		vz_input->GetValuesPtr(&vz_ptr,&dummy);
-		for(i=0;i<numvertices;i++) vz[i]=vz_ptr[i];
-	}
-	else{
-		for(i=0;i<numvertices;i++) vz[i]=0.0;
-	}
+	GetParameterListOnVertices(&vz[0],VzEnum,0);
 
 	/*Now Compute vel*/
@@ -6931,9 +6894,6 @@
 	rho_ice=matpar->GetRhoIce();
 	g=matpar->GetG();
-	inputs->GetParameterValues(&thickness[0],&gauss[0][0],3,ThicknessEnum);
-	
-	for(i=0;i<numvertices;i++){
-		pressure[i]=rho_ice*g*thickness[i];
-	}
+	GetParameterListOnVertices(&thickness[0],ThicknessEnum);
+	for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*thickness[i];
 
 	/*Now, we have to move the previous Vx and Vy inputs  to old 
Index: /issm/trunk/src/c/objects/Elements/Tria.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.h	(revision 5632)
+++ /issm/trunk/src/c/objects/Elements/Tria.h	(revision 5633)
@@ -75,5 +75,4 @@
 		void   CreateKMatrix(Mat Kgg);
 		void   CreatePVector(Vec pg);
-		void   GetBedList(double* bed_list);
 		void*  GetMatPar();
 		void   GetNodes(void** nodes);
@@ -81,5 +80,4 @@
 		bool   GetShelf(); 
 		void   GetSolutionFromInputs(Vec solution);
-		void   GetThicknessList(double* thickness_list);
 		void   GetVectorFromInputs(Vec vector,int NameEnum);
 		void   Gradj(Vec gradient,int control_type);
@@ -155,9 +153,10 @@
 		void	  GetDofList(int** pdoflist,int approximation_enum=0);
 		void	  GetDofList1(int* doflist);
+		void    GetParameterListOnVertices(double* pvalue,int enumtype);
+		void    GetParameterListOnVertices(double* pvalue,int enumtype,double defaultvalue);
 		void      GetParameterValue(double* pvalue,Node* node,int enumtype);
 		void      GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,int enumtype);
 		void      GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,Input* input_in);
 		void	  GetSolutionFromInputsDiagnosticHoriz(Vec solution);
-		void	  GetSolutionFromInputsAdjointHoriz(Vec solution);
 		void	  GetSolutionFromInputsDiagnosticHutter(Vec solution);
 		void    GetStrainRate2d(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input);
Index: /issm/trunk/src/c/objects/Loads/Icefront.cpp
===================================================================
--- /issm/trunk/src/c/objects/Loads/Icefront.cpp	(revision 5632)
+++ /issm/trunk/src/c/objects/Loads/Icefront.cpp	(revision 5633)
@@ -592,6 +592,6 @@
 	}
 
-	element->GetThicknessList(&thickness_list[0]);
-	element->GetBedList(&bed_list[0]);
+	element->GetParameterListOnVertices(&thickness_list[0],ThicknessEnum);
+	element->GetParameterListOnVertices(&bed_list[0],BedEnum);
 
 	thickness_list_quad[0]=thickness_list[grid1];
@@ -731,6 +731,6 @@
 	}
 
-	element->GetThicknessList(&thickness_list[0]);
-	element->GetBedList(&bed_list[0]);
+	element->GetParameterListOnVertices(&thickness_list[0],ThicknessEnum);
+	element->GetParameterListOnVertices(&bed_list[0],BedEnum);
 
 	thickness_list_quad[0]=thickness_list[grid1];