Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 8652)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 8653)
@@ -1773,5 +1773,5 @@
 	double    surface_normal[3];
 	double    Jdet2d,DL_scalar;
-	double    l1l6[NUMVERTICES];
+	double    basis[NUMVERTICES];
 	GaussPenta *gauss=NULL;
 	double Ke_g[numdof][numdof]; //stiffness matrix evaluated at the gaussian point.
@@ -1792,5 +1792,5 @@
 
 		GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
-		GetNodalFunctionsP1(&l1l6[0], gauss);
+		GetNodalFunctionsP1(&basis[0], gauss);
 
 		/**********************Do not forget the sign**********************************/
@@ -1798,7 +1798,7 @@
 		/******************************************************************************/
 
-		TripleMultiply( l1l6,1,numdof,1,
+		TripleMultiply( basis,1,numdof,1,
 					&DL_scalar,1,1,0,
-					l1l6,1,numdof,0,
+					basis,1,numdof,0,
 					&Ke_g[0][0],0);
 
@@ -1848,5 +1848,5 @@
 	double     Bprime_advec[3][numdof];
 	double     L[numdof];
-	double     dh1dh6[3][6];
+	double     dbasis[3][6];
 	double     D_scalar_conduct,D_scalar_advec;
 	double     D_scalar_trans,D_scalar_artdiff;
@@ -1970,5 +1970,5 @@
 		else if(artdiff==2){
 
-			GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss);
+			GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
 
 			tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,rho_ice,heatcapacity,thermalconductivity);
@@ -1977,5 +1977,5 @@
 				for(j=0;j<numdof;j++){
 					Ke_gaussian_artdiff[i][j]=tau_parameter*D_scalar_advec*
-					  ((u-um)*dh1dh6[0][i]+(v-vm)*dh1dh6[1][i]+(w-wm)*dh1dh6[2][i])*((u-um)*dh1dh6[0][j]+(v-vm)*dh1dh6[1][j]+(w-wm)*dh1dh6[2][j]);
+					  ((u-um)*dbasis[0][i]+(v-vm)*dbasis[1][i]+(w-wm)*dbasis[2][i])*((u-um)*dbasis[0][j]+(v-vm)*dbasis[1][j]+(w-wm)*dbasis[2][j]);
 				}
 			}
@@ -1983,5 +1983,5 @@
 				for(i=0;i<numdof;i++){
 					for(j=0;j<numdof;j++){
-						Ke_gaussian_artdiff[i][j]+=tau_parameter*D_scalar_trans*L[j]*((u-um)*dh1dh6[0][i]+(v-vm)*dh1dh6[1][i]+(w-wm)*dh1dh6[2][i]);
+						Ke_gaussian_artdiff[i][j]+=tau_parameter*D_scalar_trans*L[j]*((u-um)*dbasis[0][i]+(v-vm)*dbasis[1][i]+(w-wm)*dbasis[2][i]);
 					}
 				}
@@ -2013,5 +2013,5 @@
 	double    xyz_list[NUMVERTICES][3];
 	double	 xyz_list_tria[NUMVERTICES2D][3];
-	double    l1l6[NUMVERTICES];
+	double    basis[NUMVERTICES];
 	double    D_scalar;
 	double    Ke_gaussian[numdof][numdof]={0.0};
@@ -2039,12 +2039,12 @@
 		
 		GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
-		GetNodalFunctionsP1(&l1l6[0], gauss);
+		GetNodalFunctionsP1(&basis[0], gauss);
 				
 		D_scalar=gauss->weight*Jdet2d*rho_water*mixed_layer_capacity*thermal_exchange_velocity/(rho_ice*heatcapacity);
 		if(dt) D_scalar=dt*D_scalar;
 
-		TripleMultiply(&l1l6[0],numdof,1,0,
+		TripleMultiply(&basis[0],numdof,1,0,
 					&D_scalar,1,1,0,
-					&l1l6[0],1,numdof,0,
+					&basis[0],1,numdof,0,
 					&Ke_gaussian[0][0],0);
 
@@ -2138,5 +2138,5 @@
 	double     Bprime_advec[3][numdof];
 	double     L[numdof];
-	double     dh1dh6[3][6];
+	double     dbasis[3][6];
 	double     D_scalar_conduct,D_scalar_advec;
 	double     D_scalar_trans,D_scalar_artdiff;
@@ -2253,5 +2253,5 @@
 		else if(artdiff==2){
 
-			GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss);
+			GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
 
 			tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,rho_ice,heatcapacity,thermalconductivity);
@@ -2260,5 +2260,5 @@
 				for(j=0;j<numdof;j++){
 					Ke_gaussian_artdiff[i][j]=tau_parameter*D_scalar_advec*
-					  ((u-um)*dh1dh6[0][i]+(v-vm)*dh1dh6[1][i]+(w-wm)*dh1dh6[2][i])*((u-um)*dh1dh6[0][j]+(v-vm)*dh1dh6[1][j]+(w-wm)*dh1dh6[2][j]);
+					  ((u-um)*dbasis[0][i]+(v-vm)*dbasis[1][i]+(w-wm)*dbasis[2][i])*((u-um)*dbasis[0][j]+(v-vm)*dbasis[1][j]+(w-wm)*dbasis[2][j]);
 				}
 			}
@@ -2266,5 +2266,5 @@
 				for(i=0;i<numdof;i++){
 					for(j=0;j<numdof;j++){
-						Ke_gaussian_artdiff[i][j]+=tau_parameter*D_scalar_trans*L[j]*((u-um)*dh1dh6[0][i]+(v-vm)*dh1dh6[1][i]+(w-wm)*dh1dh6[2][i]);
+						Ke_gaussian_artdiff[i][j]+=tau_parameter*D_scalar_trans*L[j]*((u-um)*dbasis[0][i]+(v-vm)*dbasis[1][i]+(w-wm)*dbasis[2][i]);
 					}
 				}
@@ -2297,5 +2297,5 @@
 	double    xyz_list[NUMVERTICES][3];
 	double	 xyz_list_tria[NUMVERTICES2D][3];
-	double    l1l6[NUMVERTICES];
+	double    basis[NUMVERTICES];
 	double    D_scalar;
 	double    Ke_gaussian[numdof][numdof]={0.0};
@@ -2323,12 +2323,12 @@
 		
 		GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
-		GetNodalFunctionsP1(&l1l6[0], gauss);
+		GetNodalFunctionsP1(&basis[0], gauss);
 				
 		D_scalar=gauss->weight*Jdet2d*rho_water*mixed_layer_capacity*thermal_exchange_velocity/(heatcapacity*rho_ice);
 		if(dt) D_scalar=dt*D_scalar;
 
-		TripleMultiply(&l1l6[0],numdof,1,0,
+		TripleMultiply(&basis[0],numdof,1,0,
 					&D_scalar,1,1,0,
-					&l1l6[0],1,numdof,0,
+					&basis[0],1,numdof,0,
 					&Ke_gaussian[0][0],0);
 
@@ -2522,6 +2522,6 @@
 	double      dw[3];
 	double      xyz_list[NUMVERTICES][3];
-	double      l1l6[6]; //for the six nodes of the penta
-	double      dh1dh6[3][6]; //for the six nodes of the penta
+	double      basis[6]; //for the six nodes of the penta
+	double      dbasis[3][6]; //for the six nodes of the penta
 	GaussPenta *gauss=NULL;
 
@@ -2546,6 +2546,6 @@
 
 		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
-		GetNodalFunctionsP1(&l1l6[0], gauss);
-		GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss);
+		GetNodalFunctionsP1(&basis[0], gauss);
+		GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
 
 		vzmacayeal_input->GetParameterDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
@@ -2555,8 +2555,8 @@
 
 		for(i=0;i<NUMVERTICES;i++){
-			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];
+			pe->values[i*NDOF4+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dbasis[2][i];
+			pe->values[i*NDOF4+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dbasis[2][i];
+			pe->values[i*NDOF4+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dbasis[0][i]+dw[1]*dbasis[1][i]+2*dw[2]*dbasis[2][i]);
+			pe->values[i*NDOF4+3]+=Jdet*gauss->weight*stokesreconditioning*dw[2]*basis[i];
 		}
 	}
@@ -2584,5 +2584,5 @@
 	double	   xyz_list_tria[NUMVERTICES2D][3];
 	double      xyz_list[NUMVERTICES][3];
-	double      l1l6[6]; //for the six nodes of the penta
+	double      basis[6]; //for the six nodes of the penta
 	Tria*       tria=NULL;
 	Friction*   friction=NULL;
@@ -2616,5 +2616,5 @@
 
 		GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
-		GetNodalFunctionsP1(l1l6, gauss);
+		GetNodalFunctionsP1(basis, gauss);
 
 		vzmacayeal_input->GetParameterValue(&w, gauss);
@@ -2627,7 +2627,7 @@
 
 		for(i=0;i<NUMVERTICES2D;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];
+			pe->values[i*NDOF4+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[0])*basis[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])*basis[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])*basis[i];
 		}
 	}
@@ -2667,6 +2667,6 @@
 	double      dw[3];
 	double      xyz_list[NUMVERTICES][3];
-	double      l1l6[6]; //for the six nodes of the penta
-	double      dh1dh6[3][6]; //for the six nodes of the penta
+	double      basis[6]; //for the six nodes of the penta
+	double      dbasis[3][6]; //for the six nodes of the penta
 	GaussPenta *gauss=NULL;
 
@@ -2691,6 +2691,6 @@
 
 		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
-		GetNodalFunctionsP1(&l1l6[0], gauss);
-		GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss);
+		GetNodalFunctionsP1(&basis[0], gauss);
+		GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
 
 		vzpattyn_input->GetParameterDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
@@ -2700,8 +2700,8 @@
 
 		for(i=0;i<NUMVERTICES;i++){
-			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];
+			pe->values[i*NDOF4+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dbasis[2][i];
+			pe->values[i*NDOF4+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dbasis[2][i];
+			pe->values[i*NDOF4+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dbasis[0][i]+dw[1]*dbasis[1][i]+2*dw[2]*dbasis[2][i]);
+			pe->values[i*NDOF4+3]+=Jdet*gauss->weight*stokesreconditioning*dw[2]*basis[i];
 		}
 	}
@@ -2729,5 +2729,5 @@
 	double	   xyz_list_tria[NUMVERTICES2D][3];
 	double      xyz_list[NUMVERTICES][3];
-	double      l1l6[6]; //for the six nodes of the penta
+	double      basis[6]; //for the six nodes of the penta
 	Tria*       tria=NULL;
 	Friction*   friction=NULL;
@@ -2761,5 +2761,5 @@
 
 		GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
-		GetNodalFunctionsP1(l1l6, gauss);
+		GetNodalFunctionsP1(basis, gauss);
 
 		vzpattyn_input->GetParameterValue(&w, gauss);
@@ -2772,7 +2772,7 @@
 
 		for(i=0;i<NUMVERTICES2D;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];
+			pe->values[i*NDOF4+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[0])*basis[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])*basis[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])*basis[i];
 		}
 	}
@@ -2972,5 +2972,5 @@
 	double      driving_stress_baseline,thickness;
 	double      xyz_list[NUMVERTICES][3];
-	double      l1l6[6];
+	double      basis[6];
 	GaussPenta  *gauss=NULL;
 
@@ -2990,5 +2990,5 @@
 
 		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
-		GetNodalFunctionsP1(l1l6, gauss);
+		GetNodalFunctionsP1(basis, gauss);
 
 		thickness_input->GetParameterValue(&thickness, gauss);
@@ -2997,5 +2997,5 @@
 		driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG();
 
-		for(i=0;i<NUMVERTICES;i++) for(j=0;j<NDOF2;j++) pe->values[i*NDOF2+j]+= -driving_stress_baseline*slope[j]*Jdet*gauss->weight*l1l6[i];
+		for(i=0;i<NUMVERTICES;i++) for(j=0;j<NDOF2;j++) pe->values[i*NDOF2+j]+= -driving_stress_baseline*slope[j]*Jdet*gauss->weight*basis[i];
 	}
 
@@ -3110,5 +3110,5 @@
 	double		bed_normal[3];
 	double      dz[3];
-	double      l1l6[6]; //for the six nodes of the penta
+	double      basis[6]; //for the six nodes of the penta
 	double      Jdet2d;
 	GaussPenta  *gauss=NULL;
@@ -3139,5 +3139,5 @@
 
 		GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
-		GetNodalFunctionsP1(l1l6, gauss);
+		GetNodalFunctionsP1(basis, gauss);
 
 		BedNormal(&bed_normal[0],xyz_list_tria);
@@ -3155,5 +3155,5 @@
 		water_pressure=gravity*rho_water*bed;
 
-		for(i=0;i<NUMVERTICES;i++) for(j=0;j<3;j++) pe->values[i*NDOF4+j]+=(water_pressure+damper)*gauss->weight*Jdet2d*l1l6[i]*bed_normal[j];
+		for(i=0;i<NUMVERTICES;i++) for(j=0;j<3;j++) pe->values[i*NDOF4+j]+=(water_pressure+damper)*gauss->weight*Jdet2d*basis[i]*bed_normal[j];
 	}
 
@@ -3204,5 +3204,5 @@
 	double     dudx,dvdy,dwdz;
 	double     du[3],dv[3],dw[3];
-	double     l1l6[6];
+	double     basis[6];
 	GaussPenta *gauss=NULL;
 
@@ -3227,5 +3227,5 @@
 
 		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
-		GetNodalFunctionsP1(l1l6, gauss);
+		GetNodalFunctionsP1(basis, gauss);
 
 		vx_input->GetParameterDerivativeValue(&du[0],&xyz_list[0][0],gauss);
@@ -3239,5 +3239,5 @@
 		dvdy=dv[1];
 
-		for (i=0;i<numdof;i++) pe->values[i] += (dudx+dvdy+dwdz)*Jdet*gauss->weight*l1l6[i];
+		for (i=0;i<numdof;i++) pe->values[i] += (dudx+dvdy+dwdz)*Jdet*gauss->weight*basis[i];
 	}
 
@@ -3262,5 +3262,5 @@
 	double     vx,vy,vz,dbdx,dbdy,basalmeltingvalue;
 	double     slope[3];
-	double     l1l6[NUMVERTICES];
+	double     basis[NUMVERTICES];
 	GaussPenta* gauss=NULL;
 
@@ -3302,7 +3302,7 @@
 
 		GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
-		GetNodalFunctionsP1(&l1l6[0], gauss);
-
-		for(i=0;i<numdof;i++) pe->values[i]+=-Jdet2d*gauss->weight*(vx*dbdx+vy*dbdy-vz-basalmeltingvalue)*l1l6[i];
+		GetNodalFunctionsP1(&basis[0], gauss);
+
+		for(i=0;i<numdof;i++) pe->values[i]+=-Jdet2d*gauss->weight*(vx*dbdx+vy*dbdy-vz-basalmeltingvalue)*basis[i];
 	}
 
@@ -3347,5 +3347,5 @@
 	double xyz_list[NUMVERTICES][3];
 	double L[numdof];
-	double dh1dh6[3][6];
+	double dbasis[3][6];
 	double epsilon[6];
 	GaussPenta *gauss=NULL;
@@ -3394,5 +3394,5 @@
 
 		if(artdiff==2){
-			GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss);
+			GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
 
 			vx_input->GetParameterValue(&u, gauss);
@@ -3402,7 +3402,7 @@
 			tau_parameter=GetStabilizationParameter(u,v,w,diameter,rho_ice,heatcapacity,thermalconductivity);
 
-			for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=tau_parameter*scalar_def*(u*dh1dh6[0][i]+v*dh1dh6[1][i]+w*dh1dh6[2][i]);
+			for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=tau_parameter*scalar_def*(u*dbasis[0][i]+v*dbasis[1][i]+w*dbasis[2][i]);
 			if(dt){
-				for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=tau_parameter*scalar_transient*(u*dh1dh6[0][i]+v*dh1dh6[1][i]+w*dh1dh6[2][i]);
+				for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=tau_parameter*scalar_transient*(u*dbasis[0][i]+v*dbasis[1][i]+w*dbasis[2][i]);
 			}
 		}
@@ -3428,5 +3428,5 @@
 	double     xyz_list[NUMVERTICES][3];
 	double     xyz_list_tria[NUMVERTICES2D][3];
-	double     l1l6[NUMVERTICES];
+	double     basis[NUMVERTICES];
 	GaussPenta* gauss=NULL;
 
@@ -3455,5 +3455,5 @@
 
 		GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
-		GetNodalFunctionsP1(&l1l6[0], gauss);
+		GetNodalFunctionsP1(&basis[0], gauss);
 
 		pressure_input->GetParameterValue(&pressure,gauss);
@@ -3463,5 +3463,5 @@
 		if(dt) scalar_ocean=dt*scalar_ocean;
 
-		for(i=0;i<numdof;i++) pe->values[i]+=scalar_ocean*l1l6[i];
+		for(i=0;i<numdof;i++) pe->values[i]+=scalar_ocean*basis[i];
 	}
 
@@ -3486,5 +3486,5 @@
 	double     basalfriction,alpha2,vx,vy;
 	double     scalar;
-	double     l1l6[NUMVERTICES];
+	double     basis[NUMVERTICES];
 	Friction*  friction=NULL;
 	GaussPenta* gauss=NULL;
@@ -3520,5 +3520,5 @@
 
 		GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
-		GetNodalFunctionsP1(&l1l6[0], gauss);
+		GetNodalFunctionsP1(&basis[0], gauss);
 
 		geothermalflux_input->GetParameterValue(&geothermalflux_value,gauss);
@@ -3531,5 +3531,5 @@
 		if(dt) scalar=dt*scalar;
 
-		for(i=0;i<numdof;i++) pe->values[i]+=scalar*l1l6[i];
+		for(i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
 	}
 
@@ -3616,5 +3616,5 @@
 	double xyz_list[NUMVERTICES][3];
 	double L[numdof];
-	double dh1dh6[3][6];
+	double dbasis[3][6];
 	double epsilon[6];
 	GaussPenta *gauss=NULL;
@@ -3663,5 +3663,5 @@
 
 		if(artdiff==2){
-			GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss);
+			GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
 
 			vx_input->GetParameterValue(&u, gauss);
@@ -3671,7 +3671,7 @@
 			tau_parameter=GetStabilizationParameter(u,v,w,diameter,rho_ice,heatcapacity,thermalconductivity);
 
-			for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=tau_parameter*scalar_def*(u*dh1dh6[0][i]+v*dh1dh6[1][i]+w*dh1dh6[2][i]);
+			for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=tau_parameter*scalar_def*(u*dbasis[0][i]+v*dbasis[1][i]+w*dbasis[2][i]);
 			if(dt){
-				for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=tau_parameter*scalar_transient*(u*dh1dh6[0][i]+v*dh1dh6[1][i]+w*dh1dh6[2][i]);
+				for(i=0;i<NUMVERTICES;i++)  pe->values[i]+=tau_parameter*scalar_transient*(u*dbasis[0][i]+v*dbasis[1][i]+w*dbasis[2][i]);
 			}
 		}
@@ -3697,5 +3697,5 @@
 	double     xyz_list[NUMVERTICES][3];
 	double     xyz_list_tria[NUMVERTICES2D][3];
-	double     l1l6[NUMVERTICES];
+	double     basis[NUMVERTICES];
 	GaussPenta* gauss=NULL;
 
@@ -3724,5 +3724,5 @@
 
 		GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
-		GetNodalFunctionsP1(&l1l6[0], gauss);
+		GetNodalFunctionsP1(&basis[0], gauss);
 
 		pressure_input->GetParameterValue(&pressure,gauss);
@@ -3732,5 +3732,5 @@
 		if(dt) scalar_ocean=dt*scalar_ocean;
 
-		for(i=0;i<numdof;i++) pe->values[i]+=scalar_ocean*l1l6[i];
+		for(i=0;i<numdof;i++) pe->values[i]+=scalar_ocean*basis[i];
 	}
 
@@ -3755,5 +3755,5 @@
 	double     basalfriction,alpha2,vx,vy;
 	double     scalar;
-	double     l1l6[NUMVERTICES];
+	double     basis[NUMVERTICES];
 	Friction*  friction=NULL;
 	GaussPenta* gauss=NULL;
@@ -3789,5 +3789,5 @@
 
 		GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
-		GetNodalFunctionsP1(&l1l6[0], gauss);
+		GetNodalFunctionsP1(&basis[0], gauss);
 
 		geothermalflux_input->GetParameterValue(&geothermalflux_value,gauss);
@@ -3800,5 +3800,5 @@
 		if(dt) scalar=dt*scalar;
 
-		for(i=0;i<numdof;i++) pe->values[i]+=scalar*l1l6[i];
+		for(i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
 	}
 
Index: /issm/trunk/src/c/objects/Elements/PentaRef.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/PentaRef.cpp	(revision 8652)
+++ /issm/trunk/src/c/objects/Elements/PentaRef.cpp	(revision 8653)
@@ -70,22 +70,20 @@
 	 */
 
-	double dh1dh6[3][NUMNODESP1];
-
-	/*Get dh1dh6 in actual coordinate system: */
-	GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);
+	double dbasis[3][NUMNODESP1];
+
+	/*Get dbasis in actual coordinate system: */
+	GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss);
 
 	/*Build B: */
 	for (int i=0;i<NUMNODESP1;i++){
-		*(B+NDOF2*NUMNODESP1*0+NDOF2*i)=dh1dh6[0][i]; 
+		*(B+NDOF2*NUMNODESP1*0+NDOF2*i)=dbasis[0][i]; 
 		*(B+NDOF2*NUMNODESP1*0+NDOF2*i+1)=0.0;
 
 		*(B+NDOF2*NUMNODESP1*1+NDOF2*i)=0.0;
-		*(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=dh1dh6[1][i];
-
-		*(B+NDOF2*NUMNODESP1*2+NDOF2*i)=(float).5*dh1dh6[1][i]; 
-		*(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=(float).5*dh1dh6[0][i]; 
-
-	}
-
+		*(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=dbasis[1][i];
+
+		*(B+NDOF2*NUMNODESP1*2+NDOF2*i)=(float).5*dbasis[1][i]; 
+		*(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=(float).5*dbasis[0][i]; 
+	}
 }
 /*}}}*/
@@ -151,25 +149,25 @@
 	 */
 
-	double dh1dh6[3][NUMNODESP1];
-
-	/*Get dh1dh6 in actual coordinate system: */
-	GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss);
+	double dbasis[3][NUMNODESP1];
+
+	/*Get dbasis in actual coordinate system: */
+	GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss);
 
 	/*Build B: */
 	for (int i=0;i<NUMNODESP1;i++){
-		*(B+NDOF2*NUMNODESP1*0+NDOF2*i)=dh1dh6[0][i]; 
+		*(B+NDOF2*NUMNODESP1*0+NDOF2*i)=dbasis[0][i]; 
 		*(B+NDOF2*NUMNODESP1*0+NDOF2*i+1)=0.0;
 
 		*(B+NDOF2*NUMNODESP1*1+NDOF2*i)=0.0;
-		*(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=dh1dh6[1][i];
-
-		*(B+NDOF2*NUMNODESP1*2+NDOF2*i)=(float).5*dh1dh6[1][i]; 
-		*(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=(float).5*dh1dh6[0][i]; 
-
-		*(B+NDOF2*NUMNODESP1*3+NDOF2*i)=(float).5*dh1dh6[2][i]; 
+		*(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=dbasis[1][i];
+
+		*(B+NDOF2*NUMNODESP1*2+NDOF2*i)=(float).5*dbasis[1][i]; 
+		*(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=(float).5*dbasis[0][i]; 
+
+		*(B+NDOF2*NUMNODESP1*3+NDOF2*i)=(float).5*dbasis[2][i]; 
 		*(B+NDOF2*NUMNODESP1*3+NDOF2*i+1)=0.0;
 
 		*(B+NDOF2*NUMNODESP1*4+NDOF2*i)=0.0;
-		*(B+NDOF2*NUMNODESP1*4+NDOF2*i+1)=(float).5*dh1dh6[2][i]; 
+		*(B+NDOF2*NUMNODESP1*4+NDOF2*i+1)=(float).5*dbasis[2][i]; 
 	}
 
@@ -190,25 +188,25 @@
 	 * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1)
 	 */
-	double dh1dh6[3][NUMNODESP1];
-
-	/*Get dh1dh6 in actual coordinate system: */
-	GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list, gauss_coord);
+	double dbasis[3][NUMNODESP1];
+
+	/*Get dbasis in actual coordinate system: */
+	GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss_coord);
 
 	/*Build BPrime: */
 	for (int i=0;i<NUMNODESP1;i++){
-		*(B+NDOF2*NUMNODESP1*0+NDOF2*i)=2.0*dh1dh6[0][i]; 
-		*(B+NDOF2*NUMNODESP1*0+NDOF2*i+1)=dh1dh6[1][i];
-
-		*(B+NDOF2*NUMNODESP1*1+NDOF2*i)=dh1dh6[0][i];
-		*(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=2.0*dh1dh6[1][i];
-
-		*(B+NDOF2*NUMNODESP1*2+NDOF2*i)=dh1dh6[1][i]; 
-		*(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=dh1dh6[0][i]; 
-
-		*(B+NDOF2*NUMNODESP1*3+NDOF2*i)=dh1dh6[2][i]; 
+		*(B+NDOF2*NUMNODESP1*0+NDOF2*i)=2.0*dbasis[0][i]; 
+		*(B+NDOF2*NUMNODESP1*0+NDOF2*i+1)=dbasis[1][i];
+
+		*(B+NDOF2*NUMNODESP1*1+NDOF2*i)=dbasis[0][i];
+		*(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=2.0*dbasis[1][i];
+
+		*(B+NDOF2*NUMNODESP1*2+NDOF2*i)=dbasis[1][i]; 
+		*(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=dbasis[0][i]; 
+
+		*(B+NDOF2*NUMNODESP1*3+NDOF2*i)=dbasis[2][i]; 
 		*(B+NDOF2*NUMNODESP1*3+NDOF2*i+1)=0.0;
 
 		*(B+NDOF2*NUMNODESP1*4+NDOF2*i)=0.0;
-		*(B+NDOF2*NUMNODESP1*4+NDOF2*i+1)=dh1dh6[2][i]; 
+		*(B+NDOF2*NUMNODESP1*4+NDOF2*i+1)=dbasis[2][i]; 
 	}
 }
@@ -401,13 +399,13 @@
 
 	/*Same thing in the actual coordinate system: */
-	double dh1dh6[3][NUMNODESP1];
-
-	/*Get dh1dh6 in actual coordinates system : */
-	GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);
+	double dbasis[3][NUMNODESP1];
+
+	/*Get dbasis in actual coordinates system : */
+	GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);
 
 	/*Build B': */
 	for (int i=0;i<NUMNODESP1;i++){
-		*(B_artdiff+NDOF1*NUMNODESP1*0+NDOF1*i)=dh1dh6[0][i]; 
-		*(B_artdiff+NDOF1*NUMNODESP1*1+NDOF1*i)=dh1dh6[1][i]; 
+		*(B_artdiff+NDOF1*NUMNODESP1*0+NDOF1*i)=dbasis[0][i]; 
+		*(B_artdiff+NDOF1*NUMNODESP1*1+NDOF1*i)=dbasis[1][i]; 
 	}
 }
Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 8652)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 8653)
@@ -1538,5 +1538,5 @@
 	double         xyz_list[NUMVERTICES][3];
 	double         slope[2];
-	double         l1l2l3[3];
+	double         basis[3];
 	double         pe_g_gaussian[numdof];
 	GaussTria*     gauss=NULL;
@@ -1559,5 +1559,5 @@
 
 		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
-		GetNodalFunctions(l1l2l3, gauss);
+		GetNodalFunctions(basis, gauss);
 
 		thickness_input->GetParameterValue(&thickness,gauss);
@@ -1574,5 +1574,5 @@
 			for (i=0;i<NUMVERTICES;i++){
 				for (j=0;j<NDOF2;j++){
-					pe->values[i*NDOF2+j]+=(-driving_stress_baseline*slope[j]-plastic_stress)*Jdet*gauss->weight*l1l2l3[i]; 
+					pe->values[i*NDOF2+j]+=(-driving_stress_baseline*slope[j]-plastic_stress)*Jdet*gauss->weight*basis[i]; 
 				}
 			}
@@ -1581,5 +1581,5 @@
 			for (i=0;i<NUMVERTICES;i++){
 				for (j=0;j<NDOF2;j++){
-					pe->values[i*NDOF2+j]+=-driving_stress_baseline*slope[j]*Jdet*gauss->weight*l1l2l3[i];
+					pe->values[i*NDOF2+j]+=-driving_stress_baseline*slope[j]*Jdet*gauss->weight*basis[i];
 				}
 			}
@@ -1605,5 +1605,5 @@
 	int         num_responses;
 	double      xyz_list[NUMVERTICES][3];
-	double      l1l2l3[3];
+	double      basis[3];
 	double      dbasis[NDOF2][NUMVERTICES];
 	double      dH[2];
@@ -1628,5 +1628,5 @@
 
 		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
-		GetNodalFunctions(l1l2l3, gauss);
+		GetNodalFunctions(basis, gauss);
 		GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
 
@@ -1640,5 +1640,5 @@
 			case ThicknessAbsMisfitEnum:
 				weights_input->GetParameterValue(&weight, gauss,resp);
-				for(i=0;i<numdof;i++) pe->values[i]+=(thicknessobs-thickness)*weight*Jdet*gauss->weight*l1l2l3[i];
+				for(i=0;i<numdof;i++) pe->values[i]+=(thicknessobs-thickness)*weight*Jdet*gauss->weight*basis[i];
 				break;
 			case ThicknessAbsGradientEnum:
@@ -1673,5 +1673,5 @@
 	double     vx,vy,vxobs,vyobs,weight;
 	double     xyz_list[NUMVERTICES][3];
-	double     l1l2l3[3];
+	double     basis[3];
 	GaussTria* gauss=NULL;
 
@@ -1712,5 +1712,5 @@
 		vxobs_input->GetParameterValue(&vxobs,gauss);
 		vyobs_input->GetParameterValue(&vyobs,gauss);
-		GetNodalFunctions(l1l2l3, gauss);
+		GetNodalFunctions(basis, gauss);
 
 		/*Loop over all requested responses*/
@@ -1733,6 +1733,6 @@
 						dux=vxobs-vx;
 						duy=vyobs-vy;
-						pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 
-						pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 
+						pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
+						pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
 					}
 					break;
@@ -1754,6 +1754,6 @@
 						dux=scalex*(vxobs-vx);
 						duy=scaley*(vyobs-vy);
-						pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 
-						pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 
+						pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
+						pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
 					}
 					break;
@@ -1776,6 +1776,6 @@
 						dux=scale*vx;
 						duy=scale*vy;
-						pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 
-						pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 
+						pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
+						pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
 					}
 					break;
@@ -1794,6 +1794,6 @@
 						dux=scale*(vxobs-vx);
 						duy=scale*(vyobs-vy);
-						pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 
-						pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 
+						pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
+						pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
 					}
 					break;
@@ -1811,6 +1811,6 @@
 						dux = - pow(meanvel,2.) * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
 						duy = - pow(meanvel,2.) * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
-						pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 
-						pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 
+						pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
+						pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
 					}
 					break;
@@ -1849,5 +1849,5 @@
 	double     vx,vy,vxobs,vyobs,weight;
 	double     xyz_list[NUMVERTICES][3];
-	double     l1l2l3[3];
+	double     basis[3];
 	GaussTria* gauss=NULL;
 
@@ -1888,5 +1888,5 @@
 		vxobs_input->GetParameterValue(&vxobs,gauss);
 		vyobs_input->GetParameterValue(&vyobs,gauss);
-		GetNodalFunctions(l1l2l3, gauss);
+		GetNodalFunctions(basis, gauss);
 
 		/*Loop over all requested responses*/
@@ -1910,6 +1910,6 @@
 						dux=vxobs-vx;
 						duy=vyobs-vy;
-						pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 
-						pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 
+						pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
+						pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
 					}
 					break;
@@ -1931,6 +1931,6 @@
 						dux=scalex*(vxobs-vx);
 						duy=scaley*(vyobs-vy);
-						pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 
-						pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 
+						pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
+						pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
 					}
 					break;
@@ -1953,6 +1953,6 @@
 						dux=scale*vx;
 						duy=scale*vy;
-						pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 
-						pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 
+						pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
+						pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
 					}
 					break;
@@ -1971,6 +1971,6 @@
 						dux=scale*(vxobs-vx);
 						duy=scale*(vyobs-vy);
-						pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 
-						pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 
+						pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
+						pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
 					}
 					break;
@@ -1988,6 +1988,6 @@
 						dux = - pow(meanvel,2.) * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
 						duy = - pow(meanvel,2.) * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
-						pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*l1l2l3[i]; 
-						pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*l1l2l3[i]; 
+						pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
+						pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
 					}
 					break;
@@ -2229,5 +2229,5 @@
 	double     xyz_list[NUMVERTICES][3];
 	double     slope[2];
-	double     l1l2l3[3];
+	double     basis[3];
 	GaussTria* gauss=NULL;
 
@@ -2253,13 +2253,13 @@
 
 		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
-		GetNodalFunctions(l1l2l3, gauss);
+		GetNodalFunctions(basis, gauss);
 
 		slope_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
 
 		if ( (analysis_type==SurfaceSlopeXAnalysisEnum) || (analysis_type==BedSlopeXAnalysisEnum)){
-			for(i=0;i<numdof;i++) pe->values[i]+=Jdet*gauss->weight*slope[0]*l1l2l3[i];
+			for(i=0;i<numdof;i++) pe->values[i]+=Jdet*gauss->weight*slope[0]*basis[i];
 		}
 		if ( (analysis_type==SurfaceSlopeYAnalysisEnum) || (analysis_type==BedSlopeYAnalysisEnum)){
-			for(i=0;i<numdof;i++) pe->values[i]+=Jdet*gauss->weight*slope[1]*l1l2l3[i];
+			for(i=0;i<numdof;i++) pe->values[i]+=Jdet*gauss->weight*slope[1]*basis[i];
 		}
 	}
@@ -2949,5 +2949,5 @@
 	double     Jdet,weight;
 	double     xyz_list[NUMVERTICES][3];
-	double     dh1dh3[NDOF2][NUMVERTICES];
+	double     dbasis[NDOF2][NUMVERTICES];
 	double     dk[NDOF2]; 
 	double     grade_g[NUMVERTICES]={0.0};
@@ -2967,5 +2967,5 @@
 
 		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
-		GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss);
+		GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
 		weights_input->GetParameterValue(&weight,gauss,weight_index);
 
@@ -2974,5 +2974,5 @@
 
 		/*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
-		for (i=0;i<NUMVERTICES;i++) grade_g[i]+=-weight*Jdet*gauss->weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);
+		for (i=0;i<NUMVERTICES;i++) grade_g[i]+=-weight*Jdet*gauss->weight*(dbasis[0][i]*dk[0]+dbasis[1][i]*dk[1]);
 	}
 	VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES);
@@ -3050,9 +3050,9 @@
 	double     bed,thickness,Neff,drag;
 	double     xyz_list[NUMVERTICES][3];
-	double     dh1dh3[NDOF2][NUMVERTICES];
+	double     dbasis[NDOF2][NUMVERTICES];
 	double     dk[NDOF2]; 
 	double     grade_g[NUMVERTICES]={0.0};
 	double     grade_g_gaussian[NUMVERTICES];
-	double     l1l2l3[3];
+	double     basis[3];
 	double     epsilon[3]; /* epsilon=[exx,eyy,exy];*/
 	Friction*  friction=NULL;
@@ -3084,6 +3084,6 @@
 
 		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
-		GetNodalFunctions(l1l2l3, gauss);
-		GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss);
+		GetNodalFunctions(basis, gauss);
+		GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
 
 		/*Build alpha_complement_list: */
@@ -3100,5 +3100,5 @@
 		/*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
 		for (i=0;i<NUMVERTICES;i++){
-			grade_g_gaussian[i]=-2*drag*alpha_complement*((lambda*vx+mu*vy))*Jdet*gauss->weight*l1l2l3[i];
+			grade_g_gaussian[i]=-2*drag*alpha_complement*((lambda*vx+mu*vy))*Jdet*gauss->weight*basis[i];
 		}
 		
@@ -3128,7 +3128,7 @@
 	double     surface_normal[3],bed_normal[3];
 	double     xyz_list[NUMVERTICES][3];
-	double     dh1dh3[NDOF2][NUMVERTICES];
+	double     dbasis[NDOF2][NUMVERTICES];
 	double     dk[NDOF2]; 
-	double     l1l2l3[3];
+	double     basis[3];
 	double     epsilon[3]; /* epsilon=[exx,eyy,exy];*/
 	double     grade_g[NUMVERTICES]={0.0};
@@ -3193,8 +3193,8 @@
 
 		/* Get nodal functions value at gaussian point:*/
-		GetNodalFunctions(l1l2l3, gauss);
+		GetNodalFunctions(basis, gauss);
 
 		/*Get nodal functions derivatives*/
-		GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss);
+		GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
 
 		/*Get k derivative: dk/dx */
@@ -3208,5 +3208,5 @@
 						-mu    *(2*drag*alpha_complement*(vy - vz*bed_normal[1]*bed_normal[2]))
 						-xi    *(2*drag*alpha_complement*(-vx*bed_normal[0]*bed_normal[2]-vy*bed_normal[1]*bed_normal[2]))
-						)*Jdet*gauss->weight*l1l2l3[i]; 
+						)*Jdet*gauss->weight*basis[i]; 
 		}
 
@@ -3228,5 +3228,5 @@
 	double     Jdet,weight;
 	double     xyz_list[NUMVERTICES][3];
-	double     dh1dh3[NDOF2][NUMVERTICES];
+	double     dbasis[NDOF2][NUMVERTICES];
 	double     dk[NDOF2]; 
 	double     grade_g[NUMVERTICES]={0.0};
@@ -3247,5 +3247,5 @@
 
 		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
-		GetNodalFunctionsDerivatives(&dh1dh3[0][0],&xyz_list[0][0],gauss);
+		GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
 		weights_input->GetParameterValue(&weight,gauss,weight_index);
 
@@ -3255,5 +3255,5 @@
 		/*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
 		for (i=0;i<NUMVERTICES;i++){
-			grade_g[i]+=-weight*Jdet*gauss->weight*(dh1dh3[0][i]*dk[0]+dh1dh3[1][i]*dk[1]);
+			grade_g[i]+=-weight*Jdet*gauss->weight*(dbasis[0][i]*dk[0]+dbasis[1][i]*dk[1]);
 			_assert_(!isnan(grade_g[i]));
 		}
@@ -3289,5 +3289,5 @@
 	int        doflist1[NUMVERTICES];
 	double     thickness,Jdet;
-	double     l1l2l3[3];
+	double     basis[3];
 	double     dbasis[NDOF2][NUMVERTICES];
 	double     Dlambda[2],dp[2];
@@ -3311,5 +3311,5 @@
 
 		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
-		GetNodalFunctions(l1l2l3, gauss);
+		GetNodalFunctions(basis, gauss);
 		GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
 		
@@ -3318,5 +3318,5 @@
 		thickness_input->GetParameterDerivativeValue(&dp[0],&xyz_list[0][0],gauss);
 
-		for(i=0;i<NUMVERTICES;i++) grade_g[i]+=thickness*Dlambda[0]*Jdet*gauss->weight*l1l2l3[i];
+		for(i=0;i<NUMVERTICES;i++) grade_g[i]+=thickness*Dlambda[0]*Jdet*gauss->weight*basis[i];
 	}
 
@@ -3334,5 +3334,5 @@
 	int        doflist1[NUMVERTICES];
 	double     thickness,Jdet;
-	double     l1l2l3[3];
+	double     basis[3];
 	double     dbasis[NDOF2][NUMVERTICES];
 	double     Dlambda[2],dp[2];
@@ -3356,5 +3356,5 @@
 
 		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
-		GetNodalFunctions(l1l2l3, gauss);
+		GetNodalFunctions(basis, gauss);
 		GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
 
@@ -3363,5 +3363,5 @@
 		thickness_input->GetParameterDerivativeValue(&dp[0],&xyz_list[0][0],gauss);
 
-		for(i=0;i<NUMVERTICES;i++) grade_g[i]+=thickness*Dlambda[1]*Jdet*gauss->weight*l1l2l3[i];
+		for(i=0;i<NUMVERTICES;i++) grade_g[i]+=thickness*Dlambda[1]*Jdet*gauss->weight*basis[i];
 	}
 	VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES);
Index: /issm/trunk/src/c/objects/Elements/TriaRef.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/TriaRef.cpp	(revision 8652)
+++ /issm/trunk/src/c/objects/Elements/TriaRef.cpp	(revision 8653)
@@ -70,17 +70,17 @@
 
 	int i;
-	double dh1dh3[NDOF2][NUMNODES];
+	double dbasis[NDOF2][NUMNODES];
 
 	/*Get dh1dh2dh3 in actual coordinate system: */
-	GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list,gauss);
+	GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list,gauss);
 
 	/*Build B: */
 	for (i=0;i<NUMNODES;i++){
-		*(B+NDOF2*NUMNODES*0+NDOF2*i)=dh1dh3[0][i]; //B[0][NDOF2*i]=dh1dh3[0][i];
+		*(B+NDOF2*NUMNODES*0+NDOF2*i)=dbasis[0][i]; //B[0][NDOF2*i]=dbasis[0][i];
 		*(B+NDOF2*NUMNODES*0+NDOF2*i+1)=0;
 		*(B+NDOF2*NUMNODES*1+NDOF2*i)=0;
-		*(B+NDOF2*NUMNODES*1+NDOF2*i+1)=dh1dh3[1][i];
-		*(B+NDOF2*NUMNODES*2+NDOF2*i)=(float).5*dh1dh3[1][i]; 
-		*(B+NDOF2*NUMNODES*2+NDOF2*i+1)=(float).5*dh1dh3[0][i]; 
+		*(B+NDOF2*NUMNODES*1+NDOF2*i+1)=dbasis[1][i];
+		*(B+NDOF2*NUMNODES*2+NDOF2*i)=(float).5*dbasis[1][i]; 
+		*(B+NDOF2*NUMNODES*2+NDOF2*i+1)=(float).5*dbasis[0][i]; 
 	}
 }
@@ -101,17 +101,17 @@
 
 	/*Same thing in the actual coordinate system: */
-	double dh1dh3[NDOF2][NUMNODES];
+	double dbasis[NDOF2][NUMNODES];
 
 	/*Get dh1dh2dh3 in actual coordinates system : */
-	GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list,gauss);
+	GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list,gauss);
 
 	/*Build B': */
 	for (int i=0;i<NUMNODES;i++){
-		*(B+NDOF2*NUMNODES*0+NDOF2*i)=dh1dh3[0][i]; 
+		*(B+NDOF2*NUMNODES*0+NDOF2*i)=dbasis[0][i]; 
 		*(B+NDOF2*NUMNODES*0+NDOF2*i+1)=0; 
 		*(B+NDOF2*NUMNODES*1+NDOF2*i)=0; 
-		*(B+NDOF2*NUMNODES*1+NDOF2*i+1)=dh1dh3[1][i]; 
-		*(B+NDOF2*NUMNODES*2+NDOF2*i)=0.5*dh1dh3[1][i]; 
-		*(B+NDOF2*NUMNODES*2+NDOF2*i+1)=0.5*dh1dh3[0][i]; 
+		*(B+NDOF2*NUMNODES*1+NDOF2*i+1)=dbasis[1][i]; 
+		*(B+NDOF2*NUMNODES*2+NDOF2*i)=0.5*dbasis[1][i]; 
+		*(B+NDOF2*NUMNODES*2+NDOF2*i+1)=0.5*dbasis[0][i]; 
 	}
 }
@@ -167,13 +167,13 @@
 	 */
 
-	double l1l2l3[NUMNODES];
+	double basis[NUMNODES];
 
 	/*Get dh1dh2dh3 in actual coordinate system: */
-	GetNodalFunctions(&l1l2l3[0],gauss);
+	GetNodalFunctions(&basis[0],gauss);
 
 	/*Build B_prog: */
 	for (int i=0;i<NUMNODES;i++){
-		*(B_prog+NDOF1*NUMNODES*0+NDOF1*i)=l1l2l3[i];
-		*(B_prog+NDOF1*NUMNODES*1+NDOF1*i)=l1l2l3[i];
+		*(B_prog+NDOF1*NUMNODES*0+NDOF1*i)=basis[i];
+		*(B_prog+NDOF1*NUMNODES*1+NDOF1*i)=basis[i];
 	}
 }
@@ -194,17 +194,17 @@
 
 	/*Same thing in the actual coordinate system: */
-	double dh1dh3[NDOF2][NUMNODES];
+	double dbasis[NDOF2][NUMNODES];
 
 	/*Get dh1dh2dh3 in actual coordinates system : */
-	GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list,gauss);
+	GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list,gauss);
 
 	/*Build B': */
 	for (int i=0;i<NUMNODES;i++){
-		*(Bprime+NDOF2*NUMNODES*0+NDOF2*i)=2*dh1dh3[0][i]; 
-		*(Bprime+NDOF2*NUMNODES*0+NDOF2*i+1)=dh1dh3[1][i]; 
-		*(Bprime+NDOF2*NUMNODES*1+NDOF2*i)=dh1dh3[0][i]; 
-		*(Bprime+NDOF2*NUMNODES*1+NDOF2*i+1)=2*dh1dh3[1][i]; 
-		*(Bprime+NDOF2*NUMNODES*2+NDOF2*i)=dh1dh3[1][i]; 
-		*(Bprime+NDOF2*NUMNODES*2+NDOF2*i+1)=dh1dh3[0][i]; 
+		*(Bprime+NDOF2*NUMNODES*0+NDOF2*i)=2*dbasis[0][i]; 
+		*(Bprime+NDOF2*NUMNODES*0+NDOF2*i+1)=dbasis[1][i]; 
+		*(Bprime+NDOF2*NUMNODES*1+NDOF2*i)=dbasis[0][i]; 
+		*(Bprime+NDOF2*NUMNODES*1+NDOF2*i+1)=2*dbasis[1][i]; 
+		*(Bprime+NDOF2*NUMNODES*2+NDOF2*i)=dbasis[1][i]; 
+		*(Bprime+NDOF2*NUMNODES*2+NDOF2*i+1)=dbasis[0][i]; 
 	}
 }
@@ -226,19 +226,19 @@
 
 	/*Same thing in the actual coordinate system: */
-	double dh1dh3[NDOF2][NUMNODES];
+	double dbasis[NDOF2][NUMNODES];
 
 	/*Get dh1dh2dh3 in actual coordinates system : */
-	GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list,gauss);
+	GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list,gauss);
 
 	/*Build Bprime: */
 	for (int i=0;i<NUMNODES;i++){
-		*(Bprime+NDOF2*NUMNODES*0+NDOF2*i)=dh1dh3[0][i]; 
+		*(Bprime+NDOF2*NUMNODES*0+NDOF2*i)=dbasis[0][i]; 
 		*(Bprime+NDOF2*NUMNODES*0+NDOF2*i+1)=0; 
 		*(Bprime+NDOF2*NUMNODES*1+NDOF2*i)=0; 
-		*(Bprime+NDOF2*NUMNODES*1+NDOF2*i+1)=dh1dh3[1][i]; 
-		*(Bprime+NDOF2*NUMNODES*2+NDOF2*i)=dh1dh3[1][i]; 
-		*(Bprime+NDOF2*NUMNODES*2+NDOF2*i+1)=dh1dh3[0][i]; 
-		*(Bprime+NDOF2*NUMNODES*3+NDOF2*i)=dh1dh3[0][i]; 
-		*(Bprime+NDOF2*NUMNODES*3+NDOF2*i+1)=dh1dh3[1][i]; 
+		*(Bprime+NDOF2*NUMNODES*1+NDOF2*i+1)=dbasis[1][i]; 
+		*(Bprime+NDOF2*NUMNODES*2+NDOF2*i)=dbasis[1][i]; 
+		*(Bprime+NDOF2*NUMNODES*2+NDOF2*i+1)=dbasis[0][i]; 
+		*(Bprime+NDOF2*NUMNODES*3+NDOF2*i)=dbasis[0][i]; 
+		*(Bprime+NDOF2*NUMNODES*3+NDOF2*i+1)=dbasis[1][i]; 
 	}
 }
@@ -257,13 +257,13 @@
 
 	/*Same thing in the actual coordinate system: */
-	double dh1dh3[NDOF2][NUMNODES];
+	double dbasis[NDOF2][NUMNODES];
 
 	/*Get dh1dh2dh3 in actual coordinates system : */
-	GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list,gauss);
+	GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list,gauss);
 
 	/*Build B': */
 	for (int i=0;i<NUMNODES;i++){
-		*(Bprime_prog+NDOF1*NUMNODES*0+NDOF1*i)=dh1dh3[0][i]; 
-		*(Bprime_prog+NDOF1*NUMNODES*1+NDOF1*i)=dh1dh3[1][i]; 
+		*(Bprime_prog+NDOF1*NUMNODES*0+NDOF1*i)=dbasis[0][i]; 
+		*(Bprime_prog+NDOF1*NUMNODES*1+NDOF1*i)=dbasis[1][i]; 
 	}
 }
@@ -285,21 +285,21 @@
 
 	int i;
-	double l1l2l3[3];
-
-	/*Get l1l2l3 in actual coordinate system: */
-	GetNodalFunctions(l1l2l3,gauss);
+	double basis[3];
+
+	/*Get basis in actual coordinate system: */
+	GetNodalFunctions(basis,gauss);
 
 	/*Build L: */
 	if(numdof==1){
 		for (i=0;i<NUMNODES;i++){
-			L[i]=l1l2l3[i]; 
+			L[i]=basis[i]; 
 		}
 	}
 	else{
 		for (i=0;i<NUMNODES;i++){
-			*(L+numdof*NUMNODES*0+numdof*i)=l1l2l3[i]; //L[0][NDOF2*i]=dh1dh3[0][i];
+			*(L+numdof*NUMNODES*0+numdof*i)=basis[i]; //L[0][NDOF2*i]=dbasis[0][i];
 			*(L+numdof*NUMNODES*0+numdof*i+1)=0;
 			*(L+numdof*NUMNODES*1+numdof*i)=0;
-			*(L+numdof*NUMNODES*1+numdof*i+1)=l1l2l3[i];
+			*(L+numdof*NUMNODES*1+numdof*i+1)=basis[i];
 		}
 	}
@@ -394,15 +394,15 @@
 /*}}}*/
 /*FUNCTION TriaRef::GetNodalFunctions{{{1*/
-void TriaRef::GetNodalFunctions(double* l1l2l3,GaussTria* gauss){
+void TriaRef::GetNodalFunctions(double* basis,GaussTria* gauss){
 	/*This routine returns the values of the nodal functions  at the gaussian point.*/
 
-	l1l2l3[0]=gauss->coord1;
-	l1l2l3[1]=gauss->coord2;
-	l1l2l3[2]=gauss->coord3;
+	basis[0]=gauss->coord1;
+	basis[1]=gauss->coord2;
+	basis[2]=gauss->coord3;
 
 }
 /*}}}*/
 /*FUNCTION TriaRef::GetSegmentNodalFunctions{{{1*/
-void TriaRef::GetSegmentNodalFunctions(double* l1l2,GaussTria* gauss,int index1,int index2){
+void TriaRef::GetSegmentNodalFunctions(double* basis,GaussTria* gauss,int index1,int index2){
 	/*This routine returns the values of the nodal functions  at the gaussian point.*/
 
@@ -413,24 +413,24 @@
 	_assert_(index1>=0 && index1<3);
 	_assert_(index2>=0 && index2<3);
-	l1l2[0]=BasisFunctions[index1];
-	l1l2[1]=BasisFunctions[index2];
+	basis[0]=BasisFunctions[index1];
+	basis[1]=BasisFunctions[index2];
 }
 /*}}}*/
 /*FUNCTION TriaRef::GetNodalFunctionsDerivatives{{{1*/
-void TriaRef::GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, GaussTria* gauss){
+void TriaRef::GetNodalFunctionsDerivatives(double* dbasis,double* xyz_list, GaussTria* gauss){
 
 	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
 	 * actual coordinate system): */
 	int       i;
-	double    dh1dh3_ref[NDOF2][NUMNODES];
+	double    dbasis_ref[NDOF2][NUMNODES];
 	double    Jinv[NDOF2][NDOF2];
 
 	/*Get derivative values with respect to parametric coordinate system: */
-	GetNodalFunctionsDerivativesReference(&dh1dh3_ref[0][0], gauss); 
+	GetNodalFunctionsDerivativesReference(&dbasis_ref[0][0], gauss); 
 
 	/*Get Jacobian invert: */
 	GetJacobianInvert(&Jinv[0][0], xyz_list, gauss);
 
-	/*Build dh1dh3: 
+	/*Build dbasis: 
 	 *
 	 * [dhi/dx]= Jinv*[dhi/dr]
@@ -438,6 +438,6 @@
 	 */
 	for (i=0;i<NUMNODES;i++){
-		dh1dh3[NUMNODES*0+i]=Jinv[0][0]*dh1dh3_ref[0][i]+Jinv[0][1]*dh1dh3_ref[1][i];
-		dh1dh3[NUMNODES*1+i]=Jinv[1][0]*dh1dh3_ref[0][i]+Jinv[1][1]*dh1dh3_ref[1][i];
+		dbasis[NUMNODES*0+i]=Jinv[0][0]*dbasis_ref[0][i]+Jinv[0][1]*dbasis_ref[1][i];
+		dbasis[NUMNODES*1+i]=Jinv[1][0]*dbasis_ref[0][i]+Jinv[1][1]*dbasis_ref[1][i];
 	}
 
@@ -467,5 +467,5 @@
 
 	/*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter derivative value at gaussian 
-	 * point specified by gauss_l1l2l3:
+	 * point specified by gauss_basis:
 	 *   dp/dx=plist[0]*dh1/dx+plist[1]*dh2/dx+plist[2]*dh3/dx
 	 *   dp/dx=plist[0]*dh1/dx+plist[1]*dh2/dx+plist[2]*dh3/dx
@@ -475,12 +475,12 @@
 
 	/*Nodal Derivatives*/
-	double dh1dh3[2][3]; //nodal derivative functions in actual coordinate system.
+	double dbasis[2][3]; //nodal derivative functions in actual coordinate system.
 
 	/*Get dh1dh2dh3 in actual coordinate system: */
-	GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list, gauss);
+	GetNodalFunctionsDerivatives(&dbasis[0][0],xyz_list, gauss);
 
 	/*Assign values*/
-	*(p+0)=plist[0]*dh1dh3[0][0]+plist[1]*dh1dh3[0][1]+plist[2]*dh1dh3[0][2];
-	*(p+1)=plist[0]*dh1dh3[1][0]+plist[1]*dh1dh3[1][1]+plist[2]*dh1dh3[1][2];
+	*(p+0)=plist[0]*dbasis[0][0]+plist[1]*dbasis[0][1]+plist[2]*dbasis[0][2];
+	*(p+1)=plist[0]*dbasis[1][0]+plist[1]*dbasis[1][1]+plist[2]*dbasis[1][2];
 
 }
@@ -493,11 +493,11 @@
 
 	/*nodal functions annd output: */
-	double l1l2l3[3];
+	double basis[3];
 
 	/*Get nodal functions*/
-	GetNodalFunctions(l1l2l3, gauss);
+	GetNodalFunctions(basis, gauss);
 
 	/*Get parameter*/
-	*p=l1l2l3[0]*plist[0]+l1l2l3[1]*plist[1]+l1l2l3[2]*plist[2];
-}
-/*}}}*/
+	*p=basis[0]*plist[0]+basis[1]*plist[1]+basis[2]*plist[2];
+}
+/*}}}*/