Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 16744)
+++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp	(revision 16745)
@@ -202,4 +202,6 @@
 	IssmDouble* values        = xNew<IssmDouble>(numnodes);
 	IssmDouble* pressure      = xNew<IssmDouble>(numnodes);
+	IssmDouble* surface       = xNew<IssmDouble>(numnodes);
+	IssmDouble* B             = xNew<IssmDouble>(numnodes);
 	IssmDouble* temperature   = xNew<IssmDouble>(numnodes);
 	IssmDouble* waterfraction = xNew<IssmDouble>(numnodes);
@@ -215,5 +217,5 @@
 	/*Get all inputs and parameters*/
 	element->GetInputValue(&converged,ConvergedEnum);
-	element->GetInputListOnNodes(pressure,PressureEnum);
+	element->GetInputListOnNodes(&pressure[0],PressureEnum);
 	if(converged){
 		for(i=0;i<numnodes;i++){
@@ -229,4 +231,5 @@
 		 * otherwise the rheology could be negative*/
 		element->FindParam(&rheology_law,MaterialsRheologyLawEnum);
+		element->GetInputListOnNodes(&surface[0],SurfaceEnum);
 		switch(rheology_law){
 			case NoneEnum:
@@ -234,32 +237,18 @@
 				break;
 			case PatersonEnum:
-				for(i=0;i<numnodes;i++) T_average+=values[i]/reCast<IssmDouble>(numnodes);
-				B_average=Paterson(T_average);
-				element->AddMaterialInput(MaterialsRheologyBEnum,&B_average,P0Enum);
+				for(i=0;i<numnodes;i++) B[i]=Paterson(values[i]);
+				element->AddMaterialInput(MaterialsRheologyBEnum,&B[0],P1Enum);
 				break;
-			case ArrheniusEnum:{
-				Input* surface_input=element->GetInput(SurfaceEnum); _assert_(surface_input);
-				surface_input->GetInputAverage(&s_average);
+			case ArrheniusEnum:
 				element->GetVerticesCoordinates(&xyz_list);
-				for(i=0;i<numnodes;i++) T_average+=values[i]/reCast<IssmDouble>(numnodes);
-				//B_average=Arrhenius(T_average,
-							//s_average-((xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2]+xyz_list[3][2]+xyz_list[4][2]+xyz_list[5][2])/6.0),
-							//element->GetMaticeParameter(MaterialsRheologyNEnum));
-				element->AddMaterialInput(MaterialsRheologyBEnum,&B_average,P0Enum);
+				for(i=0;i<numnodes;i++) B[i]=Arrhenius(values[i],surface[i]-xyz_list[i*3+2],element->GetMaterialParameter(MaterialsRheologyNEnum));
+				element->AddMaterialInput(MaterialsRheologyBEnum,&B[0],P1Enum);
 				break;
-				}
 			case LliboutryDuvalEnum:
-				for(i=0;i<numnodes;i++) T_average+=values[i]/reCast<IssmDouble>(numnodes);
-				for(i=0;i<numnodes;i++) P_average+=pressure[i]/reCast<IssmDouble>(numnodes);
-				B_average=LliboutryDuval(T_average,P_average,
-							element->GetMaterialParameter(MaterialsRheologyNEnum),
-							element->GetMaterialParameter(MaterialsBetaEnum),
-							element->GetMaterialParameter(ConstantsReferencetemperatureEnum),
-							element->GetMaterialParameter(MaterialsHeatcapacityEnum),
-							element->GetMaterialParameter(MaterialsLatentheatEnum));
-				element->AddMaterialInput(MaterialsRheologyBEnum,&B_average,P0Enum);
-				break;
-			default:
-				_error_("Rheology law " << EnumToStringx(rheology_law) << " not supported yet");
+				//for(i=0;i<numnodes;i++) B[i]=LliboutryDuval(values[i],pressure[i],material->GetN(),matpar->GetBeta(),matpar->GetReferenceTemperature(),matpar->GetHeatCapacity(),matpar->GetLatentHeat()); for(i=0;i<numnodes;i++) B[i]=Paterson(values[i]); 
+				for(i=0;i<numnodes;i++) B[i]=Paterson(values[i]);
+				element->AddMaterialInput(MaterialsRheologyBEnum,&B[0],P1Enum); 
+				break; 
+			default: _error_("Rheology law " << EnumToStringx(rheology_law) << " not supported yet");
 		}
 	}
@@ -271,4 +260,6 @@
 	xDelete<IssmDouble>(values);
 	xDelete<IssmDouble>(pressure);
+	xDelete<IssmDouble>(surface);
+	xDelete<IssmDouble>(B);
 	xDelete<IssmDouble>(temperature);
 	xDelete<IssmDouble>(waterfraction);
Index: /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp	(revision 16744)
+++ /issm/trunk-jpl/src/c/analyses/ThermalAnalysis.cpp	(revision 16745)
@@ -119,5 +119,4 @@
 	bool        converged;
 	int         i,rheology_law;
-	IssmDouble  B_average,s_average,T_average=0.;
 	int        *doflist   = NULL;
 	IssmDouble *xyz_list  = NULL;
@@ -130,4 +129,6 @@
 	element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
 	IssmDouble* values    = xNew<IssmDouble>(numnodes);
+	IssmDouble* surface   = xNew<IssmDouble>(numnodes);
+	IssmDouble* B         = xNew<IssmDouble>(numnodes);
 
 	/*Use the dof list to index into the solution vector: */
@@ -144,5 +145,5 @@
 	if(hack){
 		IssmDouble* pressure = xNew<IssmDouble>(numnodes);
-		element->GetInputListOnVertices(pressure,PressureEnum);
+		element->GetInputListOnNodes(&pressure[0],PressureEnum);
 		for(i=0;i<numnodes;i++){
 			if(values[i]>element->TMeltingPoint(pressure[i]))     values[i]=element->TMeltingPoint(pressure[i]);
@@ -160,4 +161,5 @@
 		 * otherwise the rheology could be negative*/
 		element->FindParam(&rheology_law,MaterialsRheologyLawEnum);
+		element->GetInputListOnNodes(&surface[0],SurfaceEnum);
 		switch(rheology_law){
 			case NoneEnum:
@@ -165,22 +167,15 @@
 				break;
 			case PatersonEnum:
-				for(i=0;i<numnodes;i++) T_average+=values[i]/reCast<IssmDouble>(numnodes);
-				B_average=Paterson(T_average);
-				element->AddMaterialInput(MaterialsRheologyBEnum,&B_average,P0Enum);
+				for(i=0;i<numnodes;i++) B[i]=Paterson(values[i]);
+				element->AddMaterialInput(MaterialsRheologyBEnum,&B[0],P1Enum);
 				break;
 			case ArrheniusEnum:{
-				Input* surface_input=element->GetInput(SurfaceEnum); _assert_(surface_input);
-				surface_input->GetInputAverage(&s_average);
 				element->GetVerticesCoordinates(&xyz_list);
-				for(i=0;i<numnodes;i++) T_average+=values[i]/reCast<IssmDouble>(numnodes);
-				//B_average=Arrhenius(T_average,
-							//s_average-((xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2]+xyz_list[3][2]+xyz_list[4][2]+xyz_list[5][2])/6.0),
-							//element->GetMaticeParameter(MaterialsRheologyNEnum));
-				element->AddMaterialInput(MaterialsRheologyBEnum,&B_average,P0Enum);
+				for(i=0;i<numnodes;i++) B[i]=Arrhenius(values[i],surface[i]-xyz_list[i*3+2],element->GetMaterialParameter(MaterialsRheologyNEnum));
+				element->AddMaterialInput(MaterialsRheologyBEnum,&B[0],P1Enum);
 				break;
 				}
 			default:
 				_error_("Rheology law " << EnumToStringx(rheology_law) << " not supported yet");
-
 		}
 	}
@@ -191,4 +186,6 @@
 	/*Free ressources:*/
 	xDelete<IssmDouble>(values);
+	xDelete<IssmDouble>(surface);
+	xDelete<IssmDouble>(B);
 	xDelete<IssmDouble>(xyz_list);
 	xDelete<int>(doflist);
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16744)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16745)
@@ -4929,5 +4929,4 @@
 	IssmDouble waterfraction[numdof];
 	IssmDouble B[numdof];
-	IssmDouble B_average,s_average;
 	int*       doflist=NULL;
 
@@ -4964,5 +4963,5 @@
 		 * otherwise the rheology could be negative*/
 		this->parameters->FindParam(&rheology_law,MaterialsRheologyLawEnum);
-		GetInputListOnVertices(&surface[0],PressureEnum);
+		GetInputListOnVertices(&surface[0],SurfaceEnum);
 		switch(rheology_law){
 			case NoneEnum: