Index: /issm/trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp	(revision 16830)
+++ /issm/trunk-jpl/src/c/analyses/L2ProjectionBaseAnalysis.cpp	(revision 16831)
@@ -56,5 +56,57 @@
 /*Finite Element Analysis*/
 ElementMatrix* L2ProjectionBaseAnalysis::CreateKMatrix(Element* element){/*{{{*/
-	_error_("not implemented yet");
+
+	/*Intermediaries*/
+	int      meshtype;
+	Element* basalelement;
+
+	/*Get basal element*/
+	element->FindParam(&meshtype,MeshTypeEnum);
+	switch(meshtype){
+		case Mesh2DhorizontalEnum:
+			basalelement = element;
+			break;
+		case Mesh3DEnum:
+			if(!element->IsOnBed()) return NULL;
+			basalelement = element->SpawnBasalElement();
+			break;
+		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+	}
+
+	/*Intermediaries */
+	IssmDouble  D,Jdet;
+	IssmDouble *xyz_list  = NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = basalelement->GetNumberOfNodes();
+
+	/*Initialize Element vector*/
+	ElementMatrix* Ke    = basalelement->NewElementMatrix(NoneApproximationEnum);
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	basalelement->GetVerticesCoordinates(&xyz_list);
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=basalelement->NewGauss(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		basalelement->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		basalelement->NodalFunctions(basis,gauss);
+		D=gauss->weight*Jdet;
+
+		TripleMultiply(basis,1,numnodes,1,
+					&D,1,1,0,
+					basis,1,numnodes,0,
+					&Ke->values[0],1);
+	}
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(basis);
+	delete gauss;
+	if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
+	return Ke;
 }/*}}}*/
 ElementVector* L2ProjectionBaseAnalysis::CreatePVector(Element* element){/*{{{*/
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16830)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16831)
@@ -820,4 +820,6 @@
 		case HOApproximationEnum: 
 			return CreateKMatrixHO(element);
+		case FSApproximationEnum: 
+			return CreateKMatrixFS(element);
 		case NoneApproximationEnum:
 			return NULL;
@@ -1658,4 +1660,91 @@
 
 /*FS*/
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixFS(Element* element){/*{{{*/
+
+	/*compute all stiffness matrices for this element*/
+	ElementMatrix* Ke1=CreateKMatrixFSViscous(element);
+	ElementMatrix* Ke2=CreateKMatrixFSFriction(element);
+	ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
+
+	/*clean-up and return*/
+	delete Ke1;
+	delete Ke2;
+	return Ke;
+}/*}}}*/
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSViscous(Element* element){/*{{{*/
+
+	/*Intermediaries*/
+	int         i,meshtype,dim,epssize;
+	IssmDouble  viscosity,FSreconditioning,Jdet;
+	IssmDouble *xyz_list = NULL;
+
+	/*Get problem dimension*/
+	element->FindParam(&meshtype,MeshTypeEnum);
+	switch(meshtype){
+		case Mesh2DverticalEnum: dim = 2; epssize = 3; break;
+		case Mesh3DEnum:         dim = 3; epssize = 6; 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;
+	int bsize     = epssize + 2;
+
+	/*Prepare coordinate system list*/
+	int* cs_list = xNew<int>(vnumnodes+pnumnodes);
+	if(dim==2) for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
+	else       for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
+	for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
+
+	/*Initialize Element matrix and vectors*/
+	ElementMatrix* Ke     = element->NewElementMatrix(FSvelocityEnum);
+	IssmDouble*    B      = xNew<IssmDouble>(bsize*numdof);
+	IssmDouble*    Bprime = xNew<IssmDouble>(bsize*numdof);
+	IssmDouble*    D      = xNewZeroInit<IssmDouble>(bsize*bsize);
+
+	/*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;
+	if(dim==3){vz_input=element->GetInput(VzEnum); _assert_(vz_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);
+		this->GetBFS(B,element,dim,xyz_list,gauss);
+		this->GetBFSprime(Bprime,element,dim,xyz_list,gauss);
+
+		element->ViscosityFS(&viscosity,xyz_list,gauss,vx_input,vy_input,vz_input);
+		for(i=0;i<epssize;i++)     D[i*bsize+i] = + 2.*viscosity*gauss->weight*Jdet;
+		for(i=epssize;i<bsize;i++) D[i*bsize+i] = - FSreconditioning*gauss->weight*Jdet;
+
+		TripleMultiply(B,bsize,numdof,1,
+					D,bsize,bsize,0,
+					Bprime,bsize,numdof,0,
+					&Ke->values[0],1);
+	}
+
+	/*Transform Coordinate System*/
+	element->TransformStiffnessMatrixCoord(Ke,cs_list);
+
+	/*Clean up and return*/
+	delete gauss;
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(D);
+	xDelete<IssmDouble>(Bprime);
+	xDelete<IssmDouble>(B);
+	xDelete<int>(cs_list);
+	return Ke;
+}/*}}}*/
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSFriction(Element* element){/*{{{*/
+	return NULL;
+}/*}}}*/
 ElementVector* StressbalanceAnalysis::CreatePVectorFS(Element* element){/*{{{*/
 
@@ -1816,4 +1905,229 @@
 
 	return NULL;
+}/*}}}*/
+void StressbalanceAnalysis::GetBFS(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* 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   ]
+	 *
+	 *         Bpi=[  0    ]
+	 *					[  0    ]
+	 *					[  0    ]
+	 *					[ phi_p ]
+	 *					[  0    ]
+	 *
+	 *	In 3d:
+	 *     	   Bvi=[ dh/dx          0             0      ]
+	 *					[   0           dh/dy           0      ]
+	 *					[   0             0           dh/dz    ]
+	 *					[ 1/2*dh/dy    1/2*dh/dx        0      ]
+	 *					[ 1/2*dh/dz       0         1/2*dh/dx  ]
+	 *					[   0          1/2*dh/dz    1/2*dh/dy  ]
+	 *					[   0             0             0      ]
+	 *					[ dh/dx         dh/dy         dh/dz    ]
+	 *
+	 *         Bpi=[ 0 ]
+	 *					[ 0 ]
+	 *					[ 0 ]
+	 *					[ 0 ]
+	 *					[ 0 ]
+	 *					[ 0 ]
+	 *					[ h ]
+	 *					[ 0 ]
+	 *	where phi is the finiteelement 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 = element->NumberofNodesPressure();
+	int vnumnodes = element->NumberofNodesVelocity();
+
+	/*Get nodal functions derivatives*/
+	IssmDouble* vdbasis=xNew<IssmDouble>(dim*vnumnodes);
+	IssmDouble* pbasis =xNew<IssmDouble>(pnumnodes);
+	element->NodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);
+	element->NodalFunctionsPressure(pbasis,gauss);
+
+	/*Build B: */
+	if(dim==2){
+		for(int i=0;i<vnumnodes;i++){
+			B[(dim*vnumnodes+pnumnodes)*0+dim*i+0] = vdbasis[0*vnumnodes+i];
+			B[(dim*vnumnodes+pnumnodes)*0+dim*i+1] = 0.;
+			B[(dim*vnumnodes+pnumnodes)*1+dim*i+0] = 0.;
+			B[(dim*vnumnodes+pnumnodes)*1+dim*i+1] = vdbasis[1*vnumnodes+i];
+			B[(dim*vnumnodes+pnumnodes)*2+dim*i+0] = .5*vdbasis[1*vnumnodes+i];
+			B[(dim*vnumnodes+pnumnodes)*2+dim*i+1] = .5*vdbasis[0*vnumnodes+i];
+			B[(dim*vnumnodes+pnumnodes)*3+dim*i+0] = 0.;
+			B[(dim*vnumnodes+pnumnodes)*3+dim*i+1] = 0.;
+			B[(dim*vnumnodes+pnumnodes)*4+dim*i+0] = vdbasis[0*vnumnodes+i];
+			B[(dim*vnumnodes+pnumnodes)*4+dim*i+1] = vdbasis[1*vnumnodes+i];
+		}
+		for(int i=0;i<pnumnodes;i++){
+			B[(dim*vnumnodes+pnumnodes)*0+(dim*vnumnodes)+i] = 0.;
+			B[(dim*vnumnodes+pnumnodes)*1+(dim*vnumnodes)+i] = 0.;
+			B[(dim*vnumnodes+pnumnodes)*2+(dim*vnumnodes)+i] = 0.;
+			B[(dim*vnumnodes+pnumnodes)*3+(dim*vnumnodes)+i] = pbasis[i];
+			B[(dim*vnumnodes+pnumnodes)*4+(dim*vnumnodes)+i] = 0.;
+		}
+	}
+	else{
+		for(int i=0;i<vnumnodes;i++){
+			B[(dim*vnumnodes+pnumnodes)*0+dim*i+0] = vdbasis[0*vnumnodes+i];
+			B[(dim*vnumnodes+pnumnodes)*0+dim*i+1] = 0.;
+			B[(dim*vnumnodes+pnumnodes)*0+dim*i+2] = 0.;
+			B[(dim*vnumnodes+pnumnodes)*1+dim*i+0] = 0.;
+			B[(dim*vnumnodes+pnumnodes)*1+dim*i+1] = vdbasis[1*vnumnodes+i];
+			B[(dim*vnumnodes+pnumnodes)*1+dim*i+2] = 0.;
+			B[(dim*vnumnodes+pnumnodes)*2+dim*i+0] = 0.;
+			B[(dim*vnumnodes+pnumnodes)*2+dim*i+1] = 0.;
+			B[(dim*vnumnodes+pnumnodes)*2+dim*i+2] = vdbasis[2*vnumnodes+i];
+			B[(dim*vnumnodes+pnumnodes)*3+dim*i+0] = .5*vdbasis[1*vnumnodes+i];
+			B[(dim*vnumnodes+pnumnodes)*3+dim*i+1] = .5*vdbasis[0*vnumnodes+i];
+			B[(dim*vnumnodes+pnumnodes)*3+dim*i+2] = 0.;
+			B[(dim*vnumnodes+pnumnodes)*4+dim*i+0] = .5*vdbasis[2*vnumnodes+i];
+			B[(dim*vnumnodes+pnumnodes)*4+dim*i+1] = 0.;
+			B[(dim*vnumnodes+pnumnodes)*4+dim*i+2] = .5*vdbasis[0*vnumnodes+i];
+			B[(dim*vnumnodes+pnumnodes)*5+dim*i+0] = 0.;
+			B[(dim*vnumnodes+pnumnodes)*5+dim*i+1] = .5*vdbasis[2*vnumnodes+i];
+			B[(dim*vnumnodes+pnumnodes)*5+dim*i+2] = .5*vdbasis[1*vnumnodes+i];
+			B[(dim*vnumnodes+pnumnodes)*6+dim*i+0] = 0.;
+			B[(dim*vnumnodes+pnumnodes)*6+dim*i+1] = 0.;
+			B[(dim*vnumnodes+pnumnodes)*6+dim*i+2] = 0.;
+			B[(dim*vnumnodes+pnumnodes)*7+dim*i+0] = vdbasis[0*vnumnodes+i];
+			B[(dim*vnumnodes+pnumnodes)*7+dim*i+1] = vdbasis[1*vnumnodes+i];
+			B[(dim*vnumnodes+pnumnodes)*7+dim*i+2] = vdbasis[2*vnumnodes+i];
+		}
+		for(int i=0;i<pnumnodes;i++){
+			B[(dim*vnumnodes+pnumnodes)*0+(dim*vnumnodes)+i] = 0.;
+			B[(dim*vnumnodes+pnumnodes)*1+(dim*vnumnodes)+i] = 0.;
+			B[(dim*vnumnodes+pnumnodes)*2+(dim*vnumnodes)+i] = 0.;
+			B[(dim*vnumnodes+pnumnodes)*3+(dim*vnumnodes)+i] = 0.;
+			B[(dim*vnumnodes+pnumnodes)*4+(dim*vnumnodes)+i] = 0.;
+			B[(dim*vnumnodes+pnumnodes)*5+(dim*vnumnodes)+i] = 0.;
+			B[(dim*vnumnodes+pnumnodes)*6+(dim*vnumnodes)+i] = pbasis[i];
+			B[(dim*vnumnodes+pnumnodes)*7+(dim*vnumnodes)+i] = 0.;
+		}
+	}
+
+	/*Clean up*/
+	xDelete<IssmDouble>(vdbasis);
+	xDelete<IssmDouble>(pbasis);
+}/*}}}*/
+void StressbalanceAnalysis::GetBFSprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* 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 ]
+	 *
+	 *	In 3d
+	 *     	   Bvi=[ dh/dx          0             0      ]
+	 *					[   0           dh/dy           0      ]
+	 *					[   0             0           dh/dz    ]
+	 *					[ 1/2*dh/dy    1/2*dh/dx        0      ]
+	 *					[ 1/2*dh/dz       0         1/2*dh/dx  ]
+	 *					[   0          1/2*dh/dz    1/2*dh/dy  ]
+	 *					[   0             0             0      ]
+	 *					[ dh/dx         dh/dy         dh/dz    ]
+	 *
+	 *         Bpi=[ 0 ]
+	 *					[ 0 ]
+	 *					[ 0 ]
+	 *					[ 0 ]
+	 *					[ 0 ]
+	 *					[ 0 ]
+	 *					[ h ]
+	 *					[ 0 ]
+	 *	where phi is the finiteelement function for node i.
+	 *	In 3d:
+	 */
+
+	/*Fetch number of nodes for this finite element*/
+	int pnumnodes = element->NumberofNodesPressure();
+	int vnumnodes = element->NumberofNodesVelocity();
+
+	/*Get nodal functions derivatives*/
+	IssmDouble* vdbasis=xNew<IssmDouble>(dim*vnumnodes);
+	IssmDouble* pbasis =xNew<IssmDouble>(pnumnodes);
+	element->NodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);
+	element->NodalFunctionsPressure(pbasis,gauss);
+
+	/*Build B_prime: */
+	if(dim==2){
+		for(int i=0;i<vnumnodes;i++){
+			Bprime[(dim*vnumnodes+pnumnodes)*0+dim*i+0] = vdbasis[0*vnumnodes+i];
+			Bprime[(dim*vnumnodes+pnumnodes)*0+dim*i+1] = 0.;
+			Bprime[(dim*vnumnodes+pnumnodes)*1+dim*i+0] = 0.;
+			Bprime[(dim*vnumnodes+pnumnodes)*1+dim*i+1] = vdbasis[1*vnumnodes+i];
+			Bprime[(dim*vnumnodes+pnumnodes)*2+dim*i+0] = vdbasis[1*vnumnodes+i];
+			Bprime[(dim*vnumnodes+pnumnodes)*2+dim*i+1] = vdbasis[0*vnumnodes+i];
+			Bprime[(dim*vnumnodes+pnumnodes)*3+dim*i+0] = vdbasis[0*vnumnodes+i];
+			Bprime[(dim*vnumnodes+pnumnodes)*3+dim*i+1] = vdbasis[1*vnumnodes+i];
+			Bprime[(dim*vnumnodes+pnumnodes)*4+dim*i+0] = 0.;
+			Bprime[(dim*vnumnodes+pnumnodes)*4+dim*i+1] = 0.;
+		}
+		for(int i=0;i<pnumnodes;i++){
+			Bprime[(dim*vnumnodes+pnumnodes)*0+(dim*vnumnodes)+i] = 0.;
+			Bprime[(dim*vnumnodes+pnumnodes)*1+(dim*vnumnodes)+i] = 0.;
+			Bprime[(dim*vnumnodes+pnumnodes)*2+(dim*vnumnodes)+i] = 0.;
+			Bprime[(dim*vnumnodes+pnumnodes)*3+(dim*vnumnodes)+i] = 0.;
+			Bprime[(dim*vnumnodes+pnumnodes)*4+(dim*vnumnodes)+i] = pbasis[i];
+		}
+	}
+	else{
+		for(int i=0;i<vnumnodes;i++){
+			Bprime[(dim*vnumnodes+pnumnodes)*0+dim*i+0] = vdbasis[0*vnumnodes+i];
+			Bprime[(dim*vnumnodes+pnumnodes)*0+dim*i+1] = 0.;
+			Bprime[(dim*vnumnodes+pnumnodes)*0+dim*i+2] = 0.;
+			Bprime[(dim*vnumnodes+pnumnodes)*1+dim*i+0] = 0.;
+			Bprime[(dim*vnumnodes+pnumnodes)*1+dim*i+1] = vdbasis[1*vnumnodes+i];
+			Bprime[(dim*vnumnodes+pnumnodes)*1+dim*i+2] = 0.;
+			Bprime[(dim*vnumnodes+pnumnodes)*2+dim*i+0] = 0.;
+			Bprime[(dim*vnumnodes+pnumnodes)*2+dim*i+1] = 0.;
+			Bprime[(dim*vnumnodes+pnumnodes)*2+dim*i+2] = vdbasis[2*vnumnodes+i];
+			Bprime[(dim*vnumnodes+pnumnodes)*3+dim*i+0] = vdbasis[1*vnumnodes+i];
+			Bprime[(dim*vnumnodes+pnumnodes)*3+dim*i+1] = vdbasis[0*vnumnodes+i];
+			Bprime[(dim*vnumnodes+pnumnodes)*3+dim*i+2] = 0.;
+			Bprime[(dim*vnumnodes+pnumnodes)*4+dim*i+0] = vdbasis[2*vnumnodes+i];
+			Bprime[(dim*vnumnodes+pnumnodes)*4+dim*i+1] = 0.;
+			Bprime[(dim*vnumnodes+pnumnodes)*4+dim*i+2] = vdbasis[0*vnumnodes+i];
+			Bprime[(dim*vnumnodes+pnumnodes)*5+dim*i+0] = 0.;
+			Bprime[(dim*vnumnodes+pnumnodes)*5+dim*i+1] = vdbasis[2*vnumnodes+i];
+			Bprime[(dim*vnumnodes+pnumnodes)*5+dim*i+2] = vdbasis[1*vnumnodes+i];
+			Bprime[(dim*vnumnodes+pnumnodes)*6+dim*i+0] = vdbasis[0*vnumnodes+i];
+			Bprime[(dim*vnumnodes+pnumnodes)*6+dim*i+1] = vdbasis[1*vnumnodes+i];
+			Bprime[(dim*vnumnodes+pnumnodes)*6+dim*i+2] = vdbasis[2*vnumnodes+i];
+			Bprime[(dim*vnumnodes+pnumnodes)*7+dim*i+0] = 0.;
+			Bprime[(dim*vnumnodes+pnumnodes)*7+dim*i+1] = 0.;
+			Bprime[(dim*vnumnodes+pnumnodes)*7+dim*i+2] = 0.;
+		}
+		for(int i=0;i<pnumnodes;i++){
+			Bprime[(dim*vnumnodes+pnumnodes)*0+(dim*vnumnodes)+i] = 0.;
+			Bprime[(dim*vnumnodes+pnumnodes)*1+(dim*vnumnodes)+i] = 0.;
+			Bprime[(dim*vnumnodes+pnumnodes)*2+(dim*vnumnodes)+i] = 0.;
+			Bprime[(dim*vnumnodes+pnumnodes)*3+(dim*vnumnodes)+i] = 0.;
+			Bprime[(dim*vnumnodes+pnumnodes)*4+(dim*vnumnodes)+i] = 0.;
+			Bprime[(dim*vnumnodes+pnumnodes)*5+(dim*vnumnodes)+i] = 0.;
+			Bprime[(dim*vnumnodes+pnumnodes)*6+(dim*vnumnodes)+i] = 0.;
+			Bprime[(dim*vnumnodes+pnumnodes)*7+(dim*vnumnodes)+i] = pbasis[i];
+		}
+	}
+
+	/*Clean up*/
+	xDelete<IssmDouble>(vdbasis);
+	xDelete<IssmDouble>(pbasis);
 }/*}}}*/
 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 16830)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 16831)
@@ -22,34 +22,45 @@
 		/*Finite element Analysis*/
 		ElementMatrix* CreateKMatrix(Element* element);
+		ElementVector* CreatePVector(Element* element);
+		void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
+		void GetSolutionFromInputsHoriz(Vector<IssmDouble>* solution,Element* element);
+		void InputUpdateFromSolution(IssmDouble* solution,Element* element);
+
+		/*SSA*/
+		ElementMatrix* CreateKMatrixSSA(Element* element);
+		ElementMatrix* CreateKMatrixSSAViscous(Element* element);
+		ElementMatrix* CreateKMatrixSSAFriction(Element* element);
+		ElementVector* CreatePVectorSSA(Element* element);
+		ElementVector* CreatePVectorSSADrivingStress(Element* element);
+		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 InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element);
+		/*L1L2*/
+		void InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element);
+		/*HO*/
 		ElementMatrix* CreateKMatrixHO(Element* element);
 		ElementMatrix* CreateKMatrixHOViscous(Element* element);
 		ElementMatrix* CreateKMatrixHOFriction(Element* element);
-		ElementMatrix* CreateKMatrixSSA(Element* element);
-		ElementMatrix* CreateKMatrixSSAViscous(Element* element);
-		ElementMatrix* CreateKMatrixSSAFriction(Element* element);
-		ElementVector* CreatePVector(Element* element);
+		ElementVector* CreatePVectorHO(Element* element);
+		ElementVector* CreatePVectorHODrivingStress(Element* element);
+		ElementVector* CreatePVectorHOFront(Element* element);
+		void GetBHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
+		void GetBHOprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
+		void InputUpdateFromSolutionHO(IssmDouble* solution,Element* element);
+		/*FS*/
+		ElementMatrix* CreateKMatrixFS(Element* element);
+		ElementMatrix* CreateKMatrixFSViscous(Element* element);
+		ElementMatrix* CreateKMatrixFSFriction(Element* element);
 		ElementVector* CreatePVectorFS(Element* element);
 		ElementVector* CreatePVectorFSViscous(Element* element);
 		ElementVector* CreatePVectorFSShelf(Element* element);
 		ElementVector* CreatePVectorFSFront(Element* element);
-		ElementVector* CreatePVectorHO(Element* element);
-		ElementVector* CreatePVectorHODrivingStress(Element* element);
-		ElementVector* CreatePVectorHOFront(Element* element);
-		ElementVector* CreatePVectorSSA(Element* element);
-		ElementVector* CreatePVectorSSADrivingStress(Element* element);
-		ElementVector* CreatePVectorSSAFront(Element* element);
-		void GetBHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
-		void GetBHOprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
-		void GetBSSA(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
-		void GetBSSAprime(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
-		void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
+		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 GetSolutionFromInputsFS(Vector<IssmDouble>* solution,Element* element);
-		void GetSolutionFromInputsHoriz(Vector<IssmDouble>* solution,Element* element);
-		void InputUpdateFromSolution(IssmDouble* solution,Element* element);
 		void InputUpdateFromSolutionFS(IssmDouble* solution,Element* element);
-		void InputUpdateFromSolutionHO(IssmDouble* solution,Element* element);
+		/*Coupling*/
 		void InputUpdateFromSolutionHOFS(IssmDouble* solution,Element* element);
-		void InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element);
-		void InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element);
 		void InputUpdateFromSolutionSSAFS(IssmDouble* solution,Element* element);
 		void InputUpdateFromSolutionSSAHO(IssmDouble* solution,Element* element);
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 16830)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 16831)
@@ -57,7 +57,8 @@
 		virtual IssmDouble MinEdgeLength(IssmDouble* xyz_list)=0;
 		virtual void   NodalFunctions(IssmDouble* basis,Gauss* gauss)=0;
+		virtual void   NodalFunctionsVelocity(IssmDouble* basis, Gauss* gauss)=0;
+		virtual void   NodalFunctionsPressure(IssmDouble* basis, Gauss* gauss)=0;
 		virtual void   NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss)=0;
-		virtual void   NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss)=0;
-		virtual void   NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss)=0;
+		virtual void   NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss)=0;
 		virtual void   NormalBase(IssmDouble* normal,IssmDouble* xyz_list)=0;
 		virtual void   NormalTop(IssmDouble* normal,IssmDouble* xyz_list)=0;
@@ -144,4 +145,5 @@
 		virtual IssmDouble StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa)=0;
 		virtual void   ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input)=0;
+		virtual void   ViscosityFS(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input)=0;
 		virtual void   ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0;
 		virtual void   ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0;
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16830)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16831)
@@ -2568,4 +2568,12 @@
 }
 /*}}}*/
+/*FUNCTION Penta::NodalFunctionsDerivativesVelocity{{{*/
+void Penta::NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){
+
+	_assert_(gauss->Enum()==GaussPentaEnum);
+	this->GetNodalFunctionsDerivativesVelocity(dbasis,xyz_list,(GaussPenta*)gauss);
+
+}
+/*}}}*/
 /*FUNCTION Penta::NodalFunctionsVelocity{{{*/
 void Penta::NodalFunctionsVelocity(IssmDouble* basis, Gauss* gauss){
@@ -3215,4 +3223,11 @@
 
 	::TransformStiffnessMatrixCoord(Ke,this->nodes,this->NumberofNodes(),transformenum);
+
+}
+/*}}}*/
+/*FUNCTION Penta::TransformStiffnessMatrixCoord(ElementMatrix* pe,int* transformenum_list){{{*/
+void Penta::TransformStiffnessMatrixCoord(ElementMatrix* Ke,int* transformenum_list){
+
+	::TransformStiffnessMatrixCoord(Ke,this->nodes,this->NumberofNodes(),transformenum_list);
 
 }
@@ -3482,4 +3497,18 @@
 	/*Assign output pointer*/
 	*pphi = phi;
+}
+/*}}}*/
+/*FUNCTION Penta::ViscosityFS{{{*/
+void Penta::ViscosityFS(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){
+
+	/*Intermediaries*/
+	IssmDouble viscosity;
+	IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
+
+	this->GetStrainRate3d(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
+	material->GetViscosity3dFS(&viscosity, &epsilon[0]);
+
+	/*Assign output pointer*/
+	*pviscosity=viscosity;
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 16830)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 16831)
@@ -256,4 +256,5 @@
 		void           NodalFunctions(IssmDouble* basis,Gauss* gauss);
 		void           NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
+		void           NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
 		void           NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss);
 		void           NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss);
@@ -272,7 +273,8 @@
 		void           TransformSolutionCoord(IssmDouble* values,int numnodes,int* transformenum_list);/*Tiling only*/
 		void           TransformStiffnessMatrixCoord(ElementMatrix* Ke,int transformenum);
-		void           TransformStiffnessMatrixCoord(ElementMatrix* Ke,int* transformenum_list){_error_("not implemented yet");};
+		void           TransformStiffnessMatrixCoord(ElementMatrix* Ke,int* transformenum_list);
 		void           TransformStiffnessMatrixCoord(ElementMatrix* Ke,int numnodes,int transformenum){_error_("not implemented yet");};
 		void           TransformStiffnessMatrixCoord(ElementMatrix* Ke,int numnodes,int* transformenum_list){_error_("not implemented yet");};
+		void           ViscosityFS(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
 		void           ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
 		void           ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");};
Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 16830)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 16831)
@@ -120,4 +120,5 @@
 		void        NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss){_error_("not implemented yet");};
 		void        NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
+		void        NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
 		bool        NoIceInElement(){_error_("not implemented yet");};
 		void        NormalBase(IssmDouble* normal,IssmDouble* xyz_list){_error_("not implemented yet");};
@@ -145,4 +146,5 @@
 		ElementMatrix* NewElementMatrix(int approximation_enum){_error_("not implemented yet");};
 		void           ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented yet");};
+		void           ViscosityFS(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented");};
 		void           ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");};
 		void           ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented");};
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16830)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 16831)
@@ -2243,4 +2243,12 @@
 }
 /*}}}*/
+/*FUNCTION Tria::NodalFunctionsDerivativesVelocity{{{*/
+void Tria::NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){
+
+	_assert_(gauss->Enum()==GaussTriaEnum);
+	this->GetNodalFunctionsDerivativesVelocity(dbasis,xyz_list,(GaussTria*)gauss);
+
+}
+/*}}}*/
 /*FUNCTION Tria::NodalFunctionsVelocity{{{*/
 void Tria::NodalFunctionsVelocity(IssmDouble* basis, Gauss* gauss){
@@ -2763,4 +2771,11 @@
 
 	::TransformStiffnessMatrixCoord(Ke,this->nodes,this->NumberofNodes(),transformenum);
+
+}
+/*}}}*/
+/*FUNCTION Tria::TransformStiffnessMatrixCoord(ElementMatrix* pe,int* transformenum_list){{{*/
+void Tria::TransformStiffnessMatrixCoord(ElementMatrix* Ke,int* transformenum_list){
+
+	::TransformStiffnessMatrixCoord(Ke,this->nodes,this->NumberofNodes(),transformenum_list);
 
 }
@@ -2908,4 +2923,18 @@
 	delete gauss;
 
+}
+/*}}}*/
+/*FUNCTION Tria::ViscosityFS{{{*/
+void Tria::ViscosityFS(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){
+
+	/*Intermediaries*/
+	IssmDouble viscosity;
+	IssmDouble epsilon[3];/* epsilon=[exx,eyy,exy];    */
+
+	this->GetStrainRate2d(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
+	material->GetViscosity2dvertical(&viscosity, &epsilon[0]);
+
+	/*Assign output pointer*/
+	*pviscosity=viscosity;
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16830)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16831)
@@ -286,4 +286,5 @@
 		void           NodalFunctions(IssmDouble* basis,Gauss* gauss);
 		void           NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
+		void           NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
 		void           NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss);
 		void           NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss);
@@ -301,8 +302,9 @@
 		void           TransformSolutionCoord(IssmDouble* values,int numnodes,int* transformenum_list){_error_("not implemented yet");};
 		void           TransformStiffnessMatrixCoord(ElementMatrix* Ke,int transformenum);
-		void           TransformStiffnessMatrixCoord(ElementMatrix* Ke,int* transformenum_list){_error_("not implemented yet");};
+		void           TransformStiffnessMatrixCoord(ElementMatrix* Ke,int* transformenum_list);
 		void           TransformStiffnessMatrixCoord(ElementMatrix* Ke,int numnodes,int transformenum){_error_("not implemented yet");};
 		void           TransformStiffnessMatrixCoord(ElementMatrix* Ke,int numnodes,int* transformenum_list){_error_("not implemented yet");};
 		void           ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented yet");};
+		void           ViscosityFS(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
 		void           ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not implemented yet");};
 		void           ViscositySSA(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);