Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17924)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 17925)
@@ -2335,7 +2335,7 @@
 		element->NodalFunctions(basis, gauss);
 
-		x_coord=element->GetXcoord(gauss);
-		y_coord=element->GetYcoord(gauss);
-		if(dim==3) z_coord=element->GetZcoord(gauss);
+		x_coord=element->GetXcoord(xyz_list,gauss);
+		y_coord=element->GetYcoord(xyz_list,gauss);
+		if(dim==3) z_coord=element->GetZcoord(xyz_list,gauss);
 		else z_coord=0.;
 
@@ -2472,6 +2472,6 @@
 		gauss->GaussPoint(ig);
 		surface_input->GetInputValue(&surface,gauss);
-		if(dim==3) z=element->GetZcoord(gauss);
-		else       z=element->GetYcoord(gauss);
+		if(dim==3) z=element->GetZcoord(xyz_list,gauss);
+		else       z=element->GetYcoord(xyz_list,gauss);
 		element->NodalFunctions(basis,gauss);
 		element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss);
@@ -3092,7 +3092,7 @@
 		gauss->GaussPoint(ig);
 
-		x_coord=element->GetXcoord(gauss);
-		y_coord=element->GetYcoord(gauss);
-		if(dim==3) z_coord=element->GetZcoord(gauss);
+		x_coord=element->GetXcoord(xyz_list,gauss);
+		y_coord=element->GetYcoord(xyz_list,gauss);
+		if(dim==3) z_coord=element->GetZcoord(xyz_list,gauss);
 		else z_coord=0.;
 
@@ -3153,7 +3153,7 @@
 		element->NodalFunctionsVelocity(vbasis,gauss);
 
-		x_coord=element->GetXcoord(gauss);
-		y_coord=element->GetYcoord(gauss);
-		if(dim==3) z_coord=element->GetZcoord(gauss);
+		x_coord=element->GetXcoord(xyz_list,gauss);
+		y_coord=element->GetYcoord(xyz_list,gauss);
+		if(dim==3) z_coord=element->GetZcoord(xyz_list,gauss);
 		else z_coord=0.;
 
@@ -3670,6 +3670,6 @@
 		element->NodalFunctionsVelocity(vbasis,gauss);
 		surface_input->GetInputValue(&surface,gauss);
-		if(dim==3) z=element->GetZcoord(gauss);
-		else       z=element->GetYcoord(gauss);
+		if(dim==3) z=element->GetZcoord(xyz_list,gauss);
+		else       z=element->GetYcoord(xyz_list,gauss);
 		pressure = rho_water*gravity*min(0.,z);//0 if the gaussian point is above water level
 
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp	(revision 17924)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceSIAAnalysis.cpp	(revision 17925)
@@ -405,5 +405,5 @@
 			constant_part=-2.*pow(rho_ice*gravity,n)*pow(slope2,((n-1.)/2.));
 
-			z = element->GetZcoord(gauss);
+			z = element->GetZcoord(xyz_list,gauss);
 			element->JacobianDeterminantLine(&Jdet,&xyz_list_line[0][0],gauss);
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 17924)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 17925)
@@ -574,4 +574,46 @@
 }
 /*}}}*/
+IssmDouble Element::GetXcoord(IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+
+	/*output*/
+	IssmDouble x;
+
+	/*Create list of x*/
+	int         numvertices = this->GetNumberOfVertices();
+	IssmDouble* x_list      = xNew<IssmDouble>(numvertices);
+
+	for(int i=0;i<numvertices;i++) x_list[i]=xyz_list[i*3+0];
+	ValueP1OnGauss(&x,x_list,gauss);
+
+	return x;
+}/*}}}*/
+IssmDouble Element::GetYcoord(IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+
+	/*output*/
+	IssmDouble y;
+
+	/*Create list of y*/
+	int         numvertices = this->GetNumberOfVertices();
+	IssmDouble* y_list      = xNew<IssmDouble>(numvertices);
+
+	for(int i=0;i<numvertices;i++) y_list[i]=xyz_list[i*3+1];
+	ValueP1OnGauss(&y,y_list,gauss);
+
+	return y;
+}/*}}}*/
+IssmDouble Element::GetZcoord(IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+
+	/*output*/
+	IssmDouble z;
+
+	/*Create list of z*/
+	int         numvertices = this->GetNumberOfVertices();
+	IssmDouble* z_list      = xNew<IssmDouble>(numvertices);
+
+	for(int i=0;i<numvertices;i++) z_list[i]=xyz_list[i*3+2];
+	ValueP1OnGauss(&z,z_list,gauss);
+
+	return z;
+}/*}}}*/
 bool Element::HasNodeOnBase(){/*{{{*/
 	return (this->inputs->Max(MeshVertexonbaseEnum)>0.);
@@ -1100,5 +1142,5 @@
 	surface_input->GetInputValue(&s,gauss);
 	surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
-	z=GetZcoord(gauss);
+	z=this->GetZcoord(xyz_list,gauss);
 	tau_perp = matpar->GetRhoIce() * matpar->GetG() * fabs(s-z)*sqrt(slope[0]*slope[0]+slope[1]*slope[1]);
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 17924)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 17925)
@@ -91,4 +91,7 @@
 		void       GetVerticesSidList(int* sidlist);
 		void       GetVerticesConnectivityList(int* connectivitylist);
+		IssmDouble GetXcoord(IssmDouble* xyz_list,Gauss* gauss);
+		IssmDouble GetYcoord(IssmDouble* xyz_list,Gauss* gauss);
+		IssmDouble GetZcoord(IssmDouble* xyz_list,Gauss* gauss);
 		bool       HasNodeOnBase();
 		bool       HasNodeOnSurface();
@@ -194,7 +197,5 @@
 		virtual void   GetVerticesCoordinatesBase(IssmDouble** xyz_list)=0;
 		virtual void   GetVerticesCoordinatesTop(IssmDouble** xyz_list)=0;
-		virtual IssmDouble GetXcoord(Gauss* gauss)=0;
-		virtual IssmDouble GetYcoord(Gauss* gauss)=0;
-		virtual IssmDouble GetZcoord(Gauss* gauss)=0;
+
 		virtual int    GetElementType(void)=0;
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 17924)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 17925)
@@ -1054,49 +1054,4 @@
 
 	return tau_parameter;
-}
-/*}}}*/
-/*FUNCTION Penta::GetXcoord {{{*/
-IssmDouble Penta::GetXcoord(Gauss* gauss){
-
-	int    i;
-	IssmDouble x;
-	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble x_list[NUMVERTICES];
-
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	for(i=0;i<NUMVERTICES;i++) x_list[i]=xyz_list[i][0];
-	PentaRef::GetInputValue(&x,x_list,gauss,P1Enum);
-
-	return x;
-}
-/*}}}*/
-/*FUNCTION Penta::GetYcoord {{{*/
-IssmDouble Penta::GetYcoord(Gauss* gauss){
-
-	int    i;
-	IssmDouble y;
-	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble y_list[NUMVERTICES];
-
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	for(i=0;i<NUMVERTICES;i++) y_list[i]=xyz_list[i][1];
-	PentaRef::GetInputValue(&y,y_list,gauss,P1Enum);
-
-	return y;
-}
-/*}}}*/
-/*FUNCTION Penta::GetZcoord {{{*/
-IssmDouble Penta::GetZcoord(Gauss* gauss){
-
-	int    i;
-	IssmDouble z;
-	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble z_list[NUMVERTICES];
-
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	for(i=0;i<NUMVERTICES;i++) z_list[i]=xyz_list[i][2];
-	PentaRef::GetInputValue(&z,z_list,gauss,P1Enum);
-
-	return z;
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 17924)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 17925)
@@ -80,7 +80,4 @@
 		int    GetNumberOfVertices(void);
 		void   GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int enum_type);
-		IssmDouble GetXcoord(Gauss* gauss);
-		IssmDouble GetYcoord(Gauss* gauss);
-		IssmDouble GetZcoord(Gauss* gauss);
 		void   GetVerticesCoordinatesBase(IssmDouble** pxyz_list);
 		void   GetVerticesCoordinatesTop(IssmDouble** pxyz_list);
Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 17924)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 17925)
@@ -110,7 +110,4 @@
 		void        GetInputValue(IssmDouble* pvalue,Node* node,int enumtype){_error_("not implemented yet");};
 		Node*       GetNode(int node_number){_error_("Not implemented");};
-		IssmDouble  GetXcoord(Gauss* gauss){_error_("Not implemented");};
-		IssmDouble  GetYcoord(Gauss* gauss){_error_("Not implemented");};
-		IssmDouble  GetZcoord(Gauss* gauss){_error_("not implemented yet");};
 		int         GetElementType(void){_error_("not implemented yet");};
 		Gauss*      NewGauss(void);
Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp	(revision 17924)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp	(revision 17925)
@@ -272,49 +272,4 @@
 
 }/*}}}*/
-/*FUNCTION Tetra::GetXcoord {{{*/
-IssmDouble Tetra::GetXcoord(Gauss* gauss){
-
-	int    i;
-	IssmDouble z;
-	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble z_list[NUMVERTICES];
-
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	for(i=0;i<NUMVERTICES;i++) z_list[i]=xyz_list[i][0];
-	TetraRef::GetInputValue(&z,&z_list[0],(GaussTetra*)gauss,P1Enum);
-
-	return z;
-}
-/*}}}*/
-/*FUNCTION Tetra::GetYcoord {{{*/
-IssmDouble Tetra::GetYcoord(Gauss* gauss){
-
-	int    i;
-	IssmDouble z;
-	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble z_list[NUMVERTICES];
-
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	for(i=0;i<NUMVERTICES;i++) z_list[i]=xyz_list[i][1];
-	TetraRef::GetInputValue(&z,&z_list[0],(GaussTetra*)gauss,P1Enum);
-
-	return z;
-}
-/*}}}*/
-/*FUNCTION Tetra::GetZcoord {{{*/
-IssmDouble Tetra::GetZcoord(Gauss* gauss){
-
-	int    i;
-	IssmDouble z;
-	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble z_list[NUMVERTICES];
-
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	for(i=0;i<NUMVERTICES;i++) z_list[i]=xyz_list[i][2];
-	TetraRef::GetInputValue(&z,&z_list[0],(GaussTetra*)gauss,P1Enum);
-
-	return z;
-}
-/*}}}*/
 /*FUNCTION Tetra::HasFaceOnBase{{{*/
 bool Tetra::HasFaceOnBase(){
Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 17924)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h	(revision 17925)
@@ -115,7 +115,4 @@
 		void        GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);
 		Node*       GetNode(int node_number){_error_("Not implemented");};
-		IssmDouble  GetXcoord(Gauss* gauss);
-		IssmDouble  GetYcoord(Gauss* gauss);
-		IssmDouble  GetZcoord(Gauss* gauss);
 		int         GetElementType(void);
 		Gauss*      NewGauss(void);
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 17924)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 17925)
@@ -986,32 +986,4 @@
 
 }/*}}}*/
-/*FUNCTION Tria::GetXcoord {{{*/
-IssmDouble Tria::GetXcoord(Gauss* gauss){
-
-	IssmDouble x;
-	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble x_list[NUMVERTICES];
-
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	for(int i=0;i<NUMVERTICES;i++) x_list[i]=xyz_list[i][0];
-	TriaRef::GetInputValue(&x,x_list,gauss,P1Enum);
-
-	return x;
-}
-/*}}}*/
-/*FUNCTION Tria::GetYcoord {{{*/
-IssmDouble Tria::GetYcoord(Gauss* gauss){
-
-	IssmDouble y;
-	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble y_list[NUMVERTICES];
-
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	for(int i=0;i<NUMVERTICES;i++) y_list[i]=xyz_list[i][1];
-	TriaRef::GetInputValue(&y,y_list,gauss,P1Enum);
-
-	return y;
-}
-/*}}}*/
 /*FUNCTION Tria::InputDepthAverageAtBase {{{*/
 void  Tria::InputDepthAverageAtBase(int enum_type,int average_enum_type){
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 17924)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 17925)
@@ -184,7 +184,4 @@
 		void           GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints);
 		int            GetElementType(void);
-		IssmDouble     GetXcoord(Gauss* gauss);
-		IssmDouble     GetYcoord(Gauss* gauss);
-		IssmDouble     GetZcoord(Gauss* gauss){_error_("not implemented");};
 		void           NormalSection(IssmDouble* normal,IssmDouble* xyz_list);
 		void           NormalTop(IssmDouble* normal,IssmDouble* xyz_list);