Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16349)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16350)
@@ -10658,5 +10658,5 @@
 	IssmDouble* vz       = xNew<IssmDouble>(vnumnodes);
 	IssmDouble* vel      = xNew<IssmDouble>(vnumnodes);
-	IssmDouble* pressure = xNew<IssmDouble>(vnumnodes);
+	IssmDouble* pressure = xNew<IssmDouble>(pnumnodes);
 
 	/*Get dof list: */
Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 16349)
+++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 16350)
@@ -334,7 +334,7 @@
 	/*Compute B  matrix. B=[Bv1 Bv2 ... Bp1 Bp2 ...] where Bvi is of size 3*NDOF3. 
 	 * For node i, Bvi can be expressed in the actual coordinate system
-	 * by: 	   Bvi=[ dh/dx          0              0      ]
+	 * by: 	   Bvi=[ dh/dx          0             0      ]
 	 *					[   0           dh/dy           0      ]
-	 *					[   0             0           dh/dy    ]
+	 *					[   0             0           dh/dz    ]
 	 *					[ 1/2*dh/dy    1/2*dh/dx        0      ]
 	 *					[ 1/2*dh/dz       0         1/2*dh/dx  ]
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16349)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16350)
@@ -20,5 +20,7 @@
 
 /*Element macros*/
-#define NUMVERTICES 3
+#define NUMVERTICES   3
+#define NUMVERTICES1D 2
+
 /*Constructors/destructor/copy*/
 /*FUNCTION Tria::Tria(){{{*/
@@ -219,5 +221,16 @@
 		#ifdef _HAVE_STRESSBALANCE_
 		case StressbalanceAnalysisEnum:
-			return CreateKMatrixStressbalanceSSA();
+			int approximation;
+			inputs->GetInputValue(&approximation,ApproximationEnum);
+			switch(approximation){
+				case SSAApproximationEnum:
+					return CreateKMatrixStressbalanceSSA();
+				case FSApproximationEnum:
+					return CreateKMatrixStressbalanceFS();
+				case NoneApproximationEnum:
+					return NULL;
+				default:
+					_error_("Approximation " << EnumToStringx(approximation) << " not supported yet");
+			}
 			break;
 		case StressbalanceSIAAnalysisEnum:
@@ -373,5 +386,16 @@
 		#ifdef _HAVE_STRESSBALANCE_
 		case StressbalanceAnalysisEnum:
-			return CreatePVectorStressbalanceSSA();
+			int approximation;
+			inputs->GetInputValue(&approximation,ApproximationEnum);
+			switch(approximation){
+				case SSAApproximationEnum:
+					return CreatePVectorStressbalanceSSA();
+				case FSApproximationEnum:
+					return CreatePVectorStressbalanceFS();
+				case NoneApproximationEnum:
+					return NULL;
+				default:
+					_error_("Approximation " << EnumToStringx(approximation) << " not supported yet");
+			}
 			break;
 		case StressbalanceSIAAnalysisEnum:
@@ -1679,5 +1703,15 @@
 		#ifdef _HAVE_STRESSBALANCE_
 		case StressbalanceAnalysisEnum:
-			InputUpdateFromSolutionStressbalanceHoriz(solution);
+			int approximation;
+			inputs->GetInputValue(&approximation,ApproximationEnum);
+			if(approximation==FSApproximationEnum || approximation==NoneApproximationEnum){
+				InputUpdateFromSolutionStressbalanceFS(solution);
+			}
+			else if (approximation==SSAApproximationEnum || approximation==SIAApproximationEnum){
+				InputUpdateFromSolutionStressbalanceHoriz(solution);
+			}
+			else{
+				_error_("approximation not supported yet");
+			}
 			break;
 		case StressbalanceSIAAnalysisEnum:
@@ -2056,4 +2090,48 @@
 }
 /*}}}*/
+/*FUNCTION Tria::HasEdgeOnBed {{{*/
+bool Tria::HasEdgeOnBed(){
+
+	IssmDouble values[NUMVERTICES];
+	IssmDouble sum;
+
+	/*Retrieve all inputs and parameters*/
+	GetInputListOnVertices(&values[0],MeshVertexonbedEnum);
+	sum = values[0]+values[1]+values[2];
+
+	_assert_(sum==0. || sum==1. || sum==2.);
+
+	if(sum>1.){
+		return true;
+	}
+	else{
+		return false;
+	}
+}
+/*}}}*/
+/*FUNCTION Tria::EdgeOnBedIndices{{{*/
+void Tria::EdgeOnBedIndices(int* pindex1,int* pindex2){
+
+	bool found=false;
+	IssmDouble values[NUMVERTICES];
+
+	/*Retrieve all inputs and parameters*/
+	GetInputListOnVertices(&values[0],MeshVertexonbedEnum);
+
+	for(int i=0;i<NUMVERTICES;i++){
+		if(values[i]==1.){
+			if(found){
+				*pindex2 = i;
+				return;
+			}
+			else{
+				*pindex1 = i;
+			}
+		}
+	}
+
+	_error_("Could not find 2 vertices on bed");
+}
+/*}}}*/
 /*FUNCTION Tria::IsFloating {{{*/
 bool   Tria::IsFloating(){
@@ -3037,4 +3115,110 @@
 
 #ifdef _HAVE_STRESSBALANCE_
+/*FUNCTION Tria::CreateKMatrixStressbalanceFS{{{*/
+ElementMatrix* Tria::CreateKMatrixStressbalanceFS(void){
+
+	ElementMatrix* Ke1 = NULL;
+	ElementMatrix* Ke2 = NULL;
+	ElementMatrix* Ke  = NULL;
+
+	/*compute all stiffness matrices for this element*/
+	Ke1=CreateKMatrixStressbalanceFSViscous();
+	Ke2=CreateKMatrixStressbalanceFSFriction();
+	Ke =new ElementMatrix(Ke1,Ke2);
+
+	/*clean-up and return*/
+	delete Ke1;
+	delete Ke2;
+	return Ke;
+
+}
+/*}}}*/
+/*FUNCTION Tria::CreateKMatrixStressbalanceFSViscous {{{*/
+ElementMatrix* Tria::CreateKMatrixStressbalanceFSViscous(void){
+
+	/*Intermediaries */
+	int        i,approximation;
+	IssmDouble Jdet,viscosity,FSreconditioning,D_scalar;
+	IssmDouble xyz_list[NUMVERTICES][3];
+	IssmDouble epsilon[3]; /* epsilon=[exx,eyy,exy];*/
+	GaussTria *gauss=NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int vnumnodes = this->NumberofNodesVelocity();
+	int pnumnodes = this->NumberofNodesPressure();
+	int numdof    = vnumnodes*NDOF2 + pnumnodes*NDOF1;
+
+	/*Prepare coordinate system list*/
+	int* cs_list = xNew<int>(vnumnodes+pnumnodes);
+	for(i=0;i<vnumnodes;i++) cs_list[i]           = XYEnum;
+	for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
+
+	/*Initialize Element matrix and vectors*/
+	ElementMatrix* Ke     = new ElementMatrix(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
+	IssmDouble*    B      = xNew<IssmDouble>(5*numdof);
+	IssmDouble*    Bprime = xNew<IssmDouble>(5*numdof);
+	IssmDouble*    D      = xNewZeroInit<IssmDouble>(5*5);
+
+	/*Retrieve all inputs and parameters*/
+	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+	parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
+	Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
+	Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
+
+	/* Start  looping on the number of gaussian points: */
+	gauss=new GaussTria(5);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
+		GetBFS(B,&xyz_list[0][0],gauss); 
+		GetBprimeFS(Bprime,&xyz_list[0][0],gauss); 
+
+		this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
+		material->GetViscosity2dvertical(&viscosity,&epsilon[0]);
+
+		D_scalar=gauss->weight*Jdet;
+		for(i=0;i<3;i++) D[i*5+i] = +D_scalar*2.*viscosity;
+		for(i=3;i<5;i++) D[i*5+i] = -D_scalar*FSreconditioning;
+
+		TripleMultiply(B,5,numdof,1,
+					D,5,5,0,
+					Bprime,5,numdof,0,
+					Ke->values,1);
+	}
+
+	/*Transform Coordinate System*/
+	TransformStiffnessMatrixCoord(Ke,nodes,(vnumnodes+pnumnodes),cs_list);
+
+	/*Clean up and return*/
+	delete gauss;
+	xDelete<int>(cs_list);
+	xDelete<IssmDouble>(B);
+	xDelete<IssmDouble>(Bprime);
+	xDelete<IssmDouble>(D);
+	return Ke;
+}
+/*}}}*/
+/*FUNCTION Tria::CreateKMatrixStressbalanceFSFriction{{{*/
+ElementMatrix* Tria::CreateKMatrixStressbalanceFSFriction(void){
+
+	/*Intermediaries */
+	int         i,j;
+	int         analysis_type,approximation;
+	IssmDouble  alpha2,Jdet2d;
+	IssmDouble  FSreconditioning,viscosity;
+	IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
+	IssmDouble  xyz_list[NUMVERTICES][3];
+	IssmDouble  xyz_list_seg[NUMVERTICES1D][3];
+	Friction   *friction = NULL;
+	GaussTria  *gauss    = NULL;
+
+	/*Initialize Element matrix and return if necessary*/
+	if(IsFloating() || !HasEdgeOnBed()) return NULL;
+
+	_error_("STOP");
+}
+/*}}}*/
 /*FUNCTION Tria::CreateKMatrixStressbalanceSSA {{{*/
 ElementMatrix* Tria::CreateKMatrixStressbalanceSSA(void){
@@ -3226,4 +3410,181 @@
 }
 /*}}}*/
+/*FUNCTION Tria::CreatePVectorStressbalanceFS {{{*/
+ElementVector* Tria::CreatePVectorStressbalanceFS(void){
+
+	ElementVector* pe1;
+	ElementVector* pe2;
+	ElementVector* pe3;
+	ElementVector* pe;
+
+	/*compute all stiffness matrices for this element*/
+	pe1=CreatePVectorStressbalanceFSViscous();
+	pe2=CreatePVectorStressbalanceFSShelf();
+	pe3=CreatePVectorStressbalanceFSFront();
+	pe =new ElementVector(pe1,pe2,pe3);
+
+	/*clean-up and return*/
+	delete pe1;
+	delete pe2;
+	delete pe3;
+	return pe;
+}
+/*}}}*/
+/*FUNCTION Tria::CreatePVectorStressbalanceFSFront{{{*/
+ElementVector* Tria::CreatePVectorStressbalanceFSFront(void){
+
+	/*Intermediaries */
+	int         i;
+	IssmDouble  ls[NUMVERTICES];
+	IssmDouble  xyz_list[NUMVERTICES][3];
+	bool        isfront;
+
+	/*Retrieve all inputs and parameters*/
+	GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
+
+	/*If the level set is awlays <=0, there is no ice front here*/
+	isfront = false;
+	if(ls[0]>0. || ls[1]>0. || ls[2]>0.){
+		if(ls[0]*ls[1]<0. || ls[0]*ls[2]<0. || (ls[0]*ls[1]+ls[0]*ls[2]+ls[1]*ls[2]==0.)){
+			isfront = true;
+		}
+	}
+
+	/*If no front, return NULL*/
+	if(!isfront) return NULL;
+
+	_error_("STOP");
+}
+/*}}}*/
+/*FUNCTION Tria::CreatePVectorStressbalanceFSViscous {{{*/
+ElementVector* Tria::CreatePVectorStressbalanceFSViscous(void){
+
+	/*Intermediaries*/
+	int        i;
+	int        approximation;
+	IssmDouble Jdet,gravity,rho_ice;
+	IssmDouble forcex,forcey;
+	IssmDouble xyz_list[NUMVERTICES][3];
+	GaussTria *gauss=NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int vnumnodes = this->NumberofNodesVelocity();
+	int pnumnodes = this->NumberofNodesPressure();
+
+	/*Prepare coordinate system list*/
+	int* cs_list = xNew<int>(vnumnodes+pnumnodes);
+	for(i=0;i<vnumnodes;i++) cs_list[i]           = XYEnum;
+	for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
+
+	/*Initialize Element matrix and vectors*/
+	ElementVector* pe     = new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
+	IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
+
+	/*Retrieve all inputs and parameters*/
+	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+	Input* loadingforcex_input=inputs->GetInput(LoadingforceXEnum);  _assert_(loadingforcex_input);
+	Input* loadingforcey_input=inputs->GetInput(LoadingforceYEnum);  _assert_(loadingforcey_input);
+	rho_ice = matpar->GetRhoIce();
+	gravity = matpar->GetG();
+
+	/* Start  looping on the number of gaussian points: */
+	gauss=new GaussTria(5);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
+		GetNodalFunctionsVelocity(vbasis, gauss);
+
+		loadingforcex_input->GetInputValue(&forcex,gauss);
+		loadingforcey_input->GetInputValue(&forcey,gauss);
+
+		for(i=0;i<vnumnodes;i++){
+			pe->values[i*NDOF2+1] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i];
+			pe->values[i*NDOF2+0] += +rho_ice*forcex *Jdet*gauss->weight*vbasis[i];
+			pe->values[i*NDOF2+1] += +rho_ice*forcey *Jdet*gauss->weight*vbasis[i];
+		}
+	}
+
+	/*Transform coordinate system*/
+	TransformLoadVectorCoord(pe,nodes,vnumnodes+pnumnodes,cs_list);
+
+	/*Clean up and return*/
+	delete gauss;
+	xDelete<int>(cs_list);
+	xDelete<IssmDouble>(vbasis);
+	return pe;
+}
+/*}}}*/
+/*FUNCTION Tria::CreatePVectorStressbalanceFSShelf{{{*/
+ElementVector* Tria::CreatePVectorStressbalanceFSShelf(void){
+
+	/*Intermediaries*/
+	int         i,j;
+	int         indices[2];
+	IssmDouble  gravity,rho_water,bed,water_pressure;
+	IssmDouble  normal_vel,vx,vy,dt;
+	IssmDouble	xyz_list_seg[NUMVERTICES1D][3];
+	IssmDouble  xyz_list[NUMVERTICES][3];
+	IssmDouble	normal[2];
+	IssmDouble  Jdet;
+
+	/*Initialize Element vector and return if necessary*/
+	if(!HasEdgeOnBed() || !IsFloating()) return NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int vnumnodes = this->NumberofNodesVelocity();
+	int pnumnodes = this->NumberofNodesPressure();
+
+	/*Prepare coordinate system list*/
+	int* cs_list = xNew<int>(vnumnodes+pnumnodes);
+	for(i=0;i<vnumnodes;i++) cs_list[i]           = XYEnum;
+	for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
+
+	/*Initialize Element matrix and vectors*/
+	ElementVector* pe     = new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
+	IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
+
+	/*Retrieve all inputs and parameters*/
+	rho_water=matpar->GetRhoWater();
+	gravity=matpar->GetG();
+	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+	Input* bed_input=inputs->GetInput(BedEnum); _assert_(bed_input);
+
+	/*Get vertex indices that lie on bed*/
+	this->EdgeOnBedIndices(&indices[0],&indices[1]);
+
+	for(i=0;i<NUMVERTICES1D;i++) for(j=0;j<2;j++) xyz_list_seg[i][j]=xyz_list[indices[i]][j];
+
+	/* Start looping on the number of gauss 1d (nodes on the bedrock) */
+	GaussTria* gauss=new GaussTria(indices[0],indices[1],2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		GetSegmentJacobianDeterminant(&Jdet,&xyz_list_seg[0][0],gauss);
+		GetNodalFunctionsVelocity(vbasis, gauss);
+
+		GetSegmentNormal(&normal[0],xyz_list_seg);
+		_assert_(normal[1]<0.);
+		bed_input->GetInputValue(&bed, gauss);
+		water_pressure=gravity*rho_water*bed;
+
+		for(i=0;i<vnumnodes;i++){
+			pe->values[2*i+0]+=(water_pressure)*gauss->weight*Jdet*vbasis[i]*normal[0];
+			pe->values[2*i+1]+=(water_pressure)*gauss->weight*Jdet*vbasis[i]*normal[1];
+		}
+	}
+
+	/*Transform coordinate system*/
+	TransformLoadVectorCoord(pe,nodes,(vnumnodes+pnumnodes),cs_list);
+
+	/*Clean up and return*/
+	delete gauss;
+	xDelete<int>(cs_list);
+	xDelete<IssmDouble>(vbasis);
+	return pe;
+}
+/*}}}*/
 /*FUNCTION Tria::CreatePVectorStressbalanceSSA {{{*/
 ElementVector* Tria::CreatePVectorStressbalanceSSA(){
@@ -3694,4 +4055,77 @@
 	xDelete<int>(doflist);
 
+}
+/*}}}*/
+/*FUNCTION Tria::InputUpdateFromSolutionStressbalanceFS {{{*/
+void  Tria::InputUpdateFromSolutionStressbalanceFS(IssmDouble* solution){
+
+	int          i;
+	int*         vdoflist=NULL;
+	int*         pdoflist=NULL;
+	IssmDouble   FSreconditioning;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int vnumnodes = this->NumberofNodesVelocity();
+	int pnumnodes = this->NumberofNodesPressure();
+	int vnumdof   = vnumnodes*NDOF2;
+	int pnumdof   = pnumnodes*NDOF1;
+
+	/*Initialize values*/
+	IssmDouble* vvalues  = xNew<IssmDouble>(vnumdof);
+	IssmDouble* pvalues  = xNew<IssmDouble>(pnumdof);
+	IssmDouble* vx       = xNew<IssmDouble>(vnumnodes);
+	IssmDouble* vy       = xNew<IssmDouble>(vnumnodes);
+	IssmDouble* vel      = xNew<IssmDouble>(vnumnodes);
+	IssmDouble* pressure = xNew<IssmDouble>(pnumnodes);
+
+	/*Get dof list: */
+	GetDofListVelocity(&vdoflist,GsetEnum);
+	GetDofListPressure(&pdoflist,GsetEnum);
+
+	/*Use the dof list to index into the solution vector: */
+	for(i=0;i<vnumdof;i++) vvalues[i]=solution[vdoflist[i]];
+	for(i=0;i<pnumdof;i++) pvalues[i]=solution[pdoflist[i]];
+
+	/*Transform solution in Cartesian Space*/
+	TransformSolutionCoord(&vvalues[0],nodes,vnumnodes,XYEnum);
+
+	/*Ok, we have vx and vy in values, fill in all arrays: */
+	for(i=0;i<vnumnodes;i++){
+		vx[i] = vvalues[i*NDOF2+0];
+		vy[i] = vvalues[i*NDOF2+1];
+		if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
+	}
+	for(i=0;i<pnumnodes;i++){
+		pressure[i] = pvalues[i];
+		if(xIsNan<IssmDouble>(pressure[i])) _error_("NaN found in solution vector");
+	}
+
+	/*Recondition pressure and compute vel: */
+	this->parameters->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
+	for(i=0;i<pnumnodes;i++) pressure[i]=pressure[i]*FSreconditioning;
+	for(i=0;i<vnumnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i]);
+
+	/*Now, we have to move the previous inputs  to old 
+	 * status, otherwise, we'll wipe them off: */
+	this->inputs->ChangeEnum(VxEnum,VxPicardEnum);
+	this->inputs->ChangeEnum(VyEnum,VyPicardEnum);
+	this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
+
+	/*Add vx and vy as inputs to the tria element: */
+	this->inputs->AddInput(new TriaInput(VxEnum,vx,P1Enum));
+	this->inputs->AddInput(new TriaInput(VyEnum,vy,P1Enum));
+	this->inputs->AddInput(new TriaInput(VelEnum,vel,P1Enum));
+	this->inputs->AddInput(new TriaInput(PressureEnum,pressure,P1Enum));
+
+	/*Free ressources:*/
+	xDelete<IssmDouble>(pressure);
+	xDelete<IssmDouble>(vel);
+	xDelete<IssmDouble>(vy);
+	xDelete<IssmDouble>(vx);
+	xDelete<IssmDouble>(vvalues);
+	xDelete<IssmDouble>(pvalues);
+	xDelete<int>(vdoflist);
+	xDelete<int>(pdoflist);
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16349)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16350)
@@ -84,4 +84,6 @@
 		int         Sid();
 		bool        IsOnBed();
+		bool        HasEdgeOnBed();
+		void        EdgeOnBedIndices(int* pindex1,int* pindex);
 		bool        IsFloating();
 		bool        IsNodeOnShelfFromFlags(IssmDouble* flags);
@@ -233,8 +235,16 @@
 		ElementMatrix* CreateKMatrixStressbalanceSSAFriction(void);
 		ElementMatrix* CreateKMatrixStressbalanceSIA(void);
+		ElementMatrix* CreateKMatrixStressbalanceFS(void);
+		ElementMatrix* CreateKMatrixStressbalanceFSViscous(void);
+		ElementMatrix* CreateKMatrixStressbalanceFSFriction(void);
 		ElementVector* CreatePVectorStressbalanceSSA(void);
 		ElementVector* CreatePVectorStressbalanceSSADrivingStress(void);
 		ElementVector* CreatePVectorStressbalanceSSAFront(void);
 		ElementVector* CreatePVectorStressbalanceSIA(void);
+		ElementVector* CreatePVectorStressbalanceFS(void);
+		ElementVector* CreatePVectorStressbalanceFSFront(void);
+		ElementVector* CreatePVectorStressbalanceFSViscous(void);
+		void           PVectorGLSstabilization(ElementVector* pe);
+		ElementVector* CreatePVectorStressbalanceFSShelf(void);
 		ElementMatrix* CreateJacobianStressbalanceSSA(void);
 		void	  GetSolutionFromInputsStressbalanceFS(Vector<IssmDouble>* solution);
@@ -242,4 +252,5 @@
 		void	  GetSolutionFromInputsStressbalanceSIA(Vector<IssmDouble>* solution);
 		void	  InputUpdateFromSolutionStressbalanceHoriz( IssmDouble* solution);
+		void	  InputUpdateFromSolutionStressbalanceFS( IssmDouble* solution);
 		void	  InputUpdateFromSolutionStressbalanceSIA( IssmDouble* solution);
 		#endif
Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 16349)
+++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 16350)
@@ -198,4 +198,114 @@
 	/*Clean-up*/
 	xDelete<IssmDouble>(basis);
+}
+/*}}}*/
+/*FUNCTION TriaRef::GetBFS {{{*/
+void TriaRef::GetBFS(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss){
+	/*Compute B  matrix. B=[Bv1 Bv2 ... Bp1 Bp2 ...] where Bvi is of size 3*NDOF3. 
+	 * For node i, Bvi can be expressed in the actual coordinate system
+	 * by: 	   Bvi=[ dphi/dx          0        ]
+	 *					 [   0           dphi/dy     ]
+	 *					 [ 1/2*dphi/dy    1/2*dphi/dx]
+	 *					 [   0             0         ]
+	 *					 [ dphi/dx         dphi/dy   ]
+	 *
+	 * by: 	  Bpi=[  0    ]
+	 *					[  0    ]
+	 *					[  0    ]
+	 *					[ phi_p ]
+	 *					[  0    ]
+	 *	where phi is the interpolation function for node i.
+	 *	Same thing for Bb except the last column that does not exist.
+	 */
+
+	/*Fetch number of nodes for this finite element*/
+	int pnumnodes = this->NumberofNodesPressure();
+	int vnumnodes = this->NumberofNodesVelocity();
+
+	/*Get nodal functions derivatives*/
+	IssmDouble* vdbasis=xNew<IssmDouble>(2*vnumnodes);
+	IssmDouble* pbasis =xNew<IssmDouble>(pnumnodes);
+	GetNodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);
+	GetNodalFunctionsPressure(pbasis,gauss);
+
+	/*Build B: */
+	for(int i=0;i<vnumnodes;i++){
+		B[(2*vnumnodes+pnumnodes)*0+2*i+0] = vdbasis[0*vnumnodes+i];
+		B[(2*vnumnodes+pnumnodes)*0+2*i+1] = 0.;
+		B[(2*vnumnodes+pnumnodes)*1+2*i+0] = 0.;
+		B[(2*vnumnodes+pnumnodes)*1+2*i+1] = vdbasis[1*vnumnodes+i];
+		B[(2*vnumnodes+pnumnodes)*2+2*i+0] = .5*vdbasis[1*vnumnodes+i];
+		B[(2*vnumnodes+pnumnodes)*2+2*i+1] = .5*vdbasis[0*vnumnodes+i];
+		B[(2*vnumnodes+pnumnodes)*3+2*i+0] = 0.;
+		B[(2*vnumnodes+pnumnodes)*3+2*i+1] = 0.;
+		B[(2*vnumnodes+pnumnodes)*4+2*i+0] = vdbasis[0*vnumnodes+i];
+		B[(2*vnumnodes+pnumnodes)*4+2*i+1] = vdbasis[1*vnumnodes+i];
+	}
+	for(int i=0;i<pnumnodes;i++){
+		B[(2*vnumnodes+pnumnodes)*0+(2*vnumnodes)+i] = 0.;
+		B[(2*vnumnodes+pnumnodes)*1+(2*vnumnodes)+i] = 0.;
+		B[(2*vnumnodes+pnumnodes)*2+(2*vnumnodes)+i] = 0.;
+		B[(2*vnumnodes+pnumnodes)*3+(2*vnumnodes)+i] = pbasis[i];
+		B[(2*vnumnodes+pnumnodes)*4+(2*vnumnodes)+i] = 0.;
+	}
+
+	/*Clean up*/
+	xDelete<IssmDouble>(vdbasis);
+	xDelete<IssmDouble>(pbasis);
+}
+/*}}}*/
+/*FUNCTION TriaRef::GetBprimeFS {{{*/
+void TriaRef::GetBprimeFS(IssmDouble* B_prime, IssmDouble* xyz_list, GaussTria* gauss){
+	/*	Compute B'  matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2. 
+	 *	For node i, Bi' can be expressed in the actual coordinate system
+	 *	by: 
+	 *			Bvi' = [  dphi/dx     0     ]
+	 *					 [     0      dphi/dy ]
+	 *					 [  dphi/dy   dphi/dx ]
+	 *					 [  dphi/dx   dphi/dy ]
+	 *					 [     0      0       ]
+	 *
+	 * by: 	  Bpi=[  0  ]
+	 *					[  0  ]
+	 *					[  0  ]
+	 *					[  0  ]
+	 *					[ phi ]
+	 *	where phi is the interpolation function for node i.
+	 */
+
+	/*Fetch number of nodes for this finite element*/
+	int pnumnodes = this->NumberofNodesPressure();
+	int vnumnodes = this->NumberofNodesVelocity();
+
+	/*Get nodal functions derivatives*/
+	IssmDouble* vdbasis=xNew<IssmDouble>(2*vnumnodes);
+	IssmDouble* pbasis =xNew<IssmDouble>(pnumnodes);
+	GetNodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);
+	GetNodalFunctionsPressure(pbasis,gauss);
+
+	/*Build B_prime: */
+	for(int i=0;i<vnumnodes;i++){
+		B_prime[(2*vnumnodes+pnumnodes)*0+2*i+0] = vdbasis[0*vnumnodes+i];
+		B_prime[(2*vnumnodes+pnumnodes)*0+2*i+1] = 0.;
+		B_prime[(2*vnumnodes+pnumnodes)*1+2*i+0] = 0.;
+		B_prime[(2*vnumnodes+pnumnodes)*1+2*i+1] = vdbasis[1*vnumnodes+i];
+		B_prime[(2*vnumnodes+pnumnodes)*2+2*i+0] = vdbasis[1*vnumnodes+i];
+		B_prime[(2*vnumnodes+pnumnodes)*2+2*i+1] = vdbasis[0*vnumnodes+i];
+		B_prime[(2*vnumnodes+pnumnodes)*3+2*i+0] = vdbasis[0*vnumnodes+i];
+		B_prime[(2*vnumnodes+pnumnodes)*3+2*i+1] = vdbasis[1*vnumnodes+i];
+		B_prime[(2*vnumnodes+pnumnodes)*4+2*i+0] = 0.;
+		B_prime[(2*vnumnodes+pnumnodes)*4+2*i+1] = 0.;
+	}
+	for(int i=0;i<pnumnodes;i++){
+		B_prime[(2*vnumnodes+pnumnodes)*0+(2*vnumnodes)+i] = 0.;
+		B_prime[(2*vnumnodes+pnumnodes)*1+(2*vnumnodes)+i] = 0.;
+		B_prime[(2*vnumnodes+pnumnodes)*2+(2*vnumnodes)+i] = 0.;
+		B_prime[(2*vnumnodes+pnumnodes)*3+(2*vnumnodes)+i] = 0.;
+		B_prime[(2*vnumnodes+pnumnodes)*4+(2*vnumnodes)+i] = pbasis[i];
+	}
+
+	/*Clean up*/
+	xDelete<IssmDouble>(vdbasis);
+	xDelete<IssmDouble>(pbasis);
 }
 /*}}}*/
@@ -458,4 +568,74 @@
 }
 /*}}}*/
+/*FUNCTION TriaRef::GetNodalFunctionsVelocity{{{*/
+void TriaRef::GetNodalFunctionsVelocity(IssmDouble* basis,GaussTria* gauss){
+	/*This routine returns the values of the nodal functions  at the gaussian point.*/
+
+	switch(this->element_type){
+		case P1P1Enum:
+			this->element_type = P1Enum;
+			this->GetNodalFunctions(basis,gauss);
+			this->element_type = P1P1Enum;
+			return;
+		case P1P1GLSEnum:
+			this->element_type = P1Enum;
+			this->GetNodalFunctions(basis,gauss);
+			this->element_type = P1P1GLSEnum;
+			return;
+		case MINIcondensedEnum:
+			this->element_type = P1bubbleEnum;
+			this->GetNodalFunctions(basis,gauss);
+			this->element_type = MINIcondensedEnum;
+			return;
+		case MINIEnum:
+			this->element_type = P1bubbleEnum;
+			this->GetNodalFunctions(basis,gauss);
+			this->element_type = MINIEnum;
+			return;
+		case TaylorHoodEnum:
+			this->element_type = P2Enum;
+			this->GetNodalFunctions(basis,gauss);
+			this->element_type = TaylorHoodEnum;
+			return;
+		default:
+			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+	}
+}
+/*}}}*/
+/*FUNCTION TriaRef::GetNodalFunctionsPressure{{{*/
+void TriaRef::GetNodalFunctionsPressure(IssmDouble* basis,GaussTria* gauss){
+	/*This routine returns the values of the nodal functions  at the gaussian point.*/
+
+	switch(this->element_type){
+		case P1P1Enum:
+			this->element_type = P1Enum;
+			this->GetNodalFunctions(basis,gauss);
+			this->element_type = P1P1Enum;
+			return;
+		case P1P1GLSEnum:
+			this->element_type = P1Enum;
+			this->GetNodalFunctions(basis,gauss);
+			this->element_type = P1P1GLSEnum;
+			return;
+		case MINIcondensedEnum:
+			this->element_type = P1Enum;
+			this->GetNodalFunctions(basis,gauss);
+			this->element_type = MINIcondensedEnum;
+			return;
+		case MINIEnum:
+			this->element_type = P1Enum;
+			this->GetNodalFunctions(basis,gauss);
+			this->element_type = MINIEnum;
+			return;
+		case TaylorHoodEnum:
+			this->element_type = P1Enum;
+			this->GetNodalFunctions(basis,gauss);
+			this->element_type = TaylorHoodEnum;
+			return;
+		default:
+			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+	}
+}
+/*}}}*/
 /*FUNCTION TriaRef::GetSegmentNodalFunctions{{{*/
 void TriaRef::GetSegmentNodalFunctions(IssmDouble* basis,GaussTria* gauss,int index1,int index2){
@@ -527,4 +707,70 @@
 	xDelete<IssmDouble>(dbasis_ref);
 
+}
+/*}}}*/
+/*FUNCTION TriaRef::GetNodalFunctionsDerivativesPressure{{{*/
+void TriaRef::GetNodalFunctionsDerivativesPressure(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTria* gauss){
+	switch(this->element_type){
+		case P1P1Enum:
+			this->element_type = P1Enum;
+			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+			this->element_type = P1P1Enum;
+			return;
+		case P1P1GLSEnum:
+			this->element_type = P1Enum;
+			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+			this->element_type = P1P1GLSEnum;
+			return;
+		case MINIcondensedEnum:
+			this->element_type = P1Enum;
+			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+			this->element_type = MINIcondensedEnum;
+			return;
+		case MINIEnum:
+			this->element_type = P1Enum;
+			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+			this->element_type = MINIEnum;
+			return;
+		case TaylorHoodEnum:
+			this->element_type = P1Enum;
+			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+			this->element_type = TaylorHoodEnum;
+			return;
+		default:
+			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+	}
+}
+/*}}}*/
+/*FUNCTION TriaRef::GetNodalFunctionsDerivativesVelocity{{{*/
+void TriaRef::GetNodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTria* gauss){
+	switch(this->element_type){
+		case P1P1Enum:
+			this->element_type = P1Enum;
+			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+			this->element_type = P1P1Enum;
+			return;
+		case P1P1GLSEnum:
+			this->element_type = P1Enum;
+			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+			this->element_type = P1P1GLSEnum;
+			return;
+		case MINIcondensedEnum:
+			this->element_type = P1bubbleEnum;
+			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+			this->element_type = MINIcondensedEnum;
+			return;
+		case MINIEnum:
+			this->element_type = P1bubbleEnum;
+			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+			this->element_type = MINIEnum;
+			return;
+		case TaylorHoodEnum:
+			this->element_type = P2Enum;
+			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+			this->element_type = TaylorHoodEnum;
+			return;
+		default:
+			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+	}
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h	(revision 16349)
+++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h	(revision 16350)
@@ -23,6 +23,8 @@
 
 		/*Numerics*/
+		void GetBFS(IssmDouble* B_prime, IssmDouble* xyz_list, GaussTria* gauss);
 		void GetBSSA(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss);
 		void GetBSSAFS(IssmDouble* B , IssmDouble* xyz_list, GaussTria* gauss);
+		void GetBprimeFS(IssmDouble* B_prime, IssmDouble* xyz_list, GaussTria* gauss);
 		void GetBprimeSSA(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss);
 		void GetBprimeSSAFS(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss);
@@ -36,8 +38,12 @@
 		void GetJacobianInvert(IssmDouble*  Jinv, IssmDouble* xyz_list,GaussTria* gauss);
 		void GetNodalFunctions(IssmDouble* basis,GaussTria* gauss);
+		void GetNodalFunctionsVelocity(IssmDouble* basis, GaussTria* gauss);
+		void GetNodalFunctionsPressure(IssmDouble* basis, GaussTria* gauss);
 		void GetSegmentNodalFunctions(IssmDouble* basis,GaussTria* gauss, int index1,int index2);
 		void GetSegmentBFlux(IssmDouble* B,GaussTria* gauss, int index1,int index2);
 		void GetSegmentBprimeFlux(IssmDouble* Bprime,GaussTria* gauss, int index1,int index2);
 		void GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTria* gauss);
+		void GetNodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,GaussTria* gauss);
+		void GetNodalFunctionsDerivativesPressure(IssmDouble* dbasis,IssmDouble* xyz_list,GaussTria* gauss);
 		void GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussTria* gauss);
 		void GetInputValue(IssmDouble* pp, IssmDouble* plist, GaussTria* gauss);
Index: /issm/trunk-jpl/src/c/classes/Materials/Material.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Material.h	(revision 16349)
+++ /issm/trunk-jpl/src/c/classes/Materials/Material.h	(revision 16350)
@@ -27,4 +27,5 @@
 		virtual void       GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum)=0;
 		virtual void       GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon)=0;
+		virtual void       GetViscosity2dvertical(IssmDouble* pviscosity, IssmDouble* pepsilon)=0;
 		virtual void       GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon)=0;
 		virtual void       GetViscosity3dFS(IssmDouble* pviscosity3d, IssmDouble* epsilon)=0;
Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp	(revision 16349)
+++ /issm/trunk-jpl/src/c/classes/Materials/Matice.cpp	(revision 16350)
@@ -278,4 +278,70 @@
 				e=(n-1)/(2*n);
 				viscosity=(1-D)*B/(2*pow(A,e));
+			}
+		}
+	}
+
+	/*Checks in debugging mode*/
+	if(viscosity<=0) _error_("Negative viscosity");
+	_assert_(B>0);
+	_assert_(n>0);
+	_assert_(D>=0 && D<1);
+
+	/*Return: */
+	*pviscosity=viscosity;
+}
+/*}}}*/
+/*FUNCTION Matice::GetViscosity2dvertical {{{*/
+void  Matice::GetViscosity2dvertical(IssmDouble* pviscosity, IssmDouble* epsilon){
+	/*From a string tensor and a material object, return viscosity, using Glen's flow law.
+									   (1-D) B
+	  viscosity= --------------------------------------
+						  2[ exx^2+eyy^2+ 2exy^2]^[(n-1)/2n]
+
+	  where viscosity is the viscotiy, B the flow law parameter , (u,v) the velocity 
+	  vector, and n the flow law exponent.
+
+	  If epsilon is NULL, it means this is the first time SystemMatrices is being run, and we 
+	  return 10^14, initial viscosity.
+	  */
+
+	/*output: */
+	IssmDouble viscosity;
+
+	/*input strain rate: */
+	IssmDouble exx,eyy,exy;
+
+	/*Intermediary: */
+	IssmDouble A,e;
+	IssmDouble B,D,n;
+
+	/*Get B and n*/
+	B=GetB();
+	n=GetN();
+	D=GetD();
+
+	if (n==1){
+		/*Viscous behaviour! viscosity=B: */
+		viscosity=(1-D)*B/2;
+	}
+	else{
+		if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0)){
+			viscosity=0.5*pow(10.,14);
+		}
+		else{
+			/*Retrive strain rate components: */
+			exx=epsilon[0];
+			eyy=epsilon[1];
+			exy=epsilon[2];
+
+			/*Build viscosity: viscosity=B/(2*A^e) */
+			A=exx*exx+eyy*eyy+2.*exy*exy;
+			if(A==0.){
+				/*Maxiviscositym viscosity for 0 shear areas: */
+				viscosity=2.5*2.e+17;
+			}
+			else{
+				e=(n-1.)/(2.*n);
+				viscosity=(1.-D)*B/(2.*pow(A,e));
 			}
 		}
@@ -793,8 +859,8 @@
 
 		/*Intermediaries*/
-		const int num_vertices = 3; //Tria has 3 vertices
-		IssmDouble    nodeinputs[num_vertices];
-		IssmDouble    cmmininputs[num_vertices];
-		IssmDouble    cmmaxinputs[num_vertices];
+		const int  num_vertices = 3; //Tria has 3 vertices
+		IssmDouble nodeinputs[num_vertices];
+		IssmDouble cmmininputs[num_vertices];
+		IssmDouble cmmaxinputs[num_vertices];
 
 		/*Get B*/
@@ -867,5 +933,5 @@
 		if (iomodel->Data(DamageDEnum)) {
 			for(i=0;i<num_vertices;i++) nodeinputs[i]=iomodel->Data(DamageDEnum)[iomodel->elements[num_vertices*index+i]-1];
-			this->inputs->AddInput(new TriaInput(DamageDbarEnum,nodeinputs,P1Enum));
+			this->inputs->AddInput(new TriaInput(DamageDEnum,nodeinputs,P1Enum));
 		}
 	}
Index: /issm/trunk-jpl/src/c/classes/Materials/Matice.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Matice.h	(revision 16349)
+++ /issm/trunk-jpl/src/c/classes/Materials/Matice.h	(revision 16350)
@@ -54,4 +54,5 @@
 		void       SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Vertices* vertices,Materials* materials,Parameters* parameters);
 		void       GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon);
+		void       GetViscosity2dvertical(IssmDouble* pviscosity, IssmDouble* pepsilon);
 		void       GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon);
 		void       GetViscosity3dFS(IssmDouble* pviscosity3d, IssmDouble* epsilon);
Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Materials/Matpar.h	(revision 16349)
+++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.h	(revision 16350)
@@ -80,8 +80,9 @@
 		/*}}}*/
 		/*Material virtual functions resolution: {{{*/
-		void   InputDuplicate(int original_enum,int new_enum);
-		void   Configure(Elements* elements);
-		void   GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum){return;}
+		void       InputDuplicate(int original_enum,int new_enum);
+		void       Configure(Elements* elements);
+		void       GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum){return;}
 		void       GetViscosity2d(IssmDouble* pviscosity, IssmDouble* pepsilon){_error_("not supported");};
+		void       GetViscosity2dvertical(IssmDouble* pviscosity, IssmDouble* pepsilon){_error_("not supported");};
 		void       GetViscosity3d(IssmDouble* pviscosity3d, IssmDouble* pepsilon){_error_("not supported");};
 		void       GetViscosity3dFS(IssmDouble* pviscosity3d, IssmDouble* epsilon){_error_("not supported");};