Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h
===================================================================
--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 18927)
+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.h	(revision 18928)
@@ -12,12 +12,12 @@
 
   public:
 		/*Model processing*/
+		void CreateConstraints(Constraints* constraints,IoModel* iomodel);
+		void CreateLoads(Loads* loads, IoModel* iomodel);
+		void CreateNodes(Nodes* nodes,IoModel* iomodel);
 		int  DofsPerNode(int** doflist,int domaintype,int approximation);
+		void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
 		void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum);
-		void UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type);
-		void CreateNodes(Nodes* nodes,IoModel* iomodel);
-		void CreateConstraints(Constraints* constraints,IoModel* iomodel);
-		void CreateLoads(Loads* loads, IoModel* iomodel);
 
 		/*Finite element Analysis*/
 		void           Core(FemModel* femmodel);
@@ -25,92 +25,92 @@
 		ElementMatrix* CreateJacobianMatrix(Element* element);
 		ElementMatrix* CreateKMatrix(Element* element);
 		ElementVector* CreatePVector(Element* element);
-		void GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
-		void GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
-		void GetSolutionFromInputsHoriz(Vector<IssmDouble>* solution,Element* element);
-		void InputUpdateFromSolution(IssmDouble* solution,Element* element);
-		void UpdateConstraints(FemModel* femmodel);
+		void           GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element);
+		void           GetSolutionFromInputsHoriz(Vector<IssmDouble>* solution,Element* element);
+		void           GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index);
+		void           InputUpdateFromSolution(IssmDouble* solution,Element* element);
+		void           UpdateConstraints(FemModel* femmodel);
 
 		/*SSA*/
 		ElementMatrix* CreateJacobianMatrixSSA(Element* element);
 		ElementMatrix* CreateKMatrixSSA(Element* element);
-		ElementMatrix* CreateKMatrixSSAViscous(Element* element);
 		ElementMatrix* CreateKMatrixSSAFriction(Element* element);
 		ElementMatrix* CreateKMatrixSSALateralFriction(Element* element);
+		ElementMatrix* CreateKMatrixSSAViscous(Element* element);
 		ElementVector* CreatePVectorSSA(Element* element);
+		ElementVector* CreatePVectorSSAFront(Element* element);
 		ElementVector* CreatePVectorSSADrivingStress(Element* element);
-		ElementVector* CreatePVectorSSAFront(Element* element);
-		void GetBSSA(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
-		void GetBSSAprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
-		void GetBSSAFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
-		void InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element);
+		void           GetBSSA(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
+		void           GetBSSAFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
+		void           GetBSSAprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
+		void           InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element);
 		/*L1L2*/
 		ElementMatrix* CreateKMatrixL1L2(Element* element);
+		ElementMatrix* CreateKMatrixL1L2Friction(Element* element);
 		ElementMatrix* CreateKMatrixL1L2Viscous(Element* element);
-		ElementMatrix* CreateKMatrixL1L2Friction(Element* element);
 		ElementVector* CreatePVectorL1L2(Element* element);
+		ElementVector* CreatePVectorL1L2Front(Element* element);
 		ElementVector* CreatePVectorL1L2DrivingStress(Element* element);
-		ElementVector* CreatePVectorL1L2Front(Element* element);
-		void InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element);
+		void           InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element);
 		/*HO*/
 		ElementMatrix* CreateJacobianMatrixHO(Element* element);
 		ElementMatrix* CreateKMatrixHO(Element* element);
+		ElementMatrix* CreateKMatrixHOFriction(Element* element);
 		ElementMatrix* CreateKMatrixHOViscous(Element* element);
-		ElementMatrix* CreateKMatrixHOFriction(Element* element);
 		ElementVector* CreatePVectorHO(Element* element);
+		ElementVector* CreatePVectorHOFront(Element* element);
 		ElementVector* CreatePVectorHODrivingStress(Element* element);
-		ElementVector* CreatePVectorHOFront(Element* element);
-		void GetBHO(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
-		void GetBHOprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
-		void GetBHOFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
-		void InputUpdateFromSolutionHO(IssmDouble* solution,Element* element);
+		void           GetBHO(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
+		void           GetBHOFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
+		void           GetBHOprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
+		void           InputUpdateFromSolutionHO(IssmDouble* solution,Element* element);
 		/*FS*/
 		ElementVector* CreateDVectorFS(Element* element);
 		ElementMatrix* CreateJacobianMatrixFS(Element* element);
 		ElementMatrix* CreateKMatrixFS(Element* element);
+		ElementMatrix* CreateKMatrixFSFriction(Element* element);
+		ElementMatrix* CreateKMatrixFSShelf(Element* element);
+		ElementMatrix* CreateKMatrixFSViscous(Element* element);
 		ElementMatrix* CreateKMatrixFSViscousLA(Element* element);
 		ElementMatrix* CreateKMatrixFSViscousXTH(Element* element);
-		ElementMatrix* CreateKMatrixFSViscous(Element* element);
-		ElementMatrix* CreateKMatrixFSFriction(Element* element);
-		ElementMatrix* CreateKMatrixFSShelf(Element* element);
 		ElementVector* CreatePVectorFS(Element* element);
+		ElementVector* CreatePVectorFSFriction(Element* element);
+		ElementVector* CreatePVectorFSFront(Element* element);
+		ElementVector* CreatePVectorFSShelf(Element* element);
+		ElementVector* CreatePVectorFSStress(Element* element);
 		ElementVector* CreatePVectorFSViscous(Element* element);
 		ElementVector* CreatePVectorFSViscousLA(Element* element);
 		ElementVector* CreatePVectorFSViscousXTH(Element* element);
-		ElementVector* CreatePVectorFSShelf(Element* element);
-		ElementVector* CreatePVectorFSFront(Element* element);
-		ElementVector* CreatePVectorFSFriction(Element* element);
-		ElementVector* CreatePVectorFSStress(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 GetBFSvel(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
-		void GetBFSprimevel(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
-		void GetBFSFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
-		void GetBFSUzawa(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
-		void GetBFSprimeUzawa(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
-		void GetCFS(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
-		void GetCFSprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
-		void GetSolutionFromInputsFS(Vector<IssmDouble>* solution,Element* element);
-		void InputUpdateFromSolutionFS(IssmDouble* solution,Element* element);
-		void InputUpdateFromSolutionFSXTH_d(Elements* elements,Parameters* parameters);
-		void InputUpdateFromSolutionFSXTH_tau(Elements* elements,Parameters* parameters);
-		void InitializeXTH(Elements* elements,Parameters* parameters);
+		void           GetBFS(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
+		void           GetBFSFriction(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           GetBFSprimeUzawa(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
+		void           GetBFSprimevel(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
+		void           GetBFSUzawa(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
+		void           GetBFSvel(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
+		void           GetCFS(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
+		void           GetCFSprime(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss);
+		void           GetSolutionFromInputsFS(Vector<IssmDouble>* solution,Element* element);
+		void           InitializeXTH(Elements* elements,Parameters* parameters);
+		void           InputUpdateFromSolutionFS(IssmDouble* solution,Element* element);
+		void           InputUpdateFromSolutionFSXTH_d(Elements* elements,Parameters* parameters);
+		void           InputUpdateFromSolutionFSXTH_tau(Elements* elements,Parameters* parameters);
 		/*Coupling*/
-		ElementMatrix* CreateKMatrixSSA3d(Element* element);
-		ElementMatrix* CreateKMatrixSSA3dFriction(Element* element);
-		ElementMatrix* CreateKMatrixSSA3dViscous(Element* element);
-		ElementMatrix* CreateKMatrixHOFS(Element* element);
-		ElementMatrix* CreateKMatrixSSAHO(Element* element);
-		ElementMatrix* CreateKMatrixSSAFS(Element* element);
 		ElementMatrix* CreateKMatrixCouplingHOFS(Element* element);
+		ElementMatrix* CreateKMatrixCouplingSSAFS(Element* element);
+		ElementMatrix* CreateKMatrixCouplingSSAFSFriction(Element* element);
+		ElementMatrix* CreateKMatrixCouplingSSAFSViscous(Element* element);
 		ElementMatrix* CreateKMatrixCouplingSSAHO(Element* element);
 		ElementMatrix* CreateKMatrixCouplingSSAHOFriction(Element* element);
 		ElementMatrix* CreateKMatrixCouplingSSAHOViscous(Element* element);
-		ElementMatrix* CreateKMatrixCouplingSSAFS(Element* element);
-		ElementMatrix* CreateKMatrixCouplingSSAFSFriction(Element* element);
-		ElementMatrix* CreateKMatrixCouplingSSAFSViscous(Element* element);
+		ElementMatrix* CreateKMatrixHOFS(Element* element);
+		ElementMatrix* CreateKMatrixSSAFS(Element* element);
+		ElementMatrix* CreateKMatrixSSAHO(Element* element);
+		ElementMatrix* CreateKMatrixSSA3d(Element* element);
+		ElementMatrix* CreateKMatrixSSA3dFriction(Element* element);
+		ElementMatrix* CreateKMatrixSSA3dViscous(Element* element);
+		ElementVector* CreatePVectorSSAFS(Element* element);
 		ElementVector* CreatePVectorSSAHO(Element* element);
-		ElementVector* CreatePVectorSSAFS(Element* element);
 		ElementVector* CreatePVectorCouplingSSAFS(Element* element);
 		ElementVector* CreatePVectorCouplingSSAFSFriction(Element* element);
 		ElementVector* CreatePVectorCouplingSSAFSViscous(Element* element);
@@ -118,17 +118,17 @@
 		ElementVector* CreatePVectorCouplingHOFS(Element* element);
 		ElementVector* CreatePVectorCouplingHOFSFriction(Element* element);
 		ElementVector* CreatePVectorCouplingHOFSViscous(Element* element);
-		void GetBSSAHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
-		void GetBSSAFS(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
-		void GetBprimeSSAFS(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss);
-		void GetBSSAFSTria(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
-		void GetBprimeSSAFSTria(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss);
-		void GetLFSSSA(IssmDouble* L,Element* element,Gauss* gauss);
-		void GetLSSAFS(IssmDouble* L,Element* element,Gauss* gauss);
-		void GetLprimeFSSSA(IssmDouble* Lprime,Element* element,IssmDouble* xyz_list,Gauss* gauss);
-		void GetLprimeSSAFS(IssmDouble* Lprime,Element* element,IssmDouble* xyz_list,Gauss* gauss);
-		void InputUpdateFromSolutionHOFS(IssmDouble* solution,Element* element);
-		void InputUpdateFromSolutionSSAFS(IssmDouble* solution,Element* element);
-		void InputUpdateFromSolutionSSAHO(IssmDouble* solution,Element* element);
+		void           GetBprimeSSAFS(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss);
+		void           GetBprimeSSAFSTria(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss);
+		void           GetBSSAFS(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
+		void           GetBSSAFSTria(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
+		void           GetBSSAHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss);
+		void           GetLFSSSA(IssmDouble* L,Element* element,Gauss* gauss);
+		void           GetLprimeFSSSA(IssmDouble* Lprime,Element* element,IssmDouble* xyz_list,Gauss* gauss);
+		void           GetLprimeSSAFS(IssmDouble* Lprime,Element* element,IssmDouble* xyz_list,Gauss* gauss);
+		void           GetLSSAFS(IssmDouble* L,Element* element,Gauss* gauss);
+		void           InputUpdateFromSolutionHOFS(IssmDouble* solution,Element* element);
+		void           InputUpdateFromSolutionSSAFS(IssmDouble* solution,Element* element);
+		void           InputUpdateFromSolutionSSAHO(IssmDouble* solution,Element* element);
 };
 #endif
Index: ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 18927)
+++ ../trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 18928)
@@ -10,394 +10,6 @@
 //#define LATERALFRICTION 1
 
 /*Model processing*/
-int  StressbalanceAnalysis::DofsPerNode(int** pdoftype,int domaintype,int approximation){/*{{{*/
-
-	/*output*/
-	int *doftype = NULL;
-	int  numdofs;
-
-	switch(approximation){
-		case SSAApproximationEnum:
-			 switch(domaintype){
-				 case Domain3DEnum:           numdofs=2; break;
-				 case Domain2DhorizontalEnum: numdofs=2; break;
-				 case Domain2DverticalEnum:   numdofs=1; break;
-				 default: _error_("mesh type not supported yet");
-			 }
-			 break;
-		case L1L2ApproximationEnum: numdofs =2; break;
-		case HOApproximationEnum:   
-			 switch(domaintype){
-				 case Domain3DEnum:         numdofs=2; break;
-				 case Domain2DverticalEnum: numdofs=1; break;
-				 default: _error_("mesh type not supported yet");
-			 }
-			 break;
-		case SIAApproximationEnum:  numdofs =2; break;
-		case FSvelocityEnum:
-			 switch(domaintype){
-				 case Domain3DEnum:         numdofs=3; break;
-				 case Domain2DverticalEnum: numdofs=2; break;
-				 default: _error_("mesh type not supported yet");
-			}
-			break;
-		case FSpressureEnum: numdofs=1; break;
-		case NoneApproximationEnum:
-			 switch(domaintype){
-				 case Domain3DEnum:         numdofs=4; break;
-				 case Domain2DverticalEnum: numdofs=3; break;
-				 default: _error_("mesh type not supported yet");
-			}
-			break;
-		case SSAHOApproximationEnum:
-			numdofs=4;
-			doftype=xNew<int>(numdofs);
-			doftype[0]=SSAApproximationEnum;
-			doftype[1]=SSAApproximationEnum;
-			doftype[2]=HOApproximationEnum;
-			doftype[3]=HOApproximationEnum;
-			break;
-		case HOFSApproximationEnum:
-			numdofs=5;
-			doftype=xNew<int>(numdofs);
-			doftype[0]=HOApproximationEnum;
-			doftype[1]=HOApproximationEnum;
-			doftype[2]=FSvelocityEnum;
-			doftype[3]=FSvelocityEnum;
-			doftype[4]=FSvelocityEnum;
-			break;
-		case SSAFSApproximationEnum:
-			numdofs=5;
-			doftype=xNew<int>(numdofs);
-			doftype[0]=SSAApproximationEnum;
-			doftype[1]=SSAApproximationEnum;
-			doftype[2]=FSvelocityEnum;
-			doftype[3]=FSvelocityEnum;
-			doftype[4]=FSvelocityEnum;
-			break;
-		default:
-			_error_("Approximation " << EnumToStringx(approximation) << " not implemented yet");
-	}
-
-	/*Assign output pointer and return*/
-	*pdoftype = doftype;
-	return numdofs;
-}/*}}}*/
-void StressbalanceAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
-
-	/*Intermediaries*/
-	int     fe_FS;
-	int     numoutputs;
-	char**  requestedoutputs = NULL;
-	int     materials_type;
-
-	parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsSIAEnum));
-	parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsSSAEnum));
-	parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsL1L2Enum));
-	parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsHOEnum));
-	parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsFSEnum));
-	parameters->AddObject(iomodel->CopyConstantObject(FlowequationFeFSEnum));
-	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceRestolEnum));
-	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceReltolEnum));
-	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceAbstolEnum));
-	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceIsnewtonEnum));
-	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceMaxiterEnum));
-	parameters->AddObject(iomodel->CopyConstantObject(StressbalancePenaltyFactorEnum));
-	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceRiftPenaltyThresholdEnum));
-	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceFSreconditioningEnum));
-	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceShelfDampeningEnum));
-	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceViscosityOvershootEnum));
-	parameters->AddObject(iomodel->CopyConstantObject(FrictionLawEnum));
-
-	/*XTH LATH parameters*/
-	iomodel->Constant(&fe_FS,FlowequationFeFSEnum);
-	if(fe_FS==XTaylorHoodEnum || fe_FS==LATaylorHoodEnum || fe_FS==LACrouzeixRaviartEnum){
-		parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianREnum));
-		parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianRlambdaEnum));
-		parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianThetaEnum));
-	}
-
-	iomodel->Constant(&materials_type,MaterialsEnum);
-	if(materials_type==MatdamageiceEnum){
-		parameters->AddObject(iomodel->CopyConstantObject(DamageC1Enum));
-		parameters->AddObject(iomodel->CopyConstantObject(DamageStressThresholdEnum));
-	}
-
-	/*Requested outputs*/
-	iomodel->FetchData(&requestedoutputs,&numoutputs,StressbalanceRequestedOutputsEnum);
-	parameters->AddObject(new IntParam(StressbalanceNumRequestedOutputsEnum,numoutputs));
-	if(numoutputs)parameters->AddObject(new StringArrayParam(StressbalanceRequestedOutputsEnum,requestedoutputs,numoutputs));
-	iomodel->DeleteData(&requestedoutputs,numoutputs,StressbalanceRequestedOutputsEnum);
-
-	/*Deal with friction parameters*/
-	int frictionlaw;
-	iomodel->Constant(&frictionlaw,FrictionLawEnum);
-	if(frictionlaw==4 || frictionlaw==6) parameters->AddObject(iomodel->CopyConstantObject(FrictionGammaEnum));
-
-}/*}}}*/
-void StressbalanceAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
-
-	/*Intermediaries*/
-	int    materials_type,finiteelement,fe_FS;
-	int    approximation,frictionlaw;
-	int*   finiteelement_list=NULL;
-	bool   isSSA,isL1L2,isHO,isFS,iscoupling;
-	bool   control_analysis;
-	bool   dakota_analysis;
-	bool   islevelset;
-	bool   isdamage;
-
-	/*Fetch constants needed: */
-	iomodel->Constant(&isSSA,FlowequationIsSSAEnum);
-	iomodel->Constant(&isL1L2,FlowequationIsL1L2Enum);
-	iomodel->Constant(&isHO,FlowequationIsHOEnum);
-	iomodel->Constant(&isFS,FlowequationIsFSEnum);
-	iomodel->Constant(&control_analysis,InversionIscontrolEnum);
-	iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum);
-	iomodel->Constant(&materials_type,MaterialsEnum);
-	iomodel->Constant(&islevelset,TransientIslevelsetEnum);
-	iomodel->Constant(&frictionlaw,FrictionLawEnum);
-
-	/*return if no processing required*/
-	if(!isSSA & !isL1L2 & !isHO & !isFS) return;
-
-	/*Fetch data needed and allocate vectors: */
-	iomodel->FetchData(1,FlowequationElementEquationEnum);
-	finiteelement_list=xNewZeroInit<int>(iomodel->numberofelements);
-
-	/*Do we have coupling*/
-	if( (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.)
-	 iscoupling = true;
-	else
-	 iscoupling = false;
-
-	/*is damage mechanics being used?*/
-	if(materials_type==MaticeEnum) isdamage = false;
-	else if(materials_type==MatdamageiceEnum) isdamage = true;
-	else _error_("Material type not recognized");
-
-	/*Get finite element type*/
-	if(!iscoupling){
-		if(isSSA)       iomodel->Constant(&finiteelement,FlowequationFeSSAEnum);
-		else if(isL1L2) finiteelement = P1Enum;
-		else if(isHO)   iomodel->Constant(&finiteelement,FlowequationFeHOEnum);
-		else if(isFS)   iomodel->Constant(&finiteelement,FlowequationFeFSEnum);
-		for(int i=0;i<iomodel->numberofelements;i++){
-			finiteelement_list[i]=finiteelement;
-		}
-	}
-	else{
-		if(isFS){
-			for(int i=0;i<iomodel->numberofelements;i++){
-				approximation=reCast<int>(iomodel->Data(FlowequationElementEquationEnum)[i]);
-				if(approximation==FSApproximationEnum || approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){
-					finiteelement_list[i]=MINIcondensedEnum;
-				}
-				else{
-					finiteelement_list[i]=P1Enum;
-				}
-			}
-		}
-		else{
-			finiteelement = P1Enum;
-			for(int i=0;i<iomodel->numberofelements;i++){
-				finiteelement_list[i]=finiteelement;
-			}
-		}
-	}
-
-	/*Update elements: */
-	int counter=0;
-	for(int i=0;i<iomodel->numberofelements;i++){
-		if(iomodel->my_elements[i]){
-			Element* element=(Element*)elements->GetObjectByOffset(counter);
-			element->Update(i,iomodel,analysis_counter,analysis_type,finiteelement_list[i]);
-			counter++;
-		}
-	}
-
-	/*Create inputs: */
-	iomodel->FetchDataToInput(elements,ThicknessEnum);
-	iomodel->FetchDataToInput(elements,SurfaceEnum);
-	iomodel->FetchDataToInput(elements,BaseEnum);
-	iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
-	iomodel->FetchDataToInput(elements,MaskGroundediceLevelsetEnum);
-	iomodel->FetchDataToInput(elements,MaterialsRheologyBEnum);
-	iomodel->FetchDataToInput(elements,MaterialsRheologyNEnum);
-	iomodel->FetchDataToInput(elements,VxEnum,0.);
-	iomodel->FetchDataToInput(elements,VyEnum,0.);
-	iomodel->FetchDataToInput(elements,LoadingforceXEnum);
-	iomodel->FetchDataToInput(elements,LoadingforceYEnum);
-	#ifdef LATERALFRICTION
-	iomodel->FetchDataToInput(elements,MeshVertexonboundaryEnum);
-	#endif
-	if(isdamage)iomodel->FetchDataToInput(elements,DamageDEnum);
-
-	if(iomodel->domaintype!=Domain2DhorizontalEnum){
-		iomodel->FetchDataToInput(elements,MeshVertexonbaseEnum);
-		iomodel->FetchDataToInput(elements,MeshVertexonsurfaceEnum);
-	}
-	if(iomodel->domaintype==Domain3DEnum){
-		iomodel->FetchDataToInput(elements,FlowequationBorderFSEnum);
-		iomodel->FetchDataToInput(elements,LoadingforceZEnum);
-		iomodel->FetchDataToInput(elements,VzEnum,0.);
-	}
-	if(isFS){
-		iomodel->FetchDataToInput(elements,PressureEnum,0.);
-		iomodel->FetchDataToInput(elements,BasalforcingsFloatingiceMeltingRateEnum,0.);
-	}
-	if(islevelset){
-		iomodel->FetchDataToInput(elements,IceMaskNodeActivationEnum);
-	}
-	/*LATH parameters*/
-	iomodel->Constant(&fe_FS,FlowequationFeFSEnum);
-	if(fe_FS==LATaylorHoodEnum || fe_FS==LACrouzeixRaviartEnum){
-		InputUpdateFromConstantx(elements,0.,SigmaNNEnum);
-	}
-
-	/*Friction law variables*/
-	switch(frictionlaw){
-		case 1:
-			iomodel->FetchDataToInput(elements,FrictionCoefficientEnum);
-			iomodel->FetchDataToInput(elements,FrictionPEnum);
-			iomodel->FetchDataToInput(elements,FrictionQEnum);
-			break;
-		case 2:
-			iomodel->FetchDataToInput(elements,FrictionCEnum);
-			iomodel->FetchDataToInput(elements,FrictionMEnum);
-			break;
-		case 3:
-			iomodel->FetchDataToInput(elements,FrictionCEnum);
-			iomodel->FetchDataToInput(elements,FrictionAsEnum);
-			iomodel->FetchDataToInput(elements,FrictionQEnum);
-			iomodel->FetchDataToInput(elements,FrictionEffectivePressureEnum);
-			break;
-		case 4:
-			iomodel->FetchDataToInput(elements,FrictionCoefficientEnum);
-			iomodel->FetchDataToInput(elements,FrictionPEnum);
-			iomodel->FetchDataToInput(elements,FrictionQEnum);
-			iomodel->FetchDataToInput(elements,PressureEnum);
-			iomodel->FetchDataToInput(elements,TemperatureEnum);
-			break;
-		case 5:
-			iomodel->FetchDataToInput(elements,FrictionCoefficientEnum);
-			iomodel->FetchDataToInput(elements,FrictionPEnum);
-			iomodel->FetchDataToInput(elements,FrictionQEnum);
-			iomodel->FetchDataToInput(elements,FrictionWaterLayerEnum);
-			break;
-		case 6:
-			iomodel->FetchDataToInput(elements,FrictionCEnum);
-			iomodel->FetchDataToInput(elements,FrictionMEnum);
-			iomodel->FetchDataToInput(elements,PressureEnum);
-			iomodel->FetchDataToInput(elements,TemperatureEnum);
-			break;
-		default:
-			_error_("not supported");
-	}
-
-#ifdef _HAVE_ANDROID_
-	elements->InputDuplicate(FrictionCoefficientEnum,AndroidFrictionCoefficientEnum);
-#endif
-
-	/*Free data: */
-	iomodel->DeleteData(1,FlowequationElementEquationEnum);
-	xDelete<int>(finiteelement_list);
-}/*}}}*/
-void StressbalanceAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
-
-	/*Intermediary*/
-	bool isSSA,isL1L2,isHO,isFS,iscoupling;
-	int  finiteelement=-1,approximation=-1;
-
-	/*Fetch parameters: */
-	iomodel->Constant(&isSSA,FlowequationIsSSAEnum);
-	iomodel->Constant(&isL1L2,FlowequationIsL1L2Enum);
-	iomodel->Constant(&isHO,FlowequationIsHOEnum);
-	iomodel->Constant(&isFS,FlowequationIsFSEnum);
-
-	/*Now, check that we have non SIA elements */
-	if(!isSSA & !isL1L2 & !isHO & !isFS) return;
-
-	/*Do we have coupling*/
-	if( (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.)
-	 iscoupling = true;
-	else
-	 iscoupling = false;
-
-	/*If no coupling, call Regular CreateNodes, else, use P1 elements only*/
-	if(!iscoupling){
-
-		/*Get finite element type*/
-		if(isSSA){
-			approximation=SSAApproximationEnum;
-			iomodel->Constant(&finiteelement,FlowequationFeSSAEnum);
-		}
-		else if(isL1L2){
-			approximation = L1L2ApproximationEnum;
-			finiteelement = P1Enum;
-		}
-		else if(isHO){
-			approximation = HOApproximationEnum;
-			iomodel->Constant(&finiteelement,FlowequationFeHOEnum);
-		}
-		else if(isFS){
-			approximation = FSApproximationEnum;
-			iomodel->Constant(&finiteelement,FlowequationFeFSEnum);
-		}
-		iomodel->FetchData(3,FlowequationBorderSSAEnum,FlowequationVertexEquationEnum,StressbalanceReferentialEnum);
-		if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(3,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum,FlowequationBorderFSEnum);
-		::CreateNodes(nodes,iomodel,StressbalanceAnalysisEnum,finiteelement,approximation);
-		iomodel->DeleteData(6,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum,
-					FlowequationVertexEquationEnum,StressbalanceReferentialEnum);
-	}
-	else{
-		/*Coupling: we are going to create P1 Elements only*/
-
-		Node*  node  = NULL;
-		int    lid=0;
-		if(!nodes) nodes = new Nodes();
-
-		iomodel->FetchData(6,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum,
-					FlowequationVertexEquationEnum,StressbalanceReferentialEnum);
-		if(isFS){
-			/*P1+ velocity*/
-			for(int i=0;i<iomodel->numberofvertices;i++){
-				if(iomodel->my_vertices[i]){
-					approximation=reCast<int>(iomodel->Data(FlowequationVertexEquationEnum)[i]);
-					if(approximation==FSApproximationEnum)  approximation=FSvelocityEnum;
-					nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,lid++,i,iomodel,StressbalanceAnalysisEnum,approximation));
-				}
-			}
-			for(int i=0;i<iomodel->numberofelements;i++){
-				if(iomodel->my_elements[i]){
-					node = new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,iomodel,StressbalanceAnalysisEnum,FSvelocityEnum);
-					node->Deactivate();
-					nodes->AddObject(node);
-				}
-			}
-			/*P1 pressure*/
-			for(int i=0;i<iomodel->numberofvertices;i++){
-				if(iomodel->my_vertices[i]){
-					approximation=reCast<int>(iomodel->Data(FlowequationVertexEquationEnum)[i]);
-					node = new Node(iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+i+1,iomodel->numberofvertices+iomodel->numberofelements+i,lid++,i,iomodel,StressbalanceAnalysisEnum,FSpressureEnum);
-					if(approximation==HOApproximationEnum || approximation==SSAApproximationEnum){
-						node->Deactivate();
-					}
-					nodes->AddObject(node);
-				}
-			}
-		}
-		else{
-			for(int i=0;i<iomodel->numberofvertices;i++){
-				if(iomodel->my_vertices[i]){
-					nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,lid++,i,iomodel,StressbalanceAnalysisEnum,reCast<int>(iomodel->Data(FlowequationVertexEquationEnum)[i])));
-				}
-			}
-		}
-		iomodel->DeleteData(6,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum,
-					FlowequationVertexEquationEnum,StressbalanceReferentialEnum);
-	}
-}/*}}}*/
 void StressbalanceAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/
 
 	/*Intermediary*/
@@ -887,7 +499,395 @@
 		xDelete<IssmDouble>(riftinfo);
 	}
 }/*}}}*/
+void StressbalanceAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel){/*{{{*/
 
+	/*Intermediary*/
+	bool isSSA,isL1L2,isHO,isFS,iscoupling;
+	int  finiteelement=-1,approximation=-1;
+
+	/*Fetch parameters: */
+	iomodel->Constant(&isSSA,FlowequationIsSSAEnum);
+	iomodel->Constant(&isL1L2,FlowequationIsL1L2Enum);
+	iomodel->Constant(&isHO,FlowequationIsHOEnum);
+	iomodel->Constant(&isFS,FlowequationIsFSEnum);
+
+	/*Now, check that we have non SIA elements */
+	if(!isSSA & !isL1L2 & !isHO & !isFS) return;
+
+	/*Do we have coupling*/
+	if( (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.)
+	 iscoupling = true;
+	else
+	 iscoupling = false;
+
+	/*If no coupling, call Regular CreateNodes, else, use P1 elements only*/
+	if(!iscoupling){
+
+		/*Get finite element type*/
+		if(isSSA){
+			approximation=SSAApproximationEnum;
+			iomodel->Constant(&finiteelement,FlowequationFeSSAEnum);
+		}
+		else if(isL1L2){
+			approximation = L1L2ApproximationEnum;
+			finiteelement = P1Enum;
+		}
+		else if(isHO){
+			approximation = HOApproximationEnum;
+			iomodel->Constant(&finiteelement,FlowequationFeHOEnum);
+		}
+		else if(isFS){
+			approximation = FSApproximationEnum;
+			iomodel->Constant(&finiteelement,FlowequationFeFSEnum);
+		}
+		iomodel->FetchData(3,FlowequationBorderSSAEnum,FlowequationVertexEquationEnum,StressbalanceReferentialEnum);
+		if(iomodel->domaintype!=Domain2DhorizontalEnum) iomodel->FetchData(3,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum,FlowequationBorderFSEnum);
+		::CreateNodes(nodes,iomodel,StressbalanceAnalysisEnum,finiteelement,approximation);
+		iomodel->DeleteData(6,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum,
+					FlowequationVertexEquationEnum,StressbalanceReferentialEnum);
+	}
+	else{
+		/*Coupling: we are going to create P1 Elements only*/
+
+		Node*  node  = NULL;
+		int    lid=0;
+		if(!nodes) nodes = new Nodes();
+
+		iomodel->FetchData(6,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum,
+					FlowequationVertexEquationEnum,StressbalanceReferentialEnum);
+		if(isFS){
+			/*P1+ velocity*/
+			for(int i=0;i<iomodel->numberofvertices;i++){
+				if(iomodel->my_vertices[i]){
+					approximation=reCast<int>(iomodel->Data(FlowequationVertexEquationEnum)[i]);
+					if(approximation==FSApproximationEnum)  approximation=FSvelocityEnum;
+					nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,lid++,i,iomodel,StressbalanceAnalysisEnum,approximation));
+				}
+			}
+			for(int i=0;i<iomodel->numberofelements;i++){
+				if(iomodel->my_elements[i]){
+					node = new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,lid++,0,iomodel,StressbalanceAnalysisEnum,FSvelocityEnum);
+					node->Deactivate();
+					nodes->AddObject(node);
+				}
+			}
+			/*P1 pressure*/
+			for(int i=0;i<iomodel->numberofvertices;i++){
+				if(iomodel->my_vertices[i]){
+					approximation=reCast<int>(iomodel->Data(FlowequationVertexEquationEnum)[i]);
+					node = new Node(iomodel->nodecounter+iomodel->numberofvertices+iomodel->numberofelements+i+1,iomodel->numberofvertices+iomodel->numberofelements+i,lid++,i,iomodel,StressbalanceAnalysisEnum,FSpressureEnum);
+					if(approximation==HOApproximationEnum || approximation==SSAApproximationEnum){
+						node->Deactivate();
+					}
+					nodes->AddObject(node);
+				}
+			}
+		}
+		else{
+			for(int i=0;i<iomodel->numberofvertices;i++){
+				if(iomodel->my_vertices[i]){
+					nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,lid++,i,iomodel,StressbalanceAnalysisEnum,reCast<int>(iomodel->Data(FlowequationVertexEquationEnum)[i])));
+				}
+			}
+		}
+		iomodel->DeleteData(6,MeshVertexonbaseEnum,MeshVertexonsurfaceEnum,FlowequationBorderSSAEnum,FlowequationBorderFSEnum,
+					FlowequationVertexEquationEnum,StressbalanceReferentialEnum);
+	}
+}/*}}}*/
+int  StressbalanceAnalysis::DofsPerNode(int** pdoftype,int domaintype,int approximation){/*{{{*/
+
+	/*output*/
+	int *doftype = NULL;
+	int  numdofs;
+
+	switch(approximation){
+		case SSAApproximationEnum:
+			 switch(domaintype){
+				 case Domain3DEnum:           numdofs=2; break;
+				 case Domain2DhorizontalEnum: numdofs=2; break;
+				 case Domain2DverticalEnum:   numdofs=1; break;
+				 default: _error_("mesh type not supported yet");
+			 }
+			 break;
+		case L1L2ApproximationEnum: numdofs =2; break;
+		case HOApproximationEnum:   
+			 switch(domaintype){
+				 case Domain3DEnum:         numdofs=2; break;
+				 case Domain2DverticalEnum: numdofs=1; break;
+				 default: _error_("mesh type not supported yet");
+			 }
+			 break;
+		case SIAApproximationEnum:  numdofs =2; break;
+		case FSvelocityEnum:
+			 switch(domaintype){
+				 case Domain3DEnum:         numdofs=3; break;
+				 case Domain2DverticalEnum: numdofs=2; break;
+				 default: _error_("mesh type not supported yet");
+			}
+			break;
+		case FSpressureEnum: numdofs=1; break;
+		case NoneApproximationEnum:
+			 switch(domaintype){
+				 case Domain3DEnum:         numdofs=4; break;
+				 case Domain2DverticalEnum: numdofs=3; break;
+				 default: _error_("mesh type not supported yet");
+			}
+			break;
+		case SSAHOApproximationEnum:
+			numdofs=4;
+			doftype=xNew<int>(numdofs);
+			doftype[0]=SSAApproximationEnum;
+			doftype[1]=SSAApproximationEnum;
+			doftype[2]=HOApproximationEnum;
+			doftype[3]=HOApproximationEnum;
+			break;
+		case HOFSApproximationEnum:
+			numdofs=5;
+			doftype=xNew<int>(numdofs);
+			doftype[0]=HOApproximationEnum;
+			doftype[1]=HOApproximationEnum;
+			doftype[2]=FSvelocityEnum;
+			doftype[3]=FSvelocityEnum;
+			doftype[4]=FSvelocityEnum;
+			break;
+		case SSAFSApproximationEnum:
+			numdofs=5;
+			doftype=xNew<int>(numdofs);
+			doftype[0]=SSAApproximationEnum;
+			doftype[1]=SSAApproximationEnum;
+			doftype[2]=FSvelocityEnum;
+			doftype[3]=FSvelocityEnum;
+			doftype[4]=FSvelocityEnum;
+			break;
+		default:
+			_error_("Approximation " << EnumToStringx(approximation) << " not implemented yet");
+	}
+
+	/*Assign output pointer and return*/
+	*pdoftype = doftype;
+	return numdofs;
+}/*}}}*/
+void StressbalanceAnalysis::UpdateElements(Elements* elements,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/
+
+	/*Intermediaries*/
+	int    materials_type,finiteelement,fe_FS;
+	int    approximation,frictionlaw;
+	int*   finiteelement_list=NULL;
+	bool   isSSA,isL1L2,isHO,isFS,iscoupling;
+	bool   control_analysis;
+	bool   dakota_analysis;
+	bool   islevelset;
+	bool   isdamage;
+
+	/*Fetch constants needed: */
+	iomodel->Constant(&isSSA,FlowequationIsSSAEnum);
+	iomodel->Constant(&isL1L2,FlowequationIsL1L2Enum);
+	iomodel->Constant(&isHO,FlowequationIsHOEnum);
+	iomodel->Constant(&isFS,FlowequationIsFSEnum);
+	iomodel->Constant(&control_analysis,InversionIscontrolEnum);
+	iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum);
+	iomodel->Constant(&materials_type,MaterialsEnum);
+	iomodel->Constant(&islevelset,TransientIslevelsetEnum);
+	iomodel->Constant(&frictionlaw,FrictionLawEnum);
+
+	/*return if no processing required*/
+	if(!isSSA & !isL1L2 & !isHO & !isFS) return;
+
+	/*Fetch data needed and allocate vectors: */
+	iomodel->FetchData(1,FlowequationElementEquationEnum);
+	finiteelement_list=xNewZeroInit<int>(iomodel->numberofelements);
+
+	/*Do we have coupling*/
+	if( (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.)
+	 iscoupling = true;
+	else
+	 iscoupling = false;
+
+	/*is damage mechanics being used?*/
+	if(materials_type==MaticeEnum) isdamage = false;
+	else if(materials_type==MatdamageiceEnum) isdamage = true;
+	else _error_("Material type not recognized");
+
+	/*Get finite element type*/
+	if(!iscoupling){
+		if(isSSA)       iomodel->Constant(&finiteelement,FlowequationFeSSAEnum);
+		else if(isL1L2) finiteelement = P1Enum;
+		else if(isHO)   iomodel->Constant(&finiteelement,FlowequationFeHOEnum);
+		else if(isFS)   iomodel->Constant(&finiteelement,FlowequationFeFSEnum);
+		for(int i=0;i<iomodel->numberofelements;i++){
+			finiteelement_list[i]=finiteelement;
+		}
+	}
+	else{
+		if(isFS){
+			for(int i=0;i<iomodel->numberofelements;i++){
+				approximation=reCast<int>(iomodel->Data(FlowequationElementEquationEnum)[i]);
+				if(approximation==FSApproximationEnum || approximation==HOFSApproximationEnum || approximation==SSAFSApproximationEnum){
+					finiteelement_list[i]=MINIcondensedEnum;
+				}
+				else{
+					finiteelement_list[i]=P1Enum;
+				}
+			}
+		}
+		else{
+			finiteelement = P1Enum;
+			for(int i=0;i<iomodel->numberofelements;i++){
+				finiteelement_list[i]=finiteelement;
+			}
+		}
+	}
+
+	/*Update elements: */
+	int counter=0;
+	for(int i=0;i<iomodel->numberofelements;i++){
+		if(iomodel->my_elements[i]){
+			Element* element=(Element*)elements->GetObjectByOffset(counter);
+			element->Update(i,iomodel,analysis_counter,analysis_type,finiteelement_list[i]);
+			counter++;
+		}
+	}
+
+	/*Create inputs: */
+	iomodel->FetchDataToInput(elements,ThicknessEnum);
+	iomodel->FetchDataToInput(elements,SurfaceEnum);
+	iomodel->FetchDataToInput(elements,BaseEnum);
+	iomodel->FetchDataToInput(elements,MaskIceLevelsetEnum);
+	iomodel->FetchDataToInput(elements,MaskGroundediceLevelsetEnum);
+	iomodel->FetchDataToInput(elements,MaterialsRheologyBEnum);
+	iomodel->FetchDataToInput(elements,MaterialsRheologyNEnum);
+	iomodel->FetchDataToInput(elements,VxEnum,0.);
+	iomodel->FetchDataToInput(elements,VyEnum,0.);
+	iomodel->FetchDataToInput(elements,LoadingforceXEnum);
+	iomodel->FetchDataToInput(elements,LoadingforceYEnum);
+	#ifdef LATERALFRICTION
+	iomodel->FetchDataToInput(elements,MeshVertexonboundaryEnum);
+	#endif
+	if(isdamage)iomodel->FetchDataToInput(elements,DamageDEnum);
+
+	if(iomodel->domaintype!=Domain2DhorizontalEnum){
+		iomodel->FetchDataToInput(elements,MeshVertexonbaseEnum);
+		iomodel->FetchDataToInput(elements,MeshVertexonsurfaceEnum);
+	}
+	if(iomodel->domaintype==Domain3DEnum){
+		iomodel->FetchDataToInput(elements,FlowequationBorderFSEnum);
+		iomodel->FetchDataToInput(elements,LoadingforceZEnum);
+		iomodel->FetchDataToInput(elements,VzEnum,0.);
+	}
+	if(isFS){
+		iomodel->FetchDataToInput(elements,PressureEnum,0.);
+		iomodel->FetchDataToInput(elements,BasalforcingsFloatingiceMeltingRateEnum,0.);
+	}
+	if(islevelset){
+		iomodel->FetchDataToInput(elements,IceMaskNodeActivationEnum);
+	}
+	/*LATH parameters*/
+	iomodel->Constant(&fe_FS,FlowequationFeFSEnum);
+	if(fe_FS==LATaylorHoodEnum || fe_FS==LACrouzeixRaviartEnum){
+		InputUpdateFromConstantx(elements,0.,SigmaNNEnum);
+	}
+
+	/*Friction law variables*/
+	switch(frictionlaw){
+		case 1:
+			iomodel->FetchDataToInput(elements,FrictionCoefficientEnum);
+			iomodel->FetchDataToInput(elements,FrictionPEnum);
+			iomodel->FetchDataToInput(elements,FrictionQEnum);
+			break;
+		case 2:
+			iomodel->FetchDataToInput(elements,FrictionCEnum);
+			iomodel->FetchDataToInput(elements,FrictionMEnum);
+			break;
+		case 3:
+			iomodel->FetchDataToInput(elements,FrictionCEnum);
+			iomodel->FetchDataToInput(elements,FrictionAsEnum);
+			iomodel->FetchDataToInput(elements,FrictionQEnum);
+			iomodel->FetchDataToInput(elements,FrictionEffectivePressureEnum);
+			break;
+		case 4:
+			iomodel->FetchDataToInput(elements,FrictionCoefficientEnum);
+			iomodel->FetchDataToInput(elements,FrictionPEnum);
+			iomodel->FetchDataToInput(elements,FrictionQEnum);
+			iomodel->FetchDataToInput(elements,PressureEnum);
+			iomodel->FetchDataToInput(elements,TemperatureEnum);
+			break;
+		case 5:
+			iomodel->FetchDataToInput(elements,FrictionCoefficientEnum);
+			iomodel->FetchDataToInput(elements,FrictionPEnum);
+			iomodel->FetchDataToInput(elements,FrictionQEnum);
+			iomodel->FetchDataToInput(elements,FrictionWaterLayerEnum);
+			break;
+		case 6:
+			iomodel->FetchDataToInput(elements,FrictionCEnum);
+			iomodel->FetchDataToInput(elements,FrictionMEnum);
+			iomodel->FetchDataToInput(elements,PressureEnum);
+			iomodel->FetchDataToInput(elements,TemperatureEnum);
+			break;
+		default:
+			_error_("not supported");
+	}
+
+#ifdef _HAVE_ANDROID_
+	elements->InputDuplicate(FrictionCoefficientEnum,AndroidFrictionCoefficientEnum);
+#endif
+
+	/*Free data: */
+	iomodel->DeleteData(1,FlowequationElementEquationEnum);
+	xDelete<int>(finiteelement_list);
+}/*}}}*/
+void StressbalanceAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/
+
+	/*Intermediaries*/
+	int     fe_FS;
+	int     numoutputs;
+	char**  requestedoutputs = NULL;
+	int     materials_type;
+
+	parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsSIAEnum));
+	parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsSSAEnum));
+	parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsL1L2Enum));
+	parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsHOEnum));
+	parameters->AddObject(iomodel->CopyConstantObject(FlowequationIsFSEnum));
+	parameters->AddObject(iomodel->CopyConstantObject(FlowequationFeFSEnum));
+	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceRestolEnum));
+	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceReltolEnum));
+	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceAbstolEnum));
+	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceIsnewtonEnum));
+	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceMaxiterEnum));
+	parameters->AddObject(iomodel->CopyConstantObject(StressbalancePenaltyFactorEnum));
+	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceRiftPenaltyThresholdEnum));
+	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceFSreconditioningEnum));
+	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceShelfDampeningEnum));
+	parameters->AddObject(iomodel->CopyConstantObject(StressbalanceViscosityOvershootEnum));
+	parameters->AddObject(iomodel->CopyConstantObject(FrictionLawEnum));
+
+	/*XTH LATH parameters*/
+	iomodel->Constant(&fe_FS,FlowequationFeFSEnum);
+	if(fe_FS==XTaylorHoodEnum || fe_FS==LATaylorHoodEnum || fe_FS==LACrouzeixRaviartEnum){
+		parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianREnum));
+		parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianRlambdaEnum));
+		parameters->AddObject(iomodel->CopyConstantObject(AugmentedLagrangianThetaEnum));
+	}
+
+	iomodel->Constant(&materials_type,MaterialsEnum);
+	if(materials_type==MatdamageiceEnum){
+		parameters->AddObject(iomodel->CopyConstantObject(DamageC1Enum));
+		parameters->AddObject(iomodel->CopyConstantObject(DamageStressThresholdEnum));
+	}
+
+	/*Requested outputs*/
+	iomodel->FetchData(&requestedoutputs,&numoutputs,StressbalanceRequestedOutputsEnum);
+	parameters->AddObject(new IntParam(StressbalanceNumRequestedOutputsEnum,numoutputs));
+	if(numoutputs)parameters->AddObject(new StringArrayParam(StressbalanceRequestedOutputsEnum,requestedoutputs,numoutputs));
+	iomodel->DeleteData(&requestedoutputs,numoutputs,StressbalanceRequestedOutputsEnum);
+
+	/*Deal with friction parameters*/
+	int frictionlaw;
+	iomodel->Constant(&frictionlaw,FrictionLawEnum);
+	if(frictionlaw==4 || frictionlaw==6) parameters->AddObject(iomodel->CopyConstantObject(FrictionGammaEnum));
+
+}/*}}}*/
+
 /*Finite Element Analysis*/
 void           StressbalanceAnalysis::Core(FemModel* femmodel){/*{{{*/
 
@@ -1026,7 +1026,7 @@
 			_error_("Approximation "<<EnumToStringx(approximation)<<" not supported");
 	}
 }/*}}}*/
-void StressbalanceAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
+void           StressbalanceAnalysis::GetSolutionFromInputs(Vector<IssmDouble>* solution,Element* element){/*{{{*/
 
 	int approximation;
 	element->GetInputValue(&approximation,ApproximationEnum);
@@ -1047,7 +1047,7 @@
 			_error_("Approximation "<<EnumToStringx(approximation)<<" not supported");
 	}
 }/*}}}*/
-void StressbalanceAnalysis::GetSolutionFromInputsHoriz(Vector<IssmDouble>* solution,Element* element){/*{{{*/
+void           StressbalanceAnalysis::GetSolutionFromInputsHoriz(Vector<IssmDouble>* solution,Element* element){/*{{{*/
 
 	IssmDouble   vx,vy;
 	int          domaintype,dim,approximation,dofpernode;
@@ -1097,10 +1097,10 @@
 	xDelete<IssmDouble>(values);
 	xDelete<int>(doflist);
 }/*}}}*/
-void StressbalanceAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
+void           StressbalanceAnalysis::GradientJ(Vector<IssmDouble>* gradient,Element* element,int control_type,int control_index){/*{{{*/
 	_error_("Not implemented yet");
 }/*}}}*/
-void StressbalanceAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
+void           StressbalanceAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/
 
 	int approximation;
 	element->GetInputValue(&approximation,ApproximationEnum);
@@ -1132,7 +1132,7 @@
 			_error_("Approximation "<<EnumToStringx(approximation)<<" not supported");
 	}
 }/*}}}*/
-void StressbalanceAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
+void           StressbalanceAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/
 	/*Default, do nothing*/
 	bool islevelset;
 	femmodel->parameters->FindParam(&islevelset,TransientIslevelsetEnum);
@@ -1659,7 +1659,7 @@
 	delete gauss;
 	return pe;
 }/*}}}*/
-void StressbalanceAnalysis::GetBSSA(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void           StressbalanceAnalysis::GetBSSA(IssmDouble* B,Element* element,int dim,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: 
@@ -1699,7 +1699,7 @@
 	/*Clean-up*/
 	xDelete<IssmDouble>(dbasis);
 }/*}}}*/
-void StressbalanceAnalysis::GetBSSAFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void           StressbalanceAnalysis::GetBSSAFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 	/*Compute B  matrix. B=[B1 B2 B3] where Bi is square and of size 2. 
 	 * For node i, Bi can be expressed in the actual coordinate system
 	 * by: 
@@ -1736,7 +1736,7 @@
 	/*Clean-up*/
 	xDelete<IssmDouble>(basis);
 }/*}}}*/
-void StressbalanceAnalysis::GetBSSAprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void           StressbalanceAnalysis::GetBSSAprime(IssmDouble* Bprime,Element* element,int dim,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: 
@@ -1776,7 +1776,7 @@
 	/*Clean-up*/
 	xDelete<IssmDouble>(dbasis);
 }/*}}}*/
-void StressbalanceAnalysis::InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element){/*{{{*/
+void           StressbalanceAnalysis::InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element){/*{{{*/
 
 	int         i,dim,domaintype;
 	IssmDouble  rho_ice,g;
@@ -2110,7 +2110,7 @@
 	delete gauss;
 	return pe;
 }/*}}}*/
-void StressbalanceAnalysis::InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element){/*{{{*/
+void           StressbalanceAnalysis::InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element){/*{{{*/
 
 	int         i,dim,domaintype;
 	IssmDouble  rho_ice,g;
@@ -2285,78 +2285,6 @@
 	delete Ke2;
 	return Ke;
 }/*}}}*/
-ElementMatrix* StressbalanceAnalysis::CreateKMatrixHOViscous(Element* element){/*{{{*/
-
-	/* Check if ice in element */
-	if(!element->IsIceInElement()) return NULL;
-
-	/*Intermediaries*/
-	int         dim,bsize;
-	IssmDouble  viscosity,newviscosity,oldviscosity;
-	IssmDouble  viscosity_overshoot,thickness,Jdet;
-	IssmDouble  D_scalar;
-	IssmDouble *xyz_list = NULL;
-
-	/*Get problem dimension*/
-	element->FindParam(&dim,DomainDimensionEnum);
-	if(dim==2) bsize = 2;
-	else       bsize = 5;
-
-	/*Fetch number of nodes and dof for this finite element*/
-	int numnodes = element->GetNumberOfNodes();
-	int numdof   = numnodes*(dim-1);
-
-	/*Initialize Element matrix and vectors*/
-	ElementMatrix* Ke     = element->NewElementMatrix(HOApproximationEnum);
-	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);
-	Input* vx_input    = element->GetInput(VxEnum);       _assert_(vx_input);
-	Input* vxold_input = element->GetInput(VxPicardEnum); _assert_(vxold_input);
-	Input* vy_input    = NULL;
-	Input* vyold_input = NULL;
-	if(dim==3){
-		vy_input=element->GetInput(VyEnum);          _assert_(vy_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,dim,xyz_list,gauss);
-		this->GetBHOprime(Bprime,element,dim,xyz_list,gauss);
-
-		element->ViscosityHO(&viscosity,dim,xyz_list,gauss,vx_input,vy_input);
-		element->ViscosityHO(&oldviscosity,dim,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<bsize;i++) D[i*bsize+i]=D_scalar;
-
-		TripleMultiply(B,bsize,numdof,1,
-					D,bsize,bsize,0,
-					Bprime,bsize,numdof,0,
-					&Ke->values[0],1);
-	}
-
-	/*Transform Coordinate System*/
-	if(dim==3) 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){/*{{{*/
 
 	/* Check if ice in element */
@@ -2437,6 +2365,78 @@
 	xDelete<IssmDouble>(D);
 	return Ke;
 }/*}}}*/
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixHOViscous(Element* element){/*{{{*/
+
+	/* Check if ice in element */
+	if(!element->IsIceInElement()) return NULL;
+
+	/*Intermediaries*/
+	int         dim,bsize;
+	IssmDouble  viscosity,newviscosity,oldviscosity;
+	IssmDouble  viscosity_overshoot,thickness,Jdet;
+	IssmDouble  D_scalar;
+	IssmDouble *xyz_list = NULL;
+
+	/*Get problem dimension*/
+	element->FindParam(&dim,DomainDimensionEnum);
+	if(dim==2) bsize = 2;
+	else       bsize = 5;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+	int numdof   = numnodes*(dim-1);
+
+	/*Initialize Element matrix and vectors*/
+	ElementMatrix* Ke     = element->NewElementMatrix(HOApproximationEnum);
+	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);
+	Input* vx_input    = element->GetInput(VxEnum);       _assert_(vx_input);
+	Input* vxold_input = element->GetInput(VxPicardEnum); _assert_(vxold_input);
+	Input* vy_input    = NULL;
+	Input* vyold_input = NULL;
+	if(dim==3){
+		vy_input=element->GetInput(VyEnum);          _assert_(vy_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,dim,xyz_list,gauss);
+		this->GetBHOprime(Bprime,element,dim,xyz_list,gauss);
+
+		element->ViscosityHO(&viscosity,dim,xyz_list,gauss,vx_input,vy_input);
+		element->ViscosityHO(&oldviscosity,dim,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<bsize;i++) D[i*bsize+i]=D_scalar;
+
+		TripleMultiply(B,bsize,numdof,1,
+					D,bsize,bsize,0,
+					Bprime,bsize,numdof,0,
+					&Ke->values[0],1);
+	}
+
+	/*Transform Coordinate System*/
+	if(dim==3) 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;
+}/*}}}*/
 #ifdef FSANALYTICAL
 ElementVector* StressbalanceAnalysis::CreatePVectorHO(Element* element){/*{{{*/
 
@@ -2632,7 +2632,7 @@
 	delete gauss;
 	return pe;
 }/*}}}*/
-void StressbalanceAnalysis::GetBHO(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void           StressbalanceAnalysis::GetBHO(IssmDouble* B,Element* element,int dim,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: 
@@ -2681,7 +2681,44 @@
 	/*Clean-up*/
 	xDelete<IssmDouble>(dbasis);
 }/*}}}*/
-void StressbalanceAnalysis::GetBHOprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void           StressbalanceAnalysis::GetBHOFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+	/*Compute B  matrix. B=[B1 B2 B3] where Bi is square and of size 2. 
+	 * For node i, Bi can be expressed in the actual coordinate system
+	 * by: 
+	 *                       3D           2D
+	 *                 Bi=[ N   0 ]    Bi=N
+	 *                    [ 0   N ]
+	 * where N is the finiteelement function for node i.
+	 *
+	 * We assume B has been allocated already, of size: 2 x (numdof*numnodes)
+	 */
+
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = element->GetNumberOfNodes();
+
+	/*Get nodal functions derivatives*/
+	IssmDouble* basis=xNew<IssmDouble>(numnodes);
+	element->NodalFunctions(basis,gauss);
+
+	/*Build L: */
+	if(dim==3){
+		for(int i=0;i<numnodes;i++){
+			B[2*numnodes*0+2*i+0] = basis[i];
+			B[2*numnodes*0+2*i+1] = 0.;
+			B[2*numnodes*1+2*i+0] = 0.;
+			B[2*numnodes*1+2*i+1] = basis[i];
+		}
+	}
+	else{
+		for(int i=0;i<numnodes;i++){
+			B[i] = basis[i];
+		}
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(basis);
+}/*}}}*/
+void           StressbalanceAnalysis::GetBHOprime(IssmDouble* Bprime,Element* element,int dim,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: 
@@ -2726,45 +2763,8 @@
 	/*Clean-up*/
 	xDelete<IssmDouble>(dbasis);
 }/*}}}*/
-void StressbalanceAnalysis::GetBHOFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B  matrix. B=[B1 B2 B3] where Bi is square and of size 2. 
-	 * For node i, Bi can be expressed in the actual coordinate system
-	 * by: 
-	 *                       3D           2D
-	 *                 Bi=[ N   0 ]    Bi=N
-	 *                    [ 0   N ]
-	 * where N is the finiteelement function for node i.
-	 *
-	 * We assume B has been allocated already, of size: 2 x (numdof*numnodes)
-	 */
+void           StressbalanceAnalysis::InputUpdateFromSolutionHO(IssmDouble* solution,Element* element){/*{{{*/
 
-	/*Fetch number of nodes for this finite element*/
-	int numnodes = element->GetNumberOfNodes();
-
-	/*Get nodal functions derivatives*/
-	IssmDouble* basis=xNew<IssmDouble>(numnodes);
-	element->NodalFunctions(basis,gauss);
-
-	/*Build L: */
-	if(dim==3){
-		for(int i=0;i<numnodes;i++){
-			B[2*numnodes*0+2*i+0] = basis[i];
-			B[2*numnodes*0+2*i+1] = 0.;
-			B[2*numnodes*1+2*i+0] = 0.;
-			B[2*numnodes*1+2*i+1] = basis[i];
-		}
-	}
-	else{
-		for(int i=0;i<numnodes;i++){
-			B[i] = basis[i];
-		}
-	}
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(basis);
-}/*}}}*/
-void StressbalanceAnalysis::InputUpdateFromSolutionHO(IssmDouble* solution,Element* element){/*{{{*/
-
 	int         i,dim;
 	int*        doflist=NULL;
 	IssmDouble* xyz_list=NULL;
@@ -2981,6 +2981,143 @@
 	delete Ke3;
 	return Ke;
 }/*}}}*/
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSShelf(Element* element){/*{{{*/
+
+	if(!element->IsFloating() || !element->IsOnBase()) return NULL;
+
+	/*If on not water or not FS, skip stiffness: */
+	int approximation,shelf_dampening;
+	element->GetInputValue(&approximation,ApproximationEnum);
+	if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
+	element->FindParam(&shelf_dampening,StressbalanceShelfDampeningEnum);
+	if(shelf_dampening==0) return NULL;
+
+	/*Intermediaries*/
+	bool        mainlyfloating;
+	int         j,i,dim;
+	IssmDouble  Jdet,slope2,scalar,dt;
+	IssmDouble  slope[3];
+	IssmDouble *xyz_list_base = NULL;
+	IssmDouble *xyz_list      = NULL;
+	Gauss*      gauss         = NULL;
+
+	/*Get problem dimension*/
+	element->FindParam(&dim,DomainDimensionEnum);
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int vnumnodes = element->NumberofNodesVelocity();
+	int pnumnodes = element->NumberofNodesPressure();
+	int numdof    = vnumnodes*dim + pnumnodes;
+
+	/*Initialize Element matrix and vectors*/
+	ElementMatrix* Ke = element->NewElementMatrix(FSvelocityEnum);
+	IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinatesBase(&xyz_list_base);
+	element->GetVerticesCoordinates(&xyz_list);
+	element->FindParam(&dt,TimesteppingTimeStepEnum);
+	if(dt==0)   dt=1.e+5;
+	IssmDouble  rho_water     = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
+	IssmDouble  gravity       = element->GetMaterialParameter(ConstantsGEnum);
+	Input*      surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input);
+
+	/* Start  looping on the number of gaussian points: */
+	gauss=element->NewGaussBase(3);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
+		element->NodalFunctionsVelocity(vbasis,gauss);
+		element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
+		if(dim==2) slope2=slope[0]*slope[0];
+		else if(dim==3) slope2=slope[0]*slope[0]+slope[1]*slope[1];
+		scalar  = rho_water*gravity*sqrt(1+slope2)*gauss->weight*Jdet*dt; 
+		for(i=0;i<vnumnodes;i++){
+			for(j=0;j<vnumnodes;j++){
+				Ke->values[numdof*((i+1)*dim-1)+(j+1)*dim-1] += scalar*vbasis[i]*vbasis[j];
+			}
+		}
+	}
+
+	/*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/
+
+	/*Clean up and return*/
+	delete gauss;
+	xDelete<IssmDouble>(xyz_list_base);
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(vbasis);
+	return Ke;
+}/*}}}*/
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSViscous(Element* element){/*{{{*/
+
+	/*Intermediaries*/
+	int         i,dim,epssize;
+	IssmDouble  viscosity,FSreconditioning,Jdet;
+	IssmDouble *xyz_list = NULL;
+
+	/*Get problem dimension*/
+	element->FindParam(&dim,DomainDimensionEnum);
+	if(dim==2) epssize = 3;
+	else       epssize = 6;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int vnumnodes = element->NumberofNodesVelocity();
+	int pnumnodes = element->NumberofNodesPressure();
+	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,dim,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::CreateKMatrixFSViscousLA(Element* element){/*{{{*/
 
 	/*Intermediaries*/
@@ -3180,143 +3317,6 @@
 	xDelete<int>(cs_list);
 	return Ke;
 }/*}}}*/
-ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSViscous(Element* element){/*{{{*/
-
-	/*Intermediaries*/
-	int         i,dim,epssize;
-	IssmDouble  viscosity,FSreconditioning,Jdet;
-	IssmDouble *xyz_list = NULL;
-
-	/*Get problem dimension*/
-	element->FindParam(&dim,DomainDimensionEnum);
-	if(dim==2) epssize = 3;
-	else       epssize = 6;
-
-	/*Fetch number of nodes and dof for this finite element*/
-	int vnumnodes = element->NumberofNodesVelocity();
-	int pnumnodes = element->NumberofNodesPressure();
-	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,dim,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::CreateKMatrixFSShelf(Element* element){/*{{{*/
-
-	if(!element->IsFloating() || !element->IsOnBase()) return NULL;
-
-	/*If on not water or not FS, skip stiffness: */
-	int approximation,shelf_dampening;
-	element->GetInputValue(&approximation,ApproximationEnum);
-	if(approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
-	element->FindParam(&shelf_dampening,StressbalanceShelfDampeningEnum);
-	if(shelf_dampening==0) return NULL;
-
-	/*Intermediaries*/
-	bool        mainlyfloating;
-	int         j,i,dim;
-	IssmDouble  Jdet,slope2,scalar,dt;
-	IssmDouble  slope[3];
-	IssmDouble *xyz_list_base = NULL;
-	IssmDouble *xyz_list      = NULL;
-	Gauss*      gauss         = NULL;
-
-	/*Get problem dimension*/
-	element->FindParam(&dim,DomainDimensionEnum);
-
-	/*Fetch number of nodes and dof for this finite element*/
-	int vnumnodes = element->NumberofNodesVelocity();
-	int pnumnodes = element->NumberofNodesPressure();
-	int numdof    = vnumnodes*dim + pnumnodes;
-
-	/*Initialize Element matrix and vectors*/
-	ElementMatrix* Ke = element->NewElementMatrix(FSvelocityEnum);
-	IssmDouble*    vbasis = xNew<IssmDouble>(vnumnodes);
-
-	/*Retrieve all inputs and parameters*/
-	element->GetVerticesCoordinatesBase(&xyz_list_base);
-	element->GetVerticesCoordinates(&xyz_list);
-	element->FindParam(&dt,TimesteppingTimeStepEnum);
-	if(dt==0)   dt=1.e+5;
-	IssmDouble  rho_water     = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
-	IssmDouble  gravity       = element->GetMaterialParameter(ConstantsGEnum);
-	Input*      surface_input = element->GetInput(SurfaceEnum); _assert_(surface_input);
-
-	/* Start  looping on the number of gaussian points: */
-	gauss=element->NewGaussBase(3);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
-
-		surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
-		element->NodalFunctionsVelocity(vbasis,gauss);
-		element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss);
-		if(dim==2) slope2=slope[0]*slope[0];
-		else if(dim==3) slope2=slope[0]*slope[0]+slope[1]*slope[1];
-		scalar  = rho_water*gravity*sqrt(1+slope2)*gauss->weight*Jdet*dt; 
-		for(i=0;i<vnumnodes;i++){
-			for(j=0;j<vnumnodes;j++){
-				Ke->values[numdof*((i+1)*dim-1)+(j+1)*dim-1] += scalar*vbasis[i]*vbasis[j];
-			}
-		}
-	}
-
-	/*DO NOT Transform Coordinate System: this stiffness matrix is already expressed in tangential coordinates*/
-
-	/*Clean up and return*/
-	delete gauss;
-	xDelete<IssmDouble>(xyz_list_base);
-	xDelete<IssmDouble>(xyz_list);
-	xDelete<IssmDouble>(vbasis);
-	return Ke;
-}/*}}}*/
 #ifdef FSANALYTICAL
 ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSFriction(Element* element){/*{{{*/
 
@@ -3403,81 +3403,6 @@
 	/*clean-up and return*/
 	return pe;
 }/*}}}*/
-ElementVector* StressbalanceAnalysis::CreatePVectorFSViscous(Element* element){/*{{{*/
-
-	int         i,dim,fe_FS;
-	IssmDouble  x_coord,y_coord,z_coord;
-	IssmDouble  Jdet,forcex,forcey,forcez;
-	IssmDouble *xyz_list = NULL;
-
-	element->FindParam(&fe_FS,FlowequationFeFSEnum);
-	element->FindParam(&dim,DomainDimensionEnum);
-
-	/*Fetch number of nodes and dof for this finite element*/
-	int vnumnodes = element->NumberofNodesVelocity();
-	int pnumnodes = element->NumberofNodesPressure();
-
-	/*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);
-
-	/* 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);
-
-		x_coord=element->GetXcoord(xyz_list,gauss);
-		y_coord=element->GetYcoord(xyz_list,gauss);
-		if(dim==3) z_coord=element->GetZcoord(xyz_list,gauss);
-		else z_coord=0.;
-
-		forcex=fx(x_coord,y_coord,z_coord,FSANALYTICAL);
-		forcey=fy(x_coord,y_coord,z_coord,FSANALYTICAL);
-		forcez=fz(x_coord,y_coord,z_coord,FSANALYTICAL);
-
-		for(i=0;i<vnumnodes;i++){
-			pe->values[i*dim+0] += forcex *Jdet*gauss->weight*vbasis[i];
-			pe->values[i*dim+1] += forcey *Jdet*gauss->weight*vbasis[i];
-			if(dim==3) pe->values[i*dim+2] += forcez *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);
-	if(fe_FS==XTaylorHoodEnum){
-		ElementVector* pe2=CreatePVectorFSViscousXTH(element);
-		ElementVector* pe3 = new ElementVector(pe,pe2);
-		delete pe;
-		delete pe2;
-		return pe3;
-	}
-	else if(fe_FS==LATaylorHoodEnum || fe_FS==LACrouzeixRaviartEnum){
-		ElementVector* pe2=CreatePVectorFSViscousLA(element);
-		ElementVector* pe3 = new ElementVector(pe,pe2);
-		delete pe;
-		delete pe2;
-		return pe3;
-	}
-	return pe;
-}/*}}}*/
 ElementVector* StressbalanceAnalysis::CreatePVectorFSFriction(Element* element){/*{{{*/
 
 	if(!element->IsOnBase()) return NULL;
@@ -3588,6 +3513,81 @@
 	xDelete<IssmDouble>(vbasis);
 	return pe;
 }/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorFSViscous(Element* element){/*{{{*/
+
+	int         i,dim,fe_FS;
+	IssmDouble  x_coord,y_coord,z_coord;
+	IssmDouble  Jdet,forcex,forcey,forcez;
+	IssmDouble *xyz_list = NULL;
+
+	element->FindParam(&fe_FS,FlowequationFeFSEnum);
+	element->FindParam(&dim,DomainDimensionEnum);
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int vnumnodes = element->NumberofNodesVelocity();
+	int pnumnodes = element->NumberofNodesPressure();
+
+	/*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);
+
+	/* 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);
+
+		x_coord=element->GetXcoord(xyz_list,gauss);
+		y_coord=element->GetYcoord(xyz_list,gauss);
+		if(dim==3) z_coord=element->GetZcoord(xyz_list,gauss);
+		else z_coord=0.;
+
+		forcex=fx(x_coord,y_coord,z_coord,FSANALYTICAL);
+		forcey=fy(x_coord,y_coord,z_coord,FSANALYTICAL);
+		forcez=fz(x_coord,y_coord,z_coord,FSANALYTICAL);
+
+		for(i=0;i<vnumnodes;i++){
+			pe->values[i*dim+0] += forcex *Jdet*gauss->weight*vbasis[i];
+			pe->values[i*dim+1] += forcey *Jdet*gauss->weight*vbasis[i];
+			if(dim==3) pe->values[i*dim+2] += forcez *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);
+	if(fe_FS==XTaylorHoodEnum){
+		ElementVector* pe2=CreatePVectorFSViscousXTH(element);
+		ElementVector* pe3 = new ElementVector(pe,pe2);
+		delete pe;
+		delete pe2;
+		return pe3;
+	}
+	else if(fe_FS==LATaylorHoodEnum || fe_FS==LACrouzeixRaviartEnum){
+		ElementVector* pe2=CreatePVectorFSViscousLA(element);
+		ElementVector* pe3 = new ElementVector(pe,pe2);
+		delete pe;
+		delete pe2;
+		return pe3;
+	}
+	return pe;
+}/*}}}*/
 #else
 ElementMatrix* StressbalanceAnalysis::CreateKMatrixFSFriction(Element* element){/*{{{*/
 
@@ -3789,6 +3789,169 @@
 	return pe;
 }/*}}}*/
 #endif
+ElementVector* StressbalanceAnalysis::CreatePVectorFSFront(Element* element){/*{{{*/
+
+	/*If no front, return NULL*/
+	if(!element->IsIcefront()) return NULL;
+
+	/*Intermediaries*/
+	int         i,dim;
+	IssmDouble  Jdet,pressure,surface,z;
+	IssmDouble	normal[3];
+	IssmDouble *xyz_list       = NULL;
+	IssmDouble *xyz_list_front = NULL;
+	Gauss      *gauss          = NULL;
+
+	/*Make sure current element is floating*/
+	if(!element->IsFloating()) return NULL;
+
+	/*Get problem dimension*/
+	element->FindParam(&dim,DomainDimensionEnum);
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int vnumnodes   = element->NumberofNodesVelocity();
+	int pnumnodes   = element->NumberofNodesPressure();
+	int numvertices = element->GetNumberOfVertices();
+
+	/*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);
+	element->GetIcefrontCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum);
+	element->NormalSection(&normal[0],xyz_list_front);
+	Input* surface_input  = element->GetInput(SurfaceEnum); _assert_(surface_input);
+	IssmDouble  rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
+	IssmDouble  gravity   = element->GetMaterialParameter(ConstantsGEnum);
+
+	/*Initialize gauss points*/
+	IssmDouble zmax=xyz_list[0*3+(dim-1)]; for(int i=1;i<numvertices;i++) if(xyz_list[i*3+(dim-1)]>zmax) zmax=xyz_list[i*3+(dim-1)];
+	IssmDouble zmin=xyz_list[0*3+(dim-1)]; for(int i=1;i<numvertices;i++) if(xyz_list[i*3+(dim-1)]<zmin) zmin=xyz_list[i*3+(dim-1)];
+	if(zmax>0. && zmin<0.) gauss=element->NewGauss(xyz_list,xyz_list_front,3,30);//refined in vertical because of the sea level discontinuity
+	else                   gauss=element->NewGauss(xyz_list,xyz_list_front,3,3);
+
+	/* Start  looping on the number of gaussian points: */
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss);
+		element->NodalFunctionsVelocity(vbasis,gauss);
+		surface_input->GetInputValue(&surface,gauss);
+		if(dim==3) z=element->GetZcoord(xyz_list,gauss);
+		else       z=element->GetYcoord(xyz_list,gauss);
+		pressure = rho_water*gravity*min(0.,z);//0 if the gaussian point is above water level
+
+		for (int i=0;i<vnumnodes;i++){
+			pe->values[dim*i+0]+= pressure*Jdet*gauss->weight*normal[0]*vbasis[i];
+			pe->values[dim*i+1]+= pressure*Jdet*gauss->weight*normal[1]*vbasis[i];
+			if(dim==3) pe->values[dim*i+2]+= pressure*Jdet*gauss->weight*normal[2]*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);
+	xDelete<IssmDouble>(xyz_list_front);
+	return pe;
+}/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorFSShelf(Element* element){/*{{{*/
+
+	int         i,dim;
+	IssmDouble  Jdet,water_pressure,bed;
+	IssmDouble	normal[3];
+	IssmDouble *xyz_list_base = NULL;
+
+	/*Get basal element*/
+	if(!element->IsOnBase() || !element->IsFloating()) return NULL;
+
+	/*Get problem dimension*/
+	element->FindParam(&dim,DomainDimensionEnum);
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int vnumnodes = element->NumberofNodesVelocity();
+	int pnumnodes = element->NumberofNodesPressure();
+
+	/*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*      base_input=element->GetInput(BaseEnum); _assert_(base_input);
+	IssmDouble  rho_water=element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
+	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.);
+		base_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);
+
+	/* shelf dampening*/
+	int shelf_dampening;
+	element->FindParam(&shelf_dampening,StressbalanceShelfDampeningEnum);
+	if(shelf_dampening) {
+		Input*      mb_input=element->GetInput(BasalforcingsFloatingiceMeltingRateEnum); _assert_(mb_input);
+		IssmDouble dt,mb,normal_b;
+		element->FindParam(&dt,TimesteppingTimeStepEnum);
+		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);
+			if (dim==2) normal_b=normal[1];
+			else if (dim==3) normal_b=sqrt(normal[0]*normal[0]+normal[1]*normal[1]);
+			mb_input->GetInputValue(&mb, gauss);
+			for(i=0;i<vnumnodes;i++){
+				pe->values[i*dim+1] += dt*rho_water*gravity*mb*gauss->weight*Jdet*vbasis[i]*normal_b;
+			}
+		}
+	}
+
+	/*Clean up and return*/
+	delete gauss;
+	xDelete<int>(cs_list);
+	xDelete<IssmDouble>(vbasis);
+	xDelete<IssmDouble>(xyz_list_base);
+	return pe;
+}/*}}}*/
 ElementVector* StressbalanceAnalysis::CreatePVectorFSViscousLA(Element* element){/*{{{*/
 
 	int         i,dim;
@@ -4043,170 +4206,7 @@
 	xDelete<IssmDouble>(tbasis);
 	return pe;
 }/*}}}*/
-ElementVector* StressbalanceAnalysis::CreatePVectorFSShelf(Element* element){/*{{{*/
-
-	int         i,dim;
-	IssmDouble  Jdet,water_pressure,bed;
-	IssmDouble	normal[3];
-	IssmDouble *xyz_list_base = NULL;
-
-	/*Get basal element*/
-	if(!element->IsOnBase() || !element->IsFloating()) return NULL;
-
-	/*Get problem dimension*/
-	element->FindParam(&dim,DomainDimensionEnum);
-
-	/*Fetch number of nodes and dof for this finite element*/
-	int vnumnodes = element->NumberofNodesVelocity();
-	int pnumnodes = element->NumberofNodesPressure();
-
-	/*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*      base_input=element->GetInput(BaseEnum); _assert_(base_input);
-	IssmDouble  rho_water=element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
-	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.);
-		base_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);
-
-	/* shelf dampening*/
-	int shelf_dampening;
-	element->FindParam(&shelf_dampening,StressbalanceShelfDampeningEnum);
-	if(shelf_dampening) {
-		Input*      mb_input=element->GetInput(BasalforcingsFloatingiceMeltingRateEnum); _assert_(mb_input);
-		IssmDouble dt,mb,normal_b;
-		element->FindParam(&dt,TimesteppingTimeStepEnum);
-		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);
-			if (dim==2) normal_b=normal[1];
-			else if (dim==3) normal_b=sqrt(normal[0]*normal[0]+normal[1]*normal[1]);
-			mb_input->GetInputValue(&mb, gauss);
-			for(i=0;i<vnumnodes;i++){
-				pe->values[i*dim+1] += dt*rho_water*gravity*mb*gauss->weight*Jdet*vbasis[i]*normal_b;
-			}
-		}
-	}
-
-	/*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){/*{{{*/
-
-	/*If no front, return NULL*/
-	if(!element->IsIcefront()) return NULL;
-
-	/*Intermediaries*/
-	int         i,dim;
-	IssmDouble  Jdet,pressure,surface,z;
-	IssmDouble	normal[3];
-	IssmDouble *xyz_list       = NULL;
-	IssmDouble *xyz_list_front = NULL;
-	Gauss      *gauss          = NULL;
-
-	/*Make sure current element is floating*/
-	if(!element->IsFloating()) return NULL;
-
-	/*Get problem dimension*/
-	element->FindParam(&dim,DomainDimensionEnum);
-
-	/*Fetch number of nodes and dof for this finite element*/
-	int vnumnodes   = element->NumberofNodesVelocity();
-	int pnumnodes   = element->NumberofNodesPressure();
-	int numvertices = element->GetNumberOfVertices();
-
-	/*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);
-	element->GetIcefrontCoordinates(&xyz_list_front,xyz_list,MaskIceLevelsetEnum);
-	element->NormalSection(&normal[0],xyz_list_front);
-	Input* surface_input  = element->GetInput(SurfaceEnum); _assert_(surface_input);
-	IssmDouble  rho_water = element->GetMaterialParameter(MaterialsRhoSeawaterEnum);
-	IssmDouble  gravity   = element->GetMaterialParameter(ConstantsGEnum);
-
-	/*Initialize gauss points*/
-	IssmDouble zmax=xyz_list[0*3+(dim-1)]; for(int i=1;i<numvertices;i++) if(xyz_list[i*3+(dim-1)]>zmax) zmax=xyz_list[i*3+(dim-1)];
-	IssmDouble zmin=xyz_list[0*3+(dim-1)]; for(int i=1;i<numvertices;i++) if(xyz_list[i*3+(dim-1)]<zmin) zmin=xyz_list[i*3+(dim-1)];
-	if(zmax>0. && zmin<0.) gauss=element->NewGauss(xyz_list,xyz_list_front,3,30);//refined in vertical because of the sea level discontinuity
-	else                   gauss=element->NewGauss(xyz_list,xyz_list_front,3,3);
-
-	/* Start  looping on the number of gaussian points: */
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
-
-		element->JacobianDeterminantSurface(&Jdet,xyz_list_front,gauss);
-		element->NodalFunctionsVelocity(vbasis,gauss);
-		surface_input->GetInputValue(&surface,gauss);
-		if(dim==3) z=element->GetZcoord(xyz_list,gauss);
-		else       z=element->GetYcoord(xyz_list,gauss);
-		pressure = rho_water*gravity*min(0.,z);//0 if the gaussian point is above water level
-
-		for (int i=0;i<vnumnodes;i++){
-			pe->values[dim*i+0]+= pressure*Jdet*gauss->weight*normal[0]*vbasis[i];
-			pe->values[dim*i+1]+= pressure*Jdet*gauss->weight*normal[1]*vbasis[i];
-			if(dim==3) pe->values[dim*i+2]+= pressure*Jdet*gauss->weight*normal[2]*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);
-	xDelete<IssmDouble>(xyz_list_front);
-	return pe;
-}/*}}}*/
-void StressbalanceAnalysis::GetBFS(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+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        ]
@@ -4318,7 +4318,67 @@
 	xDelete<IssmDouble>(vdbasis);
 	xDelete<IssmDouble>(pbasis);
 }/*}}}*/
-void StressbalanceAnalysis::GetBFSprime(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void           StressbalanceAnalysis::GetBFSFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+	/* Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof. 
+	 * For node i, Li can be expressed in the actual coordinate system
+	 * by in 3d 
+	 *       Li=[ h 0 0 0 ]
+	 *	 	      [ 0 h 0 0 ]
+	 *	in 2d:
+	 *       Li=[ h 0 0 ]
+	 * where h is the interpolation function for node i.
+	 */
+
+	/*Fetch number of nodes for this finite element*/
+	int pnumnodes = element->NumberofNodesPressure();
+	int vnumnodes = element->NumberofNodesVelocity();
+	int pnumdof   = pnumnodes;
+	int vnumdof   = vnumnodes*dim;
+
+	/*Get nodal functions derivatives*/
+	IssmDouble* vbasis=xNew<IssmDouble>(vnumnodes);
+	element->NodalFunctionsVelocity(vbasis,gauss);
+
+	/*Build B: */
+	if(dim==3){
+		for(int i=0;i<vnumnodes;i++){
+			B[(vnumdof+pnumdof)*0+3*i+0] = vbasis[i];
+			B[(vnumdof+pnumdof)*0+3*i+1] = 0.;
+			B[(vnumdof+pnumdof)*0+3*i+2] = 0.;
+
+			B[(vnumdof+pnumdof)*1+3*i+0] = 0.;
+			B[(vnumdof+pnumdof)*1+3*i+1] = vbasis[i];
+			B[(vnumdof+pnumdof)*1+3*i+2] = 0.;
+
+			B[(vnumdof+pnumdof)*2+3*i+0] = 0.;
+			B[(vnumdof+pnumdof)*2+3*i+1] = 0.;
+			B[(vnumdof+pnumdof)*2+3*i+2] = vbasis[i];
+		}
+		for(int i=0;i<pnumnodes;i++){
+			B[(vnumdof+pnumdof)*0+i+vnumdof+0] = 0.;
+			B[(vnumdof+pnumdof)*1+i+vnumdof+0] = 0.;
+			B[(vnumdof+pnumdof)*2+i+vnumdof+0] = 0.;
+		}
+	}
+	else{
+		for(int i=0;i<vnumnodes;i++){
+			B[(vnumdof+pnumdof)*0+2*i+0] = vbasis[i];
+			B[(vnumdof+pnumdof)*0+2*i+1] = 0.;
+
+			B[(vnumdof+pnumdof)*1+2*i+0] = 0.;
+			B[(vnumdof+pnumdof)*1+2*i+1] = vbasis[i];
+		}
+
+		for(int i=0;i<pnumnodes;i++){
+			B[(vnumdof+pnumdof)*0+i+vnumdof+0] = 0.;
+			B[(vnumdof+pnumdof)*1+i+vnumdof+0] = 0.;
+		}
+	}
+
+	/*Clean-up*/
+	xDelete<IssmDouble>(vbasis);
+}/*}}}*/
+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: 
@@ -4431,24 +4491,15 @@
 	xDelete<IssmDouble>(vdbasis);
 	xDelete<IssmDouble>(pbasis);
 }/*}}}*/
-void StressbalanceAnalysis::GetBFSvel(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]
+void           StressbalanceAnalysis::GetBFSprimeUzawa(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+	/*	Compute B'  matrix. B'=[B1' B2' B3' B4' B5' B6'] where Bi' is of size 3*NDOF2. 
+	 *	For node i, Bi' can be expressed in the actual coordinate system
+	 *	by: 
+	 *			Bvi' = [  dphi/dx   dphi/dy ]
 	 *
-	 *
-	 *	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  ]
-	 *
+	 *	In 3d
+	 *     	   Bvi=[ dh/dx   dh/dy    dh/dz  ]
 	 *	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*/
@@ -4458,44 +4509,25 @@
 	IssmDouble* vdbasis=xNew<IssmDouble>(dim*vnumnodes);
 	element->NodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);
 
-	/*Build B: */
+	/*Build B_prime: */
 	if(dim==2){
 		for(int i=0;i<vnumnodes;i++){
-			B[(dim*vnumnodes)*0+dim*i+0] = vdbasis[0*vnumnodes+i];
-			B[(dim*vnumnodes)*0+dim*i+1] = 0.;
-			B[(dim*vnumnodes)*1+dim*i+0] = 0.;
-			B[(dim*vnumnodes)*1+dim*i+1] = vdbasis[1*vnumnodes+i];
-			B[(dim*vnumnodes)*2+dim*i+0] = .5*vdbasis[1*vnumnodes+i];
-			B[(dim*vnumnodes)*2+dim*i+1] = .5*vdbasis[0*vnumnodes+i];
+			Bprime[dim*i+0] = vdbasis[0*vnumnodes+i];
+			Bprime[dim*i+1] = vdbasis[1*vnumnodes+i];
 		}
 	}
 	else{
 		for(int i=0;i<vnumnodes;i++){
-			B[(dim*vnumnodes)*0+dim*i+0] = vdbasis[0*vnumnodes+i];
-			B[(dim*vnumnodes)*0+dim*i+1] = 0.;
-			B[(dim*vnumnodes)*0+dim*i+2] = 0.;
-			B[(dim*vnumnodes)*1+dim*i+0] = 0.;
-			B[(dim*vnumnodes)*1+dim*i+1] = vdbasis[1*vnumnodes+i];
-			B[(dim*vnumnodes)*1+dim*i+2] = 0.;
-			B[(dim*vnumnodes)*2+dim*i+0] = 0.;
-			B[(dim*vnumnodes)*2+dim*i+1] = 0.;
-			B[(dim*vnumnodes)*2+dim*i+2] = vdbasis[2*vnumnodes+i];
-			B[(dim*vnumnodes)*3+dim*i+0] = .5*vdbasis[1*vnumnodes+i];
-			B[(dim*vnumnodes)*3+dim*i+1] = .5*vdbasis[0*vnumnodes+i];
-			B[(dim*vnumnodes)*3+dim*i+2] = 0.;
-			B[(dim*vnumnodes)*4+dim*i+0] = .5*vdbasis[2*vnumnodes+i];
-			B[(dim*vnumnodes)*4+dim*i+1] = 0.;
-			B[(dim*vnumnodes)*4+dim*i+2] = .5*vdbasis[0*vnumnodes+i];
-			B[(dim*vnumnodes)*5+dim*i+0] = 0.;
-			B[(dim*vnumnodes)*5+dim*i+1] = .5*vdbasis[2*vnumnodes+i];
-			B[(dim*vnumnodes)*5+dim*i+2] = .5*vdbasis[1*vnumnodes+i];
+			Bprime[dim*i+0] = vdbasis[0*vnumnodes+i];
+			Bprime[dim*i+1] = vdbasis[1*vnumnodes+i];
+			Bprime[dim*i+2] = vdbasis[2*vnumnodes+i];
 		}
 	}
 
 	/*Clean up*/
 	xDelete<IssmDouble>(vdbasis);
 }/*}}}*/
-void StressbalanceAnalysis::GetBFSprimevel(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void           StressbalanceAnalysis::GetBFSprimevel(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: 
@@ -4558,7 +4590,7 @@
 	/*Clean up*/
 	xDelete<IssmDouble>(vdbasis);
 }/*}}}*/
-void StressbalanceAnalysis::GetBFSUzawa(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void           StressbalanceAnalysis::GetBFSUzawa(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 	/*Compute B  matrix. B=[Bp1 Bp2 ...] where Bpi=phi_pi. 
 	 */
 
@@ -4580,15 +4612,24 @@
 	/*Clean up*/
 	xDelete<IssmDouble>(basis);
 }/*}}}*/
-void StressbalanceAnalysis::GetBFSprimeUzawa(IssmDouble* Bprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*	Compute B'  matrix. B'=[B1' B2' B3' B4' B5' B6'] where Bi' is of size 3*NDOF2. 
-	 *	For node i, Bi' can be expressed in the actual coordinate system
-	 *	by: 
-	 *			Bvi' = [  dphi/dx   dphi/dy ]
+void           StressbalanceAnalysis::GetBFSvel(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]
 	 *
-	 *	In 3d
-	 *     	   Bvi=[ dh/dx   dh/dy    dh/dz  ]
+	 *
+	 *	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  ]
+	 *
 	 *	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*/
@@ -4598,85 +4639,44 @@
 	IssmDouble* vdbasis=xNew<IssmDouble>(dim*vnumnodes);
 	element->NodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);
 
-	/*Build B_prime: */
+	/*Build B: */
 	if(dim==2){
 		for(int i=0;i<vnumnodes;i++){
-			Bprime[dim*i+0] = vdbasis[0*vnumnodes+i];
-			Bprime[dim*i+1] = vdbasis[1*vnumnodes+i];
+			B[(dim*vnumnodes)*0+dim*i+0] = vdbasis[0*vnumnodes+i];
+			B[(dim*vnumnodes)*0+dim*i+1] = 0.;
+			B[(dim*vnumnodes)*1+dim*i+0] = 0.;
+			B[(dim*vnumnodes)*1+dim*i+1] = vdbasis[1*vnumnodes+i];
+			B[(dim*vnumnodes)*2+dim*i+0] = .5*vdbasis[1*vnumnodes+i];
+			B[(dim*vnumnodes)*2+dim*i+1] = .5*vdbasis[0*vnumnodes+i];
 		}
 	}
 	else{
 		for(int i=0;i<vnumnodes;i++){
-			Bprime[dim*i+0] = vdbasis[0*vnumnodes+i];
-			Bprime[dim*i+1] = vdbasis[1*vnumnodes+i];
-			Bprime[dim*i+2] = vdbasis[2*vnumnodes+i];
+			B[(dim*vnumnodes)*0+dim*i+0] = vdbasis[0*vnumnodes+i];
+			B[(dim*vnumnodes)*0+dim*i+1] = 0.;
+			B[(dim*vnumnodes)*0+dim*i+2] = 0.;
+			B[(dim*vnumnodes)*1+dim*i+0] = 0.;
+			B[(dim*vnumnodes)*1+dim*i+1] = vdbasis[1*vnumnodes+i];
+			B[(dim*vnumnodes)*1+dim*i+2] = 0.;
+			B[(dim*vnumnodes)*2+dim*i+0] = 0.;
+			B[(dim*vnumnodes)*2+dim*i+1] = 0.;
+			B[(dim*vnumnodes)*2+dim*i+2] = vdbasis[2*vnumnodes+i];
+			B[(dim*vnumnodes)*3+dim*i+0] = .5*vdbasis[1*vnumnodes+i];
+			B[(dim*vnumnodes)*3+dim*i+1] = .5*vdbasis[0*vnumnodes+i];
+			B[(dim*vnumnodes)*3+dim*i+2] = 0.;
+			B[(dim*vnumnodes)*4+dim*i+0] = .5*vdbasis[2*vnumnodes+i];
+			B[(dim*vnumnodes)*4+dim*i+1] = 0.;
+			B[(dim*vnumnodes)*4+dim*i+2] = .5*vdbasis[0*vnumnodes+i];
+			B[(dim*vnumnodes)*5+dim*i+0] = 0.;
+			B[(dim*vnumnodes)*5+dim*i+1] = .5*vdbasis[2*vnumnodes+i];
+			B[(dim*vnumnodes)*5+dim*i+2] = .5*vdbasis[1*vnumnodes+i];
 		}
 	}
 
 	/*Clean up*/
 	xDelete<IssmDouble>(vdbasis);
 }/*}}}*/
-void StressbalanceAnalysis::GetBFSFriction(IssmDouble* B,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/* Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof. 
-	 * For node i, Li can be expressed in the actual coordinate system
-	 * by in 3d 
-	 *       Li=[ h 0 0 0 ]
-	 *	 	      [ 0 h 0 0 ]
-	 *	in 2d:
-	 *       Li=[ h 0 0 ]
-	 * where h is the interpolation function for node i.
-	 */
-
-	/*Fetch number of nodes for this finite element*/
-	int pnumnodes = element->NumberofNodesPressure();
-	int vnumnodes = element->NumberofNodesVelocity();
-	int pnumdof   = pnumnodes;
-	int vnumdof   = vnumnodes*dim;
-
-	/*Get nodal functions derivatives*/
-	IssmDouble* vbasis=xNew<IssmDouble>(vnumnodes);
-	element->NodalFunctionsVelocity(vbasis,gauss);
-
-	/*Build B: */
-	if(dim==3){
-		for(int i=0;i<vnumnodes;i++){
-			B[(vnumdof+pnumdof)*0+3*i+0] = vbasis[i];
-			B[(vnumdof+pnumdof)*0+3*i+1] = 0.;
-			B[(vnumdof+pnumdof)*0+3*i+2] = 0.;
-
-			B[(vnumdof+pnumdof)*1+3*i+0] = 0.;
-			B[(vnumdof+pnumdof)*1+3*i+1] = vbasis[i];
-			B[(vnumdof+pnumdof)*1+3*i+2] = 0.;
-
-			B[(vnumdof+pnumdof)*2+3*i+0] = 0.;
-			B[(vnumdof+pnumdof)*2+3*i+1] = 0.;
-			B[(vnumdof+pnumdof)*2+3*i+2] = vbasis[i];
-		}
-		for(int i=0;i<pnumnodes;i++){
-			B[(vnumdof+pnumdof)*0+i+vnumdof+0] = 0.;
-			B[(vnumdof+pnumdof)*1+i+vnumdof+0] = 0.;
-			B[(vnumdof+pnumdof)*2+i+vnumdof+0] = 0.;
-		}
-	}
-	else{
-		for(int i=0;i<vnumnodes;i++){
-			B[(vnumdof+pnumdof)*0+2*i+0] = vbasis[i];
-			B[(vnumdof+pnumdof)*0+2*i+1] = 0.;
-
-			B[(vnumdof+pnumdof)*1+2*i+0] = 0.;
-			B[(vnumdof+pnumdof)*1+2*i+1] = vbasis[i];
-		}
-
-		for(int i=0;i<pnumnodes;i++){
-			B[(vnumdof+pnumdof)*0+i+vnumdof+0] = 0.;
-			B[(vnumdof+pnumdof)*1+i+vnumdof+0] = 0.;
-		}
-	}
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(vbasis);
-}/*}}}*/
-void StressbalanceAnalysis::GetCFS(IssmDouble* C,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void           StressbalanceAnalysis::GetCFS(IssmDouble* C,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 	/*Compute C  matrix. C=[Cp1 Cp2 ...] where:
 	 *     Cpi=[phi phi]. 
 	 */
@@ -4697,7 +4697,7 @@
 	/*Clean up*/
 	xDelete<IssmDouble>(basis);
 }/*}}}*/
-void StressbalanceAnalysis::GetCFSprime(IssmDouble* Cprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void           StressbalanceAnalysis::GetCFSprime(IssmDouble* Cprime,Element* element,int dim,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 	/*	Compute C'  matrix. C'=[C1' C2' ...] 
 	 *			Ci' = [  phi  0  ]
 	 *			      [   0  phi ]
@@ -4745,7 +4745,7 @@
 	/*Clean-up*/
 	xDelete<IssmDouble>(vbasis);
 }/*}}}*/
-void StressbalanceAnalysis::GetSolutionFromInputsFS(Vector<IssmDouble>* solution,Element* element){/*{{{*/
+void           StressbalanceAnalysis::GetSolutionFromInputsFS(Vector<IssmDouble>* solution,Element* element){/*{{{*/
 
 	int*         vdoflist=NULL;
 	int*         pdoflist=NULL;
@@ -4808,7 +4808,7 @@
 	xDelete<IssmDouble>(pvalues);
 	xDelete<IssmDouble>(vvalues);
 }/*}}}*/
-void StressbalanceAnalysis::InitializeXTH(Elements* elements,Parameters* parameters){/*{{{*/
+void           StressbalanceAnalysis::InitializeXTH(Elements* elements,Parameters* parameters){/*{{{*/
 
 	/*Intermediaries*/
 	int        dim;
@@ -4887,7 +4887,7 @@
 	}
 
 }/*}}}*/
-void StressbalanceAnalysis::InputUpdateFromSolutionFS(IssmDouble* solution,Element* element){/*{{{*/
+void           StressbalanceAnalysis::InputUpdateFromSolutionFS(IssmDouble* solution,Element* element){/*{{{*/
 
 	bool         results_on_nodes;
 	int          i,dim;
@@ -4977,7 +4977,7 @@
 	xDelete<int>(pdoflist);
 	xDelete<int>(cs_list);
 }/*}}}*/
-void StressbalanceAnalysis::InputUpdateFromSolutionFSXTH_d(Elements* elements,Parameters* parameters){/*{{{*/
+void           StressbalanceAnalysis::InputUpdateFromSolutionFSXTH_d(Elements* elements,Parameters* parameters){/*{{{*/
 
 	/*Intermediaries*/
 	int         dim,tausize;
@@ -5190,7 +5190,7 @@
 		xDelete<IssmDouble>(pe_yz); xDelete<IssmDouble>(d_yz);
 	}
 }/*}}}*/
-void StressbalanceAnalysis::InputUpdateFromSolutionFSXTH_tau(Elements* elements,Parameters* parameters){/*{{{*/
+void           StressbalanceAnalysis::InputUpdateFromSolutionFSXTH_tau(Elements* elements,Parameters* parameters){/*{{{*/
 
 	/*Intermediaries*/
 	int         dim,tausize;
@@ -5326,165 +5326,6 @@
 }/*}}}*/
 
 /*Coupling (Tiling)*/
-ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSA3d(Element* element){/*{{{*/
-
-	/*compute all stiffness matrices for this element*/
-	ElementMatrix* Ke1=CreateKMatrixSSA3dViscous(element);
-	ElementMatrix* Ke2=CreateKMatrixSSA3dFriction(element);
-	ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
-
-	/*clean-up and return*/
-	delete Ke1;
-	delete Ke2;
-	return Ke;
-}/*}}}*/
-ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSA3dFriction(Element* element){/*{{{*/
-
-	/*Initialize Element matrix and return if necessary*/
-	if(element->IsFloating() || !element->IsOnBase()) return NULL;
-
-	/*Build a tria element using the 3 nodes of the base of the penta. Then use 
-	 * the tria functionality to build a friction stiffness matrix on these 3
-	 * nodes: */
-	Element* basalelement = element->SpawnBasalElement();
-	ElementMatrix* Ke=CreateKMatrixSSAFriction(basalelement);
-	basalelement->DeleteMaterials(); delete basalelement;
-
-	/*clean-up and return*/
-	return Ke;
-}/*}}}*/
-ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSA3dViscous(Element* element){/*{{{*/
-
-	/*Constants*/
-	const int    numdof2d=2*3;
-
-	/*Intermediaries */
-	int         i,j,approximation;
-	int         dim=3;
-	IssmDouble  Jdet,viscosity,oldviscosity,newviscosity,viscosity_overshoot;
-	IssmDouble  epsilon[5],oldepsilon[5];       /* epsilon=[exx,eyy,exy,exz,eyz];*/
-	IssmDouble  epsilons[6];                    //6 for FS
-	IssmDouble  B[3][numdof2d];
-	IssmDouble  Bprime[3][numdof2d];
-	IssmDouble  D[3][3]= {0.0};                 // material matrix, simple scalar matrix.
-	IssmDouble  D_scalar;
-	IssmDouble  Ke_gg[numdof2d][numdof2d]={0.0};
-	IssmDouble  *xyz_list  = NULL;
-
-	/*Find penta on bed as this is a SSA elements: */
-	Element* pentabase=element->GetBasalElement();
-	Element* basaltria=pentabase->SpawnBasalElement();
-
-	/*Initialize Element matrix*/
-	ElementMatrix* Ke=basaltria->NewElementMatrix(SSAApproximationEnum);
-	element->GetInputValue(&approximation,ApproximationEnum);
-
-	/*Retrieve all inputs and parameters*/
-	element->GetVerticesCoordinates(&xyz_list);
-	element->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum);
-	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);
-	Input* vz_input   =element->GetInput(VzEnum);       _assert_(vz_input);
-
-	/* Start  looping on the number of gaussian points: */
-	Gauss* gauss=element->NewGauss(5);
-	Gauss* gauss_tria=new GaussTria();
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
-		gauss->SynchronizeGaussBase(gauss_tria);
-
-		element->JacobianDeterminant(&Jdet, xyz_list,gauss);
-		this->GetBSSA(&B[0][0],basaltria,2,xyz_list, gauss_tria);
-		this->GetBSSAprime(&Bprime[0][0],basaltria,2,xyz_list, gauss_tria);
-
-		if(approximation==SSAHOApproximationEnum){
-			element->ViscosityHO(&viscosity,dim,xyz_list,gauss,vx_input,vy_input);
-			element->ViscosityHO(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input);
-			newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
-		}
-		else if (approximation==SSAFSApproximationEnum){
-			element->ViscosityFS(&newviscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
-		}
-		else _error_("approximation " << approximation << " (" << EnumToStringx(approximation) << ") not supported yet");
-
-		D_scalar=2*newviscosity*gauss->weight*Jdet;
-		for (i=0;i<3;i++) D[i][i]=D_scalar;
-
-		TripleMultiply( &B[0][0],3,numdof2d,1,
-					&D[0][0],3,3,0,
-					&Bprime[0][0],3,numdof2d,0,
-					&Ke_gg[0][0],1);
-
-	}
-	for(i=0;i<numdof2d;i++) for(j=0;j<numdof2d;j++) Ke->values[i*numdof2d+j]+=Ke_gg[i][j];
-
-	/*Transform Coordinate System*/
-	basaltria->TransformStiffnessMatrixCoord(Ke,XYEnum);
-
-	/*Clean up and return*/
-	xDelete<IssmDouble>(xyz_list);
-	delete basaltria->material;
-	delete basaltria;
-	delete gauss_tria;
-	delete gauss;
-	return Ke;
-
-}/*}}}*/
-ElementMatrix* StressbalanceAnalysis::CreateKMatrixHOFS(Element* element){/*{{{*/
-
-	/*compute all stiffness matrices for this element*/
-	ElementMatrix* Ke1=CreateKMatrixFS(element);
-	int indices[3]={18,19,20};
-	Ke1->StaticCondensation(3,&indices[0]);
-	int init = element->FiniteElement();
-	element->SetTemporaryElementType(P1Enum); // P1 needed for HO
-	ElementMatrix* Ke2=CreateKMatrixHO(element);
-	element->SetTemporaryElementType(init); // P1 needed for HO
-	ElementMatrix* Ke3=CreateKMatrixCouplingHOFS(element);
-	ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
-
-	/*clean-up and return*/
-	delete Ke1;
-	delete Ke2;
-	delete Ke3;
-	return Ke;
-}/*}}}*/
-ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAHO(Element* element){/*{{{*/
-
-	/*compute all stiffness matrices for this element*/
-	ElementMatrix* Ke1=CreateKMatrixSSA3d(element);
-	ElementMatrix* Ke2=CreateKMatrixHO(element);
-	ElementMatrix* Ke3=CreateKMatrixCouplingSSAHO(element);
-	ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
-
-	/*clean-up and return*/
-	delete Ke1;
-	delete Ke2;
-	delete Ke3;
-	return Ke;
-}/*}}}*/
-ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAFS(Element* element){/*{{{*/
-
-	/*compute all stiffness matrices for this element*/
-	ElementMatrix* Ke1=CreateKMatrixFS(element);
-	int indices[3]={18,19,20};
-	Ke1->StaticCondensation(3,&indices[0]);
-	int init = element->FiniteElement();
-	element->SetTemporaryElementType(P1Enum); 
-	ElementMatrix* Ke2=CreateKMatrixSSA3d(element);
-	element->SetTemporaryElementType(init); 
-	ElementMatrix* Ke3=CreateKMatrixCouplingSSAFS(element);
-	ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
-
-	/*clean-up and return*/
-	delete Ke1;
-	delete Ke2;
-	delete Ke3;
-	return Ke;
-}/*}}}*/
 ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingHOFS(Element* element){/*{{{*/
 
 	/*Constants*/
@@ -5551,184 +5392,6 @@
 	delete Ke2;
 	return Ke;
 }/*}}}*/
-ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAHO(Element* element){/*{{{*/
-
-	/*compute all stiffness matrices for this element*/
-	ElementMatrix* Ke1=CreateKMatrixCouplingSSAHOViscous(element);
-	ElementMatrix* Ke2=CreateKMatrixCouplingSSAHOFriction(element);
-	ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
-
-	/*clean-up and return*/
-	delete Ke1;
-	delete Ke2;
-	return Ke;
-}/*}}}*/
-ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAHOFriction(Element* element){/*{{{*/
-
-	if(element->IsFloating() || !element->IsOnBase()) return NULL;
-
-	/*Constants*/
-	int numnodes    = element->GetNumberOfNodes();
-	int numdof      = 2*numnodes;
-	int numdoftotal = 4*numnodes;
-
-	/*Intermediaries */
-	int         i,j;
-	IssmDouble  Jdet2d,alpha2;
-	IssmDouble *xyz_list_tria = NULL;
-	IssmDouble* L             = xNewZeroInit<IssmDouble>(2*numdof);
-	IssmDouble  DL[2][2]      = {{ 0,0 },{0,0}}; //for basal drag
-	IssmDouble  DL_scalar;
-	IssmDouble* Ke_gg         = xNewZeroInit<IssmDouble>(numdof*numdof);
-	Node      **node_list     = xNew<Node*>(2*numnodes);
-	int*        cs_list       = xNew<int>(2*numnodes);
-
-	/*Initialize Element matrix and return if necessary*/
-	ElementMatrix* Ke1=element->NewElementMatrix(SSAApproximationEnum);
-	ElementMatrix* Ke2=element->NewElementMatrix(HOApproximationEnum);
-	ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
-	delete Ke1; delete Ke2;
-
-	/*Prepare node list*/
-	for(i=0;i<numnodes;i++){
-		node_list[i+0*numnodes] = element->GetNode(i);
-		node_list[i+1*numnodes] = element->GetNode(i);
-		cs_list[i+0*numnodes] = XYEnum;
-		cs_list[i+1*numnodes] = XYEnum;
-	}
-
-	/*retrieve inputs :*/
-	element->GetVerticesCoordinatesBase(&xyz_list_tria);
-
-	/*build friction object, used later on: */
-	Friction* friction=new Friction(element,2);
-
-	/* Start  looping on the number of gaussian points: */
-	Gauss* gauss=element->NewGaussBase(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
-
-		/*Friction: */
-		friction->GetAlpha2(&alpha2,gauss);
-		element->JacobianDeterminantBase(&Jdet2d, xyz_list_tria,gauss);
-		this->GetBHOFriction(L,element,3,xyz_list_tria,gauss);
-
-		DL_scalar=alpha2*gauss->weight*Jdet2d;
-		for (i=0;i<2;i++) DL[i][i]=DL_scalar; 
-
-		/*  Do the triple producte tL*D*L: */
-		TripleMultiply( L,2,numdof,1,
-					&DL[0][0],2,2,0,
-					L,2,numdof,0,
-					Ke_gg,1);
-	}
-
-	for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdoftotal+(numdof+j)]+=Ke_gg[i*numdof+j];
-	for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[(i+numdof)*numdoftotal+j]+=Ke_gg[i*numdof+j];
-
-	/*Transform Coordinate System*/
-	element->TransformStiffnessMatrixCoord(Ke,node_list,2*numnodes,cs_list);
-
-	/*Clean up and return*/
-	delete gauss;
-	delete friction;
-	xDelete<int>(cs_list);
-	xDelete<Node*>(node_list);
-	xDelete<IssmDouble>(xyz_list_tria);
-	xDelete<IssmDouble>(Ke_gg);
-	xDelete<IssmDouble>(L);
-	return Ke;
-}/*}}}*/
-ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAHOViscous(Element* element){/*{{{*/
-
-	/*Constants*/
-	int numnodes    = element->GetNumberOfNodes();
-	int numdofm     = 1 *numnodes; //*2/2
-	int numdofp     = 2 *numnodes;
-	int numdoftotal = 2 *2 *numnodes;//2 dof per nodes and 2 sets of nodes for HO and SSA
-
-	/*Intermediaries */
-	int         i,j;
-	IssmDouble  Jdet,viscosity,oldviscosity,newviscosity,viscosity_overshoot; //viscosity
-	IssmDouble  *xyz_list      = NULL;
-	IssmDouble* B              = xNew<IssmDouble>(3*numdofp);
-	IssmDouble* Bprime         = xNew<IssmDouble>(3*numdofm);
-	IssmDouble  D[3][3]={0.0}; // material matrix, simple scalar matrix.
-	IssmDouble  D_scalar;
-	IssmDouble* Ke_gg          = xNewZeroInit<IssmDouble>(numdofp*numdofm);
-	Node       **node_list     = xNew<Node*>(2*numnodes);
-	int*         cs_list= xNew<int>(2*numnodes);
-
-	/*Find penta on bed as HO must be coupled to the dofs on the bed: */
-	Element* pentabase=element->GetBasalElement();
-	Element* basaltria=pentabase->SpawnBasalElement();
-
-	/*prepare node list*/
-	for(i=0;i<numnodes;i++){
-		node_list[i+0*numnodes] = pentabase->GetNode(i);
-		node_list[i+1*numnodes] = element  ->GetNode(i);
-		cs_list[i+0*numnodes] = XYEnum;
-		cs_list[i+1*numnodes] = XYEnum;
-	}
-
-	/*Initialize Element matrix*/
-	ElementMatrix* Ke1= pentabase->NewElementMatrix(SSAApproximationEnum);
-	ElementMatrix* Ke2= element  ->NewElementMatrix(HOApproximationEnum);
-	ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
-	delete Ke1; delete Ke2;
-
-	/* Get node coordinates and dof list: */
-	element->GetVerticesCoordinates(&xyz_list);
-	element->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum);
-	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);
-
-	/* Start  looping on the number of gaussian points: */
-	Gauss* gauss=element->NewGauss(5);
-	Gauss* gauss_tria=new GaussTria();
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
-		gauss->SynchronizeGaussBase(gauss_tria);
-
-		element->JacobianDeterminant(&Jdet, xyz_list,gauss);
-		this->GetBSSAHO(B, element,xyz_list, gauss);
-		this->GetBSSAprime(Bprime,basaltria,2,xyz_list, gauss_tria); 
-		element->ViscosityHO(&viscosity,3,xyz_list,gauss,vx_input,vy_input);
-		element->ViscosityHO(&oldviscosity,3,xyz_list,gauss,vxold_input,vyold_input);
-
-		newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
-		D_scalar=2*newviscosity*gauss->weight*Jdet;
-		for (i=0;i<3;i++) D[i][i]=D_scalar;
-
-		TripleMultiply( B,3,numdofp,1,
-					&D[0][0],3,3,0,
-					Bprime,3,numdofm,0,
-					Ke_gg,1);
-	} 
-	for(i=0;i<numdofp;i++) for(j=0;j<numdofm;j++) Ke->values[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i*numdofm+j];
-	for(i=0;i<numdofm;i++) for(j=0;j<numdofp;j++) Ke->values[i*numdoftotal+(j+2*numdofm)]+=Ke_gg[j*numdofm+i];
-
-	/*Transform Coordinate System*/
-	element->TransformStiffnessMatrixCoord(Ke,node_list,2*numnodes,cs_list);
-
-	/*Clean-up and return*/
-	basaltria->DeleteMaterials(); delete basaltria;
-	
-	delete gauss;
-	delete gauss_tria;
-	xDelete<IssmDouble>(B);
-	xDelete<IssmDouble>(Bprime);
-	xDelete<IssmDouble>(Ke_gg);
-	xDelete<IssmDouble>(xyz_list);
-	xDelete<Node*>(node_list);
-	xDelete<int>(cs_list);
-	return Ke;
-
-}/*}}}*/
 ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAFS(Element* element){/*{{{*/
 
 	/*compute all stiffness matrices for this element*/
@@ -5969,64 +5632,343 @@
 	return Ke;
 
 }/*}}}*/
-ElementVector* StressbalanceAnalysis::CreatePVectorSSAHO(Element* element){/*{{{*/
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAHO(Element* element){/*{{{*/
 
-	/*compute all load vectors for this element*/
-	ElementVector* pe1=CreatePVectorSSA(element);
-	ElementVector* pe2=CreatePVectorHO(element);
-	ElementVector* pe =new ElementVector(pe1,pe2);
+	/*compute all stiffness matrices for this element*/
+	ElementMatrix* Ke1=CreateKMatrixCouplingSSAHOViscous(element);
+	ElementMatrix* Ke2=CreateKMatrixCouplingSSAHOFriction(element);
+	ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
 
 	/*clean-up and return*/
-	delete pe1;
-	delete pe2;
-	return pe;
+	delete Ke1;
+	delete Ke2;
+	return Ke;
 }/*}}}*/
-ElementVector* StressbalanceAnalysis::CreatePVectorSSAFS(Element* element){/*{{{*/
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAHOFriction(Element* element){/*{{{*/
 
-	/*compute all load vectors for this element*/
+	if(element->IsFloating() || !element->IsOnBase()) return NULL;
+
+	/*Constants*/
+	int numnodes    = element->GetNumberOfNodes();
+	int numdof      = 2*numnodes;
+	int numdoftotal = 4*numnodes;
+
+	/*Intermediaries */
+	int         i,j;
+	IssmDouble  Jdet2d,alpha2;
+	IssmDouble *xyz_list_tria = NULL;
+	IssmDouble* L             = xNewZeroInit<IssmDouble>(2*numdof);
+	IssmDouble  DL[2][2]      = {{ 0,0 },{0,0}}; //for basal drag
+	IssmDouble  DL_scalar;
+	IssmDouble* Ke_gg         = xNewZeroInit<IssmDouble>(numdof*numdof);
+	Node      **node_list     = xNew<Node*>(2*numnodes);
+	int*        cs_list       = xNew<int>(2*numnodes);
+
+	/*Initialize Element matrix and return if necessary*/
+	ElementMatrix* Ke1=element->NewElementMatrix(SSAApproximationEnum);
+	ElementMatrix* Ke2=element->NewElementMatrix(HOApproximationEnum);
+	ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
+	delete Ke1; delete Ke2;
+
+	/*Prepare node list*/
+	for(i=0;i<numnodes;i++){
+		node_list[i+0*numnodes] = element->GetNode(i);
+		node_list[i+1*numnodes] = element->GetNode(i);
+		cs_list[i+0*numnodes] = XYEnum;
+		cs_list[i+1*numnodes] = XYEnum;
+	}
+
+	/*retrieve inputs :*/
+	element->GetVerticesCoordinatesBase(&xyz_list_tria);
+
+	/*build friction object, used later on: */
+	Friction* friction=new Friction(element,2);
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=element->NewGaussBase(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+
+		/*Friction: */
+		friction->GetAlpha2(&alpha2,gauss);
+		element->JacobianDeterminantBase(&Jdet2d, xyz_list_tria,gauss);
+		this->GetBHOFriction(L,element,3,xyz_list_tria,gauss);
+
+		DL_scalar=alpha2*gauss->weight*Jdet2d;
+		for (i=0;i<2;i++) DL[i][i]=DL_scalar; 
+
+		/*  Do the triple producte tL*D*L: */
+		TripleMultiply( L,2,numdof,1,
+					&DL[0][0],2,2,0,
+					L,2,numdof,0,
+					Ke_gg,1);
+	}
+
+	for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdoftotal+(numdof+j)]+=Ke_gg[i*numdof+j];
+	for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[(i+numdof)*numdoftotal+j]+=Ke_gg[i*numdof+j];
+
+	/*Transform Coordinate System*/
+	element->TransformStiffnessMatrixCoord(Ke,node_list,2*numnodes,cs_list);
+
+	/*Clean up and return*/
+	delete gauss;
+	delete friction;
+	xDelete<int>(cs_list);
+	xDelete<Node*>(node_list);
+	xDelete<IssmDouble>(xyz_list_tria);
+	xDelete<IssmDouble>(Ke_gg);
+	xDelete<IssmDouble>(L);
+	return Ke;
+}/*}}}*/
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixCouplingSSAHOViscous(Element* element){/*{{{*/
+
+	/*Constants*/
+	int numnodes    = element->GetNumberOfNodes();
+	int numdofm     = 1 *numnodes; //*2/2
+	int numdofp     = 2 *numnodes;
+	int numdoftotal = 2 *2 *numnodes;//2 dof per nodes and 2 sets of nodes for HO and SSA
+
+	/*Intermediaries */
+	int         i,j;
+	IssmDouble  Jdet,viscosity,oldviscosity,newviscosity,viscosity_overshoot; //viscosity
+	IssmDouble  *xyz_list      = NULL;
+	IssmDouble* B              = xNew<IssmDouble>(3*numdofp);
+	IssmDouble* Bprime         = xNew<IssmDouble>(3*numdofm);
+	IssmDouble  D[3][3]={0.0}; // material matrix, simple scalar matrix.
+	IssmDouble  D_scalar;
+	IssmDouble* Ke_gg          = xNewZeroInit<IssmDouble>(numdofp*numdofm);
+	Node       **node_list     = xNew<Node*>(2*numnodes);
+	int*         cs_list= xNew<int>(2*numnodes);
+
+	/*Find penta on bed as HO must be coupled to the dofs on the bed: */
+	Element* pentabase=element->GetBasalElement();
+	Element* basaltria=pentabase->SpawnBasalElement();
+
+	/*prepare node list*/
+	for(i=0;i<numnodes;i++){
+		node_list[i+0*numnodes] = pentabase->GetNode(i);
+		node_list[i+1*numnodes] = element  ->GetNode(i);
+		cs_list[i+0*numnodes] = XYEnum;
+		cs_list[i+1*numnodes] = XYEnum;
+	}
+
+	/*Initialize Element matrix*/
+	ElementMatrix* Ke1= pentabase->NewElementMatrix(SSAApproximationEnum);
+	ElementMatrix* Ke2= element  ->NewElementMatrix(HOApproximationEnum);
+	ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
+	delete Ke1; delete Ke2;
+
+	/* Get node coordinates and dof list: */
+	element->GetVerticesCoordinates(&xyz_list);
+	element->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum);
+	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);
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=element->NewGauss(5);
+	Gauss* gauss_tria=new GaussTria();
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+		gauss->SynchronizeGaussBase(gauss_tria);
+
+		element->JacobianDeterminant(&Jdet, xyz_list,gauss);
+		this->GetBSSAHO(B, element,xyz_list, gauss);
+		this->GetBSSAprime(Bprime,basaltria,2,xyz_list, gauss_tria); 
+		element->ViscosityHO(&viscosity,3,xyz_list,gauss,vx_input,vy_input);
+		element->ViscosityHO(&oldviscosity,3,xyz_list,gauss,vxold_input,vyold_input);
+
+		newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
+		D_scalar=2*newviscosity*gauss->weight*Jdet;
+		for (i=0;i<3;i++) D[i][i]=D_scalar;
+
+		TripleMultiply( B,3,numdofp,1,
+					&D[0][0],3,3,0,
+					Bprime,3,numdofm,0,
+					Ke_gg,1);
+	} 
+	for(i=0;i<numdofp;i++) for(j=0;j<numdofm;j++) Ke->values[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i*numdofm+j];
+	for(i=0;i<numdofm;i++) for(j=0;j<numdofp;j++) Ke->values[i*numdoftotal+(j+2*numdofm)]+=Ke_gg[j*numdofm+i];
+
+	/*Transform Coordinate System*/
+	element->TransformStiffnessMatrixCoord(Ke,node_list,2*numnodes,cs_list);
+
+	/*Clean-up and return*/
+	basaltria->DeleteMaterials(); delete basaltria;
+	
+	delete gauss;
+	delete gauss_tria;
+	xDelete<IssmDouble>(B);
+	xDelete<IssmDouble>(Bprime);
+	xDelete<IssmDouble>(Ke_gg);
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<Node*>(node_list);
+	xDelete<int>(cs_list);
+	return Ke;
+
+}/*}}}*/
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixHOFS(Element* element){/*{{{*/
+
+	/*compute all stiffness matrices for this element*/
+	ElementMatrix* Ke1=CreateKMatrixFS(element);
+	int indices[3]={18,19,20};
+	Ke1->StaticCondensation(3,&indices[0]);
 	int init = element->FiniteElement();
 	element->SetTemporaryElementType(P1Enum); // P1 needed for HO
-	ElementVector* pe1=CreatePVectorSSA(element);
+	ElementMatrix* Ke2=CreateKMatrixHO(element);
 	element->SetTemporaryElementType(init); // P1 needed for HO
-	ElementVector* pe2=CreatePVectorFS(element);
-	int indices[3]={18,19,20};
-	element->SetTemporaryElementType(MINIcondensedEnum); // P1 needed for HO
-	ElementMatrix* Ke = CreateKMatrixFS(element);
-	element->SetTemporaryElementType(init); // P1 needed for HO
-	pe2->StaticCondensation(Ke,3,&indices[0]);
-	delete Ke;
-	ElementVector* pe3=CreatePVectorCouplingSSAFS(element);
-	ElementVector* pe =new ElementVector(pe1,pe2,pe3);
+	ElementMatrix* Ke3=CreateKMatrixCouplingHOFS(element);
+	ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
 
 	/*clean-up and return*/
-	delete pe1;
-	delete pe2;
-	delete pe3;
-	return pe;
+	delete Ke1;
+	delete Ke2;
+	delete Ke3;
+	return Ke;
 }/*}}}*/
-ElementVector* StressbalanceAnalysis::CreatePVectorHOFS(Element* element){/*{{{*/
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAFS(Element* element){/*{{{*/
 
-	/*compute all load vectors for this element*/
+	/*compute all stiffness matrices for this element*/
+	ElementMatrix* Ke1=CreateKMatrixFS(element);
+	int indices[3]={18,19,20};
+	Ke1->StaticCondensation(3,&indices[0]);
 	int init = element->FiniteElement();
-	element->SetTemporaryElementType(P1Enum);
-	ElementVector* pe1=CreatePVectorHO(element);
-	element->SetTemporaryElementType(init);
-	ElementVector* pe2=CreatePVectorFS(element);
-	int indices[3]={18,19,20};
-	element->SetTemporaryElementType(MINIcondensedEnum);
-	ElementMatrix* Ke = CreateKMatrixFS(element);
-	element->SetTemporaryElementType(init);
-	pe2->StaticCondensation(Ke,3,&indices[0]);
-	delete Ke;
-	ElementVector* pe3=CreatePVectorCouplingHOFS(element);
-	ElementVector* pe =new ElementVector(pe1,pe2,pe3);
+	element->SetTemporaryElementType(P1Enum); 
+	ElementMatrix* Ke2=CreateKMatrixSSA3d(element);
+	element->SetTemporaryElementType(init); 
+	ElementMatrix* Ke3=CreateKMatrixCouplingSSAFS(element);
+	ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
 
 	/*clean-up and return*/
-	delete pe1;
-	delete pe2;
-	delete pe3;
-	return pe;
+	delete Ke1;
+	delete Ke2;
+	delete Ke3;
+	return Ke;
 }/*}}}*/
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSAHO(Element* element){/*{{{*/
+
+	/*compute all stiffness matrices for this element*/
+	ElementMatrix* Ke1=CreateKMatrixSSA3d(element);
+	ElementMatrix* Ke2=CreateKMatrixHO(element);
+	ElementMatrix* Ke3=CreateKMatrixCouplingSSAHO(element);
+	ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
+
+	/*clean-up and return*/
+	delete Ke1;
+	delete Ke2;
+	delete Ke3;
+	return Ke;
+}/*}}}*/
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSA3d(Element* element){/*{{{*/
+
+	/*compute all stiffness matrices for this element*/
+	ElementMatrix* Ke1=CreateKMatrixSSA3dViscous(element);
+	ElementMatrix* Ke2=CreateKMatrixSSA3dFriction(element);
+	ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
+
+	/*clean-up and return*/
+	delete Ke1;
+	delete Ke2;
+	return Ke;
+}/*}}}*/
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSA3dFriction(Element* element){/*{{{*/
+
+	/*Initialize Element matrix and return if necessary*/
+	if(element->IsFloating() || !element->IsOnBase()) return NULL;
+
+	/*Build a tria element using the 3 nodes of the base of the penta. Then use 
+	 * the tria functionality to build a friction stiffness matrix on these 3
+	 * nodes: */
+	Element* basalelement = element->SpawnBasalElement();
+	ElementMatrix* Ke=CreateKMatrixSSAFriction(basalelement);
+	basalelement->DeleteMaterials(); delete basalelement;
+
+	/*clean-up and return*/
+	return Ke;
+}/*}}}*/
+ElementMatrix* StressbalanceAnalysis::CreateKMatrixSSA3dViscous(Element* element){/*{{{*/
+
+	/*Constants*/
+	const int    numdof2d=2*3;
+
+	/*Intermediaries */
+	int         i,j,approximation;
+	int         dim=3;
+	IssmDouble  Jdet,viscosity,oldviscosity,newviscosity,viscosity_overshoot;
+	IssmDouble  epsilon[5],oldepsilon[5];       /* epsilon=[exx,eyy,exy,exz,eyz];*/
+	IssmDouble  epsilons[6];                    //6 for FS
+	IssmDouble  B[3][numdof2d];
+	IssmDouble  Bprime[3][numdof2d];
+	IssmDouble  D[3][3]= {0.0};                 // material matrix, simple scalar matrix.
+	IssmDouble  D_scalar;
+	IssmDouble  Ke_gg[numdof2d][numdof2d]={0.0};
+	IssmDouble  *xyz_list  = NULL;
+
+	/*Find penta on bed as this is a SSA elements: */
+	Element* pentabase=element->GetBasalElement();
+	Element* basaltria=pentabase->SpawnBasalElement();
+
+	/*Initialize Element matrix*/
+	ElementMatrix* Ke=basaltria->NewElementMatrix(SSAApproximationEnum);
+	element->GetInputValue(&approximation,ApproximationEnum);
+
+	/*Retrieve all inputs and parameters*/
+	element->GetVerticesCoordinates(&xyz_list);
+	element->FindParam(&viscosity_overshoot,StressbalanceViscosityOvershootEnum);
+	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);
+	Input* vz_input   =element->GetInput(VzEnum);       _assert_(vz_input);
+
+	/* Start  looping on the number of gaussian points: */
+	Gauss* gauss=element->NewGauss(5);
+	Gauss* gauss_tria=new GaussTria();
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+
+		gauss->GaussPoint(ig);
+		gauss->SynchronizeGaussBase(gauss_tria);
+
+		element->JacobianDeterminant(&Jdet, xyz_list,gauss);
+		this->GetBSSA(&B[0][0],basaltria,2,xyz_list, gauss_tria);
+		this->GetBSSAprime(&Bprime[0][0],basaltria,2,xyz_list, gauss_tria);
+
+		if(approximation==SSAHOApproximationEnum){
+			element->ViscosityHO(&viscosity,dim,xyz_list,gauss,vx_input,vy_input);
+			element->ViscosityHO(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input);
+			newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
+		}
+		else if (approximation==SSAFSApproximationEnum){
+			element->ViscosityFS(&newviscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
+		}
+		else _error_("approximation " << approximation << " (" << EnumToStringx(approximation) << ") not supported yet");
+
+		D_scalar=2*newviscosity*gauss->weight*Jdet;
+		for (i=0;i<3;i++) D[i][i]=D_scalar;
+
+		TripleMultiply( &B[0][0],3,numdof2d,1,
+					&D[0][0],3,3,0,
+					&Bprime[0][0],3,numdof2d,0,
+					&Ke_gg[0][0],1);
+
+	}
+	for(i=0;i<numdof2d;i++) for(j=0;j<numdof2d;j++) Ke->values[i*numdof2d+j]+=Ke_gg[i][j];
+
+	/*Transform Coordinate System*/
+	basaltria->TransformStiffnessMatrixCoord(Ke,XYEnum);
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(xyz_list);
+	delete basaltria->material;
+	delete basaltria;
+	delete gauss_tria;
+	delete gauss;
+	return Ke;
+
+}/*}}}*/
 ElementVector* StressbalanceAnalysis::CreatePVectorCouplingHOFS(Element* element){/*{{{*/
 
 	/*compute all load vectors for this element*/
@@ -6350,38 +6292,149 @@
 	delete gauss;
 	return pe;
 }/*}}}*/
-void StressbalanceAnalysis::GetBSSAHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*NDOF2. 
-	 * For node i, Bi can be expressed in the actual coordinate system
+ElementVector* StressbalanceAnalysis::CreatePVectorHOFS(Element* element){/*{{{*/
+
+	/*compute all load vectors for this element*/
+	int init = element->FiniteElement();
+	element->SetTemporaryElementType(P1Enum);
+	ElementVector* pe1=CreatePVectorHO(element);
+	element->SetTemporaryElementType(init);
+	ElementVector* pe2=CreatePVectorFS(element);
+	int indices[3]={18,19,20};
+	element->SetTemporaryElementType(MINIcondensedEnum);
+	ElementMatrix* Ke = CreateKMatrixFS(element);
+	element->SetTemporaryElementType(init);
+	pe2->StaticCondensation(Ke,3,&indices[0]);
+	delete Ke;
+	ElementVector* pe3=CreatePVectorCouplingHOFS(element);
+	ElementVector* pe =new ElementVector(pe1,pe2,pe3);
+
+	/*clean-up and return*/
+	delete pe1;
+	delete pe2;
+	delete pe3;
+	return pe;
+}/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorSSAFS(Element* element){/*{{{*/
+
+	/*compute all load vectors for this element*/
+	int init = element->FiniteElement();
+	element->SetTemporaryElementType(P1Enum); // P1 needed for HO
+	ElementVector* pe1=CreatePVectorSSA(element);
+	element->SetTemporaryElementType(init); // P1 needed for HO
+	ElementVector* pe2=CreatePVectorFS(element);
+	int indices[3]={18,19,20};
+	element->SetTemporaryElementType(MINIcondensedEnum); // P1 needed for HO
+	ElementMatrix* Ke = CreateKMatrixFS(element);
+	element->SetTemporaryElementType(init); // P1 needed for HO
+	pe2->StaticCondensation(Ke,3,&indices[0]);
+	delete Ke;
+	ElementVector* pe3=CreatePVectorCouplingSSAFS(element);
+	ElementVector* pe =new ElementVector(pe1,pe2,pe3);
+
+	/*clean-up and return*/
+	delete pe1;
+	delete pe2;
+	delete pe3;
+	return pe;
+}/*}}}*/
+ElementVector* StressbalanceAnalysis::CreatePVectorSSAHO(Element* element){/*{{{*/
+
+	/*compute all load vectors for this element*/
+	ElementVector* pe1=CreatePVectorSSA(element);
+	ElementVector* pe2=CreatePVectorHO(element);
+	ElementVector* pe =new ElementVector(pe1,pe2);
+
+	/*clean-up and return*/
+	delete pe1;
+	delete pe2;
+	return pe;
+}/*}}}*/
+void           StressbalanceAnalysis::GetBprimeSSAFS(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+	/*Compute Bprime  matrix. Bprime=[Bprime1 Bprime2 Bprime3 Bprime4 Bprime5 Bprime6] where Bprimei is of size 5*NDOF2. 
+	 * For node i, Bprimei can be expressed in the actual coordinate system
 	 * by: 
-	 *       Bi=[ dh/dx          0      ]
-	 *          [   0           dh/dy   ]
-	 *          [ 1/2*dh/dy  1/2*dh/dx  ]
+	 *       Bprimei=[ 2*dh/dx    dh/dy   0   0 ]
+	 *               [  dh/dx    2*dh/dy  0   0 ]
+	 *               [  dh/dy     dh/dx   0   0 ]
 	 * where h is the interpolation function for node i.
 	 *
-	 * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1)
+	 * We assume Bprime has been allocated already, of size: 5x(NDOF2*NUMNODESP1)
 	 */
 
+	int    i;
+	IssmDouble dbasismini[3][7];
+
+	/*Get dbasis in actual coordinate system: */
+	element->NodalFunctionsMINIDerivatives(&dbasismini[0][0],xyz_list, gauss);
+
+	/*Build Bprime: */
+	for(i=0;i<6;i++){
+		Bprime[(3*7+6)*0+3*i+0] = 2.*dbasismini[0][i];
+		Bprime[(3*7+6)*0+3*i+1] = dbasismini[1][i];
+		Bprime[(3*7+6)*0+3*i+2] = 0.;
+		Bprime[(3*7+6)*1+3*i+0] = dbasismini[0][i];
+		Bprime[(3*7+6)*1+3*i+1] = 2.*dbasismini[1][i];
+		Bprime[(3*7+6)*1+3*i+2] = 0.;
+		Bprime[(3*7+6)*2+3*i+0] = dbasismini[1][i];
+		Bprime[(3*7+6)*2+3*i+1] = dbasismini[0][i];
+		Bprime[(3*7+6)*2+3*i+2] = 0.;
+	}
+
+	for(i=0;i<1;i++){ //Add zeros for the bubble function
+		Bprime[(3*7+6)*0+3*(6+i)+0] = 0.;
+		Bprime[(3*7+6)*0+3*(6+i)+1] = 0.;
+		Bprime[(3*7+6)*0+3*(6+i)+2] = 0.;
+		Bprime[(3*7+6)*1+3*(6+i)+0] = 0.;
+		Bprime[(3*7+6)*1+3*(6+i)+1] = 0.;
+		Bprime[(3*7+6)*1+3*(6+i)+2] = 0.;
+		Bprime[(3*7+6)*2+3*(6+i)+0] = 0.;
+		Bprime[(3*7+6)*2+3*(6+i)+1] = 0.;
+		Bprime[(3*7+6)*2+3*(6+i)+2] = 0.;
+	}
+
+	for(i=0;i<6;i++){ //last column not for the bubble function
+		Bprime[(3*7+6)*0+7*3+i] = 0.;
+		Bprime[(3*7+6)*1+7*3+i] = 0.;
+		Bprime[(3*7+6)*2+7*3+i] = 0.;
+	}
+}/*}}}*/
+void           StressbalanceAnalysis::GetBprimeSSAFSTria(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+	/*Compute Bprime  matrix. Bprime=[Bprime1 Bprime2 Bprime3] where Bprimei is of size 3*NDOF2. 
+	 * For node i, Bprimei can be expressed in the actual coordinate system
+	 * by: 
+	 *       Bprimei=[  dN/dx    0   ]
+	 *               [    0    dN/dy ]
+	 *               [  dN/dy  dN/dx ]
+	 N               [  dN/dx  dN/dy ]
+	 * where N is the finiteelement function for node i.
+	 *
+	 * We assume Bprime has been allocated already, of size: 3x(NDOF2*numnodes)
+	 */
+
+	/*Fetch number of nodes for this finite element*/
 	int numnodes = element->GetNumberOfNodes();
-	IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
 
-	/*Get dbasis in actual coordinate system: */
+	/*Get nodal functions*/
+	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
 	element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 
-	/*Build B: */
+	/*Build Bprime: */
 	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];
+		Bprime[2*numnodes*0+2*i+0] = dbasis[0*numnodes+i];
+		Bprime[2*numnodes*0+2*i+1] = 0.;
+		Bprime[2*numnodes*1+2*i+0] = 0.;
+		Bprime[2*numnodes*1+2*i+1] = 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[0*numnodes+i];
+		Bprime[2*numnodes*3+2*i+1] = dbasis[1*numnodes+i];
 	}
 
 	/*Clean-up*/
 	xDelete<IssmDouble>(dbasis);
 }/*}}}*/
-void StressbalanceAnalysis::GetBSSAFS(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void           StressbalanceAnalysis::GetBSSAFS(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: 
@@ -6439,56 +6492,7 @@
 		B[(3*7+6)*3+7*3+i] = basis[i];
 	}
 }/*}}}*/
-void StressbalanceAnalysis::GetBprimeSSAFS(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute Bprime  matrix. Bprime=[Bprime1 Bprime2 Bprime3 Bprime4 Bprime5 Bprime6] where Bprimei is of size 5*NDOF2. 
-	 * For node i, Bprimei can be expressed in the actual coordinate system
-	 * by: 
-	 *       Bprimei=[ 2*dh/dx    dh/dy   0   0 ]
-	 *               [  dh/dx    2*dh/dy  0   0 ]
-	 *               [  dh/dy     dh/dx   0   0 ]
-	 * where h is the interpolation function for node i.
-	 *
-	 * We assume Bprime has been allocated already, of size: 5x(NDOF2*NUMNODESP1)
-	 */
-
-	int    i;
-	IssmDouble dbasismini[3][7];
-
-	/*Get dbasis in actual coordinate system: */
-	element->NodalFunctionsMINIDerivatives(&dbasismini[0][0],xyz_list, gauss);
-
-	/*Build Bprime: */
-	for(i=0;i<6;i++){
-		Bprime[(3*7+6)*0+3*i+0] = 2.*dbasismini[0][i];
-		Bprime[(3*7+6)*0+3*i+1] = dbasismini[1][i];
-		Bprime[(3*7+6)*0+3*i+2] = 0.;
-		Bprime[(3*7+6)*1+3*i+0] = dbasismini[0][i];
-		Bprime[(3*7+6)*1+3*i+1] = 2.*dbasismini[1][i];
-		Bprime[(3*7+6)*1+3*i+2] = 0.;
-		Bprime[(3*7+6)*2+3*i+0] = dbasismini[1][i];
-		Bprime[(3*7+6)*2+3*i+1] = dbasismini[0][i];
-		Bprime[(3*7+6)*2+3*i+2] = 0.;
-	}
-
-	for(i=0;i<1;i++){ //Add zeros for the bubble function
-		Bprime[(3*7+6)*0+3*(6+i)+0] = 0.;
-		Bprime[(3*7+6)*0+3*(6+i)+1] = 0.;
-		Bprime[(3*7+6)*0+3*(6+i)+2] = 0.;
-		Bprime[(3*7+6)*1+3*(6+i)+0] = 0.;
-		Bprime[(3*7+6)*1+3*(6+i)+1] = 0.;
-		Bprime[(3*7+6)*1+3*(6+i)+2] = 0.;
-		Bprime[(3*7+6)*2+3*(6+i)+0] = 0.;
-		Bprime[(3*7+6)*2+3*(6+i)+1] = 0.;
-		Bprime[(3*7+6)*2+3*(6+i)+2] = 0.;
-	}
-
-	for(i=0;i<6;i++){ //last column not for the bubble function
-		Bprime[(3*7+6)*0+7*3+i] = 0.;
-		Bprime[(3*7+6)*1+7*3+i] = 0.;
-		Bprime[(3*7+6)*2+7*3+i] = 0.;
-	}
-}/*}}}*/
-void StressbalanceAnalysis::GetBSSAFSTria(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void           StressbalanceAnalysis::GetBSSAFSTria(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: 
@@ -6520,57 +6524,47 @@
 	/*Clean-up*/
 	xDelete<IssmDouble>(dbasis);
 }/*}}}*/
-void StressbalanceAnalysis::GetBprimeSSAFSTria(IssmDouble* Bprime,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	/*Compute Bprime  matrix. Bprime=[Bprime1 Bprime2 Bprime3] where Bprimei is of size 3*NDOF2. 
-	 * For node i, Bprimei can be expressed in the actual coordinate system
+void           StressbalanceAnalysis::GetBSSAHO(IssmDouble* B,Element* element,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*NDOF2. 
+	 * For node i, Bi can be expressed in the actual coordinate system
 	 * by: 
-	 *       Bprimei=[  dN/dx    0   ]
-	 *               [    0    dN/dy ]
-	 *               [  dN/dy  dN/dx ]
-	 N               [  dN/dx  dN/dy ]
-	 * where N is the finiteelement function for node i.
+	 *       Bi=[ dh/dx          0      ]
+	 *          [   0           dh/dy   ]
+	 *          [ 1/2*dh/dy  1/2*dh/dx  ]
+	 * where h is the interpolation function for node i.
 	 *
-	 * We assume Bprime has been allocated already, of size: 3x(NDOF2*numnodes)
+	 * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1)
 	 */
 
-	/*Fetch number of nodes for this finite element*/
 	int numnodes = element->GetNumberOfNodes();
+	IssmDouble* dbasis=xNew<IssmDouble>(3*numnodes);
 
-	/*Get nodal functions*/
-	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
+	/*Get dbasis in actual coordinate system: */
 	element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
 
-	/*Build Bprime: */
+	/*Build B: */
 	for(int i=0;i<numnodes;i++){
-		Bprime[2*numnodes*0+2*i+0] = dbasis[0*numnodes+i];
-		Bprime[2*numnodes*0+2*i+1] = 0.;
-		Bprime[2*numnodes*1+2*i+0] = 0.;
-		Bprime[2*numnodes*1+2*i+1] = 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[0*numnodes+i];
-		Bprime[2*numnodes*3+2*i+1] = dbasis[1*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::GetLSSAFS(IssmDouble* LFS,Element* element,Gauss* gauss_in){/*{{{*/
-	/*
-	 * Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof. 
-	 * For node i, Li can be expressed in the actual coordinate system
+void           StressbalanceAnalysis::GetLprimeFSSSA(IssmDouble* LprimeFS,Element* element,IssmDouble* xyz_list,Gauss* gauss_in){/*{{{*/
+	/* Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. 
+	 * For node i, Lpi can be expressed in the actual coordinate system
 	 * by: 
-	 *       Li=[ h    0 ]
-	 *	 	      [ 0    h ]
-	 *	 	      [ h    0 ]
-	 *	 	      [ 0    h ]
-	 *	 	      [ h    0 ]
-	 *	 	      [ 0    h ]
-	 *	 	      [ h    0 ]
-	 *	 	      [ 0    h ]
+	 *       Lpi=[ h    0 ]
+	 *		       [ 0    h ]
+	 *		       [ h    0 ]
+	 *		       [ 0    h ]
 	 * where h is the interpolation function for node i.
 	 */
-
 	int num_dof=2;
 	IssmDouble basis[3];
 
@@ -6583,27 +6577,19 @@
 	basis[1]=gauss->coord2*(1-gauss->coord4)/2.0;
 	basis[2]=gauss->coord3*(1-gauss->coord4)/2.0;
 
-	/*Build LFS: */
+	/*Build LprimeFS: */
 	for(int i=0;i<3;i++){
-		LFS[num_dof*3*0+num_dof*i+0] = basis[i];
-		LFS[num_dof*3*0+num_dof*i+1] = 0;
-		LFS[num_dof*3*1+num_dof*i+0] = 0;
-		LFS[num_dof*3*1+num_dof*i+1] = basis[i];
-		LFS[num_dof*3*2+num_dof*i+0] = basis[i];
-		LFS[num_dof*3*2+num_dof*i+1] = 0;
-		LFS[num_dof*3*3+num_dof*i+0] = 0;
-		LFS[num_dof*3*3+num_dof*i+1] = basis[i];
-		LFS[num_dof*3*4+num_dof*i+0] = basis[i];
-		LFS[num_dof*3*4+num_dof*i+1] = 0;
-		LFS[num_dof*3*5+num_dof*i+0] = 0;
-		LFS[num_dof*3*5+num_dof*i+1] = basis[i];
-		LFS[num_dof*3*6+num_dof*i+0] = basis[i];
-		LFS[num_dof*3*6+num_dof*i+1] = 0;
-		LFS[num_dof*3*7+num_dof*i+0] = 0;
-		LFS[num_dof*3*7+num_dof*i+1] = basis[i];
+		LprimeFS[num_dof*3*0+num_dof*i+0] = basis[i];
+		LprimeFS[num_dof*3*0+num_dof*i+1] = 0.;
+		LprimeFS[num_dof*3*1+num_dof*i+0] = 0.;
+		LprimeFS[num_dof*3*1+num_dof*i+1] = basis[i];
+		LprimeFS[num_dof*3*2+num_dof*i+0] = basis[i];
+		LprimeFS[num_dof*3*2+num_dof*i+1] = 0.;
+		LprimeFS[num_dof*3*3+num_dof*i+0] = 0.;
+		LprimeFS[num_dof*3*3+num_dof*i+1] = basis[i];
 	}
 }/*}}}*/
-void StressbalanceAnalysis::GetLprimeSSAFS(IssmDouble* LprimeFS,Element* element,IssmDouble* xyz_list,Gauss* gauss_in){/*{{{*/
+void           StressbalanceAnalysis::GetLprimeSSAFS(IssmDouble* LprimeFS,Element* element,IssmDouble* xyz_list,Gauss* gauss_in){/*{{{*/
 	/* Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. 
 	 * For node i, Lpi can be expressed in the actual coordinate system
 	 * by: 
@@ -6708,7 +6694,7 @@
 		LprimeFS[num_dof_total*7+num_dof_vel+i] = 0.;
 	}
 }/*}}}*/
-void StressbalanceAnalysis::GetLFSSSA(IssmDouble* LFS,Element* element,Gauss* gauss_in){/*{{{*/
+void           StressbalanceAnalysis::GetLFSSSA(IssmDouble* LFS,Element* element,Gauss* gauss_in){/*{{{*/
 	/* Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof. 
 	 * For node i, Li can be expressed in the actual coordinate system
 	 * by: 
@@ -6747,16 +6733,22 @@
 		LFS[num_dof*3*3+num_dof*i+2] = basis[i];
 	}
 }/*}}}*/
-void StressbalanceAnalysis::GetLprimeFSSSA(IssmDouble* LprimeFS,Element* element,IssmDouble* xyz_list,Gauss* gauss_in){/*{{{*/
-	/* Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. 
-	 * For node i, Lpi can be expressed in the actual coordinate system
+void           StressbalanceAnalysis::GetLSSAFS(IssmDouble* LFS,Element* element,Gauss* gauss_in){/*{{{*/
+	/*
+	 * Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof. 
+	 * For node i, Li can be expressed in the actual coordinate system
 	 * by: 
-	 *       Lpi=[ h    0 ]
-	 *		       [ 0    h ]
-	 *		       [ h    0 ]
-	 *		       [ 0    h ]
+	 *       Li=[ h    0 ]
+	 *	 	      [ 0    h ]
+	 *	 	      [ h    0 ]
+	 *	 	      [ 0    h ]
+	 *	 	      [ h    0 ]
+	 *	 	      [ 0    h ]
+	 *	 	      [ h    0 ]
+	 *	 	      [ 0    h ]
 	 * where h is the interpolation function for node i.
 	 */
+
 	int num_dof=2;
 	IssmDouble basis[3];
 
@@ -6769,19 +6761,27 @@
 	basis[1]=gauss->coord2*(1-gauss->coord4)/2.0;
 	basis[2]=gauss->coord3*(1-gauss->coord4)/2.0;
 
-	/*Build LprimeFS: */
+	/*Build LFS: */
 	for(int i=0;i<3;i++){
-		LprimeFS[num_dof*3*0+num_dof*i+0] = basis[i];
-		LprimeFS[num_dof*3*0+num_dof*i+1] = 0.;
-		LprimeFS[num_dof*3*1+num_dof*i+0] = 0.;
-		LprimeFS[num_dof*3*1+num_dof*i+1] = basis[i];
-		LprimeFS[num_dof*3*2+num_dof*i+0] = basis[i];
-		LprimeFS[num_dof*3*2+num_dof*i+1] = 0.;
-		LprimeFS[num_dof*3*3+num_dof*i+0] = 0.;
-		LprimeFS[num_dof*3*3+num_dof*i+1] = basis[i];
+		LFS[num_dof*3*0+num_dof*i+0] = basis[i];
+		LFS[num_dof*3*0+num_dof*i+1] = 0;
+		LFS[num_dof*3*1+num_dof*i+0] = 0;
+		LFS[num_dof*3*1+num_dof*i+1] = basis[i];
+		LFS[num_dof*3*2+num_dof*i+0] = basis[i];
+		LFS[num_dof*3*2+num_dof*i+1] = 0;
+		LFS[num_dof*3*3+num_dof*i+0] = 0;
+		LFS[num_dof*3*3+num_dof*i+1] = basis[i];
+		LFS[num_dof*3*4+num_dof*i+0] = basis[i];
+		LFS[num_dof*3*4+num_dof*i+1] = 0;
+		LFS[num_dof*3*5+num_dof*i+0] = 0;
+		LFS[num_dof*3*5+num_dof*i+1] = basis[i];
+		LFS[num_dof*3*6+num_dof*i+0] = basis[i];
+		LFS[num_dof*3*6+num_dof*i+1] = 0;
+		LFS[num_dof*3*7+num_dof*i+0] = 0;
+		LFS[num_dof*3*7+num_dof*i+1] = basis[i];
 	}
 }/*}}}*/
-void StressbalanceAnalysis::InputUpdateFromSolutionHOFS(IssmDouble* solution,Element* element){/*{{{*/
+void           StressbalanceAnalysis::InputUpdateFromSolutionHOFS(IssmDouble* solution,Element* element){/*{{{*/
 
 	int         i;
 	IssmDouble  rho_ice,g,FSreconditioning;
@@ -6878,7 +6878,7 @@
 	xDelete<int>(doflistHO);
 	xDelete<int>(cs_list);
 }/*}}}*/
-void StressbalanceAnalysis::InputUpdateFromSolutionSSAFS(IssmDouble* solution,Element* element){/*{{{*/
+void           StressbalanceAnalysis::InputUpdateFromSolutionSSAFS(IssmDouble* solution,Element* element){/*{{{*/
 
 	int         i;
 	IssmDouble  rho_ice,g,FSreconditioning;
@@ -6984,7 +6984,7 @@
 	xDelete<int>(doflistSSA);
 	xDelete<int>(cs_list);
 }/*}}}*/
-void StressbalanceAnalysis::InputUpdateFromSolutionSSAHO(IssmDouble* solution,Element* element){/*{{{*/
+void           StressbalanceAnalysis::InputUpdateFromSolutionSSAHO(IssmDouble* solution,Element* element){/*{{{*/
 
 	int         i,domaintype;
 	IssmDouble  rho_ice,g;
