Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16950)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16951)
@@ -3086,4 +3086,16 @@
 
 	return NULL;
+	/*Initialize Element matrix and return if necessary*/
+	if(element->IsFloating() || !element->IsOnBed()) return NULL;
+
+	/*Build a tria element using the 3 nodes of the base of the penta. Then use 
+	 * the tria functionality to build a friction stiffness matrix on these 3
+	 * nodes: */
+	Element* basalelement = element->SpawnBasalElement();
+	ElementMatrix* Ke=CreateKMatrixSSAFriction(element);
+	basalelement->DeleteMaterials(); delete basalelement;
+
+	/*clean-up and return*/
+	return Ke;
 }/*}}}*/
 ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSA3dViscous(Element* element){/*{{{*/
@@ -3476,9 +3488,127 @@
 }/*}}}*/
 ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAFSFriction(Element* element){/*{{{*/
-	return NULL;
-
+
+	/*Constants*/
+	const int numdofs   = (6+1)*3 + 6*1;
+	const int numdofm   = 6 *2;
+	const int numdof2d  = 3 *3;
+	const int numdof2dm = 3 *2;
+	const int numdoftot = 6*2 + (6+1)*3 +6; // HO + FS vel + FS Pressure
+
+	/*Intermediaries */
+	int        i,j,approximation;
+	int        dim=3;
+	IssmDouble FSreconditioning,viscosity,alpha2_gauss,Jdet2d;
+	IssmDouble bed_normal[3];
+	IssmDouble LSSAFS[8][numdof2dm];
+	IssmDouble LprimeSSAFS[8][numdofs];
+	IssmDouble DLSSAFS[8][8]={0.0};
+	IssmDouble LFSSSA[4][numdof2d];
+	IssmDouble LprimeFSSSA[4][numdof2dm];
+	IssmDouble DLFSSSA[4][4]={0.0};
+	IssmDouble Ke_drag_gaussian[numdof2dm][numdofs];
+	IssmDouble Ke_drag_gaussian2[numdof2d][numdof2dm];
+	IssmDouble *xyz_list      = NULL;
+	IssmDouble *xyz_list_tria = NULL;
+
+	/*If on water or not FS, skip stiffness: */
+	element->GetInputValue(&approximation,ApproximationEnum);
+	if(element->IsFloating() || !element->IsOnBed()) return NULL;
+
+	int vnumnodes = element->GetNumberOfNodesVelocity();
+	int pnumnodes = element->GetNumberOfNodesPressure();
+	int numnodes  = 2*vnumnodes-1+pnumnodes;
+
+	/*Prepare node list*/
+	int* cs_list = xNew<int>(2*vnumnodes+pnumnodes);
+	Node **node_list = xNew<Node*>(2*vnumnodes+pnumnodes);
+	for(i=0;i<vnumnodes-1;i++){
+		node_list[i] = element->GetNode(i);
+		cs_list[i]   = XYEnum;
+	}
+	for(i=0;i<vnumnodes;i++){
+		node_list[i+vnumnodes-1] = element->GetNode(i);
+		cs_list[i+vnumnodes-1]   = XYZEnum;
+	}
+	for(i=0;i<pnumnodes;i++){
+		node_list[2*vnumnodes-1+i] = element->GetNode(vnumnodes+i);
+		cs_list[2*vnumnodes-1+i]   = PressureEnum;
+	}
+
+	ElementMatrix* Ke1=element->NewElementMatrix(SSAApproximationEnum);
+	ElementMatrix* Ke2=element->NewElementMatrix(FSvelocityEnum);
+	ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
+	delete Ke1; delete Ke2;
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+	element->GetVerticesCoordinatesBase(&xyz_list_tria);
+	element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
+	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,3);
+
+	/* 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);
+
+		element->JacobianDeterminantBase(&Jdet2d,xyz_list_tria,gauss);
+		this->GetLSSAFS(&LSSAFS[0][0], element,gauss);
+		this->GetLprimeSSAFS(&LprimeSSAFS[0][0], element,xyz_list, gauss);
+		this->GetLFSSSA(&LFSSSA[0][0],element, gauss);
+		this->GetLprimeFSSSA(&LprimeFSSSA[0][0], element,xyz_list, gauss);
+
+		element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
+
+		element->NormalBase(&bed_normal[0],xyz_list_tria);
+		friction->GetAlpha2(&alpha2_gauss, gauss,vx_input,vy_input,vz_input);
+
+		DLSSAFS[0][0]=alpha2_gauss*gauss->weight*Jdet2d;
+		DLSSAFS[1][1]=alpha2_gauss*gauss->weight*Jdet2d;
+		DLSSAFS[2][2]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[0]*bed_normal[2];
+		DLSSAFS[3][3]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[1]*bed_normal[2];
+		DLSSAFS[4][4]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[0];
+		DLSSAFS[5][5]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[1];
+		DLSSAFS[6][6]=FSreconditioning*gauss->weight*Jdet2d*bed_normal[0];
+		DLSSAFS[7][7]=FSreconditioning*gauss->weight*Jdet2d*bed_normal[1];
+
+		DLFSSSA[0][0]=alpha2_gauss*gauss->weight*Jdet2d;
+		DLFSSSA[1][1]=alpha2_gauss*gauss->weight*Jdet2d;
+		DLFSSSA[2][2]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[0]*bed_normal[2];
+		DLFSSSA[3][3]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[1]*bed_normal[2];
+
+		TripleMultiply( &LSSAFS[0][0],8,numdof2dm,1,
+					&DLSSAFS[0][0],8,8,0,
+					&LprimeSSAFS[0][0],8,numdofs,0,
+					&Ke_drag_gaussian[0][0],0);
+
+		TripleMultiply( &LFSSSA[0][0],4,numdof2d,1,
+					&DLFSSSA[0][0],4,4,0,
+					&LprimeFSSSA[0][0],4,numdof2dm,0,
+					&Ke_drag_gaussian2[0][0],0);
+		for(i=0;i<numdof2dm;i++) for(j=0;j<numdofs;j++) Ke->values[i*numdoftot+j+numdofm]+=Ke_drag_gaussian[i][j];
+		for(i=0;i<numdof2d;i++) for(j=0;j<numdof2dm;j++) Ke->values[(i+numdofm)*numdoftot+j]+=Ke_drag_gaussian2[i][j];
+	}
+
+	/*Transform Coordinate System*/
+	element->TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list);
+
+	/*Clean up and return*/
+	xDelete<int>(cs_list);
+	xDelete<Node*>(node_list);
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(xyz_list_tria);
+	delete gauss;
+	delete friction;
+	return Ke;
 }/*}}}*/
 ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAFSViscous(Element* element){/*{{{*/
 
+	printf("CouplingSSAFSviscous\n");
 	/*Constants*/
 	const int numdofm     = 2 *3;
@@ -4168,4 +4298,230 @@
 	xDelete<IssmDouble>(dbasis);
 }/*}}}*/
+void StressbalanceAnalysis::GetLSSAFS(IssmDouble* LFS,Element* element,Gauss* gauss_in){/*{{{*/
+	/*
+	 * 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: 
+	 *       Li=[ h    0 ]
+	 *	 	      [ 0    h ]
+	 *	 	      [ h    0 ]
+	 *	 	      [ 0    h ]
+	 *	 	      [ h    0 ]
+	 *	 	      [ 0    h ]
+	 *	 	      [ h    0 ]
+	 *	 	      [ 0    h ]
+	 * where h is the interpolation function for node i.
+	 */
+
+	int num_dof=2;
+	IssmDouble L1L2l3[3];
+
+	/*Cast gauss to GaussPenta*/
+	_assert_(gauss_in->Enum()==GaussPentaEnum);
+	GaussPenta* gauss = dynamic_cast<GaussPenta*>(gauss_in);
+
+	/*Get L1L2l3 in actual coordinate system: */
+	L1L2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
+	L1L2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
+	L1L2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
+
+	/*Build LFS: */
+	for(int i=0;i<3;i++){
+		LFS[num_dof*3*0+num_dof*i+0] = L1L2l3[i];
+		LFS[num_dof*3*0+num_dof*i+1] = 0;
+		LFS[num_dof*3*1+num_dof*i+0] = 0;
+		LFS[num_dof*3*1+num_dof*i+1] = L1L2l3[i];
+		LFS[num_dof*3*2+num_dof*i+0] = L1L2l3[i];
+		LFS[num_dof*3*2+num_dof*i+1] = 0;
+		LFS[num_dof*3*3+num_dof*i+0] = 0;
+		LFS[num_dof*3*3+num_dof*i+1] = L1L2l3[i];
+		LFS[num_dof*3*4+num_dof*i+0] = L1L2l3[i];
+		LFS[num_dof*3*4+num_dof*i+1] = 0;
+		LFS[num_dof*3*5+num_dof*i+0] = 0;
+		LFS[num_dof*3*5+num_dof*i+1] = L1L2l3[i];
+		LFS[num_dof*3*6+num_dof*i+0] = L1L2l3[i];
+		LFS[num_dof*3*6+num_dof*i+1] = 0;
+		LFS[num_dof*3*7+num_dof*i+0] = 0;
+		LFS[num_dof*3*7+num_dof*i+1] = L1L2l3[i];
+	}
+}/*}}}*/
+void StressbalanceAnalysis::GetLprimeSSAFS(IssmDouble* LprimeFS,Element* element,IssmDouble* xyz_list,Gauss* gauss_in){/*{{{*/
+	/* Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. 
+	 * For node i, Lpi can be expressed in the actual coordinate system
+	 * by: 
+	 *       Lpi=[ h    0    0   0]
+	 *		       [ 0    h    0   0]
+	 *		       [ 0    0    h   0]
+	 *		       [ 0    0    h   0]
+	 *		       [ 0    0  dh/dz 0]
+	 *		       [ 0    0  dh/dz 0]
+	 *           [ 0    0    0   h]
+	 *           [ 0    0    0   h]
+	 * where h is the interpolation function for node i.
+	 */
+	int num_dof=3;
+	int num_dof_vel=3*7;
+	int num_dof_total=3*7+1*6;
+	IssmDouble L1L2l3[3];
+	IssmDouble dbasis[3][6];
+
+	/*Cast gauss to GaussPenta*/
+	_assert_(gauss_in->Enum()==GaussPentaEnum);
+	GaussPenta* gauss = dynamic_cast<GaussPenta*>(gauss_in);
+
+	/*Get L1L2l3 in actual coordinate system: */
+	L1L2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
+	L1L2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
+	L1L2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
+
+	element->NodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);
+
+	/*Build LprimeFS: */
+	for(int i=0;i<3;i++){
+		LprimeFS[num_dof_total*0+num_dof*i+0] = L1L2l3[i];
+		LprimeFS[num_dof_total*0+num_dof*i+1] = 0.;
+		LprimeFS[num_dof_total*0+num_dof*i+2] = 0.;
+		LprimeFS[num_dof_total*1+num_dof*i+0] = 0.;
+		LprimeFS[num_dof_total*1+num_dof*i+1] = L1L2l3[i];
+		LprimeFS[num_dof_total*1+num_dof*i+2] = 0.;
+		LprimeFS[num_dof_total*2+num_dof*i+0] = 0.;
+		LprimeFS[num_dof_total*2+num_dof*i+1] = 0.;
+		LprimeFS[num_dof_total*2+num_dof*i+2] = L1L2l3[i];
+		LprimeFS[num_dof_total*3+num_dof*i+0] = 0.;
+		LprimeFS[num_dof_total*3+num_dof*i+1] = 0.;
+		LprimeFS[num_dof_total*3+num_dof*i+2] = L1L2l3[i];
+		LprimeFS[num_dof_total*4+num_dof*i+0] = 0.;
+		LprimeFS[num_dof_total*4+num_dof*i+1] = 0.;
+		LprimeFS[num_dof_total*4+num_dof*i+2] = dbasis[2][i];
+		LprimeFS[num_dof_total*5+num_dof*i+0] = 0.;
+		LprimeFS[num_dof_total*5+num_dof*i+1] = 0.;
+		LprimeFS[num_dof_total*5+num_dof*i+2] = dbasis[2][i];
+		LprimeFS[num_dof_total*6+num_dof*i+0] = 0.;
+		LprimeFS[num_dof_total*6+num_dof*i+1] = 0.;
+		LprimeFS[num_dof_total*6+num_dof*i+2] = 0.;
+		LprimeFS[num_dof_total*7+num_dof*i+0] = 0.;
+		LprimeFS[num_dof_total*7+num_dof*i+1] = 0.;
+		LprimeFS[num_dof_total*7+num_dof*i+2] = 0.;
+	}
+	for(int i=3;i<7;i++){
+		LprimeFS[num_dof_total*0+num_dof*i+0] = 0.;
+		LprimeFS[num_dof_total*0+num_dof*i+1] = 0.;
+		LprimeFS[num_dof_total*0+num_dof*i+2] = 0.;
+		LprimeFS[num_dof_total*1+num_dof*i+0] = 0.;
+		LprimeFS[num_dof_total*1+num_dof*i+1] = 0.;
+		LprimeFS[num_dof_total*1+num_dof*i+2] = 0.;
+		LprimeFS[num_dof_total*2+num_dof*i+0] = 0.;
+		LprimeFS[num_dof_total*2+num_dof*i+1] = 0.;
+		LprimeFS[num_dof_total*2+num_dof*i+2] = 0.;
+		LprimeFS[num_dof_total*3+num_dof*i+0] = 0.;
+		LprimeFS[num_dof_total*3+num_dof*i+1] = 0.;
+		LprimeFS[num_dof_total*3+num_dof*i+2] = 0.;
+		LprimeFS[num_dof_total*4+num_dof*i+0] = 0.;
+		LprimeFS[num_dof_total*4+num_dof*i+1] = 0.;
+		LprimeFS[num_dof_total*4+num_dof*i+2] = 0.;
+		LprimeFS[num_dof_total*5+num_dof*i+0] = 0.;
+		LprimeFS[num_dof_total*5+num_dof*i+1] = 0.;
+		LprimeFS[num_dof_total*5+num_dof*i+2] = 0.;
+		LprimeFS[num_dof_total*6+num_dof*i+0] = 0.;
+		LprimeFS[num_dof_total*6+num_dof*i+1] = 0.;
+		LprimeFS[num_dof_total*6+num_dof*i+2] = 0.;
+		LprimeFS[num_dof_total*7+num_dof*i+0] = 0.;
+		LprimeFS[num_dof_total*7+num_dof*i+1] = 0.;
+		LprimeFS[num_dof_total*7+num_dof*i+2] = 0.;
+	}
+	for(int i=0;i<3;i++){
+		LprimeFS[num_dof_total*0+num_dof_vel+i] = 0.;
+		LprimeFS[num_dof_total*1+num_dof_vel+i] = 0.;
+		LprimeFS[num_dof_total*2+num_dof_vel+i] = 0.;
+		LprimeFS[num_dof_total*3+num_dof_vel+i] = 0.;
+		LprimeFS[num_dof_total*4+num_dof_vel+i] = 0.;
+		LprimeFS[num_dof_total*5+num_dof_vel+i] = 0.;
+		LprimeFS[num_dof_total*6+num_dof_vel+i] = L1L2l3[i];
+		LprimeFS[num_dof_total*7+num_dof_vel+i] = L1L2l3[i];
+	}
+	for(int i=3;i<6;i++){
+		LprimeFS[num_dof_total*0+num_dof_vel+i] = 0.;
+		LprimeFS[num_dof_total*1+num_dof_vel+i] = 0.;
+		LprimeFS[num_dof_total*2+num_dof_vel+i] = 0.;
+		LprimeFS[num_dof_total*3+num_dof_vel+i] = 0.;
+		LprimeFS[num_dof_total*4+num_dof_vel+i] = 0.;
+		LprimeFS[num_dof_total*5+num_dof_vel+i] = 0.;
+		LprimeFS[num_dof_total*6+num_dof_vel+i] = 0.;
+		LprimeFS[num_dof_total*7+num_dof_vel+i] = 0.;
+	}
+}/*}}}*/
+void StressbalanceAnalysis::GetLFSSSA(IssmDouble* LFS,Element* element,Gauss* gauss_in){/*{{{*/
+	/* 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: 
+	 *       Li=[ h    0    0 ]
+	 *	 	      [ 0    h    0 ]
+	 *		      [ 0    0    h ]
+	 *		      [ 0    0    h ]
+	 * where h is the interpolation function for node i.
+	 */
+
+	int num_dof=3;
+	IssmDouble L1L2l3[3];
+
+	/*Cast gauss to GaussPenta*/
+	_assert_(gauss_in->Enum()==GaussPentaEnum);
+	GaussPenta* gauss = dynamic_cast<GaussPenta*>(gauss_in);
+
+	/*Get L1L2l3 in actual coordinate system: */
+	L1L2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
+	L1L2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
+	L1L2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
+
+	/*Build LFS: */
+	for(int i=0;i<3;i++){
+		LFS[num_dof*3*0+num_dof*i+0] = L1L2l3[i];
+		LFS[num_dof*3*0+num_dof*i+1] = 0.;
+		LFS[num_dof*3*0+num_dof*i+2] = 0.;
+		LFS[num_dof*3*1+num_dof*i+0] = 0.;
+		LFS[num_dof*3*1+num_dof*i+1] = L1L2l3[i];
+		LFS[num_dof*3*1+num_dof*i+2] = 0.;
+		LFS[num_dof*3*2+num_dof*i+0] = 0.;
+		LFS[num_dof*3*2+num_dof*i+1] = 0.;
+		LFS[num_dof*3*2+num_dof*i+2] = L1L2l3[i];
+		LFS[num_dof*3*3+num_dof*i+0] = 0.;
+		LFS[num_dof*3*3+num_dof*i+1] = 0.;
+		LFS[num_dof*3*3+num_dof*i+2] = L1L2l3[i];
+	}
+}/*}}}*/
+void StressbalanceAnalysis::GetLprimeFSSSA(IssmDouble* LprimeFS,Element* element,IssmDouble* xyz_list,Gauss* gauss_in){/*{{{*/
+	/* Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. 
+	 * For node i, Lpi can be expressed in the actual coordinate system
+	 * by: 
+	 *       Lpi=[ h    0 ]
+	 *		       [ 0    h ]
+	 *		       [ h    0 ]
+	 *		       [ 0    h ]
+	 * where h is the interpolation function for node i.
+	 */
+	int num_dof=2;
+	IssmDouble L1L2l3[3];
+
+	/*Cast gauss to GaussPenta*/
+	_assert_(gauss_in->Enum()==GaussPentaEnum);
+	GaussPenta* gauss = dynamic_cast<GaussPenta*>(gauss_in);
+
+	/*Get L1L2l3 in actual coordinate system: */
+	L1L2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
+	L1L2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
+	L1L2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
+
+	/*Build LprimeFS: */
+	for(int i=0;i<3;i++){
+		LprimeFS[num_dof*3*0+num_dof*i+0] = L1L2l3[i];
+		LprimeFS[num_dof*3*0+num_dof*i+1] = 0.;
+		LprimeFS[num_dof*3*1+num_dof*i+0] = 0.;
+		LprimeFS[num_dof*3*1+num_dof*i+1] = L1L2l3[i];
+		LprimeFS[num_dof*3*2+num_dof*i+0] = L1L2l3[i];
+		LprimeFS[num_dof*3*2+num_dof*i+1] = 0.;
+		LprimeFS[num_dof*3*3+num_dof*i+0] = 0.;
+		LprimeFS[num_dof*3*3+num_dof*i+1] = L1L2l3[i];
+	}
+}/*}}}*/
 void StressbalanceAnalysis::InputUpdateFromSolutionHOFS(IssmDouble* solution,Element* element){/*{{{*/
 
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 16950)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 16951)
@@ -98,4 +98,8 @@
 		void GetBSSAFSTria(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
 		void GetBprimeSSAFSTria(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss);
+		void GetLFSSSA(IssmDouble* L,Element* element,Gauss* gauss);
+		void GetLSSAFS(IssmDouble* L,Element* element,Gauss* gauss);
+		void GetLprimeFSSSA(IssmDouble* Lprime,Element* element,IssmDouble* xyz_list,Gauss* gauss);
+		void GetLprimeSSAFS(IssmDouble* Lprime,Element* element,IssmDouble* xyz_list,Gauss* gauss);
 		void InputUpdateFromSolutionHOFS(IssmDouble* solution,Element* element);
 		void InputUpdateFromSolutionSSAFS(IssmDouble* solution,Element* element);