Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18920)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18921)
@@ -125,10 +125,4 @@
 
 /*Other*/
-void       Penta::AddInput(int input_enum,IssmDouble* values, int interpolation_enum){/*{{{*/
-
-	_assert_(this->inputs);
-	this->inputs->AddInput(new PentaInput(input_enum,values,interpolation_enum));
-}
-/*}}}*/
 void       Penta::AddBasalInput(int input_enum,IssmDouble* values, int interpolation_enum){/*{{{*/
 
@@ -154,4 +148,69 @@
 		else _error_("not implemented yet");
 	}
+}
+/*}}}*/
+void       Penta::AddInput(int input_enum,IssmDouble* values, int interpolation_enum){/*{{{*/
+
+	_assert_(this->inputs);
+	this->inputs->AddInput(new PentaInput(input_enum,values,interpolation_enum));
+}
+/*}}}*/
+void       Penta::AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part){/*{{{*/
+	_error_("Not supported yet!");
+}
+/*}}}*/
+void       Penta::CalvingRateLevermann(){/*{{{*/
+
+	IssmDouble  xyz_list[NUMVERTICES][3];
+	GaussPenta* gauss=NULL;
+	IssmDouble  vx,vy,vel;
+	IssmDouble  strainparallel;
+	IssmDouble  propcoeff;
+	IssmDouble  strainperpendicular;
+	IssmDouble  calvingratex[NUMVERTICES];
+	IssmDouble  calvingratey[NUMVERTICES];
+	IssmDouble  calvingrate[NUMVERTICES];
+
+
+	/* Get node coordinates and dof list: */
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+
+	/*Retrieve all inputs and parameters we will need*/
+	Input* vx_input=inputs->GetInput(VxEnum);																		_assert_(vx_input);
+	Input* vy_input=inputs->GetInput(VyEnum);																		_assert_(vy_input);
+	Input* strainparallel_input=inputs->GetInput(StrainRateparallelEnum);								_assert_(strainparallel_input);
+	Input* strainperpendicular_input=inputs->GetInput(StrainRateperpendicularEnum);              _assert_(strainperpendicular_input);
+	Input* levermanncoeff_input=inputs->GetInput(CalvinglevermannCoeffEnum);                     _assert_(levermanncoeff_input);
+
+	/* Start looping on the number of vertices: */
+	gauss=new GaussPenta();
+	for (int iv=0;iv<NUMVERTICES;iv++){
+		gauss->GaussVertex(iv);
+
+		/* Get the value we need*/
+		vx_input->GetInputValue(&vx,gauss);
+		vy_input->GetInputValue(&vy,gauss);
+		vel=vx*vx+vy*vy;
+		strainparallel_input->GetInputValue(&strainparallel,gauss);
+		strainperpendicular_input->GetInputValue(&strainperpendicular,gauss);
+		levermanncoeff_input->GetInputValue(&propcoeff,gauss);
+
+		/*Calving rate proportionnal to the positive product of the strain rate along the ice flow direction and the strain rate perpendicular to the ice flow */
+		calvingrate[iv]=propcoeff*strainparallel*strainperpendicular;	
+		if(calvingrate[iv]<0){
+			calvingrate[iv]=0;
+		}
+		calvingratex[iv]=calvingrate[iv]*vx/(vel+1.e-6);
+		calvingratey[iv]=calvingrate[iv]*vy/(vel+1.e-6);
+	}
+
+	/*Add input*/
+	this->inputs->AddInput(new PentaInput(CalvingratexEnum,&calvingratex[0],P1Enum));
+	this->inputs->AddInput(new PentaInput(CalvingrateyEnum,&calvingratey[0],P1Enum));
+	this->inputs->AddInput(new PentaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum));
+
+	/*Clean up and return*/
+	delete gauss;
+
 }
 /*}}}*/
@@ -243,4 +302,55 @@
 }
 /*}}}*/
+void       Penta::ComputeDeviatoricStressTensor(){/*{{{*/
+
+	IssmDouble  xyz_list[NUMVERTICES][3];
+	IssmDouble  viscosity;
+	IssmDouble  epsilon[6]; /* epsilon=[exx,eyy,exy];*/
+	IssmDouble  tau_xx[NUMVERTICES];
+	IssmDouble	tau_yy[NUMVERTICES];
+	IssmDouble	tau_zz[NUMVERTICES];
+	IssmDouble  tau_xy[NUMVERTICES];
+	IssmDouble	tau_xz[NUMVERTICES];
+	IssmDouble	tau_yz[NUMVERTICES];
+	GaussPenta* gauss=NULL;
+
+	/* Get node coordinates and dof list: */
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+
+	/*Retrieve all inputs we will be needing: */
+	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);
+
+	/* Start looping on the number of vertices: */
+	gauss=new GaussPenta();
+	for (int iv=0;iv<NUMVERTICES;iv++){
+		gauss->GaussVertex(iv);
+
+		/*Compute strain rate viscosity and pressure: */
+		this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
+		this->ViscosityFS(&viscosity,3,&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
+
+		/*Compute Stress*/
+		tau_xx[iv]=2*viscosity*epsilon[0]; // tau = nu eps 
+		tau_yy[iv]=2*viscosity*epsilon[1];
+		tau_zz[iv]=2*viscosity*epsilon[2];
+		tau_xy[iv]=2*viscosity*epsilon[3];
+		tau_xz[iv]=2*viscosity*epsilon[4];
+		tau_yz[iv]=2*viscosity*epsilon[5];
+	}
+
+	/*Add Stress tensor components into inputs*/
+	this->inputs->AddInput(new PentaInput(DeviatoricStressxxEnum,&tau_xx[0],P1Enum));
+	this->inputs->AddInput(new PentaInput(DeviatoricStressxyEnum,&tau_xy[0],P1Enum));
+	this->inputs->AddInput(new PentaInput(DeviatoricStressxzEnum,&tau_xz[0],P1Enum));
+	this->inputs->AddInput(new PentaInput(DeviatoricStressyyEnum,&tau_yy[0],P1Enum));
+	this->inputs->AddInput(new PentaInput(DeviatoricStressyzEnum,&tau_yz[0],P1Enum));
+	this->inputs->AddInput(new PentaInput(DeviatoricStresszzEnum,&tau_zz[0],P1Enum));
+
+	/*Clean up and return*/
+	delete gauss;
+}
+/*}}}*/
 void       Penta::ComputeStressTensor(){/*{{{*/
 
@@ -296,281 +406,4 @@
 }
 /*}}}*/
-void       Penta::ComputeDeviatoricStressTensor(){/*{{{*/
-
-	IssmDouble  xyz_list[NUMVERTICES][3];
-	IssmDouble  viscosity;
-	IssmDouble  epsilon[6]; /* epsilon=[exx,eyy,exy];*/
-	IssmDouble  tau_xx[NUMVERTICES];
-	IssmDouble	tau_yy[NUMVERTICES];
-	IssmDouble	tau_zz[NUMVERTICES];
-	IssmDouble  tau_xy[NUMVERTICES];
-	IssmDouble	tau_xz[NUMVERTICES];
-	IssmDouble	tau_yz[NUMVERTICES];
-	GaussPenta* gauss=NULL;
-
-	/* Get node coordinates and dof list: */
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-
-	/*Retrieve all inputs we will be needing: */
-	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);
-
-	/* Start looping on the number of vertices: */
-	gauss=new GaussPenta();
-	for (int iv=0;iv<NUMVERTICES;iv++){
-		gauss->GaussVertex(iv);
-
-		/*Compute strain rate viscosity and pressure: */
-		this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
-		this->ViscosityFS(&viscosity,3,&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
-
-		/*Compute Stress*/
-		tau_xx[iv]=2*viscosity*epsilon[0]; // tau = nu eps 
-		tau_yy[iv]=2*viscosity*epsilon[1];
-		tau_zz[iv]=2*viscosity*epsilon[2];
-		tau_xy[iv]=2*viscosity*epsilon[3];
-		tau_xz[iv]=2*viscosity*epsilon[4];
-		tau_yz[iv]=2*viscosity*epsilon[5];
-	}
-
-	/*Add Stress tensor components into inputs*/
-	this->inputs->AddInput(new PentaInput(DeviatoricStressxxEnum,&tau_xx[0],P1Enum));
-	this->inputs->AddInput(new PentaInput(DeviatoricStressxyEnum,&tau_xy[0],P1Enum));
-	this->inputs->AddInput(new PentaInput(DeviatoricStressxzEnum,&tau_xz[0],P1Enum));
-	this->inputs->AddInput(new PentaInput(DeviatoricStressyyEnum,&tau_yy[0],P1Enum));
-	this->inputs->AddInput(new PentaInput(DeviatoricStressyzEnum,&tau_yz[0],P1Enum));
-	this->inputs->AddInput(new PentaInput(DeviatoricStresszzEnum,&tau_zz[0],P1Enum));
-
-	/*Clean up and return*/
-	delete gauss;
-}
-/*}}}*/
-void       Penta::CalvingRateLevermann(){/*{{{*/
-
-	IssmDouble  xyz_list[NUMVERTICES][3];
-	GaussPenta* gauss=NULL;
-	IssmDouble  vx,vy,vel;
-	IssmDouble  strainparallel;
-	IssmDouble  propcoeff;
-	IssmDouble  strainperpendicular;
-	IssmDouble  calvingratex[NUMVERTICES];
-	IssmDouble  calvingratey[NUMVERTICES];
-	IssmDouble  calvingrate[NUMVERTICES];
-
-
-	/* Get node coordinates and dof list: */
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-
-	/*Retrieve all inputs and parameters we will need*/
-	Input* vx_input=inputs->GetInput(VxEnum);																		_assert_(vx_input);
-	Input* vy_input=inputs->GetInput(VyEnum);																		_assert_(vy_input);
-	Input* strainparallel_input=inputs->GetInput(StrainRateparallelEnum);								_assert_(strainparallel_input);
-	Input* strainperpendicular_input=inputs->GetInput(StrainRateperpendicularEnum);              _assert_(strainperpendicular_input);
-	Input* levermanncoeff_input=inputs->GetInput(CalvinglevermannCoeffEnum);                     _assert_(levermanncoeff_input);
-
-	/* Start looping on the number of vertices: */
-	gauss=new GaussPenta();
-	for (int iv=0;iv<NUMVERTICES;iv++){
-		gauss->GaussVertex(iv);
-
-		/* Get the value we need*/
-		vx_input->GetInputValue(&vx,gauss);
-		vy_input->GetInputValue(&vy,gauss);
-		vel=vx*vx+vy*vy;
-		strainparallel_input->GetInputValue(&strainparallel,gauss);
-		strainperpendicular_input->GetInputValue(&strainperpendicular,gauss);
-		levermanncoeff_input->GetInputValue(&propcoeff,gauss);
-
-		/*Calving rate proportionnal to the positive product of the strain rate along the ice flow direction and the strain rate perpendicular to the ice flow */
-		calvingrate[iv]=propcoeff*strainparallel*strainperpendicular;	
-		if(calvingrate[iv]<0){
-			calvingrate[iv]=0;
-		}
-		calvingratex[iv]=calvingrate[iv]*vx/(vel+1.e-6);
-		calvingratey[iv]=calvingrate[iv]*vy/(vel+1.e-6);
-	}
-
-	/*Add input*/
-	this->inputs->AddInput(new PentaInput(CalvingratexEnum,&calvingratex[0],P1Enum));
-	this->inputs->AddInput(new PentaInput(CalvingrateyEnum,&calvingratey[0],P1Enum));
-	this->inputs->AddInput(new PentaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum));
-
-	/*Clean up and return*/
-	delete gauss;
-
-}
-/*}}}*/
-void       Penta::StressIntensityFactor(){/*{{{*/
-
-	/* Check if we are on the base */
-	if(!IsOnBase()) return;
-
-	IssmDouble  ki[6]={0.};
-	IssmDouble  const_grav=9.81;
-	IssmDouble  rho_ice=900;
-	IssmDouble  rho_water=1000;
-	IssmDouble  Jdet[3];
-	IssmDouble  pressure,vx,vy,vel,deviaxx,deviaxy,deviayy,water_depth,prof,stress_xx,thickness;
-
-	Penta* penta=this;
-	for(;;){
-	
-		IssmDouble  xyz_list[NUMVERTICES][3];
-		/* Get node coordinates and dof list: */
-		::GetVerticesCoordinates(&xyz_list[0][0],penta->vertices,NUMVERTICES);
-
-		///*Compute the Jacobian for the vertical integration*/
-		Jdet[0]=(xyz_list[3][2]-xyz_list[0][2])*0.5;
-		Jdet[1]=(xyz_list[4][2]-xyz_list[1][2])*0.5;
-		Jdet[2]=(xyz_list[5][2]-xyz_list[2][2])*0.5;
-	
-		/*Retrieve all inputs we will need*/
-		Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
-		Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
-		Input* vel_input=inputs->GetInput(VelEnum);                                _assert_(vel_input);
-		Input* pressure_input=inputs->GetInput(PressureEnum);                      _assert_(pressure_input);
-		Input* deviaxx_input=inputs->GetInput(DeviatoricStressxxEnum);             _assert_(deviaxx_input);
-		Input* deviaxy_input=inputs->GetInput(DeviatoricStressxyEnum);             _assert_(deviaxy_input);
-		Input* deviayy_input=inputs->GetInput(DeviatoricStressyyEnum);             _assert_(deviayy_input);
-		Input* surface_input=inputs->GetInput(SurfaceEnum);								_assert_(surface_input);
-		Input* thickness_input=inputs->GetInput(ThicknessEnum);							_assert_(thickness_input);
-		
-		/* Start looping on the number of 2D vertices: */
-		for(int ig=0;ig<3;ig++){
-			GaussPenta* gauss=new GaussPenta(ig,3+ig,11);
-			for (int iv=gauss->begin();iv<gauss->end();iv++){
-				gauss->GaussPoint(iv);
-
-				/* Get the value we need*/
-				pressure_input->GetInputValue(&pressure,gauss);
-				vx_input->GetInputValue(&vx,gauss);
-				vy_input->GetInputValue(&vy,gauss);
-				vel_input->GetInputValue(&vel,gauss);
-				deviaxx_input->GetInputValue(&deviaxx,gauss);
-				deviaxy_input->GetInputValue(&deviaxy,gauss);
-				deviayy_input->GetInputValue(&deviayy,gauss);
-				surface_input->GetInputValue(&water_depth,gauss);
-				thickness_input->GetInputValue(&thickness,gauss);
-				prof=water_depth-penta->GetZcoord(&xyz_list[0][0],gauss);
-
-				/*stress_xx= Deviatoric stress along the ice flow direction plus cryostatic pressure */
-				stress_xx=(vx*vx*(deviaxx)+vy*vy*(deviayy)+2*vy*vx*deviaxy)/(vel*vel+1.e-6);
-
-				if(prof<water_depth&prof<thickness){
-					/* Compute the local stress intensity factor*/ 
-					ki[ig]+=Jdet[ig]*gauss->weight*stress_xx*StressIntensityIntegralWeight(prof,min(water_depth,thickness),thickness);
-				}
-			}
-			delete gauss;
-		}
-			
-		/*Stop if we have reached the surface/base*/
-		if(penta->IsOnSurface()) break;
-		
-		/*get upper Penta*/
-		penta=penta->GetUpperPenta();
-		_assert_(penta->Id()!=this->id);
-	}
-
-	/*Add input*/
-	this->inputs->AddInput(new PentaInput(StressIntensityFactorEnum,&ki[0],P1Enum));
-	this->InputExtrude(StressIntensityFactorEnum,-1);
-}
-/*}}}*/
-void       Penta::StrainRateparallel(){/*{{{*/
-
-	IssmDouble *xyz_list = NULL;
-	IssmDouble  epsilon[6];
-	GaussPenta* gauss=NULL;
-	IssmDouble  vx,vy,vel;
-	IssmDouble  strainxx;
-	IssmDouble  strainxy;
-	IssmDouble  strainyy;
-	IssmDouble  strainparallel[NUMVERTICES];
-
-	/* Get node coordinates and dof list: */
-	this->GetVerticesCoordinates(&xyz_list);
-
-	/*Retrieve all inputs we will need*/
-	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);
-
-	/* Start looping on the number of vertices: */
-	gauss=new GaussPenta();
-	for (int iv=0;iv<NUMVERTICES;iv++){
-		gauss->GaussVertex(iv);
-
-		/* Get the value we need*/
-		vx_input->GetInputValue(&vx,gauss);
-		vy_input->GetInputValue(&vy,gauss);
-		vel=vx*vx+vy*vy;
-
-		/*Compute strain rate and viscosity: */
-		this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
-		strainxx=epsilon[0];
-		strainyy=epsilon[1];
-		strainxy=epsilon[3];
-
-		/*strainparallel= Strain rate along the ice flow direction */
-		strainparallel[iv]=(vx*vx*(strainxx)+vy*vy*(strainyy)+2*vy*vx*strainxy)/(vel+1.e-6);
-	}
-
-	/*Add input*/
-	this->inputs->AddInput(new PentaInput(StrainRateparallelEnum,&strainparallel[0],P1Enum));
-
-	/*Clean up and return*/
-	delete gauss;
-	xDelete<IssmDouble>(xyz_list);
-}
-/*}}}*/
-void       Penta::StrainRateperpendicular(){/*{{{*/
-
-	IssmDouble *xyz_list = NULL;
-	IssmDouble  epsilon[6];
-	GaussPenta* gauss=NULL;
-	IssmDouble  vx,vy,vel;
-	IssmDouble  strainxx;
-	IssmDouble  strainxy;
-	IssmDouble  strainyy;
-	IssmDouble  strainperpendicular[NUMVERTICES];
-
-	/* Get node coordinates and dof list: */
-	this->GetVerticesCoordinates(&xyz_list);
-
-	/*Retrieve all inputs we will need*/
-	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);
-
-	/* Start looping on the number of vertices: */
-	gauss=new GaussPenta();
-	for (int iv=0;iv<NUMVERTICES;iv++){
-		gauss->GaussVertex(iv);
-
-		/* Get the value we need*/
-		vx_input->GetInputValue(&vx,gauss);
-		vy_input->GetInputValue(&vy,gauss);
-		vel=vx*vx+vy*vy;
-
-		/*Compute strain rate and viscosity: */
-		this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
-		strainxx=epsilon[0];
-		strainyy=epsilon[1];
-		strainxy=epsilon[3];
-
-		/*strainperpendicular= Strain rate perpendicular to the ice flow direction */
-		strainperpendicular[iv]=(vx*vx*(strainyy)+vy*vy*(strainxx)-2*vy*vx*strainxy)/(vel+1.e-6);
-	}
-
-	/*Add input*/
-	this->inputs->AddInput(new PentaInput(StrainRateperpendicularEnum,&strainperpendicular[0],P1Enum));
-
-	/*Clean up and return*/
-	delete gauss;
-	xDelete<IssmDouble>(xyz_list);
-}
-/*}}}*/
 void       Penta::Configure(Elements* elementsin, Loads* loadsin, Nodes* nodesin,Vertices* verticesin, Materials* materialsin, Parameters* parametersin){/*{{{*/
 
@@ -606,4 +439,70 @@
 }
 /*}}}*/
+void       Penta::ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){/*{{{*/
+
+	int    vertexpidlist[NUMVERTICES];
+	IssmDouble grad_list[NUMVERTICES];
+	Input* grad_input=NULL;
+	Input* input=NULL;
+
+	if(enum_type==MaterialsRheologyBbarEnum){
+		input=(Input*)inputs->GetInput(MaterialsRheologyBEnum);
+	}
+	else if(enum_type==DamageDbarEnum){
+		input=(Input*)inputs->GetInput(DamageDEnum);
+	}
+	else{
+		input=inputs->GetInput(enum_type);
+	}
+	if (!input) _error_("Input " << EnumToStringx(enum_type) << " not found");
+	if (input->ObjectEnum()!=ControlInputEnum) _error_("Input " << EnumToStringx(enum_type) << " is not a ControlInput");
+
+	GradientIndexing(&vertexpidlist[0],control_index);
+	for(int i=0;i<NUMVERTICES;i++) grad_list[i]=gradient[vertexpidlist[i]];
+	grad_input=new PentaInput(GradientEnum,grad_list,P1Enum);
+	((ControlInput*)input)->SetGradient(grad_input);
+
+}/*}}}*/
+void       Penta::ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum){/*{{{*/
+
+	Input* input=NULL;
+
+	if(control_enum==MaterialsRheologyBbarEnum){
+		input=(Input*)inputs->GetInput(MaterialsRheologyBEnum);
+	}
+	else if(control_enum==DamageDbarEnum){
+		input=(Input*)inputs->GetInput(DamageDEnum);
+	}
+	else{
+		input=inputs->GetInput(control_enum);
+	}
+	if (!input) _error_("Input " << EnumToStringx(control_enum) << " not found");
+	if (input->ObjectEnum()!=ControlInputEnum) _error_("Input " << EnumToStringx(control_enum) << " is not a ControlInput");
+
+	int         sidlist[NUMVERTICES];
+	int         connectivity[NUMVERTICES];
+	IssmPDouble values[NUMVERTICES];
+	IssmPDouble gradients[NUMVERTICES]; 
+	IssmDouble  value,gradient;
+
+	this->GetVerticesConnectivityList(&connectivity[0]);
+	this->GetVerticesSidList(&sidlist[0]);
+
+	GaussPenta* gauss=new GaussPenta();
+	for (int iv=0;iv<NUMVERTICES;iv++){
+		gauss->GaussVertex(iv);
+
+		((ControlInput*)input)->GetInputValue(&value,gauss);
+		((ControlInput*)input)->GetGradientValue(&gradient,gauss);
+
+		values[iv]    = reCast<IssmPDouble>(value)/reCast<IssmPDouble>(connectivity[iv]);
+		gradients[iv] = reCast<IssmPDouble>(gradient)/reCast<IssmPDouble>(connectivity[iv]);
+	}
+	delete gauss;
+
+	vector_control->SetValues(NUMVERTICES,&sidlist[0],&values[0],ADD_VAL);
+	vector_gradient->SetValues(NUMVERTICES,&sidlist[0],&gradients[0],ADD_VAL);
+
+}/*}}}*/
 void       Penta::Delta18oParameterization(void){/*{{{*/
 
@@ -680,4 +579,59 @@
 }
 /*}}}*/
+void       Penta::ElementResponse(IssmDouble* presponse,int response_enum){/*{{{*/
+
+	switch(response_enum){
+		case MaterialsRheologyBbarEnum:
+			*presponse=this->material->GetBbar();
+			break;
+		case DamageDbarEnum:
+			*presponse=this->material->GetDbar();
+			break;
+		case VelEnum:
+			{
+
+				/*Get input:*/
+				IssmDouble vel;
+				Input* vel_input;
+
+				vel_input=this->inputs->GetInput(VelEnum); _assert_(vel_input);
+				vel_input->GetInputAverage(&vel);
+
+				/*Assign output pointers:*/
+				*presponse=vel;
+			}
+			break;
+		default:  
+			_error_("Response type " << EnumToStringx(response_enum) << " not supported yet!");
+	}
+
+}
+/*}}}*/
+void       Penta::ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){/*{{{*/
+
+	IssmDouble xyz_list[NUMVERTICES][3];
+	IssmDouble xmin,ymin,zmin;
+	IssmDouble xmax,ymax,zmax;
+
+	/*Get xyz list: */
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+	xmin=xyz_list[0][0]; xmax=xyz_list[0][0];
+	ymin=xyz_list[0][1]; ymax=xyz_list[0][1];
+	zmin=xyz_list[0][2]; zmax=xyz_list[0][2];
+
+	for(int i=1;i<NUMVERTICES;i++){
+		if(xyz_list[i][0]<xmin) xmin=xyz_list[i][0];
+		if(xyz_list[i][0]>xmax) xmax=xyz_list[i][0];
+		if(xyz_list[i][1]<ymin) ymin=xyz_list[i][1];
+		if(xyz_list[i][1]>ymax) ymax=xyz_list[i][1];
+		if(xyz_list[i][2]<zmin) zmin=xyz_list[i][2];
+		if(xyz_list[i][2]>zmax) zmax=xyz_list[i][2];
+	}
+
+	*hx=xmax-xmin;
+	*hy=ymax-ymin;
+	*hz=zmax-zmin;
+}
+/*}}}*/
 int        Penta::FiniteElement(void){/*{{{*/
 	return this->element_type;
@@ -770,10 +724,4 @@
 }
 /*}}}*/
-int        Penta::ObjectEnum(void){/*{{{*/
-
-	return PentaEnum;
-
-}
-/*}}}*/
 void       Penta::GetAreaCoordinates(IssmDouble* area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints){/*{{{*/
 	/*Computeportion of the element that is grounded*/ 
@@ -811,40 +759,9 @@
 }
 /*}}}*/
-Penta*     Penta::GetUpperPenta(void){/*{{{*/
-
-	Penta* upper_penta=NULL;
-
-	upper_penta=(Penta*)verticalneighbors[1]; //first one (0) under, second one (1) above
-
-	return upper_penta;
-}
-/*}}}*/
-Penta*     Penta::GetLowerPenta(void){/*{{{*/
-
-	Penta* lower_penta=NULL;
-
-	lower_penta=(Penta*)verticalneighbors[0]; //first one (0) under, second one (1) above
-
-	return lower_penta;
-}
-/*}}}*/
-Penta*     Penta::GetSurfacePenta(void){/*{{{*/
+Element*   Penta::GetBasalElement(void){/*{{{*/
 
 	/*Output*/
-	Penta* penta=NULL;
-
-	/*Go through all pentas till the surface is reached*/
-	penta=this;
-	for(;;){
-		/*Stop if we have reached the surface, else, take upper penta*/
-		if (penta->IsOnSurface()) break;
-
-		/* get upper Penta*/
-		penta=penta->GetUpperPenta();
-		_assert_(penta->Id()!=this->id);
-	}
-
-	/*return output*/
-	return penta;
+	Element* element=this->GetBasalPenta();
+	return element;
 }
 /*}}}*/
@@ -869,16 +786,8 @@
 }
 /*}}}*/
-Element*   Penta::GetUpperElement(void){/*{{{*/
-
-	/*Output*/
-	Element* upper_element=this->GetUpperPenta();
-	return upper_element;
-}
-/*}}}*/
-Element*   Penta::GetBasalElement(void){/*{{{*/
-
-	/*Output*/
-	Element* element=this->GetBasalPenta();
-	return element;
+int        Penta::GetElementType(){/*{{{*/
+
+	/*return PentaRef field*/
+	return this->element_type;
 }
 /*}}}*/
@@ -1035,288 +944,4 @@
 
 	return phi;
-}
-/*}}}*/
-int        Penta::GetElementType(){/*{{{*/
-
-	/*return PentaRef field*/
-	return this->element_type;
-}
-/*}}}*/
-void       Penta::ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){/*{{{*/
-
-	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble xmin,ymin,zmin;
-	IssmDouble xmax,ymax,zmax;
-
-	/*Get xyz list: */
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	xmin=xyz_list[0][0]; xmax=xyz_list[0][0];
-	ymin=xyz_list[0][1]; ymax=xyz_list[0][1];
-	zmin=xyz_list[0][2]; zmax=xyz_list[0][2];
-
-	for(int i=1;i<NUMVERTICES;i++){
-		if(xyz_list[i][0]<xmin) xmin=xyz_list[i][0];
-		if(xyz_list[i][0]>xmax) xmax=xyz_list[i][0];
-		if(xyz_list[i][1]<ymin) ymin=xyz_list[i][1];
-		if(xyz_list[i][1]>ymax) ymax=xyz_list[i][1];
-		if(xyz_list[i][2]<zmin) zmin=xyz_list[i][2];
-		if(xyz_list[i][2]>zmax) zmax=xyz_list[i][2];
-	}
-
-	*hx=xmax-xmin;
-	*hy=ymax-ymin;
-	*hz=zmax-zmin;
-}
-/*}}}*/
-int        Penta::GetNodeIndex(Node* node){/*{{{*/
-
-	_assert_(nodes);
-	int numnodes = this->NumberofNodes(this->element_type);
-
-	for(int i=0;i<numnodes;i++){
-		if(node==nodes[i]) return i;
-	}
-	_error_("Node provided not found among element nodes");
-
-}
-/*}}}*/
-int        Penta::GetNumberOfNodes(void){/*{{{*/
-	return this->NumberofNodes(this->element_type);
-}
-/*}}}*/
-int        Penta::GetNumberOfNodes(int enum_type){/*{{{*/
-	return this->NumberofNodes(enum_type);
-}
-/*}}}*/
-int        Penta::GetNumberOfVertices(void){/*{{{*/
-	return NUMVERTICES; 
-}
-/*}}}*/
-Node*      Penta::GetNode(int node_number){/*{{{*/
-	_assert_(node_number>=0); 
-	_assert_(node_number<this->NumberofNodes(this->element_type)); 
-	return this->nodes[node_number];
-}
-/*}}}*/
-void       Penta::GetInputValue(IssmDouble* pvalue,Node* node,int enumtype){/*{{{*/
-
-	Input* input=inputs->GetInput(enumtype);
-	if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria");
-
-	GaussPenta* gauss=new GaussPenta();
-	gauss->GaussVertex(this->GetNodeIndex(node));
-
-	input->GetInputValue(pvalue,gauss);
-	delete gauss;
-}
-/*}}}*/
-void       Penta::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){/*{{{*/
-
-	IssmDouble* xyz_list = xNew<IssmDouble>(NUMVERTICES2D*3);
-	::GetVerticesCoordinates(xyz_list,this->vertices,NUMVERTICES2D);
-
-	/*Assign output pointer*/
-	*pxyz_list = xyz_list;
-
-}/*}}}*/
-void       Penta::GetVerticesCoordinatesTop(IssmDouble** pxyz_list){/*{{{*/
-
-	IssmDouble* xyz_list = xNew<IssmDouble>(NUMVERTICES2D*3);
-	::GetVerticesCoordinates(xyz_list,&this->vertices[3],NUMVERTICES2D);
-
-	/*Assign output pointer*/
-	*pxyz_list = xyz_list;
-
-}/*}}}*/
-void       Penta::NormalSection(IssmDouble* normal,IssmDouble* xyz_list){/*{{{*/
-
-	/*Build unit outward pointing vector*/
-	IssmDouble AB[3];
-	IssmDouble AC[3];
-	IssmDouble norm;
-
-	AB[0]=xyz_list[1*3+0] - xyz_list[0*3+0];
-	AB[1]=xyz_list[1*3+1] - xyz_list[0*3+1];
-	AB[2]=xyz_list[1*3+2] - xyz_list[0*3+2];
-	AC[0]=xyz_list[2*3+0] - xyz_list[0*3+0];
-	AC[1]=xyz_list[2*3+1] - xyz_list[0*3+1];
-	AC[2]=xyz_list[2*3+2] - xyz_list[0*3+2];
-
-	cross(normal,AB,AC);
-	norm=sqrt(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]);
-
-	for(int i=0;i<3;i++) normal[i]=normal[i]/norm;
-}
-/*}}}*/
-IssmDouble Penta::StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){/*{{{*/
-	/*Compute stabilization parameter*/
-	/*kappa=thermalconductivity/(rho_ice*hearcapacity) for thermal model*/
-	/*kappa=enthalpydiffusionparameter for enthalpy model*/
-
-	IssmDouble normu;
-	IssmDouble tau_parameter;
-
-	normu=pow(pow(u,2)+pow(v,2)+pow(w,2),0.5);
-	if(normu*diameter/(3*2*kappa)<1){ 
-		tau_parameter=pow(diameter,2)/(3*2*2*kappa);
-	}
-	else tau_parameter=diameter/(2*normu);
-
-	return tau_parameter;
-}
-/*}}}*/
-void       Penta::ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
-	/*Compute portion of the element that is grounded*/ 
-
-	int         normal_orientation=0;
-	IssmDouble  s1,s2;
-	IssmDouble  levelset[NUMVERTICES];
-	IssmDouble* xyz_zero = xNew<IssmDouble>(4*3);
-
-	/*Recover parameters and values*/
-	GetInputListOnVertices(&levelset[0],levelsetenum);
-
-	if(levelset[0]*levelset[1]>0.){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
-		/*Portion of the segments*/
-		s1=levelset[2]/(levelset[2]-levelset[1]);
-		s2=levelset[2]/(levelset[2]-levelset[0]);
-
-		if(levelset[2]<0.) normal_orientation=1; //orientation of quadrangle at base and top, depending on distribution of levelsetfunction
-		/*New point 1*/
-		xyz_zero[3*normal_orientation+0]=xyz_list[2*3+0]+s1*(xyz_list[1*3+0]-xyz_list[2*3+0]);
-		xyz_zero[3*normal_orientation+1]=xyz_list[2*3+1]+s1*(xyz_list[1*3+1]-xyz_list[2*3+1]);
-		xyz_zero[3*normal_orientation+2]=xyz_list[2*3+2]+s1*(xyz_list[1*3+2]-xyz_list[2*3+2]);
-
-		/*New point 0*/
-		xyz_zero[3*(1-normal_orientation)+0]=xyz_list[2*3+0]+s2*(xyz_list[0*3+0]-xyz_list[2*3+0]);
-		xyz_zero[3*(1-normal_orientation)+1]=xyz_list[2*3+1]+s2*(xyz_list[0*3+1]-xyz_list[2*3+1]);
-		xyz_zero[3*(1-normal_orientation)+2]=xyz_list[2*3+2]+s2*(xyz_list[0*3+2]-xyz_list[2*3+2]);
-
-		/*New point 3*/
-		xyz_zero[3*(2+1-normal_orientation)+0]=xyz_list[5*3+0]+s1*(xyz_list[4*3+0]-xyz_list[5*3+0]);
-		xyz_zero[3*(2+1-normal_orientation)+1]=xyz_list[5*3+1]+s1*(xyz_list[4*3+1]-xyz_list[5*3+1]);
-		xyz_zero[3*(2+1-normal_orientation)+2]=xyz_list[5*3+2]+s1*(xyz_list[4*3+2]-xyz_list[5*3+2]);
-
-		/*New point 4*/
-		xyz_zero[3*(2+normal_orientation)+0]=xyz_list[5*3+0]+s2*(xyz_list[3*3+0]-xyz_list[5*3+0]);
-		xyz_zero[3*(2+normal_orientation)+1]=xyz_list[5*3+1]+s2*(xyz_list[3*3+1]-xyz_list[5*3+1]);
-		xyz_zero[3*(2+normal_orientation)+2]=xyz_list[5*3+2]+s2*(xyz_list[3*3+2]-xyz_list[5*3+2]);
-	}
-	else if(levelset[1]*levelset[2]>0.){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
-		/*Portion of the segments*/
-		s1=levelset[0]/(levelset[0]-levelset[2]);
-		s2=levelset[0]/(levelset[0]-levelset[1]);
-
-		if(levelset[0]<0.) normal_orientation=1;
-		/*New point 1*/
-		xyz_zero[3*normal_orientation+0]=xyz_list[0*3+0]+s1*(xyz_list[2*3+0]-xyz_list[0*3+0]);
-		xyz_zero[3*normal_orientation+1]=xyz_list[0*3+1]+s1*(xyz_list[2*3+1]-xyz_list[0*3+1]);
-		xyz_zero[3*normal_orientation+2]=xyz_list[0*3+2]+s1*(xyz_list[2*3+2]-xyz_list[0*3+2]);
-
-		/*New point 2*/
-		xyz_zero[3*(1-normal_orientation)+0]=xyz_list[0*3+0]+s2*(xyz_list[1*3+0]-xyz_list[0*3+0]);
-		xyz_zero[3*(1-normal_orientation)+1]=xyz_list[0*3+1]+s2*(xyz_list[1*3+1]-xyz_list[0*3+1]);
-		xyz_zero[3*(1-normal_orientation)+2]=xyz_list[0*3+2]+s2*(xyz_list[1*3+2]-xyz_list[0*3+2]);
-
-		/*New point 3*/
-		xyz_zero[3*(2+1-normal_orientation)+0]=xyz_list[3*3+0]+s1*(xyz_list[5*3+0]-xyz_list[3*3+0]);
-		xyz_zero[3*(2+1-normal_orientation)+1]=xyz_list[3*3+1]+s1*(xyz_list[5*3+1]-xyz_list[3*3+1]);
-		xyz_zero[3*(2+1-normal_orientation)+2]=xyz_list[3*3+2]+s1*(xyz_list[5*3+2]-xyz_list[3*3+2]);
-
-		/*New point 4*/
-		xyz_zero[3*(2+normal_orientation)+0]=xyz_list[3*3+0]+s2*(xyz_list[4*3+0]-xyz_list[3*3+0]);
-		xyz_zero[3*(2+normal_orientation)+1]=xyz_list[3*3+1]+s2*(xyz_list[4*3+1]-xyz_list[3*3+1]);
-		xyz_zero[3*(2+normal_orientation)+2]=xyz_list[3*3+2]+s2*(xyz_list[4*3+2]-xyz_list[3*3+2]);
-	}
-	else if(levelset[0]*levelset[2]>0.){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
-		/*Portion of the segments*/
-		s1=levelset[1]/(levelset[1]-levelset[0]);
-		s2=levelset[1]/(levelset[1]-levelset[2]);
-
-		if(levelset[1]<0.) normal_orientation=1;
-		/*New point 0*/
-		xyz_zero[3*normal_orientation+0]=xyz_list[1*3+0]+s1*(xyz_list[0*3+0]-xyz_list[1*3+0]);
-		xyz_zero[3*normal_orientation+1]=xyz_list[1*3+1]+s1*(xyz_list[0*3+1]-xyz_list[1*3+1]);
-		xyz_zero[3*normal_orientation+2]=xyz_list[1*3+2]+s1*(xyz_list[0*3+2]-xyz_list[1*3+2]);
-
-		/*New point 2*/
-		xyz_zero[3*(1-normal_orientation)+0]=xyz_list[1*3+0]+s2*(xyz_list[2*3+0]-xyz_list[1*3+0]);
-		xyz_zero[3*(1-normal_orientation)+1]=xyz_list[1*3+1]+s2*(xyz_list[2*3+1]-xyz_list[1*3+1]);
-		xyz_zero[3*(1-normal_orientation)+2]=xyz_list[1*3+2]+s2*(xyz_list[2*3+2]-xyz_list[1*3+2]);
-
-		/*New point 3*/
-		xyz_zero[3*(2+1-normal_orientation)+0]=xyz_list[4*3+0]+s1*(xyz_list[3*3+0]-xyz_list[4*3+0]);
-		xyz_zero[3*(2+1-normal_orientation)+1]=xyz_list[4*3+1]+s1*(xyz_list[3*3+1]-xyz_list[4*3+1]);
-		xyz_zero[3*(2+1-normal_orientation)+2]=xyz_list[4*3+2]+s1*(xyz_list[3*3+2]-xyz_list[4*3+2]);
-
-		/*New point 4*/
-		xyz_zero[3*(2+normal_orientation)+0]=xyz_list[4*3+0]+s2*(xyz_list[5*3+0]-xyz_list[4*3+0]);
-		xyz_zero[3*(2+normal_orientation)+1]=xyz_list[4*3+1]+s2*(xyz_list[5*3+1]-xyz_list[4*3+1]);
-		xyz_zero[3*(2+normal_orientation)+2]=xyz_list[4*3+2]+s2*(xyz_list[5*3+2]-xyz_list[4*3+2]);
-	}
-	else if(levelset[0]==0. && levelset[1]==0.){ //front is on point 0 and 1
-		xyz_zero[3*0+0]=xyz_list[0*3+0];
-		xyz_zero[3*0+1]=xyz_list[0*3+1];
-		xyz_zero[3*0+2]=xyz_list[0*3+2];
-
-		/*New point 2*/
-		xyz_zero[3*1+0]=xyz_list[1*3+0];
-		xyz_zero[3*1+1]=xyz_list[1*3+1];
-		xyz_zero[3*1+2]=xyz_list[1*3+2];
-
-		/*New point 3*/
-		xyz_zero[3*2+0]=xyz_list[4*3+0];
-		xyz_zero[3*2+1]=xyz_list[4*3+1];
-		xyz_zero[3*2+2]=xyz_list[4*3+2];
-
-		/*New point 4*/
-		xyz_zero[3*3+0]=xyz_list[3*3+0];
-		xyz_zero[3*3+1]=xyz_list[3*3+1];
-		xyz_zero[3*3+2]=xyz_list[3*3+2];
-	}
-	else if(levelset[0]==0. && levelset[2]==0.){ //front is on point 0 and 1
-		xyz_zero[3*0+0]=xyz_list[2*3+0];
-		xyz_zero[3*0+1]=xyz_list[2*3+1];
-		xyz_zero[3*0+2]=xyz_list[2*3+2];
-
-		/*New point 2*/
-		xyz_zero[3*1+0]=xyz_list[0*3+0];
-		xyz_zero[3*1+1]=xyz_list[0*3+1];
-		xyz_zero[3*1+2]=xyz_list[0*3+2];
-
-		/*New point 3*/
-		xyz_zero[3*2+0]=xyz_list[3*3+0];
-		xyz_zero[3*2+1]=xyz_list[3*3+1];
-		xyz_zero[3*2+2]=xyz_list[3*3+2];
-
-		/*New point 4*/
-		xyz_zero[3*3+0]=xyz_list[5*3+0];
-		xyz_zero[3*3+1]=xyz_list[5*3+1];
-		xyz_zero[3*3+2]=xyz_list[5*3+2];
-	}
-	else if(levelset[1]==0. && levelset[2]==0.){ //front is on point 0 and 1
-		xyz_zero[3*0+0]=xyz_list[1*3+0];
-		xyz_zero[3*0+1]=xyz_list[1*3+1];
-		xyz_zero[3*0+2]=xyz_list[1*3+2];
-
-		/*New point 2*/
-		xyz_zero[3*1+0]=xyz_list[2*3+0];
-		xyz_zero[3*1+1]=xyz_list[2*3+1];
-		xyz_zero[3*1+2]=xyz_list[2*3+2];
-
-		/*New point 3*/
-		xyz_zero[3*2+0]=xyz_list[5*3+0];
-		xyz_zero[3*2+1]=xyz_list[5*3+1];
-		xyz_zero[3*2+2]=xyz_list[5*3+2];
-
-		/*New point 4*/
-		xyz_zero[3*3+0]=xyz_list[4*3+0];
-		xyz_zero[3*3+1]=xyz_list[4*3+1];
-		xyz_zero[3*3+2]=xyz_list[4*3+2];
-	}
-	else _error_("Case not covered");
-
-	/*Assign output pointer*/
-	*pxyz_zero= xyz_zero;
 }
 /*}}}*/
@@ -1363,4 +988,257 @@
 	xDelete<int>(indicesfront);
 }/*}}}*/
+void       Penta::GetInputValue(IssmDouble* pvalue,Node* node,int enumtype){/*{{{*/
+
+	Input* input=inputs->GetInput(enumtype);
+	if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria");
+
+	GaussPenta* gauss=new GaussPenta();
+	gauss->GaussVertex(this->GetNodeIndex(node));
+
+	input->GetInputValue(pvalue,gauss);
+	delete gauss;
+}
+/*}}}*/
+Penta*     Penta::GetLowerPenta(void){/*{{{*/
+
+	Penta* lower_penta=NULL;
+
+	lower_penta=(Penta*)verticalneighbors[0]; //first one (0) under, second one (1) above
+
+	return lower_penta;
+}
+/*}}}*/
+Node*      Penta::GetNode(int node_number){/*{{{*/
+	_assert_(node_number>=0); 
+	_assert_(node_number<this->NumberofNodes(this->element_type)); 
+	return this->nodes[node_number];
+}
+/*}}}*/
+int        Penta::GetNodeIndex(Node* node){/*{{{*/
+
+	_assert_(nodes);
+	int numnodes = this->NumberofNodes(this->element_type);
+
+	for(int i=0;i<numnodes;i++){
+		if(node==nodes[i]) return i;
+	}
+	_error_("Node provided not found among element nodes");
+
+}
+/*}}}*/
+int        Penta::GetNumberOfNodes(void){/*{{{*/
+	return this->NumberofNodes(this->element_type);
+}
+/*}}}*/
+int        Penta::GetNumberOfNodes(int enum_type){/*{{{*/
+	return this->NumberofNodes(enum_type);
+}
+/*}}}*/
+int        Penta::GetNumberOfVertices(void){/*{{{*/
+	return NUMVERTICES; 
+}
+/*}}}*/
+void       Penta::GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution, int enum_type){/*{{{*/
+
+	const int    numdof=NDOF1*NUMVERTICES;
+
+	int          i;
+	int*         doflist=NULL;
+	IssmDouble   values[numdof];
+	IssmDouble   enum_value;
+	GaussPenta   *gauss=NULL;
+
+	/*Get dof list: */
+	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
+	Input* enum_input=inputs->GetInput(enum_type); _assert_(enum_input);
+
+	gauss=new GaussPenta();
+	for(i=0;i<NUMVERTICES;i++){
+		/*Recover temperature*/
+		gauss->GaussVertex(i);
+		enum_input->GetInputValue(&enum_value,gauss);
+		values[i]=enum_value;
+	}
+
+	/*Add value to global vector*/
+	solution->SetValues(numdof,doflist,values,INS_VAL);
+
+	/*Free ressources:*/
+	delete gauss;
+	xDelete<int>(doflist);
+}
+/*}}}*/
+Penta*     Penta::GetSurfacePenta(void){/*{{{*/
+
+	/*Output*/
+	Penta* penta=NULL;
+
+	/*Go through all pentas till the surface is reached*/
+	penta=this;
+	for(;;){
+		/*Stop if we have reached the surface, else, take upper penta*/
+		if (penta->IsOnSurface()) break;
+
+		/* get upper Penta*/
+		penta=penta->GetUpperPenta();
+		_assert_(penta->Id()!=this->id);
+	}
+
+	/*return output*/
+	return penta;
+}
+/*}}}*/
+Element*   Penta::GetUpperElement(void){/*{{{*/
+
+	/*Output*/
+	Element* upper_element=this->GetUpperPenta();
+	return upper_element;
+}
+/*}}}*/
+Penta*     Penta::GetUpperPenta(void){/*{{{*/
+
+	Penta* upper_penta=NULL;
+
+	upper_penta=(Penta*)verticalneighbors[1]; //first one (0) under, second one (1) above
+
+	return upper_penta;
+}
+/*}}}*/
+void       Penta::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,const char* data,bool onsid){/*{{{*/
+
+	int vertexidlist[NUMVERTICES];
+
+	/*Get out if this is not an element input*/
+	if(!IsInput(control_enum)) return;
+
+	/*Prepare index list*/
+	GradientIndexing(&vertexidlist[0],control_index,onsid);
+
+	/*Get input (either in element or material)*/
+	if(control_enum==MaterialsRheologyBbarEnum) control_enum=MaterialsRheologyBEnum;
+	Input* input=inputs->GetInput(control_enum);
+	if(!input) _error_("Input " << EnumToStringx(control_enum) << " not found in element");
+
+	/*Check that it is a ControlInput*/
+	if (input->ObjectEnum()!=ControlInputEnum){
+		_error_("input " << EnumToStringx(control_enum) << " is not a ControlInput");
+	}
+
+	((ControlInput*)input)->GetVectorFromInputs(vector,&vertexidlist[0],data);
+}
+/*}}}*/
+void       Penta::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){/*{{{*/
+
+	IssmDouble* xyz_list = xNew<IssmDouble>(NUMVERTICES2D*3);
+	::GetVerticesCoordinates(xyz_list,this->vertices,NUMVERTICES2D);
+
+	/*Assign output pointer*/
+	*pxyz_list = xyz_list;
+
+}/*}}}*/
+void       Penta::GetVerticesCoordinatesTop(IssmDouble** pxyz_list){/*{{{*/
+
+	IssmDouble* xyz_list = xNew<IssmDouble>(NUMVERTICES2D*3);
+	::GetVerticesCoordinates(xyz_list,&this->vertices[3],NUMVERTICES2D);
+
+	/*Assign output pointer*/
+	*pxyz_list = xyz_list;
+
+}/*}}}*/
+IssmDouble Penta::IceVolume(void){/*{{{*/
+
+	/*The volume of a troncated prism is base * 1/3 sum(length of edges)*/
+	IssmDouble base,height;
+	IssmDouble xyz_list[NUMVERTICES][3];
+
+	if(!IsIceInElement())return 0;
+
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+
+	/*First calculate the area of the base (cross section triangle)
+	 * http://en.wikipedia.org/wiki/Pentangle
+	 * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
+	base = 1./2.*fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
+
+	/*Now get the average height*/
+	height = 1./3.*((xyz_list[3][2]-xyz_list[0][2])+(xyz_list[4][2]-xyz_list[1][2])+(xyz_list[5][2]-xyz_list[2][2]));
+
+	/*Return: */
+	return base*height;
+}
+/*}}}*/
+IssmDouble Penta::IceVolumeAboveFloatation(void){/*{{{*/
+
+	/*Volume above floatation: H + rho_water/rho_ice*bathymetry for nodes on the bed*/
+	IssmDouble rho_ice,rho_water;
+	IssmDouble base,bed,surface,bathymetry;
+	IssmDouble xyz_list[NUMVERTICES][3];
+
+	if(!IsIceInElement() || IsFloating() || !IsOnBase())return 0;
+
+	rho_ice=matpar->GetRhoIce();
+	rho_water=matpar->GetRhoWater();
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+
+	/*First calculate the area of the base (cross section triangle)
+	 * http://en.wikipedia.org/wiki/Pentangle
+	 * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
+	base = 1./2.*fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
+
+	/*Now get the average height above floatation*/
+	Input* surface_input    = inputs->GetInput(SurfaceEnum);    _assert_(surface_input);
+	Input* base_input        = inputs->GetInput(BaseEnum);        _assert_(base_input);
+	Input* bed_input = inputs->GetInput(BedEnum); _assert_(bed_input);
+	surface_input->GetInputAverage(&surface);
+	base_input->GetInputAverage(&bed);
+	bed_input->GetInputAverage(&bathymetry);
+
+	/*Return: */
+	return base*(surface - bed + min( rho_water/rho_ice * bathymetry, 0.) );
+}
+/*}}}*/
+void       Penta::InputControlUpdate(IssmDouble scalar,bool save_parameter){/*{{{*/
+
+	/*Intermediary*/
+	int    num_controls;
+	int*   control_type=NULL;
+	Input* input=NULL;
+
+	/*retrieve some parameters: */
+	this->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
+	this->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
+
+	for(int i=0;i<num_controls;i++){
+
+		if(control_type[i]==MaterialsRheologyBbarEnum){
+			if (!IsOnBase()) goto cleanup_and_return;
+			input=(Input*)this->inputs->GetInput(MaterialsRheologyBEnum); _assert_(input);
+		}
+		else if(control_type[i]==DamageDbarEnum){
+			if (!IsOnBase()) goto cleanup_and_return;
+			input=(Input*)this->inputs->GetInput(DamageDEnum); _assert_(input);
+		}
+		else{
+			input=(Input*)this->inputs->GetInput(control_type[i]); _assert_(input);
+		}
+		if(input->ObjectEnum()!=ControlInputEnum) _error_("input " << EnumToStringx(control_type[i]) << " is not a ControlInput");
+
+		((ControlInput*)input)->UpdateValue(scalar);
+		((ControlInput*)input)->Constrain();
+		if (save_parameter) ((ControlInput*)input)->SaveValue();
+
+		if(control_type[i]==MaterialsRheologyBbarEnum){
+			this->InputExtrude(MaterialsRheologyBEnum,-1);
+		}
+		else if(control_type[i]==DamageDbarEnum){
+			this->InputExtrude(DamageDEnum,-1);
+		}
+	}
+
+	/*Clean up and return*/
+cleanup_and_return:
+	xDelete<int>(control_type);
+}
+/*}}}*/
 void       Penta::InputDepthAverageAtBase(int enum_type,int average_enum_type){/*{{{*/
 
@@ -1732,4 +1610,37 @@
 }
 /*}}}*/
+bool       Penta::IsIcefront(void){/*{{{*/
+
+	bool isicefront;
+	int i,nrice;
+   IssmDouble ls[NUMVERTICES];
+
+	/*Retrieve all inputs and parameters*/
+	GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
+
+	/* If only one vertex has ice, there is an ice front here */
+	isicefront=false;
+	if(IsIceInElement()){
+		nrice=0;       
+		for(i=0;i<NUMVERTICES2D;i++)
+			if(ls[i]<0.) nrice++;
+		if(nrice==1) isicefront= true;
+	}
+	return isicefront;
+}/*}}}*/
+bool       Penta::IsNodeOnShelfFromFlags(IssmDouble* flags){/*{{{*/
+
+	int  i;
+	bool shelf=false;
+
+	for(i=0;i<NUMVERTICES;i++){
+		if (flags[vertices[i]->Pid()]<0.){
+			shelf=true;
+			break;
+		}
+	}
+	return shelf;
+}
+/*}}}*/
 bool       Penta::IsOnBase(void){/*{{{*/
 
@@ -1768,16 +1679,20 @@
 }
 /*}}}*/
-bool       Penta::IsNodeOnShelfFromFlags(IssmDouble* flags){/*{{{*/
-
-	int  i;
-	bool shelf=false;
-
-	for(i=0;i<NUMVERTICES;i++){
-		if (flags[vertices[i]->Pid()]<0.){
-			shelf=true;
-			break;
-		}
-	}
-	return shelf;
+bool       Penta::IsZeroLevelset(int levelset_enum){/*{{{*/
+
+	bool        iszerols;
+	IssmDouble  ls[NUMVERTICES];
+
+	/*Retrieve all inputs and parameters*/
+	GetInputListOnVertices(&ls[0],levelset_enum);
+
+	/*If the level set has always same sign, there is no ice front here*/
+	iszerols = false;
+	if(IsIceInElement()){
+		if(ls[0]*ls[1]<0. || ls[0]*ls[2]<0. || (ls[0]*ls[1]*ls[2]==0. && ls[0]*ls[1]+ls[0]*ls[2]+ls[1]*ls[2]<=0.)){
+			iszerols = true;
+		}
+	}
+	return iszerols;
 }
 /*}}}*/
@@ -1803,4 +1718,11 @@
 }
 /*}}}*/
+void       Penta::JacobianDeterminantSurface(IssmDouble* pJdet,IssmDouble* xyz_list_quad,Gauss* gauss){/*{{{*/
+
+	_assert_(gauss->Enum()==GaussPentaEnum);
+	this->GetQuadJacobianDeterminant(pJdet,xyz_list_quad,(GaussPenta*)gauss);
+
+}
+/*}}}*/
 void       Penta::JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_top,Gauss* gauss){/*{{{*/
 
@@ -1810,9 +1732,48 @@
 }
 /*}}}*/
-void       Penta::JacobianDeterminantSurface(IssmDouble* pJdet,IssmDouble* xyz_list_quad,Gauss* gauss){/*{{{*/
-
-	_assert_(gauss->Enum()==GaussPentaEnum);
-	this->GetQuadJacobianDeterminant(pJdet,xyz_list_quad,(GaussPenta*)gauss);
-
+IssmDouble Penta::MassFlux(IssmDouble* segment){/*{{{*/
+
+	IssmDouble mass_flux=0;
+
+	if(!IsOnBase()) return mass_flux;
+
+	/*Depth Averaging Vx and Vy*/
+	this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
+	this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
+
+	/*Spawn Tria element from the base of the Penta: */
+	Tria* tria=(Tria*)SpawnTria(0,1,2);
+	mass_flux=tria->MassFlux(segment);
+	delete tria->material; delete tria;
+
+	/*Delete Vx and Vy averaged*/
+	this->inputs->DeleteInput(VxAverageEnum);
+	this->inputs->DeleteInput(VyAverageEnum);
+
+	/*clean up and return*/
+	return mass_flux;
+}
+/*}}}*/
+IssmDouble Penta::MassFlux(IssmDouble x1, IssmDouble y1, IssmDouble x2, IssmDouble y2,int segment_id){/*{{{*/
+
+	IssmDouble mass_flux=0;
+
+	if(!IsOnBase()) return mass_flux;
+
+	/*Depth Averaging Vx and Vy*/
+	this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
+	this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
+
+	/*Spawn Tria element from the base of the Penta: */
+	Tria* tria=(Tria*)SpawnTria(0,1,2);
+	mass_flux=tria->MassFlux(x1,y1,x2,y2,segment_id);
+	delete tria->material; delete tria;
+
+	/*Delete Vx and Vy averaged*/
+	this->inputs->DeleteInput(VxAverageEnum);
+	this->inputs->DeleteInput(VyAverageEnum);
+
+	/*clean up and return*/
+	return mass_flux;
 }
 /*}}}*/
@@ -1836,4 +1797,103 @@
 
 	return minlength;
+}
+/*}}}*/
+Gauss*     Penta::NewGauss(void){/*{{{*/
+	return new GaussPenta();
+}
+/*}}}*/
+Gauss*     Penta::NewGauss(int order){/*{{{*/
+	return new GaussPenta(order,order);
+}
+/*}}}*/
+Gauss*     Penta::NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert){/*{{{*/
+
+	IssmDouble  area_coordinates[4][3];
+
+	GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,4);
+
+	return new GaussPenta(area_coordinates,order_horiz,order_vert);
+}
+/*}}}*/
+Gauss*     Penta::NewGaussBase(int order){/*{{{*/
+	return new GaussPenta(0,1,2,order);
+}
+/*}}}*/
+Gauss*     Penta::NewGaussLine(int vertex1,int vertex2,int order){/*{{{*/
+	return new GaussPenta(vertex1,vertex2,order);
+}
+/*}}}*/
+Gauss*     Penta::NewGaussTop(int order){/*{{{*/
+	return new GaussPenta(3,4,5,order);
+}
+/*}}}*/
+void       Penta::NodalFunctions(IssmDouble* basis, Gauss* gauss){/*{{{*/
+
+	_assert_(gauss->Enum()==GaussPentaEnum);
+	this->GetNodalFunctions(basis,(GaussPenta*)gauss,this->element_type);
+
+}
+/*}}}*/
+void       Penta::NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+
+	_assert_(gauss->Enum()==GaussPentaEnum);
+	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss,this->element_type);
+
+}
+/*}}}*/
+void       Penta::NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+
+	_assert_(gauss->Enum()==GaussPentaEnum);
+	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss,this->VelocityInterpolation());
+
+}
+/*}}}*/
+void       Penta::NodalFunctionsMINIDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+
+	_assert_(gauss->Enum()==GaussPentaEnum);
+	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss,P1bubbleEnum);
+
+}
+/*}}}*/
+void       Penta::NodalFunctionsPressure(IssmDouble* basis, Gauss* gauss){/*{{{*/
+
+	_assert_(gauss->Enum()==GaussPentaEnum);
+	this->GetNodalFunctions(basis,(GaussPenta*)gauss,this->PressureInterpolation());
+
+}
+/*}}}*/
+void       Penta::NodalFunctionsP1(IssmDouble* basis, Gauss* gauss){/*{{{*/
+
+	_assert_(gauss->Enum()==GaussPentaEnum);
+	this->GetNodalFunctions(basis,(GaussPenta*)gauss,P1Enum);
+
+}
+/*}}}*/
+void       Penta::NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+
+	_assert_(gauss->Enum()==GaussPentaEnum);
+	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss,P1Enum);
+
+}
+/*}}}*/
+void       Penta::NodalFunctionsP2(IssmDouble* basis, Gauss* gauss){/*{{{*/
+
+	_assert_(gauss->Enum()==GaussPentaEnum);
+	this->GetNodalFunctions(basis,(GaussPenta*)gauss,P2Enum);
+
+}
+/*}}}*/
+void       Penta::NodalFunctionsVelocity(IssmDouble* basis, Gauss* gauss){/*{{{*/
+
+	_assert_(gauss->Enum()==GaussPentaEnum);
+	this->GetNodalFunctions(basis,(GaussPenta*)gauss,this->VelocityInterpolation());
+
+}
+/*}}}*/
+void       Penta::NodalFunctionsTensor(IssmDouble* basis, Gauss* gauss){/*{{{*/
+
+	_assert_(gauss->Enum()==GaussPentaEnum);
+	this->GetNodalFunctions(basis,(GaussPenta*)gauss,this->TensorInterpolation());
+
 }
 /*}}}*/
@@ -1869,103 +1929,4 @@
 }
 /*}}}*/
-Gauss*     Penta::NewGauss(void){/*{{{*/
-	return new GaussPenta();
-}
-/*}}}*/
-Gauss*     Penta::NewGauss(int order){/*{{{*/
-	return new GaussPenta(order,order);
-}
-/*}}}*/
-Gauss*     Penta::NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert){/*{{{*/
-
-	IssmDouble  area_coordinates[4][3];
-
-	GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,4);
-
-	return new GaussPenta(area_coordinates,order_horiz,order_vert);
-}
-/*}}}*/
-Gauss*     Penta::NewGaussBase(int order){/*{{{*/
-	return new GaussPenta(0,1,2,order);
-}
-/*}}}*/
-Gauss*     Penta::NewGaussLine(int vertex1,int vertex2,int order){/*{{{*/
-	return new GaussPenta(vertex1,vertex2,order);
-}
-/*}}}*/
-Gauss*     Penta::NewGaussTop(int order){/*{{{*/
-	return new GaussPenta(3,4,5,order);
-}
-/*}}}*/
-void       Penta::NodalFunctions(IssmDouble* basis, Gauss* gauss){/*{{{*/
-
-	_assert_(gauss->Enum()==GaussPentaEnum);
-	this->GetNodalFunctions(basis,(GaussPenta*)gauss,this->element_type);
-
-}
-/*}}}*/
-void       Penta::NodalFunctionsP1(IssmDouble* basis, Gauss* gauss){/*{{{*/
-
-	_assert_(gauss->Enum()==GaussPentaEnum);
-	this->GetNodalFunctions(basis,(GaussPenta*)gauss,P1Enum);
-
-}
-/*}}}*/
-void       Penta::NodalFunctionsP2(IssmDouble* basis, Gauss* gauss){/*{{{*/
-
-	_assert_(gauss->Enum()==GaussPentaEnum);
-	this->GetNodalFunctions(basis,(GaussPenta*)gauss,P2Enum);
-
-}
-/*}}}*/
-void       Penta::NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-
-	_assert_(gauss->Enum()==GaussPentaEnum);
-	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss,this->element_type);
-
-}
-/*}}}*/
-void       Penta::NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-
-	_assert_(gauss->Enum()==GaussPentaEnum);
-	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss,P1Enum);
-
-}
-/*}}}*/
-void       Penta::NodalFunctionsMINIDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-
-	_assert_(gauss->Enum()==GaussPentaEnum);
-	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss,P1bubbleEnum);
-
-}
-/*}}}*/
-void       Penta::NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-
-	_assert_(gauss->Enum()==GaussPentaEnum);
-	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss,this->VelocityInterpolation());
-
-}
-/*}}}*/
-void       Penta::NodalFunctionsVelocity(IssmDouble* basis, Gauss* gauss){/*{{{*/
-
-	_assert_(gauss->Enum()==GaussPentaEnum);
-	this->GetNodalFunctions(basis,(GaussPenta*)gauss,this->VelocityInterpolation());
-
-}
-/*}}}*/
-void       Penta::NodalFunctionsPressure(IssmDouble* basis, Gauss* gauss){/*{{{*/
-
-	_assert_(gauss->Enum()==GaussPentaEnum);
-	this->GetNodalFunctions(basis,(GaussPenta*)gauss,this->PressureInterpolation());
-
-}
-/*}}}*/
-void       Penta::NodalFunctionsTensor(IssmDouble* basis, Gauss* gauss){/*{{{*/
-
-	_assert_(gauss->Enum()==GaussPentaEnum);
-	this->GetNodalFunctions(basis,(GaussPenta*)gauss,this->TensorInterpolation());
-
-}
-/*}}}*/
 void       Penta::NormalBase(IssmDouble* bed_normal,IssmDouble* xyz_list){/*{{{*/
 
@@ -1990,4 +1951,24 @@
 }
 /*}}}*/
+void       Penta::NormalSection(IssmDouble* normal,IssmDouble* xyz_list){/*{{{*/
+
+	/*Build unit outward pointing vector*/
+	IssmDouble AB[3];
+	IssmDouble AC[3];
+	IssmDouble norm;
+
+	AB[0]=xyz_list[1*3+0] - xyz_list[0*3+0];
+	AB[1]=xyz_list[1*3+1] - xyz_list[0*3+1];
+	AB[2]=xyz_list[1*3+2] - xyz_list[0*3+2];
+	AC[0]=xyz_list[2*3+0] - xyz_list[0*3+0];
+	AC[1]=xyz_list[2*3+1] - xyz_list[0*3+1];
+	AC[2]=xyz_list[2*3+2] - xyz_list[0*3+2];
+
+	cross(normal,AB,AC);
+	norm=sqrt(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]);
+
+	for(int i=0;i<3;i++) normal[i]=normal[i]/norm;
+}
+/*}}}*/
 void       Penta::NormalTop(IssmDouble* top_normal,IssmDouble* xyz_list){/*{{{*/
 
@@ -2018,4 +1999,10 @@
 int        Penta::NumberofNodesVelocity(void){/*{{{*/
 	return PentaRef::NumberofNodes(this->VelocityInterpolation());
+}
+/*}}}*/
+int        Penta::ObjectEnum(void){/*{{{*/
+
+	return PentaEnum;
+
 }
 /*}}}*/
@@ -2089,4 +2076,35 @@
 }
 /*}}}*/
+void       Penta::PotentialUngrounding(Vector<IssmDouble>* potential_ungrounding){/*{{{*/
+
+	IssmDouble  h[NUMVERTICES],r[NUMVERTICES],gl[NUMVERTICES];
+	IssmDouble  bed_hydro;
+	IssmDouble  rho_water,rho_ice,density;
+
+	/*material parameters: */
+	rho_water=matpar->GetRhoWater();
+	rho_ice=matpar->GetRhoIce();
+	density=rho_ice/rho_water;
+	GetInputListOnVertices(&h[0],ThicknessEnum);
+	GetInputListOnVertices(&r[0],BedEnum);
+	GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum);
+
+	/*go through vertices, and figure out which ones are on the ice sheet, and want to unground: */
+	for(int i=0;i<NUMVERTICES;i++){
+		/*Find if grounded vertices want to start floating*/
+		if (gl[i]>0.){
+			bed_hydro=-density*h[i];
+			if(bed_hydro>r[i]){
+				/*Vertex that could potentially unground, flag it*/
+				potential_ungrounding->SetValue(vertices[i]->Pid(),1,INS_VAL);
+			}
+		}
+	}
+}
+/*}}}*/
+int        Penta::PressureInterpolation(void){/*{{{*/
+	return PentaRef::PressureInterpolation(this->element_type);
+}
+/*}}}*/
 void       Penta::ReduceMatrices(ElementMatrix* Ke,ElementVector* pe){/*{{{*/
 
@@ -2207,24 +2225,4 @@
 }
 /*}}}*/
-void       Penta::SetClone(int* minranks){/*{{{*/
-
-	_error_("not implemented yet");
-}
-/*}}}*/
-void       Penta::SetCurrentConfiguration(Elements* elementsin, Loads* loadsin, Nodes* nodesin, Materials* materialsin, Parameters* parametersin){/*{{{*/
-
-	/*go into parameters and get the analysis_counter: */
-	int analysis_counter;
-	parametersin->FindParam(&analysis_counter,AnalysisCounterEnum);
-
-	/*Get Element type*/
-	this->element_type=this->element_type_list[analysis_counter];
-
-	/*Pick up nodes*/
-	if(this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp();
-	else this->nodes=NULL;
-
-}
-/*}}}*/
 void       Penta::ResetHooks(){/*{{{*/
 
@@ -2245,33 +2243,68 @@
 }
 /*}}}*/
+void       Penta::SetClone(int* minranks){/*{{{*/
+
+	_error_("not implemented yet");
+}
+/*}}}*/
+void       Penta::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){/*{{{*/
+
+	IssmDouble  values[NUMVERTICES];
+	int         vertexpidlist[NUMVERTICES],control_init;
+
+	/*Specific case for depth averaged quantities*/
+	control_init=control_enum;
+	if(control_enum==MaterialsRheologyBbarEnum){
+		control_enum=MaterialsRheologyBEnum;
+		if(!IsOnBase()) return;
+	}
+	if(control_enum==DamageDbarEnum){
+		control_enum=DamageDEnum;
+		if(!IsOnBase()) return;
+	}
+
+	/*Get out if this is not an element input*/
+	if(!IsInput(control_enum)) return;
+
+	/*Prepare index list*/
+	GradientIndexing(&vertexpidlist[0],control_index);
+
+	/*Get values on vertices*/
+	for(int i=0;i<NUMVERTICES;i++){
+		values[i]=vector[vertexpidlist[i]];
+	}
+	Input* new_input = new PentaInput(control_enum,values,P1Enum);
+	Input* input=(Input*)this->inputs->GetInput(control_enum);   _assert_(input);
+	if(input->ObjectEnum()!=ControlInputEnum){
+		_error_("input " << EnumToStringx(control_enum) << " is not a ControlInput");
+	}
+
+	((ControlInput*)input)->SetInput(new_input);
+
+	if(control_init==MaterialsRheologyBbarEnum){
+		this->InputExtrude(control_enum,-1);
+	}
+	if(control_init==DamageDbarEnum){
+		this->InputExtrude(control_enum,-1);
+	}
+}
+/*}}}*/
+void       Penta::SetCurrentConfiguration(Elements* elementsin, Loads* loadsin, Nodes* nodesin, Materials* materialsin, Parameters* parametersin){/*{{{*/
+
+	/*go into parameters and get the analysis_counter: */
+	int analysis_counter;
+	parametersin->FindParam(&analysis_counter,AnalysisCounterEnum);
+
+	/*Get Element type*/
+	this->element_type=this->element_type_list[analysis_counter];
+
+	/*Pick up nodes*/
+	if(this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp();
+	else this->nodes=NULL;
+
+}
+/*}}}*/
 void       Penta::SetTemporaryElementType(int element_type_in){/*{{{*/
 	this->element_type=element_type_in;
-}
-/*}}}*/
-Tria*      Penta::SpawnTria(int index1,int index2,int index3){/*{{{*/
-
-	int analysis_counter;
-
-	/*go into parameters and get the analysis_counter: */
-	this->parameters->FindParam(&analysis_counter,AnalysisCounterEnum);
-
-	/*Create Tria*/
-	Tria* tria=new Tria();
-	tria->id=this->id;
-	tria->inputs=(Inputs*)this->inputs->SpawnTriaInputs(index1,index2,index3);
-	tria->parameters=this->parameters;
-	tria->element_type=P1Enum; //Only P1 CG for now (TO BE CHANGED)
-	this->SpawnTriaHook(xDynamicCast<ElementHook*>(tria),index1,index2,index3);
-
-	/*Spawn material*/
-	tria->material=(Material*)this->material->copy2(tria);
-
-	/*recover nodes, material and matpar: */
-	tria->nodes=(Node**)tria->hnodes[analysis_counter]->deliverp();
-	tria->vertices=(Vertex**)tria->hvertices->deliverp();
-	tria->matpar=(Matpar*)tria->hmatpar->delivers();
-
-	/*Return new Tria*/
-	return tria;
 }
 /*}}}*/
@@ -2306,4 +2339,219 @@
 }
 /*}}}*/
+Tria*      Penta::SpawnTria(int index1,int index2,int index3){/*{{{*/
+
+	int analysis_counter;
+
+	/*go into parameters and get the analysis_counter: */
+	this->parameters->FindParam(&analysis_counter,AnalysisCounterEnum);
+
+	/*Create Tria*/
+	Tria* tria=new Tria();
+	tria->id=this->id;
+	tria->inputs=(Inputs*)this->inputs->SpawnTriaInputs(index1,index2,index3);
+	tria->parameters=this->parameters;
+	tria->element_type=P1Enum; //Only P1 CG for now (TO BE CHANGED)
+	this->SpawnTriaHook(xDynamicCast<ElementHook*>(tria),index1,index2,index3);
+
+	/*Spawn material*/
+	tria->material=(Material*)this->material->copy2(tria);
+
+	/*recover nodes, material and matpar: */
+	tria->nodes=(Node**)tria->hnodes[analysis_counter]->deliverp();
+	tria->vertices=(Vertex**)tria->hvertices->deliverp();
+	tria->matpar=(Matpar*)tria->hmatpar->delivers();
+
+	/*Return new Tria*/
+	return tria;
+}
+/*}}}*/
+IssmDouble Penta::StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){/*{{{*/
+	/*Compute stabilization parameter*/
+	/*kappa=thermalconductivity/(rho_ice*hearcapacity) for thermal model*/
+	/*kappa=enthalpydiffusionparameter for enthalpy model*/
+
+	IssmDouble normu;
+	IssmDouble tau_parameter;
+
+	normu=pow(pow(u,2)+pow(v,2)+pow(w,2),0.5);
+	if(normu*diameter/(3*2*kappa)<1){ 
+		tau_parameter=pow(diameter,2)/(3*2*2*kappa);
+	}
+	else tau_parameter=diameter/(2*normu);
+
+	return tau_parameter;
+}
+/*}}}*/
+void       Penta::StrainRateparallel(){/*{{{*/
+
+	IssmDouble *xyz_list = NULL;
+	IssmDouble  epsilon[6];
+	GaussPenta* gauss=NULL;
+	IssmDouble  vx,vy,vel;
+	IssmDouble  strainxx;
+	IssmDouble  strainxy;
+	IssmDouble  strainyy;
+	IssmDouble  strainparallel[NUMVERTICES];
+
+	/* Get node coordinates and dof list: */
+	this->GetVerticesCoordinates(&xyz_list);
+
+	/*Retrieve all inputs we will need*/
+	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);
+
+	/* Start looping on the number of vertices: */
+	gauss=new GaussPenta();
+	for (int iv=0;iv<NUMVERTICES;iv++){
+		gauss->GaussVertex(iv);
+
+		/* Get the value we need*/
+		vx_input->GetInputValue(&vx,gauss);
+		vy_input->GetInputValue(&vy,gauss);
+		vel=vx*vx+vy*vy;
+
+		/*Compute strain rate and viscosity: */
+		this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
+		strainxx=epsilon[0];
+		strainyy=epsilon[1];
+		strainxy=epsilon[3];
+
+		/*strainparallel= Strain rate along the ice flow direction */
+		strainparallel[iv]=(vx*vx*(strainxx)+vy*vy*(strainyy)+2*vy*vx*strainxy)/(vel+1.e-6);
+	}
+
+	/*Add input*/
+	this->inputs->AddInput(new PentaInput(StrainRateparallelEnum,&strainparallel[0],P1Enum));
+
+	/*Clean up and return*/
+	delete gauss;
+	xDelete<IssmDouble>(xyz_list);
+}
+/*}}}*/
+void       Penta::StrainRateperpendicular(){/*{{{*/
+
+	IssmDouble *xyz_list = NULL;
+	IssmDouble  epsilon[6];
+	GaussPenta* gauss=NULL;
+	IssmDouble  vx,vy,vel;
+	IssmDouble  strainxx;
+	IssmDouble  strainxy;
+	IssmDouble  strainyy;
+	IssmDouble  strainperpendicular[NUMVERTICES];
+
+	/* Get node coordinates and dof list: */
+	this->GetVerticesCoordinates(&xyz_list);
+
+	/*Retrieve all inputs we will need*/
+	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);
+
+	/* Start looping on the number of vertices: */
+	gauss=new GaussPenta();
+	for (int iv=0;iv<NUMVERTICES;iv++){
+		gauss->GaussVertex(iv);
+
+		/* Get the value we need*/
+		vx_input->GetInputValue(&vx,gauss);
+		vy_input->GetInputValue(&vy,gauss);
+		vel=vx*vx+vy*vy;
+
+		/*Compute strain rate and viscosity: */
+		this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
+		strainxx=epsilon[0];
+		strainyy=epsilon[1];
+		strainxy=epsilon[3];
+
+		/*strainperpendicular= Strain rate perpendicular to the ice flow direction */
+		strainperpendicular[iv]=(vx*vx*(strainyy)+vy*vy*(strainxx)-2*vy*vx*strainxy)/(vel+1.e-6);
+	}
+
+	/*Add input*/
+	this->inputs->AddInput(new PentaInput(StrainRateperpendicularEnum,&strainperpendicular[0],P1Enum));
+
+	/*Clean up and return*/
+	delete gauss;
+	xDelete<IssmDouble>(xyz_list);
+}
+/*}}}*/
+void       Penta::StressIntensityFactor(){/*{{{*/
+
+	/* Check if we are on the base */
+	if(!IsOnBase()) return;
+
+	IssmDouble  ki[6]={0.};
+	IssmDouble  const_grav=9.81;
+	IssmDouble  rho_ice=900;
+	IssmDouble  rho_water=1000;
+	IssmDouble  Jdet[3];
+	IssmDouble  pressure,vx,vy,vel,deviaxx,deviaxy,deviayy,water_depth,prof,stress_xx,thickness;
+
+	Penta* penta=this;
+	for(;;){
+	
+		IssmDouble  xyz_list[NUMVERTICES][3];
+		/* Get node coordinates and dof list: */
+		::GetVerticesCoordinates(&xyz_list[0][0],penta->vertices,NUMVERTICES);
+
+		///*Compute the Jacobian for the vertical integration*/
+		Jdet[0]=(xyz_list[3][2]-xyz_list[0][2])*0.5;
+		Jdet[1]=(xyz_list[4][2]-xyz_list[1][2])*0.5;
+		Jdet[2]=(xyz_list[5][2]-xyz_list[2][2])*0.5;
+	
+		/*Retrieve all inputs we will need*/
+		Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
+		Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
+		Input* vel_input=inputs->GetInput(VelEnum);                                _assert_(vel_input);
+		Input* pressure_input=inputs->GetInput(PressureEnum);                      _assert_(pressure_input);
+		Input* deviaxx_input=inputs->GetInput(DeviatoricStressxxEnum);             _assert_(deviaxx_input);
+		Input* deviaxy_input=inputs->GetInput(DeviatoricStressxyEnum);             _assert_(deviaxy_input);
+		Input* deviayy_input=inputs->GetInput(DeviatoricStressyyEnum);             _assert_(deviayy_input);
+		Input* surface_input=inputs->GetInput(SurfaceEnum);								_assert_(surface_input);
+		Input* thickness_input=inputs->GetInput(ThicknessEnum);							_assert_(thickness_input);
+		
+		/* Start looping on the number of 2D vertices: */
+		for(int ig=0;ig<3;ig++){
+			GaussPenta* gauss=new GaussPenta(ig,3+ig,11);
+			for (int iv=gauss->begin();iv<gauss->end();iv++){
+				gauss->GaussPoint(iv);
+
+				/* Get the value we need*/
+				pressure_input->GetInputValue(&pressure,gauss);
+				vx_input->GetInputValue(&vx,gauss);
+				vy_input->GetInputValue(&vy,gauss);
+				vel_input->GetInputValue(&vel,gauss);
+				deviaxx_input->GetInputValue(&deviaxx,gauss);
+				deviaxy_input->GetInputValue(&deviaxy,gauss);
+				deviayy_input->GetInputValue(&deviayy,gauss);
+				surface_input->GetInputValue(&water_depth,gauss);
+				thickness_input->GetInputValue(&thickness,gauss);
+				prof=water_depth-penta->GetZcoord(&xyz_list[0][0],gauss);
+
+				/*stress_xx= Deviatoric stress along the ice flow direction plus cryostatic pressure */
+				stress_xx=(vx*vx*(deviaxx)+vy*vy*(deviayy)+2*vy*vx*deviaxy)/(vel*vel+1.e-6);
+
+				if(prof<water_depth&prof<thickness){
+					/* Compute the local stress intensity factor*/ 
+					ki[ig]+=Jdet[ig]*gauss->weight*stress_xx*StressIntensityIntegralWeight(prof,min(water_depth,thickness),thickness);
+				}
+			}
+			delete gauss;
+		}
+			
+		/*Stop if we have reached the surface/base*/
+		if(penta->IsOnSurface()) break;
+		
+		/*get upper Penta*/
+		penta=penta->GetUpperPenta();
+		_assert_(penta->Id()!=this->id);
+	}
+
+	/*Add input*/
+	this->inputs->AddInput(new PentaInput(StressIntensityFactorEnum,&ki[0],P1Enum));
+	this->InputExtrude(StressIntensityFactorEnum,-1);
+}
+/*}}}*/
 IssmDouble Penta::SurfaceArea(void){/*{{{*/
 
@@ -2385,4 +2633,33 @@
 	return dt;
 }/*}}}*/
+IssmDouble Penta::TotalSmb(void){/*{{{*/
+
+	/*The smb[Gt yr-1] of one element is area[m2] * smb [ m ice yr^-1] * rho_ice [kg m-3] / 1e+10^12 */
+	IssmDouble base,smb,rho_ice;
+	IssmDouble Total_Smb=0;
+	IssmDouble xyz_list[NUMVERTICES][3];
+
+	/*Get material parameters :*/
+	rho_ice=matpar->GetRhoIce();
+
+	if(!IsIceInElement() || !IsOnSurface()) return 0.;
+
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+
+	/*First calculate the area of the base (cross section triangle)
+	 * http://en.wikipedia.org/wiki/Triangle
+	 * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
+	base = 1./2. * fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
+
+	/*Now get the average SMB over the element*/
+	Input* smb_input = inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(smb_input);
+
+	smb_input->GetInputAverage(&smb);
+	Total_Smb=rho_ice*base*smb;// smb on element in kg s-1
+
+	/*Return: */
+	return Total_Smb;
+}
+/*}}}*/
 void       Penta::Update(int index,IoModel* iomodel,int analysis_counter,int analysis_type,int finiteelement_type){ /*{{{*/
 
@@ -2833,10 +3110,33 @@
 }
 /*}}}*/
+int        Penta::UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf){/*{{{*/
+
+	int i;
+	int nflipped=0;
+
+	/*Go through nodes, and whoever is on the potential_ungrounding, ends up in nodes_on_iceshelf: */
+	for(i=0;i<NUMVERTICES;i++){
+		if (reCast<bool,IssmDouble>(vertices_potentially_ungrounding[vertices[i]->Pid()])){
+			vec_nodes_on_iceshelf->SetValue(vertices[i]->Pid(),-1.,INS_VAL);
+
+			/*If node was not on ice shelf, we flipped*/
+			if(nodes_on_iceshelf[vertices[i]->Pid()]>=0.){
+				nflipped++;
+			}
+		}
+	}
+	return nflipped;
+}
+/*}}}*/
+void       Penta::ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+	PentaRef::GetInputDerivativeValue(dvalue,values,xyz_list,gauss,P1Enum);
+}
+/*}}}*/
 void       Penta::ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss){/*{{{*/
 	PentaRef::GetInputValue(pvalue,values,gauss,P1Enum);
 }
 /*}}}*/
-void       Penta::ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	PentaRef::GetInputDerivativeValue(dvalue,values,xyz_list,gauss,P1Enum);
+int        Penta::VelocityInterpolation(void){/*{{{*/
+	return PentaRef::VelocityInterpolation(this->element_type);
 }
 /*}}}*/
@@ -2875,207 +3175,157 @@
 }
 /*}}}*/
-int        Penta::VelocityInterpolation(void){/*{{{*/
-	return PentaRef::VelocityInterpolation(this->element_type);
-}
-/*}}}*/
-int        Penta::PressureInterpolation(void){/*{{{*/
-	return PentaRef::PressureInterpolation(this->element_type);
-}
-/*}}}*/
-bool       Penta::IsZeroLevelset(int levelset_enum){/*{{{*/
-
-	bool        iszerols;
-	IssmDouble  ls[NUMVERTICES];
-
-	/*Retrieve all inputs and parameters*/
-	GetInputListOnVertices(&ls[0],levelset_enum);
-
-	/*If the level set has always same sign, there is no ice front here*/
-	iszerols = false;
-	if(IsIceInElement()){
-		if(ls[0]*ls[1]<0. || ls[0]*ls[2]<0. || (ls[0]*ls[1]*ls[2]==0. && ls[0]*ls[1]+ls[0]*ls[2]+ls[1]*ls[2]<=0.)){
-			iszerols = true;
-		}
-	}
-	return iszerols;
-}
-/*}}}*/
-bool       Penta::IsIcefront(void){/*{{{*/
-
-	bool isicefront;
-	int i,nrice;
-   IssmDouble ls[NUMVERTICES];
-
-	/*Retrieve all inputs and parameters*/
-	GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
-
-	/* If only one vertex has ice, there is an ice front here */
-	isicefront=false;
-	if(IsIceInElement()){
-		nrice=0;       
-		for(i=0;i<NUMVERTICES2D;i++)
-			if(ls[i]<0.) nrice++;
-		if(nrice==1) isicefront= true;
-	}
-	return isicefront;
-}/*}}}*/
-void       Penta::AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part){/*{{{*/
-	_error_("Not supported yet!");
-}
-/*}}}*/
-IssmDouble Penta::IceVolume(void){/*{{{*/
-
-	/*The volume of a troncated prism is base * 1/3 sum(length of edges)*/
-	IssmDouble base,height;
-	IssmDouble xyz_list[NUMVERTICES][3];
-
-	if(!IsIceInElement())return 0;
-
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-
-	/*First calculate the area of the base (cross section triangle)
-	 * http://en.wikipedia.org/wiki/Pentangle
-	 * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
-	base = 1./2.*fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
-
-	/*Now get the average height*/
-	height = 1./3.*((xyz_list[3][2]-xyz_list[0][2])+(xyz_list[4][2]-xyz_list[1][2])+(xyz_list[5][2]-xyz_list[2][2]));
-
-	/*Return: */
-	return base*height;
-}
-/*}}}*/
-IssmDouble Penta::IceVolumeAboveFloatation(void){/*{{{*/
-
-	/*Volume above floatation: H + rho_water/rho_ice*bathymetry for nodes on the bed*/
-	IssmDouble rho_ice,rho_water;
-	IssmDouble base,bed,surface,bathymetry;
-	IssmDouble xyz_list[NUMVERTICES][3];
-
-	if(!IsIceInElement() || IsFloating() || !IsOnBase())return 0;
-
-	rho_ice=matpar->GetRhoIce();
-	rho_water=matpar->GetRhoWater();
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-
-	/*First calculate the area of the base (cross section triangle)
-	 * http://en.wikipedia.org/wiki/Pentangle
-	 * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
-	base = 1./2.*fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
-
-	/*Now get the average height above floatation*/
-	Input* surface_input    = inputs->GetInput(SurfaceEnum);    _assert_(surface_input);
-	Input* base_input        = inputs->GetInput(BaseEnum);        _assert_(base_input);
-	Input* bed_input = inputs->GetInput(BedEnum); _assert_(bed_input);
-	surface_input->GetInputAverage(&surface);
-	base_input->GetInputAverage(&bed);
-	bed_input->GetInputAverage(&bathymetry);
-
-	/*Return: */
-	return base*(surface - bed + min( rho_water/rho_ice * bathymetry, 0.) );
-}
-/*}}}*/
-IssmDouble Penta::MassFlux( IssmDouble* segment){/*{{{*/
-
-	IssmDouble mass_flux=0;
-
-	if(!IsOnBase()) return mass_flux;
-
-	/*Depth Averaging Vx and Vy*/
-	this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
-	this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
-
-	/*Spawn Tria element from the base of the Penta: */
-	Tria* tria=(Tria*)SpawnTria(0,1,2);
-	mass_flux=tria->MassFlux(segment);
-	delete tria->material; delete tria;
-
-	/*Delete Vx and Vy averaged*/
-	this->inputs->DeleteInput(VxAverageEnum);
-	this->inputs->DeleteInput(VyAverageEnum);
-
-	/*clean up and return*/
-	return mass_flux;
-}
-/*}}}*/
-IssmDouble Penta::MassFlux( IssmDouble x1, IssmDouble y1, IssmDouble x2, IssmDouble y2,int segment_id){/*{{{*/
-
-	IssmDouble mass_flux=0;
-
-	if(!IsOnBase()) return mass_flux;
-
-	/*Depth Averaging Vx and Vy*/
-	this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
-	this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
-
-	/*Spawn Tria element from the base of the Penta: */
-	Tria* tria=(Tria*)SpawnTria(0,1,2);
-	mass_flux=tria->MassFlux(x1,y1,x2,y2,segment_id);
-	delete tria->material; delete tria;
-
-	/*Delete Vx and Vy averaged*/
-	this->inputs->DeleteInput(VxAverageEnum);
-	this->inputs->DeleteInput(VyAverageEnum);
-
-	/*clean up and return*/
-	return mass_flux;
-}
-/*}}}*/
-void       Penta::ElementResponse(IssmDouble* presponse,int response_enum){/*{{{*/
-
-	switch(response_enum){
-		case MaterialsRheologyBbarEnum:
-			*presponse=this->material->GetBbar();
-			break;
-		case DamageDbarEnum:
-			*presponse=this->material->GetDbar();
-			break;
-		case VelEnum:
-			{
-
-				/*Get input:*/
-				IssmDouble vel;
-				Input* vel_input;
-
-				vel_input=this->inputs->GetInput(VelEnum); _assert_(vel_input);
-				vel_input->GetInputAverage(&vel);
-
-				/*Assign output pointers:*/
-				*presponse=vel;
-			}
-			break;
-		default:  
-			_error_("Response type " << EnumToStringx(response_enum) << " not supported yet!");
-	}
-
-}
-/*}}}*/
-IssmDouble Penta::TotalSmb(void){/*{{{*/
-
-	/*The smb[Gt yr-1] of one element is area[m2] * smb [ m ice yr^-1] * rho_ice [kg m-3] / 1e+10^12 */
-	IssmDouble base,smb,rho_ice;
-	IssmDouble Total_Smb=0;
-	IssmDouble xyz_list[NUMVERTICES][3];
-
-	/*Get material parameters :*/
-	rho_ice=matpar->GetRhoIce();
-
-	if(!IsIceInElement() || !IsOnSurface()) return 0.;
-
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-
-	/*First calculate the area of the base (cross section triangle)
-	 * http://en.wikipedia.org/wiki/Triangle
-	 * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
-	base = 1./2. * fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));
-
-	/*Now get the average SMB over the element*/
-	Input* smb_input = inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(smb_input);
-
-	smb_input->GetInputAverage(&smb);
-	Total_Smb=rho_ice*base*smb;// smb on element in kg s-1
-
-	/*Return: */
-	return Total_Smb;
+void       Penta::ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
+	/*Compute portion of the element that is grounded*/ 
+
+	int         normal_orientation=0;
+	IssmDouble  s1,s2;
+	IssmDouble  levelset[NUMVERTICES];
+	IssmDouble* xyz_zero = xNew<IssmDouble>(4*3);
+
+	/*Recover parameters and values*/
+	GetInputListOnVertices(&levelset[0],levelsetenum);
+
+	if(levelset[0]*levelset[1]>0.){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
+		/*Portion of the segments*/
+		s1=levelset[2]/(levelset[2]-levelset[1]);
+		s2=levelset[2]/(levelset[2]-levelset[0]);
+
+		if(levelset[2]<0.) normal_orientation=1; //orientation of quadrangle at base and top, depending on distribution of levelsetfunction
+		/*New point 1*/
+		xyz_zero[3*normal_orientation+0]=xyz_list[2*3+0]+s1*(xyz_list[1*3+0]-xyz_list[2*3+0]);
+		xyz_zero[3*normal_orientation+1]=xyz_list[2*3+1]+s1*(xyz_list[1*3+1]-xyz_list[2*3+1]);
+		xyz_zero[3*normal_orientation+2]=xyz_list[2*3+2]+s1*(xyz_list[1*3+2]-xyz_list[2*3+2]);
+
+		/*New point 0*/
+		xyz_zero[3*(1-normal_orientation)+0]=xyz_list[2*3+0]+s2*(xyz_list[0*3+0]-xyz_list[2*3+0]);
+		xyz_zero[3*(1-normal_orientation)+1]=xyz_list[2*3+1]+s2*(xyz_list[0*3+1]-xyz_list[2*3+1]);
+		xyz_zero[3*(1-normal_orientation)+2]=xyz_list[2*3+2]+s2*(xyz_list[0*3+2]-xyz_list[2*3+2]);
+
+		/*New point 3*/
+		xyz_zero[3*(2+1-normal_orientation)+0]=xyz_list[5*3+0]+s1*(xyz_list[4*3+0]-xyz_list[5*3+0]);
+		xyz_zero[3*(2+1-normal_orientation)+1]=xyz_list[5*3+1]+s1*(xyz_list[4*3+1]-xyz_list[5*3+1]);
+		xyz_zero[3*(2+1-normal_orientation)+2]=xyz_list[5*3+2]+s1*(xyz_list[4*3+2]-xyz_list[5*3+2]);
+
+		/*New point 4*/
+		xyz_zero[3*(2+normal_orientation)+0]=xyz_list[5*3+0]+s2*(xyz_list[3*3+0]-xyz_list[5*3+0]);
+		xyz_zero[3*(2+normal_orientation)+1]=xyz_list[5*3+1]+s2*(xyz_list[3*3+1]-xyz_list[5*3+1]);
+		xyz_zero[3*(2+normal_orientation)+2]=xyz_list[5*3+2]+s2*(xyz_list[3*3+2]-xyz_list[5*3+2]);
+	}
+	else if(levelset[1]*levelset[2]>0.){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
+		/*Portion of the segments*/
+		s1=levelset[0]/(levelset[0]-levelset[2]);
+		s2=levelset[0]/(levelset[0]-levelset[1]);
+
+		if(levelset[0]<0.) normal_orientation=1;
+		/*New point 1*/
+		xyz_zero[3*normal_orientation+0]=xyz_list[0*3+0]+s1*(xyz_list[2*3+0]-xyz_list[0*3+0]);
+		xyz_zero[3*normal_orientation+1]=xyz_list[0*3+1]+s1*(xyz_list[2*3+1]-xyz_list[0*3+1]);
+		xyz_zero[3*normal_orientation+2]=xyz_list[0*3+2]+s1*(xyz_list[2*3+2]-xyz_list[0*3+2]);
+
+		/*New point 2*/
+		xyz_zero[3*(1-normal_orientation)+0]=xyz_list[0*3+0]+s2*(xyz_list[1*3+0]-xyz_list[0*3+0]);
+		xyz_zero[3*(1-normal_orientation)+1]=xyz_list[0*3+1]+s2*(xyz_list[1*3+1]-xyz_list[0*3+1]);
+		xyz_zero[3*(1-normal_orientation)+2]=xyz_list[0*3+2]+s2*(xyz_list[1*3+2]-xyz_list[0*3+2]);
+
+		/*New point 3*/
+		xyz_zero[3*(2+1-normal_orientation)+0]=xyz_list[3*3+0]+s1*(xyz_list[5*3+0]-xyz_list[3*3+0]);
+		xyz_zero[3*(2+1-normal_orientation)+1]=xyz_list[3*3+1]+s1*(xyz_list[5*3+1]-xyz_list[3*3+1]);
+		xyz_zero[3*(2+1-normal_orientation)+2]=xyz_list[3*3+2]+s1*(xyz_list[5*3+2]-xyz_list[3*3+2]);
+
+		/*New point 4*/
+		xyz_zero[3*(2+normal_orientation)+0]=xyz_list[3*3+0]+s2*(xyz_list[4*3+0]-xyz_list[3*3+0]);
+		xyz_zero[3*(2+normal_orientation)+1]=xyz_list[3*3+1]+s2*(xyz_list[4*3+1]-xyz_list[3*3+1]);
+		xyz_zero[3*(2+normal_orientation)+2]=xyz_list[3*3+2]+s2*(xyz_list[4*3+2]-xyz_list[3*3+2]);
+	}
+	else if(levelset[0]*levelset[2]>0.){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
+		/*Portion of the segments*/
+		s1=levelset[1]/(levelset[1]-levelset[0]);
+		s2=levelset[1]/(levelset[1]-levelset[2]);
+
+		if(levelset[1]<0.) normal_orientation=1;
+		/*New point 0*/
+		xyz_zero[3*normal_orientation+0]=xyz_list[1*3+0]+s1*(xyz_list[0*3+0]-xyz_list[1*3+0]);
+		xyz_zero[3*normal_orientation+1]=xyz_list[1*3+1]+s1*(xyz_list[0*3+1]-xyz_list[1*3+1]);
+		xyz_zero[3*normal_orientation+2]=xyz_list[1*3+2]+s1*(xyz_list[0*3+2]-xyz_list[1*3+2]);
+
+		/*New point 2*/
+		xyz_zero[3*(1-normal_orientation)+0]=xyz_list[1*3+0]+s2*(xyz_list[2*3+0]-xyz_list[1*3+0]);
+		xyz_zero[3*(1-normal_orientation)+1]=xyz_list[1*3+1]+s2*(xyz_list[2*3+1]-xyz_list[1*3+1]);
+		xyz_zero[3*(1-normal_orientation)+2]=xyz_list[1*3+2]+s2*(xyz_list[2*3+2]-xyz_list[1*3+2]);
+
+		/*New point 3*/
+		xyz_zero[3*(2+1-normal_orientation)+0]=xyz_list[4*3+0]+s1*(xyz_list[3*3+0]-xyz_list[4*3+0]);
+		xyz_zero[3*(2+1-normal_orientation)+1]=xyz_list[4*3+1]+s1*(xyz_list[3*3+1]-xyz_list[4*3+1]);
+		xyz_zero[3*(2+1-normal_orientation)+2]=xyz_list[4*3+2]+s1*(xyz_list[3*3+2]-xyz_list[4*3+2]);
+
+		/*New point 4*/
+		xyz_zero[3*(2+normal_orientation)+0]=xyz_list[4*3+0]+s2*(xyz_list[5*3+0]-xyz_list[4*3+0]);
+		xyz_zero[3*(2+normal_orientation)+1]=xyz_list[4*3+1]+s2*(xyz_list[5*3+1]-xyz_list[4*3+1]);
+		xyz_zero[3*(2+normal_orientation)+2]=xyz_list[4*3+2]+s2*(xyz_list[5*3+2]-xyz_list[4*3+2]);
+	}
+	else if(levelset[0]==0. && levelset[1]==0.){ //front is on point 0 and 1
+		xyz_zero[3*0+0]=xyz_list[0*3+0];
+		xyz_zero[3*0+1]=xyz_list[0*3+1];
+		xyz_zero[3*0+2]=xyz_list[0*3+2];
+
+		/*New point 2*/
+		xyz_zero[3*1+0]=xyz_list[1*3+0];
+		xyz_zero[3*1+1]=xyz_list[1*3+1];
+		xyz_zero[3*1+2]=xyz_list[1*3+2];
+
+		/*New point 3*/
+		xyz_zero[3*2+0]=xyz_list[4*3+0];
+		xyz_zero[3*2+1]=xyz_list[4*3+1];
+		xyz_zero[3*2+2]=xyz_list[4*3+2];
+
+		/*New point 4*/
+		xyz_zero[3*3+0]=xyz_list[3*3+0];
+		xyz_zero[3*3+1]=xyz_list[3*3+1];
+		xyz_zero[3*3+2]=xyz_list[3*3+2];
+	}
+	else if(levelset[0]==0. && levelset[2]==0.){ //front is on point 0 and 1
+		xyz_zero[3*0+0]=xyz_list[2*3+0];
+		xyz_zero[3*0+1]=xyz_list[2*3+1];
+		xyz_zero[3*0+2]=xyz_list[2*3+2];
+
+		/*New point 2*/
+		xyz_zero[3*1+0]=xyz_list[0*3+0];
+		xyz_zero[3*1+1]=xyz_list[0*3+1];
+		xyz_zero[3*1+2]=xyz_list[0*3+2];
+
+		/*New point 3*/
+		xyz_zero[3*2+0]=xyz_list[3*3+0];
+		xyz_zero[3*2+1]=xyz_list[3*3+1];
+		xyz_zero[3*2+2]=xyz_list[3*3+2];
+
+		/*New point 4*/
+		xyz_zero[3*3+0]=xyz_list[5*3+0];
+		xyz_zero[3*3+1]=xyz_list[5*3+1];
+		xyz_zero[3*3+2]=xyz_list[5*3+2];
+	}
+	else if(levelset[1]==0. && levelset[2]==0.){ //front is on point 0 and 1
+		xyz_zero[3*0+0]=xyz_list[1*3+0];
+		xyz_zero[3*0+1]=xyz_list[1*3+1];
+		xyz_zero[3*0+2]=xyz_list[1*3+2];
+
+		/*New point 2*/
+		xyz_zero[3*1+0]=xyz_list[2*3+0];
+		xyz_zero[3*1+1]=xyz_list[2*3+1];
+		xyz_zero[3*1+2]=xyz_list[2*3+2];
+
+		/*New point 3*/
+		xyz_zero[3*2+0]=xyz_list[5*3+0];
+		xyz_zero[3*2+1]=xyz_list[5*3+1];
+		xyz_zero[3*2+2]=xyz_list[5*3+2];
+
+		/*New point 4*/
+		xyz_zero[3*3+0]=xyz_list[4*3+0];
+		xyz_zero[3*3+1]=xyz_list[4*3+1];
+		xyz_zero[3*3+2]=xyz_list[4*3+2];
+	}
+	else _error_("Case not covered");
+
+	/*Assign output pointer*/
+	*pxyz_zero= xyz_zero;
 }
 /*}}}*/
@@ -3087,179 +3337,4 @@
 /*}}}*/
 #endif
-
-void       Penta::ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){/*{{{*/
-
-	int    vertexpidlist[NUMVERTICES];
-	IssmDouble grad_list[NUMVERTICES];
-	Input* grad_input=NULL;
-	Input* input=NULL;
-
-	if(enum_type==MaterialsRheologyBbarEnum){
-		input=(Input*)inputs->GetInput(MaterialsRheologyBEnum);
-	}
-	else if(enum_type==DamageDbarEnum){
-		input=(Input*)inputs->GetInput(DamageDEnum);
-	}
-	else{
-		input=inputs->GetInput(enum_type);
-	}
-	if (!input) _error_("Input " << EnumToStringx(enum_type) << " not found");
-	if (input->ObjectEnum()!=ControlInputEnum) _error_("Input " << EnumToStringx(enum_type) << " is not a ControlInput");
-
-	GradientIndexing(&vertexpidlist[0],control_index);
-	for(int i=0;i<NUMVERTICES;i++) grad_list[i]=gradient[vertexpidlist[i]];
-	grad_input=new PentaInput(GradientEnum,grad_list,P1Enum);
-	((ControlInput*)input)->SetGradient(grad_input);
-
-}/*}}}*/
-void       Penta::ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum){/*{{{*/
-
-	Input* input=NULL;
-
-	if(control_enum==MaterialsRheologyBbarEnum){
-		input=(Input*)inputs->GetInput(MaterialsRheologyBEnum);
-	}
-	else if(control_enum==DamageDbarEnum){
-		input=(Input*)inputs->GetInput(DamageDEnum);
-	}
-	else{
-		input=inputs->GetInput(control_enum);
-	}
-	if (!input) _error_("Input " << EnumToStringx(control_enum) << " not found");
-	if (input->ObjectEnum()!=ControlInputEnum) _error_("Input " << EnumToStringx(control_enum) << " is not a ControlInput");
-
-	int         sidlist[NUMVERTICES];
-	int         connectivity[NUMVERTICES];
-	IssmPDouble values[NUMVERTICES];
-	IssmPDouble gradients[NUMVERTICES]; 
-	IssmDouble  value,gradient;
-
-	this->GetVerticesConnectivityList(&connectivity[0]);
-	this->GetVerticesSidList(&sidlist[0]);
-
-	GaussPenta* gauss=new GaussPenta();
-	for (int iv=0;iv<NUMVERTICES;iv++){
-		gauss->GaussVertex(iv);
-
-		((ControlInput*)input)->GetInputValue(&value,gauss);
-		((ControlInput*)input)->GetGradientValue(&gradient,gauss);
-
-		values[iv]    = reCast<IssmPDouble>(value)/reCast<IssmPDouble>(connectivity[iv]);
-		gradients[iv] = reCast<IssmPDouble>(gradient)/reCast<IssmPDouble>(connectivity[iv]);
-	}
-	delete gauss;
-
-	vector_control->SetValues(NUMVERTICES,&sidlist[0],&values[0],ADD_VAL);
-	vector_gradient->SetValues(NUMVERTICES,&sidlist[0],&gradients[0],ADD_VAL);
-
-}/*}}}*/
-void       Penta::InputControlUpdate(IssmDouble scalar,bool save_parameter){/*{{{*/
-
-	/*Intermediary*/
-	int    num_controls;
-	int*   control_type=NULL;
-	Input* input=NULL;
-
-	/*retrieve some parameters: */
-	this->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
-	this->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
-
-	for(int i=0;i<num_controls;i++){
-
-		if(control_type[i]==MaterialsRheologyBbarEnum){
-			if (!IsOnBase()) goto cleanup_and_return;
-			input=(Input*)this->inputs->GetInput(MaterialsRheologyBEnum); _assert_(input);
-		}
-		else if(control_type[i]==DamageDbarEnum){
-			if (!IsOnBase()) goto cleanup_and_return;
-			input=(Input*)this->inputs->GetInput(DamageDEnum); _assert_(input);
-		}
-		else{
-			input=(Input*)this->inputs->GetInput(control_type[i]); _assert_(input);
-		}
-		if(input->ObjectEnum()!=ControlInputEnum) _error_("input " << EnumToStringx(control_type[i]) << " is not a ControlInput");
-
-		((ControlInput*)input)->UpdateValue(scalar);
-		((ControlInput*)input)->Constrain();
-		if (save_parameter) ((ControlInput*)input)->SaveValue();
-
-		if(control_type[i]==MaterialsRheologyBbarEnum){
-			this->InputExtrude(MaterialsRheologyBEnum,-1);
-		}
-		else if(control_type[i]==DamageDbarEnum){
-			this->InputExtrude(DamageDEnum,-1);
-		}
-	}
-
-	/*Clean up and return*/
-cleanup_and_return:
-	xDelete<int>(control_type);
-}
-/*}}}*/
-void       Penta::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,const char* data,bool onsid){/*{{{*/
-
-	int vertexidlist[NUMVERTICES];
-
-	/*Get out if this is not an element input*/
-	if(!IsInput(control_enum)) return;
-
-	/*Prepare index list*/
-	GradientIndexing(&vertexidlist[0],control_index,onsid);
-
-	/*Get input (either in element or material)*/
-	if(control_enum==MaterialsRheologyBbarEnum) control_enum=MaterialsRheologyBEnum;
-	Input* input=inputs->GetInput(control_enum);
-	if(!input) _error_("Input " << EnumToStringx(control_enum) << " not found in element");
-
-	/*Check that it is a ControlInput*/
-	if (input->ObjectEnum()!=ControlInputEnum){
-		_error_("input " << EnumToStringx(control_enum) << " is not a ControlInput");
-	}
-
-	((ControlInput*)input)->GetVectorFromInputs(vector,&vertexidlist[0],data);
-}
-/*}}}*/
-void       Penta::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){/*{{{*/
-
-	IssmDouble  values[NUMVERTICES];
-	int         vertexpidlist[NUMVERTICES],control_init;
-
-	/*Specific case for depth averaged quantities*/
-	control_init=control_enum;
-	if(control_enum==MaterialsRheologyBbarEnum){
-		control_enum=MaterialsRheologyBEnum;
-		if(!IsOnBase()) return;
-	}
-	if(control_enum==DamageDbarEnum){
-		control_enum=DamageDEnum;
-		if(!IsOnBase()) return;
-	}
-
-	/*Get out if this is not an element input*/
-	if(!IsInput(control_enum)) return;
-
-	/*Prepare index list*/
-	GradientIndexing(&vertexpidlist[0],control_index);
-
-	/*Get values on vertices*/
-	for(int i=0;i<NUMVERTICES;i++){
-		values[i]=vector[vertexpidlist[i]];
-	}
-	Input* new_input = new PentaInput(control_enum,values,P1Enum);
-	Input* input=(Input*)this->inputs->GetInput(control_enum);   _assert_(input);
-	if(input->ObjectEnum()!=ControlInputEnum){
-		_error_("input " << EnumToStringx(control_enum) << " is not a ControlInput");
-	}
-
-	((ControlInput*)input)->SetInput(new_input);
-
-	if(control_init==MaterialsRheologyBbarEnum){
-		this->InputExtrude(control_enum,-1);
-	}
-	if(control_init==DamageDbarEnum){
-		this->InputExtrude(control_enum,-1);
-	}
-}
-/*}}}*/
 
 #ifdef _HAVE_DAKOTA_
@@ -3403,78 +3478,2 @@
 /*}}}*/
 #endif
-void       Penta::GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution, int enum_type){/*{{{*/
-
-	const int    numdof=NDOF1*NUMVERTICES;
-
-	int          i;
-	int*         doflist=NULL;
-	IssmDouble   values[numdof];
-	IssmDouble   enum_value;
-	GaussPenta   *gauss=NULL;
-
-	/*Get dof list: */
-	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-	Input* enum_input=inputs->GetInput(enum_type); _assert_(enum_input);
-
-	gauss=new GaussPenta();
-	for(i=0;i<NUMVERTICES;i++){
-		/*Recover temperature*/
-		gauss->GaussVertex(i);
-		enum_input->GetInputValue(&enum_value,gauss);
-		values[i]=enum_value;
-	}
-
-	/*Add value to global vector*/
-	solution->SetValues(numdof,doflist,values,INS_VAL);
-
-	/*Free ressources:*/
-	delete gauss;
-	xDelete<int>(doflist);
-}
-/*}}}*/
-void       Penta::PotentialUngrounding(Vector<IssmDouble>* potential_ungrounding){/*{{{*/
-
-	IssmDouble  h[NUMVERTICES],r[NUMVERTICES],gl[NUMVERTICES];
-	IssmDouble  bed_hydro;
-	IssmDouble  rho_water,rho_ice,density;
-
-	/*material parameters: */
-	rho_water=matpar->GetRhoWater();
-	rho_ice=matpar->GetRhoIce();
-	density=rho_ice/rho_water;
-	GetInputListOnVertices(&h[0],ThicknessEnum);
-	GetInputListOnVertices(&r[0],BedEnum);
-	GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum);
-
-	/*go through vertices, and figure out which ones are on the ice sheet, and want to unground: */
-	for(int i=0;i<NUMVERTICES;i++){
-		/*Find if grounded vertices want to start floating*/
-		if (gl[i]>0.){
-			bed_hydro=-density*h[i];
-			if(bed_hydro>r[i]){
-				/*Vertex that could potentially unground, flag it*/
-				potential_ungrounding->SetValue(vertices[i]->Pid(),1,INS_VAL);
-			}
-		}
-	}
-}
-/*}}}*/
-int        Penta::UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf){/*{{{*/
-
-	int i;
-	int nflipped=0;
-
-	/*Go through nodes, and whoever is on the potential_ungrounding, ends up in nodes_on_iceshelf: */
-	for(i=0;i<NUMVERTICES;i++){
-		if (reCast<bool,IssmDouble>(vertices_potentially_ungrounding[vertices[i]->Pid()])){
-			vec_nodes_on_iceshelf->SetValue(vertices[i]->Pid(),-1.,INS_VAL);
-
-			/*If node was not on ice shelf, we flipped*/
-			if(nodes_on_iceshelf[vertices[i]->Pid()]>=0.){
-				nflipped++;
-			}
-		}
-	}
-	return nflipped;
-}
-/*}}}*/