Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16881)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16882)
@@ -980,5 +980,5 @@
 
 	/*Intermediaries*/
-	bool       mainlyfloating;
+	bool        mainlyfloating;
 	int         migration_style,point1;
 	IssmDouble  alpha2,Jdet,fraction1,fraction2;
@@ -1755,5 +1755,74 @@
 	if(element->IsFloating() || !element->IsOnBed()) return NULL;
 
-	return NULL;
+	/*Intermediaries*/
+	bool        mainlyfloating;
+	int         migration_style,point1;
+	IssmDouble  alpha2,Jdet,fraction1,fraction2;
+	IssmDouble  gllevelset,phi=1.;
+	IssmDouble *xyz_list_base = NULL;
+	Gauss*      gauss         = NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+	int numdof   = numnodes*2;
+
+	/*Initialize Element matrix and vectors*/
+	ElementMatrix* Ke      = element->NewElementMatrix(HOApproximationEnum);
+	IssmDouble*    B       = xNew<IssmDouble>(2*numdof);
+	IssmDouble     D[2][2] = {0.};
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinatesBase(&xyz_list_base);
+	element->FindParam(&migration_style,GroundinglineMigrationEnum);
+	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);
+	Input* gllevelset_input = NULL;
+
+	/*build friction object, used later on: */
+	Friction* friction=new Friction(element,2);
+
+	/*Recover portion of element that is grounded*/
+	if(migration_style==SubelementMigrationEnum) phi=element->GetGroundedPortion(xyz_list_base);
+	if(migration_style==SubelementMigration2Enum){
+		gllevelset_input=element->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
+		element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
+		//gauss = element->NewGauss(point1,fraction1,fraction2,mainlyfloating,2);
+		gauss=element->NewGaussBase(2);
+	}
+	else{
+		gauss=element->NewGaussBase(2);
+	}
+
+	/* Start  looping on the number of gaussian points: */
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input);
+		if(migration_style==SubelementMigrationEnum) alpha2=phi*alpha2;
+		if(migration_style==SubelementMigration2Enum){
+			gllevelset_input->GetInputValue(&gllevelset, gauss);
+			if(gllevelset<0.) alpha2=0.;
+		}
+
+		this->GetBHOFriction(B,element,xyz_list_base,gauss);
+		element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
+		for(int i=0;i<2;i++) D[i][i]=alpha2*gauss->weight*Jdet;
+
+		TripleMultiply(B,2,numdof,1,
+					&D[0][0],2,2,0,
+					B,2,numdof,0,
+					&Ke->values[0],1);
+	}
+
+	/*Transform Coordinate System*/
+	element->TransformStiffnessMatrixCoord(Ke,XYEnum);
+
+	/*Clean up and return*/
+	delete gauss;
+	delete friction;
+	xDelete<IssmDouble>(xyz_list_base);
+	xDelete<IssmDouble>(B);
+	return Ke;
 }/*}}}*/
 ElementVector* StressbalanceAnalysis::CreatePVectorHO(Element* element){/*{{{*/
@@ -1950,4 +2019,33 @@
 	xDelete<IssmDouble>(dbasis);
 }/*}}}*/
+void StressbalanceAnalysis::GetBHOFriction(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* 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 finiteelement function for node i.
+	 *
+	 * We assume B has been allocated already, of size: 2 x (numdof*numnodes)
+	 */
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Get nodal functions derivatives*/
+	IssmDouble* basis=xNew<IssmDouble>(numnodes);
+	element->NodalFunctions(basis,gauss);
+
+	/*Build L: */
+	for(int i=0;i<numnodes;i++){
+		B[2*numnodes*0+2*i+0] = basis[i];
+		B[2*numnodes*0+2*i+1] = 0.;
+		B[2*numnodes*1+2*i+0] = 0.;
+		B[2*numnodes*1+2*i+1] = basis[i];
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(basis);
+}/*}}}*/
 void StressbalanceAnalysis::InputUpdateFromSolutionHO(IssmDouble* solution,Element* element){/*{{{*/
 
@@ -2108,5 +2206,93 @@
 }/*}}}*/
 ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSFriction(Element* element){/*{{{*/
-	return NULL;
+
+	if(element->IsFloating() || !element->IsOnBed()) return NULL;
+
+	/*If on water or not FS, skip stiffness: */
+	int approximation;
+	element->GetInputValue(&approximation,ApproximationEnum);
+	if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
+
+	/*Intermediaries*/
+	bool        mainlyfloating;
+	int         i,meshtype,dim,epssize;
+	int         migration_style,point1;
+	IssmDouble  alpha2,Jdet,fraction1,fraction2;
+	IssmDouble  gllevelset,phi=1.;
+	IssmDouble *xyz_list_base = NULL;
+	Gauss*      gauss         = NULL;
+
+	/*Get problem dimension*/
+	element->FindParam(&meshtype,MeshTypeEnum);
+	switch(meshtype){
+		case Mesh2DverticalEnum: dim = 2;break;
+		case Mesh3DEnum:         dim = 3;break;
+		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+	}
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int vnumnodes = element->GetNumberOfNodesVelocity();
+	int pnumnodes = element->GetNumberOfNodesPressure();
+	int numdof    = vnumnodes*dim + pnumnodes;
+
+	/*Initialize Element matrix and vectors*/
+	ElementMatrix* Ke      = element->NewElementMatrix(FSvelocityEnum);
+	IssmDouble*    B       = xNew<IssmDouble>((dim-1)*numdof);
+	IssmDouble*    D       = xNewZeroInit<IssmDouble>((dim-1)*(dim-1));
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinatesBase(&xyz_list_base);
+	element->FindParam(&migration_style,GroundinglineMigrationEnum);
+	Input* vx_input         = element->GetInput(VxEnum);      _assert_(vx_input);
+	Input* vy_input         = element->GetInput(VyEnum);      _assert_(vy_input);
+	Input* vz_input         = NULL;
+	if(dim==3){    vz_input = element->GetInput(VzEnum);      _assert_(vz_input);}
+	Input* gllevelset_input = NULL;
+
+	/*build friction object, used later on: */
+	Friction* friction=new Friction(element,dim==3?3:1);
+
+	/*Recover portion of element that is grounded*/
+	if(migration_style==SubelementMigrationEnum) phi=element->GetGroundedPortion(xyz_list_base);
+	if(migration_style==SubelementMigration2Enum){
+		gllevelset_input=element->GetInput(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
+		element->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
+		//gauss = element->NewGauss(point1,fraction1,fraction2,mainlyfloating,2);
+		gauss=element->NewGaussBase(3);
+	}
+	else{
+		gauss=element->NewGaussBase(3);
+	}
+
+	/* Start  looping on the number of gaussian points: */
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		friction->GetAlpha2(&alpha2,gauss,vx_input,vy_input,vz_input);
+		if(migration_style==SubelementMigrationEnum) alpha2=phi*alpha2;
+		if(migration_style==SubelementMigration2Enum){
+			gllevelset_input->GetInputValue(&gllevelset, gauss);
+			if(gllevelset<0.) alpha2=0.;
+		}
+
+		this->GetBFSFriction(B,element,dim,xyz_list_base,gauss);
+		element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
+		for(int i=0;i<2;i++) D[i*(dim-1)+i] = alpha2*gauss->weight*Jdet; //taub_x = -alpha2 v_x (same for y)
+
+		TripleMultiply(B,dim-1,numdof,1,
+					D,dim-1,dim-1,0,
+					B,dim-1,numdof,0,
+					&Ke->values[0],1);
+	}
+
+	/*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/
+
+	/*Clean up and return*/
+	delete gauss;
+	delete friction;
+	xDelete<IssmDouble>(xyz_list_base);
+	xDelete<IssmDouble>(B);
+	xDelete<IssmDouble>(D);
+	return Ke;
 }/*}}}*/
 ElementVector* StressbalanceAnalysis::CreatePVectorFS(Element* element){/*{{{*/
@@ -2575,4 +2761,55 @@
 	xDelete<IssmDouble>(pbasis);
 }/*}}}*/
+void StressbalanceAnalysis::GetBFSFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+	/* 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 in 3d 
+	 *       Li=[ h 0 0 0 ]
+	 *	 	      [ 0 h 0 0 ]
+	 *	in 2d:
+	 *       Li=[ h 0 0 ]
+	 * where h is the interpolation function for node i.
+	 */
+
+	/*Fetch number of nodes for this finite element*/
+	int pnumnodes = element->GetNumberOfNodesPressure();
+	int vnumnodes = element->GetNumberOfNodesVelocity();
+	int pnumdof   = pnumnodes;
+	int vnumdof   = vnumnodes*dim;
+
+	/*Get nodal functions derivatives*/
+	IssmDouble* vbasis=xNew<IssmDouble>(vnumnodes);
+	element->NodalFunctionsVelocity(vbasis,gauss);
+
+	/*Build B: */
+	if(dim==3){
+		for(int i=0;i<vnumnodes;i++){
+			B[(vnumdof+pnumdof)*0+3*i+0] = vbasis[i];
+			B[(vnumdof+pnumdof)*0+3*i+1] = 0.;
+			B[(vnumdof+pnumdof)*0+3*i+2] = 0.;
+
+			B[(vnumdof+pnumdof)*1+3*i+0] = 0.;
+			B[(vnumdof+pnumdof)*1+3*i+1] = vbasis[i];
+			B[(vnumdof+pnumdof)*1+3*i+2] = 0.;
+		}
+		for(int i=0;i<pnumnodes;i++){
+			B[(vnumdof+pnumdof)*0+i+vnumdof+0] = 0.;
+			B[(vnumdof+pnumdof)*1+i+vnumdof+0] = 0.;
+		}
+	}
+	else{
+		for(int i=0;i<vnumnodes;i++){
+			B[2*i+0] = vbasis[i];
+			B[2*i+1] = 0.;
+		}
+
+		for(int i=0;i<pnumnodes;i++){
+			B[i+vnumdof+0] = 0.;
+		}
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(vbasis);
+}/*}}}*/
 void StressbalanceAnalysis::GetSolutionFromInputsFS(Vector<IssmDouble>* solution,Element* element){/*{{{*/
 
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 16881)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 16882)
@@ -35,6 +35,6 @@
 		ElementVector* CreatePVectorSSAFront(Element* element);
 		void GetBSSA(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
+		void GetBSSAprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
 		void GetBSSAFriction(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
-		void GetBSSAprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
 		void InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element);
 		/*L1L2*/
@@ -52,4 +52,5 @@
 		void GetBHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
 		void GetBHOprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
+		void GetBHOFriction(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
 		void InputUpdateFromSolutionHO(IssmDouble* solution,Element* element);
 		/*FS*/
@@ -63,4 +64,5 @@
 		void GetBFS(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
 		void GetBFSprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
+		void GetBFSFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
 		void GetSolutionFromInputsFS(Vector<IssmDouble>* solution,Element* element);
 		void InputUpdateFromSolutionFS(IssmDouble* solution,Element* element);