Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16908)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16909)
@@ -3090,6 +3090,80 @@
 }/*}}}*/
 ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAHOFriction(Element* element){/*{{{*/
-	return NULL;
-
+
+	if(element->IsFloating() || !element->IsOnBed()) return NULL;
+
+	/*Constants*/
+	int numnodes    = element->GetNumberOfNodes();
+	int numdof      = 2*numnodes;
+	int numdoftotal = 4*numnodes;
+
+	/*Intermediaries */
+	int         i,j;
+	IssmDouble  Jdet2d,alpha2;
+	IssmDouble *xyz_list_tria=NULL;
+	IssmDouble* L           = xNewZeroInit<IssmDouble>(2*numdof);
+	IssmDouble  DL[2][2]     ={{ 0,0 },{0,0}}; //for basal drag
+	IssmDouble  DL_scalar;
+	IssmDouble* Ke_gg       = xNewZeroInit<IssmDouble>(numdof*numdof);
+	Node       *node_list[2*numnodes];
+	int*        cs_list= xNew<int>(2*numnodes);
+
+	/*Initialize Element matrix and return if necessary*/
+	ElementMatrix* Ke1=element->NewElementMatrix(SSAApproximationEnum);
+	ElementMatrix* Ke2=element->NewElementMatrix(HOApproximationEnum);
+	ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
+	delete Ke1; delete Ke2;
+
+	/*Prepare node list*/
+	for(i=0;i<numnodes;i++){
+		node_list[i+0*numnodes] = element->GetNode(i);
+		node_list[i+1*numnodes] = element->GetNode(i);
+		cs_list[i+0*numnodes] = XYEnum;
+		cs_list[i+1*numnodes] = XYEnum;
+	}
+
+	/*retrieve inputs :*/
+	element->GetVerticesCoordinatesBase(&xyz_list_tria);
+	Input* vx_input=element->GetInput(VxEnum);  _assert_(vx_input);
+	Input* vy_input=element->GetInput(VyEnum);  _assert_(vy_input);
+	Input* vz_input=element->GetInput(VzEnum);  _assert_(vz_input);
+
+	/*build friction object, used later on: */
+	Friction* friction=new Friction(element,2);
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=element->NewGaussBase(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		/*Friction: */
+		friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input);
+		element->JacobianDeterminantBase(&Jdet2d, xyz_list_tria,gauss);
+		this->GetBHOFriction(L,element,xyz_list_tria,gauss);
+
+		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,2,numdof,1,
+					&DL[0][0],2,2,0,
+					L,2,numdof,0,
+					Ke_gg,1);
+	}
+
+	for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdoftotal+(numdof+j)]+=Ke_gg[i*numdof+j];
+	for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[(i+numdof)*numdoftotal+j]+=Ke_gg[i*numdof+j];
+
+	/*Transform Coordinate System*/
+	element->TransformStiffnessMatrixCoord(Ke,node_list,2*numnodes,cs_list);
+
+	/*Clean up and return*/
+	delete gauss;
+	delete friction;
+	xDelete<int>(cs_list);
+	xDelete<IssmDouble>(xyz_list_tria);
+	xDelete<IssmDouble>(Ke_gg);
+	return Ke;
 }/*}}}*/
 ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAHOViscous(Element* element){/*{{{*/
@@ -3110,5 +3184,4 @@
 	IssmDouble  D_scalar;
 	IssmDouble* Ke_gg          = xNewZeroInit<IssmDouble>(numdofp*numdofm);
-	IssmDouble* Ke_gg_gaussian = xNew<IssmDouble>(numdofp*numdofm);
 	Node        *node_list[2*numnodes];
 	int*         cs_list= xNew<int>(2*numnodes);
@@ -3161,7 +3234,5 @@
 					&D[0][0],3,3,0,
 					Bprime,3,numdofm,0,
-					Ke_gg_gaussian,0);
-
-		for( i=0; i<numdofp; i++) for(j=0;j<numdofm; j++) Ke_gg[i*numdofm+j]+=Ke_gg_gaussian[i*numdofm+j];
+					Ke_gg,1);
 	} 
 	for(i=0;i<numdofp;i++) for(j=0;j<numdofm;j++) Ke->values[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i*numdofm+j];
@@ -3177,5 +3248,5 @@
 	delete gauss_tria;
 	xDelete<IssmDouble>(Ke_gg);
-	xDelete<IssmDouble>(Ke_gg_gaussian);
+	xDelete<IssmDouble>(xyz_list);
 	xDelete<int>(cs_list);
 	return Ke;