Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16742)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16743)
@@ -4851,5 +4851,4 @@
 	IssmDouble  values[numdof];
 	IssmDouble  B[numdof],surface[numdof];
-	IssmDouble  B_average;
 	int        *doflist = NULL;
 	bool        hack    = false;
@@ -4890,4 +4889,5 @@
 		 * otherwise the rheology could be negative*/
 		this->parameters->FindParam(&rheology_law,MaterialsRheologyLawEnum);
+		GetInputListOnVertices(&surface[0],SurfaceEnum);
 		switch(rheology_law){
 			case NoneEnum:
@@ -4895,10 +4895,8 @@
 				break;
 			case PatersonEnum:
-				B_average=Paterson((values[0]+values[1]+values[2]+values[3]+values[4]+values[5])/6.0);
-				for(i=0;i<numdof;i++) B[i]=B_average;
+				for(i=0;i<numdof;i++) B[i]=Paterson(values[i]);
 				this->material->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum));
 				break;
 			case ArrheniusEnum:
-				GetInputListOnVertices(&surface[0],SurfaceEnum);
 				for(i=0;i<numdof;i++) B[i]=Arrhenius(values[i],surface[i]-xyz_list[i][2],material->GetN());
 				this->material->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum));
@@ -4926,5 +4924,6 @@
 	IssmDouble xyz_list[NUMVERTICES][3];
 	IssmDouble values[numdof];
-	IssmDouble pressure[NUMVERTICES];
+	IssmDouble pressure[numdof];
+	IssmDouble surface[numdof];
 	IssmDouble temperatures[numdof];
 	IssmDouble waterfraction[numdof];
@@ -4965,4 +4964,5 @@
 		 * otherwise the rheology could be negative*/
 		this->parameters->FindParam(&rheology_law,MaterialsRheologyLawEnum);
+		GetInputListOnVertices(&surface[0],PressureEnum);
 		switch(rheology_law){
 			case NoneEnum:
@@ -4970,22 +4970,13 @@
 				break;
 			case PatersonEnum:
-				B_average=Paterson((temperatures[0]+temperatures[1]+temperatures[2]+temperatures[3]+temperatures[4]+temperatures[5])/6.0);
-				for(i=0;i<numdof;i++) B[i]=B_average;
+				for(i=0;i<numdof;i++) B[i]=Paterson(values[i]);
 				this->material->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum));
 				break;
 			case ArrheniusEnum:
-				surface_input->GetInputAverage(&s_average);
-				B_average=Arrhenius((temperatures[0]+temperatures[1]+temperatures[2]+temperatures[3]+temperatures[4]+temperatures[5])/6.0,
-							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),
-							material->GetN());
-				for(i=0;i<numdof;i++) B[i]=B_average;
+				for(i=0;i<numdof;i++) B[i]=Arrhenius(values[i],surface[i]-xyz_list[i][2],material->GetN());
 				this->material->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum));
 				break;
 			case LliboutryDuvalEnum:
-				B_average=LliboutryDuval((values[0]+values[1]+values[2]+values[3]+values[4]+values[5])/6.0,
-																 (pressure[0]+pressure[1]+pressure[2]+pressure[3]+pressure[4]+pressure[5])/6.0,
-																 material->GetN(), 
-																 matpar->GetBeta(), matpar->GetReferenceTemperature(), matpar->GetHeatCapacity(), matpar->GetLatentHeat());
-				for(i=0;i<numdof;i++) B[i]=B_average;
+				for(i=0;i<numdof;i++) B[i]=LliboutryDuval(values[i],pressure[i],material->GetN(),matpar->GetBeta(),matpar->GetReferenceTemperature(),matpar->GetHeatCapacity(),matpar->GetLatentHeat());
 				this->material->inputs->AddInput(new PentaInput(MaterialsRheologyBEnum,B,P1Enum));
 				break;