Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16828)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16829)
@@ -818,4 +818,6 @@
 		case SSAApproximationEnum: 
 			return CreateKMatrixSSA(element);
+		case HOApproximationEnum: 
+			return CreateKMatrixHO(element);
 		case NoneApproximationEnum:
 			return NULL;
@@ -823,98 +825,4 @@
 			_error_("Approximation "<<EnumToStringx(approximation)<<" not supported");
 	}
-}/*}}}*/
-ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSA(Element* element){/*{{{*/
-
-	/*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");
-	}
-
-	/*compute all stiffness matrices for this element*/
-	ElementMatrix* Ke1=CreateKMatrixSSAViscous(basalelement);
-	ElementMatrix* Ke2=CreateKMatrixSSAFriction(basalelement);
-	ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
-
-	/*clean-up and return*/
-	if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
-	delete Ke1;
-	delete Ke2;
-	return Ke;
-}/*}}}*/
-ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAViscous(Element* element){/*{{{*/
-
-	/*Intermediaries*/
-	IssmDouble  viscosity,newviscosity,oldviscosity;
-	IssmDouble  viscosity_overshoot,thickness,Jdet;
-	IssmDouble  D_scalar;
-	IssmDouble *xyz_list = 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>(3*numdof);
-	IssmDouble*    Bprime = xNew<IssmDouble>(3*numdof);
-	IssmDouble*    D      = xNewZeroInit<IssmDouble>(3*3);
-
-	/*Retrieve all inputs and parameters*/
-	element->GetVerticesCoordinates(&xyz_list);
-	Input* thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input);
-	Input* vx_input=element->GetInput(VxEnum);               _assert_(vx_input);
-	Input* vy_input=element->GetInput(VyEnum);               _assert_(vy_input);
-	Input* vxold_input=element->GetInput(VxPicardEnum);      _assert_(vxold_input);
-	Input* vyold_input=element->GetInput(VyPicardEnum);      _assert_(vyold_input);
-	element->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum);
-
-	/* Start  looping on the number of gaussian points: */
-	Gauss* gauss = element->NewGauss(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
-
-		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-		this->GetBSSA(B,element,xyz_list,gauss);
-		this->GetBSSAprime(Bprime,element,xyz_list,gauss);
-
-		element->ViscositySSA(&viscosity,xyz_list,gauss,vx_input,vy_input);
-		element->ViscositySSA(&oldviscosity,xyz_list,gauss,vxold_input,vyold_input);
-		thickness_input->GetInputValue(&thickness, gauss);
-
-		newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
-		D_scalar=2.*newviscosity*thickness*gauss->weight*Jdet;
-		for(int i=0;i<3;i++) D[i*3+i]=D_scalar;
-
-		TripleMultiply(B,3,numdof,1,
-					D,3,3,0,
-					Bprime,3,numdof,0,
-					&Ke->values[0],1);
-	}
-
-	/*Transform Coordinate System*/
-	element->TransformStiffnessMatrixCoord(Ke,XYEnum);
-
-	/*Clean up and return*/
-	delete gauss;
-	xDelete<IssmDouble>(xyz_list);
-	xDelete<IssmDouble>(D);
-	xDelete<IssmDouble>(Bprime);
-	xDelete<IssmDouble>(B);
-	return Ke;
-}/*}}}*/
-ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAFriction(Element* element){/*{{{*/
-	return NULL;
 }/*}}}*/
 ElementVector* StressbalanceAnalysis::CreatePVector(Element* element){/*{{{*/
@@ -934,365 +842,4 @@
 			_error_("Approximation "<<EnumToStringx(approximation)<<" not supported");
 	}
-}/*}}}*/
-ElementVector* StressbalanceAnalysis::CreatePVectorFS(Element* element){/*{{{*/
-
-	/*compute all load vectors for this element*/
-	ElementVector* pe1=CreatePVectorFSViscous(element);
-	ElementVector* pe2=CreatePVectorFSShelf(element);
-	ElementVector* pe3=CreatePVectorFSFront(element);
-	ElementVector* pe =new ElementVector(pe1,pe2,pe3);
-
-	/*clean-up and return*/
-	delete pe1;
-	delete pe2;
-	delete pe3;
-	return pe;
-}/*}}}*/
-ElementVector* StressbalanceAnalysis::CreatePVectorFSViscous(Element* element){/*{{{*/
-
-	int         i,meshtype,dim;
-	IssmDouble  Jdet,forcex,forcey,forcez;
-	IssmDouble *xyz_list = NULL;
-
-	/*Get problem dimension*/
-	element->FindParam(&meshtype,MeshTypeEnum);
-	switch(meshtype){
-		case Mesh2DverticalEnum: dim = 2; break;
-		case Mesh3DEnum:         dim = 3; 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();
-
-	/*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 vectors*/
-	ElementVector* pe     = element->NewElementVector(FSvelocityEnum);
-	IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
-
-	/*Retrieve all inputs and parameters*/
-	element->GetVerticesCoordinates(&xyz_list);
-	Input*      loadingforcex_input=element->GetInput(LoadingforceXEnum);  _assert_(loadingforcex_input);
-	Input*      loadingforcey_input=element->GetInput(LoadingforceYEnum);  _assert_(loadingforcey_input);
-	Input*      loadingforcez_input=element->GetInput(LoadingforceZEnum);  _assert_(loadingforcez_input);
-	IssmDouble  rho_ice =element->GetMaterialParameter(MaterialsRhoIceEnum);
-	IssmDouble  gravity =element->GetMaterialParameter(ConstantsGEnum);
-
-	/* 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->NodalFunctionsVelocity(vbasis,gauss);
-
-		loadingforcex_input->GetInputValue(&forcex,gauss);
-		loadingforcey_input->GetInputValue(&forcey,gauss);
-		if(dim==3) loadingforcez_input->GetInputValue(&forcez,gauss);
-
-		for(i=0;i<vnumnodes;i++){
-			pe->values[i*dim+0] += +rho_ice*forcex *Jdet*gauss->weight*vbasis[i];
-			pe->values[i*dim+1] += +rho_ice*forcey *Jdet*gauss->weight*vbasis[i];
-			if(dim==3){
-				pe->values[i*dim+2] += +rho_ice*forcez*Jdet*gauss->weight*vbasis[i];
-				pe->values[i*dim+2] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i];
-			}
-			else{
-				pe->values[i*dim+1] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i];
-			}
-		}
-	}
-
-	/*Transform coordinate system*/
-	element->TransformLoadVectorCoord(pe,cs_list);
-
-	/*Clean up and return*/
-	delete gauss;
-	xDelete<int>(cs_list);
-	xDelete<IssmDouble>(vbasis);
-	xDelete<IssmDouble>(xyz_list);
-	return pe;
-}/*}}}*/
-ElementVector* StressbalanceAnalysis::CreatePVectorFSShelf(Element* element){/*{{{*/
-
-	int         i,meshtype,dim;
-	IssmDouble  Jdet,water_pressure,bed;
-	IssmDouble	normal[3];
-	IssmDouble *xyz_list_base = NULL;
-
-	/*Get basal element*/
-	if(!element->IsOnBed() || !element->IsFloating()) return NULL;
-
-	/*Get problem dimension*/
-	element->FindParam(&meshtype,MeshTypeEnum);
-	switch(meshtype){
-		case Mesh2DverticalEnum: dim = 2; break;
-		case Mesh3DEnum:         dim = 3; 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();
-
-	/*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 vectors*/
-	ElementVector* pe     = element->NewElementVector(FSvelocityEnum);
-	IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
-
-	/*Retrieve all inputs and parameters*/
-	element->GetVerticesCoordinatesBase(&xyz_list_base);
-	Input*      bed_input=element->GetInput(BedEnum); _assert_(bed_input);
-	IssmDouble  rho_water=element->GetMaterialParameter(MaterialsRhoWaterEnum);
-	IssmDouble  gravity  =element->GetMaterialParameter(ConstantsGEnum);
-
-	/* Start  looping on the number of gaussian points: */
-	Gauss* gauss=element->NewGaussBase(5);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
-
-		element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
-		element->NodalFunctionsVelocity(vbasis,gauss);
-
-		element->NormalBase(&normal[0],xyz_list_base);
-		_assert_(normal[dim-1]<0.);
-		bed_input->GetInputValue(&bed, gauss);
-		water_pressure=gravity*rho_water*bed;
-
-		for(i=0;i<vnumnodes;i++){
-			pe->values[i*dim+0] += water_pressure*gauss->weight*Jdet*vbasis[i]*normal[0];
-			pe->values[i*dim+1] += water_pressure*gauss->weight*Jdet*vbasis[i]*normal[1];
-			if(dim==3){
-				pe->values[i*dim+2]+=water_pressure*gauss->weight*Jdet*vbasis[i]*normal[2];
-			}
-		}
-	}
-
-	/*Transform coordinate system*/
-	element->TransformLoadVectorCoord(pe,cs_list);
-
-	/*Clean up and return*/
-	delete gauss;
-	xDelete<int>(cs_list);
-	xDelete<IssmDouble>(vbasis);
-	xDelete<IssmDouble>(xyz_list_base);
-	return pe;
-}/*}}}*/
-ElementVector* StressbalanceAnalysis::CreatePVectorFSFront(Element* element){/*{{{*/
-
-	return NULL;
-}/*}}}*/
-ElementVector* StressbalanceAnalysis::CreatePVectorHO(Element* element){/*{{{*/
-
-	/*compute all load vectors for this element*/
-	ElementVector* pe1=CreatePVectorHODrivingStress(element);
-	ElementVector* pe2=CreatePVectorHOFront(element);
-	ElementVector* pe =new ElementVector(pe1,pe2);
-
-	/*clean-up and return*/
-	delete pe1;
-	delete pe2;
-	return pe;
-}/*}}}*/
-ElementVector* StressbalanceAnalysis::CreatePVectorHODrivingStress(Element* element){/*{{{*/
-
-	/*Intermediaries */
-	IssmDouble  Jdet,slope[3];
-	IssmDouble* xyz_list = NULL;
-
-	/*Fetch number of nodes and dof for this finite element*/
-	int numnodes = element->GetNumberOfNodes();
-
-	/*Initialize Element vector and vectors*/
-	ElementVector* pe=element->NewElementVector(HOApproximationEnum);
-	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
-
-	/*Retrieve all inputs and parameters*/
-	element->GetVerticesCoordinates(&xyz_list);
-	Input*     surface_input = element->GetInput(SurfaceEnum);   _assert_(surface_input);
-	IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum);
-
-	/* Start  looping on the number of gaussian points: */
-	Gauss* gauss=element->NewGauss(3);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
-
-		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-		element->NodalFunctions(basis, gauss);
-		surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
-
-		for(int i=0;i<numnodes;i++){
-			pe->values[i*2+0]+=-rhog*slope[0]*Jdet*gauss->weight*basis[i];
-			pe->values[i*2+1]+=-rhog*slope[1]*Jdet*gauss->weight*basis[i];
-		}
-	}
-
-	/*Transform coordinate system*/
-	element->TransformLoadVectorCoord(pe,XYEnum);
-
-	/*Clean up and return*/
-	xDelete<IssmDouble>(basis);
-	xDelete<IssmDouble>(xyz_list);
-	delete gauss;
-	return pe;
-}/*}}}*/
-ElementVector* StressbalanceAnalysis::CreatePVectorHOFront(Element* element){/*{{{*/
-
-	return NULL;
-}/*}}}*/
-ElementVector* StressbalanceAnalysis::CreatePVectorSSA(Element* element){/*{{{*/
-
-	/*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");
-	}
-
-	/*compute all load vectors for this element*/
-	ElementVector* pe1=CreatePVectorSSADrivingStress(basalelement);
-	ElementVector* pe2=CreatePVectorSSAFront(basalelement);
-	ElementVector* pe =new ElementVector(pe1,pe2);
-
-	/*clean-up and return*/
-	if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
-	delete pe1;
-	delete pe2;
-	return pe;
-}/*}}}*/
-ElementVector* StressbalanceAnalysis::CreatePVectorSSADrivingStress(Element* element){/*{{{*/
-
-	/*Intermediaries */
-	IssmDouble  thickness,Jdet,slope[2];
-	IssmDouble* xyz_list = NULL;
-
-	/*Fetch number of nodes and dof for this finite element*/
-	int numnodes = element->GetNumberOfNodes();
-
-	/*Initialize Element vector and vectors*/
-	ElementVector* pe    = element->NewElementVector(SSAApproximationEnum);
-	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
-
-	/*Retrieve all inputs and parameters*/
-	element->GetVerticesCoordinates(&xyz_list);
-	Input*     thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input); 
-	Input*     surface_input  =element->GetInput(SurfaceEnum);   _assert_(surface_input);
-	IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum);
-
-	/* Start  looping on the number of gaussian points: */
-	Gauss* gauss=element->NewGauss(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
-
-		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
-		element->NodalFunctions(basis, gauss);
-
-		thickness_input->GetInputValue(&thickness,gauss);
-		surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
-
-		for(int i=0;i<numnodes;i++){
-			pe->values[i*2+0]+=-rhog*thickness*slope[0]*Jdet*gauss->weight*basis[i];
-			pe->values[i*2+1]+=-rhog*thickness*slope[1]*Jdet*gauss->weight*basis[i];
-		}
-	}
-
-	/*Transform coordinate system*/
-	element->TransformLoadVectorCoord(pe,XYEnum);
-
-	/*Clean up and return*/
-	xDelete<IssmDouble>(xyz_list);
-	xDelete<IssmDouble>(basis);
-	delete gauss;
-	return pe;
-}/*}}}*/
-ElementVector* StressbalanceAnalysis::CreatePVectorSSAFront(Element* element){/*{{{*/
-
-	return NULL;
-}/*}}}*/
-void StressbalanceAnalysis::GetBSSA(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
-	 * For node i, Bi can be expressed in the actual coordinate system
-	 * by: 
-	 *       Bi=[ dN/dx           0    ]
-	 *          [   0           dN/dy  ]
-	 *          [ 1/2*dN/dy  1/2*dN/dx ]
-	 * where N is the finiteelement function for node i.
-	 *
-	 * We assume B has been allocated already, of size: 3x(NDOF2*numnodes)
-	 */
-
-	/*Fetch number of nodes for this finite element*/
-	int numnodes = element->GetNumberOfNodes();
-
-	/*Get nodal functions derivatives*/
-	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
-	element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-
-	/*Build B: */
-	for(int i=0;i<numnodes;i++){
-		B[2*numnodes*0+2*i+0] = dbasis[0*numnodes+i];
-		B[2*numnodes*0+2*i+1] = 0.;
-		B[2*numnodes*1+2*i+0] = 0.;
-		B[2*numnodes*1+2*i+1] = dbasis[1*numnodes+i];
-		B[2*numnodes*2+2*i+0] = .5*dbasis[1*numnodes+i];
-		B[2*numnodes*2+2*i+1] = .5*dbasis[0*numnodes+i];
-	}
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(dbasis);
-}/*}}}*/
-void StressbalanceAnalysis::GetBSSAprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
-	 * For node i, Bi' can be expressed in the actual coordinate system
-	 * by: 
-	 *       Bi_prime=[ 2*dN/dx    dN/dy ]
-	 *                [   dN/dx  2*dN/dy ]
-	 *                [   dN/dy    dN/dx ]
-	 * where hNis the finiteelement function for node i.
-	 *
-	 * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
-	 */
-
-	/*Fetch number of nodes for this finite element*/
-	int numnodes = element->GetNumberOfNodes();
-
-	/*Get nodal functions derivatives*/
-	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
-	element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-
-	/*Build B': */
-	for(int i=0;i<numnodes;i++){
-		Bprime[2*numnodes*0+2*i+0] = 2.*dbasis[0*numnodes+i];
-		Bprime[2*numnodes*0+2*i+1] =    dbasis[1*numnodes+i];
-		Bprime[2*numnodes*1+2*i+0] =    dbasis[0*numnodes+i];
-		Bprime[2*numnodes*1+2*i+1] = 2.*dbasis[1*numnodes+i];
-		Bprime[2*numnodes*2+2*i+0] =    dbasis[1*numnodes+i];
-		Bprime[2*numnodes*2+2*i+1] =    dbasis[0*numnodes+i];
-	}
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(dbasis);
 }/*}}}*/
 void StressbalanceAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
@@ -1317,72 +864,4 @@
 	}
 }/*}}}*/
-void StressbalanceAnalysis::GetSolutionFromInputsFS(Vector<IssmDouble>* solution,Element* element){/*{{{*/
-
-	int*         vdoflist=NULL;
-	int*         pdoflist=NULL;
-	Input*       vz_input=NULL;
-	int          meshtype,dim;
-	IssmDouble   vx,vy,vz,p;
-	IssmDouble   FSreconditioning;
-
-	/*Get some parameters*/
-	element->FindParam(&meshtype,MeshTypeEnum);
-	element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
-	switch(meshtype){
-		case Mesh2DverticalEnum: dim = 2; break;
-		case Mesh3DEnum:         dim = 3; break;
-		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
-	}
-
-	/*Fetch number of nodes and dof for this finite element*/
-	int vnumnodes = element->NumberofNodesVelocity();
-	int pnumnodes = element->NumberofNodesPressure();
-	int vnumdof   = vnumnodes*dim;
-	int pnumdof   = pnumnodes*1;
-
-	/*Initialize values*/
-	IssmDouble* vvalues = xNew<IssmDouble>(vnumdof);
-	IssmDouble* pvalues = xNew<IssmDouble>(pnumdof);
-
-	/*Get dof list: */
-	element->GetDofListVelocity(&vdoflist,GsetEnum);
-	element->GetDofListPressure(&pdoflist,GsetEnum);
-	Input*     vx_input=element->GetInput(VxEnum);       _assert_(vx_input);
-	Input*     vy_input=element->GetInput(VyEnum);       _assert_(vy_input);
-	if(dim==3){vz_input=element->GetInput(VzEnum);       _assert_(vz_input);}
-	Input*     p_input =element->GetInput(PressureEnum); _assert_(p_input);
-
-	element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
-
-	/*Ok, we have the velocities in inputs, fill in solution */
-	Gauss* gauss = element->NewGauss();
-	for(int i=0;i<vnumnodes;i++){
-		gauss->GaussNode(element->VelocityInterpolation(),i);
-		vx_input->GetInputValue(&vx,gauss);
-		vy_input->GetInputValue(&vy,gauss);
-		vvalues[i*dim+0]=vx;
-		vvalues[i*dim+1]=vy;
-		if(dim==3){
-			vz_input->GetInputValue(&vz,gauss);
-			vvalues[i*dim+2]=vz;
-		}
-	}
-	for(int i=0;i<pnumnodes;i++){
-		gauss->GaussNode(element->PressureInterpolation(),i);
-		p_input->GetInputValue(&p ,gauss);
-		pvalues[i]=p/FSreconditioning;
-	}
-
-	/*Add value to global vector*/
-	solution->SetValues(vnumdof,vdoflist,vvalues,INS_VAL);
-	solution->SetValues(pnumdof,pdoflist,pvalues,INS_VAL);
-
-	/*Free ressources:*/
-	delete gauss;
-	xDelete<int>(pdoflist);
-	xDelete<int>(vdoflist);
-	xDelete<IssmDouble>(pvalues);
-	xDelete<IssmDouble>(vvalues);
-}/*}}}*/
 void StressbalanceAnalysis::GetSolutionFromInputsHoriz(Vector<IssmDouble>* solution,Element* element){/*{{{*/
 
@@ -1412,6 +891,6 @@
 		vx_input->GetInputValue(&vx,gauss);
 		vy_input->GetInputValue(&vy,gauss);
-		values[i*NDOF2+0]=vx;
-		values[i*NDOF2+1]=vy;
+		values[i*2+0]=vx;
+		values[i*2+1]=vy;
 	}
 
@@ -1455,88 +934,334 @@
 	}
 }/*}}}*/
-void StressbalanceAnalysis::InputUpdateFromSolutionFS(IssmDouble* solution,Element* element){/*{{{*/
-
-	int          i,dim,meshtype;
-	int*         vdoflist=NULL;
-	int*         pdoflist=NULL;
-	IssmDouble   FSreconditioning;
-
+
+/*SSA*/
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSA(Element* element){/*{{{*/
+
+	/*Intermediaries*/
+	int      meshtype;
+	Element* basalelement;
+
+	/*Get basal element*/
 	element->FindParam(&meshtype,MeshTypeEnum);
-	element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
 	switch(meshtype){
-		case Mesh2DverticalEnum: dim = 2; break;
-		case Mesh3DEnum:         dim = 3; break;
+		case Mesh2DhorizontalEnum:
+			basalelement = element;
+			break;
+		case Mesh3DEnum:
+			if(!element->IsOnBed()) return NULL;
+			basalelement = element->SpawnBasalElement();
+			break;
 		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
 	}
 
+	/*compute all stiffness matrices for this element*/
+	ElementMatrix* Ke1=CreateKMatrixSSAViscous(basalelement);
+	ElementMatrix* Ke2=CreateKMatrixSSAFriction(basalelement);
+	ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
+
+	/*clean-up and return*/
+	if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
+	delete Ke1;
+	delete Ke2;
+	return Ke;
+}/*}}}*/
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAViscous(Element* element){/*{{{*/
+
+	/*Intermediaries*/
+	IssmDouble  viscosity,newviscosity,oldviscosity;
+	IssmDouble  viscosity_overshoot,thickness,Jdet;
+	IssmDouble  D_scalar;
+	IssmDouble *xyz_list = NULL;
+
 	/*Fetch number of nodes and dof for this finite element*/
-	int vnumnodes = element->GetNumberOfNodesVelocity();
-	int pnumnodes = element->GetNumberOfNodesPressure();
-	int vnumdof   = vnumnodes*dim;
-	int pnumdof   = pnumnodes*1;
-
-	/*Initialize values*/
-	IssmDouble* values   = xNew<IssmDouble>(vnumdof+pnumdof);
-	IssmDouble* vx       = xNew<IssmDouble>(vnumnodes);
-	IssmDouble* vy       = xNew<IssmDouble>(vnumnodes);
-	IssmDouble* vz       = xNew<IssmDouble>(vnumnodes);
-	IssmDouble* vel      = xNew<IssmDouble>(vnumnodes);
-	IssmDouble* pressure = xNew<IssmDouble>(pnumnodes);
-
-	/*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;
-
-	/*Get dof list: */
-	element->GetDofListVelocity(&vdoflist,GsetEnum);
-	element->GetDofListPressure(&pdoflist,GsetEnum);
+	int numnodes = element->GetNumberOfNodes();
+	int numdof   = numnodes*2;
+
+	/*Initialize Element matrix and vectors*/
+	ElementMatrix* Ke     = element->NewElementMatrix(SSAApproximationEnum);
+	IssmDouble*    B      = xNew<IssmDouble>(3*numdof);
+	IssmDouble*    Bprime = xNew<IssmDouble>(3*numdof);
+	IssmDouble*    D      = xNewZeroInit<IssmDouble>(3*3);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+	Input* thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input);
+	Input* vx_input=element->GetInput(VxEnum);               _assert_(vx_input);
+	Input* vy_input=element->GetInput(VyEnum);               _assert_(vy_input);
+	Input* vxold_input=element->GetInput(VxPicardEnum);      _assert_(vxold_input);
+	Input* vyold_input=element->GetInput(VyPicardEnum);      _assert_(vyold_input);
+	element->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum);
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss = element->NewGauss(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		this->GetBSSA(B,element,xyz_list,gauss);
+		this->GetBSSAprime(Bprime,element,xyz_list,gauss);
+
+		element->ViscositySSA(&viscosity,xyz_list,gauss,vx_input,vy_input);
+		element->ViscositySSA(&oldviscosity,xyz_list,gauss,vxold_input,vyold_input);
+		thickness_input->GetInputValue(&thickness, gauss);
+
+		newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
+		D_scalar=2.*newviscosity*thickness*gauss->weight*Jdet;
+		for(int i=0;i<3;i++) D[i*3+i]=D_scalar;
+
+		TripleMultiply(B,3,numdof,1,
+					D,3,3,0,
+					Bprime,3,numdof,0,
+					&Ke->values[0],1);
+	}
+
+	/*Transform Coordinate System*/
+	element->TransformStiffnessMatrixCoord(Ke,XYEnum);
+
+	/*Clean up and return*/
+	delete gauss;
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(D);
+	xDelete<IssmDouble>(Bprime);
+	xDelete<IssmDouble>(B);
+	return Ke;
+}/*}}}*/
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAFriction(Element* element){/*{{{*/
+	return NULL;
+}/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorSSA(Element* element){/*{{{*/
+
+	/*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");
+	}
+
+	/*compute all load vectors for this element*/
+	ElementVector* pe1=CreatePVectorSSADrivingStress(basalelement);
+	ElementVector* pe2=CreatePVectorSSAFront(basalelement);
+	ElementVector* pe =new ElementVector(pe1,pe2);
+
+	/*clean-up and return*/
+	if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
+	delete pe1;
+	delete pe2;
+	return pe;
+}/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorSSADrivingStress(Element* element){/*{{{*/
+
+	/*Intermediaries */
+	IssmDouble  thickness,Jdet,slope[2];
+	IssmDouble* xyz_list = NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Initialize Element vector and vectors*/
+	ElementVector* pe    = element->NewElementVector(SSAApproximationEnum);
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+	Input*     thickness_input=element->GetInput(ThicknessEnum); _assert_(thickness_input); 
+	Input*     surface_input  =element->GetInput(SurfaceEnum);   _assert_(surface_input);
+	IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum);
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=element->NewGauss(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		element->NodalFunctions(basis, gauss);
+
+		thickness_input->GetInputValue(&thickness,gauss);
+		surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
+
+		for(int i=0;i<numnodes;i++){
+			pe->values[i*2+0]+=-rhog*thickness*slope[0]*Jdet*gauss->weight*basis[i];
+			pe->values[i*2+1]+=-rhog*thickness*slope[1]*Jdet*gauss->weight*basis[i];
+		}
+	}
+
+	/*Transform coordinate system*/
+	element->TransformLoadVectorCoord(pe,XYEnum);
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(basis);
+	delete gauss;
+	return pe;
+}/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorSSAFront(Element* element){/*{{{*/
+
+	return NULL;
+}/*}}}*/
+void StressbalanceAnalysis::GetBSSA(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
+	 * For node i, Bi can be expressed in the actual coordinate system
+	 * by: 
+	 *       Bi=[ dN/dx           0    ]
+	 *          [   0           dN/dy  ]
+	 *          [ 1/2*dN/dy  1/2*dN/dx ]
+	 * where N is the finiteelement function for node i.
+	 *
+	 * We assume B has been allocated already, of size: 3x(NDOF2*numnodes)
+	 */
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Get nodal functions derivatives*/
+	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
+	element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+	/*Build B: */
+	for(int i=0;i<numnodes;i++){
+		B[2*numnodes*0+2*i+0] = dbasis[0*numnodes+i];
+		B[2*numnodes*0+2*i+1] = 0.;
+		B[2*numnodes*1+2*i+0] = 0.;
+		B[2*numnodes*1+2*i+1] = dbasis[1*numnodes+i];
+		B[2*numnodes*2+2*i+0] = .5*dbasis[1*numnodes+i];
+		B[2*numnodes*2+2*i+1] = .5*dbasis[0*numnodes+i];
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(dbasis);
+}/*}}}*/
+void StressbalanceAnalysis::GetBSSAprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
+	 * For node i, Bi' can be expressed in the actual coordinate system
+	 * by: 
+	 *       Bi_prime=[ 2*dN/dx    dN/dy ]
+	 *                [   dN/dx  2*dN/dy ]
+	 *                [   dN/dy    dN/dx ]
+	 * where hNis the finiteelement function for node i.
+	 *
+	 * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
+	 */
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Get nodal functions derivatives*/
+	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
+	element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+	/*Build B': */
+	for(int i=0;i<numnodes;i++){
+		Bprime[2*numnodes*0+2*i+0] = 2.*dbasis[0*numnodes+i];
+		Bprime[2*numnodes*0+2*i+1] =    dbasis[1*numnodes+i];
+		Bprime[2*numnodes*1+2*i+0] =    dbasis[0*numnodes+i];
+		Bprime[2*numnodes*1+2*i+1] = 2.*dbasis[1*numnodes+i];
+		Bprime[2*numnodes*2+2*i+0] =    dbasis[1*numnodes+i];
+		Bprime[2*numnodes*2+2*i+1] =    dbasis[0*numnodes+i];
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(dbasis);
+}/*}}}*/
+void StressbalanceAnalysis::InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element){/*{{{*/
+
+	int         i,meshtype;
+	IssmDouble  rho_ice,g;
+	int*        doflist=NULL;
+	IssmDouble* xyz_list=NULL;
+	Element*    basalelement=NULL;
+
+	/*Deal with pressure first*/
+	int numvertices = element->GetNumberOfVertices();
+	IssmDouble* pressure  = xNew<IssmDouble>(numvertices);
+	IssmDouble* thickness = xNew<IssmDouble>(numvertices);
+	IssmDouble* surface   = xNew<IssmDouble>(numvertices);
+
+	element->FindParam(&meshtype,MeshTypeEnum);
+	rho_ice =element->GetMaterialParameter(MaterialsRhoIceEnum);
+	g       =element->GetMaterialParameter(ConstantsGEnum);
+	switch(meshtype){
+		case Mesh2DhorizontalEnum:
+			element->GetInputListOnVertices(thickness,ThicknessEnum);
+			for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*thickness[i];
+			break;
+		case Mesh3DEnum:
+			element->GetVerticesCoordinates(&xyz_list);
+			element->GetInputListOnVertices(surface,SurfaceEnum);
+			for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]);
+			break;
+		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+	}
+	element->AddInput(PressureEnum,pressure,P1Enum);
+	xDelete<IssmDouble>(pressure);
+	xDelete<IssmDouble>(thickness);
+	xDelete<IssmDouble>(surface);
+
+	/*Get basal element*/
+	switch(meshtype){
+		case Mesh2DhorizontalEnum:
+			basalelement = element;
+			break;
+		case Mesh3DEnum:
+			if(!element->IsOnBed()){xDelete<IssmDouble>(xyz_list); return;}
+			basalelement=element->SpawnBasalElement();
+			break;
+		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+	}
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = basalelement->GetNumberOfNodes();
+	int numdof   = numnodes*2;
+
+	/*Fetch dof list and allocate solution vectors*/
+	basalelement->GetDofList(&doflist,SSAApproximationEnum,GsetEnum);
+	IssmDouble* values    = xNew<IssmDouble>(numdof);
+	IssmDouble* vx        = xNew<IssmDouble>(numnodes);
+	IssmDouble* vy        = xNew<IssmDouble>(numnodes);
+	IssmDouble* vz        = xNew<IssmDouble>(numnodes);
+	IssmDouble* vel       = xNew<IssmDouble>(numnodes);
 
 	/*Use the dof list to index into the solution vector: */
-	for(i=0;i<vnumdof;i++) values[i]        =solution[vdoflist[i]];
-	for(i=0;i<pnumdof;i++) values[vnumdof+i]=solution[pdoflist[i]];
+	for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
 
 	/*Transform solution in Cartesian Space*/
-	element->TransformSolutionCoord(values,cs_list);
-
-	/*Ok, we have vx and vy in values, fill in all arrays: */
-	for(i=0;i<vnumnodes;i++){
-		vx[i] = values[i*dim+0];
-		vy[i] = values[i*dim+1];
+	basalelement->TransformSolutionCoord(&values[0],XYEnum);
+	basalelement->FindParam(&meshtype,MeshTypeEnum);
+
+	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
+	for(i=0;i<numnodes;i++){
+		vx[i]=values[i*2+0];
+		vy[i]=values[i*2+1];
+
+		/*Check solution*/
 		if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
 		if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
-
-		if(dim==3){
-			vz[i] = values[i*dim+2];
-			if(xIsNan<IssmDouble>(vz[i])) _error_("NaN found in solution vector");
-		}
-	}
-	for(i=0;i<pnumnodes;i++){
-		pressure[i] = values[vnumdof+i];
-		if(xIsNan<IssmDouble>(pressure[i])) _error_("NaN found in solution vector");
-	}
-
-	/*Recondition pressure and compute vel: */
-	for(i=0;i<pnumnodes;i++) pressure[i] = pressure[i]*FSreconditioning;
-	if(dim==3) for(i=0;i<vnumnodes;i++) vel[i] = sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
-	else       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 
+	}
+
+	/*Get Vz and compute vel*/
+	basalelement->GetInputListOnNodes(&vz[0],VzEnum,0.);
+	for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
+
+	/*Now, we have to move the previous Vx and Vy inputs  to old 
 	 * status, otherwise, we'll wipe them off: */
 	element->InputChangeName(VxEnum,VxPicardEnum);
 	element->InputChangeName(VyEnum,VyPicardEnum);
-	element->InputChangeName(PressureEnum,PressurePicardEnum);
-	if(dim==3) element->InputChangeName(VzEnum,VzPicardEnum);
 
 	/*Add vx and vy as inputs to the tria element: */
-	element->AddInput(VxEnum,vx,P1Enum);
-	element->AddInput(VyEnum,vy,P1Enum);
-	element->AddInput(VelEnum,vel,P1Enum);
-	element->AddInput(PressureEnum,pressure,P1Enum);
-	if(dim==3) element->AddInput(VzEnum,vz,P1Enum);
+	element->AddBasalInput(VxEnum,vx,P1Enum);
+	element->AddBasalInput(VyEnum,vy,P1Enum);
+	element->AddBasalInput(VelEnum,vel,P1Enum);
 
 	/*Free ressources:*/
-	xDelete<IssmDouble>(pressure);
 	xDelete<IssmDouble>(vel);
 	xDelete<IssmDouble>(vz);
@@ -1544,7 +1269,319 @@
 	xDelete<IssmDouble>(vx);
 	xDelete<IssmDouble>(values);
-	xDelete<int>(vdoflist);
-	xDelete<int>(pdoflist);
-	xDelete<int>(cs_list);
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<int>(doflist);
+	if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
+}/*}}}*/
+
+/*L1L2*/
+void StressbalanceAnalysis::InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element){/*{{{*/
+
+	int         i,meshtype;
+	IssmDouble  rho_ice,g;
+	int*        doflist=NULL;
+	IssmDouble* xyz_list=NULL;
+	Element*    basalelement=NULL;
+
+	/*Deal with pressure first*/
+	int numvertices = element->GetNumberOfVertices();
+	IssmDouble* pressure  = xNew<IssmDouble>(numvertices);
+	IssmDouble* thickness = xNew<IssmDouble>(numvertices);
+	IssmDouble* surface   = xNew<IssmDouble>(numvertices);
+
+	element->FindParam(&meshtype,MeshTypeEnum);
+	rho_ice =element->GetMaterialParameter(MaterialsRhoIceEnum);
+	g       =element->GetMaterialParameter(ConstantsGEnum);
+	switch(meshtype){
+		case Mesh2DhorizontalEnum:
+			element->GetInputListOnVertices(thickness,ThicknessEnum);
+			for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*thickness[i];
+			break;
+		case Mesh3DEnum:
+			element->GetVerticesCoordinates(&xyz_list);
+			element->GetInputListOnVertices(surface,SurfaceEnum);
+			for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]);
+			break;
+		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+	}
+	element->AddInput(PressureEnum,pressure,P1Enum);
+	xDelete<IssmDouble>(pressure);
+	xDelete<IssmDouble>(thickness);
+	xDelete<IssmDouble>(surface);
+
+	/*Get basal element*/
+	switch(meshtype){
+		case Mesh2DhorizontalEnum:
+			basalelement = element;
+			break;
+		case Mesh3DEnum:
+			if(!element->IsOnBed()){xDelete<IssmDouble>(xyz_list); return;}
+			basalelement=element->SpawnBasalElement();
+			break;
+		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+	}
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = basalelement->GetNumberOfNodes();
+	int numdof   = numnodes*2;
+
+	/*Fetch dof list and allocate solution vectors*/
+	basalelement->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
+	IssmDouble* values    = xNew<IssmDouble>(numdof);
+	IssmDouble* vx        = xNew<IssmDouble>(numnodes);
+	IssmDouble* vy        = xNew<IssmDouble>(numnodes);
+	IssmDouble* vz        = xNew<IssmDouble>(numnodes);
+	IssmDouble* vel       = xNew<IssmDouble>(numnodes);
+
+	/*Use the dof list to index into the solution vector: */
+	for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
+
+	/*Transform solution in Cartesian Space*/
+	basalelement->TransformSolutionCoord(&values[0],XYEnum);
+	basalelement->FindParam(&meshtype,MeshTypeEnum);
+
+	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
+	for(i=0;i<numnodes;i++){
+		vx[i]=values[i*2+0];
+		vy[i]=values[i*2+1];
+
+		/*Check solution*/
+		if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
+	}
+
+	/*Get Vz and compute vel*/
+	basalelement->GetInputListOnNodes(&vz[0],VzEnum,0.);
+	for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
+
+	/*Now, we have to move the previous Vx and Vy inputs  to old 
+	 * status, otherwise, we'll wipe them off: */
+	element->InputChangeName(VxEnum,VxPicardEnum);
+	element->InputChangeName(VyEnum,VyPicardEnum);
+
+	/*Add vx and vy as inputs to the tria element: */
+	element->AddBasalInput(VxEnum,vx,P1Enum);
+	element->AddBasalInput(VyEnum,vy,P1Enum);
+	element->AddBasalInput(VelEnum,vel,P1Enum);
+
+	/*Free ressources:*/
+	xDelete<IssmDouble>(vel);
+	xDelete<IssmDouble>(vz);
+	xDelete<IssmDouble>(vy);
+	xDelete<IssmDouble>(vx);
+	xDelete<IssmDouble>(values);
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<int>(doflist);
+	if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
+}/*}}}*/
+
+/*HO*/
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixHO(Element* element){/*{{{*/
+
+	/*compute all stiffness matrices for this element*/
+	ElementMatrix* Ke1=CreateKMatrixHOViscous(element);
+	ElementMatrix* Ke2=CreateKMatrixHOFriction(element);
+	ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
+
+	/*clean-up and return*/
+	delete Ke1;
+	delete Ke2;
+	return Ke;
+}/*}}}*/
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixHOViscous(Element* element){/*{{{*/
+
+	/*Intermediaries*/
+	IssmDouble  viscosity,newviscosity,oldviscosity;
+	IssmDouble  viscosity_overshoot,thickness,Jdet;
+	IssmDouble  D_scalar;
+	IssmDouble *xyz_list = 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(HOApproximationEnum);
+	IssmDouble*    B      = xNew<IssmDouble>(5*numdof);
+	IssmDouble*    Bprime = xNew<IssmDouble>(5*numdof);
+	IssmDouble*    D      = xNewZeroInit<IssmDouble>(5*5);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+	Input* vx_input=element->GetInput(VxEnum);               _assert_(vx_input);
+	Input* vy_input=element->GetInput(VyEnum);               _assert_(vy_input);
+	Input* vxold_input=element->GetInput(VxPicardEnum);      _assert_(vxold_input);
+	Input* vyold_input=element->GetInput(VyPicardEnum);      _assert_(vyold_input);
+	element->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum);
+
+	/* 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->GetBHO(B,element,xyz_list,gauss);
+		this->GetBHOprime(Bprime,element,xyz_list,gauss);
+
+		element->ViscosityHO(&viscosity,xyz_list,gauss,vx_input,vy_input);
+		element->ViscosityHO(&oldviscosity,xyz_list,gauss,vxold_input,vyold_input);
+
+		newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
+		D_scalar=2.*newviscosity*gauss->weight*Jdet;
+		for(int i=0;i<5;i++) D[i*5+i]=D_scalar;
+
+		TripleMultiply(B,5,numdof,1,
+					D,5,5,0,
+					Bprime,5,numdof,0,
+					&Ke->values[0],1);
+	}
+
+	/*Transform Coordinate System*/
+	element->TransformStiffnessMatrixCoord(Ke,XYEnum);
+
+	/*Clean up and return*/
+	delete gauss;
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(D);
+	xDelete<IssmDouble>(Bprime);
+	xDelete<IssmDouble>(B);
+	return Ke;
+}/*}}}*/
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixHOFriction(Element* element){/*{{{*/
+
+	if(element->IsFloating() || !element->IsOnBed()) return NULL;
+
+	return NULL;
+}/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorHO(Element* element){/*{{{*/
+
+	/*compute all load vectors for this element*/
+	ElementVector* pe1=CreatePVectorHODrivingStress(element);
+	ElementVector* pe2=CreatePVectorHOFront(element);
+	ElementVector* pe =new ElementVector(pe1,pe2);
+
+	/*clean-up and return*/
+	delete pe1;
+	delete pe2;
+	return pe;
+}/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorHODrivingStress(Element* element){/*{{{*/
+
+	/*Intermediaries */
+	IssmDouble  Jdet,slope[3];
+	IssmDouble* xyz_list = NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Initialize Element vector and vectors*/
+	ElementVector* pe=element->NewElementVector(HOApproximationEnum);
+	IssmDouble*    basis = xNew<IssmDouble>(numnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+	Input*     surface_input = element->GetInput(SurfaceEnum);   _assert_(surface_input);
+	IssmDouble rhog = element->GetMaterialParameter(MaterialsRhoIceEnum)*element->GetMaterialParameter(ConstantsGEnum);
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=element->NewGauss(3);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		element->NodalFunctions(basis, gauss);
+		surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
+
+		for(int i=0;i<numnodes;i++){
+			pe->values[i*2+0]+=-rhog*slope[0]*Jdet*gauss->weight*basis[i];
+			pe->values[i*2+1]+=-rhog*slope[1]*Jdet*gauss->weight*basis[i];
+		}
+	}
+
+	/*Transform coordinate system*/
+	element->TransformLoadVectorCoord(pe,XYEnum);
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(basis);
+	xDelete<IssmDouble>(xyz_list);
+	delete gauss;
+	return pe;
+}/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorHOFront(Element* element){/*{{{*/
+
+	return NULL;
+}/*}}}*/
+void StressbalanceAnalysis::GetBHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. 
+	 * For node i, Bi can be expressed in the actual coordinate system
+	 * by: 
+	 *       Bi=[ dh/dx          0      ]
+	 *          [   0           dh/dy   ]
+	 *          [ 1/2*dh/dy  1/2*dh/dx  ]
+	 *          [ 1/2*dh/dz      0      ]
+	 *          [  0         1/2*dh/dz  ]
+	 * where h is the interpolation function for node i.
+	 *
+	 * We assume B has been allocated already, of size: 5x(NDOF2*numnodes)
+	 */
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Get nodal functions derivatives*/
+	IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
+	element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+	/*Build B: */
+	for(int i=0;i<numnodes;i++){
+		B[2*numnodes*0+2*i+0] = dbasis[0*numnodes+i];
+		B[2*numnodes*0+2*i+1] = 0.;
+		B[2*numnodes*1+2*i+0] = 0.;
+		B[2*numnodes*1+2*i+1] = dbasis[1*numnodes+i];
+		B[2*numnodes*2+2*i+0] = .5*dbasis[1*numnodes+i];
+		B[2*numnodes*2+2*i+1] = .5*dbasis[0*numnodes+i];
+		B[2*numnodes*3+2*i+0] = .5*dbasis[2*numnodes+i];
+		B[2*numnodes*3+2*i+1] = 0.;
+		B[2*numnodes*4+2*i+0] = 0.;
+		B[2*numnodes*4+2*i+1] = .5*dbasis[2*numnodes+i];
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(dbasis);
+}/*}}}*/
+void StressbalanceAnalysis::GetBHOprime(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
+	 * For node i, Bi' can be expressed in the actual coordinate system
+	 * by: 
+	 *       Bi_prime=[ 2*dN/dx    dN/dy ]
+	 *                [   dN/dx  2*dN/dy ]
+	 *                [   dN/dy    dN/dx ]
+	 * where hNis the finiteelement function for node i.
+	 *
+	 * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
+	 */
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Get nodal functions derivatives*/
+	IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
+	element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
+
+	/*Build B': */
+	for(int i=0;i<numnodes;i++){
+		Bprime[2*numnodes*0+2*i+0] = 2.*dbasis[0*numnodes+i];
+		Bprime[2*numnodes*0+2*i+1] = dbasis[1*numnodes+i];
+		Bprime[2*numnodes*1+2*i+0] = dbasis[0*numnodes+i];
+		Bprime[2*numnodes*1+2*i+1] = 2.*dbasis[1*numnodes+i];
+		Bprime[2*numnodes*2+2*i+0] = dbasis[1*numnodes+i];
+		Bprime[2*numnodes*2+2*i+1] = dbasis[0*numnodes+i];
+		Bprime[2*numnodes*3+2*i+0] = dbasis[2*numnodes+i];
+		Bprime[2*numnodes*3+2*i+1] = 0.;
+		Bprime[2*numnodes*4+2*i+0] = 0.;
+		Bprime[2*numnodes*4+2*i+1] = dbasis[2*numnodes+i];
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(dbasis);
 }/*}}}*/
 void StressbalanceAnalysis::InputUpdateFromSolutionHO(IssmDouble* solution,Element* element){/*{{{*/
@@ -1619,4 +1656,328 @@
 	xDelete<int>(doflist);
 }/*}}}*/
+
+/*FS*/
+ElementVector* StressbalanceAnalysis::CreatePVectorFS(Element* element){/*{{{*/
+
+	/*compute all load vectors for this element*/
+	ElementVector* pe1=CreatePVectorFSViscous(element);
+	ElementVector* pe2=CreatePVectorFSShelf(element);
+	ElementVector* pe3=CreatePVectorFSFront(element);
+	ElementVector* pe =new ElementVector(pe1,pe2,pe3);
+
+	/*clean-up and return*/
+	delete pe1;
+	delete pe2;
+	delete pe3;
+	return pe;
+}/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorFSViscous(Element* element){/*{{{*/
+
+	int         i,meshtype,dim;
+	IssmDouble  Jdet,forcex,forcey,forcez;
+	IssmDouble *xyz_list = NULL;
+
+	/*Get problem dimension*/
+	element->FindParam(&meshtype,MeshTypeEnum);
+	switch(meshtype){
+		case Mesh2DverticalEnum: dim = 2; break;
+		case Mesh3DEnum:         dim = 3; 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();
+
+	/*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 vectors*/
+	ElementVector* pe     = element->NewElementVector(FSvelocityEnum);
+	IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+	Input*      loadingforcex_input=element->GetInput(LoadingforceXEnum);  _assert_(loadingforcex_input);
+	Input*      loadingforcey_input=element->GetInput(LoadingforceYEnum);  _assert_(loadingforcey_input);
+	Input*      loadingforcez_input=element->GetInput(LoadingforceZEnum);  _assert_(loadingforcez_input);
+	IssmDouble  rho_ice =element->GetMaterialParameter(MaterialsRhoIceEnum);
+	IssmDouble  gravity =element->GetMaterialParameter(ConstantsGEnum);
+
+	/* 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->NodalFunctionsVelocity(vbasis,gauss);
+
+		loadingforcex_input->GetInputValue(&forcex,gauss);
+		loadingforcey_input->GetInputValue(&forcey,gauss);
+		if(dim==3) loadingforcez_input->GetInputValue(&forcez,gauss);
+
+		for(i=0;i<vnumnodes;i++){
+			pe->values[i*dim+0] += +rho_ice*forcex *Jdet*gauss->weight*vbasis[i];
+			pe->values[i*dim+1] += +rho_ice*forcey *Jdet*gauss->weight*vbasis[i];
+			if(dim==3){
+				pe->values[i*dim+2] += +rho_ice*forcez*Jdet*gauss->weight*vbasis[i];
+				pe->values[i*dim+2] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i];
+			}
+			else{
+				pe->values[i*dim+1] += -rho_ice*gravity*Jdet*gauss->weight*vbasis[i];
+			}
+		}
+	}
+
+	/*Transform coordinate system*/
+	element->TransformLoadVectorCoord(pe,cs_list);
+
+	/*Clean up and return*/
+	delete gauss;
+	xDelete<int>(cs_list);
+	xDelete<IssmDouble>(vbasis);
+	xDelete<IssmDouble>(xyz_list);
+	return pe;
+}/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorFSShelf(Element* element){/*{{{*/
+
+	int         i,meshtype,dim;
+	IssmDouble  Jdet,water_pressure,bed;
+	IssmDouble	normal[3];
+	IssmDouble *xyz_list_base = NULL;
+
+	/*Get basal element*/
+	if(!element->IsOnBed() || !element->IsFloating()) return NULL;
+
+	/*Get problem dimension*/
+	element->FindParam(&meshtype,MeshTypeEnum);
+	switch(meshtype){
+		case Mesh2DverticalEnum: dim = 2; break;
+		case Mesh3DEnum:         dim = 3; 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();
+
+	/*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 vectors*/
+	ElementVector* pe     = element->NewElementVector(FSvelocityEnum);
+	IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinatesBase(&xyz_list_base);
+	Input*      bed_input=element->GetInput(BedEnum); _assert_(bed_input);
+	IssmDouble  rho_water=element->GetMaterialParameter(MaterialsRhoWaterEnum);
+	IssmDouble  gravity  =element->GetMaterialParameter(ConstantsGEnum);
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=element->NewGaussBase(5);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
+		element->NodalFunctionsVelocity(vbasis,gauss);
+
+		element->NormalBase(&normal[0],xyz_list_base);
+		_assert_(normal[dim-1]<0.);
+		bed_input->GetInputValue(&bed, gauss);
+		water_pressure=gravity*rho_water*bed;
+
+		for(i=0;i<vnumnodes;i++){
+			pe->values[i*dim+0] += water_pressure*gauss->weight*Jdet*vbasis[i]*normal[0];
+			pe->values[i*dim+1] += water_pressure*gauss->weight*Jdet*vbasis[i]*normal[1];
+			if(dim==3){
+				pe->values[i*dim+2]+=water_pressure*gauss->weight*Jdet*vbasis[i]*normal[2];
+			}
+		}
+	}
+
+	/*Transform coordinate system*/
+	element->TransformLoadVectorCoord(pe,cs_list);
+
+	/*Clean up and return*/
+	delete gauss;
+	xDelete<int>(cs_list);
+	xDelete<IssmDouble>(vbasis);
+	xDelete<IssmDouble>(xyz_list_base);
+	return pe;
+}/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorFSFront(Element* element){/*{{{*/
+
+	return NULL;
+}/*}}}*/
+void StressbalanceAnalysis::GetSolutionFromInputsFS(Vector<IssmDouble>* solution,Element* element){/*{{{*/
+
+	int*         vdoflist=NULL;
+	int*         pdoflist=NULL;
+	Input*       vz_input=NULL;
+	int          meshtype,dim;
+	IssmDouble   vx,vy,vz,p;
+	IssmDouble   FSreconditioning;
+
+	/*Get some parameters*/
+	element->FindParam(&meshtype,MeshTypeEnum);
+	element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
+	switch(meshtype){
+		case Mesh2DverticalEnum: dim = 2; break;
+		case Mesh3DEnum:         dim = 3; break;
+		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
+	}
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int vnumnodes = element->NumberofNodesVelocity();
+	int pnumnodes = element->NumberofNodesPressure();
+	int vnumdof   = vnumnodes*dim;
+	int pnumdof   = pnumnodes*1;
+
+	/*Initialize values*/
+	IssmDouble* vvalues = xNew<IssmDouble>(vnumdof);
+	IssmDouble* pvalues = xNew<IssmDouble>(pnumdof);
+
+	/*Get dof list: */
+	element->GetDofListVelocity(&vdoflist,GsetEnum);
+	element->GetDofListPressure(&pdoflist,GsetEnum);
+	Input*     vx_input=element->GetInput(VxEnum);       _assert_(vx_input);
+	Input*     vy_input=element->GetInput(VyEnum);       _assert_(vy_input);
+	if(dim==3){vz_input=element->GetInput(VzEnum);       _assert_(vz_input);}
+	Input*     p_input =element->GetInput(PressureEnum); _assert_(p_input);
+
+	element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
+
+	/*Ok, we have the velocities in inputs, fill in solution */
+	Gauss* gauss = element->NewGauss();
+	for(int i=0;i<vnumnodes;i++){
+		gauss->GaussNode(element->VelocityInterpolation(),i);
+		vx_input->GetInputValue(&vx,gauss);
+		vy_input->GetInputValue(&vy,gauss);
+		vvalues[i*dim+0]=vx;
+		vvalues[i*dim+1]=vy;
+		if(dim==3){
+			vz_input->GetInputValue(&vz,gauss);
+			vvalues[i*dim+2]=vz;
+		}
+	}
+	for(int i=0;i<pnumnodes;i++){
+		gauss->GaussNode(element->PressureInterpolation(),i);
+		p_input->GetInputValue(&p ,gauss);
+		pvalues[i]=p/FSreconditioning;
+	}
+
+	/*Add value to global vector*/
+	solution->SetValues(vnumdof,vdoflist,vvalues,INS_VAL);
+	solution->SetValues(pnumdof,pdoflist,pvalues,INS_VAL);
+
+	/*Free ressources:*/
+	delete gauss;
+	xDelete<int>(pdoflist);
+	xDelete<int>(vdoflist);
+	xDelete<IssmDouble>(pvalues);
+	xDelete<IssmDouble>(vvalues);
+}/*}}}*/
+void StressbalanceAnalysis::InputUpdateFromSolutionFS(IssmDouble* solution,Element* element){/*{{{*/
+
+	int          i,dim,meshtype;
+	int*         vdoflist=NULL;
+	int*         pdoflist=NULL;
+	IssmDouble   FSreconditioning;
+
+	element->FindParam(&meshtype,MeshTypeEnum);
+	element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
+	switch(meshtype){
+		case Mesh2DverticalEnum: dim = 2; break;
+		case Mesh3DEnum:         dim = 3; 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 vnumdof   = vnumnodes*dim;
+	int pnumdof   = pnumnodes*1;
+
+	/*Initialize values*/
+	IssmDouble* values   = xNew<IssmDouble>(vnumdof+pnumdof);
+	IssmDouble* vx       = xNew<IssmDouble>(vnumnodes);
+	IssmDouble* vy       = xNew<IssmDouble>(vnumnodes);
+	IssmDouble* vz       = xNew<IssmDouble>(vnumnodes);
+	IssmDouble* vel      = xNew<IssmDouble>(vnumnodes);
+	IssmDouble* pressure = xNew<IssmDouble>(pnumnodes);
+
+	/*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;
+
+	/*Get dof list: */
+	element->GetDofListVelocity(&vdoflist,GsetEnum);
+	element->GetDofListPressure(&pdoflist,GsetEnum);
+
+	/*Use the dof list to index into the solution vector: */
+	for(i=0;i<vnumdof;i++) values[i]        =solution[vdoflist[i]];
+	for(i=0;i<pnumdof;i++) values[vnumdof+i]=solution[pdoflist[i]];
+
+	/*Transform solution in Cartesian Space*/
+	element->TransformSolutionCoord(values,cs_list);
+
+	/*Ok, we have vx and vy in values, fill in all arrays: */
+	for(i=0;i<vnumnodes;i++){
+		vx[i] = values[i*dim+0];
+		vy[i] = values[i*dim+1];
+		if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
+
+		if(dim==3){
+			vz[i] = values[i*dim+2];
+			if(xIsNan<IssmDouble>(vz[i])) _error_("NaN found in solution vector");
+		}
+	}
+	for(i=0;i<pnumnodes;i++){
+		pressure[i] = values[vnumdof+i];
+		if(xIsNan<IssmDouble>(pressure[i])) _error_("NaN found in solution vector");
+	}
+
+	/*Recondition pressure and compute vel: */
+	for(i=0;i<pnumnodes;i++) pressure[i] = pressure[i]*FSreconditioning;
+	if(dim==3) for(i=0;i<vnumnodes;i++) vel[i] = sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
+	else       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: */
+	element->InputChangeName(VxEnum,VxPicardEnum);
+	element->InputChangeName(VyEnum,VyPicardEnum);
+	element->InputChangeName(PressureEnum,PressurePicardEnum);
+	if(dim==3) element->InputChangeName(VzEnum,VzPicardEnum);
+
+	/*Add vx and vy as inputs to the tria element: */
+	element->AddInput(VxEnum,vx,P1Enum);
+	element->AddInput(VyEnum,vy,P1Enum);
+	element->AddInput(VelEnum,vel,P1Enum);
+	element->AddInput(PressureEnum,pressure,P1Enum);
+	if(dim==3) element->AddInput(VzEnum,vz,P1Enum);
+
+	/*Free ressources:*/
+	xDelete<IssmDouble>(pressure);
+	xDelete<IssmDouble>(vel);
+	xDelete<IssmDouble>(vz);
+	xDelete<IssmDouble>(vy);
+	xDelete<IssmDouble>(vx);
+	xDelete<IssmDouble>(values);
+	xDelete<int>(vdoflist);
+	xDelete<int>(pdoflist);
+	xDelete<int>(cs_list);
+}/*}}}*/
+
+/*Coupling (Tiling)*/
 void StressbalanceAnalysis::InputUpdateFromSolutionHOFS(IssmDouble* solution,Element* element){/*{{{*/
 
@@ -1716,202 +2077,4 @@
 	xDelete<int>(cs_list);
 }/*}}}*/
-void StressbalanceAnalysis::InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element){/*{{{*/
-
-	int         i,meshtype;
-	IssmDouble  rho_ice,g;
-	int*        doflist=NULL;
-	IssmDouble* xyz_list=NULL;
-	Element*    basalelement=NULL;
-
-	/*Deal with pressure first*/
-	int numvertices = element->GetNumberOfVertices();
-	IssmDouble* pressure  = xNew<IssmDouble>(numvertices);
-	IssmDouble* thickness = xNew<IssmDouble>(numvertices);
-	IssmDouble* surface   = xNew<IssmDouble>(numvertices);
-
-	element->FindParam(&meshtype,MeshTypeEnum);
-	rho_ice =element->GetMaterialParameter(MaterialsRhoIceEnum);
-	g       =element->GetMaterialParameter(ConstantsGEnum);
-	switch(meshtype){
-		case Mesh2DhorizontalEnum:
-			element->GetInputListOnVertices(thickness,ThicknessEnum);
-			for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*thickness[i];
-			break;
-		case Mesh3DEnum:
-			element->GetVerticesCoordinates(&xyz_list);
-			element->GetInputListOnVertices(surface,SurfaceEnum);
-			for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]);
-			break;
-		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
-	}
-	element->AddInput(PressureEnum,pressure,P1Enum);
-	xDelete<IssmDouble>(pressure);
-	xDelete<IssmDouble>(thickness);
-	xDelete<IssmDouble>(surface);
-
-	/*Get basal element*/
-	switch(meshtype){
-		case Mesh2DhorizontalEnum:
-			basalelement = element;
-			break;
-		case Mesh3DEnum:
-			if(!element->IsOnBed()){xDelete<IssmDouble>(xyz_list); return;}
-			basalelement=element->SpawnBasalElement();
-			break;
-		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
-	}
-
-	/*Fetch number of nodes and dof for this finite element*/
-	int numnodes = basalelement->GetNumberOfNodes();
-	int numdof   = numnodes*2;
-
-	/*Fetch dof list and allocate solution vectors*/
-	basalelement->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-	IssmDouble* values    = xNew<IssmDouble>(numdof);
-	IssmDouble* vx        = xNew<IssmDouble>(numnodes);
-	IssmDouble* vy        = xNew<IssmDouble>(numnodes);
-	IssmDouble* vz        = xNew<IssmDouble>(numnodes);
-	IssmDouble* vel       = xNew<IssmDouble>(numnodes);
-
-	/*Use the dof list to index into the solution vector: */
-	for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
-
-	/*Transform solution in Cartesian Space*/
-	basalelement->TransformSolutionCoord(&values[0],XYEnum);
-	basalelement->FindParam(&meshtype,MeshTypeEnum);
-
-	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
-	for(i=0;i<numnodes;i++){
-		vx[i]=values[i*2+0];
-		vy[i]=values[i*2+1];
-
-		/*Check solution*/
-		if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
-		if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
-	}
-
-	/*Get Vz and compute vel*/
-	basalelement->GetInputListOnNodes(&vz[0],VzEnum,0.);
-	for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
-
-	/*Now, we have to move the previous Vx and Vy inputs  to old 
-	 * status, otherwise, we'll wipe them off: */
-	element->InputChangeName(VxEnum,VxPicardEnum);
-	element->InputChangeName(VyEnum,VyPicardEnum);
-
-	/*Add vx and vy as inputs to the tria element: */
-	element->AddBasalInput(VxEnum,vx,P1Enum);
-	element->AddBasalInput(VyEnum,vy,P1Enum);
-	element->AddBasalInput(VelEnum,vel,P1Enum);
-
-	/*Free ressources:*/
-	xDelete<IssmDouble>(vel);
-	xDelete<IssmDouble>(vz);
-	xDelete<IssmDouble>(vy);
-	xDelete<IssmDouble>(vx);
-	xDelete<IssmDouble>(values);
-	xDelete<IssmDouble>(xyz_list);
-	xDelete<int>(doflist);
-	if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
-}/*}}}*/
-void StressbalanceAnalysis::InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element){/*{{{*/
-
-	int         i,meshtype;
-	IssmDouble  rho_ice,g;
-	int*        doflist=NULL;
-	IssmDouble* xyz_list=NULL;
-	Element*    basalelement=NULL;
-
-	/*Deal with pressure first*/
-	int numvertices = element->GetNumberOfVertices();
-	IssmDouble* pressure  = xNew<IssmDouble>(numvertices);
-	IssmDouble* thickness = xNew<IssmDouble>(numvertices);
-	IssmDouble* surface   = xNew<IssmDouble>(numvertices);
-
-	element->FindParam(&meshtype,MeshTypeEnum);
-	rho_ice =element->GetMaterialParameter(MaterialsRhoIceEnum);
-	g       =element->GetMaterialParameter(ConstantsGEnum);
-	switch(meshtype){
-		case Mesh2DhorizontalEnum:
-			element->GetInputListOnVertices(thickness,ThicknessEnum);
-			for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*thickness[i];
-			break;
-		case Mesh3DEnum:
-			element->GetVerticesCoordinates(&xyz_list);
-			element->GetInputListOnVertices(surface,SurfaceEnum);
-			for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]);
-			break;
-		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
-	}
-	element->AddInput(PressureEnum,pressure,P1Enum);
-	xDelete<IssmDouble>(pressure);
-	xDelete<IssmDouble>(thickness);
-	xDelete<IssmDouble>(surface);
-
-	/*Get basal element*/
-	switch(meshtype){
-		case Mesh2DhorizontalEnum:
-			basalelement = element;
-			break;
-		case Mesh3DEnum:
-			if(!element->IsOnBed()){xDelete<IssmDouble>(xyz_list); return;}
-			basalelement=element->SpawnBasalElement();
-			break;
-		default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
-	}
-
-	/*Fetch number of nodes and dof for this finite element*/
-	int numnodes = basalelement->GetNumberOfNodes();
-	int numdof   = numnodes*2;
-
-	/*Fetch dof list and allocate solution vectors*/
-	basalelement->GetDofList(&doflist,SSAApproximationEnum,GsetEnum);
-	IssmDouble* values    = xNew<IssmDouble>(numdof);
-	IssmDouble* vx        = xNew<IssmDouble>(numnodes);
-	IssmDouble* vy        = xNew<IssmDouble>(numnodes);
-	IssmDouble* vz        = xNew<IssmDouble>(numnodes);
-	IssmDouble* vel       = xNew<IssmDouble>(numnodes);
-
-	/*Use the dof list to index into the solution vector: */
-	for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
-
-	/*Transform solution in Cartesian Space*/
-	basalelement->TransformSolutionCoord(&values[0],XYEnum);
-	basalelement->FindParam(&meshtype,MeshTypeEnum);
-
-	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
-	for(i=0;i<numnodes;i++){
-		vx[i]=values[i*2+0];
-		vy[i]=values[i*2+1];
-
-		/*Check solution*/
-		if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
-		if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
-	}
-
-	/*Get Vz and compute vel*/
-	basalelement->GetInputListOnNodes(&vz[0],VzEnum,0.);
-	for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
-
-	/*Now, we have to move the previous Vx and Vy inputs  to old 
-	 * status, otherwise, we'll wipe them off: */
-	element->InputChangeName(VxEnum,VxPicardEnum);
-	element->InputChangeName(VyEnum,VyPicardEnum);
-
-	/*Add vx and vy as inputs to the tria element: */
-	element->AddBasalInput(VxEnum,vx,P1Enum);
-	element->AddBasalInput(VyEnum,vy,P1Enum);
-	element->AddBasalInput(VelEnum,vel,P1Enum);
-
-	/*Free ressources:*/
-	xDelete<IssmDouble>(vel);
-	xDelete<IssmDouble>(vz);
-	xDelete<IssmDouble>(vy);
-	xDelete<IssmDouble>(vx);
-	xDelete<IssmDouble>(values);
-	xDelete<IssmDouble>(xyz_list);
-	xDelete<int>(doflist);
-	if(meshtype!=Mesh2DhorizontalEnum){basalelement->DeleteMaterials(); delete basalelement;};
-}/*}}}*/
 void StressbalanceAnalysis::InputUpdateFromSolutionSSAFS(IssmDouble* solution,Element* element){/*{{{*/
 
Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 16828)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 16829)
@@ -22,4 +22,7 @@
 		/*Finite element Analysis*/
 		ElementMatrix* CreateKMatrix(Element* element);
+		ElementMatrix* CreateKMatrixHO(Element* element);
+		ElementMatrix* CreateKMatrixHOViscous(Element* element);
+		ElementMatrix* CreateKMatrixHOFriction(Element* element);
 		ElementMatrix* CreateKMatrixSSA(Element* element);
 		ElementMatrix* CreateKMatrixSSAViscous(Element* element);
@@ -36,4 +39,6 @@
 		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);
Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 16828)
+++ /issm/trunk-jpl/src/c/classes/Elements/Element.h	(revision 16829)
@@ -144,4 +144,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   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;
 		virtual int    VelocityInterpolation()=0;
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16828)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16829)
@@ -1543,5 +1543,5 @@
 /*}}}*/
 /*FUNCTION Penta::GetStrainRate3dHO{{{*/
-void Penta::GetStrainRate3dHO(IssmDouble* epsilon,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input){
+void Penta::GetStrainRate3dHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss, Input* vx_input, Input* vy_input){
 	/*Compute the 3d Blatter/HOStrain Rate (5 components):
 	 *
@@ -1572,5 +1572,5 @@
 /*}}}*/
 /*FUNCTION Penta::GetStrainRate3d{{{*/
-void Penta::GetStrainRate3d(IssmDouble* epsilon,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input, Input* vz_input){
+void Penta::GetStrainRate3d(IssmDouble* epsilon,IssmDouble* xyz_list, Gauss* gauss, Input* vx_input, Input* vy_input, Input* vz_input){
 	/*Compute the 3d Strain Rate (6 components):
 	 *
@@ -3482,4 +3482,18 @@
 	/*Assign output pointer*/
 	*pphi = phi;
+}
+/*}}}*/
+/*FUNCTION Penta::ViscosityHO{{{*/
+void Penta::ViscosityHO(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){
+
+	/*Intermediaries*/
+	IssmDouble viscosity;
+	IssmDouble epsilon[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/
+
+	this->GetStrainRate3dHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
+	material->GetViscosity3d(&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 16828)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 16829)
@@ -230,6 +230,6 @@
 		void	         GetPhi(IssmDouble* phi, IssmDouble*  epsilon, IssmDouble viscosity);
 		void           GetQuadNormal(IssmDouble* normal,IssmDouble xyz_list[4][3]);
-		void           GetStrainRate3dHO(IssmDouble* epsilon,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input);
-		void           GetStrainRate3d(IssmDouble* epsilon,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input, Input* vz_input);
+		void           GetStrainRate3dHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss, Input* vx_input, Input* vy_input);
+		void           GetStrainRate3d(IssmDouble* epsilon,IssmDouble* xyz_list, Gauss* gauss, Input* vx_input, Input* vy_input, Input* vz_input);
 		Penta*         GetUpperElement(void);
 		void           GetZeroLevelsetCoordinates(IssmDouble* xyz_zero,IssmDouble xyz_list[6][3],int levelsetenum);
@@ -275,4 +275,5 @@
 		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           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 16828)
+++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h	(revision 16829)
@@ -145,4 +145,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           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");};
 		#ifdef _HAVE_THERMAL_
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16828)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 16829)
@@ -305,4 +305,5 @@
 		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           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);