Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5933)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5934)
@@ -2371,32 +2371,20 @@
 ElementMatrix* Penta::CreateKMatrixDiagnosticPattynViscous(void){
 
-	/* local declarations */
-	const int    numdof=2*NUMVERTICES;
-	int             i,j,ig;
-	double       xyz_list[NUMVERTICES][3];
+	/*Constants*/
+	const int    numdof=NDOF2*NUMVERTICES;
+
+	/*Intermediaries */
+	int        i,j,ig;
+	double     xyz_list[NUMVERTICES][3];
+	double     Jdet;
+	double     viscosity,oldviscosity,newviscosity,viscosity_overshoot; //viscosity
+	double     epsilon[5],oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
+	double     D_scalar;
+	double     D[5][5]={0.0};            // material matrix, simple scalar matrix.
+	double     B[5][numdof];
+	double     Bprime[5][numdof];
+	double     Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
+	Tria*      tria=NULL;
 	GaussPenta *gauss=NULL;
-	double viscosity; //viscosity
-	double oldviscosity; //viscosity
-	double newviscosity; //viscosity
-	double viscosity_overshoot;
-	double epsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
-	double oldepsilon[5]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
-	double B[5][numdof];
-	double Bprime[5][numdof];
-	double L[2][numdof];
-	double D[5][5]={0.0};            // material matrix, simple scalar matrix.
-	double D_scalar;
-	double DL[2][2]={0.0}; //for basal drag
-	double DL_scalar;
-	double Ke_gg[numdof][numdof]={0.0}; //local element stiffness matrix 
-	double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
-	double Jdet;
-	double  slope[2]={0.0};
-	double  slope_magnitude;
-	double  alpha2_list[3];
-	double  alpha2;
-	double MAXSLOPE=.06; // 6 %
-	double MOUNTAINKEXPONENT=10;
-	Tria*  tria=NULL;
 
 	/*Initialize Element matrix and return if necessary*/
@@ -2418,13 +2406,12 @@
 		gauss->GaussPoint(ig);
 
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
+		GetBPattyn(&B[0][0], &xyz_list[0][0], gauss);
+		GetBprimePattyn(&Bprime[0][0], &xyz_list[0][0], gauss);
+
 		this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
 		this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
 		matice->GetViscosity3d(&viscosity, &epsilon[0]);
 		matice->GetViscosity3d(&oldviscosity, &oldepsilon[0]);
-
-		GetBPattyn(&B[0][0], &xyz_list[0][0], gauss);
-		GetBprimePattyn(&Bprime[0][0], &xyz_list[0][0], gauss);
-
-		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
 
 		newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
@@ -2495,27 +2482,16 @@
 ElementMatrix* Penta::CreateKMatrixDiagnosticStokesViscous(void){
 
-	const int numdof=NUMVERTICES*NDOF4;
-	int i,j;
-	int     ig;
+	/*Intermediaries */
+	int        i,j,ig,approximation;
+	double     Jdet,viscosity,stokesreconditioning;
 	double     xyz_list[NUMVERTICES][3];
-	double	  xyz_list_tria[NUMVERTICES2D][3];
-	double	  bed_normal[3];
-	double     Ke_temp[27][27]={0.0}; //for the six nodes and the bubble 
-	double     Ke_reduced[numdof][numdof]; //for the six nodes only
-	double     Ke_gaussian[27][27];
+	double     epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
 	double     B[8][27];
 	double     B_prime[8][27];
-	double     Jdet;
+	double     D_scalar;
 	double     D[8][8]={0.0};
-	double     D_scalar;
-	double     DLStokes[14][14]={0.0};
+	double     Ke_temp[27][27]={0.0}; //for the six nodes and the bubble 
+	double     Ke_gaussian[27][27];
 	GaussPenta *gauss=NULL;
-	double  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
-	double  viscosity;
-	double  alpha2_gauss;
-	Friction* friction=NULL;
-	double stokesreconditioning;
-	int analysis_type;
-	int approximation;
 
 	/*If on water or not Stokes, skip stiffness: */
@@ -2526,5 +2502,4 @@
 	/*Retrieve all inputs and parameters*/
 	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
-	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
 	parameters->FindParam(&stokesreconditioning,StokesReconditioningEnum);
 	Input* vx_input=inputs->GetInput(VxEnum); ISSMASSERT(vx_input);
@@ -2538,14 +2513,11 @@
 		gauss->GaussPoint(ig);
 
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
+		GetBStokes(&B[0][0],&xyz_list[0][0],gauss); 
+		GetBprimeStokes(&B_prime[0][0],&xyz_list[0][0],gauss); 
+
 		this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
 		matice->GetViscosity3dStokes(&viscosity,&epsilon[0]);
 
-		GetBStokes(&B[0][0],&xyz_list[0][0],gauss); 
-		GetBprimeStokes(&B_prime[0][0],&xyz_list[0][0],gauss); 
-
-		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
-
-		/* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant 
-		 * onto this scalar matrix, so that we win some computational time: */
 		D_scalar=gauss->weight*Jdet;
 		for (i=0;i<6;i++) D[i][i]=D_scalar*2*viscosity;
@@ -2571,32 +2543,24 @@
 ElementMatrix* Penta::CreateKMatrixDiagnosticStokesFriction(void){
 
+	/*Constants*/
 	const int numdof=NUMVERTICES*NDOF4;
 	const int numdof2d=NUMVERTICES2D*NDOF4;
-	int i,j;
-	int     ig;
+
+	/*Intermediaries */
+	int        i,j,ig;
+	int        analysis_type,approximation;
+	double     stokesreconditioning;
+	double     viscosity,alpha2_gauss,Jdet2d;
+	double	  bed_normal[3];
+	double     epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
 	double     xyz_list[NUMVERTICES][3];
 	double	  xyz_list_tria[NUMVERTICES2D][3];
-	double	  bed_normal[3];
-	double     Ke_temp[27][27]={0.0}; //for the six nodes and the bubble 
-	double     Ke_reduced[numdof][numdof]; //for the six nodes only
-	double     Ke_gaussian[27][27];
-	double     Ke_drag_gaussian[numdof2d][numdof2d];
-	double     B[8][27];
-	double     B_prime[8][27];
 	double     LStokes[14][numdof2d];
 	double     LprimeStokes[14][numdof2d];
-	double     Jdet;
-	double     Jdet2d;
-	double     D[8][8]={0.0};
-	double     D_scalar;
 	double     DLStokes[14][14]={0.0};
+	double     Ke_temp[27][27]={0.0}; //for the six nodes and the bubble 
+	double     Ke_drag_gaussian[numdof2d][numdof2d];
+	Friction*  friction=NULL;
 	GaussPenta *gauss=NULL;
-	double  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
-	double  viscosity;
-	double  alpha2_gauss;
-	Friction* friction=NULL;
-	double stokesreconditioning;
-	int analysis_type;
-	int approximation;
 
 	/*If on water or not Stokes, skip stiffness: */
@@ -2623,4 +2587,5 @@
 		gauss->GaussPoint(ig);
 
+		GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
 		GetLStokes(&LStokes[0][0], gauss);
 		GetLprimeStokes(&LprimeStokes[0][0], &xyz_list[0][0], gauss);
@@ -2630,5 +2595,4 @@
 
 		BedNormal(&bed_normal[0],xyz_list_tria);
-		GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
 		friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum);
 
Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 5933)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 5934)
@@ -2401,5 +2401,4 @@
 ElementMatrix* Tria::CreateKMatrixBalancedthickness_CG(void){
 
-
 	/*Constants*/
 	const int    numdof=NDOF1*NUMVERTICES;