Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5784)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5785)
@@ -3577,8 +3577,13 @@
 	double du[3];
 	double dv[3];
-	double dudx,dvdy;
+	double dw[3];
+	double dudx,dvdy,dwdz;
+
+	/*Get some parameters*/
+	int approximation;
 
 	/*If on water, skip: */
 	if(IsOnWater())return;
+	inputs->GetParameterValue(&approximation,ApproximationEnum);
 
 	/*If we are on the bedrock, spawn a tria on the bedrock, and use it to build the 
@@ -3598,4 +3603,8 @@
 	Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input);
 	Input* vy_input=inputs->GetInput(VyEnum); ISSMASSERT(vy_input);
+	Input* vz_input=NULL;
+	if(approximation==PattynStokesApproximationEnum){
+		vz_input=inputs->GetInput(VzEnum); ISSMASSERT(vz_input);
+	}
 
 	/* Start  looping on the number of gaussian points: */
@@ -3608,4 +3617,9 @@
 		vx_input->GetParameterDerivativeValue(&du[0],&xyz_list[0][0],gauss);
 		vy_input->GetParameterDerivativeValue(&dv[0],&xyz_list[0][0],gauss);
+		if(approximation==PattynStokesApproximationEnum){
+			vz_input->GetParameterDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
+			dwdz=dw[2];
+		}
+		else dwdz=0;
 		dudx=du[0];
 		dvdy=dv[1];
@@ -3619,5 +3633,5 @@
 		/*Build pe_g_gaussian vector: */
 		for (i=0;i<NUMVERTICES;i++){
-			pe_g_gaussian[i]=(dudx+dvdy)*Jdet*gauss->weight*l1l6[i];
+			pe_g_gaussian[i]=(dudx+dvdy+dwdz)*Jdet*gauss->weight*l1l6[i];
 		}
 
Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 5784)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 5785)
@@ -4137,8 +4137,9 @@
 
 	/*input parameters for structural analysis (diagnostic): */
-	double  vx,vy;
+	double  vx,vy,vz;
 	double  meltingvalue;
 	double  slope[2];
 	double  dbdx,dbdy;
+	int     approximation;
 
 	/* Get node coordinates and dof list: */
@@ -4147,7 +4148,12 @@
 
 	/*Retrieve all inputs we will be needing: */
+	inputs->GetParameterValue(&approximation,ApproximationEnum);
 	Input* melting_input=inputs->GetInput(MeltingRateEnum); ISSMASSERT(melting_input);
 	Input* vx_input=inputs->GetInput(VxEnum);               ISSMASSERT(vx_input);
 	Input* vy_input=inputs->GetInput(VyEnum);               ISSMASSERT(vy_input);
+	Input* vz_input=NULL;
+	if(approximation==PattynStokesApproximationEnum){
+		vz_input=inputs->GetInput(VzEnum);            ISSMASSERT(vz_input);
+	}
 	Input* bed_input=inputs->GetInput(BedEnum);             ISSMASSERT(bed_input);
 
@@ -4165,4 +4171,8 @@
 		vx_input->GetParameterValue(&vx, gauss);
 		vy_input->GetParameterValue(&vy, gauss);
+		if(approximation==PattynStokesApproximationEnum){
+			vz_input->GetParameterValue(&vz, gauss);
+		}
+		else vz=0;
 
 		/*Get bed slope: */
@@ -4179,5 +4189,5 @@
 		/*Build gaussian vector: */
 		for(i=0;i<NUMVERTICES;i++){
-			pe_g_gaussian[i]=-Jdet*gauss->weight*(vx*dbdx+vy*dbdy-meltingvalue)*L[i];
+			pe_g_gaussian[i]=-Jdet*gauss->weight*(vx*dbdx+vy*dbdy-vz-meltingvalue)*L[i];
 		}