Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 6076)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 6077)
@@ -2082,7 +2082,7 @@
 
 	/*Constants*/
-	const int    numdofp=2*NUMVERTICES;
-	const int    numdofs=4*NUMVERTICES;
-	const int    numdoftotal=(4+2)*NUMVERTICES;
+	const int    numdofp=NDOF2*NUMVERTICES;
+	const int    numdofs=NDOF4*NUMVERTICES;
+	const int    numdoftotal=(NDOF2+NDOF4)*NUMVERTICES;
 	int          i,j;
 
@@ -2359,4 +2359,5 @@
 	/*Intermediaries */
 	int        i,j,ig;
+	int        approximation;
 	double     xyz_list[NUMVERTICES][3];
 	double     Jdet;
@@ -2376,4 +2377,5 @@
 
 	/*Retrieve all inputs and parameters*/
+	inputs->GetParameterValue(&approximation,ApproximationEnum);
 	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
 	this->parameters->FindParam(&viscosity_overshoot,ViscosityOvershootEnum);
@@ -3045,7 +3047,7 @@
 
 		for(i=0;i<NUMVERTICES;i++){
-			pe->values[i*NDOF4+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dh1dh6[2][i]/2;
-			pe->values[i*NDOF4+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dh1dh6[2][i]/2;
-			pe->values[i*NDOF4+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dh1dh6[0][i]+dw[1]*dh1dh6[1][i]+dw[2]*dh1dh6[2][i])/2;
+			pe->values[i*NDOF4+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dh1dh6[2][i];
+			pe->values[i*NDOF4+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dh1dh6[2][i];
+			pe->values[i*NDOF4+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dh1dh6[0][i]+dw[1]*dh1dh6[1][i]+2*dw[2]*dh1dh6[2][i]);
 			pe->values[i*NDOF4+3]+=Jdet*gauss->weight*stokesreconditioning*dw[2]*l1l6[i];
 		}
@@ -3117,7 +3119,7 @@
 
 		for(i=0;i<NUMVERTICES2D;i++){
-			pe->values[i*NDOF4+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+viscosity*dw[2]*bed_normal[0])*l1l6[i];
-			pe->values[i*NDOF4+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+viscosity*dw[2]*bed_normal[1])*l1l6[i];
-			pe->values[i*NDOF4+2]+=Jdet2d*gauss->weight*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*l1l6[i];
+			pe->values[i*NDOF4+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[0])*l1l6[i];
+			pe->values[i*NDOF4+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[1])*l1l6[i];
+			pe->values[i*NDOF4+2]+=Jdet2d*gauss->weight*2*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*l1l6[i];
 		}
 	}