Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 24433)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 24434)
@@ -1108,39 +1108,14 @@
 /*}}}*/
 void       Element::GetInputListOnNodes(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue){/*{{{*/
-
-	_assert_(pvalue);
-
 	Input2 *input    = this->GetInput2(enumtype);
-	int    numnodes = this->GetNumberOfNodes();
-
-	/* Start looping on the number of vertices: */
-	if(input){
-		Gauss* gauss=this->NewGauss();
-		for(int iv=0;iv<numnodes;iv++){
-			gauss->GaussNode(this->FiniteElement(),iv);
-			input->GetInputValue(&pvalue[iv],gauss);
-		}
-		delete gauss;
-	}
-	else{
-		for(int iv=0;iv<numnodes;iv++) pvalue[iv]=defaultvalue;
-	}
+	this->GetInputListOnNodes(pvalue,input,defaultvalue);
 }
 /*}}}*/
 void       Element::GetInputListOnNodes(IssmDouble* pvalue,int enumtype){/*{{{*/
 
-	_assert_(pvalue);
-
-	int     numnodes = this->GetNumberOfNodes();
 	Input2 *input    = this->GetInput2(enumtype);
 	if(!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element");
-
-	/* Start looping on the number of vertices: */
-	Gauss* gauss=this->NewGauss();
-	for(int iv=0;iv<numnodes;iv++){
-		gauss->GaussNode(this->FiniteElement(),iv);
-		input->GetInputValue(&pvalue[iv],gauss);
-	}
-	delete gauss;
+	this->GetInputListOnNodes(pvalue,input,0.);
+
 }
 /*}}}*/
@@ -1162,40 +1137,10 @@
 }
 /*}}}*/
-void       Element::GetInputListOnVertices(IssmDouble* pvalue,ElementInput2* input){/*{{{*/
-
-	/*Fetch number vertices for this element*/
-	const int NUM_VERTICES = this->GetNumberOfVertices();
-
-	/*Checks in debugging mode*/
-	_assert_(input);
-	_assert_(pvalue);
-
-	/* Start looping on the number of vertices: */
-	Gauss*gauss=this->NewGauss();
-	for(int iv=0;iv<NUM_VERTICES;iv++){
-		gauss->GaussVertex(iv);
-		input->GetInputValue(&pvalue[iv],gauss);
-	}
-	delete gauss;
-}
-/*}}}*/
 void       Element::GetInputListOnVertices(IssmDouble* pvalue,int enumtype){/*{{{*/
-
-	/*Fetch number vertices for this element*/
-	const int NUM_VERTICES = this->GetNumberOfVertices();
-
-	/*Checks in debugging mode*/
-	_assert_(pvalue);
 
 	/*Recover input*/
 	Input2* input2=this->GetInput2(enumtype);
 	if(!input2) _error_("input "<<EnumToStringx(enumtype)<<" not found in element");
-	/* Start looping on the number of vertices: */
-	Gauss*gauss=this->NewGauss();
-	for(int iv=0;iv<NUM_VERTICES;iv++){
-		gauss->GaussVertex(iv);
-		input2->GetInputValue(&pvalue[iv],gauss);
-	}
-	delete gauss;
+	this->GetInputListOnVertices(pvalue,input2,0.);
 }
 /*}}}*/
@@ -1205,45 +1150,10 @@
 	Input2* input=this->GetInput2(enumtype,time);
 	if (!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element");
-
-	/*Fetch number vertices for this element*/
-	const int NUM_VERTICES = this->GetNumberOfVertices();
-
-	/*Checks in debugging mode*/
-	_assert_(pvalue);
-
-	/* Start looping on the number of vertices: */
-	Gauss*gauss=this->NewGauss();
-	for(int iv=0;iv<NUM_VERTICES;iv++){
-		gauss->GaussVertex(iv);
-		input->GetInputValue(&pvalue[iv],gauss);
-	}
-
-	/*clean-up*/
-	delete gauss;
+	this->GetInputListOnVertices(pvalue,input,0.);
 }
 /*}}}*/
 void       Element::GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue){/*{{{*/
-
-	/*Recover input*/
 	Input2* input=this->GetInput2(enumtype);
-
-	/*Checks in debugging mode*/
-	_assert_(pvalue);
-
-	/*Fetch number vertices for this element*/
-	const int NUM_VERTICES = this->GetNumberOfVertices();
-
-	/* Start looping on the number of vertices: */
-	if (input){
-		Gauss* gauss=this->NewGauss();
-		for (int iv=0;iv<NUM_VERTICES;iv++){
-			gauss->GaussVertex(iv);
-			input->GetInputValue(&pvalue[iv],gauss);
-		}
-		delete gauss;
-	}
-	else{
-		for(int iv=0;iv<NUM_VERTICES;iv++) pvalue[iv]=defaultvalue;
-	}
+	this->GetInputListOnVertices(pvalue,input,defaultvalue);
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 24433)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 24434)
@@ -94,5 +94,4 @@
 		void               GetInputListOnNodes(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue);
 		void               GetInputListOnNodesVelocity(IssmDouble* pvalue,int enumtype);
-		void               GetInputListOnVertices(IssmDouble* pvalue,ElementInput2* input);
 		void               GetInputListOnVertices(IssmDouble* pvalue,int enumtype);
 		void               GetInputListOnVerticesAtTime(IssmDouble* pvalue,int enumtype,IssmDouble time);
@@ -253,4 +252,6 @@
 		virtual void       GetInputValue(IssmDouble* pvalue,Vertex* vertex,int enumtype){_error_("not implemented yet");};
 		virtual void       GetInputValue(IssmDouble* pvalue,Node* node,int enumtype){_error_("not implemented yet");};
+		virtual void       GetInputListOnVertices(IssmDouble* pvalue,Input2* input,IssmDouble default_value)=0;
+		virtual void       GetInputListOnNodes(IssmDouble* pvalue,Input2* input,IssmDouble default_value)=0;
 		virtual void       GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level)=0;
 		virtual void       GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* levelsetvalues)=0;
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 24433)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 24434)
@@ -432,5 +432,5 @@
 
 		/*Recover parameters and values*/
-		GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
+		Element::GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
 
 		/*Be sure that values are not zero*/
@@ -554,5 +554,5 @@
 
 		/*Recover parameters and values*/
-		GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
+		Element::GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
 
 		/*Be sure that values are not zero*/
@@ -985,5 +985,5 @@
 	/*Get current field and vertex coordinates*/
 	IssmDouble ls[NUMVERTICES],distance;
-	GetInputListOnVertices(&ls[0],distanceenum);
+	Element::GetInputListOnVertices(&ls[0],distanceenum);
 
 	/*Get distance from list of segments and reset ls*/
@@ -1110,9 +1110,9 @@
 	IssmDouble  sigmayz[NUMVERTICES],sigmaxz[NUMVERTICES],sigma_nn[NUMVERTICES];
 	IssmDouble  viscosity,epsilon[NUMVERTICES];
-	GetInputListOnVertices(&base[0],BaseEnum);
-	GetInputListOnVertices(&bed[0],BedEnum);
-	GetInputListOnVertices(&surface[0],SurfaceEnum);
-	GetInputListOnVertices(&pressure[0],PressureEnum);
-	GetInputListOnVertices(&phi[0],MaskGroundediceLevelsetEnum);
+	Element::GetInputListOnVertices(&base[0],BaseEnum);
+	Element::GetInputListOnVertices(&bed[0],BedEnum);
+	Element::GetInputListOnVertices(&surface[0],SurfaceEnum);
+	Element::GetInputListOnVertices(&pressure[0],PressureEnum);
+	Element::GetInputListOnVertices(&phi[0],MaskGroundediceLevelsetEnum);
 	IssmDouble rho_ice   = FindParam(MaterialsRhoIceEnum);
 	IssmDouble rho_water = FindParam(MaterialsRhoSeawaterEnum);
@@ -1253,5 +1253,5 @@
 
 	/*Recover parameters and values*/
-	GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum);
+	Element::GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum);
 
 	/*Be sure that values are not zero*/
@@ -1307,5 +1307,5 @@
 
 	/*Recover parameters and values*/
-	GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum);
+	Element::GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum);
 
 	/*Be sure that values are not zero*/
@@ -1413,8 +1413,8 @@
 
 	/*Retrieve all inputs and parameters*/
-	GetInputListOnVertices(&bed[0],BedEnum);
-	GetInputListOnVertices(&surfaces[0],SurfaceEnum);
-	GetInputListOnVertices(&bases[0],BaseEnum);
-	GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
+	Element::GetInputListOnVertices(&bed[0],BedEnum);
+	Element::GetInputListOnVertices(&surfaces[0],SurfaceEnum);
+	Element::GetInputListOnVertices(&bases[0],BaseEnum);
+	Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
 
 	nrfrontbed=0;
@@ -1508,5 +1508,5 @@
 
 	/*Recover parameters and values*/
-	GetInputListOnVertices(&levelset[0],levelsetenum);
+	Element::GetInputListOnVertices(&levelset[0],levelsetenum);
 
 	int* indicesfront = xNew<int>(NUMVERTICES);
@@ -1596,4 +1596,44 @@
 	return input;
 }/*}}}*/
+void       Penta::GetInputListOnVertices(IssmDouble* pvalue,Input2* input,IssmDouble default_value){/*{{{*/
+
+	/*Checks in debugging mode*/
+	_assert_(pvalue);
+
+	/* Start looping on the number of vertices: */
+	if(input){
+		GaussPenta gauss;
+		for(int iv=0;iv<NUMVERTICES;iv++){
+			gauss.GaussVertex(iv);
+			input->GetInputValue(&pvalue[iv],&gauss);
+		}
+	}
+	else{
+		for(int iv=0;iv<NUMVERTICES;iv++) pvalue[iv] = default_value;
+	}
+}
+/*}}}*/
+void       Penta::GetInputListOnNodes(IssmDouble* pvalue,Input2* input,IssmDouble default_value){/*{{{*/
+
+	/*Checks in debugging mode*/
+	_assert_(pvalue);
+
+	/*What type of finite element are we dealing with?*/
+	int fe       = this->FiniteElement();
+	int numnodes = this->GetNumberOfNodes();
+
+	/* Start looping on the number of vertices: */
+	if(input){
+		GaussPenta gauss;
+		for(int iv=0;iv<numnodes;iv++){
+			gauss.GaussNode(fe,iv);
+			input->GetInputValue(&pvalue[iv],&gauss);
+		}
+	}
+	else{
+		for(int iv=0;iv<numnodes;iv++) pvalue[iv] = default_value;
+	}
+}
+/*}}}*/
 DatasetInput2* Penta::GetDatasetInput2(int inputenum){/*{{{*/
 
@@ -1675,5 +1715,5 @@
 
 	/*Retrieve all inputs and parameters*/
-	GetInputListOnVertices(&lsf[0],levelset_enum);
+	Element::GetInputListOnVertices(&lsf[0],levelset_enum);
 
 	/* Determine distribution of ice over element.
@@ -1911,5 +1951,5 @@
 
 	/*Recover parameters and values*/
-	GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum);
+	Element::GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum);
 
 	/*Be sure that values are not zero*/
@@ -2030,5 +2070,5 @@
 
 	/*Recover parameters and values*/
-	GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
+	Element::GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
 
 	/*Be sure that values are not zero*/
@@ -2291,5 +2331,5 @@
 			/*Extrude values first*/
 			IssmDouble extrudedvalues[NUMVERTICES];
-			this->GetInputListOnVertices(&extrudedvalues[0],pentainput);
+			this->GetInputListOnVertices(&extrudedvalues[0],pentainput,0.);
 
 			if(start==-1){
@@ -2359,7 +2399,7 @@
 		IssmDouble extrudedvalues3[NUMVERTICES];
 
-		this->GetInputListOnVertices(&extrudedvalues[0],pentainput);
-		this->GetInputListOnVertices(&extrudedvalues2[0],pentainput2);
-		this->GetInputListOnVertices(&extrudedvalues3[0],pentainput3);
+		this->GetInputListOnVertices(&extrudedvalues[0],pentainput,0.);
+		this->GetInputListOnVertices(&extrudedvalues2[0],pentainput2,0.);
+		this->GetInputListOnVertices(&extrudedvalues3[0],pentainput3,0.);
 
 		if(start==-1){
@@ -2425,5 +2465,5 @@
 		IssmDouble extrudedvalues[NUMVERTICES];
 
-		this->GetInputListOnVertices(&extrudedvalues[0],enum_type);
+		Element::GetInputListOnVertices(&extrudedvalues[0],enum_type);
 		if(start==-1){
 			for(int i=0;i<NUMVERTICES2D;i++) extrudedvalues[i+NUMVERTICES2D]=extrudedvalues[i];
@@ -2628,5 +2668,5 @@
 
 	/*Retrieve all inputs and parameters*/
-	GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
+	Element::GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
 
 	/* If only one vertex has ice, there is an ice front here */
@@ -2660,5 +2700,5 @@
 
 	/*Retrieve all inputs and parameters*/
-	GetInputListOnVertices(&ls[0],levelset_enum);
+	Element::GetInputListOnVertices(&ls[0],levelset_enum);
 
 	/*If the level set has always same sign, there is no ice front here*/
@@ -3011,7 +3051,7 @@
 	rho_ice=FindParam(MaterialsRhoIceEnum);
 	density=rho_ice/rho_water;
-	GetInputListOnVertices(&h[0],ThicknessEnum);
-	GetInputListOnVertices(&r[0],BedEnum);
-	GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum);
+	Element::GetInputListOnVertices(&h[0],ThicknessEnum);
+	Element::GetInputListOnVertices(&r[0],BedEnum);
+	Element::GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum);
 
 	/*go through vertices, and figure out which ones are on the ice sheet, and want to unground: */
@@ -3187,5 +3227,5 @@
 	Input2* qsg_input = this->GetInput2(FrontalForcingsSubglacialDischargeEnum);		 _assert_(qsg_input);
 	Input2* TF_input  = this->GetInput2(FrontalForcingsThermalForcingEnum);          _assert_(TF_input);
-	GetInputListOnVertices(&basinid[0],FrontalForcingsBasinIdEnum);
+	Element::GetInputListOnVertices(&basinid[0],FrontalForcingsBasinIdEnum);
 
 	this->FindParam(&yts, ConstantsYtsEnum);
@@ -3750,5 +3790,5 @@
 
 	/*Recover parameters and values*/
-	GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
+	Element::GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
 
 	/*Be sure that values are not zero*/
@@ -3868,5 +3908,5 @@
 
 	/*Recover parameters and values*/
-	GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
+	Element::GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
 
 	/*Be sure that values are not zero*/
@@ -4735,8 +4775,8 @@
 
 					/*retrieve inputs: */
-					GetInputListOnVertices(&thickness_init[0],ThicknessEnum);
-					GetInputListOnVertices(&hydrostatic_ratio[0],GeometryHydrostaticRatioEnum);
-					GetInputListOnVertices(&bed[0],BaseEnum);
-					GetInputListOnVertices(&surface[0],SurfaceEnum);
+					Element::GetInputListOnVertices(&thickness_init[0],ThicknessEnum);
+					Element::GetInputListOnVertices(&hydrostatic_ratio[0],GeometryHydrostaticRatioEnum);
+					Element::GetInputListOnVertices(&bed[0],BaseEnum);
+					Element::GetInputListOnVertices(&surface[0],SurfaceEnum);
 
 					/*build new thickness: */
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 24433)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 24434)
@@ -85,4 +85,6 @@
 		Input2*        GetInput2(int enumtype);
 		Input2*        GetInput2(int enumtype,IssmDouble time);
+		void        GetInputListOnVertices(IssmDouble* pvalue,Input2* input,IssmDouble default_value);
+		void        GetInputListOnNodes(IssmDouble* pvalue,Input2* input,IssmDouble default_value);
 		DatasetInput2* GetDatasetInput2(int inputenum);
 		void           GetInputValue(IssmDouble* pvalue,Vertex* vertex,int enumtype);
Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp	(revision 24433)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.cpp	(revision 24434)
@@ -164,5 +164,6 @@
 
 	/*Recover parameters and values*/
-	GetInputListOnVertices(&levelset[0],levelsetenum);
+	Element::GetInputListOnVertices(&levelset[0],levelsetenum);
+
 	/* Get nodes where there is no ice */
 	nrfrontnodes=0;
@@ -263,5 +264,5 @@
 
 	/*Recover parameters and values*/
-	GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum);
+	Element::GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum);
 
 	/*Be sure that values are not zero*/
@@ -300,4 +301,44 @@
 
 }/*}}}*/
+void       Seg::GetInputListOnVertices(IssmDouble* pvalue,Input2* input,IssmDouble default_value){/*{{{*/
+
+	/*Checks in debugging mode*/
+	_assert_(pvalue);
+
+	/* Start looping on the number of vertices: */
+	if(input){
+		GaussSeg gauss;
+		for(int iv=0;iv<NUMVERTICES;iv++){
+			gauss.GaussVertex(iv);
+			input->GetInputValue(&pvalue[iv],&gauss);
+		}
+	}
+	else{
+		for(int iv=0;iv<NUMVERTICES;iv++) pvalue[iv] = default_value;
+	}
+}
+/*}}}*/
+void       Seg::GetInputListOnNodes(IssmDouble* pvalue,Input2* input,IssmDouble default_value){/*{{{*/
+
+	/*Checks in debugging mode*/
+	_assert_(pvalue);
+
+	/*What type of finite element are we dealing with?*/
+	int fe       = this->FiniteElement();
+	int numnodes = this->GetNumberOfNodes();
+
+	/* Start looping on the number of vertices: */
+	if(input){
+		GaussSeg gauss;
+		for(int iv=0;iv<numnodes;iv++){
+			gauss.GaussNode(fe,iv);
+			input->GetInputValue(&pvalue[iv],&gauss);
+		}
+	}
+	else{
+		for(int iv=0;iv<numnodes;iv++) pvalue[iv] = default_value;
+	}
+}
+/*}}}*/
 bool       Seg::IsIcefront(void){/*{{{*/
 
@@ -307,5 +348,5 @@
 
 	/*Retrieve all inputs and parameters*/
-	GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
+	Element::GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
 
 	/* If only one vertex has ice, there is an ice front here */
Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 24433)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 24434)
@@ -65,4 +65,6 @@
 		Input2*     GetInput2(int enumtype);
 		Input2*     GetInput2(int enumtype,IssmDouble time);
+		void        GetInputListOnVertices(IssmDouble* pvalue,Input2* input,IssmDouble default_value);
+		void        GetInputListOnNodes(IssmDouble* pvalue,Input2* input,IssmDouble default_value);
 		void        GetInputValue(IssmDouble* pvalue,Vertex* vertex,int enumtype){_error_("not implemented yet");};
 		void		   GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level){_error_("not implemented");};
Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp	(revision 24433)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp	(revision 24434)
@@ -194,5 +194,5 @@
 
 	/*Retrieve all inputs and parameters*/
-	GetInputListOnVertices(&values[0],MeshVertexonbaseEnum);
+	Element::GetInputListOnVertices(&values[0],MeshVertexonbaseEnum);
 
 	for(int i=0;i<4;i++){
@@ -214,5 +214,5 @@
 
 	/*Retrieve all inputs and parameters*/
-	GetInputListOnVertices(&values[0],MaskIceLevelsetEnum);
+	Element::GetInputListOnVertices(&values[0],MaskIceLevelsetEnum);
 
 	for(int i=0;i<4;i++){
@@ -234,5 +234,5 @@
 
 	/*Retrieve all inputs and parameters*/
-	GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
+	Element::GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
 
 	for(int i=0;i<4;i++){
@@ -263,4 +263,44 @@
 	_error_("not implemented yet");
 }/*}}}*/
+void       Tetra::GetInputListOnVertices(IssmDouble* pvalue,Input2* input,IssmDouble default_value){/*{{{*/
+
+	/*Checks in debugging mode*/
+	_assert_(pvalue);
+
+	/* Start looping on the number of vertices: */
+	if(input){
+		GaussTetra gauss;
+		for(int iv=0;iv<NUMVERTICES;iv++){
+			gauss.GaussVertex(iv);
+			input->GetInputValue(&pvalue[iv],&gauss);
+		}
+	}
+	else{
+		for(int iv=0;iv<NUMVERTICES;iv++) pvalue[iv] = default_value;
+	}
+}
+/*}}}*/
+void       Tetra::GetInputListOnNodes(IssmDouble* pvalue,Input2* input,IssmDouble default_value){/*{{{*/
+
+	/*Checks in debugging mode*/
+	_assert_(pvalue);
+
+	/*What type of finite element are we dealing with?*/
+	int fe       = this->FiniteElement();
+	int numnodes = this->GetNumberOfNodes();
+
+	/* Start looping on the number of vertices: */
+	if(input){
+		GaussTetra gauss;
+		for(int iv=0;iv<numnodes;iv++){
+			gauss.GaussNode(fe,iv);
+			input->GetInputValue(&pvalue[iv],&gauss);
+		}
+	}
+	else{
+		for(int iv=0;iv<numnodes;iv++) pvalue[iv] = default_value;
+	}
+}
+/*}}}*/
 void     Tetra::GetInputValue(IssmDouble* pvalue,Node* node,int enumtype){/*{{{*/
 
@@ -323,5 +363,5 @@
 
 	/*Retrieve all inputs and parameters*/
-	GetInputListOnVertices(&values[0],MeshVertexonbaseEnum);
+	Element::GetInputListOnVertices(&values[0],MeshVertexonbaseEnum);
 	sum = values[0]+values[1]+values[2]+values[3];
 
@@ -342,5 +382,5 @@
 
 	/*Retrieve all inputs and parameters*/
-	GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
+	Element::GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
 	sum = values[0]+values[1]+values[2]+values[3];
 
@@ -413,5 +453,5 @@
 	/*Retrieve all inputs and parameters*/
 	IssmDouble ls[NUMVERTICES];
-	GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
+	Element::GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
 
 	/* If only one vertex has ice, there is an ice front here */
Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 24433)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 24434)
@@ -69,4 +69,6 @@
 		Input2*     GetInput2(int enumtype);
 		Input2*     GetInput2(int enumtype,IssmDouble time);
+		void        GetInputListOnVertices(IssmDouble* pvalue,Input2* input,IssmDouble default_value);
+		void        GetInputListOnNodes(IssmDouble* pvalue,Input2* input,IssmDouble default_value);
 		void        GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);
 		void		   GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 24433)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 24434)
@@ -561,5 +561,5 @@
 		/*Recover parameters and values*/
 		parameters->FindParam(&domaintype,DomainTypeEnum);
-		GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
+		Element::GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
 
 		/*Be sure that values are not zero*/
@@ -698,5 +698,5 @@
 		/*Recover parameters and values*/
 		parameters->FindParam(&domaintype,DomainTypeEnum);
-		GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
+		Element::GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
 
 		/*Be sure that values are not zero*/
@@ -1169,5 +1169,5 @@
 	/*Get current field and vertex coordinates*/
 	IssmDouble ls[NUMVERTICES],distance;
-	GetInputListOnVertices(&ls[0],distanceenum);
+	Element::GetInputListOnVertices(&ls[0],distanceenum);
 
 	/*Get distance from list of segments and reset ls*/
@@ -1197,5 +1197,5 @@
 
 	/*Retrieve all inputs and parameters*/
-	GetInputListOnVertices(&values[0],MeshVertexonbaseEnum);
+	Element::GetInputListOnVertices(&values[0],MeshVertexonbaseEnum);
 
 	for(int i=0;i<3;i++){
@@ -1215,5 +1215,5 @@
 
 	/*Retrieve all inputs and parameters*/
-	GetInputListOnVertices(&values[0],MeshVertexonbaseEnum);
+	Element::GetInputListOnVertices(&values[0],MeshVertexonbaseEnum);
 
 	for(int i=0;i<3;i++){
@@ -1235,5 +1235,5 @@
 
 	/*Retrieve all inputs and parameters*/
-	GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
+	Element::GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
 
 	for(int i=0;i<3;i++){
@@ -1253,5 +1253,5 @@
 
 	/*Retrieve all inputs and parameters*/
-	GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
+	Element::GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
 
 	for(int i=0;i<3;i++){
@@ -1360,9 +1360,9 @@
 	IssmDouble  sigmaxx[NUMVERTICES],sigmayy[NUMVERTICES],sigmaxy[NUMVERTICES],sigma_nn[NUMVERTICES];
 	IssmDouble  viscosity,epsilon[NUMVERTICES];
-	GetInputListOnVertices(&base[0],BaseEnum);
-	GetInputListOnVertices(&bed[0],BedEnum);
-	GetInputListOnVertices(&surface[0],SurfaceEnum);
-	GetInputListOnVertices(&pressure[0],PressureEnum);
-	GetInputListOnVertices(&phi[0],MaskGroundediceLevelsetEnum);
+	Element::GetInputListOnVertices(&base[0],BaseEnum);
+	Element::GetInputListOnVertices(&bed[0],BedEnum);
+	Element::GetInputListOnVertices(&surface[0],SurfaceEnum);
+	Element::GetInputListOnVertices(&pressure[0],PressureEnum);
+	Element::GetInputListOnVertices(&phi[0],MaskGroundediceLevelsetEnum);
 	IssmDouble rho_ice   = FindParam(MaterialsRhoIceEnum);
 	IssmDouble rho_water = FindParam(MaterialsRhoSeawaterEnum);
@@ -1570,5 +1570,5 @@
 
 	/*Recover parameters and values*/
-	GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum);
+	Element::GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum);
 
 	/*Be sure that values are not zero*/
@@ -1626,5 +1626,5 @@
 	/*Recover parameters and values*/
 	parameters->FindParam(&domaintype,DomainTypeEnum);
-	GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum);
+	Element::GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum);
 
 	/*Be sure that values are not zero*/
@@ -1752,8 +1752,8 @@
 
 	/*Retrieve all inputs and parameters*/
-	GetInputListOnVertices(&bed[0],BedEnum);
-	GetInputListOnVertices(&surfaces[0],SurfaceEnum);
-	GetInputListOnVertices(&bases[0],BaseEnum);
-	GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
+	Element::GetInputListOnVertices(&bed[0],BedEnum);
+	Element::GetInputListOnVertices(&surfaces[0],SurfaceEnum);
+	Element::GetInputListOnVertices(&bases[0],BaseEnum);
+	Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);
 
 	nrfrontbed=0;
@@ -1846,5 +1846,5 @@
 
 	/*Recover parameters and values*/
-	GetInputListOnVertices(&levelset[0],levelsetenum);
+	Element::GetInputListOnVertices(&levelset[0],levelsetenum);
 
 	/* Get nodes where there is no ice */
@@ -1911,4 +1911,44 @@
 	}
 }/*}}}*/
+void       Tria::GetInputListOnVertices(IssmDouble* pvalue,Input2* input,IssmDouble default_value){/*{{{*/
+
+	/*Checks in debugging mode*/
+	_assert_(pvalue);
+
+	/* Start looping on the number of vertices: */
+	if(input){
+		GaussTria gauss;
+		for(int iv=0;iv<NUMVERTICES;iv++){
+			gauss.GaussVertex(iv);
+			input->GetInputValue(&pvalue[iv],&gauss);
+		}
+	}
+	else{
+		for(int iv=0;iv<NUMVERTICES;iv++) pvalue[iv] = default_value;
+	}
+}
+/*}}}*/
+void       Tria::GetInputListOnNodes(IssmDouble* pvalue,Input2* input,IssmDouble default_value){/*{{{*/
+
+	/*Checks in debugging mode*/
+	_assert_(pvalue);
+
+	/*What type of finite element are we dealing with?*/
+	int fe       = this->FiniteElement();
+	int numnodes = this->GetNumberOfNodes();
+
+	/* Start looping on the number of vertices: */
+	if(input){
+		GaussTria gauss;
+		for(int iv=0;iv<numnodes;iv++){
+			gauss.GaussNode(fe,iv);
+			input->GetInputValue(&pvalue[iv],&gauss);
+		}
+	}
+	else{
+		for(int iv=0;iv<numnodes;iv++) pvalue[iv] = default_value;
+	}
+}
+/*}}}*/
 void       Tria::InputServe(Input2* input_in){/*{{{*/
 
@@ -2252,5 +2292,5 @@
 
 	/*Recover parameters and values*/
-	GetInputListOnVertices(&levelset[0],levelsetenum);
+	Element::GetInputListOnVertices(&levelset[0],levelsetenum);
 
 	/* Get nodes where there is no ice */
@@ -2297,5 +2337,5 @@
 
 	/*Retrieve all inputs and parameters*/
-	GetInputListOnVertices(&lsf[0],levelset_enum);
+	Element::GetInputListOnVertices(&lsf[0],levelset_enum);
 
 	/* Determine distribution of ice over element.
@@ -2602,5 +2642,5 @@
 
 	/*Retrieve all inputs and parameters*/
-	GetInputListOnVertices(&values[0],MeshVertexonbaseEnum);
+	Element::GetInputListOnVertices(&values[0],MeshVertexonbaseEnum);
 	sum = values[0]+values[1]+values[2];
 
@@ -2623,5 +2663,5 @@
 
 	/*Retrieve all inputs and parameters*/
-	GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
+	Element::GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
 	sum = values[0]+values[1]+values[2];
 
@@ -2714,5 +2754,5 @@
 	/*Recover parameters and values*/
 	parameters->FindParam(&domaintype,DomainTypeEnum);
-	GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
+	Element::GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
 
 	/*Be sure that values are not zero*/
@@ -2843,5 +2883,5 @@
 	/*Recover parameters and values*/
 	parameters->FindParam(&domaintype,DomainTypeEnum);
-	GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum);
+	Element::GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum);
 
 	/*Be sure that values are not zero*/
@@ -2979,5 +3019,5 @@
 		area_base=this->GetAreaIce();
 		if(scaled==true){
-			GetInputListOnVertices(&scalefactors[0],MeshScaleFactorEnum);
+			Element::GetInputListOnVertices(&scalefactors[0],MeshScaleFactorEnum);
 			for(i=0;i<NUMVERTICES;i++) SFaux[i]= scalefactors[indices[i]]; //sort thicknesses in ice/noice
 			switch(numiceverts){
@@ -3002,6 +3042,6 @@
 			area_base=area_base*scalefactor;
 		}
-		GetInputListOnVertices(&surfaces[0],SurfaceEnum);
-		GetInputListOnVertices(&bases[0],BaseEnum);
+		Element::GetInputListOnVertices(&surfaces[0],SurfaceEnum);
+		Element::GetInputListOnVertices(&bases[0],BaseEnum);
 		for(i=0;i<NUMVERTICES;i++) Haux[i]= surfaces[indices[i]]-bases[indices[i]]; //sort thicknesses in ice/noice
 		switch(numiceverts){
@@ -3270,5 +3310,5 @@
 
 	/*Retrieve all inputs and parameters*/
-	GetInputListOnVertices(&values[0],MeshVertexonboundaryEnum);
+	Element::GetInputListOnVertices(&values[0],MeshVertexonboundaryEnum);
 	sum = values[0]+values[1]+values[2];
 
@@ -3291,5 +3331,5 @@
 
 	/*Retrieve all inputs and parameters*/
-	GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
+	Element::GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
 
 	/* If only one vertex has ice, there is an ice front here */
@@ -3323,5 +3363,5 @@
 
 	/*Retrieve all inputs and parameters*/
-	GetInputListOnVertices(&ls[0],levelset_enum);
+	Element::GetInputListOnVertices(&ls[0],levelset_enum);
 
 	/*If the level set is awlays <0, there is no ice front here*/
@@ -3773,7 +3813,7 @@
 	rho_ice=FindParam(MaterialsRhoIceEnum);
 	density=rho_ice/rho_water;
-	GetInputListOnVertices(&h[0],ThicknessEnum);
-	GetInputListOnVertices(&r[0],BedEnum);
-	GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum);
+	Element::GetInputListOnVertices(&h[0],ThicknessEnum);
+	Element::GetInputListOnVertices(&r[0],BedEnum);
+	Element::GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum);
 
 	/*go through vertices, and figure out which ones are grounded and want to unground: */
@@ -3915,5 +3955,5 @@
 	Input2* qsg_input = this->GetInput2(FrontalForcingsSubglacialDischargeEnum);		 _assert_(qsg_input);
 	Input2* TF_input  = this->GetInput2(FrontalForcingsThermalForcingEnum);          _assert_(TF_input);
-	GetInputListOnVertices(&basinid[0],FrontalForcingsBasinIdEnum);
+	Element::GetInputListOnVertices(&basinid[0],FrontalForcingsBasinIdEnum);
 
 	this->FindParam(&yts, ConstantsYtsEnum);
@@ -4326,5 +4366,5 @@
 	/*Recover parameters and values*/
 	parameters->FindParam(&domaintype,DomainTypeEnum);
-	GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
+	Element::GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
 
 	/*Be sure that values are not zero*/
@@ -4455,5 +4495,5 @@
 	/*Recover parameters and values*/
 	parameters->FindParam(&domaintype,DomainTypeEnum);
-	GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
+	Element::GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);
 
 	/*Be sure that values are not zero*/
@@ -4936,5 +4976,5 @@
 	/*Get field on vertices (we do not allow for higher order elements!!)*/
 	IssmDouble lsf[NUMVERTICES];
-	this->GetInputListOnVertices(&lsf[0],fieldenum);
+	Element::GetInputListOnVertices(&lsf[0],fieldenum);
 
 	/*1. check that we do cross fieldvalue in this element*/
@@ -6144,8 +6184,8 @@
 
 					/*retrieve inputs: */
-					GetInputListOnVertices(&thickness_init[0],ThicknessEnum);
-					GetInputListOnVertices(&hydrostatic_ratio[0],GeometryHydrostaticRatioEnum);
-					GetInputListOnVertices(&bed[0],BaseEnum);
-					GetInputListOnVertices(&surface[0],SurfaceEnum);
+					Element::GetInputListOnVertices(&thickness_init[0],ThicknessEnum);
+					Element::GetInputListOnVertices(&hydrostatic_ratio[0],GeometryHydrostaticRatioEnum);
+					Element::GetInputListOnVertices(&bed[0],BaseEnum);
+					Element::GetInputListOnVertices(&surface[0],SurfaceEnum);
 
 					/*build new bed and surface: */
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 24433)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 24434)
@@ -84,4 +84,6 @@
 		IssmDouble  GetIcefrontArea();
 		void	      GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum);
+		void        GetInputListOnVertices(IssmDouble* pvalue,Input2* input,IssmDouble default_value);
+		void        GetInputListOnNodes(IssmDouble* pvalue,Input2* input,IssmDouble default_value);
 		void	      GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level);
 		int         GetVertexIndex(Vertex* vertex);