Index: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18920)
+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 18921)
@@ -124,12 +124,6 @@
 /*}}}*/
 
 /*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){/*{{{*/
 
 	_assert_(this->inputs);
@@ -155,6 +149,71 @@
 	}
 }
 /*}}}*/
+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;
+
+}
+/*}}}*/
 void       Penta::ComputeBasalStress(Vector<IssmDouble>* sigma_b){/*{{{*/
 
 	int         i,j;
@@ -242,59 +301,6 @@
 	sigma_b->SetValue(id-1,value,INS_VAL);
 }
 /*}}}*/
-void       Penta::ComputeStressTensor(){/*{{{*/
-
-	IssmDouble  xyz_list[NUMVERTICES][3];
-	IssmDouble  pressure,viscosity;
-	IssmDouble  epsilon[6]; /* epsilon=[exx,eyy,exy];*/
-	IssmDouble  sigma_xx[NUMVERTICES];
-	IssmDouble	sigma_yy[NUMVERTICES];
-	IssmDouble	sigma_zz[NUMVERTICES];
-	IssmDouble  sigma_xy[NUMVERTICES];
-	IssmDouble	sigma_xz[NUMVERTICES];
-	IssmDouble	sigma_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* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input);
-	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);
-		pressure_input->GetInputValue(&pressure,gauss);
-
-		/*Compute Stress*/
-		sigma_xx[iv]=2*viscosity*epsilon[0]-pressure; // sigma = nu eps - pressure
-		sigma_yy[iv]=2*viscosity*epsilon[1]-pressure;
-		sigma_zz[iv]=2*viscosity*epsilon[2]-pressure;
-		sigma_xy[iv]=2*viscosity*epsilon[3];
-		sigma_xz[iv]=2*viscosity*epsilon[4];
-		sigma_yz[iv]=2*viscosity*epsilon[5];
-	}
-
-	/*Add Stress tensor components into inputs*/
-	this->inputs->AddInput(new PentaInput(StressTensorxxEnum,&sigma_xx[0],P1Enum));
-	this->inputs->AddInput(new PentaInput(StressTensorxyEnum,&sigma_xy[0],P1Enum));
-	this->inputs->AddInput(new PentaInput(StressTensorxzEnum,&sigma_xz[0],P1Enum));
-	this->inputs->AddInput(new PentaInput(StressTensoryyEnum,&sigma_yy[0],P1Enum));
-	this->inputs->AddInput(new PentaInput(StressTensoryzEnum,&sigma_yz[0],P1Enum));
-	this->inputs->AddInput(new PentaInput(StressTensorzzEnum,&sigma_zz[0],P1Enum));
-
-	/*Clean up and return*/
-	delete gauss;
-}
-/*}}}*/
 void       Penta::ComputeDeviatoricStressTensor(){/*{{{*/
 
 	IssmDouble  xyz_list[NUMVERTICES][3];
@@ -346,232 +352,59 @@
 	delete gauss;
 }
 /*}}}*/
-void       Penta::CalvingRateLevermann(){/*{{{*/
+void       Penta::ComputeStressTensor(){/*{{{*/
 
 	IssmDouble  xyz_list[NUMVERTICES][3];
+	IssmDouble  pressure,viscosity;
+	IssmDouble  epsilon[6]; /* epsilon=[exx,eyy,exy];*/
+	IssmDouble  sigma_xx[NUMVERTICES];
+	IssmDouble	sigma_yy[NUMVERTICES];
+	IssmDouble	sigma_zz[NUMVERTICES];
+	IssmDouble  sigma_xy[NUMVERTICES];
+	IssmDouble	sigma_xz[NUMVERTICES];
+	IssmDouble	sigma_yz[NUMVERTICES];
 	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);
+	/*Retrieve all inputs we will be needing: */
+	Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input);
+	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;
-		strainparallel_input->GetInputValue(&strainparallel,gauss);
-		strainperpendicular_input->GetInputValue(&strainperpendicular,gauss);
-		levermanncoeff_input->GetInputValue(&propcoeff,gauss);
+		/*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);
+		pressure_input->GetInputValue(&pressure,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);
+		/*Compute Stress*/
+		sigma_xx[iv]=2*viscosity*epsilon[0]-pressure; // sigma = nu eps - pressure
+		sigma_yy[iv]=2*viscosity*epsilon[1]-pressure;
+		sigma_zz[iv]=2*viscosity*epsilon[2]-pressure;
+		sigma_xy[iv]=2*viscosity*epsilon[3];
+		sigma_xz[iv]=2*viscosity*epsilon[4];
+		sigma_yz[iv]=2*viscosity*epsilon[5];
 	}
 
-	/*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));
+	/*Add Stress tensor components into inputs*/
+	this->inputs->AddInput(new PentaInput(StressTensorxxEnum,&sigma_xx[0],P1Enum));
+	this->inputs->AddInput(new PentaInput(StressTensorxyEnum,&sigma_xy[0],P1Enum));
+	this->inputs->AddInput(new PentaInput(StressTensorxzEnum,&sigma_xz[0],P1Enum));
+	this->inputs->AddInput(new PentaInput(StressTensoryyEnum,&sigma_yy[0],P1Enum));
+	this->inputs->AddInput(new PentaInput(StressTensoryzEnum,&sigma_yz[0],P1Enum));
+	this->inputs->AddInput(new PentaInput(StressTensorzzEnum,&sigma_zz[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){/*{{{*/
 
 	int analysis_counter;
@@ -605,6 +438,72 @@
 	this->inputs->Configure(parameters);
 }
 /*}}}*/
+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){/*{{{*/
 
 	/*Are we on the base? If not, return*/
@@ -679,6 +578,61 @@
 	delete gauss;
 }
 /*}}}*/
+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;
 }
@@ -769,12 +723,6 @@
 	}
 }
 /*}}}*/
-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*/ 
 
@@ -810,42 +758,11 @@
 	}
 }
 /*}}}*/
-Penta*     Penta::GetUpperPenta(void){/*{{{*/
+Element*   Penta::GetBasalElement(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){/*{{{*/
-
 	/*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;
 }
 /*}}}*/
 Penta*     Penta::GetBasalPenta(void){/*{{{*/
@@ -868,20 +785,12 @@
 	return penta;
 }
 /*}}}*/
-Element*   Penta::GetUpperElement(void){/*{{{*/
+int        Penta::GetElementType(){/*{{{*/
 
-	/*Output*/
-	Element* upper_element=this->GetUpperPenta();
-	return upper_element;
+	/*return PentaRef field*/
+	return this->element_type;
 }
 /*}}}*/
-Element*   Penta::GetBasalElement(void){/*{{{*/
-
-	/*Output*/
-	Element* element=this->GetBasalPenta();
-	return element;
-}
-/*}}}*/
 void       Penta::GetGroundedPart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlyfloating){/*{{{*/
 	/*Computeportion of the element that is grounded*/ 
 
@@ -1036,38 +945,75 @@
 	return phi;
 }
 /*}}}*/
-int        Penta::GetElementType(){/*{{{*/
+void       Penta::GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
+	
+	/* Intermediaries */
+	const int dim=3;
+	int i, dir,nrfrontnodes;
+	IssmDouble  levelset[NUMVERTICES];
 
-	/*return PentaRef field*/
-	return this->element_type;
-}
-/*}}}*/
-void       Penta::ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){/*{{{*/
+	/*Recover parameters and values*/
+	GetInputListOnVertices(&levelset[0],levelsetenum);
 
-	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble xmin,ymin,zmin;
-	IssmDouble xmax,ymax,zmax;
+	int* indicesfront = xNew<int>(NUMVERTICES);
+	/* Get basal nodes where there is no ice */
+	nrfrontnodes=0;
+	for(i=0;i<NUMVERTICES2D;i++){
+		if(levelset[i]>=0.){
+			indicesfront[nrfrontnodes]=i;
+			nrfrontnodes++;
+		}
+	}
+	_assert_(nrfrontnodes==2);
 
-	/*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];
+	/* arrange order of basal frontnodes such that they are oriented counterclockwise */
+	if((NUMVERTICES2D+indicesfront[0]-indicesfront[1])%NUMVERTICES2D!=NUMVERTICES2D-1){
+		int index=indicesfront[0];
+		indicesfront[0]=indicesfront[1];
+		indicesfront[1]=index;
+	}	
 
-	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];
+	IssmDouble* xyz_front = xNew<IssmDouble>(2*dim*nrfrontnodes);
+	/* Return basal and top front nodes */
+	for(i=0;i<nrfrontnodes;i++){
+		for(dir=0;dir<dim;dir++){
+			int ind1=i*dim+dir, ind2=(2*nrfrontnodes-1-i)*dim+dir; // vertex structure front segment: base0, base1, top1, top0
+			xyz_front[ind1]=xyz_list[dim*indicesfront[i]+dir];
+			xyz_front[ind2]=xyz_list[dim*(indicesfront[i]+NUMVERTICES2D)+dir];
+		}
 	}
 
-	*hx=xmax-xmin;
-	*hy=ymax-ymin;
-	*hz=zmax-zmin;
+	*pxyz_front=xyz_front;
+
+	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);
@@ -1092,276 +1038,208 @@
 	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){/*{{{*/
+void       Penta::GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution, int enum_type){/*{{{*/
 
-	Input* input=inputs->GetInput(enumtype);
-	if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria");
+	const int    numdof=NDOF1*NUMVERTICES;
 
-	GaussPenta* gauss=new GaussPenta();
-	gauss->GaussVertex(this->GetNodeIndex(node));
+	int          i;
+	int*         doflist=NULL;
+	IssmDouble   values[numdof];
+	IssmDouble   enum_value;
+	GaussPenta   *gauss=NULL;
 
-	input->GetInputValue(pvalue,gauss);
+	/*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::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){/*{{{*/
+Penta*     Penta::GetSurfacePenta(void){/*{{{*/
 
-	IssmDouble* xyz_list = xNew<IssmDouble>(NUMVERTICES2D*3);
-	::GetVerticesCoordinates(xyz_list,this->vertices,NUMVERTICES2D);
+	/*Output*/
+	Penta* penta=NULL;
 
-	/*Assign output pointer*/
-	*pxyz_list = xyz_list;
+	/*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;
 
-}/*}}}*/
-void       Penta::GetVerticesCoordinatesTop(IssmDouble** pxyz_list){/*{{{*/
+		/* get upper Penta*/
+		penta=penta->GetUpperPenta();
+		_assert_(penta->Id()!=this->id);
+	}
 
-	IssmDouble* xyz_list = xNew<IssmDouble>(NUMVERTICES2D*3);
-	::GetVerticesCoordinates(xyz_list,&this->vertices[3],NUMVERTICES2D);
+	/*return output*/
+	return penta;
+}
+/*}}}*/
+Element*   Penta::GetUpperElement(void){/*{{{*/
 
-	/*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;
+	/*Output*/
+	Element* upper_element=this->GetUpperPenta();
+	return upper_element;
 }
 /*}}}*/
-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*/
+Penta*     Penta::GetUpperPenta(void){/*{{{*/
 
-	IssmDouble normu;
-	IssmDouble tau_parameter;
+	Penta* upper_penta=NULL;
 
-	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);
+	upper_penta=(Penta*)verticalneighbors[1]; //first one (0) under, second one (1) above
 
-	return tau_parameter;
+	return upper_penta;
 }
 /*}}}*/
-void       Penta::ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
-	/*Compute portion of the element that is grounded*/ 
+void       Penta::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,const char* data,bool onsid){/*{{{*/
 
-	int         normal_orientation=0;
-	IssmDouble  s1,s2;
-	IssmDouble  levelset[NUMVERTICES];
-	IssmDouble* xyz_zero = xNew<IssmDouble>(4*3);
+	int vertexidlist[NUMVERTICES];
 
-	/*Recover parameters and values*/
-	GetInputListOnVertices(&levelset[0],levelsetenum);
+	/*Get out if this is not an element input*/
+	if(!IsInput(control_enum)) return;
 
-	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]);
+	/*Prepare index list*/
+	GradientIndexing(&vertexidlist[0],control_index,onsid);
 
-		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]);
+	/*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");
 
-		/*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]);
+	/*Check that it is a ControlInput*/
+	if (input->ObjectEnum()!=ControlInputEnum){
+		_error_("input " << EnumToStringx(control_enum) << " is not a ControlInput");
+	}
 
-		/*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]);
+	((ControlInput*)input)->GetVectorFromInputs(vector,&vertexidlist[0],data);
+}
+/*}}}*/
+void       Penta::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){/*{{{*/
 
-		/*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]);
+	IssmDouble* xyz_list = xNew<IssmDouble>(NUMVERTICES2D*3);
+	::GetVerticesCoordinates(xyz_list,this->vertices,NUMVERTICES2D);
 
-		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]);
+	/*Assign output pointer*/
+	*pxyz_list = xyz_list;
 
-		/*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]);
+}/*}}}*/
+void       Penta::GetVerticesCoordinatesTop(IssmDouble** pxyz_list){/*{{{*/
 
-		/*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]);
+	IssmDouble* xyz_list = xNew<IssmDouble>(NUMVERTICES2D*3);
+	::GetVerticesCoordinates(xyz_list,&this->vertices[3],NUMVERTICES2D);
 
-		/*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]);
+	/*Assign output pointer*/
+	*pxyz_list = xyz_list;
 
-		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]);
+}/*}}}*/
+IssmDouble Penta::IceVolume(void){/*{{{*/
 
-		/*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]);
+	/*The volume of a troncated prism is base * 1/3 sum(length of edges)*/
+	IssmDouble base,height;
+	IssmDouble xyz_list[NUMVERTICES][3];
 
-		/*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]);
+	if(!IsIceInElement())return 0;
 
-		/*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];
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 
-		/*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];
+	/*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]));
 
-		/*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];
+	/*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]));
 
-		/*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];
+	/*Return: */
+	return base*height;
+}
+/*}}}*/
+IssmDouble Penta::IceVolumeAboveFloatation(void){/*{{{*/
 
-		/*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];
+	/*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];
 
-		/*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];
+	if(!IsIceInElement() || IsFloating() || !IsOnBase())return 0;
 
-		/*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];
+	rho_ice=matpar->GetRhoIce();
+	rho_water=matpar->GetRhoWater();
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 
-		/*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];
+	/*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]));
 
-		/*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];
+	/*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);
 
-		/*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;
+	/*Return: */
+	return base*(surface - bed + min( rho_water/rho_ice * bathymetry, 0.) );
 }
 /*}}}*/
-void       Penta::GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
-	
-	/* Intermediaries */
-	const int dim=3;
-	int i, dir,nrfrontnodes;
-	IssmDouble  levelset[NUMVERTICES];
+void       Penta::InputControlUpdate(IssmDouble scalar,bool save_parameter){/*{{{*/
 
-	/*Recover parameters and values*/
-	GetInputListOnVertices(&levelset[0],levelsetenum);
+	/*Intermediary*/
+	int    num_controls;
+	int*   control_type=NULL;
+	Input* input=NULL;
 
-	int* indicesfront = xNew<int>(NUMVERTICES);
-	/* Get basal nodes where there is no ice */
-	nrfrontnodes=0;
-	for(i=0;i<NUMVERTICES2D;i++){
-		if(levelset[i]>=0.){
-			indicesfront[nrfrontnodes]=i;
-			nrfrontnodes++;
+	/*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);
 		}
-	}
-	_assert_(nrfrontnodes==2);
+		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");
 
-	/* arrange order of basal frontnodes such that they are oriented counterclockwise */
-	if((NUMVERTICES2D+indicesfront[0]-indicesfront[1])%NUMVERTICES2D!=NUMVERTICES2D-1){
-		int index=indicesfront[0];
-		indicesfront[0]=indicesfront[1];
-		indicesfront[1]=index;
-	}	
+		((ControlInput*)input)->UpdateValue(scalar);
+		((ControlInput*)input)->Constrain();
+		if (save_parameter) ((ControlInput*)input)->SaveValue();
 
-	IssmDouble* xyz_front = xNew<IssmDouble>(2*dim*nrfrontnodes);
-	/* Return basal and top front nodes */
-	for(i=0;i<nrfrontnodes;i++){
-		for(dir=0;dir<dim;dir++){
-			int ind1=i*dim+dir, ind2=(2*nrfrontnodes-1-i)*dim+dir; // vertex structure front segment: base0, base1, top1, top0
-			xyz_front[ind1]=xyz_list[dim*indicesfront[i]+dir];
-			xyz_front[ind2]=xyz_list[dim*(indicesfront[i]+NUMVERTICES2D)+dir];
+		if(control_type[i]==MaterialsRheologyBbarEnum){
+			this->InputExtrude(MaterialsRheologyBEnum,-1);
 		}
+		else if(control_type[i]==DamageDbarEnum){
+			this->InputExtrude(DamageDEnum,-1);
+		}
 	}
 
-	*pxyz_front=xyz_front;
-
-	xDelete<int>(indicesfront);
-}/*}}}*/
+	/*Clean up and return*/
+cleanup_and_return:
+	xDelete<int>(control_type);
+}
+/*}}}*/
 void       Penta::InputDepthAverageAtBase(int enum_type,int average_enum_type){/*{{{*/
 
 	int  step,i;
@@ -1731,6 +1609,39 @@
 	}
 }
 /*}}}*/
+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){/*{{{*/
 
 	IssmDouble values[NUMVERTICES];
@@ -1767,18 +1678,22 @@
 	}
 }
 /*}}}*/
-bool       Penta::IsNodeOnShelfFromFlags(IssmDouble* flags){/*{{{*/
+bool       Penta::IsZeroLevelset(int levelset_enum){/*{{{*/
 
-	int  i;
-	bool shelf=false;
+	bool        iszerols;
+	IssmDouble  ls[NUMVERTICES];
 
-	for(i=0;i<NUMVERTICES;i++){
-		if (flags[vertices[i]->Pid()]<0.){
-			shelf=true;
-			break;
+	/*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 shelf;
+	return iszerols;
 }
 /*}}}*/
 void       Penta::JacobianDeterminant(IssmDouble* pJdet,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
@@ -1802,6 +1717,13 @@
 
 }
 /*}}}*/
+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){/*{{{*/
 
 	_assert_(gauss->Enum()==GaussPentaEnum);
@@ -1809,13 +1731,52 @@
 
 }
 /*}}}*/
-void       Penta::JacobianDeterminantSurface(IssmDouble* pJdet,IssmDouble* xyz_list_quad,Gauss* gauss){/*{{{*/
+IssmDouble Penta::MassFlux(IssmDouble* segment){/*{{{*/
 
-	_assert_(gauss->Enum()==GaussPentaEnum);
-	this->GetQuadJacobianDeterminant(pJdet,xyz_list_quad,(GaussPenta*)gauss);
+	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;
+}
+/*}}}*/
 IssmDouble Penta::MinEdgeLength(IssmDouble* xyz_list){/*{{{*/
 	/*Return the minimum lenght of the nine egdes of the penta*/
 
@@ -1837,37 +1798,6 @@
 	return minlength;
 }
 /*}}}*/
-int        Penta::NodalValue(IssmDouble* pvalue, int index, int natureofdataenum){/*{{{*/
-
-	int i;
-	int found=0;
-	IssmDouble value;
-	Input* data=NULL;
-	GaussPenta* gauss=NULL;
-
-	/*First, serarch the input: */
-	data=inputs->GetInput(natureofdataenum); 
-
-	/*figure out if we have the vertex id: */
-	found=0;
-	for(i=0;i<NUMVERTICES;i++){
-		if(index==vertices[i]->Id()){
-			/*Do we have natureofdataenum in our inputs? :*/
-			if(data){
-				/*ok, we are good. retrieve value of input at vertex :*/
-				gauss=new GaussPenta(); gauss->GaussVertex(i);
-				data->GetInputValue(&value,gauss);
-				found=1;
-				break;
-			}
-		}
-	}
-
-	delete gauss;
-	if(found)*pvalue=value;
-	return found;
-}
-/*}}}*/
 Gauss*     Penta::NewGauss(void){/*{{{*/
 	return new GaussPenta();
 }
@@ -1904,59 +1834,59 @@
 
 }
 /*}}}*/
-void       Penta::NodalFunctionsP1(IssmDouble* basis, Gauss* gauss){/*{{{*/
+void       Penta::NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 
 	_assert_(gauss->Enum()==GaussPentaEnum);
-	this->GetNodalFunctions(basis,(GaussPenta*)gauss,P1Enum);
+	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss,this->element_type);
 
 }
 /*}}}*/
-void       Penta::NodalFunctionsP2(IssmDouble* basis, Gauss* gauss){/*{{{*/
+void       Penta::NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 
 	_assert_(gauss->Enum()==GaussPentaEnum);
-	this->GetNodalFunctions(basis,(GaussPenta*)gauss,P2Enum);
+	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss,this->VelocityInterpolation());
 
 }
 /*}}}*/
-void       Penta::NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void       Penta::NodalFunctionsMINIDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 
 	_assert_(gauss->Enum()==GaussPentaEnum);
-	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss,this->element_type);
+	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss,P1bubbleEnum);
 
 }
 /*}}}*/
-void       Penta::NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void       Penta::NodalFunctionsPressure(IssmDouble* basis, Gauss* gauss){/*{{{*/
 
 	_assert_(gauss->Enum()==GaussPentaEnum);
-	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss,P1Enum);
+	this->GetNodalFunctions(basis,(GaussPenta*)gauss,this->PressureInterpolation());
 
 }
 /*}}}*/
-void       Penta::NodalFunctionsMINIDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void       Penta::NodalFunctionsP1(IssmDouble* basis, Gauss* gauss){/*{{{*/
 
 	_assert_(gauss->Enum()==GaussPentaEnum);
-	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss,P1bubbleEnum);
+	this->GetNodalFunctions(basis,(GaussPenta*)gauss,P1Enum);
 
 }
 /*}}}*/
-void       Penta::NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void       Penta::NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 
 	_assert_(gauss->Enum()==GaussPentaEnum);
-	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss,this->VelocityInterpolation());
+	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussPenta*)gauss,P1Enum);
 
 }
 /*}}}*/
-void       Penta::NodalFunctionsVelocity(IssmDouble* basis, Gauss* gauss){/*{{{*/
+void       Penta::NodalFunctionsP2(IssmDouble* basis, Gauss* gauss){/*{{{*/
 
 	_assert_(gauss->Enum()==GaussPentaEnum);
-	this->GetNodalFunctions(basis,(GaussPenta*)gauss,this->VelocityInterpolation());
+	this->GetNodalFunctions(basis,(GaussPenta*)gauss,P2Enum);
 
 }
 /*}}}*/
-void       Penta::NodalFunctionsPressure(IssmDouble* basis, Gauss* gauss){/*{{{*/
+void       Penta::NodalFunctionsVelocity(IssmDouble* basis, Gauss* gauss){/*{{{*/
 
 	_assert_(gauss->Enum()==GaussPentaEnum);
-	this->GetNodalFunctions(basis,(GaussPenta*)gauss,this->PressureInterpolation());
+	this->GetNodalFunctions(basis,(GaussPenta*)gauss,this->VelocityInterpolation());
 
 }
 /*}}}*/
@@ -1967,6 +1897,37 @@
 
 }
 /*}}}*/
+int        Penta::NodalValue(IssmDouble* pvalue, int index, int natureofdataenum){/*{{{*/
+
+	int i;
+	int found=0;
+	IssmDouble value;
+	Input* data=NULL;
+	GaussPenta* gauss=NULL;
+
+	/*First, serarch the input: */
+	data=inputs->GetInput(natureofdataenum); 
+
+	/*figure out if we have the vertex id: */
+	found=0;
+	for(i=0;i<NUMVERTICES;i++){
+		if(index==vertices[i]->Id()){
+			/*Do we have natureofdataenum in our inputs? :*/
+			if(data){
+				/*ok, we are good. retrieve value of input at vertex :*/
+				gauss=new GaussPenta(); gauss->GaussVertex(i);
+				data->GetInputValue(&value,gauss);
+				found=1;
+				break;
+			}
+		}
+	}
+
+	delete gauss;
+	if(found)*pvalue=value;
+	return found;
+}
+/*}}}*/
 void       Penta::NormalBase(IssmDouble* bed_normal,IssmDouble* xyz_list){/*{{{*/
 
 	IssmDouble v13[3],v23[3];
@@ -1989,6 +1950,26 @@
 	bed_normal[2]=-normal[2]/normal_norm;
 }
 /*}}}*/
+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){/*{{{*/
 
 	int i;
@@ -2019,6 +2000,12 @@
 	return PentaRef::NumberofNodes(this->VelocityInterpolation());
 }
 /*}}}*/
+int        Penta::ObjectEnum(void){/*{{{*/
+
+	return PentaEnum;
+
+}
+/*}}}*/
 void       Penta::PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm){/*{{{*/
 
    IssmDouble agd[NUMVERTICES];             // surface mass balance
@@ -2088,6 +2075,37 @@
 	delete gauss;
 }
 /*}}}*/
+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){/*{{{*/
 
 	int analysis_type;
@@ -2206,11 +2224,71 @@
 	delete gauss;
 }
 /*}}}*/
+void       Penta::ResetHooks(){/*{{{*/
+
+	this->nodes=NULL;
+	this->vertices=NULL;
+	this->material=NULL;
+	this->matpar=NULL;
+	this->verticalneighbors=NULL;
+	this->parameters=NULL;
+
+	//deal with ElementHook mother class
+	for(int i=0;i<this->numanalyses;i++) if(this->hnodes[i]) this->hnodes[i]->reset();
+	this->hvertices->reset();
+	this->hmaterial->reset();
+	this->hmatpar->reset();
+	if(this->hneighbors) this->hneighbors->reset();
+
+}
+/*}}}*/
 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: */
@@ -2226,26 +2304,38 @@
 
 }
 /*}}}*/
-void       Penta::ResetHooks(){/*{{{*/
+void       Penta::SetTemporaryElementType(int element_type_in){/*{{{*/
+	this->element_type=element_type_in;
+}
+/*}}}*/
+Element*   Penta::SpawnBasalElement(void){/*{{{*/
 
-	this->nodes=NULL;
-	this->vertices=NULL;
-	this->material=NULL;
-	this->matpar=NULL;
-	this->verticalneighbors=NULL;
-	this->parameters=NULL;
+	_assert_(this->IsOnBase());
 
-	//deal with ElementHook mother class
-	for(int i=0;i<this->numanalyses;i++) if(this->hnodes[i]) this->hnodes[i]->reset();
-	this->hvertices->reset();
-	this->hmaterial->reset();
-	this->hmatpar->reset();
-	if(this->hneighbors) this->hneighbors->reset();
+	this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum);
+	if(this->material->IsDamage())this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum);
+	if(this->inputs->GetInput(VxEnum)) this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
+	if(this->inputs->GetInput(VyEnum)) this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
+	if(this->inputs->GetInput(CalvingratexEnum)) this->InputDepthAverageAtBase(CalvingratexEnum,CalvingratexAverageEnum);
+	if(this->inputs->GetInput(CalvingrateyEnum)) this->InputDepthAverageAtBase(CalvingrateyEnum,CalvingrateyAverageEnum);
+	Tria* tria=(Tria*)SpawnTria(0,1,2);
+	this->inputs->DeleteInput(MaterialsRheologyBbarEnum);
+	this->inputs->DeleteInput(DamageDbarEnum);
+	this->inputs->DeleteInput(VxAverageEnum);
+	this->inputs->DeleteInput(VyAverageEnum);
+	this->inputs->DeleteInput(CalvingratexAverageEnum);
+	this->inputs->DeleteInput(CalvingrateyAverageEnum);
 
+	return tria;
 }
 /*}}}*/
-void       Penta::SetTemporaryElementType(int element_type_in){/*{{{*/
-	this->element_type=element_type_in;
+Element*   Penta::SpawnTopElement(void){/*{{{*/
+
+	_assert_(this->IsOnSurface());
+
+	Tria* tria=(Tria*)SpawnTria(3,4,5);
+
+	return tria;
 }
 /*}}}*/
 Tria*      Penta::SpawnTria(int index1,int index2,int index3){/*{{{*/
@@ -2275,36 +2365,194 @@
 	return tria;
 }
 /*}}}*/
-Element*   Penta::SpawnBasalElement(void){/*{{{*/
+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*/
 
-	_assert_(this->IsOnBase());
+	IssmDouble normu;
+	IssmDouble tau_parameter;
 
-	this->InputDepthAverageAtBase(MaterialsRheologyBEnum,MaterialsRheologyBbarEnum);
-	if(this->material->IsDamage())this->InputDepthAverageAtBase(DamageDEnum,DamageDbarEnum);
-	if(this->inputs->GetInput(VxEnum)) this->InputDepthAverageAtBase(VxEnum,VxAverageEnum);
-	if(this->inputs->GetInput(VyEnum)) this->InputDepthAverageAtBase(VyEnum,VyAverageEnum);
-	if(this->inputs->GetInput(CalvingratexEnum)) this->InputDepthAverageAtBase(CalvingratexEnum,CalvingratexAverageEnum);
-	if(this->inputs->GetInput(CalvingrateyEnum)) this->InputDepthAverageAtBase(CalvingrateyEnum,CalvingrateyAverageEnum);
-	Tria* tria=(Tria*)SpawnTria(0,1,2);
-	this->inputs->DeleteInput(MaterialsRheologyBbarEnum);
-	this->inputs->DeleteInput(DamageDbarEnum);
-	this->inputs->DeleteInput(VxAverageEnum);
-	this->inputs->DeleteInput(VyAverageEnum);
-	this->inputs->DeleteInput(CalvingratexAverageEnum);
-	this->inputs->DeleteInput(CalvingrateyAverageEnum);
+	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 tria;
+	return tau_parameter;
 }
 /*}}}*/
-Element*   Penta::SpawnTopElement(void){/*{{{*/
+void       Penta::StrainRateparallel(){/*{{{*/
 
-	_assert_(this->IsOnSurface());
+	IssmDouble *xyz_list = NULL;
+	IssmDouble  epsilon[6];
+	GaussPenta* gauss=NULL;
+	IssmDouble  vx,vy,vel;
+	IssmDouble  strainxx;
+	IssmDouble  strainxy;
+	IssmDouble  strainyy;
+	IssmDouble  strainparallel[NUMVERTICES];
 
-	Tria* tria=(Tria*)SpawnTria(3,4,5);
+	/* Get node coordinates and dof list: */
+	this->GetVerticesCoordinates(&xyz_list);
 
-	return tria;
+	/*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){/*{{{*/
 
 	int    approximation;
@@ -2384,6 +2632,35 @@
 
 	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){ /*{{{*/
 
 	/*Intermediaries*/
@@ -2832,14 +3109,37 @@
 
 }
 /*}}}*/
-void       Penta::ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss){/*{{{*/
-	PentaRef::GetInputValue(pvalue,values,gauss,P1Enum);
+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);
+}
+/*}}}*/
+int        Penta::VelocityInterpolation(void){/*{{{*/
+	return PentaRef::VelocityInterpolation(this->element_type);
+}
+/*}}}*/
 int        Penta::VertexConnectivity(int vertexindex){/*{{{*/
 	_assert_(this->vertices);
 	return this->vertices[vertexindex]->Connectivity();
@@ -2874,393 +3174,168 @@
 	*pphi = phi;
 }
 /*}}}*/
-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){/*{{{*/
+void       Penta::ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
+	/*Compute portion of the element that is grounded*/ 
 
-	bool        iszerols;
-	IssmDouble  ls[NUMVERTICES];
+	int         normal_orientation=0;
+	IssmDouble  s1,s2;
+	IssmDouble  levelset[NUMVERTICES];
+	IssmDouble* xyz_zero = xNew<IssmDouble>(4*3);
 
-	/*Retrieve all inputs and parameters*/
-	GetInputListOnVertices(&ls[0],levelset_enum);
+	/*Recover parameters and values*/
+	GetInputListOnVertices(&levelset[0],levelsetenum);
 
-	/*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){/*{{{*/
+	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]);
 
-	bool isicefront;
-	int i,nrice;
-   IssmDouble ls[NUMVERTICES];
+		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]);
 
-	/*Retrieve all inputs and parameters*/
-	GetInputListOnVertices(&ls[0],MaskIceLevelsetEnum);
+		/*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]);
 
-	/* 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){/*{{{*/
+		/*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]);
 
-	/*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!");
+		/*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]);
 
-}
-/*}}}*/
-IssmDouble Penta::TotalSmb(void){/*{{{*/
+		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]);
 
-	/*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];
+		/*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]);
 
-	/*Get material parameters :*/
-	rho_ice=matpar->GetRhoIce();
+		/*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]);
 
-	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;
-}
-/*}}}*/
-
-#ifdef _HAVE_GIA_
-void       Penta::GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,IssmDouble* x,IssmDouble* y){/*{{{*/
-	_error_("GIA deflection not implemented yet!");
-}
-/*}}}*/
-#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);
+		/*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(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");
+	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]);
 
-	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);
+		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]);
 
-}/*}}}*/
-void       Penta::ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum){/*{{{*/
+		/*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]);
 
-	Input* input=NULL;
+		/*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]);
 
-	if(control_enum==MaterialsRheologyBbarEnum){
-		input=(Input*)inputs->GetInput(MaterialsRheologyBEnum);
+		/*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(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");
+	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];
 
-	int         sidlist[NUMVERTICES];
-	int         connectivity[NUMVERTICES];
-	IssmPDouble values[NUMVERTICES];
-	IssmPDouble gradients[NUMVERTICES]; 
-	IssmDouble  value,gradient;
+		/*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];
 
-	this->GetVerticesConnectivityList(&connectivity[0]);
-	this->GetVerticesSidList(&sidlist[0]);
+		/*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];
 
-	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]);
+		/*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];
 	}
-	delete gauss;
+	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];
 
-	vector_control->SetValues(NUMVERTICES,&sidlist[0],&values[0],ADD_VAL);
-	vector_gradient->SetValues(NUMVERTICES,&sidlist[0],&gradients[0],ADD_VAL);
+		/*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];
 
-}/*}}}*/
-void       Penta::InputControlUpdate(IssmDouble scalar,bool save_parameter){/*{{{*/
+		/*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];
 
-	/*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);
-		}
+		/*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];
 
-	/*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){/*{{{*/
+		/*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];
 
-	int vertexidlist[NUMVERTICES];
+		/*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];
 
-	/*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");
+		/*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");
 
-	((ControlInput*)input)->GetVectorFromInputs(vector,&vertexidlist[0],data);
+	/*Assign output pointer*/
+	*pxyz_zero= xyz_zero;
 }
 /*}}}*/
-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_GIA_
+void       Penta::GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,IssmDouble* x,IssmDouble* y){/*{{{*/
+	_error_("GIA deflection not implemented yet!");
 }
 /*}}}*/
+#endif
 
 #ifdef _HAVE_DAKOTA_
 void       Penta::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){/*{{{*/
@@ -3402,79 +3477,3 @@
 }
 /*}}}*/
 #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;
-}
-/*}}}*/
