Index: ../trunk-jpl/src/c/classes/Elements/Element.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 18919)
+++ ../trunk-jpl/src/c/classes/Elements/Element.cpp	(revision 18920)
@@ -305,9 +305,6 @@
 	*ptransform=transform;
 }
 /*}}}*/
-void       Element::DeleteMaterials(void){/*{{{*/
-	delete this->material;
-}/*}}}*/
 void       Element::DeepEcho(void){/*{{{*/
 
 	_printf_(EnumToStringx(this->ObjectEnum())<<" element:\n");
@@ -342,40 +339,15 @@
 	return;
 }
 /*}}}*/
-void       Element::Echo(void){/*{{{*/
-	_printf_(EnumToStringx(this->ObjectEnum())<<" element:\n");
-	_printf_("   id : "<<this->id <<"\n");
-	_printf_("   sid: "<<this->sid<<"\n");
-	if(vertices){
-		int numvertices = this->GetNumberOfVertices();
-		for(int i=0;i<numvertices;i++) vertices[i]->Echo();
-	}
-	else _printf_("vertices = NULL\n");
+void       Element::DeleteInput(int input_enum){/*{{{*/
 
-	if(nodes){
-		int numnodes = this->GetNumberOfNodes();
-		for(int i=0;i<numnodes;i++) {
-			_printf_("nodes[" << i << "] = " << nodes[i]);	
-			nodes[i]->Echo();
-		}
-	}
-	else _printf_("nodes = NULL\n");
+	inputs->DeleteInput(input_enum);
 
-	if (material) material->Echo();
-	else _printf_("material = NULL\n");
-
-	if (matpar) matpar->Echo();
-	else _printf_("matpar = NULL\n");
-
-	_printf_("   parameters\n");
-	if (parameters) parameters->Echo();
-	else _printf_("parameters = NULL\n");
-
-	_printf_("   inputs\n");
-	if (inputs) inputs->Echo();
-	else _printf_("inputs=NULL\n");
 }
 /*}}}*/
+void       Element::DeleteMaterials(void){/*{{{*/
+	delete this->material;
+}/*}}}*/
 IssmDouble Element::Divergence(void){/*{{{*/
 	/*Compute element divergence*/
 
@@ -418,25 +390,25 @@
 	delete gauss;
 	return divergence;
 }/*}}}*/
-void       Element::dViscositydBSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
+void       Element::dViscositydBFS(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/
 
 	/*Intermediaries*/
 	IssmDouble dmudB;
+	IssmDouble epsilon3d[6];/* epsilon=[exx,eyy,exy,exy,exz,eyz];    */
 	IssmDouble epsilon2d[3];/* epsilon=[exx,eyy,exy];    */
-	IssmDouble epsilon1d;   /* epsilon=[exx];    */
 	IssmDouble eps_eff;
+	IssmDouble eps0=1.e-27;
 
-	 if(dim==2){
-		 /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy*/
-		 this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
-		 eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + epsilon2d[2]*epsilon2d[2] + epsilon2d[0]*epsilon2d[1]);
-	 }
-	 else{
-		 /* eps_eff^2 = 1/2 exx^2*/
-		 this->StrainRateSSA1d(&epsilon1d,xyz_list,gauss,vx_input);
-		 eps_eff = sqrt(epsilon1d*epsilon1d/2.);
-	 }
-
+	if(dim==3){
+		/* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */
+		this->StrainRateFS(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input,vz_input);
+		eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[3]*epsilon3d[3] +  epsilon3d[4]*epsilon3d[4] + epsilon3d[5]*epsilon3d[5] + epsilon3d[0]*epsilon3d[1]+eps0*eps0);
+	}
+	else{
+		/* eps_eff^2 = 1/2 ( exx^2 + eyy^2 + 2*exy^2 )*/
+		this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
+		eps_eff = 1./sqrt(2.)*sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + 2.*epsilon2d[2]*epsilon2d[2]);
+	}
 	/*Get viscosity*/
 	material->GetViscosity_B(&dmudB,eps_eff);
 
@@ -472,25 +444,25 @@
 
 }
 /*}}}*/
-void       Element::dViscositydBFS(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/
+void       Element::dViscositydBSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
 
 	/*Intermediaries*/
 	IssmDouble dmudB;
-	IssmDouble epsilon3d[6];/* epsilon=[exx,eyy,exy,exy,exz,eyz];    */
 	IssmDouble epsilon2d[3];/* epsilon=[exx,eyy,exy];    */
+	IssmDouble epsilon1d;   /* epsilon=[exx];    */
 	IssmDouble eps_eff;
-	IssmDouble eps0=1.e-27;
 
-	if(dim==3){
-		/* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */
-		this->StrainRateFS(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input,vz_input);
-		eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[3]*epsilon3d[3] +  epsilon3d[4]*epsilon3d[4] + epsilon3d[5]*epsilon3d[5] + epsilon3d[0]*epsilon3d[1]+eps0*eps0);
-	}
-	else{
-		/* eps_eff^2 = 1/2 ( exx^2 + eyy^2 + 2*exy^2 )*/
-		this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
-		eps_eff = 1./sqrt(2.)*sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + 2.*epsilon2d[2]*epsilon2d[2]);
-	}
+	 if(dim==2){
+		 /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy*/
+		 this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
+		 eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + epsilon2d[2]*epsilon2d[2] + epsilon2d[0]*epsilon2d[1]);
+	 }
+	 else{
+		 /* eps_eff^2 = 1/2 exx^2*/
+		 this->StrainRateSSA1d(&epsilon1d,xyz_list,gauss,vx_input);
+		 eps_eff = sqrt(epsilon1d*epsilon1d/2.);
+	 }
+
 	/*Get viscosity*/
 	material->GetViscosity_B(&dmudB,eps_eff);
 
@@ -526,18 +498,49 @@
 
 }
 /*}}}*/
-void       Element::ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){/*{{{*/
-	matpar->ThermalToEnthalpy(penthalpy,temperature,waterfraction,pressure);
-}/*}}}*/
-void       Element::EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){/*{{{*/
-	matpar->EnthalpyToThermal(ptemperature,pwaterfraction,enthalpy,pressure);
-}/*}}}*/
+void       Element::Echo(void){/*{{{*/
+	_printf_(EnumToStringx(this->ObjectEnum())<<" element:\n");
+	_printf_("   id : "<<this->id <<"\n");
+	_printf_("   sid: "<<this->sid<<"\n");
+	if(vertices){
+		int numvertices = this->GetNumberOfVertices();
+		for(int i=0;i<numvertices;i++) vertices[i]->Echo();
+	}
+	else _printf_("vertices = NULL\n");
+
+	if(nodes){
+		int numnodes = this->GetNumberOfNodes();
+		for(int i=0;i<numnodes;i++) {
+			_printf_("nodes[" << i << "] = " << nodes[i]);	
+			nodes[i]->Echo();
+		}
+	}
+	else _printf_("nodes = NULL\n");
+
+	if (material) material->Echo();
+	else _printf_("material = NULL\n");
+
+	if (matpar) matpar->Echo();
+	else _printf_("matpar = NULL\n");
+
+	_printf_("   parameters\n");
+	if (parameters) parameters->Echo();
+	else _printf_("parameters = NULL\n");
+
+	_printf_("   inputs\n");
+	if (inputs) inputs->Echo();
+	else _printf_("inputs=NULL\n");
+}
+/*}}}*/
 IssmDouble Element::EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){/*{{{*/
 	return matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure);
 }/*}}}*/
 IssmDouble Element::EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){/*{{{*/
 	return matpar->GetEnthalpyDiffusionParameterVolume(numvertices,enthalpy,pressure);
 }/*}}}*/
+void       Element::EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){/*{{{*/
+	matpar->EnthalpyToThermal(ptemperature,pwaterfraction,enthalpy,pressure);
+}/*}}}*/
 void       Element::FindParam(bool* pvalue,int paramenum){/*{{{*/
 	this->parameters->FindParam(pvalue,paramenum);
 }/*}}}*/
@@ -573,155 +576,56 @@
 	*pdoflist=doflist;
 }
 /*}}}*/
-void       Element::GetDofListVelocity(int** pdoflist,int setenum){/*{{{*/
+void       Element::GetDofListPressure(int** pdoflist,int setenum){/*{{{*/
 
 	/*Fetch number of nodes and dof for this finite element*/
-	int numnodes = this->NumberofNodesVelocity();
+	int vnumnodes = this->NumberofNodesVelocity();
+	int pnumnodes = this->NumberofNodesPressure();
 
 	/*First, figure out size of doflist and create it: */
 	int numberofdofs=0;
-	for(int i=0;i<numnodes;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(FSvelocityEnum,setenum);
+	for(int i=vnumnodes;i<vnumnodes+pnumnodes;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(FSApproximationEnum,setenum);
 
 	/*Allocate output*/
 	int* doflist=xNew<int>(numberofdofs);
 
 	/*Populate: */
 	int count=0;
-	for(int i=0;i<numnodes;i++){
-		nodes[i]->GetDofList(doflist+count,FSvelocityEnum,setenum);
-		count+=nodes[i]->GetNumberOfDofs(FSvelocityEnum,setenum);
+	for(int i=vnumnodes;i<vnumnodes+pnumnodes;i++){
+		nodes[i]->GetDofList(doflist+count,FSApproximationEnum,setenum);
+		count+=nodes[i]->GetNumberOfDofs(FSApproximationEnum,setenum);
 	}
 
 	/*Assign output pointers:*/
 	*pdoflist=doflist;
 }
 /*}}}*/
-void       Element::GetDofListPressure(int** pdoflist,int setenum){/*{{{*/
+void       Element::GetDofListVelocity(int** pdoflist,int setenum){/*{{{*/
 
 	/*Fetch number of nodes and dof for this finite element*/
-	int vnumnodes = this->NumberofNodesVelocity();
-	int pnumnodes = this->NumberofNodesPressure();
+	int numnodes = this->NumberofNodesVelocity();
 
 	/*First, figure out size of doflist and create it: */
 	int numberofdofs=0;
-	for(int i=vnumnodes;i<vnumnodes+pnumnodes;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(FSApproximationEnum,setenum);
+	for(int i=0;i<numnodes;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(FSvelocityEnum,setenum);
 
 	/*Allocate output*/
 	int* doflist=xNew<int>(numberofdofs);
 
 	/*Populate: */
 	int count=0;
-	for(int i=vnumnodes;i<vnumnodes+pnumnodes;i++){
-		nodes[i]->GetDofList(doflist+count,FSApproximationEnum,setenum);
-		count+=nodes[i]->GetNumberOfDofs(FSApproximationEnum,setenum);
+	for(int i=0;i<numnodes;i++){
+		nodes[i]->GetDofList(doflist+count,FSvelocityEnum,setenum);
+		count+=nodes[i]->GetNumberOfDofs(FSvelocityEnum,setenum);
 	}
 
 	/*Assign output pointers:*/
 	*pdoflist=doflist;
 }
 /*}}}*/
-IssmDouble Element::GetMaterialParameter(int enum_in){/*{{{*/
-
-	_assert_(this->matpar);
-	switch(enum_in){ // FIXME: change this to material
-		case MaterialsRheologyNEnum:
-			return this->material->GetN();
-		case MaterialsRheologyBEnum:
-			return this->material->GetB();
-		case MaterialsRheologyBbarEnum:
-			return this->material->GetBbar();
-		default:
-			return this->matpar->GetMaterialParameter(enum_in);
-	}
-}/*}}}*/
-void       Element::GetPhi(IssmDouble* phi, IssmDouble*  epsilon, IssmDouble viscosity){/*{{{*/
-	/*Compute deformational heating from epsilon and viscosity */
-
-	IssmDouble epsilon_matrix[3][3];
-	IssmDouble epsilon_eff;
-	IssmDouble epsilon_sqr[3][3];
-
-	/* Build epsilon matrix */
-	epsilon_matrix[0][0]=epsilon[0];
-	epsilon_matrix[1][0]=epsilon[3];
-	epsilon_matrix[2][0]=epsilon[4];
-	epsilon_matrix[0][1]=epsilon[3];
-	epsilon_matrix[1][1]=epsilon[1];
-	epsilon_matrix[2][1]=epsilon[5];
-	epsilon_matrix[0][2]=epsilon[4];
-	epsilon_matrix[1][2]=epsilon[5];
-	epsilon_matrix[2][2]=epsilon[2];
-
-	/* Effective value of epsilon_matrix */
-	epsilon_sqr[0][0]=epsilon_matrix[0][0]*epsilon_matrix[0][0];
-	epsilon_sqr[1][0]=epsilon_matrix[1][0]*epsilon_matrix[1][0];
-	epsilon_sqr[2][0]=epsilon_matrix[2][0]*epsilon_matrix[2][0];
-	epsilon_sqr[0][1]=epsilon_matrix[0][1]*epsilon_matrix[0][1];
-	epsilon_sqr[1][1]=epsilon_matrix[1][1]*epsilon_matrix[1][1];
-	epsilon_sqr[2][1]=epsilon_matrix[2][1]*epsilon_matrix[2][1];
-	epsilon_sqr[0][2]=epsilon_matrix[0][2]*epsilon_matrix[0][2];
-	epsilon_sqr[1][2]=epsilon_matrix[1][2]*epsilon_matrix[1][2];
-	epsilon_sqr[2][2]=epsilon_matrix[2][2]*epsilon_matrix[2][2];
-	epsilon_eff=1/sqrt(2.)*sqrt(epsilon_sqr[0][0]+epsilon_sqr[0][1]+ epsilon_sqr[0][2]+ epsilon_sqr[1][0]+ epsilon_sqr[1][1]+ epsilon_sqr[1][2]+ epsilon_sqr[2][0]+ epsilon_sqr[2][1]+ epsilon_sqr[2][2]);
-
-	/*Phi = Tr(sigma * eps) 
-	 *    = Tr(sigma'* eps)
-	 *    = 2 * eps_eff * sigma'_eff
-	 *    = 4 * mu * eps_eff ^2*/
-	*phi=4.*epsilon_eff*epsilon_eff*viscosity;
-}
-/*}}}*/
 Input*     Element::GetInput(int inputenum){/*{{{*/
 	return inputs->GetInput(inputenum);
 }/*}}}*/
-void       Element::GetInputListOnVertices(IssmDouble* pvalue,int enumtype){/*{{{*/
-
-	/*Recover input*/
-	Input* input=this->GetInput(enumtype);
-	if (!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element");
-
-	/*Fetch number vertices for this element*/
-	int numvertices = this->GetNumberOfVertices();
-
-	/*Checks in debugging mode*/
-	_assert_(pvalue);
-
-	/* Start looping on the number of vertices: */
-	Gauss*gauss=this->NewGauss();
-	for(int iv=0;iv<numvertices;iv++){
-		gauss->GaussVertex(iv);
-		input->GetInputValue(&pvalue[iv],gauss);
-	}
-
-	/*clean-up*/
-	delete gauss;
-}
-/*}}}*/
-void       Element::GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue){/*{{{*/
-
-	/*Recover input*/
-	Input* input=this->GetInput(enumtype);
-
-	/*Checks in debugging mode*/
-	_assert_(pvalue);
-
-	/*Fetch number vertices for this element*/
-	int numvertices = this->GetNumberOfVertices();
-
-	/* Start looping on the number of vertices: */
-	if (input){
-		Gauss* gauss=this->NewGauss();
-		for (int iv=0;iv<numvertices;iv++){
-			gauss->GaussVertex(iv);
-			input->GetInputValue(&pvalue[iv],gauss);
-		}
-		delete gauss;
-	}
-	else{
-		for(int iv=0;iv<numvertices;iv++) pvalue[iv]=defaultvalue;
-	}
-}
-/*}}}*/
 void       Element::GetInputListOnNodes(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue){/*{{{*/
 
 	_assert_(pvalue);
@@ -777,6 +681,54 @@
 	delete gauss;
 }
 /*}}}*/
+void       Element::GetInputListOnVertices(IssmDouble* pvalue,int enumtype){/*{{{*/
+
+	/*Recover input*/
+	Input* input=this->GetInput(enumtype);
+	if (!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element");
+
+	/*Fetch number vertices for this element*/
+	int numvertices = this->GetNumberOfVertices();
+
+	/*Checks in debugging mode*/
+	_assert_(pvalue);
+
+	/* Start looping on the number of vertices: */
+	Gauss*gauss=this->NewGauss();
+	for(int iv=0;iv<numvertices;iv++){
+		gauss->GaussVertex(iv);
+		input->GetInputValue(&pvalue[iv],gauss);
+	}
+
+	/*clean-up*/
+	delete gauss;
+}
+/*}}}*/
+void       Element::GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue){/*{{{*/
+
+	/*Recover input*/
+	Input* input=this->GetInput(enumtype);
+
+	/*Checks in debugging mode*/
+	_assert_(pvalue);
+
+	/*Fetch number vertices for this element*/
+	int numvertices = this->GetNumberOfVertices();
+
+	/* Start looping on the number of vertices: */
+	if (input){
+		Gauss* gauss=this->NewGauss();
+		for (int iv=0;iv<numvertices;iv++){
+			gauss->GaussVertex(iv);
+			input->GetInputValue(&pvalue[iv],gauss);
+		}
+		delete gauss;
+	}
+	else{
+		for(int iv=0;iv<numvertices;iv++) pvalue[iv]=defaultvalue;
+	}
+}
+/*}}}*/
 void       Element::GetInputLocalMinMaxOnNodes(IssmDouble* min,IssmDouble* max,IssmDouble* ug){/*{{{*/
 
 
@@ -832,26 +784,77 @@
 	input->GetInputValue(pvalue,gauss);
 
 }/*}}}*/
-void       Element::GetNodesSidList(int* sidlist){/*{{{*/
+IssmDouble Element::GetMaterialParameter(int enum_in){/*{{{*/
 
-	_assert_(sidlist);
+	_assert_(this->matpar);
+	switch(enum_in){ // FIXME: change this to material
+		case MaterialsRheologyNEnum:
+			return this->material->GetN();
+		case MaterialsRheologyBEnum:
+			return this->material->GetB();
+		case MaterialsRheologyBbarEnum:
+			return this->material->GetBbar();
+		default:
+			return this->matpar->GetMaterialParameter(enum_in);
+	}
+}/*}}}*/
+void       Element::GetNodesLidList(int* lidlist){/*{{{*/
+
+	_assert_(lidlist);
 	_assert_(nodes);
 	int numnodes = this->GetNumberOfNodes();
 	for(int i=0;i<numnodes;i++){
-		sidlist[i]=nodes[i]->Sid();
+		lidlist[i]=nodes[i]->Lid();
 	}
 }
 /*}}}*/
-void       Element::GetNodesLidList(int* lidlist){/*{{{*/
+void       Element::GetNodesSidList(int* sidlist){/*{{{*/
 
-	_assert_(lidlist);
+	_assert_(sidlist);
 	_assert_(nodes);
 	int numnodes = this->GetNumberOfNodes();
 	for(int i=0;i<numnodes;i++){
-		lidlist[i]=nodes[i]->Lid();
+		sidlist[i]=nodes[i]->Sid();
 	}
 }
 /*}}}*/
+void       Element::GetPhi(IssmDouble* phi, IssmDouble*  epsilon, IssmDouble viscosity){/*{{{*/
+	/*Compute deformational heating from epsilon and viscosity */
+
+	IssmDouble epsilon_matrix[3][3];
+	IssmDouble epsilon_eff;
+	IssmDouble epsilon_sqr[3][3];
+
+	/* Build epsilon matrix */
+	epsilon_matrix[0][0]=epsilon[0];
+	epsilon_matrix[1][0]=epsilon[3];
+	epsilon_matrix[2][0]=epsilon[4];
+	epsilon_matrix[0][1]=epsilon[3];
+	epsilon_matrix[1][1]=epsilon[1];
+	epsilon_matrix[2][1]=epsilon[5];
+	epsilon_matrix[0][2]=epsilon[4];
+	epsilon_matrix[1][2]=epsilon[5];
+	epsilon_matrix[2][2]=epsilon[2];
+
+	/* Effective value of epsilon_matrix */
+	epsilon_sqr[0][0]=epsilon_matrix[0][0]*epsilon_matrix[0][0];
+	epsilon_sqr[1][0]=epsilon_matrix[1][0]*epsilon_matrix[1][0];
+	epsilon_sqr[2][0]=epsilon_matrix[2][0]*epsilon_matrix[2][0];
+	epsilon_sqr[0][1]=epsilon_matrix[0][1]*epsilon_matrix[0][1];
+	epsilon_sqr[1][1]=epsilon_matrix[1][1]*epsilon_matrix[1][1];
+	epsilon_sqr[2][1]=epsilon_matrix[2][1]*epsilon_matrix[2][1];
+	epsilon_sqr[0][2]=epsilon_matrix[0][2]*epsilon_matrix[0][2];
+	epsilon_sqr[1][2]=epsilon_matrix[1][2]*epsilon_matrix[1][2];
+	epsilon_sqr[2][2]=epsilon_matrix[2][2]*epsilon_matrix[2][2];
+	epsilon_eff=1/sqrt(2.)*sqrt(epsilon_sqr[0][0]+epsilon_sqr[0][1]+ epsilon_sqr[0][2]+ epsilon_sqr[1][0]+ epsilon_sqr[1][1]+ epsilon_sqr[1][2]+ epsilon_sqr[2][0]+ epsilon_sqr[2][1]+ epsilon_sqr[2][2]);
+
+	/*Phi = Tr(sigma * eps) 
+	 *    = Tr(sigma'* eps)
+	 *    = 2 * eps_eff * sigma'_eff
+	 *    = 4 * mu * eps_eff ^2*/
+	*phi=4.*epsilon_eff*epsilon_eff*viscosity;
+}
+/*}}}*/
 void       Element::GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum){/*{{{*/
 
 	/*Fetch number vertices for this element and allocate arrays*/
@@ -877,6 +880,12 @@
 
 }
 /*}}}*/
+void       Element::GetVerticesConnectivityList(int* connectivity){/*{{{*/
+
+	int numvertices = this->GetNumberOfVertices();
+	for(int i=0;i<numvertices;i++) connectivity[i]=this->vertices[i]->Connectivity();
+}
+/*}}}*/
 void       Element::GetVerticesCoordinates(IssmDouble** pxyz_list){/*{{{*/
 
 	int         numvertices = this->GetNumberOfVertices();
@@ -892,12 +901,6 @@
 	for(int i=0;i<numvertices;i++) sidlist[i]=this->vertices[i]->Sid();
 }
 /*}}}*/
-void       Element::GetVerticesConnectivityList(int* connectivity){/*{{{*/
-
-	int numvertices = this->GetNumberOfVertices();
-	for(int i=0;i<numvertices;i++) connectivity[i]=this->vertices[i]->Connectivity();
-}
-/*}}}*/
 IssmDouble Element::GetXcoord(IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 
 	/*output*/
@@ -1057,12 +1060,6 @@
 
 }
 /*}}}*/
-void       Element::DeleteInput(int input_enum){/*{{{*/
-
-	inputs->DeleteInput(input_enum);
-
-}
-/*}}}*/
 void       Element::InputUpdateFromConstant(int constant, int name){/*{{{*/
 
 	/*Check that name is an element input*/
@@ -1297,10 +1294,6 @@
 
 }
 /*}}}*/
-ElementVector* Element::NewElementVector(int approximation_enum){/*{{{*/
-	return new ElementVector(nodes,this->GetNumberOfNodes(),this->parameters,approximation_enum);
-}
-/*}}}*/
 ElementMatrix* Element::NewElementMatrix(int approximation_enum){/*{{{*/
 	return new ElementMatrix(nodes,this->GetNumberOfNodes(),this->parameters,approximation_enum);
 }
@@ -1309,6 +1302,10 @@
 	return new ElementMatrix(nodes,number_nodes,this->parameters,approximation_enum);
 }
 /*}}}*/
+ElementVector* Element::NewElementVector(int approximation_enum){/*{{{*/
+	return new ElementVector(nodes,this->GetNumberOfNodes(),this->parameters,approximation_enum);
+}
+/*}}}*/
 IssmDouble Element::PureIceEnthalpy(IssmDouble pressure){/*{{{*/
 	return this->matpar->PureIceEnthalpy(pressure);
 }/*}}}*/
@@ -1392,6 +1389,14 @@
 	*pinterpolation   = input->GetResultInterpolation();
 	*pnodesperelement = input->GetResultNumberOfNodes();
 }/*}}}*/
+void       Element::ResultToPatch(IssmDouble* values,int nodesperelement,int output_enum){/*{{{*/
+
+	Input* input=this->inputs->GetInput(output_enum);
+	if(!input) _error_("input "<<EnumToStringx(output_enum)<<" not found in element");
+
+	input->ResultToPatch(values,nodesperelement,this->Sid());
+
+} /*}}}*/
 void       Element::ResultToVector(Vector<IssmDouble>* vector,int output_enum){/*{{{*/
 
 	Input* input=this->inputs->GetInput(output_enum);
@@ -1440,14 +1445,6 @@
 					 _error_("interpolation "<<EnumToStringx(input->GetResultInterpolation())<<" not supported yet");
 	}
 } /*}}}*/
-void       Element::ResultToPatch(IssmDouble* values,int nodesperelement,int output_enum){/*{{{*/
-
-	Input* input=this->inputs->GetInput(output_enum);
-	if(!input) _error_("input "<<EnumToStringx(output_enum)<<" not found in element");
-
-	input->ResultToPatch(values,nodesperelement,this->Sid());
-
-} /*}}}*/
 void       Element::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){/*{{{*/
 
 	/*Intermediaries*/
@@ -1524,295 +1521,6 @@
 
 }
 /*}}}*/
-IssmDouble Element::TMeltingPoint(IssmDouble pressure){/*{{{*/
-	_assert_(matpar);
-	return this->matpar->TMeltingPoint(pressure);
-}/*}}}*/
-void       Element::StressMaxPrincipalCreateInput(void){/*{{{*/
-
-	/*Intermediaries*/
-	IssmDouble *xyz_list = NULL;
-	IssmDouble  sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz;
-	IssmDouble  a,b,c,d,x[3],max;
-	int         dim,numroots;
-
-	/*First: get stress tensor*/
-	this->ComputeStressTensor();
-
-	/*Get domain dimension*/
-	this->FindParam(&dim,DomainDimensionEnum);
-
-	/*Fetch number vertices and allocate memory*/
-	int         numvertices  = this->GetNumberOfVertices();
-	IssmDouble* maxprincipal = xNew<IssmDouble>(numvertices);
-
-	/*Retrieve all inputs and parameters*/
-	this->GetVerticesCoordinatesBase(&xyz_list);
-	Input* sigma_xx_input  = this->GetInput(StressTensorxxEnum); _assert_(sigma_xx_input);
-	Input* sigma_yy_input  = this->GetInput(StressTensoryyEnum); _assert_(sigma_yy_input);
-	Input* sigma_xy_input  = this->GetInput(StressTensorxyEnum); _assert_(sigma_xy_input);
-	Input* sigma_xz_input  = NULL;
-	Input* sigma_yz_input  = NULL;
-	Input* sigma_zz_input  = NULL;
-	if(dim==3){
-		sigma_xz_input  = this->GetInput(StressTensorxzEnum); _assert_(sigma_xz_input);
-		sigma_yz_input  = this->GetInput(StressTensoryzEnum); _assert_(sigma_yz_input);
-		sigma_zz_input  = this->GetInput(StressTensorzzEnum); _assert_(sigma_zz_input);
-	}
-
-	/*loop over vertices: */
-	Gauss* gauss=this->NewGauss();
-	for (int iv=0;iv<numvertices;iv++){
-		gauss->GaussVertex(iv);
-
-		sigma_xx_input->GetInputValue(&sigma_xx,gauss);
-		sigma_yy_input->GetInputValue(&sigma_yy,gauss);
-		sigma_xy_input->GetInputValue(&sigma_xy,gauss);
-		if(dim==3){
-			sigma_xz_input->GetInputValue(&sigma_xz,gauss);
-			sigma_yz_input->GetInputValue(&sigma_yz,gauss);
-			sigma_zz_input->GetInputValue(&sigma_zz,gauss);
-		}
-
-		if(dim==2){
-			a = 0.;
-			b = 1.;
-			c = -sigma_yy -sigma_xx;
-			d = sigma_xx*sigma_yy - sigma_xy*sigma_xy;
-		}
-		else{
-			a = -1.;
-			b = sigma_xx+sigma_yy+sigma_zz;
-			c = -sigma_xx*sigma_yy -sigma_xx*sigma_zz -sigma_yy*sigma_zz + sigma_xy*sigma_xy +sigma_xz*sigma_xz +sigma_yz*sigma_yz;
-			d = sigma_xx*sigma_yy*sigma_zz - sigma_xx*sigma_yz*sigma_yz -sigma_yy*sigma_xz*sigma_xz - sigma_zz*sigma_xy*sigma_xy + 2.*sigma_xy*sigma_xz*sigma_yz;
-		}
-
-		/*Get roots of polynomials*/
-		cubic(a,b,c,d,x,&numroots);
-
-		/*Initialize maximum eigne value*/
-		if(numroots>0){
-			max = fabs(x[0]);
-		}
-		else{
-			_error_("No eigen value found");
-		}
-
-		/*Get max*/
-		for(int i=1;i<numroots;i++){
-			if(fabs(x[i])>max) max = fabs(x[i]);
-		}
-
-		maxprincipal[iv]=max;
-	}
-
-	/*Create input*/
-	this->AddInput(StressMaxPrincipalEnum,maxprincipal,P1Enum);
-
-	/*Clean up and return*/
-	xDelete<IssmDouble>(maxprincipal);
-	delete gauss;
-}
-/*}}}*/
-void       Element::ViscousHeatingCreateInput(void){/*{{{*/
-
-	/*Intermediaries*/
-	IssmDouble phi;
-	IssmDouble viscosity;
-	IssmDouble epsilon[6];
-	IssmDouble thickness;
-	IssmDouble *xyz_list = NULL;
-
-	/*Fetch number vertices and allocate memory*/
-	int         numvertices    = this->GetNumberOfVertices();
-	IssmDouble* viscousheating = xNew<IssmDouble>(numvertices);
-
-	/*Retrieve all inputs and parameters*/
-	this->GetVerticesCoordinatesBase(&xyz_list);
-	Input* vx_input        = this->GetInput(VxEnum); _assert_(vx_input);
-	Input* vy_input        = this->GetInput(VyEnum); _assert_(vy_input);
-	Input* vz_input        = this->GetInput(VzEnum); _assert_(vz_input);
-	Input* thickness_input = this->GetInput(ThicknessEnum); _assert_(thickness_input);
-
-	/*loop over vertices: */
-	Gauss* gauss=this->NewGauss();
-	for (int iv=0;iv<numvertices;iv++){
-		gauss->GaussVertex(iv);
-
-		thickness_input->GetInputValue(&thickness,gauss);
-
-		this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
-		this->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input);
-		this->GetPhi(&phi,&epsilon[0],viscosity);
-
-		viscousheating[iv]=phi*thickness;
-	}
-
-	/*Create PentaVertex input, which will hold the basal friction:*/
-	this->AddInput(ViscousHeatingEnum,viscousheating,P1Enum);
-
-	/*Clean up and return*/
-	xDelete<IssmDouble>(viscousheating);
-	delete gauss;
-}
-/*}}}*/
-void       Element::ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/
-	/*The effective strain rate is defined in Paterson 3d Ed p 91 eq 9,
-	 * and Cuffey p 303 eq 8.18:
-	 *
-	 *  2 eps_eff^2 = eps_xx^2 + eps_yy^2 + eps_zz^2 + 2(eps_xy^2 + eps_xz^2 + eps_yz^2)
-	 *
-	 *  or
-	 *
-	 *  eps_eff = 1/sqrt(2) sqrt( \sum_ij eps_ij^2 )
-	 *
-	 *          = 1/sqrt(2) ||eps||_F
-	 *
-	 *  where ||.||_F is the Frobenius norm */
-
-	/*Intermediaries*/
-	IssmDouble viscosity;
-	IssmDouble epsilon3d[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
-	IssmDouble epsilon2d[3]; /* epsilon=[exx,eyy,exy];            */
-	IssmDouble eps_eff;
-	IssmDouble eps0=1.e-27;
-
-	if(dim==3){
-		/* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */
-		this->StrainRateFS(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input,vz_input);
-		eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[3]*epsilon3d[3] +  epsilon3d[4]*epsilon3d[4] + epsilon3d[5]*epsilon3d[5] + epsilon3d[0]*epsilon3d[1]+eps0*eps0);
-	}
-	else{
-		/* eps_eff^2 = 1/2 ( exx^2 + eyy^2 + 2*exy^2 )*/
-		this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
-		eps_eff = 1./sqrt(2.)*sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + 2.*epsilon2d[2]*epsilon2d[2]);
-	}
-
-	/*Get viscosity*/
-	material->GetViscosity(&viscosity,eps_eff);
-
-	/*Assign output pointer*/
-	*pviscosity=viscosity;
-}
-/*}}}*/
-void       Element::ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surface_input){/*{{{*/
-	/*Compute the L1L2 viscosity
-	 *
-	 *      1
-	 * mu = - A^-1 (sigma'_e)^(1-n)
-	 *      2
-	 *
-	 * sigma'_e^2 = |sigma'_//|^2 + |sigma'_perp|^2 (see Perego 2012 eq. 17,18)
-	 *
-	 * L1L2 assumptions:
-	 *
-	 * (1) |eps_b|_// = A (|sigma'_//|^2 + |sigma'_perp|^2)^((n-1)/2) |sigma'_//|
-	 * (2) |sigma'_perp|^2 = |rho g (s-z) grad(s)|^2
-	 *
-	 * Assuming that n = 3, we have a polynom of degree 3 to solve (the only unkown is X=|sigma'_//|)
-	 *
-	 * A X^3 + A |rho g (s-z) grad(s)|^2 X - |eps_b|_// = 0     */
-
-	IssmDouble z,s,viscosity,p,q,delta;
-	IssmDouble tau_perp,tau_par,eps_b,A;
-	IssmDouble epsilon[5];   /*exx eyy exy exz eyz*/
-	IssmDouble slope[3];
-
-	/*Check that both inputs have been found*/
-	if (!vx_input || !vy_input || !surface_input) _error_("Input missing");
-
-	/*Get tau_perp*/
-	surface_input->GetInputValue(&s,gauss);
-	surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
-	z=this->GetZcoord(xyz_list,gauss);
-	tau_perp = matpar->GetRhoIce() * matpar->GetG() * fabs(s-z)*sqrt(slope[0]*slope[0]+slope[1]*slope[1]);
-
-	/* Get eps_b*/
-	this->StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
-	eps_b = sqrt(epsilon[0]*epsilon[0] + epsilon[1]*epsilon[1] + epsilon[0]*epsilon[1] + epsilon[2]*epsilon[2]);
-	if(eps_b==0.){
-		*pviscosity = 2.5e+17;
-		return;
-	}
-
-	/*Get A*/
-	_assert_(material->GetN()==3.0);
-	A=material->GetA();
-
-	/*Solve for tau_perp (http://fr.wikipedia.org/wiki/Méthode_de_Cardan)*/
-	p     = tau_perp *tau_perp;
-	q     = - eps_b/A;
-	delta = q *q + p*p*p*4./27.;
-	_assert_(delta>0);
-	tau_par = pow(0.5*(-q+sqrt(delta)),1./3.) - pow(0.5*(q+sqrt(delta)),1./3.);
-
-	/*Viscosity*/
-	viscosity = 1./(2.*A*(tau_par*tau_par + tau_perp*tau_perp));
-	_assert_(!xIsNan(viscosity));
-	_assert_(viscosity > 0.);
-
-	/*Assign output pointer*/
-	*pviscosity = viscosity;
-}/*}}}*/
-void       Element::ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
-
-	/*Intermediaries*/
-	IssmDouble viscosity;
-	IssmDouble epsilon3d[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/
-	IssmDouble epsilon2d[2];/* epsilon=[exx,exy];            */
-	IssmDouble eps_eff;
-
-	if(dim==3){
-		/* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */
-		this->StrainRateHO(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input);
-		eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[2]*epsilon3d[2] +  epsilon3d[3]*epsilon3d[3] + epsilon3d[4]*epsilon3d[4] + epsilon3d[0]*epsilon3d[1]);
-	}
-	else{
-		/* eps_eff^2 = 1/2 (exx^2 + 2*exy^2 )*/
-		this->StrainRateHO2dvertical(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
-		eps_eff = 1./sqrt(2.)*sqrt(epsilon2d[0]*epsilon2d[0] + 2.*epsilon2d[1]*epsilon2d[1]);
-	}
-
-	/*Get viscosity*/
-	material->GetViscosity(&viscosity,eps_eff);
-
-	/*Assign output pointer*/
-	*pviscosity=viscosity;
-}/*}}}*/
-void       Element::ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
-
-	/*Intermediaries*/
-	IssmDouble viscosity;
-	IssmDouble epsilon2d[3];/* epsilon=[exx,eyy,exy];    */
-	IssmDouble epsilon1d;   /* epsilon=[exx];    */
-	IssmDouble eps_eff;
-
-	 if(dim==2){
-		 /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy*/
-		 this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
-		 eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + epsilon2d[2]*epsilon2d[2] + epsilon2d[0]*epsilon2d[1]);
-	 }
-	 else{
-		 /* eps_eff^2 = 1/2 exx^2*/
-		 this->StrainRateSSA1d(&epsilon1d,xyz_list,gauss,vx_input);
-		 eps_eff = sqrt(epsilon1d*epsilon1d/2.);
-	 }
-
-	/*Get viscosity*/
-	material->GetViscosityBar(&viscosity,eps_eff);
-
-	/*Assign output pointer*/
-	*pviscosity=viscosity;
-}/*}}}*/
-void       Element::ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
-	this->material->GetViscosity2dDerivativeEpsSquare(pmu_prime,epsilon);
-}/*}}}*/
-void       Element::ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
-	this->material->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon);
-}/*}}}*/
-void       Element::ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
-	this->material->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon);
-}/*}}}*/
 void       Element::StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/
 	/*Compute the 3d Strain Rate (6 components):
 	 *
@@ -1929,6 +1637,99 @@
 	*epsilon = dvx[0];
 
 }/*}}}*/
+void       Element::StressMaxPrincipalCreateInput(void){/*{{{*/
+
+	/*Intermediaries*/
+	IssmDouble *xyz_list = NULL;
+	IssmDouble  sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz;
+	IssmDouble  a,b,c,d,x[3],max;
+	int         dim,numroots;
+
+	/*First: get stress tensor*/
+	this->ComputeStressTensor();
+
+	/*Get domain dimension*/
+	this->FindParam(&dim,DomainDimensionEnum);
+
+	/*Fetch number vertices and allocate memory*/
+	int         numvertices  = this->GetNumberOfVertices();
+	IssmDouble* maxprincipal = xNew<IssmDouble>(numvertices);
+
+	/*Retrieve all inputs and parameters*/
+	this->GetVerticesCoordinatesBase(&xyz_list);
+	Input* sigma_xx_input  = this->GetInput(StressTensorxxEnum); _assert_(sigma_xx_input);
+	Input* sigma_yy_input  = this->GetInput(StressTensoryyEnum); _assert_(sigma_yy_input);
+	Input* sigma_xy_input  = this->GetInput(StressTensorxyEnum); _assert_(sigma_xy_input);
+	Input* sigma_xz_input  = NULL;
+	Input* sigma_yz_input  = NULL;
+	Input* sigma_zz_input  = NULL;
+	if(dim==3){
+		sigma_xz_input  = this->GetInput(StressTensorxzEnum); _assert_(sigma_xz_input);
+		sigma_yz_input  = this->GetInput(StressTensoryzEnum); _assert_(sigma_yz_input);
+		sigma_zz_input  = this->GetInput(StressTensorzzEnum); _assert_(sigma_zz_input);
+	}
+
+	/*loop over vertices: */
+	Gauss* gauss=this->NewGauss();
+	for (int iv=0;iv<numvertices;iv++){
+		gauss->GaussVertex(iv);
+
+		sigma_xx_input->GetInputValue(&sigma_xx,gauss);
+		sigma_yy_input->GetInputValue(&sigma_yy,gauss);
+		sigma_xy_input->GetInputValue(&sigma_xy,gauss);
+		if(dim==3){
+			sigma_xz_input->GetInputValue(&sigma_xz,gauss);
+			sigma_yz_input->GetInputValue(&sigma_yz,gauss);
+			sigma_zz_input->GetInputValue(&sigma_zz,gauss);
+		}
+
+		if(dim==2){
+			a = 0.;
+			b = 1.;
+			c = -sigma_yy -sigma_xx;
+			d = sigma_xx*sigma_yy - sigma_xy*sigma_xy;
+		}
+		else{
+			a = -1.;
+			b = sigma_xx+sigma_yy+sigma_zz;
+			c = -sigma_xx*sigma_yy -sigma_xx*sigma_zz -sigma_yy*sigma_zz + sigma_xy*sigma_xy +sigma_xz*sigma_xz +sigma_yz*sigma_yz;
+			d = sigma_xx*sigma_yy*sigma_zz - sigma_xx*sigma_yz*sigma_yz -sigma_yy*sigma_xz*sigma_xz - sigma_zz*sigma_xy*sigma_xy + 2.*sigma_xy*sigma_xz*sigma_yz;
+		}
+
+		/*Get roots of polynomials*/
+		cubic(a,b,c,d,x,&numroots);
+
+		/*Initialize maximum eigne value*/
+		if(numroots>0){
+			max = fabs(x[0]);
+		}
+		else{
+			_error_("No eigen value found");
+		}
+
+		/*Get max*/
+		for(int i=1;i<numroots;i++){
+			if(fabs(x[i])>max) max = fabs(x[i]);
+		}
+
+		maxprincipal[iv]=max;
+	}
+
+	/*Create input*/
+	this->AddInput(StressMaxPrincipalEnum,maxprincipal,P1Enum);
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(maxprincipal);
+	delete gauss;
+}
+/*}}}*/
+void       Element::ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){/*{{{*/
+	matpar->ThermalToEnthalpy(penthalpy,temperature,waterfraction,pressure);
+}/*}}}*/
+IssmDouble Element::TMeltingPoint(IssmDouble pressure){/*{{{*/
+	_assert_(matpar);
+	return this->matpar->TMeltingPoint(pressure);
+}/*}}}*/
 void       Element::TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,int transformenum){/*{{{*/
 
 	/*All nodes have the same Coordinate System*/
@@ -2152,3 +1953,202 @@
 	xDelete<IssmDouble>(transform);
 	xDelete<IssmDouble>(values);
 }/*}}}*/
+void       Element::ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/
+	/*The effective strain rate is defined in Paterson 3d Ed p 91 eq 9,
+	 * and Cuffey p 303 eq 8.18:
+	 *
+	 *  2 eps_eff^2 = eps_xx^2 + eps_yy^2 + eps_zz^2 + 2(eps_xy^2 + eps_xz^2 + eps_yz^2)
+	 *
+	 *  or
+	 *
+	 *  eps_eff = 1/sqrt(2) sqrt( \sum_ij eps_ij^2 )
+	 *
+	 *          = 1/sqrt(2) ||eps||_F
+	 *
+	 *  where ||.||_F is the Frobenius norm */
+
+	/*Intermediaries*/
+	IssmDouble viscosity;
+	IssmDouble epsilon3d[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
+	IssmDouble epsilon2d[3]; /* epsilon=[exx,eyy,exy];            */
+	IssmDouble eps_eff;
+	IssmDouble eps0=1.e-27;
+
+	if(dim==3){
+		/* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */
+		this->StrainRateFS(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input,vz_input);
+		eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[3]*epsilon3d[3] +  epsilon3d[4]*epsilon3d[4] + epsilon3d[5]*epsilon3d[5] + epsilon3d[0]*epsilon3d[1]+eps0*eps0);
+	}
+	else{
+		/* eps_eff^2 = 1/2 ( exx^2 + eyy^2 + 2*exy^2 )*/
+		this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
+		eps_eff = 1./sqrt(2.)*sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + 2.*epsilon2d[2]*epsilon2d[2]);
+	}
+
+	/*Get viscosity*/
+	material->GetViscosity(&viscosity,eps_eff);
+
+	/*Assign output pointer*/
+	*pviscosity=viscosity;
+}
+/*}}}*/
+void       Element::ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
+	this->material->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon);
+}/*}}}*/
+void       Element::ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
+
+	/*Intermediaries*/
+	IssmDouble viscosity;
+	IssmDouble epsilon3d[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/
+	IssmDouble epsilon2d[2];/* epsilon=[exx,exy];            */
+	IssmDouble eps_eff;
+
+	if(dim==3){
+		/* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */
+		this->StrainRateHO(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input);
+		eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[2]*epsilon3d[2] +  epsilon3d[3]*epsilon3d[3] + epsilon3d[4]*epsilon3d[4] + epsilon3d[0]*epsilon3d[1]);
+	}
+	else{
+		/* eps_eff^2 = 1/2 (exx^2 + 2*exy^2 )*/
+		this->StrainRateHO2dvertical(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
+		eps_eff = 1./sqrt(2.)*sqrt(epsilon2d[0]*epsilon2d[0] + 2.*epsilon2d[1]*epsilon2d[1]);
+	}
+
+	/*Get viscosity*/
+	material->GetViscosity(&viscosity,eps_eff);
+
+	/*Assign output pointer*/
+	*pviscosity=viscosity;
+}/*}}}*/
+void       Element::ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
+	this->material->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon);
+}/*}}}*/
+void       Element::ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surface_input){/*{{{*/
+	/*Compute the L1L2 viscosity
+	 *
+	 *      1
+	 * mu = - A^-1 (sigma'_e)^(1-n)
+	 *      2
+	 *
+	 * sigma'_e^2 = |sigma'_//|^2 + |sigma'_perp|^2 (see Perego 2012 eq. 17,18)
+	 *
+	 * L1L2 assumptions:
+	 *
+	 * (1) |eps_b|_// = A (|sigma'_//|^2 + |sigma'_perp|^2)^((n-1)/2) |sigma'_//|
+	 * (2) |sigma'_perp|^2 = |rho g (s-z) grad(s)|^2
+	 *
+	 * Assuming that n = 3, we have a polynom of degree 3 to solve (the only unkown is X=|sigma'_//|)
+	 *
+	 * A X^3 + A |rho g (s-z) grad(s)|^2 X - |eps_b|_// = 0     */
+
+	IssmDouble z,s,viscosity,p,q,delta;
+	IssmDouble tau_perp,tau_par,eps_b,A;
+	IssmDouble epsilon[5];   /*exx eyy exy exz eyz*/
+	IssmDouble slope[3];
+
+	/*Check that both inputs have been found*/
+	if (!vx_input || !vy_input || !surface_input) _error_("Input missing");
+
+	/*Get tau_perp*/
+	surface_input->GetInputValue(&s,gauss);
+	surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
+	z=this->GetZcoord(xyz_list,gauss);
+	tau_perp = matpar->GetRhoIce() * matpar->GetG() * fabs(s-z)*sqrt(slope[0]*slope[0]+slope[1]*slope[1]);
+
+	/* Get eps_b*/
+	this->StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
+	eps_b = sqrt(epsilon[0]*epsilon[0] + epsilon[1]*epsilon[1] + epsilon[0]*epsilon[1] + epsilon[2]*epsilon[2]);
+	if(eps_b==0.){
+		*pviscosity = 2.5e+17;
+		return;
+	}
+
+	/*Get A*/
+	_assert_(material->GetN()==3.0);
+	A=material->GetA();
+
+	/*Solve for tau_perp (http://fr.wikipedia.org/wiki/Méthode_de_Cardan)*/
+	p     = tau_perp *tau_perp;
+	q     = - eps_b/A;
+	delta = q *q + p*p*p*4./27.;
+	_assert_(delta>0);
+	tau_par = pow(0.5*(-q+sqrt(delta)),1./3.) - pow(0.5*(q+sqrt(delta)),1./3.);
+
+	/*Viscosity*/
+	viscosity = 1./(2.*A*(tau_par*tau_par + tau_perp*tau_perp));
+	_assert_(!xIsNan(viscosity));
+	_assert_(viscosity > 0.);
+
+	/*Assign output pointer*/
+	*pviscosity = viscosity;
+}/*}}}*/
+void       Element::ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
+
+	/*Intermediaries*/
+	IssmDouble viscosity;
+	IssmDouble epsilon2d[3];/* epsilon=[exx,eyy,exy];    */
+	IssmDouble epsilon1d;   /* epsilon=[exx];    */
+	IssmDouble eps_eff;
+
+	 if(dim==2){
+		 /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy*/
+		 this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
+		 eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + epsilon2d[2]*epsilon2d[2] + epsilon2d[0]*epsilon2d[1]);
+	 }
+	 else{
+		 /* eps_eff^2 = 1/2 exx^2*/
+		 this->StrainRateSSA1d(&epsilon1d,xyz_list,gauss,vx_input);
+		 eps_eff = sqrt(epsilon1d*epsilon1d/2.);
+	 }
+
+	/*Get viscosity*/
+	material->GetViscosityBar(&viscosity,eps_eff);
+
+	/*Assign output pointer*/
+	*pviscosity=viscosity;
+}/*}}}*/
+void       Element::ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
+	this->material->GetViscosity2dDerivativeEpsSquare(pmu_prime,epsilon);
+}/*}}}*/
+void       Element::ViscousHeatingCreateInput(void){/*{{{*/
+
+	/*Intermediaries*/
+	IssmDouble phi;
+	IssmDouble viscosity;
+	IssmDouble epsilon[6];
+	IssmDouble thickness;
+	IssmDouble *xyz_list = NULL;
+
+	/*Fetch number vertices and allocate memory*/
+	int         numvertices    = this->GetNumberOfVertices();
+	IssmDouble* viscousheating = xNew<IssmDouble>(numvertices);
+
+	/*Retrieve all inputs and parameters*/
+	this->GetVerticesCoordinatesBase(&xyz_list);
+	Input* vx_input        = this->GetInput(VxEnum); _assert_(vx_input);
+	Input* vy_input        = this->GetInput(VyEnum); _assert_(vy_input);
+	Input* vz_input        = this->GetInput(VzEnum); _assert_(vz_input);
+	Input* thickness_input = this->GetInput(ThicknessEnum); _assert_(thickness_input);
+
+	/*loop over vertices: */
+	Gauss* gauss=this->NewGauss();
+	for (int iv=0;iv<numvertices;iv++){
+		gauss->GaussVertex(iv);
+
+		thickness_input->GetInputValue(&thickness,gauss);
+
+		this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
+		this->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input);
+		this->GetPhi(&phi,&epsilon[0],viscosity);
+
+		viscousheating[iv]=phi*thickness;
+	}
+
+	/*Create PentaVertex input, which will hold the basal friction:*/
+	this->AddInput(ViscousHeatingEnum,viscousheating,P1Enum);
+
+	/*Clean up and return*/
+	xDelete<IssmDouble>(viscousheating);
+	delete gauss;
+}
+/*}}}*/
