Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15474)
+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15475)
@@ -6917,7 +6917,7 @@
 	GaussPenta *gauss=NULL;
 
 	/*Stabilization*/
-	bool       stabilization = false;
+	bool       stabilization = true;
 	IssmDouble dbasis[3][6];
 	IssmDouble dmu[3];
 	IssmDouble mu;
@@ -6943,6 +6943,20 @@
 	rigidity=material->GetB();
 	diameter=MinEdgeLength(xyz_list);
 
+	if(stabilization){
+		gauss=new GaussPenta();
+		for(int iv=0;iv<6;iv++){
+			gauss->GaussVertex(iv);
+			GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
+			for(i=0;i<6;i++){
+				for(j=0;j<3;j++){
+					dnodalbasis[i][iv][j] = dbasis[j][i];
+				}
+			}
+		}
+		delete gauss;
+	}
+
 	/* Start  looping on the number of gaussian points: */
 	gauss=new GaussPenta(3,2);
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
@@ -6971,14 +6985,7 @@
 		if(stabilization){
 			GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
 			dmu[0]=0.; dmu[1]=0.; dmu[2]=0.;
-			mu = viscosity;
-			for(i=0;i<6;i++){
-				for(p=0;p<6;p++){
-					for(j=0;j<3;j++){
-						dnodalbasis[i][p][j] = dbasis[j][i];
-					}
-				}
-			}
+			mu = 2.*viscosity;
 			for(p=0;p<6;p++) for(i=0;i<4;i++) for(j=0;j<4;j++){
 				SU[p][i][j]=0.;
 				SW[p][i][j]=0.;
@@ -7017,34 +7024,8 @@
 			}
 
 		}
-
-
-		//viscosity = 1000*rigidity;
-		//D_scalar_stab=-gauss->weight*Jdet*1./3.*pow(diameter,2)/(8.*2.*viscosity)*stokesreconditioning;
-		//if(id==24){
-		//printf("tau %g\n",1./3.*pow(diameter,2)/(8.*2.*viscosity));
-		//}
-		//GetBConduct(&B_stab[0][0],&xyz_list[0][0],gauss); 
-
-		//D_stab[0][0]=D_scalar_stab; D_stab[0][1]=0;             D_stab[0][2]=0;
-		//D_stab[1][0]=0;             D_stab[1][1]=D_scalar_stab; D_stab[1][2]=0;
-		//D_stab[2][0]=0;             D_stab[2][1]=0;             D_stab[2][2]=D_scalar_stab;
-
-		//TripleMultiply(&B_stab[0][0],3,6,1,
-		//			&D_stab[0][0],3,3,0,
-		//			&B_stab[0][0],3,6,0,
-		//			&Ke_temp_stab[0][0],0);
-
-		//for(i=0;i<NUMVERTICES;i++) for(j=0;j<NUMVERTICES;j++) Ke->values[numdof*(i*4+3)+j*4+3]+=Ke_temp_stab[i][j];
-
 	}
 
-//	if(id==24){
-//		_error_("");
-//	}
-
-	//Ke->Echo();
-	//_error_("stop");
 	/*Transform Coordinate System*/
 	TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYZPEnum);
 	if(id==24){
@@ -7943,8 +7924,22 @@
 	/*Find minimal length*/
 	diameter=MinEdgeLength(xyz_list);
 
+	if(stabilization){
+		gauss=new GaussPenta();
+		for(int iv=0;iv<6;iv++){
+			gauss->GaussVertex(iv);
+			GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
+			for(i=0;i<6;i++){
+				for(j=0;j<3;j++){
+					dnodalbasis[i][iv][j] = dbasis[j][i];
+				}
+			}
+		}
+		delete gauss;
+	}
+
 	/* Start  looping on the number of gaussian points: */
-	gauss=new GaussPenta(2,3);
+	gauss=new GaussPenta(3,2);
 	for(int ig=gauss->begin();ig<gauss->end();ig++){
 
 		gauss->GaussPoint(ig);
@@ -7969,15 +7964,7 @@
 		if(stabilization){
 			GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
 			dmu[0]=0.; dmu[1]=0.; dmu[2]=0.;
-			mu = viscosity;
-			for(i=0;i<6;i++){
-				for(p=0;p<6;p++){
-					for(j=0;j<3;j++){
-						dnodalbasis[i][p][j] = dbasis[j][i];
-					}
-				}
-			}
-			//dNodalBasisdx(1:n,p,:) = dBasisdx(1:n,:)
+			mu = 2.*viscosity;
 			for(p=0;p<6;p++) for(i=0;i<4;i++) for(j=0;j<4;j++){
 				SW[p][i][j]=0.;
 			}
@@ -7988,10 +7975,8 @@
 						SW[p][i][i] += -dmu[j]*dbasis[j][p];
 						SW[p][j][i] += -dmu[j]*dbasis[i][p];
 						for(ii=0;ii<6;ii++){
-							//SW[p][i][i] += -mu*dnodalbasis[p][ii][j]*dbasis[j][ii];
-							//SW[p][j][i] += -mu*dnodalbasis[p][ii][i]*dbasis[j][ii];
-							SW[p][i][i] += -mu*dnodalbasis[p][ii][j];
-							SW[p][j][i] += -mu*dnodalbasis[p][ii][i];
+							SW[p][i][i] += -mu*dnodalbasis[p][ii][j]*dbasis[j][ii];
+							SW[p][j][i] += -mu*dnodalbasis[p][ii][i]*dbasis[j][ii];
 						}
 					}
 				}
@@ -8006,26 +7991,11 @@
 			}
 
 		}
-		///*Add stabilization*/
-		//D_scalar_stab=-gauss->weight*Jdet*1./3.*pow(diameter,2)/(8.*2.*viscosity)*stokesreconditioning;
-		//GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0],gauss);
-
-		//for(i=0;i<NUMVERTICES;i++){
-		//	Pe_gaussian[i*NDOF4+3]+=-rho_ice*gravity*D_scalar_stab*dh1dh6[2][i];
-		//	Pe_gaussian[i*NDOF4+3]+=forcex*D_scalar_stab*dh1dh6[0][i];
-		//	Pe_gaussian[i*NDOF4+3]+=forcey*D_scalar_stab*dh1dh6[1][i];
-		//	Pe_gaussian[i*NDOF4+3]+=forcez*D_scalar_stab*dh1dh6[2][i];
-		//}
 	}
 
 	/*Transform coordinate system*/
 	TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZPEnum);
 
-	if(id==24){
-		printarray(pe->values,24,1);
-		_error_("STOP");
-	}
-
 	/*Clean up and return*/
 	delete gauss;
 	return pe;