Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 8416)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 8417)
@@ -1371,15 +1371,75 @@
 ElementMatrix* Penta::CreateKMatrixDiagnosticPattynFriction(void){
 
+	/*Constants*/
+	const int numdof   = NDOF2*NUMVERTICES;
+	
+	/*Intermediaries */
+	int       i,j,ig;
+	int       analysis_type,drag_type;
+	double    xyz_list[NUMVERTICES][3];
+	double    xyz_list_tria[NUMVERTICES2D][3]={0.0};
+	double    slope_magnitude,alpha2,Jdet;
+	double    slope[3]={0.0,0.0,0.0};
+	double    MAXSLOPE=.06; // 6 %
+	double    MOUNTAINKEXPONENT=10;
+	double    L[2][numdof];
+	double    DL[2][2]={{ 0,0 },{0,0}}; //for basal drag
+	double    DL_scalar;
+	double    Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
+	Friction  *friction = NULL;
+	GaussPenta *gauss=NULL;
+
 	/*Initialize Element matrix and return if necessary*/
 	if(IsOnShelf() || !IsOnBed()) return NULL;
 
-	/*Build a tria element using the 3 nodes of the base of the penta. Then use 
-	 * the tria functionality to build a friction stiffness matrix on these 3
-	 * nodes: */
-	Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
-	ElementMatrix* Ke=tria->CreateKMatrixDiagnosticPattynFriction();
-	delete tria->matice; delete tria;
-
-	/*clean-up and return*/
+	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,PattynApproximationEnum);
+
+	/*Retrieve all inputs and parameters*/
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
+	for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<2;j++) xyz_list_tria[i][j]=xyz_list[i][j];
+	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+	inputs->GetParameterValue(&drag_type,DragTypeEnum);
+	Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input);
+	Input* vx_input=inputs->GetInput(VxEnum);           _assert_(vx_input);
+	Input* vy_input=inputs->GetInput(VyEnum);           _assert_(vy_input);
+	Input* vz_input=inputs->GetInput(VzEnum);           _assert_(vz_input);
+
+	/*build friction object, used later on: */
+	if (drag_type!=2)_error_(" non-viscous friction not supported yet!");
+	friction=new Friction("2d",inputs,matpar,analysis_type);
+
+	/* Start  looping on the number of gaussian points: */
+	gauss=new GaussPenta(0,1,2,2);
+	for (ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		GetTriaJacobianDeterminant(&Jdet, &xyz_list_tria[0][0],gauss);
+		GetL(&L[0][0], gauss,NDOF2);
+
+		surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
+		friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); 
+		slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2));
+
+		// If we have a slope > 6% for this element,  it means  we are on a mountain. In this particular case, 
+		//velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */
+		if (slope_magnitude>MAXSLOPE){
+			alpha2=pow((double)10,MOUNTAINKEXPONENT);
+		}
+		
+		DL_scalar=alpha2*gauss->weight*Jdet;
+		for (i=0;i<2;i++) DL[i][i]=DL_scalar;
+		
+		TripleMultiply( &L[0][0],2,numdof,1,
+					&DL[0][0],2,2,0,
+					&L[0][0],2,numdof,0,
+					&Ke_gg_gaussian[0][0],0);
+
+		for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_gaussian[i][j];
+	}
+
+	/*Clean up and return*/
+	delete gauss;
+	delete friction;
 	return Ke;
 }
Index: /issm/trunk/src/c/objects/Elements/PentaRef.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/PentaRef.cpp	(revision 8416)
+++ /issm/trunk/src/c/objects/Elements/PentaRef.cpp	(revision 8417)
@@ -519,4 +519,41 @@
 
 }
+/*}}}*/
+/*FUNCTION PentaRef::GetL{{{1*/
+void PentaRef::GetL(double* 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
+	 ** 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)
+	 **/
+
+	int i;
+	double l1l6[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];
+		}
+	}
+} 
 /*}}}*/
 /*FUNCTION PentaRef::GetLStokes {{{1*/
Index: /issm/trunk/src/c/objects/Elements/PentaRef.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/PentaRef.h	(revision 8416)
+++ /issm/trunk/src/c/objects/Elements/PentaRef.h	(revision 8417)
@@ -49,4 +49,5 @@
 		void GetBVert(double* B, double* xyz_list, GaussPenta* gauss);
 		void GetBprimeAdvec(double* Bprime_advec, double* xyz_list, GaussPenta* gauss);
+		void GetL(double* L, GaussPenta* gauss,int numdof);
 		void GetLStokes(double* LStokes, GaussPenta* gauss);
 		void GetLprimeStokes(double* LprimeStokes, double* xyz_list, GaussPenta* gauss);
Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 8416)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 8417)
@@ -988,78 +988,4 @@
 	for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdoftotal+(numdof+j)]+=Ke_gg[i][j];
 	for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[(i+numdof)*numdoftotal+j]+=Ke_gg[i][j];
-
-	/*Clean up and return*/
-	delete gauss;
-	delete friction;
-	return Ke;
-}	
-/*}}}*/
-/*FUNCTION Tria::CreateKMatrixDiagnosticPattynFriction {{{1*/
-ElementMatrix* Tria::CreateKMatrixDiagnosticPattynFriction(void){
-
-	/*Constants*/
-	const int numdof   = NDOF2*NUMVERTICES;
-	
-	/*Intermediaries */
-	int       i,j,ig;
-	int       analysis_type,drag_type;
-	double    xyz_list[NUMVERTICES][3];
-	double    slope_magnitude,alpha2,Jdet;
-	double    slope[2]={0.0,0.0};
-	double    MAXSLOPE=.06; // 6 %
-	double    MOUNTAINKEXPONENT=10;
-	double    L[2][numdof];
-	double    DL[2][2]={{ 0,0 },{0,0}}; //for basal drag
-	double    DL_scalar;
-	double    Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
-	Friction  *friction = NULL;
-	GaussTria *gauss=NULL;
-
-	/*Initialize Element matrix and return if necessary*/
-	if(IsOnShelf()) return NULL;
-	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,PattynApproximationEnum);
-
-	/*Retrieve all inputs and parameters*/
-	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
-	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
-	inputs->GetParameterValue(&drag_type,DragTypeEnum);
-	Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input);
-	Input* vx_input=inputs->GetInput(VxEnum);           _assert_(vx_input);
-	Input* vy_input=inputs->GetInput(VyEnum);           _assert_(vy_input);
-	Input* vz_input=inputs->GetInput(VzEnum);           _assert_(vz_input);
-
-	/*build friction object, used later on: */
-	if (drag_type!=2)_error_(" non-viscous friction not supported yet!");
-	friction=new Friction("2d",inputs,matpar,analysis_type);
-
-	/* Start  looping on the number of gaussian points: */
-	gauss=new GaussTria(2);
-	for (ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
-
-		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
-		GetL(&L[0][0], &xyz_list[0][0], gauss,NDOF2);
-
-		surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
-		friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); // TO UNCOMMENT
-		slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2));
-
-		// If we have a slope > 6% for this element,  it means  we are on a mountain. In this particular case, 
-		//velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */
-		if (slope_magnitude>MAXSLOPE){
-			alpha2=pow((double)10,MOUNTAINKEXPONENT);
-		}
-		
-		DL_scalar=alpha2*gauss->weight*Jdet;
-		for (i=0;i<2;i++) DL[i][i]=DL_scalar;
-		
-		TripleMultiply( &L[0][0],2,numdof,1,
-					&DL[0][0],2,2,0,
-					&L[0][0],2,numdof,0,
-					&Ke_gg_gaussian[0][0],0);
-
-		for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_gaussian[i][j];
-	}
 
 	/*Clean up and return*/
Index: /issm/trunk/src/c/objects/Elements/Tria.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.h	(revision 8416)
+++ /issm/trunk/src/c/objects/Elements/Tria.h	(revision 8417)
@@ -151,5 +151,4 @@
 		ElementMatrix* CreateKMatrixDiagnosticMacAyealFriction(void);
 		ElementMatrix* CreateKMatrixCouplingMacAyealPattynFriction(void);
-		ElementMatrix* CreateKMatrixDiagnosticPattynFriction(void);
 		ElementMatrix* CreateKMatrixDiagnosticHutter(void);
 		ElementMatrix* CreateKMatrixMelting(void);