Index: issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15556)
+++ issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15564)
@@ -217,5 +217,5 @@
 	IssmDouble  xyz_list[NUMVERTICES][3];
 	IssmDouble  xyz_list_tria[3][3];
-	IssmDouble  rho_ice,gravity,stokesreconditioning;
+	IssmDouble  rho_ice,gravity,FSreconditioning;
 	IssmDouble  pressure,viscosity,Jdet2d;
 	IssmDouble  bed_normal[3];
@@ -234,8 +234,8 @@
 	/*Check analysis_types*/
 	if (analysis_type!=DiagnosticHorizAnalysisEnum) _error_("Not supported yet!");
-	if (approximation!=StokesApproximationEnum) _error_("Not supported yet!");
+	if (approximation!=FSApproximationEnum) _error_("Not supported yet!");
 
 	/*retrieve some parameters: */
-	this->parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
+	this->parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
 
 	if(!IsOnBed()){
@@ -267,11 +267,11 @@
 		/*Compute strain rate viscosity and pressure: */
 		this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
-		material->GetViscosity3dStokes(&viscosity,&epsilon[0]);
+		material->GetViscosity3dFS(&viscosity,&epsilon[0]);
 		pressure_input->GetInputValue(&pressure,gauss);
 
 		/*Compute Stress*/
-		sigma_xx=2*viscosity*epsilon[0]-pressure*stokesreconditioning; // sigma = nu eps - pressure
-		sigma_yy=2*viscosity*epsilon[1]-pressure*stokesreconditioning;
-		sigma_zz=2*viscosity*epsilon[2]-pressure*stokesreconditioning;
+		sigma_xx=2*viscosity*epsilon[0]-pressure*FSreconditioning; // sigma = nu eps - pressure
+		sigma_yy=2*viscosity*epsilon[1]-pressure*FSreconditioning;
+		sigma_zz=2*viscosity*epsilon[2]-pressure*FSreconditioning;
 		sigma_xy=2*viscosity*epsilon[3];
 		sigma_xz=2*viscosity*epsilon[4];
@@ -415,6 +415,6 @@
 			Ke=CreateKMatrixAdjointHoriz();
 			break;
-		case DiagnosticHutterAnalysisEnum:
-			Ke=CreateKMatrixDiagnosticHutter();
+		case DiagnosticSIAAnalysisEnum:
+			Ke=CreateKMatrixDiagnosticSIA();
 			break;
 		case DiagnosticVertAnalysisEnum:
@@ -546,6 +546,6 @@
 			pe=CreatePVectorDiagnosticHoriz();
 			break;
-		case DiagnosticHutterAnalysisEnum:
-			pe=CreatePVectorDiagnosticHutter();
+		case DiagnosticSIAAnalysisEnum:
+			pe=CreatePVectorDiagnosticSIA();
 			break;
 		case DiagnosticVertAnalysisEnum:
@@ -1190,16 +1190,16 @@
 		int approximation;
 		inputs->GetInputValue(&approximation,ApproximationEnum);
-		if(approximation==StokesApproximationEnum || approximation==NoneApproximationEnum){
-			GetSolutionFromInputsDiagnosticStokes(solution);
-		}
-		else if (approximation==MacAyealApproximationEnum || approximation==PattynApproximationEnum || approximation==HutterApproximationEnum){
+		if(approximation==FSApproximationEnum || approximation==NoneApproximationEnum){
+			GetSolutionFromInputsDiagnosticFS(solution);
+		}
+		else if (approximation==MacAyealApproximationEnum || approximation==HOApproximationEnum || approximation==SIAApproximationEnum){
 			GetSolutionFromInputsDiagnosticHoriz(solution);
 		}
-		else if (approximation==MacAyealPattynApproximationEnum || approximation==PattynStokesApproximationEnum || approximation==MacAyealStokesApproximationEnum){
+		else if (approximation==MacAyealHOApproximationEnum || approximation==HOFSApproximationEnum || approximation==MacAyealFSApproximationEnum){
 			return; //the elements around will create the solution
 		}
 		break;
-	case DiagnosticHutterAnalysisEnum:
-		GetSolutionFromInputsDiagnosticHutter(solution);
+	case DiagnosticSIAAnalysisEnum:
+		GetSolutionFromInputsDiagnosticSIA(solution);
 		break;
 	case DiagnosticVertAnalysisEnum:
@@ -1249,7 +1249,7 @@
 }
 /*}}}*/
-/*FUNCTION Penta::GetStrainRate3dPattyn{{{*/
-void Penta::GetStrainRate3dPattyn(IssmDouble* epsilon,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input){
-	/*Compute the 3d Blatter/PattynStrain Rate (5 components):
+/*FUNCTION Penta::GetStrainRate3dHO{{{*/
+void Penta::GetStrainRate3dHO(IssmDouble* epsilon,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input){
+	/*Compute the 3d Blatter/HOStrain Rate (5 components):
 	 *
 	 * epsilon=[exx eyy exy exz eyz]
@@ -1271,6 +1271,6 @@
 
 	/*Get strain rate assuming that epsilon has been allocated*/
-	vx_input->GetVxStrainRate3dPattyn(epsilonvx,xyz_list,gauss);
-	vy_input->GetVyStrainRate3dPattyn(epsilonvy,xyz_list,gauss);
+	vx_input->GetVxStrainRate3dHO(epsilonvx,xyz_list,gauss);
+	vy_input->GetVyStrainRate3dHO(epsilonvy,xyz_list,gauss);
 
 	/*Sum all contributions*/
@@ -1973,24 +1973,24 @@
 			this->inputs->AddInput(new IntInput(ApproximationEnum,MacAyealApproximationEnum));
 		}
-		else if (iomodel->Data(FlowequationElementEquationEnum)[index]==PattynApproximationEnum){
-			this->inputs->AddInput(new IntInput(ApproximationEnum,PattynApproximationEnum));
-		}
-		else if (iomodel->Data(FlowequationElementEquationEnum)[index]==MacAyealPattynApproximationEnum){
-			this->inputs->AddInput(new IntInput(ApproximationEnum,MacAyealPattynApproximationEnum));
-		}
-		else if (iomodel->Data(FlowequationElementEquationEnum)[index]==HutterApproximationEnum){
-			this->inputs->AddInput(new IntInput(ApproximationEnum,HutterApproximationEnum));
+		else if (iomodel->Data(FlowequationElementEquationEnum)[index]==HOApproximationEnum){
+			this->inputs->AddInput(new IntInput(ApproximationEnum,HOApproximationEnum));
+		}
+		else if (iomodel->Data(FlowequationElementEquationEnum)[index]==MacAyealHOApproximationEnum){
+			this->inputs->AddInput(new IntInput(ApproximationEnum,MacAyealHOApproximationEnum));
+		}
+		else if (iomodel->Data(FlowequationElementEquationEnum)[index]==SIAApproximationEnum){
+			this->inputs->AddInput(new IntInput(ApproximationEnum,SIAApproximationEnum));
 		}
 		else if (iomodel->Data(FlowequationElementEquationEnum)[index]==L1L2ApproximationEnum){
 			this->inputs->AddInput(new IntInput(ApproximationEnum,L1L2ApproximationEnum));
 		}
-		else if (iomodel->Data(FlowequationElementEquationEnum)[index]==StokesApproximationEnum){
-			this->inputs->AddInput(new IntInput(ApproximationEnum,StokesApproximationEnum));
-		}
-		else if (iomodel->Data(FlowequationElementEquationEnum)[index]==MacAyealStokesApproximationEnum){
-			this->inputs->AddInput(new IntInput(ApproximationEnum,MacAyealStokesApproximationEnum));
-		}
-		else if (iomodel->Data(FlowequationElementEquationEnum)[index]==PattynStokesApproximationEnum){
-			this->inputs->AddInput(new IntInput(ApproximationEnum,PattynStokesApproximationEnum));
+		else if (iomodel->Data(FlowequationElementEquationEnum)[index]==FSApproximationEnum){
+			this->inputs->AddInput(new IntInput(ApproximationEnum,FSApproximationEnum));
+		}
+		else if (iomodel->Data(FlowequationElementEquationEnum)[index]==MacAyealFSApproximationEnum){
+			this->inputs->AddInput(new IntInput(ApproximationEnum,MacAyealFSApproximationEnum));
+		}
+		else if (iomodel->Data(FlowequationElementEquationEnum)[index]==HOFSApproximationEnum){
+			this->inputs->AddInput(new IntInput(ApproximationEnum,HOFSApproximationEnum));
 		}
 		else if (iomodel->Data(FlowequationElementEquationEnum)[index]==NoneApproximationEnum){
@@ -2031,6 +2031,6 @@
 		InputUpdateFromSolutionDiagnosticHoriz( solution);
 		break;
-	case DiagnosticHutterAnalysisEnum:
-		InputUpdateFromSolutionDiagnosticHutter( solution);
+	case DiagnosticSIAAnalysisEnum:
+		InputUpdateFromSolutionDiagnosticSIA( solution);
 		break;
 	case DiagnosticVertAnalysisEnum:
@@ -2042,6 +2042,6 @@
 		int approximation;
 		inputs->GetInputValue(&approximation,ApproximationEnum);
-		if(approximation==StokesApproximationEnum || approximation==NoneApproximationEnum){
-			InputUpdateFromSolutionAdjointStokes( solution);
+		if(approximation==FSApproximationEnum || approximation==NoneApproximationEnum){
+			InputUpdateFromSolutionAdjointFS( solution);
 		}
 		else{
@@ -2651,6 +2651,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::ReduceMatrixStokes {{{*/
-void Penta::ReduceMatrixStokes(IssmDouble* Ke_reduced, IssmDouble* Ke_temp){
+/*FUNCTION Penta::ReduceMatrixFS {{{*/
+void Penta::ReduceMatrixFS(IssmDouble* Ke_reduced, IssmDouble* Ke_temp){
 
 	int    i,j;
@@ -2693,6 +2693,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::ReduceVectorStokes {{{*/
-void Penta::ReduceVectorStokes(IssmDouble* Pe_reduced, IssmDouble* Ke_temp, IssmDouble* Pe_temp){
+/*FUNCTION Penta::ReduceVectorFS {{{*/
+void Penta::ReduceVectorFS(IssmDouble* Pe_reduced, IssmDouble* Ke_temp, IssmDouble* Pe_temp){
 
 	int    i,j;
@@ -2786,7 +2786,7 @@
 	IssmDouble xz_plane[6];
 
-	/*For Stokes only: we want the CS to be tangential to the bedrock*/
+	/*For FS only: we want the CS to be tangential to the bedrock*/
 	inputs->GetInputValue(&approximation,ApproximationEnum);
-	if(IsFloating() || !IsOnBed() || (approximation!=StokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum &&  approximation!=PattynStokesApproximationEnum)) return;
+	if(IsFloating() || !IsOnBed() || (approximation!=FSApproximationEnum && approximation!=MacAyealFSApproximationEnum &&  approximation!=HOFSApproximationEnum)) return;
 
 	/*Get slope on each node*/
@@ -2883,5 +2883,5 @@
 	if(analysis_type==DiagnosticHorizAnalysisEnum){
 		inputs->GetInputValue(&approximation,ApproximationEnum);
-		if(approximation==MacAyealPattynApproximationEnum || approximation==MacAyealStokesApproximationEnum){
+		if(approximation==MacAyealHOApproximationEnum || approximation==MacAyealFSApproximationEnum){
 			parameters->FindParam(&numlayers,MeshNumberoflayersEnum);
 			o_nz += numlayers*3;
@@ -3102,5 +3102,5 @@
 	int        stabilization;
 	bool       dakota_analysis;
-	bool       isstokes;
+	bool       isFS;
 	IssmDouble beta,heatcapacity,referencetemperature,meltingpoint,latentheat;
 
@@ -3109,5 +3109,5 @@
 	iomodel->Constant(&stabilization,PrognosticStabilizationEnum);
 	iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum);
-	iomodel->Constant(&isstokes,FlowequationIsstokesEnum);
+	iomodel->Constant(&isFS,FlowequationIsFSEnum);
 	iomodel->Constant(&beta,MaterialsBetaEnum);
 	iomodel->Constant(&heatcapacity,MaterialsHeatcapacityEnum);
@@ -3175,34 +3175,34 @@
 					this->inputs->AddInput(new PentaInput(QmuPressureEnum,nodeinputs,P1Enum));
 				}
-				if(isstokes){
+				if(isFS){
 					this->inputs->AddInput(new PentaInput(PressureEnum,nodeinputs,P1Enum));
 					this->inputs->AddInput(new PentaInput(PressurePicardEnum,nodeinputs,P1Enum));
 				}
 			}
-			if(*(iomodel->Data(FlowequationElementEquationEnum)+index)==PattynStokesApproximationEnum){
-				/*Create VzPattyn and VzStokes Enums*/
-				if(iomodel->Data(VzEnum) && iomodel->Data(FlowequationBorderstokesEnum)){
-					for(i=0;i<6;i++) nodeinputs[i]=iomodel->Data(VzEnum)[penta_vertex_ids[i]-1]*iomodel->Data(FlowequationBorderstokesEnum)[penta_vertex_ids[i]-1];
-					this->inputs->AddInput(new PentaInput(VzStokesEnum,nodeinputs,P1Enum));
-					for(i=0;i<6;i++) nodeinputs[i]=iomodel->Data(VzEnum)[penta_vertex_ids[i]-1]*(1-iomodel->Data(FlowequationBorderstokesEnum)[penta_vertex_ids[i]-1]);
-					this->inputs->AddInput(new PentaInput(VzPattynEnum,nodeinputs,P1Enum));
+			if(*(iomodel->Data(FlowequationElementEquationEnum)+index)==HOFSApproximationEnum){
+				/*Create VzHO and VzFS Enums*/
+				if(iomodel->Data(VzEnum) && iomodel->Data(FlowequationBorderFSEnum)){
+					for(i=0;i<6;i++) nodeinputs[i]=iomodel->Data(VzEnum)[penta_vertex_ids[i]-1]*iomodel->Data(FlowequationBorderFSEnum)[penta_vertex_ids[i]-1];
+					this->inputs->AddInput(new PentaInput(VzFSEnum,nodeinputs,P1Enum));
+					for(i=0;i<6;i++) nodeinputs[i]=iomodel->Data(VzEnum)[penta_vertex_ids[i]-1]*(1-iomodel->Data(FlowequationBorderFSEnum)[penta_vertex_ids[i]-1]);
+					this->inputs->AddInput(new PentaInput(VzHOEnum,nodeinputs,P1Enum));
 				}
 				else{
 					for(i=0;i<6;i++)nodeinputs[i]=0;
-					this->inputs->AddInput(new PentaInput(VzStokesEnum,nodeinputs,P1Enum));
-					this->inputs->AddInput(new PentaInput(VzPattynEnum,nodeinputs,P1Enum));
+					this->inputs->AddInput(new PentaInput(VzFSEnum,nodeinputs,P1Enum));
+					this->inputs->AddInput(new PentaInput(VzHOEnum,nodeinputs,P1Enum));
 				}
 			}
-			if(*(iomodel->Data(FlowequationElementEquationEnum)+index)==MacAyealStokesApproximationEnum){
-				/*Create VzMacAyeal and VzStokes Enums*/
-				if(iomodel->Data(VzEnum) && iomodel->Data(FlowequationBorderstokesEnum)){
-					for(i=0;i<6;i++) nodeinputs[i]=iomodel->Data(VzEnum)[penta_vertex_ids[i]-1]*iomodel->Data(FlowequationBorderstokesEnum)[penta_vertex_ids[i]-1];
-					this->inputs->AddInput(new PentaInput(VzStokesEnum,nodeinputs,P1Enum));
-					for(i=0;i<6;i++) nodeinputs[i]=iomodel->Data(VzEnum)[penta_vertex_ids[i]-1]*(1-iomodel->Data(FlowequationBorderstokesEnum)[penta_vertex_ids[i]-1]);
+			if(*(iomodel->Data(FlowequationElementEquationEnum)+index)==MacAyealFSApproximationEnum){
+				/*Create VzMacAyeal and VzFS Enums*/
+				if(iomodel->Data(VzEnum) && iomodel->Data(FlowequationBorderFSEnum)){
+					for(i=0;i<6;i++) nodeinputs[i]=iomodel->Data(VzEnum)[penta_vertex_ids[i]-1]*iomodel->Data(FlowequationBorderFSEnum)[penta_vertex_ids[i]-1];
+					this->inputs->AddInput(new PentaInput(VzFSEnum,nodeinputs,P1Enum));
+					for(i=0;i<6;i++) nodeinputs[i]=iomodel->Data(VzEnum)[penta_vertex_ids[i]-1]*(1-iomodel->Data(FlowequationBorderFSEnum)[penta_vertex_ids[i]-1]);
 					this->inputs->AddInput(new PentaInput(VzMacAyealEnum,nodeinputs,P1Enum));
 				}
 				else{
 					for(i=0;i<6;i++)nodeinputs[i]=0;
-					this->inputs->AddInput(new PentaInput(VzStokesEnum,nodeinputs,P1Enum));
+					this->inputs->AddInput(new PentaInput(VzFSEnum,nodeinputs,P1Enum));
 					this->inputs->AddInput(new PentaInput(VzMacAyealEnum,nodeinputs,P1Enum));
 				}
@@ -3289,5 +3289,5 @@
 
 		this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
-		material->GetViscosity3dStokes(&viscosity,&epsilon[0]);
+		material->GetViscosity3dFS(&viscosity,&epsilon[0]);
 		GetPhi(&phi, &epsilon[0], viscosity);
 
@@ -4067,5 +4067,5 @@
 
 		this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
-		material->GetViscosity3dStokes(&viscosity,&epsilon[0]);
+		material->GetViscosity3dFS(&viscosity,&epsilon[0]);
 		GetPhi(&phi, &epsilon[0], viscosity);
 
@@ -4323,5 +4323,5 @@
 
 		this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
-		material->GetViscosity3dStokes(&viscosity,&epsilon[0]);
+		material->GetViscosity3dFS(&viscosity,&epsilon[0]);
 		GetPhi(&phi, &epsilon[0], viscosity);
 
@@ -4776,8 +4776,8 @@
 		case MacAyealApproximationEnum:
 			return CreateKMatrixAdjointMacAyeal2d();
-		case PattynApproximationEnum:
-			return CreateKMatrixAdjointPattyn();
-		case StokesApproximationEnum:
-			return CreateKMatrixAdjointStokes();
+		case HOApproximationEnum:
+			return CreateKMatrixAdjointHO();
+		case FSApproximationEnum:
+			return CreateKMatrixAdjointFS();
 		case NoneApproximationEnum:
 			return NULL;
@@ -4821,6 +4821,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreateKMatrixAdjointPattyn{{{*/
-ElementMatrix* Penta::CreateKMatrixAdjointPattyn(void){
+/*FUNCTION Penta::CreateKMatrixAdjointHO{{{*/
+ElementMatrix* Penta::CreateKMatrixAdjointHO(void){
 
 	/*Intermediaries */
@@ -4838,7 +4838,7 @@
 	GaussPenta *gauss=NULL;
 
-	/*Initialize Jacobian with regular Pattyn (first part of the Gateau derivative)*/
+	/*Initialize Jacobian with regular HO (first part of the Gateau derivative)*/
 	parameters->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum);
-	ElementMatrix* Ke=CreateKMatrixDiagnosticPattyn();
+	ElementMatrix* Ke=CreateKMatrixDiagnosticHO();
 	if(incomplete_adjoint) return Ke;
 
@@ -4857,5 +4857,5 @@
 		GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss);
 
-		this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
+		this->GetStrainRate3dHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
 		material->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);
 		eps1[0]=2*epsilon[0]+epsilon[1];   eps2[0]=epsilon[2];
@@ -4886,6 +4886,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreateKMatrixAdjointStokes{{{*/
-ElementMatrix* Penta::CreateKMatrixAdjointStokes(void){
+/*FUNCTION Penta::CreateKMatrixAdjointFS{{{*/
+ElementMatrix* Penta::CreateKMatrixAdjointFS(void){
 
 	/*Constants*/
@@ -4906,7 +4906,7 @@
 	GaussPenta *gauss=NULL;
 
-	/*Initialize Jacobian with regular Stokes (first part of the Gateau derivative)*/
+	/*Initialize Jacobian with regular FS (first part of the Gateau derivative)*/
 	parameters->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum);
-	ElementMatrix* Ke=CreateKMatrixDiagnosticStokes();
+	ElementMatrix* Ke=CreateKMatrixDiagnosticFS();
 	if(incomplete_adjoint) return Ke;
 
@@ -4926,5 +4926,5 @@
 		GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss);
 
-		this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
+		this->GetStrainRate3dHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
 		material->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);
 		eps1[0]=epsilon[0];   eps2[0]=epsilon[2];   eps3[0]=epsilon[3];
@@ -4973,10 +4973,10 @@
 		case MacAyealApproximationEnum:
 			return CreatePVectorAdjointMacAyeal();
-		case PattynApproximationEnum:
-			return CreatePVectorAdjointPattyn();
+		case HOApproximationEnum:
+			return CreatePVectorAdjointHO();
 		case NoneApproximationEnum:
 			return NULL;
-		case StokesApproximationEnum:
-			return CreatePVectorAdjointStokes();
+		case FSApproximationEnum:
+			return CreatePVectorAdjointFS();
 		default:
 			_error_("Approximation " << EnumToStringx(approximation) << " not supported yet");
@@ -4998,6 +4998,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreatePVectorAdjointPattyn{{{*/
-ElementVector* Penta::CreatePVectorAdjointPattyn(void){
+/*FUNCTION Penta::CreatePVectorAdjointHO{{{*/
+ElementVector* Penta::CreatePVectorAdjointHO(void){
 
 	if (!IsOnSurface()) return NULL;
@@ -5012,6 +5012,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreatePVectorAdjointStokes{{{*/
-ElementVector* Penta::CreatePVectorAdjointStokes(void){
+/*FUNCTION Penta::CreatePVectorAdjointFS{{{*/
+ElementVector* Penta::CreatePVectorAdjointFS(void){
 
 	if (!IsOnSurface()) return NULL;
@@ -5019,5 +5019,5 @@
 	/*Call Tria function*/
 	Tria* tria=(Tria*)SpawnTria(3,4,5); //nodes 3, 4 and 5 make the new tria (upper face).
-	ElementVector* pe=tria->CreatePVectorAdjointStokes();
+	ElementVector* pe=tria->CreatePVectorAdjointFS();
 	delete tria->material; delete tria;
 
@@ -5059,9 +5059,9 @@
 					GradjDragMacAyeal(gradient,control_index);
 					break;
-				case PattynApproximationEnum:
-					GradjDragPattyn(gradient,control_index);
+				case HOApproximationEnum:
+					GradjDragHO(gradient,control_index);
 					break;
-				case StokesApproximationEnum:
-					GradjDragStokes(gradient,control_index);
+				case FSApproximationEnum:
+					GradjDragFS(gradient,control_index);
 					break;
 				case NoneApproximationEnum:
@@ -5079,9 +5079,9 @@
 					GradjBbarMacAyeal(gradient,control_index);
 					break;
-				case PattynApproximationEnum:
-					GradjBbarPattyn(gradient,control_index);
+				case HOApproximationEnum:
+					GradjBbarHO(gradient,control_index);
 					break;
-				case StokesApproximationEnum:
-					GradjBbarStokes(gradient,control_index);
+				case FSApproximationEnum:
+					GradjBbarFS(gradient,control_index);
 					break;
 				case NoneApproximationEnum:
@@ -5143,6 +5143,6 @@
 
 } /*}}}*/
-/*FUNCTION Penta::GradjDragPattyn {{{*/
-void  Penta::GradjDragPattyn(Vector<IssmDouble>* gradient,int control_index){
+/*FUNCTION Penta::GradjDragHO {{{*/
+void  Penta::GradjDragHO(Vector<IssmDouble>* gradient,int control_index){
 
 	int        i,j;
@@ -5214,6 +5214,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::GradjDragStokes {{{*/
-void  Penta::GradjDragStokes(Vector<IssmDouble>* gradient,int control_index){
+/*FUNCTION Penta::GradjDragFS {{{*/
+void  Penta::GradjDragFS(Vector<IssmDouble>* gradient,int control_index){
 
 	int        i,j;
@@ -5324,6 +5324,6 @@
 
 } /*}}}*/
-/*FUNCTION Penta::GradjBbarPattyn {{{*/
-void  Penta::GradjBbarPattyn(Vector<IssmDouble>* gradient,int control_index){
+/*FUNCTION Penta::GradjBbarHO {{{*/
+void  Penta::GradjBbarHO(Vector<IssmDouble>* gradient,int control_index){
 
 	/*Gradient is computed on bed only (Bbar)*/
@@ -5341,6 +5341,6 @@
 	this->material->inputs->DeleteInput(MaterialsRheologyBbarEnum);
 } /*}}}*/
-/*FUNCTION Penta::GradjBbarStokes {{{*/
-void  Penta::GradjBbarStokes(Vector<IssmDouble>* gradient,int control_index){
+/*FUNCTION Penta::GradjBbarFS {{{*/
+void  Penta::GradjBbarFS(Vector<IssmDouble>* gradient,int control_index){
 
 	/*Gradient is computed on bed only (Bbar)*/
@@ -5403,6 +5403,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::InputUpdateFromSolutionAdjointStokes {{{*/
-void  Penta::InputUpdateFromSolutionAdjointStokes(IssmDouble* solution){
+/*FUNCTION Penta::InputUpdateFromSolutionAdjointFS {{{*/
+void  Penta::InputUpdateFromSolutionAdjointFS(IssmDouble* solution){
 
 	const int    numdof=NDOF4*NUMVERTICES;
@@ -5951,13 +5951,13 @@
 
 	switch(approximation){
-		case StokesApproximationEnum:
-			return CreateDVectorDiagnosticStokes();
+		case FSApproximationEnum:
+			return CreateDVectorDiagnosticFS();
 		default:
-			return NULL; //no need for doftypes outside of stokes approximation
-	}
-}
-/*}}}*/
-/*FUNCTION Penta::CreateDVectorDiagnosticStokes{{{*/
-ElementVector* Penta::CreateDVectorDiagnosticStokes(void){
+			return NULL; //no need for doftypes outside of FS approximation
+	}
+}
+/*}}}*/
+/*FUNCTION Penta::CreateDVectorDiagnosticFS{{{*/
+ElementVector* Penta::CreateDVectorDiagnosticFS(void){
 
 	/*output: */
@@ -5969,7 +5969,7 @@
 	/*Initialize Element vector and return if necessary*/
 	inputs->GetInputValue(&approximation,ApproximationEnum);
-	if(approximation!=StokesApproximationEnum) return NULL;
-
-	De=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
+	if(approximation!=FSApproximationEnum) return NULL;
+
+	De=new ElementVector(nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
 
 	for (i=0;i<NUMVERTICES;i++){
@@ -5983,10 +5983,10 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreateKMatrixCouplingMacAyealPattyn{{{*/
-ElementMatrix* Penta::CreateKMatrixCouplingMacAyealPattyn(void){
+/*FUNCTION Penta::CreateKMatrixCouplingMacAyealHO{{{*/
+ElementMatrix* Penta::CreateKMatrixCouplingMacAyealHO(void){
 
 	/*compute all stiffness matrices for this element*/
-	ElementMatrix* Ke1=CreateKMatrixCouplingMacAyealPattynViscous();
-	ElementMatrix* Ke2=CreateKMatrixCouplingMacAyealPattynFriction();
+	ElementMatrix* Ke1=CreateKMatrixCouplingMacAyealHOViscous();
+	ElementMatrix* Ke2=CreateKMatrixCouplingMacAyealHOFriction();
 	ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
 
@@ -5997,6 +5997,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreateKMatrixCouplingMacAyealPattynViscous{{{*/
-ElementMatrix* Penta::CreateKMatrixCouplingMacAyealPattynViscous(void){
+/*FUNCTION Penta::CreateKMatrixCouplingMacAyealHOViscous{{{*/
+ElementMatrix* Penta::CreateKMatrixCouplingMacAyealHOViscous(void){
 
 	/*Constants*/
@@ -6023,5 +6023,5 @@
 	int         cs_list[numnodes];
 
-	/*Find penta on bed as pattyn must be coupled to the dofs on the bed: */
+	/*Find penta on bed as HO must be coupled to the dofs on the bed: */
 	Penta* pentabase=GetBasalElement();
 	Tria*  tria=pentabase->SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
@@ -6037,5 +6037,5 @@
 	/*Initialize Element matrix*/
 	ElementMatrix* Ke1=new ElementMatrix(pentabase->nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
-	ElementMatrix* Ke2=new ElementMatrix(this->nodes     ,NUMVERTICES,this->parameters,PattynApproximationEnum);
+	ElementMatrix* Ke2=new ElementMatrix(this->nodes     ,NUMVERTICES,this->parameters,HOApproximationEnum);
 	ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
 	delete Ke1; delete Ke2;
@@ -6058,9 +6058,9 @@
 
 		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
-		GetBMacAyealPattyn(&B[0][0], &xyz_list[0][0], gauss);
+		GetBMacAyealHO(&B[0][0], &xyz_list[0][0], gauss);
 		tria->GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss_tria);
 
-		this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
-		this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
+		this->GetStrainRate3dHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
+		this->GetStrainRate3dHO(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
 		material->GetViscosity3d(&viscosity, &epsilon[0]);
 		material->GetViscosity3d(&oldviscosity, &oldepsilon[0]);
@@ -6090,6 +6090,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreateKMatrixCouplingMacAyealPattynFriction{{{*/
-ElementMatrix* Penta::CreateKMatrixCouplingMacAyealPattynFriction(void){
+/*FUNCTION Penta::CreateKMatrixCouplingMacAyealHOFriction{{{*/
+ElementMatrix* Penta::CreateKMatrixCouplingMacAyealHOFriction(void){
 
 	/*Constants*/
@@ -6116,5 +6116,5 @@
 	if(IsFloating() || !IsOnBed()) return NULL;
 	ElementMatrix* Ke1=new ElementMatrix(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
-	ElementMatrix* Ke2=new ElementMatrix(nodes,NUMVERTICES,this->parameters,PattynApproximationEnum);
+	ElementMatrix* Ke2=new ElementMatrix(nodes,NUMVERTICES,this->parameters,HOApproximationEnum);
 	ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
 	delete Ke1; delete Ke2;
@@ -6149,5 +6149,5 @@
 
 		GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
-		GetBPattynFriction(&L[0][0],gauss);
+		GetBHOFriction(&L[0][0],gauss);
 
 		DL_scalar=alpha2*gauss->weight*Jdet2d;
@@ -6175,10 +6175,10 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreateKMatrixCouplingMacAyealStokes{{{*/
-ElementMatrix* Penta::CreateKMatrixCouplingMacAyealStokes(void){
+/*FUNCTION Penta::CreateKMatrixCouplingMacAyealFS{{{*/
+ElementMatrix* Penta::CreateKMatrixCouplingMacAyealFS(void){
 
 	/*compute all stiffness matrices for this element*/
-	ElementMatrix* Ke1=CreateKMatrixCouplingMacAyealStokesViscous();
-	ElementMatrix* Ke2=CreateKMatrixCouplingMacAyealStokesFriction();
+	ElementMatrix* Ke1=CreateKMatrixCouplingMacAyealFSViscous();
+	ElementMatrix* Ke2=CreateKMatrixCouplingMacAyealFSFriction();
 	ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
 
@@ -6189,6 +6189,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreateKMatrixCouplingMacAyealStokesViscous{{{*/
-ElementMatrix* Penta::CreateKMatrixCouplingMacAyealStokesViscous(void){
+/*FUNCTION Penta::CreateKMatrixCouplingMacAyealFSViscous{{{*/
+ElementMatrix* Penta::CreateKMatrixCouplingMacAyealFSViscous(void){
 
 	/*Constants*/
@@ -6201,5 +6201,5 @@
 	int         i,j;
 	IssmDouble Jdet;
-	IssmDouble viscosity,stokesreconditioning; //viscosity
+	IssmDouble viscosity,FSreconditioning; //viscosity
 	IssmDouble epsilon[6]; /* epsilon=[exx,eyy,exy,exz,eyz];*/
 	IssmDouble xyz_list[NUMVERTICES][3];
@@ -6220,5 +6220,5 @@
 	int         cs_list[numnodes];
 
-	/*Find penta on bed as stokes must be coupled to the dofs on the bed: */
+	/*Find penta on bed as FS must be coupled to the dofs on the bed: */
 	Penta* pentabase=GetBasalElement();
 	Tria* tria=pentabase->SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
@@ -6234,5 +6234,5 @@
 	/*Initialize Element matrix and return if necessary*/
 	ElementMatrix* Ke1=new ElementMatrix(pentabase->nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
-	ElementMatrix* Ke2=new ElementMatrix(this->nodes     ,NUMVERTICES,this->parameters,StokesApproximationEnum);
+	ElementMatrix* Ke2=new ElementMatrix(this->nodes     ,NUMVERTICES,this->parameters,FSApproximationEnum);
 	ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
 	delete Ke1; delete Ke2;
@@ -6240,5 +6240,5 @@
 	/* Get node coordinates and dof list: */
 	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
+	parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
 	Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
 	Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
@@ -6254,15 +6254,15 @@
 
 		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
-		GetBMacAyealStokes(&B[0][0], &xyz_list[0][0], gauss);
-		tria->GetBprimeMacAyealStokes(&Bprime[0][0], &xyz_list[0][0], gauss_tria);
-		tria->GetBMacAyealStokes(&B2[0][0], &xyz_list[0][0], gauss_tria);
-		GetBprimeMacAyealStokes(&Bprime2[0][0], &xyz_list[0][0], gauss);
+		GetBMacAyealFS(&B[0][0], &xyz_list[0][0], gauss);
+		tria->GetBprimeMacAyealFS(&Bprime[0][0], &xyz_list[0][0], gauss_tria);
+		tria->GetBMacAyealFS(&B2[0][0], &xyz_list[0][0], gauss_tria);
+		GetBprimeMacAyealFS(&Bprime2[0][0], &xyz_list[0][0], gauss);
 
 		this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
-		material->GetViscosity3dStokes(&viscosity, &epsilon[0]);
+		material->GetViscosity3dFS(&viscosity, &epsilon[0]);
 
 		D_scalar=2*viscosity*gauss->weight*Jdet;
 		for (i=0;i<3;i++) D[i][i]=D_scalar;
-		D[3][3]=-gauss->weight*Jdet*stokesreconditioning;
+		D[3][3]=-gauss->weight*Jdet*FSreconditioning;
 		for (i=0;i<3;i++) D2[i][i]=D_scalar;
 
@@ -6293,6 +6293,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreateKMatrixCouplingMacAyealStokesFriction {{{*/
-ElementMatrix* Penta::CreateKMatrixCouplingMacAyealStokesFriction(void){
+/*FUNCTION Penta::CreateKMatrixCouplingMacAyealFSFriction {{{*/
+ElementMatrix* Penta::CreateKMatrixCouplingMacAyealFSFriction(void){
 
 	/*Constants*/
@@ -6307,5 +6307,5 @@
 	int        i,j;
 	int        analysis_type,approximation;
-	IssmDouble stokesreconditioning;
+	IssmDouble FSreconditioning;
 	IssmDouble viscosity,alpha2_gauss,Jdet2d;
 	IssmDouble bed_normal[3];
@@ -6313,10 +6313,10 @@
 	IssmDouble xyz_list[NUMVERTICES][3];
 	IssmDouble xyz_list_tria[NUMVERTICES2D][3];
-	IssmDouble LMacAyealStokes[8][numdof2dm];
-	IssmDouble LprimeMacAyealStokes[8][numdof2d];
-	IssmDouble DLMacAyealStokes[8][8]={0.0};
-	IssmDouble LStokesMacAyeal[4][numdof2d];
-	IssmDouble LprimeStokesMacAyeal[4][numdof2dm];
-	IssmDouble DLStokesMacAyeal[4][4]={0.0};
+	IssmDouble LMacAyealFS[8][numdof2dm];
+	IssmDouble LprimeMacAyealFS[8][numdof2d];
+	IssmDouble DLMacAyealFS[8][8]={0.0};
+	IssmDouble LFSMacAyeal[4][numdof2d];
+	IssmDouble LprimeFSMacAyeal[4][numdof2dm];
+	IssmDouble DLFSMacAyeal[4][4]={0.0};
 	IssmDouble Ke_drag_gaussian[numdof2dm][numdof2d];
 	IssmDouble Ke_drag_gaussian2[numdof2d][numdof2dm];
@@ -6326,9 +6326,9 @@
 	int         cs_list[numnodes];
 
-	/*If on water or not Stokes, skip stiffness: */
+	/*If on water or not FS, skip stiffness: */
 	inputs->GetInputValue(&approximation,ApproximationEnum);
 	if(IsFloating() || !IsOnBed()) return NULL;
 	ElementMatrix* Ke1=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);
-	ElementMatrix* Ke2=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
+	ElementMatrix* Ke2=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
 	ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
 	delete Ke1; delete Ke2;
@@ -6345,5 +6345,5 @@
 	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
-	parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
+	parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
 	Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
 	Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
@@ -6361,37 +6361,37 @@
 
 		GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
-		GetLMacAyealStokes(&LMacAyealStokes[0][0], gauss);
-		GetLprimeMacAyealStokes(&LprimeMacAyealStokes[0][0], &xyz_list[0][0], gauss);
-		GetLStokesMacAyeal(&LStokesMacAyeal[0][0], gauss);
-		GetLprimeStokesMacAyeal(&LprimeStokesMacAyeal[0][0], &xyz_list[0][0], gauss);
+		GetLMacAyealFS(&LMacAyealFS[0][0], gauss);
+		GetLprimeMacAyealFS(&LprimeMacAyealFS[0][0], &xyz_list[0][0], gauss);
+		GetLFSMacAyeal(&LFSMacAyeal[0][0], gauss);
+		GetLprimeFSMacAyeal(&LprimeFSMacAyeal[0][0], &xyz_list[0][0], gauss);
 
 		this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
-		material->GetViscosity3dStokes(&viscosity,&epsilon[0]);
+		material->GetViscosity3dFS(&viscosity,&epsilon[0]);
 
 		BedNormal(&bed_normal[0],xyz_list_tria);
 		friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum);
 
-		DLMacAyealStokes[0][0]=alpha2_gauss*gauss->weight*Jdet2d;
-		DLMacAyealStokes[1][1]=alpha2_gauss*gauss->weight*Jdet2d;
-		DLMacAyealStokes[2][2]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[0]*bed_normal[2];
-		DLMacAyealStokes[3][3]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[1]*bed_normal[2];
-		DLMacAyealStokes[4][4]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[0];
-		DLMacAyealStokes[5][5]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[1];
-		DLMacAyealStokes[6][6]=stokesreconditioning*gauss->weight*Jdet2d*bed_normal[0];
-		DLMacAyealStokes[7][7]=stokesreconditioning*gauss->weight*Jdet2d*bed_normal[1];
-
-		DLStokesMacAyeal[0][0]=alpha2_gauss*gauss->weight*Jdet2d;
-		DLStokesMacAyeal[1][1]=alpha2_gauss*gauss->weight*Jdet2d;
-		DLStokesMacAyeal[2][2]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[0]*bed_normal[2];
-		DLStokesMacAyeal[3][3]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[1]*bed_normal[2];
-
-		TripleMultiply( &LMacAyealStokes[0][0],8,numdof2dm,1,
-					&DLMacAyealStokes[0][0],8,8,0,
-					&LprimeMacAyealStokes[0][0],8,numdof2d,0,
+		DLMacAyealFS[0][0]=alpha2_gauss*gauss->weight*Jdet2d;
+		DLMacAyealFS[1][1]=alpha2_gauss*gauss->weight*Jdet2d;
+		DLMacAyealFS[2][2]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[0]*bed_normal[2];
+		DLMacAyealFS[3][3]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[1]*bed_normal[2];
+		DLMacAyealFS[4][4]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[0];
+		DLMacAyealFS[5][5]=-2*viscosity*gauss->weight*Jdet2d*bed_normal[1];
+		DLMacAyealFS[6][6]=FSreconditioning*gauss->weight*Jdet2d*bed_normal[0];
+		DLMacAyealFS[7][7]=FSreconditioning*gauss->weight*Jdet2d*bed_normal[1];
+
+		DLFSMacAyeal[0][0]=alpha2_gauss*gauss->weight*Jdet2d;
+		DLFSMacAyeal[1][1]=alpha2_gauss*gauss->weight*Jdet2d;
+		DLFSMacAyeal[2][2]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[0]*bed_normal[2];
+		DLFSMacAyeal[3][3]=-alpha2_gauss*gauss->weight*Jdet2d*bed_normal[1]*bed_normal[2];
+
+		TripleMultiply( &LMacAyealFS[0][0],8,numdof2dm,1,
+					&DLMacAyealFS[0][0],8,8,0,
+					&LprimeMacAyealFS[0][0],8,numdof2d,0,
 					&Ke_drag_gaussian[0][0],0);
 
-		TripleMultiply( &LStokesMacAyeal[0][0],4,numdof2d,1,
-					&DLStokesMacAyeal[0][0],4,4,0,
-					&LprimeStokesMacAyeal[0][0],4,numdof2dm,0,
+		TripleMultiply( &LFSMacAyeal[0][0],4,numdof2d,1,
+					&DLFSMacAyeal[0][0],4,4,0,
+					&LprimeFSMacAyeal[0][0],4,numdof2dm,0,
 					&Ke_drag_gaussian2[0][0],0);
 
@@ -6409,6 +6409,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreateKMatrixCouplingPattynStokes{{{*/
-ElementMatrix* Penta::CreateKMatrixCouplingPattynStokes(void){
+/*FUNCTION Penta::CreateKMatrixCouplingHOFS{{{*/
+ElementMatrix* Penta::CreateKMatrixCouplingHOFS(void){
 
 	/*Constants*/
@@ -6432,11 +6432,11 @@
 
 	/*compute all stiffness matrices for this element*/
-	ElementMatrix* Ke1=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,PattynApproximationEnum);
-	ElementMatrix* Ke2=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
+	ElementMatrix* Ke1=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,HOApproximationEnum);
+	ElementMatrix* Ke2=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
 	ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
 	delete Ke1;
 	delete Ke2;
-	Ke1=CreateKMatrixDiagnosticPattyn(); TransformInvStiffnessMatrixCoord(Ke1,this->nodes,NUMVERTICES,XYEnum);
-	Ke2=CreateKMatrixDiagnosticStokes(); TransformInvStiffnessMatrixCoord(Ke2,this->nodes,NUMVERTICES,XYZPEnum);
+	Ke1=CreateKMatrixDiagnosticHO(); TransformInvStiffnessMatrixCoord(Ke1,this->nodes,NUMVERTICES,XYEnum);
+	Ke2=CreateKMatrixDiagnosticFS(); TransformInvStiffnessMatrixCoord(Ke2,this->nodes,NUMVERTICES,XYZPEnum);
 
 	for(i=0;i<numdofs;i++) for(j=0;j<NUMVERTICES;j++){
@@ -6469,18 +6469,18 @@
 		case L1L2ApproximationEnum:
 			return CreateKMatrixDiagnosticL1L2();
-		case PattynApproximationEnum:
-			return CreateKMatrixDiagnosticPattyn();
-		case StokesApproximationEnum:
-			return CreateKMatrixDiagnosticStokes();
-		case HutterApproximationEnum:
+		case HOApproximationEnum:
+			return CreateKMatrixDiagnosticHO();
+		case FSApproximationEnum:
+			return CreateKMatrixDiagnosticFS();
+		case SIAApproximationEnum:
 			return NULL;
 		case NoneApproximationEnum:
 			return NULL;
-		case MacAyealPattynApproximationEnum:
-			return CreateKMatrixDiagnosticMacAyealPattyn();
-		case MacAyealStokesApproximationEnum:
-			return CreateKMatrixDiagnosticMacAyealStokes();
-		case PattynStokesApproximationEnum:
-			return CreateKMatrixDiagnosticPattynStokes();
+		case MacAyealHOApproximationEnum:
+			return CreateKMatrixDiagnosticMacAyealHO();
+		case MacAyealFSApproximationEnum:
+			return CreateKMatrixDiagnosticMacAyealFS();
+		case HOFSApproximationEnum:
+			return CreateKMatrixDiagnosticHOFS();
 		default:
 			_error_("Approximation " << EnumToStringx(approximation) << " not supported yet");
@@ -6488,6 +6488,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreateKMatrixDiagnosticHutter{{{*/
-ElementMatrix* Penta::CreateKMatrixDiagnosticHutter(void){
+/*FUNCTION Penta::CreateKMatrixDiagnosticSIA{{{*/
+ElementMatrix* Penta::CreateKMatrixDiagnosticSIA(void){
 
 	/*Constants*/
@@ -6610,5 +6610,5 @@
 	IssmDouble  viscosity , oldviscosity, newviscosity, viscosity_overshoot;
 	IssmDouble  epsilon[5],oldepsilon[5];       /* epsilon=[exx,eyy,exy,exz,eyz];*/
-	IssmDouble  epsilons[6];                    //6 for stokes
+	IssmDouble  epsilons[6];                    //6 for FS
 	IssmDouble  xyz_list[NUMVERTICES][3];
 	IssmDouble  B[3][numdof2d];
@@ -6622,5 +6622,5 @@
 	GaussTria  *gauss_tria = NULL;
 
-	/*Find penta on bed as this is a macayeal elements: */
+	/*Find penta on bed as this is a SSA elements: */
 	pentabase=GetBasalElement();
 	tria=pentabase->SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
@@ -6651,7 +6651,7 @@
 		tria->GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss_tria);
 
-		if(approximation==MacAyealPattynApproximationEnum){
-			this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
-			this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
+		if(approximation==MacAyealHOApproximationEnum){
+			this->GetStrainRate3dHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
+			this->GetStrainRate3dHO(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
 			material->GetViscosity3d(&viscosity, &epsilon[0]);
 			material->GetViscosity3d(&oldviscosity, &oldepsilon[0]);
@@ -6659,7 +6659,7 @@
 			newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
 		}
-		else if (approximation==MacAyealStokesApproximationEnum){
+		else if (approximation==MacAyealFSApproximationEnum){
 			this->GetStrainRate3d(&epsilons[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
-			material->GetViscosity3dStokes(&newviscosity,&epsilons[0]);
+			material->GetViscosity3dFS(&newviscosity,&epsilons[0]);
 		}
 		else _error_("approximation " << approximation << " (" << EnumToStringx(approximation) << ") not supported yet");
@@ -6704,11 +6704,11 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreateKMatrixDiagnosticMacAyealPattyn{{{*/
-ElementMatrix* Penta::CreateKMatrixDiagnosticMacAyealPattyn(void){
+/*FUNCTION Penta::CreateKMatrixDiagnosticMacAyealHO{{{*/
+ElementMatrix* Penta::CreateKMatrixDiagnosticMacAyealHO(void){
 
 	/*compute all stiffness matrices for this element*/
 	ElementMatrix* Ke1=CreateKMatrixDiagnosticMacAyeal3d();
-	ElementMatrix* Ke2=CreateKMatrixDiagnosticPattyn();
-	ElementMatrix* Ke3=CreateKMatrixCouplingMacAyealPattyn();
+	ElementMatrix* Ke2=CreateKMatrixDiagnosticHO();
+	ElementMatrix* Ke3=CreateKMatrixCouplingMacAyealHO();
 	ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
 
@@ -6720,11 +6720,11 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreateKMatrixDiagnosticMacAyealStokes{{{*/
-ElementMatrix* Penta::CreateKMatrixDiagnosticMacAyealStokes(void){
+/*FUNCTION Penta::CreateKMatrixDiagnosticMacAyealFS{{{*/
+ElementMatrix* Penta::CreateKMatrixDiagnosticMacAyealFS(void){
 
 	/*compute all stiffness matrices for this element*/
 	ElementMatrix* Ke1=CreateKMatrixDiagnosticMacAyeal3d();
-	ElementMatrix* Ke2=CreateKMatrixDiagnosticStokes();
-	ElementMatrix* Ke3=CreateKMatrixCouplingMacAyealStokes();
+	ElementMatrix* Ke2=CreateKMatrixDiagnosticFS();
+	ElementMatrix* Ke3=CreateKMatrixCouplingMacAyealFS();
 	ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
 
@@ -6770,5 +6770,5 @@
 	GaussTria  *gauss_tria = NULL;
 
-	/*Find penta on bed as this is a macayeal elements: */
+	/*Find penta on bed as this is a SSA elements: */
 	pentabase=GetBasalElement();
 	tria=pentabase->SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
@@ -6835,10 +6835,10 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreateKMatrixDiagnosticPattyn{{{*/
-ElementMatrix* Penta::CreateKMatrixDiagnosticPattyn(void){
+/*FUNCTION Penta::CreateKMatrixDiagnosticHO{{{*/
+ElementMatrix* Penta::CreateKMatrixDiagnosticHO(void){
 
 	/*compute all stiffness matrices for this element*/
-	ElementMatrix* Ke1=CreateKMatrixDiagnosticPattynViscous();
-	ElementMatrix* Ke2=CreateKMatrixDiagnosticPattynFriction();
+	ElementMatrix* Ke1=CreateKMatrixDiagnosticHOViscous();
+	ElementMatrix* Ke2=CreateKMatrixDiagnosticHOFriction();
 	ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
 
@@ -6850,6 +6850,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreateKMatrixDiagnosticPattynViscous{{{*/
-ElementMatrix* Penta::CreateKMatrixDiagnosticPattynViscous(void){
+/*FUNCTION Penta::CreateKMatrixDiagnosticHOViscous{{{*/
+ElementMatrix* Penta::CreateKMatrixDiagnosticHOViscous(void){
 
 	/*Constants*/
@@ -6871,5 +6871,5 @@
 
 	/*Initialize Element matrix*/
-	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,PattynApproximationEnum);
+	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,HOApproximationEnum);
 
 	/*Retrieve all inputs and parameters*/
@@ -6889,9 +6889,9 @@
 
 		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
-		GetBPattyn(&B[0][0], &xyz_list[0][0], gauss);
-		GetBprimePattyn(&Bprime[0][0], &xyz_list[0][0], gauss);
-
-		this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
-		this->GetStrainRate3dPattyn(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
+		GetBHO(&B[0][0], &xyz_list[0][0], gauss);
+		GetBprimeHO(&Bprime[0][0], &xyz_list[0][0], gauss);
+
+		this->GetStrainRate3dHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
+		this->GetStrainRate3dHO(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
 		material->GetViscosity3d(&viscosity, &epsilon[0]);
 		material->GetViscosity3d(&oldviscosity, &oldepsilon[0]);
@@ -6915,6 +6915,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreateKMatrixDiagnosticPattynFriction{{{*/
-ElementMatrix* Penta::CreateKMatrixDiagnosticPattynFriction(void){
+/*FUNCTION Penta::CreateKMatrixDiagnosticHOFriction{{{*/
+ElementMatrix* Penta::CreateKMatrixDiagnosticHOFriction(void){
 
 	/*Constants*/
@@ -6937,5 +6937,5 @@
 	if(IsFloating() || !IsOnBed()) return NULL;
 
-	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,PattynApproximationEnum);
+	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,HOApproximationEnum);
 
 	/*Retrieve all inputs and parameters*/
@@ -6962,5 +6962,5 @@
 
 		GetTriaJacobianDeterminant(&Jdet, &xyz_list_tria[0][0],gauss);
-		GetBPattynFriction(&L[0][0],gauss);
+		GetBHOFriction(&L[0][0],gauss);
 
 		friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum); 
@@ -6985,11 +6985,11 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreateKMatrixDiagnosticPattynStokes{{{*/
-ElementMatrix* Penta::CreateKMatrixDiagnosticPattynStokes(void){
+/*FUNCTION Penta::CreateKMatrixDiagnosticHOFS{{{*/
+ElementMatrix* Penta::CreateKMatrixDiagnosticHOFS(void){
 
 	/*compute all stiffness matrices for this element*/
-	ElementMatrix* Ke1=CreateKMatrixDiagnosticPattyn();
-	ElementMatrix* Ke2=CreateKMatrixDiagnosticStokes();
-	ElementMatrix* Ke3=CreateKMatrixCouplingPattynStokes();
+	ElementMatrix* Ke1=CreateKMatrixDiagnosticHO();
+	ElementMatrix* Ke2=CreateKMatrixDiagnosticFS();
+	ElementMatrix* Ke3=CreateKMatrixCouplingHOFS();
 	ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
 
@@ -7001,28 +7001,28 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreateKMatrixDiagnosticStokes{{{*/
-ElementMatrix* Penta::CreateKMatrixDiagnosticStokes(void){
-
-	int fe_stokes;
+/*FUNCTION Penta::CreateKMatrixDiagnosticFS{{{*/
+ElementMatrix* Penta::CreateKMatrixDiagnosticFS(void){
+
+	int fe_FS;
 	ElementMatrix* Ke1;
 	ElementMatrix* Ke2;
 	ElementMatrix* Ke;
-	parameters->FindParam(&fe_stokes,FlowequationFeStokesEnum);
-
-	switch(fe_stokes){
+	parameters->FindParam(&fe_FS,FlowequationFeFSEnum);
+
+	switch(fe_FS){
 		case 0:
 			/*compute all stiffness matrices for this element*/
-			Ke1=CreateKMatrixDiagnosticStokesViscous();
-			Ke2=CreateKMatrixDiagnosticStokesFriction();
+			Ke1=CreateKMatrixDiagnosticFSViscous();
+			Ke2=CreateKMatrixDiagnosticFSFriction();
 			Ke =new ElementMatrix(Ke1,Ke2);
 			break;
 		case 1:
 			/*compute all stiffness matrices for this element*/
-			Ke1=CreateKMatrixDiagnosticStokesGLSViscous();
-			Ke2=CreateKMatrixDiagnosticStokesFriction();
+			Ke1=CreateKMatrixDiagnosticFSGLSViscous();
+			Ke2=CreateKMatrixDiagnosticFSFriction();
 			Ke =new ElementMatrix(Ke1,Ke2);
 			break;
 		default:
-			_error_("Finite element" << fe_stokes << " not supported yet");
+			_error_("Finite element" << fe_FS << " not supported yet");
 	}
 
@@ -7034,10 +7034,10 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreateKMatrixDiagnosticStokesViscous {{{*/
-ElementMatrix* Penta::CreateKMatrixDiagnosticStokesViscous(void){
+/*FUNCTION Penta::CreateKMatrixDiagnosticFSViscous {{{*/
+ElementMatrix* Penta::CreateKMatrixDiagnosticFSViscous(void){
 
 	/*Intermediaries */
 	int        i,approximation;
-	IssmDouble Jdet,viscosity,stokesreconditioning;
+	IssmDouble Jdet,viscosity,FSreconditioning;
 	IssmDouble xyz_list[NUMVERTICES][3];
 	IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
@@ -7049,12 +7049,12 @@
 	GaussPenta *gauss=NULL;
 
-	/*If on water or not Stokes, skip stiffness: */
+	/*If on water or not FS, skip stiffness: */
 	inputs->GetInputValue(&approximation,ApproximationEnum);
-	if(approximation!=StokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum && approximation!=PattynStokesApproximationEnum) return NULL;
-	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
+	if(approximation!=FSApproximationEnum && approximation!=MacAyealFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
+	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
 
 	/*Retrieve all inputs and parameters*/
 	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
+	parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
 	Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
 	Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
@@ -7068,13 +7068,13 @@
 
 		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
-		GetBStokes(&B[0][0],&xyz_list[0][0],gauss); 
-		GetBprimeStokes(&B_prime[0][0],&xyz_list[0][0],gauss); 
+		GetBFS(&B[0][0],&xyz_list[0][0],gauss); 
+		GetBprimeFS(&B_prime[0][0],&xyz_list[0][0],gauss); 
 
 		this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
-		material->GetViscosity3dStokes(&viscosity,&epsilon[0]);
+		material->GetViscosity3dFS(&viscosity,&epsilon[0]);
 
 		D_scalar=gauss->weight*Jdet;
 		for (i=0;i<6;i++) D[i][i]=D_scalar*2*viscosity;
-		for (i=6;i<8;i++) D[i][i]=-D_scalar*stokesreconditioning;
+		for (i=6;i<8;i++) D[i][i]=-D_scalar*FSreconditioning;
 
 		TripleMultiply( &B[0][0],8,27,1,
@@ -7085,5 +7085,5 @@
 
 	/*Condensation*/
-	ReduceMatrixStokes(Ke->values, &Ke_temp[0][0]);
+	ReduceMatrixFS(Ke->values, &Ke_temp[0][0]);
 	//Ke->Echo();
 	//_error_("stop");
@@ -7097,6 +7097,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreateKMatrixDiagnosticStokesGLSViscous {{{*/
-ElementMatrix* Penta::CreateKMatrixDiagnosticStokesGLSViscous(void){
+/*FUNCTION Penta::CreateKMatrixDiagnosticFSGLSViscous {{{*/
+ElementMatrix* Penta::CreateKMatrixDiagnosticFSGLSViscous(void){
 
 	int        numdof  = NUMVERTICES*NDOF4;
@@ -7104,5 +7104,5 @@
 	/*Intermediaries */
 	int        i,j,approximation;
-	IssmDouble Jdet,viscosity,stokesreconditioning,diameter,rigidity;
+	IssmDouble Jdet,viscosity,FSreconditioning,diameter,rigidity;
 	IssmDouble xyz_list[NUMVERTICES][3];
 	IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
@@ -7128,12 +7128,12 @@
 	int c=3; //index of pressure
 
-	/*If on water or not Stokes, skip stiffness: */
+	/*If on water or not FS, skip stiffness: */
 	inputs->GetInputValue(&approximation,ApproximationEnum);
-	if(approximation!=StokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum && approximation!=PattynStokesApproximationEnum) return NULL;
-	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
+	if(approximation!=FSApproximationEnum && approximation!=MacAyealFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
+	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
 
 	/*Retrieve all inputs and parameters*/
 	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
+	parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
 	Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
 	Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
@@ -7165,13 +7165,13 @@
 
 		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
-		GetBStokesGLS(&B[0][0],&xyz_list[0][0],gauss); 
-		GetBprimeStokesGLS(&B_prime[0][0],&xyz_list[0][0],gauss); 
+		GetBFSGLS(&B[0][0],&xyz_list[0][0],gauss); 
+		GetBprimeFSGLS(&B_prime[0][0],&xyz_list[0][0],gauss); 
 
 		this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
-		material->GetViscosity3dStokes(&viscosity,&epsilon[0]);
+		material->GetViscosity3dFS(&viscosity,&epsilon[0]);
 
 		D_scalar=gauss->weight*Jdet;
 		for (i=0;i<6;i++) D[i][i]=D_scalar*2.*2.*viscosity;
-		for (i=6;i<8;i++) D[i][i]=-D_scalar*stokesreconditioning;
+		for (i=6;i<8;i++) D[i][i]=-D_scalar*FSreconditioning;
 
 		TripleMultiply( &B[0][0],8,24,1,
@@ -7235,6 +7235,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreateKMatrixDiagnosticStokesFriction{{{*/
-ElementMatrix* Penta::CreateKMatrixDiagnosticStokesFriction(void){
+/*FUNCTION Penta::CreateKMatrixDiagnosticFSFriction{{{*/
+ElementMatrix* Penta::CreateKMatrixDiagnosticFSFriction(void){
 
 	/*Constants*/
@@ -7246,23 +7246,23 @@
 	int        analysis_type,approximation;
 	IssmDouble alpha2,Jdet2d;
-	IssmDouble stokesreconditioning,viscosity;
+	IssmDouble FSreconditioning,viscosity;
 	IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
 	IssmDouble xyz_list[NUMVERTICES][3];
 	IssmDouble xyz_list_tria[NUMVERTICES2D][3];
-	IssmDouble LStokes[2][numdof2d];
-	IssmDouble DLStokes[2][2]={0.0};
+	IssmDouble LFS[2][numdof2d];
+	IssmDouble DLFS[2][2]={0.0};
 	IssmDouble Ke_drag_gaussian[numdof2d][numdof2d];
 	Friction*  friction=NULL;
 	GaussPenta *gauss=NULL;
 
-	/*If on water or not Stokes, skip stiffness: */
+	/*If on water or not FS, skip stiffness: */
 	inputs->GetInputValue(&approximation,ApproximationEnum);
-	if(IsFloating() || !IsOnBed() || (approximation!=StokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum &&  approximation!=PattynStokesApproximationEnum)) return NULL;
-	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
+	if(IsFloating() || !IsOnBed() || (approximation!=FSApproximationEnum && approximation!=MacAyealFSApproximationEnum &&  approximation!=HOFSApproximationEnum)) return NULL;
+	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
 
 	/*Retrieve all inputs and parameters*/
 	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
-	parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
+	parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
 	Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
 	Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
@@ -7280,17 +7280,17 @@
 
 		GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0],gauss);
-		GetLStokes(&LStokes[0][0], gauss);
+		GetLFS(&LFS[0][0], gauss);
 
 		this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
-		material->GetViscosity3dStokes(&viscosity,&epsilon[0]);
+		material->GetViscosity3dFS(&viscosity,&epsilon[0]);
 
 		friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);
 
-		DLStokes[0][0] = +alpha2*gauss->weight*Jdet2d; //taub_x = -alpha2 vx
-		DLStokes[1][1] = +alpha2*gauss->weight*Jdet2d; //taub_y = -alpha2 vy
-
-		TripleMultiply( &LStokes[0][0],2,numdof2d,1,
-					&DLStokes[0][0],2,2,0,
-					&LStokes[0][0],2,numdof2d,0,
+		DLFS[0][0] = +alpha2*gauss->weight*Jdet2d; //taub_x = -alpha2 vx
+		DLFS[1][1] = +alpha2*gauss->weight*Jdet2d; //taub_y = -alpha2 vy
+
+		TripleMultiply( &LFS[0][0],2,numdof2d,1,
+					&DLFS[0][0],2,2,0,
+					&LFS[0][0],2,numdof2d,0,
 					&Ke_drag_gaussian[0][0],0);
 
@@ -7412,10 +7412,10 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreatePVectorCouplingMacAyealStokes {{{*/
-ElementVector* Penta::CreatePVectorCouplingMacAyealStokes(void){
+/*FUNCTION Penta::CreatePVectorCouplingMacAyealFS {{{*/
+ElementVector* Penta::CreatePVectorCouplingMacAyealFS(void){
 
 	/*compute all load vectors for this element*/
-	ElementVector* pe1=CreatePVectorCouplingMacAyealStokesViscous();
-	ElementVector* pe2=CreatePVectorCouplingMacAyealStokesFriction();
+	ElementVector* pe1=CreatePVectorCouplingMacAyealFSViscous();
+	ElementVector* pe2=CreatePVectorCouplingMacAyealFSFriction();
 	ElementVector* pe =new ElementVector(pe1,pe2);
 
@@ -7426,6 +7426,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreatePVectorCouplingMacAyealStokesViscous {{{*/
-ElementVector* Penta::CreatePVectorCouplingMacAyealStokesViscous(void){
+/*FUNCTION Penta::CreatePVectorCouplingMacAyealFSViscous {{{*/
+ElementVector* Penta::CreatePVectorCouplingMacAyealFSViscous(void){
 
 	/*Constants*/
@@ -7436,5 +7436,5 @@
 	int         approximation;
 	IssmDouble  viscosity,Jdet;
-	IssmDouble  stokesreconditioning;
+	IssmDouble  FSreconditioning;
 	IssmDouble  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
 	IssmDouble  dw[3];
@@ -7446,14 +7446,14 @@
 	/*Initialize Element vector and return if necessary*/
 	inputs->GetInputValue(&approximation,ApproximationEnum);
-	if(approximation!=MacAyealStokesApproximationEnum) return NULL;
-	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
+	if(approximation!=MacAyealFSApproximationEnum) return NULL;
+	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
 
 	/*Retrieve all inputs and parameters*/
 	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	this->parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
+	this->parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
 	Input* vx_input=inputs->GetInput(VxEnum);               _assert_(vx_input);
 	Input* vy_input=inputs->GetInput(VyEnum);               _assert_(vy_input);
 	Input* vz_input=inputs->GetInput(VzEnum);               _assert_(vz_input);
-	Input* vzmacayeal_input=inputs->GetInput(VzMacAyealEnum);   _assert_(vzmacayeal_input);
+	Input* vzSSA_input=inputs->GetInput(VzMacAyealEnum);   _assert_(vzSSA_input);
 
 	/* Start  looping on the number of gaussian points: */
@@ -7467,8 +7467,8 @@
 		GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
 
-		vzmacayeal_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
+		vzSSA_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
 
 		this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
-		material->GetViscosity3dStokes(&viscosity,&epsilon[0]);
+		material->GetViscosity3dFS(&viscosity,&epsilon[0]);
 
 		for(i=0;i<NUMVERTICES;i++){
@@ -7476,5 +7476,5 @@
 			pe->values[i*NDOF4+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dbasis[2][i];
 			pe->values[i*NDOF4+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dbasis[0][i]+dw[1]*dbasis[1][i]+2*dw[2]*dbasis[2][i]);
-			pe->values[i*NDOF4+3]+=Jdet*gauss->weight*stokesreconditioning*dw[2]*basis[i];
+			pe->values[i*NDOF4+3]+=Jdet*gauss->weight*FSreconditioning*dw[2]*basis[i];
 		}
 	}
@@ -7488,6 +7488,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreatePVectorCouplingMacAyealStokesFriction{{{*/
-ElementVector* Penta::CreatePVectorCouplingMacAyealStokesFriction(void){
+/*FUNCTION Penta::CreatePVectorCouplingMacAyealFSFriction{{{*/
+ElementVector* Penta::CreatePVectorCouplingMacAyealFSFriction(void){
 
 	/*Constants*/
@@ -7498,5 +7498,5 @@
 	int         approximation,analysis_type;
 	IssmDouble  Jdet,Jdet2d;
-	IssmDouble  stokesreconditioning;
+	IssmDouble  FSreconditioning;
 	IssmDouble	bed_normal[3];
 	IssmDouble  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
@@ -7513,15 +7513,15 @@
 	if(!IsOnBed() || IsFloating()) return NULL;
 	inputs->GetInputValue(&approximation,ApproximationEnum);
-	if(approximation!=MacAyealStokesApproximationEnum) return NULL;
-	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
+	if(approximation!=MacAyealFSApproximationEnum) return NULL;
+	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
 
 	/*Retrieve all inputs and parameters*/
 	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
-	this->parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
+	this->parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
 	Input* vx_input=inputs->GetInput(VxEnum);               _assert_(vx_input);
 	Input* vy_input=inputs->GetInput(VyEnum);               _assert_(vy_input);
 	Input* vz_input=inputs->GetInput(VzEnum);               _assert_(vz_input);
-	Input* vzmacayeal_input=inputs->GetInput(VzMacAyealEnum);   _assert_(vzmacayeal_input);
+	Input* vzSSA_input=inputs->GetInput(VzMacAyealEnum);   _assert_(vzSSA_input);
 
 	for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
@@ -7539,10 +7539,10 @@
 		GetNodalFunctionsP1(basis, gauss);
 
-		vzmacayeal_input->GetInputValue(&w, gauss);
-		vzmacayeal_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
+		vzSSA_input->GetInputValue(&w, gauss);
+		vzSSA_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
 
 		BedNormal(&bed_normal[0],xyz_list_tria);
 		this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
-		material->GetViscosity3dStokes(&viscosity,&epsilon[0]);
+		material->GetViscosity3dFS(&viscosity,&epsilon[0]);
 		friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum);
 
@@ -7563,10 +7563,10 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreatePVectorCouplingPattynStokes {{{*/
-ElementVector* Penta::CreatePVectorCouplingPattynStokes(void){
+/*FUNCTION Penta::CreatePVectorCouplingHOFS {{{*/
+ElementVector* Penta::CreatePVectorCouplingHOFS(void){
 
 	/*compute all load vectors for this element*/
-	ElementVector* pe1=CreatePVectorCouplingPattynStokesViscous();
-	ElementVector* pe2=CreatePVectorCouplingPattynStokesFriction();
+	ElementVector* pe1=CreatePVectorCouplingHOFSViscous();
+	ElementVector* pe2=CreatePVectorCouplingHOFSFriction();
 	ElementVector* pe =new ElementVector(pe1,pe2);
 
@@ -7577,6 +7577,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreatePVectorCouplingPattynStokesViscous {{{*/
-ElementVector* Penta::CreatePVectorCouplingPattynStokesViscous(void){
+/*FUNCTION Penta::CreatePVectorCouplingHOFSViscous {{{*/
+ElementVector* Penta::CreatePVectorCouplingHOFSViscous(void){
 
 	/*Constants*/
@@ -7587,5 +7587,5 @@
 	int         approximation;
 	IssmDouble  viscosity,Jdet;
-	IssmDouble  stokesreconditioning;
+	IssmDouble  FSreconditioning;
 	IssmDouble  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
 	IssmDouble  dw[3];
@@ -7597,14 +7597,14 @@
 	/*Initialize Element vector and return if necessary*/
 	inputs->GetInputValue(&approximation,ApproximationEnum);
-	if(approximation!=PattynStokesApproximationEnum) return NULL;
-	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
+	if(approximation!=HOFSApproximationEnum) return NULL;
+	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
 
 	/*Retrieve all inputs and parameters*/
 	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	this->parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
+	this->parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
 	Input* vx_input=inputs->GetInput(VxEnum);               _assert_(vx_input);
 	Input* vy_input=inputs->GetInput(VyEnum);               _assert_(vy_input);
 	Input* vz_input=inputs->GetInput(VzEnum);               _assert_(vz_input);
-	Input* vzpattyn_input=inputs->GetInput(VzPattynEnum);   _assert_(vzpattyn_input);
+	Input* vzHO_input=inputs->GetInput(VzHOEnum);   _assert_(vzHO_input);
 
 	/* Start  looping on the number of gaussian points: */
@@ -7618,8 +7618,8 @@
 		GetNodalFunctionsP1Derivatives(&dbasis[0][0],&xyz_list[0][0], gauss);
 
-		vzpattyn_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
+		vzHO_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
 
 		this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
-		material->GetViscosity3dStokes(&viscosity,&epsilon[0]);
+		material->GetViscosity3dFS(&viscosity,&epsilon[0]);
 
 		for(i=0;i<NUMVERTICES;i++){
@@ -7627,5 +7627,5 @@
 			pe->values[i*NDOF4+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dbasis[2][i];
 			pe->values[i*NDOF4+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dbasis[0][i]+dw[1]*dbasis[1][i]+2*dw[2]*dbasis[2][i]);
-			pe->values[i*NDOF4+3]+=Jdet*gauss->weight*stokesreconditioning*dw[2]*basis[i];
+			pe->values[i*NDOF4+3]+=Jdet*gauss->weight*FSreconditioning*dw[2]*basis[i];
 		}
 	}
@@ -7639,6 +7639,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreatePVectorCouplingPattynStokesFriction{{{*/
-ElementVector* Penta::CreatePVectorCouplingPattynStokesFriction(void){
+/*FUNCTION Penta::CreatePVectorCouplingHOFSFriction{{{*/
+ElementVector* Penta::CreatePVectorCouplingHOFSFriction(void){
 
 	/*Constants*/
@@ -7649,5 +7649,5 @@
 	int         approximation,analysis_type;
 	IssmDouble  Jdet,Jdet2d;
-	IssmDouble  stokesreconditioning;
+	IssmDouble  FSreconditioning;
 	IssmDouble	bed_normal[3];
 	IssmDouble  epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
@@ -7664,15 +7664,15 @@
 	if(!IsOnBed() || IsFloating()) return NULL;
 	inputs->GetInputValue(&approximation,ApproximationEnum);
-	if(approximation!=PattynStokesApproximationEnum) return NULL;
-	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
+	if(approximation!=HOFSApproximationEnum) return NULL;
+	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
 
 	/*Retrieve all inputs and parameters*/
 	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
-	this->parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
+	this->parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
 	Input* vx_input=inputs->GetInput(VxEnum);               _assert_(vx_input);
 	Input* vy_input=inputs->GetInput(VyEnum);               _assert_(vy_input);
 	Input* vz_input=inputs->GetInput(VzEnum);               _assert_(vz_input);
-	Input* vzpattyn_input=inputs->GetInput(VzPattynEnum);   _assert_(vzpattyn_input);
+	Input* vzHO_input=inputs->GetInput(VzHOEnum);   _assert_(vzHO_input);
 
 	for(i=0;i<NUMVERTICES2D;i++) for(j=0;j<3;j++) xyz_list_tria[i][j]=xyz_list[i][j];
@@ -7690,10 +7690,10 @@
 		GetNodalFunctionsP1(basis, gauss);
 
-		vzpattyn_input->GetInputValue(&w, gauss);
-		vzpattyn_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
+		vzHO_input->GetInputValue(&w, gauss);
+		vzHO_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
 
 		BedNormal(&bed_normal[0],xyz_list_tria);
 		this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
-		material->GetViscosity3dStokes(&viscosity,&epsilon[0]);
+		material->GetViscosity3dFS(&viscosity,&epsilon[0]);
 		friction->GetAlpha2(&alpha2_gauss, gauss,VxEnum,VyEnum,VzEnum);
 
@@ -7723,20 +7723,20 @@
 		case MacAyealApproximationEnum:
 			return CreatePVectorDiagnosticMacAyeal();
-		case PattynApproximationEnum:
-			return CreatePVectorDiagnosticPattyn();
+		case HOApproximationEnum:
+			return CreatePVectorDiagnosticHO();
 		case L1L2ApproximationEnum:
 			return CreatePVectorDiagnosticL1L2();
-		case HutterApproximationEnum:
+		case SIAApproximationEnum:
 			return NULL;
 		case NoneApproximationEnum:
 			return NULL;
-		case StokesApproximationEnum:
-			return CreatePVectorDiagnosticStokes();
-		case MacAyealPattynApproximationEnum:
-			return CreatePVectorDiagnosticMacAyealPattyn();
-		case MacAyealStokesApproximationEnum:
-			return CreatePVectorDiagnosticMacAyealStokes();
-		case PattynStokesApproximationEnum:
-			return CreatePVectorDiagnosticPattynStokes();
+		case FSApproximationEnum:
+			return CreatePVectorDiagnosticFS();
+		case MacAyealHOApproximationEnum:
+			return CreatePVectorDiagnosticMacAyealHO();
+		case MacAyealFSApproximationEnum:
+			return CreatePVectorDiagnosticMacAyealFS();
+		case HOFSApproximationEnum:
+			return CreatePVectorDiagnosticHOFS();
 		default:
 			_error_("Approximation " << EnumToStringx(approximation) << " not supported yet");
@@ -7744,10 +7744,10 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreatePVectorDiagnosticMacAyealPattyn{{{*/
-ElementVector* Penta::CreatePVectorDiagnosticMacAyealPattyn(void){
+/*FUNCTION Penta::CreatePVectorDiagnosticMacAyealHO{{{*/
+ElementVector* Penta::CreatePVectorDiagnosticMacAyealHO(void){
 
 	/*compute all load vectors for this element*/
 	ElementVector* pe1=CreatePVectorDiagnosticMacAyeal();
-	ElementVector* pe2=CreatePVectorDiagnosticPattyn();
+	ElementVector* pe2=CreatePVectorDiagnosticHO();
 	ElementVector* pe =new ElementVector(pe1,pe2);
 
@@ -7758,11 +7758,11 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreatePVectorDiagnosticMacAyealStokes{{{*/
-ElementVector* Penta::CreatePVectorDiagnosticMacAyealStokes(void){
+/*FUNCTION Penta::CreatePVectorDiagnosticMacAyealFS{{{*/
+ElementVector* Penta::CreatePVectorDiagnosticMacAyealFS(void){
 
 	/*compute all load vectors for this element*/
 	ElementVector* pe1=CreatePVectorDiagnosticMacAyeal();
-	ElementVector* pe2=CreatePVectorDiagnosticStokes();
-	ElementVector* pe3=CreatePVectorCouplingMacAyealStokes();
+	ElementVector* pe2=CreatePVectorDiagnosticFS();
+	ElementVector* pe3=CreatePVectorCouplingMacAyealFS();
 	ElementVector* pe =new ElementVector(pe1,pe2,pe3);
 
@@ -7774,11 +7774,11 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreatePVectorDiagnosticPattynStokes{{{*/
-ElementVector* Penta::CreatePVectorDiagnosticPattynStokes(void){
+/*FUNCTION Penta::CreatePVectorDiagnosticHOFS{{{*/
+ElementVector* Penta::CreatePVectorDiagnosticHOFS(void){
 
 	/*compute all load vectors for this element*/
-	ElementVector* pe1=CreatePVectorDiagnosticPattyn();
-	ElementVector* pe2=CreatePVectorDiagnosticStokes();
-	ElementVector* pe3=CreatePVectorCouplingPattynStokes();
+	ElementVector* pe1=CreatePVectorDiagnosticHO();
+	ElementVector* pe2=CreatePVectorDiagnosticFS();
+	ElementVector* pe3=CreatePVectorCouplingHOFS();
 	ElementVector* pe =new ElementVector(pe1,pe2,pe3);
 
@@ -7790,6 +7790,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreatePVectorDiagnosticHutter{{{*/
-ElementVector* Penta::CreatePVectorDiagnosticHutter(void){
+/*FUNCTION Penta::CreatePVectorDiagnosticSIA{{{*/
+ElementVector* Penta::CreatePVectorDiagnosticSIA(void){
 
 	/*Intermediaries*/
@@ -7903,10 +7903,10 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreatePVectorDiagnosticPattyn{{{*/
-ElementVector* Penta::CreatePVectorDiagnosticPattyn(void){
+/*FUNCTION Penta::CreatePVectorDiagnosticHO{{{*/
+ElementVector* Penta::CreatePVectorDiagnosticHO(void){
 
 	/*compute all load vectors for this element*/
-	ElementVector* pe1=CreatePVectorDiagnosticPattynDrivingStress();
-	ElementVector* pe2=CreatePVectorDiagnosticPattynFront();
+	ElementVector* pe1=CreatePVectorDiagnosticHODrivingStress();
+	ElementVector* pe2=CreatePVectorDiagnosticHOFront();
 	ElementVector* pe =new ElementVector(pe1,pe2);
 
@@ -7917,6 +7917,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreatePVectorDiagnosticPattynDrivingStress{{{*/
-ElementVector* Penta::CreatePVectorDiagnosticPattynDrivingStress(void){
+/*FUNCTION Penta::CreatePVectorDiagnosticHODrivingStress{{{*/
+ElementVector* Penta::CreatePVectorDiagnosticHODrivingStress(void){
 
 	/*Constants*/
@@ -7933,5 +7933,5 @@
 
 	/*Initialize Element vector*/
-	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,PattynApproximationEnum);
+	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,HOApproximationEnum);
 
 	/*Retrieve all inputs and parameters*/
@@ -7965,6 +7965,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreatePVectorDiagnosticPattynFront{{{*/
-ElementVector* Penta::CreatePVectorDiagnosticPattynFront(void){
+/*FUNCTION Penta::CreatePVectorDiagnosticHOFront{{{*/
+ElementVector* Penta::CreatePVectorDiagnosticHOFront(void){
 
 	/*Intermediaries */
@@ -8000,5 +8000,5 @@
 
 	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	ElementVector* pe=new ElementVector(nodes,numnodes,this->parameters,PattynApproximationEnum);
+	ElementVector* pe=new ElementVector(nodes,numnodes,this->parameters,HOApproximationEnum);
 	Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input);
 	rho_water=matpar->GetRhoWater();
@@ -8043,31 +8043,31 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreatePVectorDiagnosticStokes {{{*/
-ElementVector* Penta::CreatePVectorDiagnosticStokes(void){
-
-	int fe_stokes;
+/*FUNCTION Penta::CreatePVectorDiagnosticFS {{{*/
+ElementVector* Penta::CreatePVectorDiagnosticFS(void){
+
+	int fe_FS;
 	ElementVector* pe1;
 	ElementVector* pe2;
 	ElementVector* pe3;
 	ElementVector* pe;
-	parameters->FindParam(&fe_stokes,FlowequationFeStokesEnum);
-
-	switch(fe_stokes){
+	parameters->FindParam(&fe_FS,FlowequationFeFSEnum);
+
+	switch(fe_FS){
 		case 0:
 			/*compute all stiffness matrices for this element*/
-			pe1=CreatePVectorDiagnosticStokesViscous();
-			pe2=CreatePVectorDiagnosticStokesShelf();
-			pe3=CreatePVectorDiagnosticStokesFront();
+			pe1=CreatePVectorDiagnosticFSViscous();
+			pe2=CreatePVectorDiagnosticFSShelf();
+			pe3=CreatePVectorDiagnosticFSFront();
 			pe =new ElementVector(pe1,pe2,pe3);
 			break;
 		case 1:
 			/*compute all stiffness matrices for this element*/
-			pe1=CreatePVectorDiagnosticStokesGLSViscous();
-			pe2=CreatePVectorDiagnosticStokesShelf();
-			pe3=CreatePVectorDiagnosticStokesFront();
+			pe1=CreatePVectorDiagnosticFSGLSViscous();
+			pe2=CreatePVectorDiagnosticFSShelf();
+			pe3=CreatePVectorDiagnosticFSFront();
 			pe =new ElementVector(pe1,pe2,pe3);
 			break;
 		default:
-			_error_("Finite element" << fe_stokes << " not supported yet");
+			_error_("Finite element" << fe_FS << " not supported yet");
 	}
 
@@ -8080,6 +8080,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreatePVectorDiagnosticStokesViscous {{{*/
-ElementVector* Penta::CreatePVectorDiagnosticStokesViscous(void){
+/*FUNCTION Penta::CreatePVectorDiagnosticFSViscous {{{*/
+ElementVector* Penta::CreatePVectorDiagnosticFSViscous(void){
 
 	/*Constants*/
@@ -8090,5 +8090,5 @@
 	int        approximation;
 	IssmDouble Jdet,viscosity;
-	IssmDouble gravity,rho_ice,stokesreconditioning;
+	IssmDouble gravity,rho_ice,FSreconditioning;
 	IssmDouble forcex,forcey,forcez;
 	IssmDouble xyz_list[NUMVERTICES][3];
@@ -8106,9 +8106,9 @@
 	/*Initialize Element vector and return if necessary*/
 	inputs->GetInputValue(&approximation,ApproximationEnum);
-	if(approximation!=StokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum && approximation!=PattynStokesApproximationEnum) return NULL;
-	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
+	if(approximation!=FSApproximationEnum && approximation!=MacAyealFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
+	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
 
 	/*Retrieve all inputs and parameters*/
-	this->parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
+	this->parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
 	rho_ice=matpar->GetRhoIce();
 	gravity=matpar->GetG();
@@ -8128,10 +8128,10 @@
 
 		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
-		GetBStokes(&B[0][0],&xyz_list[0][0],gauss); 
-		GetBprimeStokes(&B_prime[0][0],&xyz_list[0][0], gauss); 
+		GetBFS(&B[0][0],&xyz_list[0][0],gauss); 
+		GetBprimeFS(&B_prime[0][0],&xyz_list[0][0], gauss); 
 		GetNodalFunctionsMINI(&l1l7[0], gauss);
 
 		this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
-		material->GetViscosity3dStokes(&viscosity,&epsilon[0]);
+		material->GetViscosity3dFS(&viscosity,&epsilon[0]);
 
 		loadingforcex_input->GetInputValue(&forcex, gauss);
@@ -8151,5 +8151,5 @@
 		D_scalar=gauss->weight*Jdet;
 		for (i=0;i<6;i++) D[i][i]=D_scalar*2*viscosity;
-		for (i=6;i<8;i++) D[i][i]=-D_scalar*stokesreconditioning;
+		for (i=6;i<8;i++) D[i][i]=-D_scalar*FSreconditioning;
 
 		TripleMultiply(&B[0][0],8,numdofbubble,1,
@@ -8160,5 +8160,5 @@
 
 	/*Condensation*/
-	ReduceVectorStokes(pe->values, &Ke_temp[0][0], &Pe_gaussian[0]);
+	ReduceVectorFS(pe->values, &Ke_temp[0][0], &Pe_gaussian[0]);
 
 	/*Transform coordinate system*/
@@ -8170,6 +8170,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreatePVectorDiagnosticStokesFront{{{*/
-ElementVector* Penta::CreatePVectorDiagnosticStokesFront(void){
+/*FUNCTION Penta::CreatePVectorDiagnosticFSFront{{{*/
+ElementVector* Penta::CreatePVectorDiagnosticFSFront(void){
 
 	/*Intermediaries */
@@ -8205,5 +8205,5 @@
 
 	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	ElementVector* pe=new ElementVector(nodes,numnodes,this->parameters,StokesApproximationEnum);
+	ElementVector* pe=new ElementVector(nodes,numnodes,this->parameters,FSApproximationEnum);
 	rho_water=matpar->GetRhoWater();
 	rho_ice  =matpar->GetRhoIce();
@@ -8247,6 +8247,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreatePVectorDiagnosticStokesGLSViscous {{{*/
-ElementVector* Penta::CreatePVectorDiagnosticStokesGLSViscous(void){
+/*FUNCTION Penta::CreatePVectorDiagnosticFSGLSViscous {{{*/
+ElementVector* Penta::CreatePVectorDiagnosticFSGLSViscous(void){
 
 	/*Constants*/
@@ -8257,5 +8257,5 @@
 	int        approximation;
 	IssmDouble Jdet,gravity,rho_ice,B,D_scalar_stab,viscosity;
-	IssmDouble forcex,forcey,forcez,diameter,stokesreconditioning;
+	IssmDouble forcex,forcey,forcez,diameter,FSreconditioning;
 	IssmDouble xyz_list[NUMVERTICES][3];
 	IssmDouble epsilon[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
@@ -8276,7 +8276,7 @@
 	/*Initialize Element vector and return if necessary*/
 	inputs->GetInputValue(&approximation,ApproximationEnum);
-	if(approximation!=StokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum && approximation!=PattynStokesApproximationEnum) return NULL;
-	parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
-	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
+	if(approximation!=FSApproximationEnum && approximation!=MacAyealFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
+	parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
+	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
 
 	/*Retrieve all inputs and parameters*/
@@ -8318,5 +8318,5 @@
 		GetNodalFunctionsP1(&l1l6[0], gauss);
 		this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
-		material->GetViscosity3dStokes(&viscosity,&epsilon[0]);
+		material->GetViscosity3dFS(&viscosity,&epsilon[0]);
 
 		loadingforcex_input->GetInputValue(&forcex, gauss);
@@ -8372,6 +8372,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreatePVectorDiagnosticStokesShelf{{{*/
-ElementVector* Penta::CreatePVectorDiagnosticStokesShelf(void){
+/*FUNCTION Penta::CreatePVectorDiagnosticFSShelf{{{*/
+ElementVector* Penta::CreatePVectorDiagnosticFSShelf(void){
 
 	/*Intermediaries*/
@@ -8392,6 +8392,6 @@
 	inputs->GetInputValue(&approximation,ApproximationEnum);
 	this->parameters->FindParam(&shelf_dampening,DiagnosticShelfDampeningEnum);
-	if(approximation!=StokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum && approximation!=PattynStokesApproximationEnum) return NULL;
-	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
+	if(approximation!=FSApproximationEnum && approximation!=MacAyealFSApproximationEnum && approximation!=HOFSApproximationEnum) return NULL;
+	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
 
 	/*Retrieve all inputs and parameters*/
@@ -8477,7 +8477,7 @@
 	Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
 	Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
-	Input* vzstokes_input=NULL;
-	if(approximation==PattynStokesApproximationEnum || approximation==MacAyealStokesApproximationEnum){
-		vzstokes_input=inputs->GetInput(VzStokesEnum); _assert_(vzstokes_input);
+	Input* vzFS_input=NULL;
+	if(approximation==HOFSApproximationEnum || approximation==MacAyealFSApproximationEnum){
+		vzFS_input=inputs->GetInput(VzFSEnum); _assert_(vzFS_input);
 	}
 
@@ -8493,6 +8493,6 @@
 		vx_input->GetInputDerivativeValue(&du[0],&xyz_list[0][0],gauss);
 		vy_input->GetInputDerivativeValue(&dv[0],&xyz_list[0][0],gauss);
-		if(approximation==PattynStokesApproximationEnum || approximation==MacAyealStokesApproximationEnum){
-			vzstokes_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
+		if(approximation==HOFSApproximationEnum || approximation==MacAyealFSApproximationEnum){
+			vzFS_input->GetInputDerivativeValue(&dw[0],&xyz_list[0][0],gauss);
 			dwdz=dw[2];
 		}
@@ -8539,7 +8539,7 @@
 	Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
 	Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
-	Input* vzstokes_input=NULL;
-	if(approximation==PattynStokesApproximationEnum || approximation==MacAyealStokesApproximationEnum){
-		vzstokes_input=inputs->GetInput(VzStokesEnum);       _assert_(vzstokes_input);
+	Input* vzFS_input=NULL;
+	if(approximation==HOFSApproximationEnum || approximation==MacAyealFSApproximationEnum){
+		vzFS_input=inputs->GetInput(VzFSEnum);       _assert_(vzFS_input);
 	}
 
@@ -8554,6 +8554,6 @@
 		vx_input->GetInputValue(&vx, gauss);
 		vy_input->GetInputValue(&vy, gauss);
-		if(approximation==PattynStokesApproximationEnum || approximation==MacAyealStokesApproximationEnum){
-			vzstokes_input->GetInputValue(&vz, gauss);
+		if(approximation==HOFSApproximationEnum || approximation==MacAyealFSApproximationEnum){
+			vzFS_input->GetInputValue(&vz, gauss);
 		}
 		else vz=0;
@@ -8581,9 +8581,9 @@
 	switch(approximation){
 		case MacAyealApproximationEnum:
-			return CreateJacobianDiagnosticMacayeal2d();
-		case PattynApproximationEnum:
-			return CreateJacobianDiagnosticPattyn();
-		case StokesApproximationEnum:
-			return CreateJacobianDiagnosticStokes();
+			return CreateJacobianDiagnosticSSA2d();
+		case HOApproximationEnum:
+			return CreateJacobianDiagnosticHO();
+		case FSApproximationEnum:
+			return CreateJacobianDiagnosticFS();
 		case NoneApproximationEnum:
 			return NULL;
@@ -8593,6 +8593,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreateJacobianDiagnosticMacayeal2d{{{*/
-ElementMatrix* Penta::CreateJacobianDiagnosticMacayeal2d(void){
+/*FUNCTION Penta::CreateJacobianDiagnosticSSA2d{{{*/
+ElementMatrix* Penta::CreateJacobianDiagnosticSSA2d(void){
 
 	/*Figure out if this penta is collapsed. If so, then bailout, except if it is at the 
@@ -8616,5 +8616,5 @@
 	/*Call Tria function*/
 	Tria* tria=(Tria*)SpawnTria(0,1,2); //nodes 0, 1 and 2 make the new tria.
-	ElementMatrix* Ke=tria->CreateJacobianDiagnosticMacayeal();
+	ElementMatrix* Ke=tria->CreateJacobianDiagnosticSSA();
 	delete tria->material; delete tria;
 
@@ -8627,6 +8627,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreateJacobianDiagnosticPattyn{{{*/
-ElementMatrix* Penta::CreateJacobianDiagnosticPattyn(void){
+/*FUNCTION Penta::CreateJacobianDiagnosticHO{{{*/
+ElementMatrix* Penta::CreateJacobianDiagnosticHO(void){
 
 	/*Constants*/
@@ -8645,6 +8645,6 @@
 	GaussPenta *gauss=NULL;
 
-	/*Initialize Jacobian with regular Pattyn (first part of the Gateau derivative)*/
-	ElementMatrix* Ke=CreateKMatrixDiagnosticPattyn();
+	/*Initialize Jacobian with regular HO (first part of the Gateau derivative)*/
+	ElementMatrix* Ke=CreateKMatrixDiagnosticHO();
 
 	/*Retrieve all inputs and parameters*/
@@ -8662,5 +8662,5 @@
 		GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss);
 
-		this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
+		this->GetStrainRate3dHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
 		material->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);
 		eps1[0]=2*epsilon[0]+epsilon[1];   eps2[0]=epsilon[2];
@@ -8691,6 +8691,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::CreateJacobianDiagnosticStokes{{{*/
-ElementMatrix* Penta::CreateJacobianDiagnosticStokes(void){
+/*FUNCTION Penta::CreateJacobianDiagnosticFS{{{*/
+ElementMatrix* Penta::CreateJacobianDiagnosticFS(void){
 
 	/*Constants*/
@@ -8710,6 +8710,6 @@
 	GaussPenta *gauss=NULL;
 
-	/*Initialize Jacobian with regular Stokes (first part of the Gateau derivative)*/
-	ElementMatrix* Ke=CreateKMatrixDiagnosticStokes();
+	/*Initialize Jacobian with regular FS (first part of the Gateau derivative)*/
+	ElementMatrix* Ke=CreateKMatrixDiagnosticFS();
 
 	/*Retrieve all inputs and parameters*/
@@ -8728,5 +8728,5 @@
 		GetNodalFunctionsP1Derivatives(&dphi[0][0],&xyz_list[0][0],gauss);
 
-		this->GetStrainRate3dPattyn(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
+		this->GetStrainRate3dHO(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
 		material->GetViscosityDerivativeEpsSquare(&mu_prime,&epsilon[0]);
 		eps1[0]=epsilon[0];   eps2[0]=epsilon[2];   eps3[0]=epsilon[3];
@@ -8784,5 +8784,5 @@
 
 	/*If the element is a coupling, do nothing: every node is also on an other elements 
-	 * (as coupling is between MacAyeal and Pattyn) so the other element will take care of it*/
+	 * (as coupling is between MacAyeal and HO) so the other element will take care of it*/
 	GetDofList(&doflist,approximation,GsetEnum);
 
@@ -8808,6 +8808,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::GetSolutionFromInputsDiagnosticHutter{{{*/
-void  Penta::GetSolutionFromInputsDiagnosticHutter(Vector<IssmDouble>* solution){
+/*FUNCTION Penta::GetSolutionFromInputsDiagnosticSIA{{{*/
+void  Penta::GetSolutionFromInputsDiagnosticSIA(Vector<IssmDouble>* solution){
 
 	const int    numdof=NDOF2*NUMVERTICES;
@@ -8877,6 +8877,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::GetSolutionFromInputsDiagnosticStokes{{{*/
-void  Penta::GetSolutionFromInputsDiagnosticStokes(Vector<IssmDouble>* solution){
+/*FUNCTION Penta::GetSolutionFromInputsDiagnosticFS{{{*/
+void  Penta::GetSolutionFromInputsDiagnosticFS(Vector<IssmDouble>* solution){
 
 	const int    numdof=NDOF4*NUMVERTICES;
@@ -8885,10 +8885,10 @@
 	int*         doflist=NULL;
 	IssmDouble       vx,vy,vz,p;
-	IssmDouble       stokesreconditioning;
+	IssmDouble       FSreconditioning;
 	IssmDouble       values[numdof];
 	GaussPenta   *gauss;
 
 	/*Get dof list: */
-	GetDofList(&doflist,StokesApproximationEnum,GsetEnum);
+	GetDofList(&doflist,FSApproximationEnum,GsetEnum);
 	Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
 	Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
@@ -8897,5 +8897,5 @@
 
 	/*Recondition pressure: */
-	this->parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
+	this->parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
 
 	/*Ok, we have vx vy vz and P in values, fill in vx vy vz P arrays: */
@@ -8911,5 +8911,5 @@
 		values[i*NDOF4+1]=vy;
 		values[i*NDOF4+2]=vz;
-		values[i*NDOF4+3]=p/stokesreconditioning;
+		values[i*NDOF4+3]=p/FSreconditioning;
 	}
 
@@ -8961,6 +8961,6 @@
 
 	/* Get eps_b*/
-	vx_input->GetVxStrainRate3dPattyn(epsilonvx,xyz_list,gauss);
-	vy_input->GetVyStrainRate3dPattyn(epsilonvy,xyz_list,gauss);
+	vx_input->GetVxStrainRate3dHO(epsilonvx,xyz_list,gauss);
+	vy_input->GetVyStrainRate3dHO(epsilonvy,xyz_list,gauss);
 	for(i=0;i<5;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i];
 	eps_b = sqrt(epsilon[0]*epsilon[0] + epsilon[1]*epsilon[1] + epsilon[0]*epsilon[1] + epsilon[2]*epsilon[2]);
@@ -9015,18 +9015,18 @@
 		return;
 	}
-	else if (approximation==PattynApproximationEnum){
-		InputUpdateFromSolutionDiagnosticPattyn(solution);
-	}
-	else if (approximation==PattynStokesApproximationEnum){
-		InputUpdateFromSolutionDiagnosticPattynStokes(solution);
-	}
-	else if (approximation==MacAyealStokesApproximationEnum){
-		InputUpdateFromSolutionDiagnosticMacAyealStokes(solution);
-	}
-	else if (approximation==StokesApproximationEnum || approximation==NoneApproximationEnum){
-		InputUpdateFromSolutionDiagnosticStokes(solution);
-	}
-	else if (approximation==MacAyealPattynApproximationEnum){
-		InputUpdateFromSolutionDiagnosticMacAyealPattyn(solution);
+	else if (approximation==HOApproximationEnum){
+		InputUpdateFromSolutionDiagnosticHO(solution);
+	}
+	else if (approximation==HOFSApproximationEnum){
+		InputUpdateFromSolutionDiagnosticHOFS(solution);
+	}
+	else if (approximation==MacAyealFSApproximationEnum){
+		InputUpdateFromSolutionDiagnosticMacAyealFS(solution);
+	}
+	else if (approximation==FSApproximationEnum || approximation==NoneApproximationEnum){
+		InputUpdateFromSolutionDiagnosticFS(solution);
+	}
+	else if (approximation==MacAyealHOApproximationEnum){
+		InputUpdateFromSolutionDiagnosticMacAyealHO(solution);
 	}
 }
@@ -9113,6 +9113,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticMacAyealPattyn {{{*/
-void  Penta::InputUpdateFromSolutionDiagnosticMacAyealPattyn(IssmDouble* solution){
+/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticMacAyealHO {{{*/
+void  Penta::InputUpdateFromSolutionDiagnosticMacAyealHO(IssmDouble* solution){
 
 	const int    numdof=NDOF2*NUMVERTICES;
@@ -9121,6 +9121,6 @@
 	int     i;
 	IssmDouble  rho_ice,g;
-	IssmDouble  macayeal_values[numdof];
-	IssmDouble  pattyn_values[numdof];
+	IssmDouble  SSA_values[numdof];
+	IssmDouble  HO_values[numdof];
 	IssmDouble  vx[NUMVERTICES];
 	IssmDouble  vy[NUMVERTICES];
@@ -9134,10 +9134,10 @@
 	Penta   *penta   = NULL;
 
-	/*OK, we have to add results of this element for pattyn 
-	 * and results from the penta at base for macayeal. Now recover results*/
+	/*OK, we have to add results of this element for HO 
+	 * and results from the penta at base for SSA. Now recover results*/
 	penta=GetBasalElement();
 
-	/*Get dof listof this element (pattyn dofs) and of the penta at base (macayeal dofs): */
-	GetDofList(&doflistp,PattynApproximationEnum,GsetEnum);
+	/*Get dof listof this element (HO dofs) and of the penta at base (SSA dofs): */
+	GetDofList(&doflistp,HOApproximationEnum,GsetEnum);
 	penta->GetDofList(&doflistm,MacAyealApproximationEnum,GsetEnum);
 
@@ -9147,20 +9147,20 @@
 	/*Use the dof list to index into the solution vector: */
 	for(i=0;i<numdof2d;i++){
-		pattyn_values[i]=solution[doflistp[i]];
-		macayeal_values[i]=solution[doflistm[i]];
+		HO_values[i]=solution[doflistp[i]];
+		SSA_values[i]=solution[doflistm[i]];
 	}
 	for(i=numdof2d;i<numdof;i++){
-		pattyn_values[i]=solution[doflistp[i]];
-		macayeal_values[i]=macayeal_values[i-numdof2d];
+		HO_values[i]=solution[doflistp[i]];
+		SSA_values[i]=SSA_values[i-numdof2d];
 	}
 
 	/*Transform solution in Cartesian Space*/
-	TransformSolutionCoord(&macayeal_values[0],penta->nodes,NUMVERTICES,XYEnum);
-	TransformSolutionCoord(&pattyn_values[0],   this->nodes,NUMVERTICES,XYEnum);
+	TransformSolutionCoord(&SSA_values[0],penta->nodes,NUMVERTICES,XYEnum);
+	TransformSolutionCoord(&HO_values[0],   this->nodes,NUMVERTICES,XYEnum);
 
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
 	for(i=0;i<NUMVERTICES;i++){
-		vx[i]=macayeal_values[i*NDOF2+0]+pattyn_values[i*NDOF2+0];
-		vy[i]=macayeal_values[i*NDOF2+1]+pattyn_values[i*NDOF2+1];
+		vx[i]=SSA_values[i*NDOF2+0]+HO_values[i*NDOF2+0];
+		vy[i]=SSA_values[i*NDOF2+1]+HO_values[i*NDOF2+1];
 
 		/*Check solution*/
@@ -9197,6 +9197,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticMacAyealStokes {{{*/
-void  Penta::InputUpdateFromSolutionDiagnosticMacAyealStokes(IssmDouble* solution){
+/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticMacAyealFS {{{*/
+void  Penta::InputUpdateFromSolutionDiagnosticMacAyealFS(IssmDouble* solution){
 
 	const int    numdofm=NDOF2*NUMVERTICES;
@@ -9205,12 +9205,12 @@
 
 	int     i;
-	IssmDouble  stokesreconditioning;
-	IssmDouble  macayeal_values[numdofm];
-	IssmDouble  stokes_values[numdofs];
+	IssmDouble  FSreconditioning;
+	IssmDouble  SSA_values[numdofm];
+	IssmDouble  FS_values[numdofs];
 	IssmDouble  vx[NUMVERTICES];
 	IssmDouble  vy[NUMVERTICES];
 	IssmDouble  vz[NUMVERTICES];
-	IssmDouble  vzmacayeal[NUMVERTICES];
-	IssmDouble  vzstokes[NUMVERTICES];
+	IssmDouble  vzSSA[NUMVERTICES];
+	IssmDouble  vzFS[NUMVERTICES];
 	IssmDouble  vel[NUMVERTICES];
 	IssmDouble  pressure[NUMVERTICES];
@@ -9220,12 +9220,12 @@
 	Penta   *penta          = NULL;
 
-	/*OK, we have to add results of this element for macayeal 
-	 * and results from the penta at base for macayeal. Now recover results*/
+	/*OK, we have to add results of this element for SSA 
+	 * and results from the penta at base for SSA. Now recover results*/
 	penta=GetBasalElement();
 
-	/*Get dof listof this element (macayeal dofs) and of the penta at base (macayeal dofs): */
+	/*Get dof listof this element (SSA dofs) and of the penta at base (SSA dofs): */
 	penta->GetDofList(&doflistm,MacAyealApproximationEnum,GsetEnum);
-	GetDofList(&doflists,StokesApproximationEnum,GsetEnum);
-	this->parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
+	GetDofList(&doflists,FSApproximationEnum,GsetEnum);
+	this->parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
 
 	/*Get node data: */
@@ -9234,36 +9234,36 @@
 	/*Use the dof list to index into the solution vector: */
 	for(i=0;i<numdof2d;i++){
-		macayeal_values[i]=solution[doflistm[i]];
-		macayeal_values[i+numdof2d]=solution[doflistm[i]];
+		SSA_values[i]=solution[doflistm[i]];
+		SSA_values[i+numdof2d]=solution[doflistm[i]];
 	}
 	for(i=0;i<numdofs;i++){
-		stokes_values[i]=solution[doflists[i]];
+		FS_values[i]=solution[doflists[i]];
 	}
 
 	/*Transform solution in Cartesian Space*/
-	TransformSolutionCoord(&macayeal_values[0],this->nodes,NUMVERTICES,XYEnum);
-	TransformSolutionCoord(&stokes_values[0],this->nodes,NUMVERTICES,XYZPEnum);
+	TransformSolutionCoord(&SSA_values[0],this->nodes,NUMVERTICES,XYEnum);
+	TransformSolutionCoord(&FS_values[0],this->nodes,NUMVERTICES,XYZPEnum);
 
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
 	for(i=0;i<NUMVERTICES;i++){
-		vx[i]=stokes_values[i*NDOF4+0]+macayeal_values[i*NDOF2+0];
-		vy[i]=stokes_values[i*NDOF4+1]+macayeal_values[i*NDOF2+1];
-		vzstokes[i]=stokes_values[i*NDOF4+2];
-		pressure[i]=stokes_values[i*NDOF4+3]*stokesreconditioning;
+		vx[i]=FS_values[i*NDOF4+0]+SSA_values[i*NDOF2+0];
+		vy[i]=FS_values[i*NDOF4+1]+SSA_values[i*NDOF2+1];
+		vzFS[i]=FS_values[i*NDOF4+2];
+		pressure[i]=FS_values[i*NDOF4+3]*FSreconditioning;
 
 		/*Check solution*/
 		if(xIsNan<IssmDouble>(vx[i]))       _error_("NaN found in solution vector");
 		if(xIsNan<IssmDouble>(vy[i]))       _error_("NaN found in solution vector");
-		if(xIsNan<IssmDouble>(vzstokes[i])) _error_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vzFS[i])) _error_("NaN found in solution vector");
 		if(xIsNan<IssmDouble>(pressure[i])) _error_("NaN found in solution vector");
 	}
 
 	/*Get Vz*/
-	Input* vzmacayeal_input=inputs->GetInput(VzMacAyealEnum);
-	if (vzmacayeal_input){
-		if (vzmacayeal_input->ObjectEnum()!=PentaInputEnum){
-			_error_("Cannot compute Vel as VzMacAyeal is of type " << EnumToStringx(vzmacayeal_input->ObjectEnum()));
-		}
-		GetInputListOnVertices(&vzmacayeal[0],VzMacAyealEnum);
+	Input* vzSSA_input=inputs->GetInput(VzMacAyealEnum);
+	if (vzSSA_input){
+		if (vzSSA_input->ObjectEnum()!=PentaInputEnum){
+			_error_("Cannot compute Vel as VzMacAyeal is of type " << EnumToStringx(vzSSA_input->ObjectEnum()));
+		}
+		GetInputListOnVertices(&vzSSA[0],VzMacAyealEnum);
 	}
 	else{
@@ -9273,5 +9273,5 @@
 	/*Now Compute vel*/
 	for(i=0;i<NUMVERTICES;i++) {
-		vz[i]=vzmacayeal[i]+vzstokes[i];
+		vz[i]=vzSSA[i]+vzFS[i];
 		vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
 	}
@@ -9288,5 +9288,5 @@
 	this->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum));
 	this->inputs->AddInput(new PentaInput(VzEnum,vz,P1Enum));
-	this->inputs->AddInput(new PentaInput(VzStokesEnum,vzstokes,P1Enum));
+	this->inputs->AddInput(new PentaInput(VzFSEnum,vzFS,P1Enum));
 	this->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum));
 	this->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum));
@@ -9378,6 +9378,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticPattyn {{{*/
-void  Penta::InputUpdateFromSolutionDiagnosticPattyn(IssmDouble* solution){
+/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticHO {{{*/
+void  Penta::InputUpdateFromSolutionDiagnosticHO(IssmDouble* solution){
 
 	const int    numdof=NDOF2*NUMVERTICES;
@@ -9396,5 +9396,5 @@
 
 	/*Get dof list: */
-	GetDofList(&doflist,PattynApproximationEnum,GsetEnum);
+	GetDofList(&doflist,HOApproximationEnum,GsetEnum);
 
 	/*Get node data: */
@@ -9452,6 +9452,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticPattynStokes {{{*/
-void  Penta::InputUpdateFromSolutionDiagnosticPattynStokes(IssmDouble* solution){
+/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticHOFS {{{*/
+void  Penta::InputUpdateFromSolutionDiagnosticHOFS(IssmDouble* solution){
 
 	const int    numdofp=NDOF2*NUMVERTICES;
@@ -9459,27 +9459,27 @@
 
 	int    i;
-	IssmDouble pattyn_values[numdofp];
-	IssmDouble stokes_values[numdofs];
+	IssmDouble HO_values[numdofp];
+	IssmDouble FS_values[numdofs];
 	IssmDouble vx[NUMVERTICES];
 	IssmDouble vy[NUMVERTICES];
 	IssmDouble vz[NUMVERTICES];
-	IssmDouble vzpattyn[NUMVERTICES];
-	IssmDouble vzstokes[NUMVERTICES];
+	IssmDouble vzHO[NUMVERTICES];
+	IssmDouble vzFS[NUMVERTICES];
 	IssmDouble vel[NUMVERTICES];
 	IssmDouble pressure[NUMVERTICES];
 	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble stokesreconditioning;
+	IssmDouble FSreconditioning;
 	int*   doflistp      = NULL;
 	int*   doflists      = NULL;
 	Penta  *penta        = NULL;
 
-	/*OK, we have to add results of this element for pattyn 
-	 * and results from the penta at base for macayeal. Now recover results*/
+	/*OK, we have to add results of this element for HO 
+	 * and results from the penta at base for SSA. Now recover results*/
 	penta=GetBasalElement();
 
-	/*Get dof listof this element (pattyn dofs) and of the penta at base (macayeal dofs): */
-	GetDofList(&doflistp,PattynApproximationEnum,GsetEnum);
-	GetDofList(&doflists,StokesApproximationEnum,GsetEnum);
-	this->parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
+	/*Get dof listof this element (HO dofs) and of the penta at base (SSA dofs): */
+	GetDofList(&doflistp,HOApproximationEnum,GsetEnum);
+	GetDofList(&doflists,FSApproximationEnum,GsetEnum);
+	this->parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
 
 	/*Get node data: */
@@ -9487,40 +9487,40 @@
 
 	/*Use the dof list to index into the solution vector: */
-	for(i=0;i<numdofp;i++) pattyn_values[i]=solution[doflistp[i]];
-	for(i=0;i<numdofs;i++) stokes_values[i]=solution[doflists[i]];
+	for(i=0;i<numdofp;i++) HO_values[i]=solution[doflistp[i]];
+	for(i=0;i<numdofs;i++) FS_values[i]=solution[doflists[i]];
 
 	/*Transform solution in Cartesian Space*/
-	TransformSolutionCoord(&pattyn_values[0],this->nodes,NUMVERTICES,XYEnum);
-	TransformSolutionCoord(&stokes_values[0],this->nodes,NUMVERTICES,XYZPEnum);
+	TransformSolutionCoord(&HO_values[0],this->nodes,NUMVERTICES,XYEnum);
+	TransformSolutionCoord(&FS_values[0],this->nodes,NUMVERTICES,XYZPEnum);
 
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
 	for(i=0;i<NUMVERTICES;i++){
-		vx[i]=stokes_values[i*NDOF4+0]+pattyn_values[i*NDOF2+0];
-		vy[i]=stokes_values[i*NDOF4+1]+pattyn_values[i*NDOF2+1];
-		vzstokes[i]=stokes_values[i*NDOF4+2];
-		pressure[i]=stokes_values[i*NDOF4+3]*stokesreconditioning;
+		vx[i]=FS_values[i*NDOF4+0]+HO_values[i*NDOF2+0];
+		vy[i]=FS_values[i*NDOF4+1]+HO_values[i*NDOF2+1];
+		vzFS[i]=FS_values[i*NDOF4+2];
+		pressure[i]=FS_values[i*NDOF4+3]*FSreconditioning;
 
 		/*Check solution*/
 		if(xIsNan<IssmDouble>(vx[i]))       _error_("NaN found in solution vector");
 		if(xIsNan<IssmDouble>(vy[i]))       _error_("NaN found in solution vector");
-		if(xIsNan<IssmDouble>(vzstokes[i])) _error_("NaN found in solution vector");
+		if(xIsNan<IssmDouble>(vzFS[i])) _error_("NaN found in solution vector");
 		if(xIsNan<IssmDouble>(pressure[i])) _error_("NaN found in solution vector");
 	}
 
 	/*Get Vz*/
-	Input* vzpattyn_input=inputs->GetInput(VzPattynEnum);
-	if (vzpattyn_input){
-		if (vzpattyn_input->ObjectEnum()!=PentaInputEnum){
-			_error_("Cannot compute Vel as VzPattyn is of type " << EnumToStringx(vzpattyn_input->ObjectEnum()));
-		}
-		GetInputListOnVertices(&vzpattyn[0],VzPattynEnum);
+	Input* vzHO_input=inputs->GetInput(VzHOEnum);
+	if (vzHO_input){
+		if (vzHO_input->ObjectEnum()!=PentaInputEnum){
+			_error_("Cannot compute Vel as VzHO is of type " << EnumToStringx(vzHO_input->ObjectEnum()));
+		}
+		GetInputListOnVertices(&vzHO[0],VzHOEnum);
 	}
 	else{
-		_error_("Cannot update solution as VzPattyn is not present");
+		_error_("Cannot update solution as VzHO is not present");
 	}
 
 	/*Now Compute vel*/
 	for(i=0;i<NUMVERTICES;i++) {
-		vz[i]=vzpattyn[i]+vzstokes[i];
+		vz[i]=vzHO[i]+vzFS[i];
 		vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
 	}
@@ -9537,5 +9537,5 @@
 	this->inputs->AddInput(new PentaInput(VyEnum,vy,P1Enum));
 	this->inputs->AddInput(new PentaInput(VzEnum,vz,P1Enum));
-	this->inputs->AddInput(new PentaInput(VzStokesEnum,vzstokes,P1Enum));
+	this->inputs->AddInput(new PentaInput(VzFSEnum,vzFS,P1Enum));
 	this->inputs->AddInput(new PentaInput(VelEnum,vel,P1Enum));
 	this->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum));
@@ -9546,6 +9546,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticHutter {{{*/
-void  Penta::InputUpdateFromSolutionDiagnosticHutter(IssmDouble* solution){
+/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticSIA {{{*/
+void  Penta::InputUpdateFromSolutionDiagnosticSIA(IssmDouble* solution){
 
 	const int    numdof=NDOF2*NUMVERTICES;
@@ -9621,7 +9621,7 @@
 	IssmDouble   vy[NUMVERTICES];
 	IssmDouble   vz[NUMVERTICES];
-	IssmDouble   vzmacayeal[NUMVERTICES];
-	IssmDouble   vzpattyn[NUMVERTICES];
-	IssmDouble   vzstokes[NUMVERTICES];
+	IssmDouble   vzSSA[NUMVERTICES];
+	IssmDouble   vzHO[NUMVERTICES];
+	IssmDouble   vzFS[NUMVERTICES];
 	IssmDouble   vel[NUMVERTICES];
 	IssmDouble   pressure[NUMVERTICES];
@@ -9630,7 +9630,7 @@
 	int*     doflist      = NULL;
 
-	/*Get the approximation and do nothing if the element in Stokes or None*/
+	/*Get the approximation and do nothing if the element in FS or None*/
 	inputs->GetInputValue(&approximation,ApproximationEnum);
-	if(approximation==StokesApproximationEnum || approximation==NoneApproximationEnum){
+	if(approximation==FSApproximationEnum || approximation==NoneApproximationEnum){
 		return;
 	}
@@ -9653,27 +9653,27 @@
 	GetInputListOnVertices(&vy[0],VyEnum,0.0); //default is 0
 
-	/*Do some modifications if we actually have a PattynStokes or MacAyealStokes element*/
-	if(approximation==PattynStokesApproximationEnum){
-		Input* vzstokes_input=inputs->GetInput(VzStokesEnum);
-		if (vzstokes_input){
-			if (vzstokes_input->ObjectEnum()!=PentaInputEnum) _error_("Cannot compute Vel as VzStokes is of type " << EnumToStringx(vzstokes_input->ObjectEnum()));
-			GetInputListOnVertices(&vzstokes[0],VzStokesEnum);
-		}
-		else _error_("Cannot compute Vz as VzStokes in not present in PattynStokes element");
+	/*Do some modifications if we actually have a HOFS or MacAyealFS element*/
+	if(approximation==HOFSApproximationEnum){
+		Input* vzFS_input=inputs->GetInput(VzFSEnum);
+		if (vzFS_input){
+			if (vzFS_input->ObjectEnum()!=PentaInputEnum) _error_("Cannot compute Vel as VzFS is of type " << EnumToStringx(vzFS_input->ObjectEnum()));
+			GetInputListOnVertices(&vzFS[0],VzFSEnum);
+		}
+		else _error_("Cannot compute Vz as VzFS in not present in HOFS element");
 		for(i=0;i<NUMVERTICES;i++){
-			vzpattyn[i]=vz[i];
-			vz[i]=vzpattyn[i]+vzstokes[i];
-		}
-	}
-	else if(approximation==MacAyealStokesApproximationEnum){
-		Input* vzstokes_input=inputs->GetInput(VzStokesEnum);
-		if (vzstokes_input){
-			if (vzstokes_input->ObjectEnum()!=PentaInputEnum) _error_("Cannot compute Vel as VzStokes is of type " << EnumToStringx(vzstokes_input->ObjectEnum()));
-			GetInputListOnVertices(&vzstokes[0],VzStokesEnum);
-		}
-		else _error_("Cannot compute Vz as VzStokes in not present in MacAyealStokes element");
+			vzHO[i]=vz[i];
+			vz[i]=vzHO[i]+vzFS[i];
+		}
+	}
+	else if(approximation==MacAyealFSApproximationEnum){
+		Input* vzFS_input=inputs->GetInput(VzFSEnum);
+		if (vzFS_input){
+			if (vzFS_input->ObjectEnum()!=PentaInputEnum) _error_("Cannot compute Vel as VzFS is of type " << EnumToStringx(vzFS_input->ObjectEnum()));
+			GetInputListOnVertices(&vzFS[0],VzFSEnum);
+		}
+		else _error_("Cannot compute Vz as VzFS in not present in MacAyealFS element");
 		for(i=0;i<NUMVERTICES;i++){
-			vzmacayeal[i]=vz[i];
-			vz[i]=vzmacayeal[i]+vzstokes[i];
+			vzSSA[i]=vz[i];
+			vz[i]=vzSSA[i]+vzFS[i];
 		}
 	}
@@ -9683,6 +9683,6 @@
 
 	/*For pressure: we have not computed pressure in this analysis, for this element. We are in 3D, 
-	 *so the pressure is just the pressure at the z elevation: except it this is a PattynStokes element */
-	if(approximation!=PattynStokesApproximationEnum &&  approximation!=MacAyealStokesApproximationEnum){
+	 *so the pressure is just the pressure at the z elevation: except it this is a HOFS element */
+	if(approximation!=HOFSApproximationEnum &&  approximation!=MacAyealFSApproximationEnum){
 		rho_ice=matpar->GetRhoIce();
 		g=matpar->GetG();
@@ -9695,13 +9695,13 @@
 	this->inputs->ChangeEnum(VzEnum,VzPicardEnum);
 
-	if(approximation!=PattynStokesApproximationEnum && approximation!=MacAyealStokesApproximationEnum){
+	if(approximation!=HOFSApproximationEnum && approximation!=MacAyealFSApproximationEnum){
 		this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);
 		this->inputs->AddInput(new PentaInput(PressureEnum,pressure,P1Enum));
 	}
-	else if(approximation==PattynStokesApproximationEnum){
-		this->inputs->AddInput(new PentaInput(VzPattynEnum,vzpattyn,P1Enum));
-	}
-	else if(approximation==MacAyealStokesApproximationEnum){
-		this->inputs->AddInput(new PentaInput(VzMacAyealEnum,vzmacayeal,P1Enum));
+	else if(approximation==HOFSApproximationEnum){
+		this->inputs->AddInput(new PentaInput(VzHOEnum,vzHO,P1Enum));
+	}
+	else if(approximation==MacAyealFSApproximationEnum){
+		this->inputs->AddInput(new PentaInput(VzMacAyealEnum,vzSSA,P1Enum));
 	}
 	this->inputs->AddInput(new PentaInput(VzEnum,vz,P1Enum));
@@ -9712,6 +9712,6 @@
 }
 /*}}}*/
-/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticStokes {{{*/
-void  Penta::InputUpdateFromSolutionDiagnosticStokes(IssmDouble* solution){
+/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticFS {{{*/
+void  Penta::InputUpdateFromSolutionDiagnosticFS(IssmDouble* solution){
 
 	const int numdof=NDOF4*NUMVERTICES;
@@ -9724,9 +9724,9 @@
 	IssmDouble  vel[NUMVERTICES];
 	IssmDouble  pressure[NUMVERTICES];
-	IssmDouble  stokesreconditioning;
+	IssmDouble  FSreconditioning;
 	int*    doflist=NULL;
 
 	/*Get dof list: */
-	GetDofList(&doflist,StokesApproximationEnum,GsetEnum);
+	GetDofList(&doflist,FSApproximationEnum,GsetEnum);
 
 	/*Use the dof list to index into the solution vector: */
@@ -9751,6 +9751,6 @@
 
 	/*Recondition pressure and compute vel: */
-	this->parameters->FindParam(&stokesreconditioning,DiagnosticStokesreconditioningEnum);
-	for(i=0;i<NUMVERTICES;i++) pressure[i]=pressure[i]*stokesreconditioning;
+	this->parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
+	for(i=0;i<NUMVERTICES;i++) pressure[i]=pressure[i]*FSreconditioning;
 	for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
 
Index: issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 15556)
+++ issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 15564)
@@ -145,9 +145,9 @@
 		void   Gradj(Vector<IssmDouble>* gradient,int control_type,int control_index);
 		void   GradjDragMacAyeal(Vector<IssmDouble>* gradient,int control_index);
-		void   GradjDragPattyn(Vector<IssmDouble>* gradient,int control_index);
-		void   GradjDragStokes(Vector<IssmDouble>* gradient,int control_index);
+		void   GradjDragHO(Vector<IssmDouble>* gradient,int control_index);
+		void   GradjDragFS(Vector<IssmDouble>* gradient,int control_index);
 		void   GradjBbarMacAyeal(Vector<IssmDouble>* gradient,int control_index);
-		void   GradjBbarPattyn(Vector<IssmDouble>* gradient,int control_index);
-		void   GradjBbarStokes(Vector<IssmDouble>* gradient,int control_index);
+		void   GradjBbarHO(Vector<IssmDouble>* gradient,int control_index);
+		void   GradjBbarFS(Vector<IssmDouble>* gradient,int control_index);
 		void   GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data);
 		void   SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index);
@@ -197,5 +197,5 @@
 		void	         GetSolutionFromInputsEnthalpy(Vector<IssmDouble>* solutiong);
 		IssmDouble     GetStabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa);
-		void    GetStrainRate3dPattyn(IssmDouble* epsilon,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input);
+		void    GetStrainRate3dHO(IssmDouble* epsilon,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input);
 		void    GetStrainRate3d(IssmDouble* epsilon,IssmDouble* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input, Input* vz_input);
 		Penta*  GetUpperElement(void);
@@ -215,6 +215,6 @@
 		bool    IsOnWater(void); 
 		IssmDouble  MinEdgeLength(IssmDouble xyz_list[6][3]);
-		void	  ReduceMatrixStokes(IssmDouble* Ke_reduced, IssmDouble* Ke_temp);
-		void	  ReduceVectorStokes(IssmDouble* Pe_reduced, IssmDouble* Ke_temp, IssmDouble* Pe_temp);
+		void	  ReduceMatrixFS(IssmDouble* Ke_reduced, IssmDouble* Ke_temp);
+		void	  ReduceVectorFS(IssmDouble* Pe_reduced, IssmDouble* Ke_temp, IssmDouble* Pe_temp);
 		void	  SetClone(int* minranks);
 		Tria*	  SpawnTria(int g0, int g1, int g2);
@@ -223,76 +223,76 @@
 
 		#ifdef _HAVE_DIAGNOSTIC_
-		ElementMatrix* CreateKMatrixCouplingMacAyealPattyn(void);
-		ElementMatrix* CreateKMatrixCouplingMacAyealPattynViscous(void);
-		ElementMatrix* CreateKMatrixCouplingMacAyealPattynFriction(void);
-		ElementMatrix* CreateKMatrixCouplingMacAyealStokes(void);
-		ElementMatrix* CreateKMatrixCouplingMacAyealStokesViscous(void);
-		ElementMatrix* CreateKMatrixCouplingMacAyealStokesFriction(void);
-		ElementMatrix* CreateKMatrixCouplingPattynStokes(void);
+		ElementMatrix* CreateKMatrixCouplingMacAyealHO(void);
+		ElementMatrix* CreateKMatrixCouplingMacAyealHOViscous(void);
+		ElementMatrix* CreateKMatrixCouplingMacAyealHOFriction(void);
+		ElementMatrix* CreateKMatrixCouplingMacAyealFS(void);
+		ElementMatrix* CreateKMatrixCouplingMacAyealFSViscous(void);
+		ElementMatrix* CreateKMatrixCouplingMacAyealFSFriction(void);
+		ElementMatrix* CreateKMatrixCouplingHOFS(void);
 		ElementMatrix* CreateKMatrixDiagnosticHoriz(void);
 		ElementMatrix* CreateKMatrixAdjointHoriz(void);
 		ElementVector* CreateDVectorDiagnosticHoriz(void);
-		ElementVector* CreateDVectorDiagnosticStokes(void);
-		ElementMatrix* CreateKMatrixDiagnosticHutter(void);
+		ElementVector* CreateDVectorDiagnosticFS(void);
+		ElementMatrix* CreateKMatrixDiagnosticSIA(void);
 		ElementMatrix* CreateKMatrixDiagnosticMacAyeal2d(void);
 		ElementMatrix* CreateKMatrixDiagnosticMacAyeal3d(void);
 		ElementMatrix* CreateKMatrixDiagnosticMacAyeal3dViscous(void);
 		ElementMatrix* CreateKMatrixDiagnosticMacAyeal3dFriction(void);
-		ElementMatrix* CreateKMatrixDiagnosticMacAyealPattyn(void);
-		ElementMatrix* CreateKMatrixDiagnosticMacAyealStokes(void);
+		ElementMatrix* CreateKMatrixDiagnosticMacAyealHO(void);
+		ElementMatrix* CreateKMatrixDiagnosticMacAyealFS(void);
 		ElementMatrix* CreateKMatrixDiagnosticL1L2(void);
 		ElementMatrix* CreateKMatrixDiagnosticL1L2Viscous(void);
 		ElementMatrix* CreateKMatrixDiagnosticL1L2Friction(void);
-		ElementMatrix* CreateKMatrixDiagnosticPattyn(void);
-		ElementMatrix* CreateKMatrixDiagnosticPattynViscous(void);
-		ElementMatrix* CreateKMatrixDiagnosticPattynFriction(void);
-		ElementMatrix* CreateKMatrixDiagnosticPattynStokes(void);
-		ElementMatrix* CreateKMatrixDiagnosticStokes(void);
-		ElementMatrix* CreateKMatrixDiagnosticStokesViscous(void);
-		ElementMatrix* CreateKMatrixDiagnosticStokesGLSViscous(void);
-		ElementMatrix* CreateKMatrixDiagnosticStokesFriction(void);
-		ElementMatrix* CreateKMatrixDiagnosticStokesGLSFriction(void);
+		ElementMatrix* CreateKMatrixDiagnosticHO(void);
+		ElementMatrix* CreateKMatrixDiagnosticHOViscous(void);
+		ElementMatrix* CreateKMatrixDiagnosticHOFriction(void);
+		ElementMatrix* CreateKMatrixDiagnosticHOFS(void);
+		ElementMatrix* CreateKMatrixDiagnosticFS(void);
+		ElementMatrix* CreateKMatrixDiagnosticFSViscous(void);
+		ElementMatrix* CreateKMatrixDiagnosticFSGLSViscous(void);
+		ElementMatrix* CreateKMatrixDiagnosticFSFriction(void);
+		ElementMatrix* CreateKMatrixDiagnosticFSGLSFriction(void);
 		ElementMatrix* CreateKMatrixDiagnosticVert(void);
 		ElementMatrix* CreateKMatrixDiagnosticVertVolume(void);
 		ElementMatrix* CreateKMatrixDiagnosticVertSurface(void);
 		ElementMatrix* CreateJacobianDiagnosticHoriz(void);
-		ElementMatrix* CreateJacobianDiagnosticMacayeal2d(void);
-		ElementMatrix* CreateJacobianDiagnosticPattyn(void);
-		ElementMatrix* CreateJacobianDiagnosticStokes(void);
+		ElementMatrix* CreateJacobianDiagnosticSSA2d(void);
+		ElementMatrix* CreateJacobianDiagnosticHO(void);
+		ElementMatrix* CreateJacobianDiagnosticFS(void);
 		void           InputUpdateFromSolutionDiagnosticHoriz( IssmDouble* solutiong);
 		void           InputUpdateFromSolutionDiagnosticMacAyeal( IssmDouble* solutiong);
-		void           InputUpdateFromSolutionDiagnosticMacAyealPattyn( IssmDouble* solutiong);
-		void           InputUpdateFromSolutionDiagnosticMacAyealStokes( IssmDouble* solutiong);
+		void           InputUpdateFromSolutionDiagnosticMacAyealHO( IssmDouble* solutiong);
+		void           InputUpdateFromSolutionDiagnosticMacAyealFS( IssmDouble* solutiong);
 		void           InputUpdateFromSolutionDiagnosticL1L2( IssmDouble* solutiong);
-		void           InputUpdateFromSolutionDiagnosticPattyn( IssmDouble* solutiong);
-		void           InputUpdateFromSolutionDiagnosticPattynStokes( IssmDouble* solutiong);
-		void           InputUpdateFromSolutionDiagnosticHutter( IssmDouble* solutiong);
+		void           InputUpdateFromSolutionDiagnosticHO( IssmDouble* solutiong);
+		void           InputUpdateFromSolutionDiagnosticHOFS( IssmDouble* solutiong);
+		void           InputUpdateFromSolutionDiagnosticSIA( IssmDouble* solutiong);
 		void           InputUpdateFromSolutionDiagnosticVert( IssmDouble* solutiong);
-		void           InputUpdateFromSolutionDiagnosticStokes( IssmDouble* solutiong);
+		void           InputUpdateFromSolutionDiagnosticFS( IssmDouble* solutiong);
 		void	         GetSolutionFromInputsDiagnosticHoriz(Vector<IssmDouble>* solutiong);
-		void	         GetSolutionFromInputsDiagnosticHutter(Vector<IssmDouble>* solutiong);
-		void	         GetSolutionFromInputsDiagnosticStokes(Vector<IssmDouble>* solutiong);
+		void	         GetSolutionFromInputsDiagnosticSIA(Vector<IssmDouble>* solutiong);
+		void	         GetSolutionFromInputsDiagnosticFS(Vector<IssmDouble>* solutiong);
 		void	         GetSolutionFromInputsDiagnosticVert(Vector<IssmDouble>* solutiong);
-		ElementVector* CreatePVectorCouplingMacAyealStokes(void);
-		ElementVector* CreatePVectorCouplingMacAyealStokesViscous(void);
-		ElementVector* CreatePVectorCouplingMacAyealStokesFriction(void);
-		ElementVector* CreatePVectorCouplingPattynStokes(void);
-		ElementVector* CreatePVectorCouplingPattynStokesViscous(void);
-		ElementVector* CreatePVectorCouplingPattynStokesFriction(void);
+		ElementVector* CreatePVectorCouplingMacAyealFS(void);
+		ElementVector* CreatePVectorCouplingMacAyealFSViscous(void);
+		ElementVector* CreatePVectorCouplingMacAyealFSFriction(void);
+		ElementVector* CreatePVectorCouplingHOFS(void);
+		ElementVector* CreatePVectorCouplingHOFSViscous(void);
+		ElementVector* CreatePVectorCouplingHOFSFriction(void);
 		ElementVector* CreatePVectorDiagnosticHoriz(void);
-		ElementVector* CreatePVectorDiagnosticHutter(void);
+		ElementVector* CreatePVectorDiagnosticSIA(void);
 		ElementVector* CreatePVectorDiagnosticMacAyeal(void);
-		ElementVector* CreatePVectorDiagnosticMacAyealPattyn(void);
-		ElementVector* CreatePVectorDiagnosticMacAyealStokes(void);
+		ElementVector* CreatePVectorDiagnosticMacAyealHO(void);
+		ElementVector* CreatePVectorDiagnosticMacAyealFS(void);
 		ElementVector* CreatePVectorDiagnosticL1L2(void);
-		ElementVector* CreatePVectorDiagnosticPattyn(void);
-		ElementVector* CreatePVectorDiagnosticPattynDrivingStress(void);
-		ElementVector* CreatePVectorDiagnosticPattynFront(void);
-		ElementVector* CreatePVectorDiagnosticPattynStokes(void);
-		ElementVector* CreatePVectorDiagnosticStokes(void);
-		ElementVector* CreatePVectorDiagnosticStokesFront(void);
-		ElementVector* CreatePVectorDiagnosticStokesViscous(void);
-		ElementVector* CreatePVectorDiagnosticStokesGLSViscous(void);
-		ElementVector* CreatePVectorDiagnosticStokesShelf(void);
+		ElementVector* CreatePVectorDiagnosticHO(void);
+		ElementVector* CreatePVectorDiagnosticHODrivingStress(void);
+		ElementVector* CreatePVectorDiagnosticHOFront(void);
+		ElementVector* CreatePVectorDiagnosticHOFS(void);
+		ElementVector* CreatePVectorDiagnosticFS(void);
+		ElementVector* CreatePVectorDiagnosticFSFront(void);
+		ElementVector* CreatePVectorDiagnosticFSViscous(void);
+		ElementVector* CreatePVectorDiagnosticFSGLSViscous(void);
+		ElementVector* CreatePVectorDiagnosticFSShelf(void);
 		ElementVector* CreatePVectorDiagnosticVert(void);
 		ElementVector* CreatePVectorDiagnosticVertVolume(void);
@@ -304,11 +304,11 @@
 		ElementVector* CreatePVectorAdjointHoriz(void);
 		ElementMatrix* CreateKMatrixAdjointMacAyeal2d(void);
-		ElementMatrix* CreateKMatrixAdjointPattyn(void);
-		ElementMatrix* CreateKMatrixAdjointStokes(void);
+		ElementMatrix* CreateKMatrixAdjointHO(void);
+		ElementMatrix* CreateKMatrixAdjointFS(void);
 		ElementVector* CreatePVectorAdjointMacAyeal(void);
-		ElementVector* CreatePVectorAdjointPattyn(void);
-		ElementVector* CreatePVectorAdjointStokes(void);
+		ElementVector* CreatePVectorAdjointHO(void);
+		ElementVector* CreatePVectorAdjointFS(void);
 		void    InputUpdateFromSolutionAdjointHoriz( IssmDouble* solutiong);
-		void    InputUpdateFromSolutionAdjointStokes( IssmDouble* solutiong);
+		void    InputUpdateFromSolutionAdjointFS( IssmDouble* solutiong);
 		#endif
 
Index: issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 15556)
+++ issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 15564)
@@ -54,6 +54,6 @@
 
 /*Reference Element numerics*/
-/*FUNCTION PentaRef::GetBMacAyealPattyn {{{*/
-void PentaRef::GetBMacAyealPattyn(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){
+/*FUNCTION PentaRef::GetBMacAyealHO {{{*/
+void PentaRef::GetBMacAyealHO(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* 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
@@ -85,6 +85,6 @@
 }
 /*}}}*/
-/*FUNCTION PentaRef::GetBMacAyealStokes{{{*/
-void PentaRef::GetBMacAyealStokes(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){
+/*FUNCTION PentaRef::GetBMacAyealFS{{{*/
+void PentaRef::GetBMacAyealFS(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* 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
@@ -131,6 +131,6 @@
 }
 /*}}}*/
-/*FUNCTION PentaRef::GetBPattyn {{{*/
-void PentaRef::GetBPattyn(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){
+/*FUNCTION PentaRef::GetBHO {{{*/
+void PentaRef::GetBHO(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* 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
@@ -171,6 +171,6 @@
 }
 /*}}}*/
-/*FUNCTION PentaRef::GetBprimePattyn {{{*/
-void PentaRef::GetBprimePattyn(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss_coord){
+/*FUNCTION PentaRef::GetBprimeHO {{{*/
+void PentaRef::GetBprimeHO(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss_coord){
 	/*Compute B  prime 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
@@ -209,6 +209,6 @@
 }
 /*}}}*/
-/*FUNCTION PentaRef::GetBprimeMacAyealStokes{{{*/
-void PentaRef::GetBprimeMacAyealStokes(IssmDouble* Bprime, IssmDouble* xyz_list, GaussPenta* gauss){
+/*FUNCTION PentaRef::GetBprimeMacAyealFS{{{*/
+void PentaRef::GetBprimeMacAyealFS(IssmDouble* Bprime, IssmDouble* xyz_list, GaussPenta* 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
@@ -249,6 +249,6 @@
 }
 /*}}}*/
-/*FUNCTION PentaRef::GetBStokes {{{*/
-void PentaRef::GetBStokes(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){
+/*FUNCTION PentaRef::GetBFS {{{*/
+void PentaRef::GetBFS(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){
 
 	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*NDOF4. 
@@ -316,6 +316,6 @@
 }
 /*}}}*/
-/*FUNCTION PentaRef::GetBStokesGLS {{{*/
-void PentaRef::GetBStokesGLS(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){
+/*FUNCTION PentaRef::GetBFSGLS {{{*/
+void PentaRef::GetBFSGLS(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){
 
 	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*NDOF4. 
@@ -382,6 +382,6 @@
 }
 /*}}}*/
-/*FUNCTION PentaRef::GetBprimeStokes {{{*/
-void PentaRef::GetBprimeStokes(IssmDouble* B_prime, IssmDouble* xyz_list, GaussPenta* gauss){
+/*FUNCTION PentaRef::GetBprimeFS {{{*/
+void PentaRef::GetBprimeFS(IssmDouble* B_prime, IssmDouble* xyz_list, GaussPenta* 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
@@ -449,6 +449,6 @@
 }
 /*}}}*/
-/*FUNCTION PentaRef::GetBprimeStokesGLS {{{*/
-void PentaRef::GetBprimeStokesGLS(IssmDouble* B_prime, IssmDouble* xyz_list, GaussPenta* gauss){
+/*FUNCTION PentaRef::GetBprimeFSGLS {{{*/
+void PentaRef::GetBprimeFSGLS(IssmDouble* B_prime, IssmDouble* xyz_list, GaussPenta* 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
@@ -618,6 +618,6 @@
 }
 /*}}}*/
-/*FUNCTION PentaRef::GetBPattynFriction{{{*/
-void PentaRef::GetBPattynFriction(IssmDouble* B, GaussPenta* gauss){
+/*FUNCTION PentaRef::GetBHOFriction{{{*/
+void PentaRef::GetBHOFriction(IssmDouble* B, GaussPenta* gauss){
 	/*Compute B  matrix. B=[B1 B2 B3] where Bi is square and of size 2x2. 
 	 ** For node i, Bi can be expressed in the actual coordinate system
@@ -642,6 +642,6 @@
 } 
 /*}}}*/
-/*FUNCTION PentaRef::GetLStokes{{{*/
-void PentaRef::GetLStokes(IssmDouble* LStokes, GaussPenta* gauss){
+/*FUNCTION PentaRef::GetLFS{{{*/
+void PentaRef::GetLFS(IssmDouble* LFS, GaussPenta* 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
@@ -655,27 +655,27 @@
 
 	const int num_dof=4;
-	IssmDouble l1l2l3[NUMNODESP1_2d];
-
-	/*Get l1l2l3 in actual coordinate system: */
-	l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
-	l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
-	l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
-
-	/*Build LStokes: */
+	IssmDouble L1L2l3[NUMNODESP1_2d];
+
+	/*Get L1L2l3 in actual coordinate system: */
+	L1L2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
+	L1L2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
+	L1L2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
+
+	/*Build LFS: */
 	for(int i=0;i<3;i++){
-		LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = l1l2l3[i];
-		LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;
-		LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+2] = 0.;
-		LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+3] = 0.;
-
-		LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;
-		LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = l1l2l3[i];
-		LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+2] = 0.;
-		LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+3] = 0.;
-	}
-}
-/*}}}*/
-/*FUNCTION PentaRef::GetLprimeStokes {{{*/
-void PentaRef::GetLprimeStokes(IssmDouble* LprimeStokes, IssmDouble* xyz_list, GaussPenta* gauss){
+		LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = L1L2l3[i];
+		LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;
+		LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+2] = 0.;
+		LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+3] = 0.;
+
+		LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;
+		LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = L1L2l3[i];
+		LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+2] = 0.;
+		LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+3] = 0.;
+	}
+}
+/*}}}*/
+/*FUNCTION PentaRef::GetLprimeFS {{{*/
+void PentaRef::GetLprimeFS(IssmDouble* LprimeFS, IssmDouble* xyz_list, GaussPenta* gauss){
 	/* 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
@@ -715,77 +715,77 @@
 	int num_dof=4;
 
-	IssmDouble l1l2l3[NUMNODESP1_2d];
+	IssmDouble L1L2l3[NUMNODESP1_2d];
 	IssmDouble dbasis[3][NUMNODESP1];
 
-	/*Get l1l2l3 in actual coordinate system: */
-	l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
-	l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
-	l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
+	/*Get L1L2l3 in actual coordinate system: */
+	L1L2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
+	L1L2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
+	L1L2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
 
 	GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);
 
-	/*Build LprimeStokes: */
+	/*Build LprimeFS: */
 	for(int i=0;i<3;i++){
-		LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+0]  = l1l2l3[i];
-		LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+1]  = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+2]  = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+3]  = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+0]  = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+1]  = l1l2l3[i];
-		LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+2]  = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+3]  = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+0]  = l1l2l3[i];
-		LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+1]  = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+2]  = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+3]  = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+0]  = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+1]  = l1l2l3[i];
-		LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+2]  = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+3]  = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+0]  = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+1]  = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+2]  = l1l2l3[i];
-		LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+3]  = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+0]  = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+1]  = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+2]  = l1l2l3[i];
-		LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+3]  = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+0]  = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+1]  = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+2]  = dbasis[2][i];
-		LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+3]  = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+0]  = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+1]  = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+2]  = dbasis[2][i];
-		LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+3]  = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*8+num_dof*i+0]  = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*8+num_dof*i+1]  = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*8+num_dof*i+2]  = dbasis[2][i];
-		LprimeStokes[num_dof*NUMNODESP1_2d*8+num_dof*i+3]  = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*9+num_dof*i+0]  = dbasis[2][i];
-		LprimeStokes[num_dof*NUMNODESP1_2d*9+num_dof*i+1]  = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*9+num_dof*i+2]  = dbasis[0][i];
-		LprimeStokes[num_dof*NUMNODESP1_2d*9+num_dof*i+3]  = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*10+num_dof*i+0] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*10+num_dof*i+1] = dbasis[2][i];
-		LprimeStokes[num_dof*NUMNODESP1_2d*10+num_dof*i+2] = dbasis[1][i];
-		LprimeStokes[num_dof*NUMNODESP1_2d*10+num_dof*i+3] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*11+num_dof*i+0] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*11+num_dof*i+1] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*11+num_dof*i+2] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*11+num_dof*i+3] = l1l2l3[i];
-		LprimeStokes[num_dof*NUMNODESP1_2d*12+num_dof*i+0] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*12+num_dof*i+1] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*12+num_dof*i+2] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*12+num_dof*i+3] = l1l2l3[i];
-		LprimeStokes[num_dof*NUMNODESP1_2d*13+num_dof*i+0] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*13+num_dof*i+1] = 0;
-		LprimeStokes[num_dof*NUMNODESP1_2d*13+num_dof*i+2] = 0;
-		LprimeStokes[num_dof*NUMNODESP1_2d*13+num_dof*i+3] = l1l2l3[i];
-	}
-}
-/*}}}*/
-/*FUNCTION PentaRef::GetLMacAyealStokes {{{*/
-void PentaRef::GetLMacAyealStokes(IssmDouble* LStokes, GaussPenta* gauss){
+		LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+0]  = L1L2l3[i];
+		LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+1]  = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+2]  = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+3]  = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+0]  = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+1]  = L1L2l3[i];
+		LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+2]  = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+3]  = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+0]  = L1L2l3[i];
+		LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+1]  = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+2]  = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+3]  = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+0]  = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+1]  = L1L2l3[i];
+		LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+2]  = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+3]  = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*4+num_dof*i+0]  = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*4+num_dof*i+1]  = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*4+num_dof*i+2]  = L1L2l3[i];
+		LprimeFS[num_dof*NUMNODESP1_2d*4+num_dof*i+3]  = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*5+num_dof*i+0]  = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*5+num_dof*i+1]  = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*5+num_dof*i+2]  = L1L2l3[i];
+		LprimeFS[num_dof*NUMNODESP1_2d*5+num_dof*i+3]  = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*6+num_dof*i+0]  = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*6+num_dof*i+1]  = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*6+num_dof*i+2]  = dbasis[2][i];
+		LprimeFS[num_dof*NUMNODESP1_2d*6+num_dof*i+3]  = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*7+num_dof*i+0]  = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*7+num_dof*i+1]  = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*7+num_dof*i+2]  = dbasis[2][i];
+		LprimeFS[num_dof*NUMNODESP1_2d*7+num_dof*i+3]  = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*8+num_dof*i+0]  = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*8+num_dof*i+1]  = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*8+num_dof*i+2]  = dbasis[2][i];
+		LprimeFS[num_dof*NUMNODESP1_2d*8+num_dof*i+3]  = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*9+num_dof*i+0]  = dbasis[2][i];
+		LprimeFS[num_dof*NUMNODESP1_2d*9+num_dof*i+1]  = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*9+num_dof*i+2]  = dbasis[0][i];
+		LprimeFS[num_dof*NUMNODESP1_2d*9+num_dof*i+3]  = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*10+num_dof*i+0] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*10+num_dof*i+1] = dbasis[2][i];
+		LprimeFS[num_dof*NUMNODESP1_2d*10+num_dof*i+2] = dbasis[1][i];
+		LprimeFS[num_dof*NUMNODESP1_2d*10+num_dof*i+3] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*11+num_dof*i+0] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*11+num_dof*i+1] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*11+num_dof*i+2] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*11+num_dof*i+3] = L1L2l3[i];
+		LprimeFS[num_dof*NUMNODESP1_2d*12+num_dof*i+0] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*12+num_dof*i+1] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*12+num_dof*i+2] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*12+num_dof*i+3] = L1L2l3[i];
+		LprimeFS[num_dof*NUMNODESP1_2d*13+num_dof*i+0] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*13+num_dof*i+1] = 0;
+		LprimeFS[num_dof*NUMNODESP1_2d*13+num_dof*i+2] = 0;
+		LprimeFS[num_dof*NUMNODESP1_2d*13+num_dof*i+3] = L1L2l3[i];
+	}
+}
+/*}}}*/
+/*FUNCTION PentaRef::GetLMacAyealFS {{{*/
+void PentaRef::GetLMacAyealFS(IssmDouble* LFS, GaussPenta* gauss){
 	/*
 	 * Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof. 
@@ -804,34 +804,34 @@
 
 	int num_dof=2;
-	IssmDouble l1l2l3[NUMNODESP1_2d];
-
-	/*Get l1l2l3 in actual coordinate system: */
-	l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
-	l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
-	l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
-
-	/*Build LStokes: */
+	IssmDouble L1L2l3[NUMNODESP1_2d];
+
+	/*Get L1L2l3 in actual coordinate system: */
+	L1L2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
+	L1L2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
+	L1L2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
+
+	/*Build LFS: */
 	for(int i=0;i<3;i++){
-		LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = l1l2l3[i];
-		LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0;
-		LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0;
-		LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = l1l2l3[i];
-		LStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = l1l2l3[i];
-		LStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0;
-		LStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0;
-		LStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = l1l2l3[i];
-		LStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+0] = l1l2l3[i];
-		LStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+1] = 0;
-		LStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+0] = 0;
-		LStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+1] = l1l2l3[i];
-		LStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+0] = l1l2l3[i];
-		LStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+1] = 0;
-		LStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+0] = 0;
-		LStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+1] = l1l2l3[i];
-	}
-}
-/*}}}*/
-/*FUNCTION PentaRef::GetLprimeMacAyealStokes {{{*/
-void PentaRef::GetLprimeMacAyealStokes(IssmDouble* LprimeStokes, IssmDouble* xyz_list, GaussPenta* gauss){
+		LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = L1L2l3[i];
+		LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0;
+		LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0;
+		LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = L1L2l3[i];
+		LFS[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = L1L2l3[i];
+		LFS[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0;
+		LFS[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0;
+		LFS[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = L1L2l3[i];
+		LFS[num_dof*NUMNODESP1_2d*4+num_dof*i+0] = L1L2l3[i];
+		LFS[num_dof*NUMNODESP1_2d*4+num_dof*i+1] = 0;
+		LFS[num_dof*NUMNODESP1_2d*5+num_dof*i+0] = 0;
+		LFS[num_dof*NUMNODESP1_2d*5+num_dof*i+1] = L1L2l3[i];
+		LFS[num_dof*NUMNODESP1_2d*6+num_dof*i+0] = L1L2l3[i];
+		LFS[num_dof*NUMNODESP1_2d*6+num_dof*i+1] = 0;
+		LFS[num_dof*NUMNODESP1_2d*7+num_dof*i+0] = 0;
+		LFS[num_dof*NUMNODESP1_2d*7+num_dof*i+1] = L1L2l3[i];
+	}
+}
+/*}}}*/
+/*FUNCTION PentaRef::GetLprimeMacAyealFS {{{*/
+void PentaRef::GetLprimeMacAyealFS(IssmDouble* LprimeFS, IssmDouble* xyz_list, GaussPenta* gauss){
 	/* 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
@@ -848,53 +848,53 @@
 	 */
 	int num_dof=4;
-	IssmDouble l1l2l3[NUMNODESP1_2d];
+	IssmDouble L1L2l3[NUMNODESP1_2d];
 	IssmDouble dbasis[3][NUMNODESP1];
 
-	/*Get l1l2l3 in actual coordinate system: */
-	l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
-	l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
-	l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
+	/*Get L1L2l3 in actual coordinate system: */
+	L1L2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
+	L1L2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
+	L1L2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
 
 	GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);
 
-	/*Build LprimeStokes: */
+	/*Build LprimeFS: */
 	for(int i=0;i<3;i++){
-		LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = l1l2l3[i];
-		LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+2] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+3] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = l1l2l3[i];
-		LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+2] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+3] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+2] = l1l2l3[i];
-		LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+3] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+2] = l1l2l3[i];
-		LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+3] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+0] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+1] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+2] = dbasis[2][i];
-		LprimeStokes[num_dof*NUMNODESP1_2d*4+num_dof*i+3] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+0] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+1] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+2] = dbasis[2][i];
-		LprimeStokes[num_dof*NUMNODESP1_2d*5+num_dof*i+3] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+0] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+1] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+2] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*6+num_dof*i+3] = l1l2l3[i];
-		LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+0] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+1] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+2] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*7+num_dof*i+3] = l1l2l3[i];
-	}
-}
-/*}}}*/
-/*FUNCTION PentaRef::GetLStokesMacAyeal {{{*/
-void PentaRef::GetLStokesMacAyeal(IssmDouble* LStokes, GaussPenta* gauss){
+		LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = L1L2l3[i];
+		LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+2] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+3] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = L1L2l3[i];
+		LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+2] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+3] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+2] = L1L2l3[i];
+		LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+3] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+2] = L1L2l3[i];
+		LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+3] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*4+num_dof*i+0] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*4+num_dof*i+1] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*4+num_dof*i+2] = dbasis[2][i];
+		LprimeFS[num_dof*NUMNODESP1_2d*4+num_dof*i+3] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*5+num_dof*i+0] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*5+num_dof*i+1] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*5+num_dof*i+2] = dbasis[2][i];
+		LprimeFS[num_dof*NUMNODESP1_2d*5+num_dof*i+3] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*6+num_dof*i+0] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*6+num_dof*i+1] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*6+num_dof*i+2] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*6+num_dof*i+3] = L1L2l3[i];
+		LprimeFS[num_dof*NUMNODESP1_2d*7+num_dof*i+0] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*7+num_dof*i+1] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*7+num_dof*i+2] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*7+num_dof*i+3] = L1L2l3[i];
+	}
+}
+/*}}}*/
+/*FUNCTION PentaRef::GetLFSMacAyeal {{{*/
+void PentaRef::GetLFSMacAyeal(IssmDouble* LFS, GaussPenta* 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
@@ -908,34 +908,34 @@
 
 	int num_dof=4;
-	IssmDouble l1l2l3[NUMNODESP1_2d];
-
-	/*Get l1l2l3 in actual coordinate system: */
-	l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
-	l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
-	l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
-
-	/*Build LStokes: */
+	IssmDouble L1L2l3[NUMNODESP1_2d];
+
+	/*Get L1L2l3 in actual coordinate system: */
+	L1L2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
+	L1L2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
+	L1L2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
+
+	/*Build LFS: */
 	for(int i=0;i<3;i++){
-		LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = l1l2l3[i];
-		LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;
-		LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+2] = 0.;
-		LStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+3] = 0.;
-		LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;
-		LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = l1l2l3[i];
-		LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+2] = 0.;
-		LStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+3] = 0.;
-		LStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = 0.;
-		LStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0.;
-		LStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+2] = l1l2l3[i];
-		LStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+3] = 0.;
-		LStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0.;
-		LStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = 0.;
-		LStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+2] = l1l2l3[i];
-		LStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+3] = 0.;
-	}
-}
-/*}}}*/
-/*FUNCTION PentaRef::GetLprimeStokesMacAyeal {{{*/
-void PentaRef::GetLprimeStokesMacAyeal(IssmDouble* LprimeStokes, IssmDouble* xyz_list, GaussPenta* gauss){
+		LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = L1L2l3[i];
+		LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;
+		LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+2] = 0.;
+		LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+3] = 0.;
+		LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;
+		LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = L1L2l3[i];
+		LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+2] = 0.;
+		LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+3] = 0.;
+		LFS[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = 0.;
+		LFS[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0.;
+		LFS[num_dof*NUMNODESP1_2d*2+num_dof*i+2] = L1L2l3[i];
+		LFS[num_dof*NUMNODESP1_2d*2+num_dof*i+3] = 0.;
+		LFS[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0.;
+		LFS[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = 0.;
+		LFS[num_dof*NUMNODESP1_2d*3+num_dof*i+2] = L1L2l3[i];
+		LFS[num_dof*NUMNODESP1_2d*3+num_dof*i+3] = 0.;
+	}
+}
+/*}}}*/
+/*FUNCTION PentaRef::GetLprimeFSMacAyeal {{{*/
+void PentaRef::GetLprimeFSMacAyeal(IssmDouble* LprimeFS, IssmDouble* xyz_list, GaussPenta* gauss){
 	/* 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
@@ -948,23 +948,23 @@
 	 */
 	int num_dof=2;
-	IssmDouble l1l2l3[NUMNODESP1_2d];
+	IssmDouble L1L2l3[NUMNODESP1_2d];
 	IssmDouble dbasis[3][NUMNODESP1];
 
-	/*Get l1l2l3 in actual coordinate system: */
-	l1l2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
-	l1l2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
-	l1l2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
+	/*Get L1L2l3 in actual coordinate system: */
+	L1L2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
+	L1L2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
+	L1L2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
 	GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);
 
-	/*Build LprimeStokes: */
+	/*Build LprimeFS: */
 	for(int i=0;i<3;i++){
-		LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = l1l2l3[i];
-		LprimeStokes[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = l1l2l3[i];
-		LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = l1l2l3[i];
-		LprimeStokes[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0.;
-		LprimeStokes[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = l1l2l3[i];
+		LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = L1L2l3[i];
+		LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = L1L2l3[i];
+		LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = L1L2l3[i];
+		LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0.;
+		LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = L1L2l3[i];
 	}
 }
Index: issm/trunk-jpl/src/c/classes/Elements/PentaRef.h
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/PentaRef.h	(revision 15556)
+++ issm/trunk-jpl/src/c/classes/Elements/PentaRef.h	(revision 15564)
@@ -38,13 +38,13 @@
 		void GetSegmentJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussPenta* gauss);
 		void GetJacobianInvert(IssmDouble*  Jinv, IssmDouble* xyz_list,GaussPenta* gauss);
-		void GetBMacAyealPattyn(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
-		void GetBMacAyealStokes(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
-		void GetBPattyn(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
-		void GetBStokes(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
-		void GetBStokesGLS(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
-		void GetBprimeMacAyealStokes(IssmDouble* Bprime, IssmDouble* xyz_list, GaussPenta* gauss);
-		void GetBprimePattyn(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
-		void GetBprimeStokes(IssmDouble* B_prime, IssmDouble* xyz_list, GaussPenta* gauss);
-		void GetBprimeStokesGLS(IssmDouble* B_prime, IssmDouble* xyz_list, GaussPenta* gauss);
+		void GetBMacAyealHO(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
+		void GetBMacAyealFS(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
+		void GetBHO(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
+		void GetBFS(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
+		void GetBFSGLS(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
+		void GetBprimeMacAyealFS(IssmDouble* Bprime, IssmDouble* xyz_list, GaussPenta* gauss);
+		void GetBprimeHO(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
+		void GetBprimeFS(IssmDouble* B_prime, IssmDouble* xyz_list, GaussPenta* gauss);
+		void GetBprimeFSGLS(IssmDouble* B_prime, IssmDouble* xyz_list, GaussPenta* gauss);
 		void GetBprimeVert(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
 		void GetBAdvec(IssmDouble* B_advec, IssmDouble* xyz_list, GaussPenta* gauss);
@@ -52,11 +52,11 @@
 		void GetBVert(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss);
 		void GetBprimeAdvec(IssmDouble* Bprime_advec, IssmDouble* xyz_list, GaussPenta* gauss);
-		void GetBPattynFriction(IssmDouble* L, GaussPenta* gauss);
-		void GetLStokes(IssmDouble* LStokes, GaussPenta* gauss);
-		void GetLprimeStokes(IssmDouble* LprimeStokes, IssmDouble* xyz_list, GaussPenta* gauss);
-		void GetLMacAyealStokes(IssmDouble* LMacAyealStokes, GaussPenta* gauss);
-		void GetLprimeMacAyealStokes(IssmDouble* LprimeMacAyealStokes, IssmDouble* xyz_list, GaussPenta* gauss);
-		void GetLStokesMacAyeal(IssmDouble* LStokesMacAyeal, GaussPenta* gauss);
-		void GetLprimeStokesMacAyeal(IssmDouble* LprimeStokesMacAyeal, IssmDouble* xyz_list, GaussPenta* gauss);
+		void GetBHOFriction(IssmDouble* L, GaussPenta* gauss);
+		void GetLFS(IssmDouble* LFS, GaussPenta* gauss);
+		void GetLprimeFS(IssmDouble* LprimeFS, IssmDouble* xyz_list, GaussPenta* gauss);
+		void GetLMacAyealFS(IssmDouble* LMacAyealFS, GaussPenta* gauss);
+		void GetLprimeMacAyealFS(IssmDouble* LprimeMacAyealFS, IssmDouble* xyz_list, GaussPenta* gauss);
+		void GetLFSMacAyeal(IssmDouble* LFSMacAyeal, GaussPenta* gauss);
+		void GetLprimeFSMacAyeal(IssmDouble* LprimeFSMacAyeal, IssmDouble* xyz_list, GaussPenta* gauss);
 		void GetInputValue(IssmDouble* pvalue,IssmDouble* plist, GaussPenta* gauss);
 		void GetInputValue(IssmDouble* pvalue,IssmDouble* plist,GaussTria* gauss){_error_("only PentaGauss are supported");};
Index: issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15556)
+++ issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15564)
@@ -193,6 +193,6 @@
 			Ke=CreateKMatrixDiagnosticMacAyeal();
 			break;
-		case DiagnosticHutterAnalysisEnum:
-			Ke=CreateKMatrixDiagnosticHutter();
+		case DiagnosticSIAAnalysisEnum:
+			Ke=CreateKMatrixDiagnosticSIA();
 			break;
 		 #endif
@@ -306,6 +306,6 @@
 			pe=CreatePVectorDiagnosticMacAyeal();
 			break;
-		case DiagnosticHutterAnalysisEnum:
-			pe=CreatePVectorDiagnosticHutter();
+		case DiagnosticSIAAnalysisEnum:
+			pe=CreatePVectorDiagnosticSIA();
 			break;
 		#endif
@@ -423,5 +423,5 @@
 #ifdef _HAVE_DIAGNOSTIC_
 		case DiagnosticHorizAnalysisEnum:
-			Ke=CreateJacobianDiagnosticMacayeal();
+			Ke=CreateJacobianDiagnosticSSA();
 			break;
 #endif
@@ -1201,6 +1201,6 @@
 		GetSolutionFromInputsDiagnosticHoriz(solution);
 		break;
-	case DiagnosticHutterAnalysisEnum:
-		GetSolutionFromInputsDiagnosticHutter(solution);
+	case DiagnosticSIAAnalysisEnum:
+		GetSolutionFromInputsDiagnosticSIA(solution);
 		break;
 	#endif
@@ -1522,5 +1522,5 @@
 			InputUpdateFromSolutionDiagnosticHoriz(solution);
 			break;
-		case DiagnosticHutterAnalysisEnum:
+		case DiagnosticSIAAnalysisEnum:
 			InputUpdateFromSolutionDiagnosticHoriz(solution);
 			break;
@@ -2343,5 +2343,5 @@
 	iomodel->Constant(&progstabilization,PrognosticStabilizationEnum);
 	iomodel->Constant(&balancestabilization,BalancethicknessStabilizationEnum);
-	iomodel->Constant(&fe_ssa,FlowequationFeSsaEnum);
+	iomodel->Constant(&fe_ssa,FlowequationFeSSAEnum);
 	iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum);
 
@@ -3055,6 +3055,6 @@
 }
 /*}}}*/
-/*FUNCTION Tria::CreateKMatrixDiagnosticHutter{{{*/
-ElementMatrix* Tria::CreateKMatrixDiagnosticHutter(void){
+/*FUNCTION Tria::CreateKMatrixDiagnosticSIA{{{*/
+ElementMatrix* Tria::CreateKMatrixDiagnosticSIA(void){
 
 	/*Intermediaries*/
@@ -3224,6 +3224,6 @@
 }
 /*}}}*/
-/*FUNCTION Tria::CreatePVectorDiagnosticHutter{{{*/
-ElementVector* Tria::CreatePVectorDiagnosticHutter(void){
+/*FUNCTION Tria::CreatePVectorDiagnosticSIA{{{*/
+ElementVector* Tria::CreatePVectorDiagnosticSIA(void){
 
 	/*Intermediaries */
@@ -3274,6 +3274,6 @@
 }
 /*}}}*/
-/*FUNCTION Tria::CreateJacobianDiagnosticMacayeal{{{*/
-ElementMatrix* Tria::CreateJacobianDiagnosticMacayeal(void){
+/*FUNCTION Tria::CreateJacobianDiagnosticSSA{{{*/
+ElementMatrix* Tria::CreateJacobianDiagnosticSSA(void){
 
 	/*Intermediaries */
@@ -3382,6 +3382,6 @@
 }
 /*}}}*/
-/*FUNCTION Tria::GetSolutionFromInputsDiagnosticHutter{{{*/
-void  Tria::GetSolutionFromInputsDiagnosticHutter(Vector<IssmDouble>* solution){
+/*FUNCTION Tria::GetSolutionFromInputsDiagnosticSIA{{{*/
+void  Tria::GetSolutionFromInputsDiagnosticSIA(Vector<IssmDouble>* solution){
 
 	const int    numdof=NDOF2*NUMVERTICES;
@@ -3494,6 +3494,6 @@
 }
 /*}}}*/
-/*FUNCTION Tria::InputUpdateFromSolutionDiagnosticHutter {{{*/
-void  Tria::InputUpdateFromSolutionDiagnosticHutter(IssmDouble* solution){
+/*FUNCTION Tria::InputUpdateFromSolutionDiagnosticSIA {{{*/
+void  Tria::InputUpdateFromSolutionDiagnosticSIA(IssmDouble* solution){
 
 	int        i;
@@ -5137,6 +5137,6 @@
 }
 /*}}}*/
-/*FUNCTION Tria::CreatePVectorAdjointStokes{{{*/
-ElementVector* Tria::CreatePVectorAdjointStokes(void){
+/*FUNCTION Tria::CreatePVectorAdjointFS{{{*/
+ElementVector* Tria::CreatePVectorAdjointFS(void){
 
 	/*Intermediaries */
@@ -5156,5 +5156,5 @@
 
 	/*Initialize Element vector*/
-	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);
+	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,FSApproximationEnum);
 
 	/*Retrieve all inputs and parameters*/
Index: issm/trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 15556)
+++ issm/trunk-jpl/src/c/classes/Elements/Tria.h	(revision 15564)
@@ -147,5 +147,5 @@
 		void   GradjZMacAyeal(Vector<IssmDouble>* gradient,int control_index);
 		void   GradjDragMacAyeal(Vector<IssmDouble>* gradient,int control_index);
-		void   GradjDragStokes(Vector<IssmDouble>* gradient,int control_index);
+		void   GradjDragFS(Vector<IssmDouble>* gradient,int control_index);
 		void   GradjDragGradient(Vector<IssmDouble>* gradient,int weight_index,int control_index);
 		void   GradjDhDtBalancedthickness(Vector<IssmDouble>* gradient,int control_index);
@@ -223,14 +223,14 @@
 		ElementMatrix* CreateKMatrixDiagnosticMacAyealViscous(void);
 		ElementMatrix* CreateKMatrixDiagnosticMacAyealFriction(void);
-		ElementMatrix* CreateKMatrixDiagnosticHutter(void);
+		ElementMatrix* CreateKMatrixDiagnosticSIA(void);
 		ElementVector* CreatePVectorDiagnosticMacAyeal(void);
 		ElementVector* CreatePVectorDiagnosticMacAyealDrivingStress(void);
 		ElementVector* CreatePVectorDiagnosticMacAyealFront(void);
-		ElementVector* CreatePVectorDiagnosticHutter(void);
-		ElementMatrix* CreateJacobianDiagnosticMacayeal(void);
+		ElementVector* CreatePVectorDiagnosticSIA(void);
+		ElementMatrix* CreateJacobianDiagnosticSSA(void);
 		void	  GetSolutionFromInputsDiagnosticHoriz(Vector<IssmDouble>* solution);
-		void	  GetSolutionFromInputsDiagnosticHutter(Vector<IssmDouble>* solution);
+		void	  GetSolutionFromInputsDiagnosticSIA(Vector<IssmDouble>* solution);
 		void	  InputUpdateFromSolutionDiagnosticHoriz( IssmDouble* solution);
-		void	  InputUpdateFromSolutionDiagnosticHutter( IssmDouble* solution);
+		void	  InputUpdateFromSolutionDiagnosticSIA( IssmDouble* solution);
 		#endif
 
@@ -239,5 +239,5 @@
 		ElementMatrix* CreateKMatrixAdjointMacAyeal(void);
 		ElementVector* CreatePVectorAdjointHoriz(void);
-		ElementVector* CreatePVectorAdjointStokes(void);
+		ElementVector* CreatePVectorAdjointFS(void);
 		ElementVector* CreatePVectorAdjointBalancethickness(void);
 		void	  InputUpdateFromSolutionAdjointBalancethickness( IssmDouble* solution);
Index: issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 15556)
+++ issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 15564)
@@ -114,6 +114,6 @@
 }
 /*}}}*/
-/*FUNCTION TriaRef::GetBMacAyealStokes {{{*/
-void TriaRef::GetBMacAyealStokes(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss){
+/*FUNCTION TriaRef::GetBMacAyealFS {{{*/
+void TriaRef::GetBMacAyealFS(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss){
 
 	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
@@ -265,6 +265,6 @@
 }
 /*}}}*/
-/*FUNCTION TriaRef::GetBprimeMacAyealStokes {{{*/
-void TriaRef::GetBprimeMacAyealStokes(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss){
+/*FUNCTION TriaRef::GetBprimeMacAyealFS {{{*/
+void TriaRef::GetBprimeMacAyealFS(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* 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
Index: issm/trunk-jpl/src/c/classes/Elements/TriaRef.h
===================================================================
--- issm/trunk-jpl/src/c/classes/Elements/TriaRef.h	(revision 15556)
+++ issm/trunk-jpl/src/c/classes/Elements/TriaRef.h	(revision 15564)
@@ -24,7 +24,7 @@
 		/*Numerics*/
 		void GetBMacAyeal(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss);
-		void GetBMacAyealStokes(IssmDouble* B , IssmDouble* xyz_list, GaussTria* gauss);
+		void GetBMacAyealFS(IssmDouble* B , IssmDouble* xyz_list, GaussTria* gauss);
 		void GetBprimeMacAyeal(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss);
-		void GetBprimeMacAyealStokes(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss);
+		void GetBprimeMacAyealFS(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss);
 		void GetBprimePrognostic(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss);
 		void GetBPrognostic(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss);
