Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 8417)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 8418)
@@ -777,11 +777,84 @@
 ElementMatrix* Penta::CreateKMatrixCouplingMacAyealPattynFriction(void){
 
+	/*Constants*/
+	const int numdof        = NDOF2 *NUMVERTICES;
+	const int numdoftotal   = NDOF4 *NUMVERTICES;
+	
+	/*Intermediaries */
+	int       i,j,ig,analysis_type,drag_type;
+	double    Jdet2d,slope_magnitude,alpha2;
+	double    xyz_list[NUMVERTICES][3];
+	double    xyz_list_tria[NUMVERTICES2D][3]={0.0};
+	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[numdof][numdof]     ={0.0};
+	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;
-
-	Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
-	ElementMatrix* Ke=tria->CreateKMatrixCouplingMacAyealPattynFriction();
-	delete tria->matice; delete tria;
-
+	ElementMatrix* Ke1=new ElementMatrix(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
+	ElementMatrix* Ke2=new ElementMatrix(nodes,NUMVERTICES,this->parameters,PattynApproximationEnum);
+	ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
+	delete Ke1; delete Ke2;
+
+	/*retrieve inputs :*/
+	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);
+
+		/*Friction: */
+		friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);
+
+		// 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: */
+		surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
+		slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2));
+
+		if (slope_magnitude>MAXSLOPE){
+			alpha2=pow((double)10,MOUNTAINKEXPONENT);
+		}
+
+		GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
+		GetL(&L[0][0], gauss,NDOF2);
+
+		DL_scalar=alpha2*gauss->weight*Jdet2d;
+		for (i=0;i<2;i++) DL[i][i]=DL_scalar; 
+		
+		/*  Do the triple producte tL*D*L: */
+		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_gg[i][j]+=Ke_gg_gaussian[i][j];
+	}
+
+	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;
 }
Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 8417)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 8418)
@@ -904,94 +904,4 @@
 	return Ke;
 }
-/*}}}*/
-/*FUNCTION Tria::CreateKMatrixCouplingMacAyealPattynFriction {{{1*/
-ElementMatrix* Tria::CreateKMatrixCouplingMacAyealPattynFriction(void){
-
-	/*Constants*/
-	const int numdof        = NDOF2 *NUMVERTICES;
-	const int numdoftotal   = NDOF4 *NUMVERTICES;
-	
-	/*Intermediaries */
-	int       i,j,ig,analysis_type,drag_type;
-	double    Jdet,slope_magnitude,alpha2;
-	double    xyz_list[NUMVERTICES][3];
-	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[numdof][numdof]     ={0.0};
-	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* Ke1=new ElementMatrix(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
-	ElementMatrix* Ke2=new ElementMatrix(nodes,NUMVERTICES,this->parameters,PattynApproximationEnum);
-	ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
-	delete Ke1; delete Ke2;
-
-	/*retrieve inputs :*/
-	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);
-
-		/*Friction: */
-		friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);
-
-		// 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: */
-		surface_input->GetParameterDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
-		slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2));
-
-		if (slope_magnitude>MAXSLOPE){
-			alpha2=pow((double)10,MOUNTAINKEXPONENT);
-		}
-
-		/* Get Jacobian determinant: */
-		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
-
-		/*Get L matrix: */
-		GetL(&L[0][0], &xyz_list[0][0], gauss,NDOF2);
-
-		
-		DL_scalar=alpha2*gauss->weight*Jdet;
-		for (i=0;i<2;i++){
-			DL[i][i]=DL_scalar;
-		}
-		
-		/*  Do the triple producte tL*D*L: */
-		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_gg[i][j]+=Ke_gg_gaussian[i][j];
-	}
-
-	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::CreateKMatrixDiagnosticHutter{{{1*/
Index: /issm/trunk/src/c/objects/Elements/Tria.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.h	(revision 8417)
+++ /issm/trunk/src/c/objects/Elements/Tria.h	(revision 8418)
@@ -150,5 +150,4 @@
 		ElementMatrix* CreateKMatrixDiagnosticMacAyealViscous(void);
 		ElementMatrix* CreateKMatrixDiagnosticMacAyealFriction(void);
-		ElementMatrix* CreateKMatrixCouplingMacAyealPattynFriction(void);
 		ElementMatrix* CreateKMatrixDiagnosticHutter(void);
 		ElementMatrix* CreateKMatrixMelting(void);