Index: /issm/trunk/src/c/Container/Inputs.cpp
===================================================================
--- /issm/trunk/src/c/Container/Inputs.cpp	(revision 4994)
+++ /issm/trunk/src/c/Container/Inputs.cpp	(revision 4995)
@@ -301,4 +301,7 @@
 	Input* input=NULL;
 
+	/*In debugging mode, check that the input is not a NULL pointer*/
+	ISSMASSERT(in_input);
+
 	for ( object=objects.begin() ; object < objects.end(); object++ ){
 
Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 4994)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 4995)
@@ -712,6 +712,10 @@
 	else{
 
-		/*Depth Average B*/
-		this->InputDepthAverageAtBase(RheologyBEnum,RheologyB2dEnum,MaterialsEnum);
+		/*Depth Average B and put it in inputs*/
+		Penta* penta=GetBasalElement();
+		penta->InputDepthAverageAtBase(RheologyBEnum,RheologyB2dEnum,MaterialsEnum);
+		Input* B_input=penta->matice->inputs->GetInput(RheologyB2dEnum);
+		Input* B_copy=(Input*)B_input->copy();
+		this->matice->inputs->AddInput((Input*)B_copy);
 
 		tria=(Tria*)SpawnTria(3,4,5); //grids 3, 4 and 5 make the new tria (upper face).
@@ -721,4 +725,5 @@
 		/*delete B average*/
 		this->matice->inputs->DeleteInput(RheologyB2dEnum);
+		penta->matice->inputs->DeleteInput(RheologyB2dEnum);
 
 		return J;
@@ -4661,4 +4666,35 @@
 }
 /*}}}*/
+/*FUNCTION Penta::GetBasalElement{{{1*/
+Penta* Penta::GetBasalElement(void){
+
+	/*Output*/
+	Penta* penta=NULL;
+
+	/*Go through all elements till the bed is reached*/
+	penta=this;
+	for(;;){
+
+		/*Stop if we have reached the surface, else, take lower penta*/
+		if (penta->IsOnBed()) break;
+
+		/* get lower Penta*/
+		penta=penta->GetLowerElement();
+		ISSMASSERT(penta->Id()!=this->id);
+	}
+
+	/*return output*/
+	return penta;
+}
+/*}}}*/
+/*FUNCTION Penta::GetLowerElement{{{1*/
+Penta* Penta::GetLowerElement(void){
+
+	Penta* upper_penta=NULL;
+	upper_penta=(Penta*)neighbors[0]; //first one (0) under, second one (1) above
+	return upper_penta;
+
+}
+/*}}}*/
 /*FUNCTION Penta::GetUpperElement{{{1*/
 Penta* Penta::GetUpperElement(void){
@@ -4708,4 +4744,7 @@
 	}
 	else{
+		/*Gradient is computed on bed only (B2d)*/
+		if (!onbed) return;
+
 		/*Depth Average B*/
 		this->InputDepthAverageAtBase(RheologyBEnum,RheologyB2dEnum,MaterialsEnum);
@@ -5471,4 +5510,13 @@
 }
 /*}}}*/
+/*FUNCTION Penta::IsOnBed{{{1*/
+bool Penta::IsOnBed(void){
+
+	bool onbed;
+	inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+	return onbed;
+
+}
+/*}}}*/
 /*FUNCTION Penta::ReduceMatrixStokes {{{1*/
 void Penta::ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp){
Index: /issm/trunk/src/c/objects/Elements/Penta.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.h	(revision 4994)
+++ /issm/trunk/src/c/objects/Elements/Penta.h	(revision 4995)
@@ -156,4 +156,6 @@
 		void    GetStrainRate3d(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input, Input* vz_input);
 		Penta*  GetUpperElement(void);
+		Penta*  GetLowerElement(void);
+		Penta*  GetBasalElement(void);
 		void	  InputExtrude(int enum_type,int object_type);
 		void    InputUpdateFromSolutionAdjointHoriz( double* solutiong);
@@ -170,4 +172,5 @@
 		bool	  IsInput(int name);
 		bool	  IsOnSurface(void);
+		bool	  IsOnBed(void);
 		void	  ReduceMatrixStokes(double* Ke_reduced, double* Ke_temp);
 		void	  ReduceVectorStokes(double* Pe_reduced, double* Ke_temp, double* Pe_temp);