Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15312)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15313)
@@ -6060,5 +6060,5 @@
 
 		GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
-		GetL(&L[0][0], gauss,NDOF2);
+		GetBPattynFriction(&L[0][0],gauss);
 
 		DL_scalar=alpha2*gauss->weight*Jdet2d;
@@ -6873,5 +6873,5 @@
 
 		GetTriaJacobianDeterminant(&Jdet, &xyz_list_tria[0][0],gauss);
-		GetL(&L[0][0], gauss,NDOF2);
+		GetBPattynFriction(&L[0][0],gauss);
 
 		friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); 
Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 15312)
+++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 15313)
@@ -488,38 +488,26 @@
 }
 /*}}}*/
-/*FUNCTION PentaRef::GetL{{{*/
-void PentaRef::GetL(IssmDouble* L, GaussPenta* gauss, int numdof){
-	/*Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof. 
-	 ** For node i, Li can be expressed in the actual coordinate system
+/*FUNCTION PentaRef::GetBPattynFriction{{{*/
+void PentaRef::GetBPattynFriction(IssmDouble* B, GaussPenta* gauss){
+	/*Compute B  matrix. B=[B1 B2 B3] where Bi is square and of size 2x2. 
+	 ** For node i, Bi can be expressed in the actual coordinate system
 	 ** by: 
-	 **       numdof=1: 
-	 **                 Li=h;
-	 **       numdof=2:
-	 **                 Li=[ h   0 ]
-	 **                    [ 0   h ]
-	 ** where h is the interpolation function for node i.
+	 **                 Bi=[ N   0 ]
+	 **                    [ 0   N ]
+	 ** where N is the interpolation function for node i.
 	 **
-	 ** We assume L has been allocated already, of size: NUMNODESP1 (numdof=1), or numdofx(numdof*NUMNODESP1) (numdof=2)
+	 ** We assume B has been allocated already, of size: 2 (2 x numnodes)
 	 **/
 
-	int i;
-	IssmDouble l1l6[6];
+	IssmDouble basis[6];
 
 	/*Get l1l6 in actual coordinate system: */
-	GetNodalFunctionsP1(l1l6,gauss);
-
-	/*Build L: */
-	if(numdof==1){
-		for (i=0;i<NUMNODESP1;i++){
-			L[i]=l1l6[i]; 
-		}
-	}
-	else{
-		for (i=0;i<NUMNODESP1;i++){
-			*(L+numdof*NUMNODESP1*0+numdof*i)=l1l6[i]; 
-			*(L+numdof*NUMNODESP1*0+numdof*i+1)=0;
-			*(L+numdof*NUMNODESP1*1+numdof*i)=0;
-			*(L+numdof*NUMNODESP1*1+numdof*i+1)=l1l6[i];
-		}
+	GetNodalFunctionsP1(&basis[0],gauss);
+
+	for(int i=0;i<NUMNODESP1;i++){
+		B[2*NUMNODESP1*0+2*i+0]=basis[i]; 
+		B[2*NUMNODESP1*0+2*i+1]=0.;
+		B[2*NUMNODESP1*1+2*i+0]=0.;
+		B[2*NUMNODESP1*1+2*i+1]=basis[i];
 	}
 } 
Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h	(revision 15312)
+++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h	(revision 15313)
@@ -47,5 +47,5 @@
 		void GetBVert(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
 		void GetBprimeAdvec(IssmDouble* Bprime_advec, IssmDouble* xyz_list, GaussPenta* gauss);
-		void GetL(IssmDouble* L, GaussPenta* gauss,int numdof);
+		void GetBPattynFriction(IssmDouble* L, GaussPenta* gauss);
 		void GetLStokes(IssmDouble* LStokes, GaussPenta* gauss);
 		void GetLprimeStokes(IssmDouble* LprimeStokes, IssmDouble* xyz_list, GaussPenta* gauss);
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15312)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15313)
@@ -292,5 +292,5 @@
 	IssmDouble Jdet,D_scalar;
 	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble L[3];
+	IssmDouble basis[3];
 	GaussTria *gauss=NULL;
 
@@ -309,12 +309,12 @@
 		gauss->GaussPoint(ig);
 
-		GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
+		GetNodalFunctions(&basis[0],gauss);
 		GetJacobianDeterminant(&Jdet, &xyz_list[0][0], gauss);
 
 		D_scalar=latentheat/heatcapacity*gauss->weight*Jdet;
 
-		TripleMultiply(&L[0],numdof,1,0,
+		TripleMultiply(&basis[0],numdof,1,0,
 					&D_scalar,1,1,0,
-					&L[0],1,numdof,0,
+					&basis[0],1,numdof,0,
 					&Ke->values[0],1);
 	}
@@ -353,5 +353,5 @@
 	IssmDouble v_gauss[2]={0.0};
 	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble L[NUMVERTICES];
+	IssmDouble basis[NUMVERTICES];
 	IssmDouble B[2][NUMVERTICES];
 	IssmDouble Bprime[2][NUMVERTICES];
@@ -388,5 +388,5 @@
 
 		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
-		GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
+		GetNodalFunctions(&basis[0],gauss);
 
 		vxaverage_input->GetInputValue(&vx,gauss);
@@ -397,7 +397,7 @@
 		DL_scalar=gauss->weight*Jdettria;
 
-		TripleMultiply( &L[0],1,numdof,1,
+		TripleMultiply(&basis[0],1,numdof,1,
 					&DL_scalar,1,1,0,
-					&L[0],1,numdof,0,
+					&basis[0],1,numdof,0,
 					&Ke->values[0],1);
 
@@ -468,5 +468,5 @@
 	IssmDouble xyz_list[NUMVERTICES][3];
 	IssmDouble Jdettria,dt,vx,vy;
-	IssmDouble L[NUMVERTICES];
+	IssmDouble basis[NUMVERTICES];
 	IssmDouble B[2][NUMVERTICES];
 	IssmDouble Bprime[2][NUMVERTICES];
@@ -504,11 +504,11 @@
 
 		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
-		GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
+		GetNodalFunctions(&basis[0],gauss);
 
 		DL_scalar=gauss->weight*Jdettria;
 
-		TripleMultiply( &L[0],1,numdof,1,
+		TripleMultiply(&basis[0],1,numdof,1,
 					&DL_scalar,1,1,0,
-					&L[0],1,numdof,0,
+					&basis[0],1,numdof,0,
 					&Ke->values[0],1);
 
@@ -542,5 +542,5 @@
 	IssmDouble  D,Jdet;
 	IssmDouble  xyz_list[NUMVERTICES][3];
-	IssmDouble  L[1][numdof];
+	IssmDouble  basis[numdof];
 	GaussTria  *gauss = NULL;
 
@@ -556,12 +556,11 @@
 		gauss->GaussPoint(ig);
 
+		GetNodalFunctions(&basis[0],gauss);
 		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
 		D=gauss->weight*Jdet;
 
-		GetL(&L[0][0], &xyz_list[0][0], gauss,NDOF1);
-
-		TripleMultiply(&L[0][0],1,3,1,
+		TripleMultiply(&basis[0],1,3,1,
 					&D,1,1,0,
-					&L[0][0],1,3,0,
+					&basis[0],1,3,0,
 					&Ke->values[0],1);
 	}
@@ -668,5 +667,5 @@
 	IssmDouble surface_mass_balance_g,basal_melting_g,basal_melting_correction_g,thickness_g;
 	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble L[NUMVERTICES];
+	IssmDouble basis[NUMVERTICES];
 	GaussTria* gauss=NULL;
 
@@ -690,5 +689,5 @@
 
 		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
-		GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
+		GetNodalFunctions(&basis[0],gauss);
 
 		surface_mass_balance_input->GetInputValue(&surface_mass_balance_g,gauss);
@@ -700,5 +699,5 @@
 		 basal_melting_correction_g=0.;
 
-		for(int i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(thickness_g+dt*(surface_mass_balance_g-basal_melting_g-basal_melting_correction_g))*L[i];
+		for(int i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(thickness_g+dt*(surface_mass_balance_g-basal_melting_g-basal_melting_correction_g))*basis[i];
 	}
 
@@ -718,5 +717,5 @@
 	IssmDouble surface_mass_balance_g,basal_melting_g,thickness_g;
 	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble L[NUMVERTICES];
+	IssmDouble basis[NUMVERTICES];
 	GaussTria* gauss=NULL;
 
@@ -738,5 +737,5 @@
 
 		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
-		GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
+		GetNodalFunctions(&basis[0],gauss);
 
 		surface_mass_balance_input->GetInputValue(&surface_mass_balance_g,gauss);
@@ -744,5 +743,5 @@
 		thickness_input->GetInputValue(&thickness_g,gauss);
 
-		for(int i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(thickness_g+dt*(surface_mass_balance_g-basal_melting_g))*L[i];
+		for(int i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(thickness_g+dt*(surface_mass_balance_g-basal_melting_g))*basis[i];
 	}
 
@@ -3401,5 +3400,5 @@
 		}
 
-		GetL(&B[0], &xyz_list[0][0], gauss,NDOF2);
+		GetBMacAyealFriction(&B[0], &xyz_list[0][0], gauss);
 		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
 		D_scalar=alpha2*gauss->weight*Jdet;
@@ -5947,5 +5946,5 @@
 	IssmDouble v_gauss[2]={0.0};
 	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble L[NUMVERTICES];
+	IssmDouble basis[NUMVERTICES];
 	IssmDouble B[2][NUMVERTICES];
 	IssmDouble Bprime[2][NUMVERTICES];
@@ -5980,5 +5979,5 @@
 
 		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
-		GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
+		GetNodalFunctions(&basis[0],gauss);
 
 		vx_input->GetInputValue(&vx,gauss);
@@ -5989,7 +5988,7 @@
 		DL_scalar=gauss->weight*Jdettria;
 
-		TripleMultiply( &L[0],1,numdof,1,
+		TripleMultiply(&basis[0],1,numdof,1,
 					&DL_scalar,1,1,0,
-					&L[0],1,numdof,0,
+					&basis[0],1,numdof,0,
 					&Ke->values[0],1);
 
@@ -6050,5 +6049,5 @@
 	IssmDouble  xyz_list[NUMVERTICES][3];
 	IssmDouble  B[2][numdof];
-	IssmDouble  L[NUMVERTICES];
+	IssmDouble  basis[NUMVERTICES];
 	IssmDouble  D[2][2];
 	GaussTria   *gauss = NULL;
@@ -6078,17 +6077,17 @@
 		GetBHydro(&B[0][0],&xyz_list[0][0],gauss); 
 		TripleMultiply(&B[0][0],2,numdof,1,
-									 &D[0][0],2,2,0,
-									 &B[0][0],2,numdof,0,
-									 &Ke->values[0],1);
+					&D[0][0],2,2,0,
+					&B[0][0],2,numdof,0,
+					&Ke->values[0],1);
 
 		/*Transient*/
 		if(reCast<bool,IssmDouble>(dt)){
-			GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
+			GetNodalFunctions(&basis[0],gauss);
 			D_scalar=sediment_storing*gauss->weight*Jdet;
 
-			TripleMultiply(&L[0],numdof,1,0,
-										 &D_scalar,1,1,0,
-										 &L[0],1,numdof,0,
-										 &Ke->values[0],1);
+			TripleMultiply(&basis[0],numdof,1,0,
+						&D_scalar,1,1,0,
+						&basis[0],1,numdof,0,
+						&Ke->values[0],1);
 		}
 	}
@@ -6110,5 +6109,5 @@
 	IssmDouble  xyz_list[NUMVERTICES][3];
 	IssmDouble  B[2][numdof];
-	IssmDouble  L[NUMVERTICES];
+	IssmDouble  basis[NUMVERTICES];
 	IssmDouble  D[2][2];
 	GaussTria   *gauss = NULL;
@@ -6149,11 +6148,11 @@
 		/*Transient*/
 		if(reCast<bool,IssmDouble>(dt)){
-			GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
+			GetNodalFunctions(&basis[0],gauss);
 			D_scalar=epl_storing*gauss->weight*Jdet;
 
-			TripleMultiply(&L[0],numdof,1,0,
-										 &D_scalar,1,1,0,
-										 &L[0],1,numdof,0,
-										 &Ke->values[0],1);
+			TripleMultiply(&basis[0],numdof,1,0,
+						&D_scalar,1,1,0,
+						&basis[0],1,numdof,0,
+						&Ke->values[0],1);
 		}
 	}
@@ -7058,5 +7057,5 @@
 	IssmDouble xyz_list[NUMVERTICES][3];
 	IssmDouble dhdt_g,basal_melting_g,surface_mass_balance_g,Jdettria;
-	IssmDouble L[NUMVERTICES];
+	IssmDouble basis[NUMVERTICES];
 	GaussTria* gauss=NULL;
 
@@ -7081,7 +7080,7 @@
 
 		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
-		GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
-
-		for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*L[i];
+		GetNodalFunctions(&basis[0],gauss);
+
+		for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*basis[i];
 	}
 
@@ -7101,5 +7100,5 @@
 	IssmDouble xyz_list[NUMVERTICES][3];
 	IssmDouble basal_melting_g,surface_mass_balance_g,dhdt_g,Jdettria;
-	IssmDouble L[NUMVERTICES];
+	IssmDouble basis[NUMVERTICES];
 	GaussTria* gauss=NULL;
 
@@ -7124,7 +7123,7 @@
 
 		GetJacobianDeterminant(&Jdettria, &xyz_list[0][0],gauss);
-		GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
-
-		for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*L[i];
+		GetNodalFunctions(&basis[0],gauss);
+
+		for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*basis[i];
 	}
 
Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 15312)
+++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 15313)
@@ -298,16 +298,13 @@
 /*}}}*/
 /*FUNCTION TriaRef::GetL{{{*/
-void TriaRef::GetL(IssmDouble* L, IssmDouble* xyz_list,GaussTria* gauss,int numdof){
-	/*Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof. 
-	 * For node i, Li can be expressed in the actual coordinate system
-	 * by: 
-	 *       numdof=1: 
-	 *                 Li=h;
-	 *       numdof=2:
-	 *                 Li=[ h   0 ]
-	 *                    [ 0   h ]
-	 * where h is the interpolation function for node i.
-	 *
-	 * We assume L has been allocated already, of size: NUMNODESP1 (numdof=1), or numdofx(numdof*NUMNODESP1) (numdof=2)
+void TriaRef::GetBMacAyealFriction(IssmDouble* B, IssmDouble* xyz_list,GaussTria* gauss){
+	/*Compute B  matrix. B=[B1 B2 B3] where Bi is square and of size 2. 
+	 * For node i, Bi can be expressed in the actual coordinate system
+	 * by: 
+	 *                 Bi=[ N   0 ]
+	 *                    [ 0   N ]
+	 * where N is the interpolation function for node i.
+	 *
+	 * We assume B has been allocated already, of size: 2 x (numdof*NUMNODESP1)
 	 */
 
@@ -319,16 +316,9 @@
 
 	/*Build L: */
-	if(numdof==1){
-		for (i=0;i<NUMNODESP1;i++){
-			L[i]=basis[i]; 
-		}
-	}
-	else{
-		for (i=0;i<NUMNODESP1;i++){
-			*(L+numdof*NUMNODESP1*0+numdof*i)=basis[i]; //L[0][NDOF2*i]=dbasis[0][i];
-			*(L+numdof*NUMNODESP1*0+numdof*i+1)=0;
-			*(L+numdof*NUMNODESP1*1+numdof*i)=0;
-			*(L+numdof*NUMNODESP1*1+numdof*i+1)=basis[i];
-		}
+	for (i=0;i<NUMNODESP1;i++){
+		*(B+2*NUMNODESP1*0+2*i)=basis[i]; //L[0][NDOF2*i]=dbasis[0][i];
+		*(B+2*NUMNODESP1*0+2*i+1)=0;
+		*(B+2*NUMNODESP1*1+2*i)=0;
+		*(B+2*NUMNODESP1*1+2*i+1)=basis[i];
 	}
 }
Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h	(revision 15312)
+++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h	(revision 15313)
@@ -30,5 +30,5 @@
 		void GetBPrognostic(IssmDouble* B_prog, IssmDouble* xyz_list, GaussTria* gauss);
 		void GetBHydro(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss);
-		void GetL(IssmDouble* L, IssmDouble* xyz_list,GaussTria* gauss,int numdof);
+		void GetBMacAyealFriction(IssmDouble* L, IssmDouble* xyz_list,GaussTria* gauss);
 		void GetJacobian(IssmDouble* J, IssmDouble* xyz_list,GaussTria* gauss);
 		void GetSegmentJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,GaussTria* gauss);