Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5386)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5387)
@@ -2502,5 +2502,5 @@
 		 * grids: */
 
-		tria->CreateKMatrixDiagnosticHorizFriction(Kgg);
+		tria->CreateKMatrixCouplingMacAyealPattynFriction(Kgg);
 	}
 
@@ -2835,5 +2835,5 @@
 
 		tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
-		tria->CreateKMatrixDiagnosticHorizFriction(Kgg);
+		tria->CreateKMatrixDiagnosticPattynFriction(Kgg);
 		delete tria->matice; delete tria;
 	}
Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 5386)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 5387)
@@ -3591,5 +3591,5 @@
 	/*Do not forget to include friction: */
 	if(!shelf){
-		CreateKMatrixDiagnosticHorizFriction(Kgg);
+		CreateKMatrixDiagnosticMacAyealFriction(Kgg);
 	}
 
@@ -3601,6 +3601,6 @@
 }
 /*}}}*/
-/*FUNCTION Tria::CreateKMatrixDiagnosticHorizFriction {{{1*/
-void  Tria::CreateKMatrixDiagnosticHorizFriction(Mat Kgg){
+/*FUNCTION Tria::CreateKMatrixCouplingMacAyealPattynFriction {{{1*/
+void  Tria::CreateKMatrixCouplingMacAyealPattynFriction(Mat Kgg){
 
 	/* local declarations */
@@ -3612,5 +3612,6 @@
 	const int numdof   = 2 *numvertices;
 	double    xyz_list[numvertices][3];
-	int*      doflist=NULL;
+	int*      doflistm=NULL;
+	int*      doflistp=NULL;
 	int       numberofdofspernode=2;
 
@@ -3669,5 +3670,290 @@
 	/* Get node coordinates and dof list: */
 	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
-	GetDofList(&doflist);
+	GetDofList(&doflistm,MacAyealApproximationEnum);
+	GetDofList(&doflistp,PattynApproximationEnum);
+
+	if (shelf){
+		/*no friction, do nothing*/
+		return;
+	}
+
+	/*build friction object, used later on: */
+	if (drag_type!=2)ISSMERROR(" non-viscous friction not supported yet!");
+	friction=new Friction("2d",inputs,matpar,analysis_type);
+
+	/*COMPUT alpha2_list (TO BE DELETED)*/
+	for(i=0;i<numvertices;i++){
+		friction->GetAlpha2(&alpha2_list[i],&gauss[i][0],VxEnum,VyEnum,VzEnum);
+	}
+
+	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
+	GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
+
+	/* Start  looping on the number of gaussian points: */
+	for (ig=0; ig<num_gauss; ig++){
+		/*Pick up the gaussian point: */
+		gauss_weight=*(gauss_weights+ig);
+		gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 
+		gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
+		gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
+
+		/*Friction: */
+		// friction->GetAlpha2(&alpha2, gauss_l1l2l3,VxEnum,VyEnum,VzEnum); // TO UNCOMMENT
+		TriaRef::GetParameterValue(&alpha2,&alpha2_list[0],gauss_l1l2l3); // TO BE DELETED
+
+		// 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_l1l2l3[0]);
+		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_l1l2l3);
+
+		/*Get L matrix: */
+		GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode);
+
+		
+		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 (ig=0; ig<num_gauss; ig++)
+
+	/*Add Ke_gg to global matrix Kgg: */
+	MatSetValues(Kgg,numdof,doflistm,numdof,doflistp,(const double*)Ke_gg,ADD_VALUES);
+	MatSetValues(Kgg,numdof,doflistp,numdof,doflistm,(const double*)Ke_gg,ADD_VALUES);
+
+	xfree((void**)&first_gauss_area_coord);
+	xfree((void**)&second_gauss_area_coord);
+	xfree((void**)&third_gauss_area_coord);
+	xfree((void**)&gauss_weights);
+	xfree((void**)&doflistm);
+	xfree((void**)&doflistp);
+	delete friction;
+}	
+/*}}}*/
+/*FUNCTION Tria::CreateKMatrixDiagnosticMacAyealFriction {{{1*/
+void  Tria::CreateKMatrixDiagnosticMacAyealFriction(Mat Kgg){
+
+	/* local declarations */
+	int       i,j;
+	int       analysis_type;
+
+	/* node data: */
+	const int numvertices = 3;
+	const int numdof   = 2 *numvertices;
+	double    xyz_list[numvertices][3];
+	int*      doflist=NULL;
+	int       numberofdofspernode=2;
+
+	/* gaussian points: */
+	int     num_gauss,ig;
+	double *first_gauss_area_coord  = NULL;
+	double *second_gauss_area_coord = NULL;
+	double *third_gauss_area_coord  = NULL;
+	double *gauss_weights           = NULL;
+	double  gauss_weight;
+	double  gauss_l1l2l3[3];
+
+	/* matrices: */
+	double L[2][numdof];
+	double DL[2][2]={{ 0,0 },{0,0}}; //for basal drag
+	double DL_scalar;
+
+	/* local element matrices: */
+	double Ke_gg[numdof][numdof]={0.0};
+	double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
+	
+	double Jdet;
+	
+	/*slope: */
+	double  slope[2]={0.0,0.0};
+	double  slope_magnitude;
+
+	/*friction: */
+	Friction *friction = NULL;
+	double    alpha2;
+	double    alpha2_list[numvertices];                                       //TO BE DELETED
+	double    gauss[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED
+
+	double MAXSLOPE=.06; // 6 %
+	double MOUNTAINKEXPONENT=10;
+
+	/*inputs: */
+	bool shelf;
+	int  drag_type;
+	Input* surface_input=NULL;
+	Input* vx_input=NULL;
+	Input* vy_input=NULL;
+	Input* vz_input=NULL;
+
+	/*retrive parameters: */
+	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+
+	/*retrieve inputs :*/
+	inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
+	inputs->GetParameterValue(&drag_type,DragTypeEnum);
+	surface_input=inputs->GetInput(SurfaceEnum);
+	vx_input=inputs->GetInput(VxEnum);
+	vy_input=inputs->GetInput(VyEnum);
+	vz_input=inputs->GetInput(VzEnum);
+	
+	/* Get node coordinates and dof list: */
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetDofList(&doflist,MacAyealApproximationEnum);
+
+	if (shelf){
+		/*no friction, do nothing*/
+		return;
+	}
+
+	/*build friction object, used later on: */
+	if (drag_type!=2)ISSMERROR(" non-viscous friction not supported yet!");
+	friction=new Friction("2d",inputs,matpar,analysis_type);
+
+	/*COMPUT alpha2_list (TO BE DELETED)*/
+	for(i=0;i<numvertices;i++){
+		friction->GetAlpha2(&alpha2_list[i],&gauss[i][0],VxEnum,VyEnum,VzEnum);
+	}
+
+	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
+	GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
+
+	/* Start  looping on the number of gaussian points: */
+	for (ig=0; ig<num_gauss; ig++){
+		/*Pick up the gaussian point: */
+		gauss_weight=*(gauss_weights+ig);
+		gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 
+		gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
+		gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);
+
+		/*Friction: */
+		// friction->GetAlpha2(&alpha2, gauss_l1l2l3,VxEnum,VyEnum,VzEnum); // TO UNCOMMENT
+		TriaRef::GetParameterValue(&alpha2,&alpha2_list[0],gauss_l1l2l3); // TO BE DELETED
+
+		// 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_l1l2l3[0]);
+		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_l1l2l3);
+
+		/*Get L matrix: */
+		GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode);
+
+		
+		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 (ig=0; ig<num_gauss; ig++)
+
+	/*Add Ke_gg to global matrix Kgg: */
+	MatSetValues(Kgg,numdof,doflist,numdof,doflist,(const double*)Ke_gg,ADD_VALUES);
+
+	xfree((void**)&first_gauss_area_coord);
+	xfree((void**)&second_gauss_area_coord);
+	xfree((void**)&third_gauss_area_coord);
+	xfree((void**)&gauss_weights);
+	xfree((void**)&doflist);
+	delete friction;
+}	
+/*}}}*/
+/*FUNCTION Tria::CreateKMatrixDiagnosticPattynFriction {{{1*/
+void  Tria::CreateKMatrixDiagnosticPattynFriction(Mat Kgg){
+
+	/* local declarations */
+	int       i,j;
+	int       analysis_type;
+
+	/* node data: */
+	const int numvertices = 3;
+	const int numdof   = 2 *numvertices;
+	double    xyz_list[numvertices][3];
+	int*      doflist=NULL;
+	int       numberofdofspernode=2;
+
+	/* gaussian points: */
+	int     num_gauss,ig;
+	double *first_gauss_area_coord  = NULL;
+	double *second_gauss_area_coord = NULL;
+	double *third_gauss_area_coord  = NULL;
+	double *gauss_weights           = NULL;
+	double  gauss_weight;
+	double  gauss_l1l2l3[3];
+
+	/* matrices: */
+	double L[2][numdof];
+	double DL[2][2]={{ 0,0 },{0,0}}; //for basal drag
+	double DL_scalar;
+
+	/* local element matrices: */
+	double Ke_gg[numdof][numdof]={0.0};
+	double Ke_gg_gaussian[numdof][numdof]; //stiffness matrix contribution from drag
+	
+	double Jdet;
+	
+	/*slope: */
+	double  slope[2]={0.0,0.0};
+	double  slope_magnitude;
+
+	/*friction: */
+	Friction *friction = NULL;
+	double    alpha2;
+	double    alpha2_list[numvertices];                                       //TO BE DELETED
+	double    gauss[numvertices][numvertices] = {{1,0,0},{0,1,0},{0,0,1}}; //TO BE DELETED
+
+	double MAXSLOPE=.06; // 6 %
+	double MOUNTAINKEXPONENT=10;
+
+	/*inputs: */
+	bool shelf;
+	int  drag_type;
+	Input* surface_input=NULL;
+	Input* vx_input=NULL;
+	Input* vy_input=NULL;
+	Input* vz_input=NULL;
+
+	/*retrive parameters: */
+	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
+
+	/*retrieve inputs :*/
+	inputs->GetParameterValue(&shelf,ElementOnIceShelfEnum);
+	inputs->GetParameterValue(&drag_type,DragTypeEnum);
+	surface_input=inputs->GetInput(SurfaceEnum);
+	vx_input=inputs->GetInput(VxEnum);
+	vy_input=inputs->GetInput(VyEnum);
+	vz_input=inputs->GetInput(VzEnum);
+	
+	/* Get node coordinates and dof list: */
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+	GetDofList(&doflist,PattynApproximationEnum);
 
 	if (shelf){
Index: /issm/trunk/src/c/objects/Elements/Tria.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.h	(revision 5386)
+++ /issm/trunk/src/c/objects/Elements/Tria.h	(revision 5387)
@@ -124,5 +124,7 @@
 		void	  CreateKMatrixBalancedvelocities(Mat Kgg);
 		void	  CreateKMatrixDiagnosticMacAyeal(Mat Kgg);
-		void	  CreateKMatrixDiagnosticHorizFriction(Mat Kgg);
+		void	  CreateKMatrixCouplingMacAyealPattynFriction(Mat Kgg);
+		void	  CreateKMatrixDiagnosticMacAyealFriction(Mat Kgg);
+		void	  CreateKMatrixDiagnosticPattynFriction(Mat Kgg);
 		void	  CreateKMatrixDiagnosticHutter(Mat Kgg);
 		void	  CreateKMatrixDiagnosticSurfaceVert(Mat Kgg);