Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16871)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16872)
@@ -974,4 +974,85 @@
 	return Ke;
 }/*}}}*/
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAFriction(Element* element){/*{{{*/
+	/*Intermediaries*/
+	bool       mainlyfloating;
+	int        analysis_type,migration_style,point1;
+	IssmDouble alpha2,Jdet,fraction1,fraction2;
+	IssmDouble phi=1.0;
+	IssmDouble D_scalar,gllevelset;
+	IssmDouble *xyz_list = NULL;
+	Friction   *friction = NULL;
+	Gauss*     gauss     = NULL;
+
+	/*Return if element is inactive*/
+	if(element->IsFloating()) return 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(SSAApproximationEnum);
+	IssmDouble*    B      = xNew<IssmDouble>(2*numdof);
+	IssmDouble*    D      = xNewZeroInit<IssmDouble>(2*2);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+	element->FindParam(&migration_style,GroundinglineMigrationEnum);
+	Input* surface_input=element->GetInput(SurfaceEnum); _assert_(surface_input);
+	Input* vx_input=element->GetInput(VxEnum);           _assert_(vx_input);
+	Input* vy_input=element->GetInput(VyEnum);           _assert_(vy_input);
+	Input* gllevelset_input=NULL;
+	element->FindParam(&analysis_type,AnalysisTypeEnum);
+
+	/*build friction object, used later on: */
+	//friction=new Friction("2d",inputs,matpar,analysis_type);
+
+	/*Recover portion of element that is grounded*/
+	if(migration_style==SubelementMigrationEnum) phi=element->GetGroundedPortion(xyz_list);
+	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);
+	}
+	else{
+		gauss = element->NewGauss(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,VxEnum,VyEnum,VzEnum);
+		alpha2=1.0;
+		if(migration_style==SubelementMigrationEnum) alpha2=phi*alpha2;
+		if(migration_style==SubelementMigration2Enum){
+			gllevelset_input->GetInputValue(&gllevelset, gauss);
+			if(gllevelset<0) alpha2=0;
+		}
+
+		this->GetBSSAFriction(B, element, xyz_list, gauss);
+		element->JacobianDeterminant(&Jdet, xyz_list,gauss);
+		D_scalar=alpha2*gauss->weight*Jdet;
+		for(int i=0;i<2;i++) D[i*2+i]=D_scalar;
+
+		TripleMultiply(B,2,numdof,1,
+					D,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);
+	xDelete<IssmDouble>(D);
+	xDelete<IssmDouble>(B);
+	return Ke;
+}/*}}}*/
 ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAViscous(Element* element){/*{{{*/
 
@@ -1034,7 +1115,4 @@
 	xDelete<IssmDouble>(B);
 	return Ke;
-}/*}}}*/
-ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAFriction(Element* element){/*{{{*/
-	return NULL;
 }/*}}}*/
 ElementVector* StressbalanceAnalysis::CreatePVectorSSA(Element* element){/*{{{*/
@@ -1207,4 +1285,33 @@
 	/*Clean-up*/
 	xDelete<IssmDouble>(dbasis);
+}/*}}}*/
+void StressbalanceAnalysis::GetBSSAFriction(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::GetBSSAprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 16871)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 16872)
@@ -35,4 +35,5 @@
 		ElementVector* CreatePVectorSSAFront(Element* element);
 		void GetBSSA(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);
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 16871)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 16872)
@@ -106,4 +106,6 @@
 		virtual bool   IsOnBed()=0;
 		virtual bool   IsOnSurface()=0;
+		virtual void   GetGroundedPart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlyfloating)=0;
+		virtual IssmDouble GetGroundedPortion(IssmDouble* xyz_list)=0;
 		virtual void   GetInputListOnNodes(IssmDouble* pvalue,int enumtype)=0;
 		virtual void   GetInputListOnNodes(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue)=0;
@@ -144,4 +146,5 @@
       virtual Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order)=0;
       virtual Gauss* NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert)=0;
+      virtual Gauss* NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,bool mainlyfloating,int order)=0;
 		virtual Gauss* NewGaussBase(int order)=0;
 		virtual Gauss* NewGaussLine(int vertex1,int vertex2,int order)=0;
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 16871)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 16872)
@@ -92,4 +92,6 @@
 		void	 GetDofListVelocity(int** pdoflist,int setenum);
 		void	 GetDofListPressure(int** pdoflist,int setenum);
+		void   GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating);
+		IssmDouble GetGroundedPortion(IssmDouble* xyz_list);
 		int    GetNodeIndex(Node* node);
 		void   GetNodesSidList(int* sidlist);
@@ -228,6 +230,4 @@
 		void           GetVertexSidList(int* sidlist);
 		void           GetConnectivityList(int* connectivity);
-		void           GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating);
-		IssmDouble     GetGroundedPortion(IssmDouble* xyz_list);
 		int            GetElementType(void);
 		Input*         GetInput(int inputenum);
@@ -266,4 +266,5 @@
       Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order){_error_("not implemented yet");};
       Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert);
+      Gauss*         NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,bool mainlyfloating,int order){_error_("not implemented yet");};
 		Gauss*         NewGaussBase(int order);
 		Gauss*         NewGaussLine(int vertex1,int vertex2,int order);
Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 16871)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 16872)
@@ -142,4 +142,6 @@
 		void        ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
 		int         VelocityInterpolation(void){_error_("not implemented yet");};
+		void        GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating){_error_("not implemented yet");};
+		IssmDouble  GetGroundedPortion(IssmDouble* xyz_list){_error_("not implemented yet");};
 		Input*      GetInput(int inputenum);
 		Input*      GetMaterialInput(int inputenum){_error_("not implemented yet");};
@@ -157,4 +159,5 @@
       Gauss*      NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order){_error_("not implemented yet");};
       Gauss*      NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert){_error_("not implemented yet");};
+      Gauss*      NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,bool mainlyfloating,int order){_error_("not implemented yet");};
 		Gauss*      NewGaussBase(int order){_error_("not implemented yet");};
 		Gauss*      NewGaussLine(int vertex1,int vertex2,int order){_error_("not implemented yet");};
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16871)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16872)
@@ -2298,4 +2298,10 @@
 
 	return new GaussTria(area_coordinates,order);
+}
+/*}}}*/
+/*FUNCTION Tria::NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,bool mainlyfloating){{{*/
+Gauss* Tria::NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,bool mainlyfloating,int order){
+
+	return new GaussTria(point1,fraction1,fraction2,mainlyfloating,order);
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16871)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16872)
@@ -91,4 +91,6 @@
 		void	      GetDofListVelocity(int** pdoflist,int setenum);
 		void	      GetDofListPressure(int** pdoflist,int setenum);
+		void        GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating);
+		IssmDouble  GetGroundedPortion(IssmDouble* xyz_list);
 		int         GetNodeIndex(Node* node);
 		void        GetNodesSidList(int* sidlist);
@@ -264,6 +266,4 @@
 		void           GetVertexSidList(int* sidlist);
 		void           GetConnectivityList(int* connectivity);
-		void           GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating);
-		IssmDouble     GetGroundedPortion(IssmDouble* xyz_list);
 		IssmDouble     GetYcoord(Gauss* gauss);
 		IssmDouble     GetZcoord(Gauss* gauss){_error_("not implemented");};
@@ -298,4 +298,5 @@
 		Gauss*         NewGauss(int order);
       Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order);
+      Gauss*         NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,bool mainlyfloating,int order);
       Gauss*         NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert);
 		Gauss*         NewGaussBase(int order);