Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16939)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16940)
@@ -855,4 +855,6 @@
 		case HOFSApproximationEnum: 
 			return CreatePVectorHOFS(element);
+		case SSAFSApproximationEnum: 
+			return CreatePVectorSSAFS(element);
 		case NoneApproximationEnum:
 			return NULL;
@@ -3356,4 +3358,27 @@
 	return pe;
 }/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorSSAFS(Element* element){/*{{{*/
+
+	/*compute all load vectors for this element*/
+	int init = element->FiniteElement();
+	element->SetTemporaryElementType(P1Enum); // P1 needed for HO
+	ElementVector* pe1=CreatePVectorSSA(element);
+	element->SetTemporaryElementType(init); // P1 needed for HO
+	ElementVector* pe2=CreatePVectorFS(element);
+	int indices[3]={18,19,20};
+	element->SetTemporaryElementType(MINIcondensedEnum); // P1 needed for HO
+	ElementMatrix* Ke = CreateKMatrixFS(element);
+	element->SetTemporaryElementType(init); // P1 needed for HO
+	pe2->StaticCondensation(Ke,3,&indices[0]);
+	delete Ke;
+	ElementVector* pe3=CreatePVectorCouplingSSAFS(element);
+	ElementVector* pe =new ElementVector(pe1,pe2,pe3);
+
+	/*clean-up and return*/
+	delete pe1;
+	delete pe2;
+	delete pe3;
+	return pe;
+}/*}}}*/
 ElementVector* StressbalanceAnalysis::CreatePVectorHOFS(Element* element){/*{{{*/
 
@@ -3543,4 +3568,161 @@
 	return pe;
 }/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorCouplingSSAFS(Element* element){/*{{{*/
+
+	/*compute all load vectors for this element*/
+	ElementVector* pe1=CreatePVectorCouplingSSAFSViscous(element);
+	ElementVector* pe2=CreatePVectorCouplingSSAFSFriction(element);
+	ElementVector* pe =new ElementVector(pe1,pe2);
+
+	/*clean-up and return*/
+	delete pe1;
+	delete pe2;
+	return pe;
+}/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorCouplingSSAFSFriction(Element* element){/*{{{*/
+
+	/*Intermediaries*/
+	int         i,j,approximation;
+	int         dim=3;
+	IssmDouble  Jdet,Jdet2d,FSreconditioning;
+	IssmDouble	bed_normal[3];
+	IssmDouble  viscosity, w, alpha2_gauss;
+	IssmDouble  dw[3];
+	IssmDouble  basis[6]; //for the six nodes of the penta
+	IssmDouble	*xyz_list_tria = NULL;
+	IssmDouble  *xyz_list      = NULL;
+
+	/*Initialize Element vector and return if necessary*/
+	if(!element->IsOnBed() || element->IsFloating()) return NULL;
+	element->GetInputValue(&approximation,ApproximationEnum);
+	if(approximation!=SSAFSApproximationEnum) return NULL;
+	int vnumnodes = element->GetNumberOfNodesVelocity();
+	int pnumnodes = element->GetNumberOfNodesPressure();
+
+	/*Prepare coordinate system list*/
+	int* cs_list     = xNew<int>(vnumnodes+pnumnodes);
+	Node **node_list = xNew<Node*>(vnumnodes+pnumnodes);
+	for(i=0;i<vnumnodes;i++){
+		cs_list[i]             = XYZEnum;
+		node_list[i]           = element->GetNode(i);
+	}
+	for(i=0;i<pnumnodes;i++){
+		cs_list[vnumnodes+i]   = PressureEnum;
+		node_list[vnumnodes+i] = element->GetNode(vnumnodes+i);
+	}
+	ElementVector* pe=element->NewElementVector(FSvelocityEnum);
+
+	/*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);
+	Input* vzSSA_input=element->GetInput(VzSSAEnum); _assert_(vzSSA_input);
+
+	/*build friction object, used later on: */
+	Friction* friction=new Friction(element,3);
+
+	/* Start looping on the number of gauss 2d (nodes on the bedrock) */
+	Gauss* gauss=element->NewGaussBase(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminantBase(&Jdet2d,xyz_list_tria,gauss);
+		element->NodalFunctionsP1(basis, gauss);
+
+		vzSSA_input->GetInputValue(&w, gauss);
+		vzSSA_input->GetInputDerivativeValue(&dw[0],xyz_list,gauss);
+
+		element->NormalBase(&bed_normal[0],xyz_list_tria);
+		element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
+		friction->GetAlpha2(&alpha2_gauss, gauss,vx_input,vy_input,vz_input);
+
+		for(i=0;i<3;i++){
+			pe->values[i*3+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[0])*basis[i];
+			pe->values[i*3+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[1])*basis[i];
+			pe->values[i*3+2]+=Jdet2d*gauss->weight*2*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*basis[i];
+		}
+	}
+
+	/*Transform coordinate system*/
+	element->TransformLoadVectorCoord(pe,node_list,vnumnodes+pnumnodes,cs_list);
+
+	/*Clean up and return*/
+	xDelete<int>(cs_list);
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(xyz_list_tria);
+	xDelete<Node*>(node_list);
+	delete gauss;
+	delete friction;
+	return pe;
+}/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorCouplingSSAFSViscous(Element* element){/*{{{*/
+
+	/*Intermediaries */
+	int         i,approximation;
+	IssmDouble  viscosity,Jdet,FSreconditioning;
+	IssmDouble  dw[3];
+	IssmDouble  *xyz_list = NULL;
+	IssmDouble  basis[6]; //for the six nodes of the penta
+	IssmDouble  dbasis[3][6]; //for the six nodes of the penta
+
+	/*Initialize Element vector and return if necessary*/
+	element->GetInputValue(&approximation,ApproximationEnum);
+	if(approximation!=SSAFSApproximationEnum) return NULL;
+	int vnumnodes = element->GetNumberOfNodesVelocity();
+	int pnumnodes = element->GetNumberOfNodesPressure();
+
+	/*Prepare coordinate system list*/
+	int* cs_list = xNew<int>(vnumnodes+pnumnodes);
+	Node **node_list = xNew<Node*>(vnumnodes+pnumnodes);
+	for(i=0;i<vnumnodes;i++){
+		cs_list[i]             = XYZEnum;
+		node_list[i]           = element->GetNode(i);
+	}
+	for(i=0;i<pnumnodes;i++){
+		cs_list[vnumnodes+i]   = PressureEnum;
+		node_list[vnumnodes+i] = element->GetNode(vnumnodes+i);
+	}
+	ElementVector* pe=element->NewElementVector(FSvelocityEnum);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+	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);
+	Input* vzSSA_input=element->GetInput(VzSSAEnum);   _assert_(vzSSA_input);
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=element->NewGauss(5);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		element->NodalFunctionsP1(&basis[0], gauss);
+		element->NodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss);
+
+		vzSSA_input->GetInputDerivativeValue(&dw[0],xyz_list,gauss);
+		element->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input);
+
+		for(i=0;i<6;i++){
+			pe->values[i*3+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dbasis[2][i];
+			pe->values[i*3+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dbasis[2][i];
+			pe->values[i*3+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dbasis[0][i]+dw[1]*dbasis[1][i]+2*dw[2]*dbasis[2][i]);
+			pe->values[3*vnumnodes+i]+=Jdet*gauss->weight*FSreconditioning*dw[2]*basis[i];
+		}
+	}
+
+	/*Transform coordinate system*/
+	element->TransformLoadVectorCoord(pe,node_list,vnumnodes+pnumnodes,cs_list);
+
+	/*Clean up and return*/
+	xDelete<int>(cs_list);
+	delete gauss;
+	return pe;
+}/*}}}*/
 void StressbalanceAnalysis::GetBSSAHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*NDOF2. 
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 16939)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 16940)
@@ -78,4 +78,8 @@
 		ElementMatrix* CreateKMatrixCouplingSSAHOViscous(Element* element);
 		ElementVector* CreatePVectorSSAHO(Element* element);
+		ElementVector* CreatePVectorSSAFS(Element* element);
+		ElementVector* CreatePVectorCouplingSSAFS(Element* element);
+		ElementVector* CreatePVectorCouplingSSAFSFriction(Element* element);
+		ElementVector* CreatePVectorCouplingSSAFSViscous(Element* element);
 		ElementVector* CreatePVectorHOFS(Element* element);
 		ElementVector* CreatePVectorCouplingHOFS(Element* element);