Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18910)
+++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 18911)
@@ -154,6 +154,99 @@
 	this->inputs->AddInput(new TriaInput(input_enum,values,interpolation_enum));
 }
 /*}}}*/
+void       Tria::AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part){/*{{{*/
+
+	bool       already = false;
+	int        i,j;
+	int        partition[NUMVERTICES];
+	int        offsetsid[NUMVERTICES];
+	int        offsetdof[NUMVERTICES];
+	IssmDouble area;
+	IssmDouble mean;
+
+	/*First, get the area: */
+	area=this->GetArea();
+
+	/*Figure out the average for this element: */
+	this->GetVerticesSidList(&offsetsid[0]);
+	this->GetVertexPidList(&offsetdof[0]);
+	mean=0;
+	for(i=0;i<NUMVERTICES;i++){
+		partition[i]=reCast<int>(qmu_part[offsetsid[i]]);
+		mean=mean+1.0/NUMVERTICES*vertex_response[offsetdof[i]];
+	}
+
+	/*Add contribution: */
+	for(i=0;i<NUMVERTICES;i++){
+		already=false;
+		for(j=0;j<i;j++){
+			if (partition[i]==partition[j]){
+				already=true;
+				break;
+			}
+		}
+		if(!already){
+			partition_contributions->SetValue(partition[i],mean*area,ADD_VAL);
+			partition_areas->SetValue(partition[i],area,ADD_VAL);
+		};
+	}
+}
+/*}}}*/
+void       Tria::CalvingRateLevermann(){/*{{{*/
+
+	IssmDouble  xyz_list[NUMVERTICES][3];
+	GaussTria* 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 GaussTria();
+	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 TriaInput(CalvingratexEnum,&calvingratex[0],P1Enum));
+	this->inputs->AddInput(new TriaInput(CalvingrateyEnum,&calvingratey[0],P1Enum));
+	this->inputs->AddInput(new TriaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum));
+
+	/*Clean up and return*/
+	delete gauss;
+
+}
+/*}}}*/
 IssmDouble Tria::CharacteristicLength(void){/*{{{*/
 
 	return sqrt(2*this->GetArea());
@@ -163,6 +256,56 @@
 	_error_("Not Implemented yet");
 }
 /*}}}*/
+void       Tria::ComputeDeviatoricStressTensor(){/*{{{*/
+
+	IssmDouble  xyz_list[NUMVERTICES][3];
+	IssmDouble  viscosity;
+	IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
+	IssmDouble  tau_xx[NUMVERTICES];
+	IssmDouble	tau_yy[NUMVERTICES];
+	IssmDouble	tau_zz[NUMVERTICES]={0,0,0};
+	IssmDouble  tau_xy[NUMVERTICES];
+	IssmDouble	tau_xz[NUMVERTICES]={0,0,0};
+	IssmDouble	tau_yz[NUMVERTICES]={0,0,0};
+	GaussTria*  gauss=NULL;
+	int domaintype,dim=2;
+
+	/* Get node coordinates and dof list: */
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+
+	/*Retrieve all inputs we will be needing: */
+	this->FindParam(&domaintype,DomainTypeEnum);
+	if(domaintype!=Domain2DhorizontalEnum) _error_("deviatoric stress tensor calculation not implemented for mesh of type " <<EnumToStringx(domaintype));
+	Input* vx_input=inputs->GetInput(VxEnum);             _assert_(vx_input);
+	Input* vy_input=inputs->GetInput(VyEnum);             _assert_(vy_input);
+
+	/* Start looping on the number of vertices: */
+	gauss=new GaussTria();
+	for (int iv=0;iv<NUMVERTICES;iv++){
+		gauss->GaussVertex(iv);
+
+		/*Compute strain rate and viscosity: */
+		this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
+		this->ViscositySSA(&viscosity,dim,&xyz_list[0][0],gauss,vx_input,vy_input);
+
+		/*Compute Stress*/
+		tau_xx[iv]=2*viscosity*epsilon[0]; // tau = nu eps
+		tau_yy[iv]=2*viscosity*epsilon[1];
+		tau_xy[iv]=2*viscosity*epsilon[2];
+	}
+
+	/*Add Stress tensor components into inputs*/
+	this->inputs->AddInput(new TriaInput(DeviatoricStressxxEnum,&tau_xx[0],P1Enum));
+	this->inputs->AddInput(new TriaInput(DeviatoricStressxyEnum,&tau_xy[0],P1Enum));
+	this->inputs->AddInput(new TriaInput(DeviatoricStressxzEnum,&tau_xz[0],P1Enum));
+	this->inputs->AddInput(new TriaInput(DeviatoricStressyyEnum,&tau_yy[0],P1Enum));
+	this->inputs->AddInput(new TriaInput(DeviatoricStressyzEnum,&tau_yz[0],P1Enum));
+	this->inputs->AddInput(new TriaInput(DeviatoricStresszzEnum,&tau_zz[0],P1Enum));
+
+	/*Clean up and return*/
+	delete gauss;
+}
+/*}}}*/
 void       Tria::ComputeSigmaNN(){/*{{{*/
 
 	if(!IsOnBase()){
@@ -275,203 +418,6 @@
 	delete gauss;
 }
 /*}}}*/
-void       Tria::ComputeDeviatoricStressTensor(){/*{{{*/
-
-	IssmDouble  xyz_list[NUMVERTICES][3];
-	IssmDouble  viscosity;
-	IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
-	IssmDouble  tau_xx[NUMVERTICES];
-	IssmDouble	tau_yy[NUMVERTICES];
-	IssmDouble	tau_zz[NUMVERTICES]={0,0,0};
-	IssmDouble  tau_xy[NUMVERTICES];
-	IssmDouble	tau_xz[NUMVERTICES]={0,0,0};
-	IssmDouble	tau_yz[NUMVERTICES]={0,0,0};
-	GaussTria*  gauss=NULL;
-	int domaintype,dim=2;
-
-	/* Get node coordinates and dof list: */
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-
-	/*Retrieve all inputs we will be needing: */
-	this->FindParam(&domaintype,DomainTypeEnum);
-	if(domaintype!=Domain2DhorizontalEnum) _error_("deviatoric stress tensor calculation not implemented for mesh of type " <<EnumToStringx(domaintype));
-	Input* vx_input=inputs->GetInput(VxEnum);             _assert_(vx_input);
-	Input* vy_input=inputs->GetInput(VyEnum);             _assert_(vy_input);
-
-	/* Start looping on the number of vertices: */
-	gauss=new GaussTria();
-	for (int iv=0;iv<NUMVERTICES;iv++){
-		gauss->GaussVertex(iv);
-
-		/*Compute strain rate and viscosity: */
-		this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
-		this->ViscositySSA(&viscosity,dim,&xyz_list[0][0],gauss,vx_input,vy_input);
-
-		/*Compute Stress*/
-		tau_xx[iv]=2*viscosity*epsilon[0]; // tau = nu eps
-		tau_yy[iv]=2*viscosity*epsilon[1];
-		tau_xy[iv]=2*viscosity*epsilon[2];
-	}
-
-	/*Add Stress tensor components into inputs*/
-	this->inputs->AddInput(new TriaInput(DeviatoricStressxxEnum,&tau_xx[0],P1Enum));
-	this->inputs->AddInput(new TriaInput(DeviatoricStressxyEnum,&tau_xy[0],P1Enum));
-	this->inputs->AddInput(new TriaInput(DeviatoricStressxzEnum,&tau_xz[0],P1Enum));
-	this->inputs->AddInput(new TriaInput(DeviatoricStressyyEnum,&tau_yy[0],P1Enum));
-	this->inputs->AddInput(new TriaInput(DeviatoricStressyzEnum,&tau_yz[0],P1Enum));
-	this->inputs->AddInput(new TriaInput(DeviatoricStresszzEnum,&tau_zz[0],P1Enum));
-
-	/*Clean up and return*/
-	delete gauss;
-}
-/*}}}*/
-void       Tria::StrainRateparallel(){/*{{{*/
-
-	IssmDouble *xyz_list = NULL;
-	IssmDouble  epsilon[3];
-	GaussTria* 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);
-
-	/* Start looping on the number of vertices: */
-	gauss=new GaussTria();
-	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 viscosity and pressure: */
-		this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
-		strainxx=epsilon[0];
-		strainyy=epsilon[1];
-		strainxy=epsilon[2];
-
-		/*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 TriaInput(StrainRateparallelEnum,&strainparallel[0],P1Enum));
-
-	/*Clean up and return*/
-	delete gauss;
-	xDelete<IssmDouble>(xyz_list);
-}
-/*}}}*/
-void       Tria::StrainRateperpendicular(){/*{{{*/
-
-	IssmDouble *xyz_list = NULL;
-	GaussTria* gauss=NULL;
-	IssmDouble  epsilon[3];
-	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);
-
-	/* Start looping on the number of vertices: */
-	gauss=new GaussTria();
-	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 viscosity and pressure: */
-		this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
-		strainxx=epsilon[0];
-		strainyy=epsilon[1];
-		strainxy=epsilon[2];
-
-		/*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 TriaInput(StrainRateperpendicularEnum,&strainperpendicular[0],P1Enum));
-
-	/*Clean up and return*/
-	delete gauss;
-	xDelete<IssmDouble>(xyz_list);
-}
-/*}}}*/
-void       Tria::CalvingRateLevermann(){/*{{{*/
-
-	IssmDouble  xyz_list[NUMVERTICES][3];
-	GaussTria* 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 GaussTria();
-	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 TriaInput(CalvingratexEnum,&calvingratex[0],P1Enum));
-	this->inputs->AddInput(new TriaInput(CalvingrateyEnum,&calvingratey[0],P1Enum));
-	this->inputs->AddInput(new TriaInput(CalvingCalvingrateEnum,&calvingrate[0],P1Enum));
-
-	/*Clean up and return*/
-	delete gauss;
-
-}
-/*}}}*/
 void       Tria::Configure(Elements* elementsin, Loads* loadsin,Nodes* nodesin,Vertices *verticesin,Materials* materialsin, Parameters* parametersin){/*{{{*/
 
 	/*go into parameters and get the analysis_counter: */
@@ -507,6 +453,54 @@
 
 }
 /*}}}*/
+void       Tria::ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){/*{{{*/
+
+	int    vertexpidlist[NUMVERTICES];
+	IssmDouble grad_list[NUMVERTICES];
+	Input* grad_input=NULL;
+
+	Input* 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 TriaInput(GradientEnum,grad_list,P1Enum);
+
+	((ControlInput*)input)->SetGradient(grad_input);
+
+}/*}}}*/
+void       Tria::ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum){/*{{{*/
+
+	Input* 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]);
+
+	GaussTria* gauss=new GaussTria();
+	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       Tria::Delta18oParameterization(void){/*{{{*/
 
 	int        i;
@@ -577,6 +571,107 @@
 	delete gauss;
 }
 /*}}}*/
+int        Tria::EdgeOnBaseIndex(void){/*{{{*/
+
+	IssmDouble values[NUMVERTICES];
+	int        indices[3][2] = {{1,2},{2,0},{0,1}};
+
+	/*Retrieve all inputs and parameters*/
+	GetInputListOnVertices(&values[0],MeshVertexonbaseEnum);
+
+	for(int i=0;i<3;i++){
+		if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){
+			return i;
+		}
+	}
+
+	_printf_("list of vertices on bed: "<<values[0]<<" "<<values[1]<<" "<<values[2]);
+	_error_("Could not find 2 vertices on bed");
+}
+/*}}}*/
+void       Tria::EdgeOnBaseIndices(int* pindex1,int* pindex2){/*{{{*/
+
+	IssmDouble values[NUMVERTICES];
+	int        indices[3][2] = {{1,2},{2,0},{0,1}};
+
+	/*Retrieve all inputs and parameters*/
+	GetInputListOnVertices(&values[0],MeshVertexonbaseEnum);
+
+	for(int i=0;i<3;i++){
+		if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){
+			*pindex1 = indices[i][0];
+			*pindex2 = indices[i][1];
+			return;
+		}
+	}
+
+	_printf_("list of vertices on bed: "<<values[0]<<" "<<values[1]<<" "<<values[2]);
+	_error_("Could not find 2 vertices on bed");
+}
+/*}}}*/
+int        Tria::EdgeOnSurfaceIndex(void){/*{{{*/
+
+	IssmDouble values[NUMVERTICES];
+	int        indices[3][2] = {{1,2},{2,0},{0,1}};
+
+	/*Retrieve all inputs and parameters*/
+	GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
+
+	for(int i=0;i<3;i++){
+		if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){
+			return i;
+		}
+	}
+
+	_printf_("list of vertices on surface: "<<values[0]<<" "<<values[1]<<" "<<values[2]);
+	_error_("Could not find 2 vertices on surface");
+}
+/*}}}*/
+void       Tria::EdgeOnSurfaceIndices(int* pindex1,int* pindex2){/*{{{*/
+
+	IssmDouble values[NUMVERTICES];
+	int        indices[3][2] = {{1,2},{2,0},{0,1}};
+
+	/*Retrieve all inputs and parameters*/
+	GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
+
+	for(int i=0;i<3;i++){
+		if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){
+			*pindex1 = indices[i][0];
+			*pindex2 = indices[i][1];
+			return;
+		}
+	}
+
+	_printf_("list of vertices on surface: "<<values[0]<<" "<<values[1]<<" "<<values[2]);
+	_error_("Could not find 2 vertices on surface");
+}
+/*}}}*/
+void       Tria::ElementResponse(IssmDouble* presponse,int response_enum){/*{{{*/
+
+	switch(response_enum){
+		case MaterialsRheologyBbarEnum:
+			*presponse=this->material->GetBbar();
+			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       Tria::ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz){/*{{{*/
 
 	IssmDouble xyz_list[NUMVERTICES][3];
@@ -604,10 +699,90 @@
 	return this->element_type;
 }
 /*}}}*/
-int        Tria::ObjectEnum(void){/*{{{*/
+void       Tria::FSContactMigration(Vector<IssmDouble>* vertexgrounded,Vector<IssmDouble>* vertexfloating){/*{{{*/
 
-	return TriaEnum;
+	if(!IsOnBase()) return;
 
+	int approximation;
+	inputs->GetInputValue(&approximation,ApproximationEnum);
+
+	if(approximation==HOApproximationEnum || approximation==SSAApproximationEnum || approximation==SSAHOApproximationEnum){
+		for(int i=0;i<NUMVERTICES;i++){
+			vertexgrounded->SetValue(vertices[i]->Pid(),+9999.,INS_VAL);
+			vertexfloating->SetValue(vertices[i]->Pid(),+9999.,INS_VAL);
+		}
+	}
+	else{
+		/*Intermediaries*/
+		IssmDouble* xyz_list = NULL;
+		IssmDouble* xyz_list_base = NULL;
+		IssmDouble  pressure,water_pressure,sigma_nn,viscosity,bed,base;
+		IssmDouble  bed_normal[2];
+		IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
+		IssmDouble  surface=0,value=0;
+		bool grounded;
+
+		/* Get node coordinates and dof list: */
+		GetVerticesCoordinates(&xyz_list);
+		GetVerticesCoordinatesBase(&xyz_list_base);
+
+		/*Retrieve all inputs we will be needing: */
+		Input* pressure_input = inputs->GetInput(PressureEnum); _assert_(pressure_input);
+		Input* base_input     = inputs->GetInput(BaseEnum);     _assert_(base_input);
+		Input* bed_input      = inputs->GetInput(BedEnum);      _assert_(bed_input);
+		Input* vx_input       = inputs->GetInput(VxEnum);       _assert_(vx_input);
+		Input* vy_input       = inputs->GetInput(VyEnum);       _assert_(vy_input);
+
+		/*Create gauss point in the middle of the basal edge*/
+		Gauss* gauss=NewGaussBase(1);
+		gauss->GaussPoint(0);
+
+		if(!IsFloating()){ 
+			/*Check for basal force only if grounded and touching GL*/
+			//		if(this->inputs->Min(MaskGroundediceLevelsetEnum)==0.){
+			this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
+			this->ViscosityFS(&viscosity,2,xyz_list,gauss,vx_input,vy_input,NULL);
+			pressure_input->GetInputValue(&pressure, gauss);
+			base_input->GetInputValue(&base, gauss); 
+
+			/*Compute Stress*/
+			IssmDouble sigma_xx=2.*viscosity*epsilon[0]-pressure;
+			IssmDouble sigma_yy=2.*viscosity*epsilon[1]-pressure;
+			IssmDouble sigma_xy=2.*viscosity*epsilon[2];
+
+			/*Get normal vector to the bed */
+			NormalBase(&bed_normal[0],xyz_list_base);
+
+			/*basalforce*/
+			sigma_nn = sigma_xx*bed_normal[0]*bed_normal[0] + sigma_yy*bed_normal[1]*bed_normal[1] + 2.*sigma_xy*bed_normal[0]*bed_normal[1];
+
+			/*Compute water pressure*/
+			IssmDouble rho_ice   = matpar->GetRhoIce();
+			IssmDouble rho_water = matpar->GetRhoWater();
+			IssmDouble gravity   = matpar->GetG();
+			water_pressure=gravity*rho_water*base;
+
+			/*Compare basal stress to water pressure and determine whether it should ground*/
+			if (sigma_nn<water_pressure) grounded=true;
+			else                         grounded=false;
+		}
+		else{
+			/*Check for basal elevation if floating*/
+			base_input->GetInputValue(&base, gauss);
+			bed_input->GetInputValue(&bed, gauss);
+			if(base<bed) grounded=true;
+			else         grounded=false;
+		}
+		for(int i=0;i<NUMVERTICES;i++){
+			if(grounded) vertexgrounded->SetValue(vertices[i]->Pid(),+1.,INS_VAL);
+			else         vertexfloating->SetValue(vertices[i]->Pid(),+1.,INS_VAL);
+		}
+
+		/*clean up*/
+		delete gauss;
+		xDelete<IssmDouble>(xyz_list);
+		xDelete<IssmDouble>(xyz_list_base);
+	}
 }
 /*}}}*/
 IssmDouble Tria::GetArea(void){/*{{{*/
@@ -667,56 +842,6 @@
 
 }
 /*}}}*/
-void        Tria::GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* gl){/*{{{*/
-	
-	/*Computeportion of the element that has a positive levelset*/ 
-
-	bool               negative=true;
-	int                point;
-	const IssmPDouble  epsilon= 1.e-15;
-	IssmDouble         f1,f2;
-
-	/*Be sure that values are not zero*/
-	if(gl[0]==0.) gl[0]=gl[0]+epsilon;
-	if(gl[1]==0.) gl[1]=gl[1]+epsilon;
-	if(gl[2]==0.) gl[2]=gl[2]+epsilon;
-
-	/*Check that not all nodes are positive or negative*/
-	if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All positive
-		point=0;
-		f1=1.;
-		f2=1.;
-	}
-	else if(gl[0]<0 && gl[1]<0 && gl[2]<0){ //All negative
-		point=0;
-		f1=0.;
-		f2=0.;
-	}
-	else{
-		if(gl[0]*gl[1]*gl[2]<0) negative=false;
-
-		if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
-			point=2;
-			f1=gl[2]/(gl[2]-gl[0]);
-			f2=gl[2]/(gl[2]-gl[1]);
-		}
-		else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
-			point=0;
-			f1=gl[0]/(gl[0]-gl[1]);
-			f2=gl[0]/(gl[0]-gl[2]);
-		}
-		else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
-			point=1;
-			f1=gl[1]/(gl[1]-gl[2]);
-			f2=gl[1]/(gl[1]-gl[0]);
-		}
-	}
-	*point1=point;
-	*fraction1=f1;
-	*fraction2=f2;
-	*mainlynegative=negative;
-}
-/*}}}*/
 void       Tria::GetGroundedPart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlyfloating){/*{{{*/
 	/*Computeportion of the element that is grounded*/ 
 
@@ -888,149 +1013,6 @@
 	return phi;
 }
 /*}}}*/
-void       Tria::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){/*{{{*/
-
-	int        indices[2];
-	IssmDouble xyz_list[NUMVERTICES][3];
-
-	/*Element XYZ list*/
-	::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
-
-	/*Allocate Output*/
-	IssmDouble* xyz_list_edge = xNew<IssmDouble>(2*3);
-	this->EdgeOnBaseIndices(&indices[0],&indices[1]);
-	for(int i=0;i<2;i++) for(int j=0;j<2;j++) xyz_list_edge[i*3+j]=xyz_list[indices[i]][j];
-
-	/*Assign output pointer*/
-	*pxyz_list = xyz_list_edge;
-
-}/*}}}*/
-void       Tria::GetVerticesCoordinatesTop(IssmDouble** pxyz_list){/*{{{*/
-
-	int        indices[2];
-	IssmDouble xyz_list[NUMVERTICES][3];
-
-	/*Element XYZ list*/
-	::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
-
-	/*Allocate Output*/
-	IssmDouble* xyz_list_edge = xNew<IssmDouble>(2*3);
-	this->EdgeOnSurfaceIndices(&indices[0],&indices[1]);
-	for(int i=0;i<2;i++) for(int j=0;j<2;j++) xyz_list_edge[i*3+j]=xyz_list[indices[i]][j];
-
-	/*Assign output pointer*/
-	*pxyz_list = xyz_list_edge;
-
-}/*}}}*/
-void       Tria::NormalSection(IssmDouble* normal,IssmDouble* xyz_list){/*{{{*/
-
-	/*Build unit outward pointing vector*/
-	IssmDouble vector[2];
-	IssmDouble norm;
-
-	vector[0]=xyz_list[1*3+0] - xyz_list[0*3+0];
-	vector[1]=xyz_list[1*3+1] - xyz_list[0*3+1];
-
-	norm=sqrt(vector[0]*vector[0] + vector[1]*vector[1]);
-
-	normal[0]= + vector[1]/norm;
-	normal[1]= - vector[0]/norm;
-}
-/*}}}*/
-void       Tria::ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
-
-	int         normal_orientation=0;
-	IssmDouble  s1,s2;
-	IssmDouble  levelset[NUMVERTICES];
-
-	/*Recover parameters and values*/
-	IssmDouble* xyz_zero = xNew<IssmDouble>(2*3);
-	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 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]);
-	}
-	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]);
-	}
-	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]);
-	}
-	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];
-	}
-	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];
-	}
-	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];
-	}
-	else _error_("Case not covered");
-
-	/*Assign output pointer*/
-	*pxyz_zero= xyz_zero;
-}
-/*}}}*/
 void       Tria::GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
 	
 	/* Intermediaries */
@@ -1071,6 +1053,18 @@
 
 	xDelete<int>(indicesfront);
 }/*}}}*/
+void       Tria::GetInputValue(IssmDouble* pvalue,Node* node,int enumtype){/*{{{*/
+
+	Input* input=inputs->GetInput(enumtype);
+	if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria");
+
+	GaussTria* gauss=new GaussTria();
+	gauss->GaussVertex(this->GetNodeIndex(node));
+
+	input->GetInputValue(pvalue,gauss);
+	delete gauss;
+}
+/*}}}*/
 void       Tria::GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level){/*{{{*/
 
 	/* Intermediaries */
@@ -1111,6 +1105,62 @@
 
 	xDelete<int>(indicesfront);
 }/*}}}*/
+void       Tria::GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* gl){/*{{{*/
+	
+	/*Computeportion of the element that has a positive levelset*/ 
+
+	bool               negative=true;
+	int                point;
+	const IssmPDouble  epsilon= 1.e-15;
+	IssmDouble         f1,f2;
+
+	/*Be sure that values are not zero*/
+	if(gl[0]==0.) gl[0]=gl[0]+epsilon;
+	if(gl[1]==0.) gl[1]=gl[1]+epsilon;
+	if(gl[2]==0.) gl[2]=gl[2]+epsilon;
+
+	/*Check that not all nodes are positive or negative*/
+	if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All positive
+		point=0;
+		f1=1.;
+		f2=1.;
+	}
+	else if(gl[0]<0 && gl[1]<0 && gl[2]<0){ //All negative
+		point=0;
+		f1=0.;
+		f2=0.;
+	}
+	else{
+		if(gl[0]*gl[1]*gl[2]<0) negative=false;
+
+		if(gl[0]*gl[1]>0){ //Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2
+			point=2;
+			f1=gl[2]/(gl[2]-gl[0]);
+			f2=gl[2]/(gl[2]-gl[1]);
+		}
+		else if(gl[1]*gl[2]>0){ //Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2
+			point=0;
+			f1=gl[0]/(gl[0]-gl[1]);
+			f2=gl[0]/(gl[0]-gl[2]);
+		}
+		else if(gl[0]*gl[2]>0){ //Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2
+			point=1;
+			f1=gl[1]/(gl[1]-gl[2]);
+			f2=gl[1]/(gl[1]-gl[0]);
+		}
+	}
+	*point1=point;
+	*fraction1=f1;
+	*fraction2=f2;
+	*mainlynegative=negative;
+}
+/*}}}*/
+Node*      Tria::GetNode(int node_number){/*{{{*/
+	_assert_(node_number>=0); 
+	_assert_(node_number<this->NumberofNodes(this->element_type)); 
+	return this->nodes[node_number];
+
+}/*}}}*/
 int        Tria::GetNodeIndex(Node* node){/*{{{*/
 
 	_assert_(nodes);
@@ -1134,24 +1184,221 @@
 	return NUMVERTICES;
 }
 /*}}}*/
-void       Tria::GetInputValue(IssmDouble* pvalue,Node* node,int enumtype){/*{{{*/
+void       Tria::GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution, int enum_type){/*{{{*/
 
-	Input* input=inputs->GetInput(enumtype);
-	if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria");
+	int        *doflist = NULL;
+	IssmDouble  value;
 
+	/*Fetch number of nodes for this finite element*/
+	int numnodes = this->NumberofNodes(this->element_type);
+
+	/*Fetch dof list and allocate solution vector*/
+	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
+	IssmDouble* values = xNew<IssmDouble>(numnodes);
+
+	/*Get inputs*/
+	Input* enum_input=inputs->GetInput(enum_type); _assert_(enum_input);
+
+	/*Ok, we have the values, fill in the array: */
 	GaussTria* gauss=new GaussTria();
-	gauss->GaussVertex(this->GetNodeIndex(node));
+	for(int i=0;i<numnodes;i++){
+		gauss->GaussNode(this->element_type,i);
 
-	input->GetInputValue(pvalue,gauss);
+		enum_input->GetInputValue(&value,gauss);
+		values[i]=value;
+	}
+
+	solution->SetValues(numnodes,doflist,values,INS_VAL);
+
+	/*Free ressources:*/
+	xDelete<int>(doflist);
+	xDelete<IssmDouble>(values);
 	delete gauss;
 }
 /*}}}*/
-Node*      Tria::GetNode(int node_number){/*{{{*/
-	_assert_(node_number>=0); 
-	_assert_(node_number<this->NumberofNodes(this->element_type)); 
-	return this->nodes[node_number];
+void       Tria::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,const char* data,bool onsid){/*{{{*/
 
+	int vertexidlist[NUMVERTICES];
+	Input *input=NULL;
+
+	/*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)*/
+	input=(Input*)this->inputs->GetInput(control_enum);   _assert_(input);
+
+	/*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       Tria::GetVerticesCoordinatesBase(IssmDouble** pxyz_list){/*{{{*/
+
+	int        indices[2];
+	IssmDouble xyz_list[NUMVERTICES][3];
+
+	/*Element XYZ list*/
+	::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
+
+	/*Allocate Output*/
+	IssmDouble* xyz_list_edge = xNew<IssmDouble>(2*3);
+	this->EdgeOnBaseIndices(&indices[0],&indices[1]);
+	for(int i=0;i<2;i++) for(int j=0;j<2;j++) xyz_list_edge[i*3+j]=xyz_list[indices[i]][j];
+
+	/*Assign output pointer*/
+	*pxyz_list = xyz_list_edge;
+
 }/*}}}*/
+void       Tria::GetVerticesCoordinatesTop(IssmDouble** pxyz_list){/*{{{*/
+
+	int        indices[2];
+	IssmDouble xyz_list[NUMVERTICES][3];
+
+	/*Element XYZ list*/
+	::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
+
+	/*Allocate Output*/
+	IssmDouble* xyz_list_edge = xNew<IssmDouble>(2*3);
+	this->EdgeOnSurfaceIndices(&indices[0],&indices[1]);
+	for(int i=0;i<2;i++) for(int j=0;j<2;j++) xyz_list_edge[i*3+j]=xyz_list[indices[i]][j];
+
+	/*Assign output pointer*/
+	*pxyz_list = xyz_list_edge;
+
+}/*}}}*/
+bool       Tria::HasEdgeOnBase(){/*{{{*/
+
+	IssmDouble values[NUMVERTICES];
+	IssmDouble sum;
+
+	/*Retrieve all inputs and parameters*/
+	GetInputListOnVertices(&values[0],MeshVertexonbaseEnum);
+	sum = values[0]+values[1]+values[2];
+
+	_assert_(sum==0. || sum==1. || sum==2.);
+
+	if(sum==3.)  _error_("Two edges on bed not supported yet...");
+
+	if(sum>1.){
+		return true;
+	}
+	else{
+		return false;
+	}
+}
+/*}}}*/
+bool       Tria::HasEdgeOnSurface(){/*{{{*/
+
+	IssmDouble values[NUMVERTICES];
+	IssmDouble sum;
+
+	/*Retrieve all inputs and parameters*/
+	GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
+	sum = values[0]+values[1]+values[2];
+
+	_assert_(sum==0. || sum==1. || sum==2.);
+
+	if(sum==3.)  _error_("Two edges on surface not supported yet...");
+
+	if(sum>1.){
+		return true;
+	}
+	else{
+		return false;
+	}
+}
+/*}}}*/
+IssmDouble Tria::IceVolume(void){/*{{{*/
+
+	/*The volume of a troncated prism is base * 1/3 sum(length of edges)*/
+	IssmDouble base,surface,bed;
+	IssmDouble xyz_list[NUMVERTICES][3];
+
+	if(!IsIceInElement())return 0;
+
+	/*First get back the area of the base*/
+	base=this->GetArea();
+
+	/*Now get the average height*/
+	Input* surface_input = inputs->GetInput(SurfaceEnum); _assert_(surface_input);
+	Input* base_input     = inputs->GetInput(BaseEnum);     _assert_(base_input);
+	surface_input->GetInputAverage(&surface);
+	base_input->GetInputAverage(&bed);
+
+	/*Return: */
+	int domaintype;
+	parameters->FindParam(&domaintype,DomainTypeEnum);
+	if(domaintype==Domain2DverticalEnum){
+	  return base;
+	}
+	else{
+	  return base*(surface-bed);
+	}
+}
+/*}}}*/
+IssmDouble Tria::IceVolumeAboveFloatation(void){/*{{{*/
+
+	/*The volume above floatation: H + rho_water/rho_ice * bathymetry */
+	IssmDouble rho_ice,rho_water;
+	IssmDouble base,surface,bed,bathymetry;
+	IssmDouble xyz_list[NUMVERTICES][3];
+
+	if(!IsIceInElement() || IsFloating())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/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 height and bathymetry*/
+	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       Tria::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++){
+		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();
+
+	}
+
+	/*Clean up and return*/
+	xDelete<int>(control_type);
+}
+/*}}}*/
 void       Tria::InputDepthAverageAtBase(int enum_type,int average_enum_type){/*{{{*/
 
 	/*New input*/
@@ -1370,6 +1617,59 @@
 
 }
 /*}}}*/
+bool       Tria::IsFaceOnBoundary(void){/*{{{*/
+
+	IssmDouble values[NUMVERTICES];
+	IssmDouble sum;
+
+	/*Retrieve all inputs and parameters*/
+	GetInputListOnVertices(&values[0],MeshVertexonboundaryEnum);
+	sum = values[0]+values[1]+values[2];
+
+	_assert_(sum==0. || sum==1. || sum==2.);
+
+	if(sum==3.)  _error_("Two edges on boundary not supported yet...");
+
+	if(sum>1.){
+		return true;
+	}
+	else{
+		return false;
+	}
+}/*}}}*/
+bool       Tria::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<NUMVERTICES;i++)
+			if(ls[i]<0.) nrice++;
+		if(nrice==1) isicefront= true;
+	}
+	return isicefront;
+}/*}}}*/
+bool       Tria::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       Tria::IsOnBase(){/*{{{*/
 
 	int domaintype;
@@ -1396,6 +1696,25 @@
 	}
 }
 /*}}}*/
+bool       Tria::IsZeroLevelset(int levelset_enum){/*{{{*/
+
+	bool iszerols;
+	IssmDouble ls[NUMVERTICES];
+
+	/*Retrieve all inputs and parameters*/
+	GetInputListOnVertices(&ls[0],levelset_enum);
+
+	/*If the level set is awlays <0, 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;
+}
+/*}}}*/
 void       Tria::JacobianDeterminant(IssmDouble* pJdet,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 
 	_assert_(gauss->Enum()==GaussTriaEnum);
@@ -1424,222 +1743,267 @@
 
 }
 /*}}}*/
-bool       Tria::HasEdgeOnBase(){/*{{{*/
+IssmDouble Tria::Masscon(IssmDouble* levelset){ /*{{{*/
 
-	IssmDouble values[NUMVERTICES];
-	IssmDouble sum;
 
-	/*Retrieve all inputs and parameters*/
-	GetInputListOnVertices(&values[0],MeshVertexonbaseEnum);
-	sum = values[0]+values[1]+values[2];
+	/*intermediary: */
+	IssmDouble* values=NULL;
+	Input*      thickness_input=NULL;
+	IssmDouble  thickness;
+	IssmDouble  weight;
+	IssmDouble  Jdet;
+	IssmDouble  volume;
+	IssmDouble  rho_ice;
+	IssmDouble* xyz_list=NULL;
+	int         point1;
+	IssmDouble  fraction1,fraction2;
+	bool        mainlynegative=true;
+	
+	/*Output:*/
+	volume=0;
 
-	_assert_(sum==0. || sum==1. || sum==2.);
+	/* Get node coordinates and dof list: */
+	GetVerticesCoordinates(&xyz_list);
 
-	if(sum==3.)  _error_("Two edges on bed not supported yet...");
+	/*Retrieve inputs required:*/
+	thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input);
+	
+	/*Retrieve material parameters: */
+	rho_ice=matpar->GetRhoIce();
 
-	if(sum>1.){
-		return true;
+	/*Retrieve values of the levelset defining the masscon: */
+	values = xNew<IssmDouble>(NUMVERTICES);
+	for(int i=0;i<NUMVERTICES;i++){
+		values[i]=levelset[this->vertices[i]->Sid()];
 	}
-	else{
-		return false;
+		
+	/*Ok, use the level set values to figure out where we put our gaussian points:*/
+	this->GetLevelsetPositivePart(&point1,&fraction1,&fraction2,&mainlynegative,values);
+	Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,mainlynegative,4);
+
+	volume=0;
+
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
+		gauss->GaussPoint(ig);
+
+		this->JacobianDeterminant(&Jdet,xyz_list,gauss);
+		thickness_input->GetInputValue(&thickness, gauss);
+
+		volume+=thickness*gauss->weight*Jdet;
 	}
+
+	/* clean up and Return: */
+	xDelete<IssmDouble>(xyz_list);
+	xDelete<IssmDouble>(values);
+	delete gauss;
+	return rho_ice*volume;
 }
 /*}}}*/
-bool       Tria::HasEdgeOnSurface(){/*{{{*/
+IssmDouble Tria::MassFlux(IssmDouble x1, IssmDouble y1, IssmDouble x2, IssmDouble y2,int segment_id){/*{{{*/
 
-	IssmDouble values[NUMVERTICES];
-	IssmDouble sum;
+	int        domaintype;
+	IssmDouble mass_flux=0.;
+	IssmDouble xyz_list[NUMVERTICES][3];
+	IssmDouble normal[2];
+	IssmDouble length,rho_ice;
+	IssmDouble h1,h2;
+	IssmDouble vx1,vx2,vy1,vy2;
+	GaussTria* gauss_1=NULL;
+	GaussTria* gauss_2=NULL;
 
-	/*Retrieve all inputs and parameters*/
-	GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
-	sum = values[0]+values[1]+values[2];
+	/*Get material parameters :*/
+	rho_ice=matpar->GetRhoIce();
 
-	_assert_(sum==0. || sum==1. || sum==2.);
+	/*First off, check that this segment belongs to this element: */
+	if (segment_id!=this->id)_error_("error message: segment with id " << segment_id << " does not belong to element with id:" << this->id);
 
-	if(sum==3.)  _error_("Two edges on surface not supported yet...");
+	/*Get xyz list: */
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 
-	if(sum>1.){
-		return true;
+	/*get area coordinates of 0 and 1 locations: */
+	gauss_1=new GaussTria();
+	gauss_1->GaussFromCoords(x1,y1,&xyz_list[0][0]);
+	gauss_2=new GaussTria();
+	gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]);
+
+	normal[0]=cos(atan2(x1-x2,y2-y1));
+	normal[1]=sin(atan2(x1-x2,y2-y1));
+
+	length=sqrt(pow(x2-x1,2)+pow(y2-y1,2));
+
+	Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
+	this->parameters->FindParam(&domaintype,DomainTypeEnum);
+	Input* vx_input=NULL;
+	Input* vy_input=NULL;
+	if(domaintype==Domain2DhorizontalEnum){
+		vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
+		vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
 	}
 	else{
-		return false;
+		vx_input=inputs->GetInput(VxAverageEnum); _assert_(vx_input);
+		vy_input=inputs->GetInput(VyAverageEnum); _assert_(vy_input);
 	}
-}
-/*}}}*/
-void       Tria::EdgeOnBaseIndices(int* pindex1,int* pindex2){/*{{{*/
 
-	IssmDouble values[NUMVERTICES];
-	int        indices[3][2] = {{1,2},{2,0},{0,1}};
+	thickness_input->GetInputValue(&h1, gauss_1);
+	thickness_input->GetInputValue(&h2, gauss_2);
+	vx_input->GetInputValue(&vx1,gauss_1);
+	vx_input->GetInputValue(&vx2,gauss_2);
+	vy_input->GetInputValue(&vy1,gauss_1);
+	vy_input->GetInputValue(&vy2,gauss_2);
 
-	/*Retrieve all inputs and parameters*/
-	GetInputListOnVertices(&values[0],MeshVertexonbaseEnum);
+	mass_flux= rho_ice*length*(  
+				(ONETHIRD*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*normal[0]+
+				(ONETHIRD*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*normal[1]
+				);
 
-	for(int i=0;i<3;i++){
-		if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){
-			*pindex1 = indices[i][0];
-			*pindex2 = indices[i][1];
-			return;
-		}
-	}
-
-	_printf_("list of vertices on bed: "<<values[0]<<" "<<values[1]<<" "<<values[2]);
-	_error_("Could not find 2 vertices on bed");
+	/*clean up and return:*/
+	delete gauss_1;
+	delete gauss_2;
+	return mass_flux;
 }
 /*}}}*/
-void       Tria::EdgeOnSurfaceIndices(int* pindex1,int* pindex2){/*{{{*/
+IssmDouble Tria::MassFlux(IssmDouble* segment){/*{{{*/
 
-	IssmDouble values[NUMVERTICES];
-	int        indices[3][2] = {{1,2},{2,0},{0,1}};
+	int        domaintype;
+	IssmDouble mass_flux=0.;
+	IssmDouble xyz_list[NUMVERTICES][3];
+	IssmDouble normal[2];
+	IssmDouble length,rho_ice;
+	IssmDouble x1,y1,x2,y2,h1,h2;
+	IssmDouble vx1,vx2,vy1,vy2;
+	GaussTria* gauss_1=NULL;
+	GaussTria* gauss_2=NULL;
 
-	/*Retrieve all inputs and parameters*/
-	GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
+	/*Get material parameters :*/
+	rho_ice=matpar->GetRhoIce();
 
-	for(int i=0;i<3;i++){
-		if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){
-			*pindex1 = indices[i][0];
-			*pindex2 = indices[i][1];
-			return;
-		}
-	}
+	/*First off, check that this segment belongs to this element: */
+	if (reCast<int>(*(segment+4))!=this->id)_error_("error message: segment with id " << reCast<int>(*(segment+4)) << " does not belong to element with id:" << this->id);
 
-	_printf_("list of vertices on surface: "<<values[0]<<" "<<values[1]<<" "<<values[2]);
-	_error_("Could not find 2 vertices on surface");
-}
-/*}}}*/
-int        Tria::EdgeOnBaseIndex(void){/*{{{*/
+	/*Recover segment node locations: */
+	x1=*(segment+0); y1=*(segment+1); x2=*(segment+2); y2=*(segment+3);
 
-	IssmDouble values[NUMVERTICES];
-	int        indices[3][2] = {{1,2},{2,0},{0,1}};
+	/*Get xyz list: */
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 
-	/*Retrieve all inputs and parameters*/
-	GetInputListOnVertices(&values[0],MeshVertexonbaseEnum);
+	/*get area coordinates of 0 and 1 locations: */
+	gauss_1=new GaussTria();
+	gauss_1->GaussFromCoords(x1,y1,&xyz_list[0][0]);
+	gauss_2=new GaussTria();
+	gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]);
 
-	for(int i=0;i<3;i++){
-		if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){
-			return i;
-		}
+	normal[0]=cos(atan2(x1-x2,y2-y1));
+	normal[1]=sin(atan2(x1-x2,y2-y1));
+
+	length=sqrt(pow(x2-x1,2)+pow(y2-y1,2));
+
+	Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
+	this->parameters->FindParam(&domaintype,DomainTypeEnum);
+	Input* vx_input=NULL;
+	Input* vy_input=NULL;
+	if(domaintype==Domain2DhorizontalEnum){
+		vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
+		vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
 	}
+	else{
+		vx_input=inputs->GetInput(VxAverageEnum); _assert_(vx_input);
+		vy_input=inputs->GetInput(VyAverageEnum); _assert_(vy_input);
+	}
 
-	_printf_("list of vertices on bed: "<<values[0]<<" "<<values[1]<<" "<<values[2]);
-	_error_("Could not find 2 vertices on bed");
+	thickness_input->GetInputValue(&h1, gauss_1);
+	thickness_input->GetInputValue(&h2, gauss_2);
+	vx_input->GetInputValue(&vx1,gauss_1);
+	vx_input->GetInputValue(&vx2,gauss_2);
+	vy_input->GetInputValue(&vy1,gauss_1);
+	vy_input->GetInputValue(&vy2,gauss_2);
+
+	mass_flux= rho_ice*length*(  
+				(ONETHIRD*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*normal[0]+
+				(ONETHIRD*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*normal[1]
+				);
+
+	/*clean up and return:*/
+	delete gauss_1;
+	delete gauss_2;
+	return mass_flux;
 }
 /*}}}*/
-int        Tria::EdgeOnSurfaceIndex(void){/*{{{*/
+IssmDouble Tria::Misfit(int modelenum,int observationenum,int weightsenum){/*{{{*/
 
-	IssmDouble values[NUMVERTICES];
-	int        indices[3][2] = {{1,2},{2,0},{0,1}};
+	/*Intermediaries*/
+	IssmDouble model,observation,weight;
+	IssmDouble Jdet;
+	IssmDouble Jelem = 0;
+	IssmDouble xyz_list[NUMVERTICES][3];
+	GaussTria *gauss = NULL;
 
-	/*Retrieve all inputs and parameters*/
-	GetInputListOnVertices(&values[0],MeshVertexonsurfaceEnum);
+	/*If on water, return 0: */
+	if(!IsIceInElement())return 0;
 
-	for(int i=0;i<3;i++){
-		if(values[indices[i][0]] == 1. && values[indices[i][1]] == 1.){
-			return i;
-		}
-	}
+	/*Retrieve all inputs we will be needing: */
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+	Input* model_input=inputs->GetInput(modelenum);   _assert_(model_input);
+	Input* observation_input=inputs->GetInput(observationenum);_assert_(observation_input);
+	Input* weights_input     =inputs->GetInput(weightsenum);     _assert_(weights_input);
 
-	_printf_("list of vertices on surface: "<<values[0]<<" "<<values[1]<<" "<<values[2]);
-	_error_("Could not find 2 vertices on surface");
-}
-/*}}}*/
-void       Tria::FSContactMigration(Vector<IssmDouble>* vertexgrounded,Vector<IssmDouble>* vertexfloating){/*{{{*/
+	/* Start  looping on the number of gaussian points: */
+	gauss=new GaussTria(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
 
-	if(!IsOnBase()) return;
+		gauss->GaussPoint(ig);
 
-	int approximation;
-	inputs->GetInputValue(&approximation,ApproximationEnum);
+		/* Get Jacobian determinant: */
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
 
-	if(approximation==HOApproximationEnum || approximation==SSAApproximationEnum || approximation==SSAHOApproximationEnum){
-		for(int i=0;i<NUMVERTICES;i++){
-			vertexgrounded->SetValue(vertices[i]->Pid(),+9999.,INS_VAL);
-			vertexfloating->SetValue(vertices[i]->Pid(),+9999.,INS_VAL);
-		}
+		/*Get parameters at gauss point*/
+		model_input->GetInputValue(&model,gauss);
+		observation_input->GetInputValue(&observation,gauss);
+		weights_input->GetInputValue(&weight,gauss);
+
+		/*compute misfit between model and observation */
+		Jelem+=0.5*(model-observation)*(model-observation)*Jdet*weight*gauss->weight;
 	}
-	else{
-		/*Intermediaries*/
-		IssmDouble* xyz_list = NULL;
-		IssmDouble* xyz_list_base = NULL;
-		IssmDouble  pressure,water_pressure,sigma_nn,viscosity,bed,base;
-		IssmDouble  bed_normal[2];
-		IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
-		IssmDouble  surface=0,value=0;
-		bool grounded;
 
-		/* Get node coordinates and dof list: */
-		GetVerticesCoordinates(&xyz_list);
-		GetVerticesCoordinatesBase(&xyz_list_base);
+	/* clean up and Return: */
+	delete gauss;
+	return Jelem;
+}
+/*}}}*/
+IssmDouble Tria::MisfitArea(int weightsenum){/*{{{*/
 
-		/*Retrieve all inputs we will be needing: */
-		Input* pressure_input = inputs->GetInput(PressureEnum); _assert_(pressure_input);
-		Input* base_input     = inputs->GetInput(BaseEnum);     _assert_(base_input);
-		Input* bed_input      = inputs->GetInput(BedEnum);      _assert_(bed_input);
-		Input* vx_input       = inputs->GetInput(VxEnum);       _assert_(vx_input);
-		Input* vy_input       = inputs->GetInput(VyEnum);       _assert_(vy_input);
+	/*Intermediaries*/
+	IssmDouble weight;
+	IssmDouble Jdet;
+	IssmDouble Jelem = 0;
+	IssmDouble xyz_list[NUMVERTICES][3];
+	GaussTria *gauss = NULL;
 
-		/*Create gauss point in the middle of the basal edge*/
-		Gauss* gauss=NewGaussBase(1);
-		gauss->GaussPoint(0);
+	/*If on water, return 0: */
+	if(!IsIceInElement())return 0;
 
-		if(!IsFloating()){ 
-			/*Check for basal force only if grounded and touching GL*/
-			//		if(this->inputs->Min(MaskGroundediceLevelsetEnum)==0.){
-			this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
-			this->ViscosityFS(&viscosity,2,xyz_list,gauss,vx_input,vy_input,NULL);
-			pressure_input->GetInputValue(&pressure, gauss);
-			base_input->GetInputValue(&base, gauss); 
+	/*Retrieve all inputs we will be needing: */
+	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+	Input* weights_input     =inputs->GetInput(weightsenum);     _assert_(weights_input);
 
-			/*Compute Stress*/
-			IssmDouble sigma_xx=2.*viscosity*epsilon[0]-pressure;
-			IssmDouble sigma_yy=2.*viscosity*epsilon[1]-pressure;
-			IssmDouble sigma_xy=2.*viscosity*epsilon[2];
+	/* Start  looping on the number of gaussian points: */
+	gauss=new GaussTria(2);
+	for(int ig=gauss->begin();ig<gauss->end();ig++){
 
-			/*Get normal vector to the bed */
-			NormalBase(&bed_normal[0],xyz_list_base);
+		gauss->GaussPoint(ig);
 
-			/*basalforce*/
-			sigma_nn = sigma_xx*bed_normal[0]*bed_normal[0] + sigma_yy*bed_normal[1]*bed_normal[1] + 2.*sigma_xy*bed_normal[0]*bed_normal[1];
+		/* Get Jacobian determinant: */
+		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
 
-			/*Compute water pressure*/
-			IssmDouble rho_ice   = matpar->GetRhoIce();
-			IssmDouble rho_water = matpar->GetRhoWater();
-			IssmDouble gravity   = matpar->GetG();
-			water_pressure=gravity*rho_water*base;
+		/*Get parameters at gauss point*/
+		weights_input->GetInputValue(&weight,gauss);
 
-			/*Compare basal stress to water pressure and determine whether it should ground*/
-			if (sigma_nn<water_pressure) grounded=true;
-			else                         grounded=false;
-		}
-		else{
-			/*Check for basal elevation if floating*/
-			base_input->GetInputValue(&base, gauss);
-			bed_input->GetInputValue(&bed, gauss);
-			if(base<bed) grounded=true;
-			else         grounded=false;
-		}
-		for(int i=0;i<NUMVERTICES;i++){
-			if(grounded) vertexgrounded->SetValue(vertices[i]->Pid(),+1.,INS_VAL);
-			else         vertexfloating->SetValue(vertices[i]->Pid(),+1.,INS_VAL);
-		}
-
-		/*clean up*/
-		delete gauss;
-		xDelete<IssmDouble>(xyz_list);
-		xDelete<IssmDouble>(xyz_list_base);
+		/*compute misfit between model and observation */
+		Jelem+=Jdet*weight*gauss->weight;
 	}
-}
-/*}}}*/
-bool       Tria::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;
+	/* clean up and Return: */
+	delete gauss;
+	return Jelem;
 }
 /*}}}*/
 Gauss*     Tria::NewGauss(void){/*{{{*/
@@ -1690,59 +2054,59 @@
 
 }
 /*}}}*/
-void       Tria::NodalFunctionsP1(IssmDouble* basis, Gauss* gauss){/*{{{*/
+void       Tria::NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 
 	_assert_(gauss->Enum()==GaussTriaEnum);
-	this->GetNodalFunctions(basis,(GaussTria*)gauss,P1Enum);
+	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss,this->element_type);
 
 }
 /*}}}*/
-void       Tria::NodalFunctionsP2(IssmDouble* basis, Gauss* gauss){/*{{{*/
+void       Tria::NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 
 	_assert_(gauss->Enum()==GaussTriaEnum);
-	this->GetNodalFunctions(basis,(GaussTria*)gauss,P2Enum);
+	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss,this->VelocityInterpolation());
 
 }
 /*}}}*/
-void       Tria::NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void       Tria::NodalFunctionsPressure(IssmDouble* basis, Gauss* gauss){/*{{{*/
 
 	_assert_(gauss->Enum()==GaussTriaEnum);
-	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss,this->element_type);
+	this->GetNodalFunctions(basis,(GaussTria*)gauss,this->PressureInterpolation());
 
 }
 /*}}}*/
-void       Tria::NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void       Tria::NodalFunctionsP1(IssmDouble* basis, Gauss* gauss){/*{{{*/
 
 	_assert_(gauss->Enum()==GaussTriaEnum);
-	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss,P1Enum);
+	this->GetNodalFunctions(basis,(GaussTria*)gauss,P1Enum);
 
 }
 /*}}}*/
-void       Tria::NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+void       Tria::NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
 
 	_assert_(gauss->Enum()==GaussTriaEnum);
-	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss,this->VelocityInterpolation());
+	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss,P1Enum);
 
 }
 /*}}}*/
-void       Tria::NodalFunctionsVelocity(IssmDouble* basis, Gauss* gauss){/*{{{*/
+void       Tria::NodalFunctionsP2(IssmDouble* basis, Gauss* gauss){/*{{{*/
 
 	_assert_(gauss->Enum()==GaussTriaEnum);
-	this->GetNodalFunctions(basis,(GaussTria*)gauss,this->VelocityInterpolation());
+	this->GetNodalFunctions(basis,(GaussTria*)gauss,P2Enum);
 
 }
 /*}}}*/
-void       Tria::NodalFunctionsPressure(IssmDouble* basis, Gauss* gauss){/*{{{*/
+void       Tria::NodalFunctionsTensor(IssmDouble* basis, Gauss* gauss){/*{{{*/
 
 	_assert_(gauss->Enum()==GaussTriaEnum);
-	this->GetNodalFunctions(basis,(GaussTria*)gauss,this->PressureInterpolation());
+	this->GetNodalFunctions(basis,(GaussTria*)gauss,this->TensorInterpolation());
 
 }
 /*}}}*/
-void       Tria::NodalFunctionsTensor(IssmDouble* basis, Gauss* gauss){/*{{{*/
+void       Tria::NodalFunctionsVelocity(IssmDouble* basis, Gauss* gauss){/*{{{*/
 
 	_assert_(gauss->Enum()==GaussTriaEnum);
-	this->GetNodalFunctions(basis,(GaussTria*)gauss,this->TensorInterpolation());
+	this->GetNodalFunctions(basis,(GaussTria*)gauss,this->VelocityInterpolation());
 
 }
 /*}}}*/
@@ -1794,6 +2158,21 @@
 	_assert_(bed_normal[1]<0); 
 }
 /*}}}*/
+void       Tria::NormalSection(IssmDouble* normal,IssmDouble* xyz_list){/*{{{*/
+
+	/*Build unit outward pointing vector*/
+	IssmDouble vector[2];
+	IssmDouble norm;
+
+	vector[0]=xyz_list[1*3+0] - xyz_list[0*3+0];
+	vector[1]=xyz_list[1*3+1] - xyz_list[0*3+1];
+
+	norm=sqrt(vector[0]*vector[0] + vector[1]*vector[1]);
+
+	normal[0]= + vector[1]/norm;
+	normal[1]= - vector[0]/norm;
+}
+/*}}}*/
 void       Tria::NormalTop(IssmDouble* top_normal,IssmDouble* xyz_list){/*{{{*/
 
 	/*Build unit outward pointing vector*/
@@ -1812,18 +2191,16 @@
 	_assert_(top_normal[1]>0); 
 }
 /*}}}*/
-int        Tria::VelocityInterpolation(void){/*{{{*/
-	return TriaRef::VelocityInterpolation(this->element_type);
+int        Tria::ObjectEnum(void){/*{{{*/
+
+	return TriaEnum;
+
 }
 /*}}}*/
 int        Tria::PressureInterpolation(void){/*{{{*/
 	return TriaRef::PressureInterpolation(this->element_type);
 }
 /*}}}*/
-int        Tria::TensorInterpolation(void){/*{{{*/
-	return TriaRef::TensorInterpolation(this->element_type);
-}
-/*}}}*/
 int        Tria::NumberofNodesPressure(void){/*{{{*/
 	return TriaRef::NumberofNodes(this->PressureInterpolation());
 }
@@ -1900,6 +2277,33 @@
 	delete gauss;
 }
 /*}}}*/
+void       Tria::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 grounded 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);
+			}
+		}
+	}
+}
+/*}}}*/
 void       Tria::ReduceMatrices(ElementMatrix* Ke,ElementVector* pe){/*{{{*/
 
 	/*Static condensation if requested*/
@@ -2010,6 +2414,82 @@
 	_error_("not implemented yet");
 }
 /*}}}*/
+void       Tria::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){/*{{{*/
+
+	IssmDouble  values[NUMVERTICES];
+	int         vertexpidlist[NUMVERTICES],control_init;
+
+
+	/*Get Domain type*/
+	int domaintype;
+	parameters->FindParam(&domaintype,DomainTypeEnum);
+
+	/*Specific case for depth averaged quantities*/
+	control_init=control_enum;
+	if(domaintype==Domain2DverticalEnum){
+		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 TriaInput(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);
+}
+/*}}}*/
+void       Tria::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*/
+	if(this->element_type_list) this->element_type=this->element_type_list[analysis_counter];
+
+	/*Pick up nodes*/
+	if(this->hnodes && this->hnodes[analysis_counter]){
+		this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp();
+	}
+
+}
+/*}}}*/
+Element*   Tria::SpawnBasalElement(void){/*{{{*/
+
+	int index1,index2;
+	int domaintype;
+
+	this->parameters->FindParam(&domaintype,DomainTypeEnum);
+	switch(domaintype){
+		case Domain2DhorizontalEnum:
+			return this;
+		case Domain2DverticalEnum:
+			_assert_(HasEdgeOnBase());
+			this->EdgeOnBaseIndices(&index1,&index2);
+			return SpawnSeg(index1,index2);
+		default:
+			_error_("not implemented yet");
+	}
+}
+/*}}}*/
 Seg*       Tria::SpawnSeg(int index1,int index2){/*{{{*/
 
 	int analysis_counter;
@@ -2037,7 +2517,7 @@
 	return seg;
 }
 /*}}}*/
-Element*   Tria::SpawnBasalElement(void){/*{{{*/
+Element*   Tria::SpawnTopElement(void){/*{{{*/
 
 	int index1,index2;
 	int domaintype;
@@ -2047,46 +2527,104 @@
 		case Domain2DhorizontalEnum:
 			return this;
 		case Domain2DverticalEnum:
-			_assert_(HasEdgeOnBase());
-			this->EdgeOnBaseIndices(&index1,&index2);
-			return SpawnSeg(index1,index2);
+			_assert_(HasEdgeOnSurface());
+			this->EdgeOnSurfaceIndices(&index1,&index2);
+			return SpawnSeg(index2,index1); //reverse order
 		default:
 			_error_("not implemented yet");
 	}
 }
 /*}}}*/
-Element*   Tria::SpawnTopElement(void){/*{{{*/
+void       Tria::StrainRateparallel(){/*{{{*/
 
-	int index1,index2;
-	int domaintype;
+	IssmDouble *xyz_list = NULL;
+	IssmDouble  epsilon[3];
+	GaussTria* gauss=NULL;
+	IssmDouble  vx,vy,vel;
+	IssmDouble  strainxx;
+	IssmDouble  strainxy;
+	IssmDouble  strainyy;
+	IssmDouble  strainparallel[NUMVERTICES];
 
-	this->parameters->FindParam(&domaintype,DomainTypeEnum);
-	switch(domaintype){
-		case Domain2DhorizontalEnum:
-			return this;
-		case Domain2DverticalEnum:
-			_assert_(HasEdgeOnSurface());
-			this->EdgeOnSurfaceIndices(&index1,&index2);
-			return SpawnSeg(index2,index1); //reverse order
-		default:
-			_error_("not implemented yet");
+	/* 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);
+
+	/* Start looping on the number of vertices: */
+	gauss=new GaussTria();
+	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 viscosity and pressure: */
+		this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
+		strainxx=epsilon[0];
+		strainyy=epsilon[1];
+		strainxy=epsilon[2];
+
+		/*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 TriaInput(StrainRateparallelEnum,&strainparallel[0],P1Enum));
+
+	/*Clean up and return*/
+	delete gauss;
+	xDelete<IssmDouble>(xyz_list);
 }
 /*}}}*/
-void       Tria::SetCurrentConfiguration(Elements* elementsin, Loads* loadsin, Nodes* nodesin, Materials* materialsin, Parameters* parametersin){/*{{{*/
+void       Tria::StrainRateperpendicular(){/*{{{*/
 
-	/*go into parameters and get the analysis_counter: */
-	int analysis_counter;
-	parametersin->FindParam(&analysis_counter,AnalysisCounterEnum);
+	IssmDouble *xyz_list = NULL;
+	GaussTria* gauss=NULL;
+	IssmDouble  epsilon[3];
+	IssmDouble  vx,vy,vel;
+	IssmDouble  strainxx;
+	IssmDouble  strainxy;
+	IssmDouble  strainyy;
+	IssmDouble  strainperpendicular[NUMVERTICES];
 
-	/*Get Element type*/
-	if(this->element_type_list) this->element_type=this->element_type_list[analysis_counter];
+	/* Get node coordinates and dof list: */
+	this->GetVerticesCoordinates(&xyz_list);
 
-	/*Pick up nodes*/
-	if(this->hnodes && this->hnodes[analysis_counter]){
-		this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp();
+	/*Retrieve all inputs we will need*/
+	Input* vx_input=inputs->GetInput(VxEnum);                                  _assert_(vx_input);
+	Input* vy_input=inputs->GetInput(VyEnum);                                  _assert_(vy_input);
+
+	/* Start looping on the number of vertices: */
+	gauss=new GaussTria();
+	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 viscosity and pressure: */
+		this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
+		strainxx=epsilon[0];
+		strainyy=epsilon[1];
+		strainxy=epsilon[2];
+
+		/*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 TriaInput(StrainRateperpendicularEnum,&strainperpendicular[0],P1Enum));
+
+	/*Clean up and return*/
+	delete gauss;
+	xDelete<IssmDouble>(xyz_list);
 }
 /*}}}*/
 IssmDouble Tria::SurfaceArea(void){/*{{{*/
@@ -2116,6 +2654,10 @@
 	return S;
 }
 /*}}}*/
+int        Tria::TensorInterpolation(void){/*{{{*/
+	return TriaRef::TensorInterpolation(this->element_type);
+}
+/*}}}*/
 IssmDouble Tria::TimeAdapt(void){/*{{{*/
 
 	/*intermediary: */
@@ -2157,6 +2699,34 @@
 	return dt;
 }
 /*}}}*/
+IssmDouble Tria::TotalSmb(void){/*{{{*/
+
+	/*The smb[kg yr-1] of one element is area[m2] * smb [kg m^-2 yr^-1]*/
+	IssmDouble base,smb,rho_ice;
+	IssmDouble Total_Smb=0;
+	IssmDouble xyz_list[NUMVERTICES][3];
+
+	/*Get material parameters :*/
+	rho_ice=matpar->GetRhoIce();
+
+   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/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]));	// area of element in m2
+
+	/*Now get the average SMB over the element*/
+	Input* smb_input = inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(smb_input);
+	smb_input->GetInputAverage(&smb);																								// average smb on element in m ice s-1
+   Total_Smb=rho_ice*base*smb;																											// smb on element in kg s-1
+
+	/*Return: */
+	return Total_Smb;
+}
+/*}}}*/
 void       Tria::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finiteelement_type){/*{{{*/
 
 	/*Intermediaries*/
@@ -2346,487 +2916,134 @@
 
 }
 /*}}}*/
-void       Tria::ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss){/*{{{*/
-	TriaRef::GetInputValue(pvalue,values,gauss,P1Enum);
-}
-/*}}}*/
-void       Tria::ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
-	TriaRef::GetInputDerivativeValue(dvalue,values,xyz_list,gauss,P1Enum);
-}
-/*}}}*/
-int        Tria::VertexConnectivity(int vertexindex){/*{{{*/
-	_assert_(this->vertices);
-	return this->vertices[vertexindex]->Connectivity();
-}
-/*}}}*/
-bool       Tria::IsZeroLevelset(int levelset_enum){/*{{{*/
+int        Tria::UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf){/*{{{*/
 
-	bool iszerols;
-	IssmDouble ls[NUMVERTICES];
+	int i;
+	int nflipped=0;
 
-	/*Retrieve all inputs and parameters*/
-	GetInputListOnVertices(&ls[0],levelset_enum);
+	/*Go through nodes, and whoever is on the potential_ungrounding, ends up in nodes_on_iceshelf: */
+	for(i=0;i<3;i++){
+		if (reCast<bool>(vertices_potentially_ungrounding[vertices[i]->Pid()])){
+			vec_nodes_on_iceshelf->SetValue(vertices[i]->Pid(),-1.,INS_VAL);
 
-	/*If the level set is awlays <0, 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;
+			/*If node was not on ice shelf, we flipped*/
+			if(nodes_on_iceshelf[vertices[i]->Pid()]>=0.){
+				nflipped++;
+			}
 		}
 	}
-
-	return iszerols;
+	return nflipped;
 }
 /*}}}*/
-bool       Tria::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<NUMVERTICES;i++)
-			if(ls[i]<0.) nrice++;
-		if(nrice==1) isicefront= true;
-	}
-	return isicefront;
-}/*}}}*/
-bool       Tria::IsFaceOnBoundary(void){/*{{{*/
-
-	IssmDouble values[NUMVERTICES];
-	IssmDouble sum;
-
-	/*Retrieve all inputs and parameters*/
-	GetInputListOnVertices(&values[0],MeshVertexonboundaryEnum);
-	sum = values[0]+values[1]+values[2];
-
-	_assert_(sum==0. || sum==1. || sum==2.);
-
-	if(sum==3.)  _error_("Two edges on boundary not supported yet...");
-
-	if(sum>1.){
-		return true;
-	}
-	else{
-		return false;
-	}
-}/*}}}*/
-void       Tria::AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part){/*{{{*/
-
-	bool       already = false;
-	int        i,j;
-	int        partition[NUMVERTICES];
-	int        offsetsid[NUMVERTICES];
-	int        offsetdof[NUMVERTICES];
-	IssmDouble area;
-	IssmDouble mean;
-
-	/*First, get the area: */
-	area=this->GetArea();
-
-	/*Figure out the average for this element: */
-	this->GetVerticesSidList(&offsetsid[0]);
-	this->GetVertexPidList(&offsetdof[0]);
-	mean=0;
-	for(i=0;i<NUMVERTICES;i++){
-		partition[i]=reCast<int>(qmu_part[offsetsid[i]]);
-		mean=mean+1.0/NUMVERTICES*vertex_response[offsetdof[i]];
-	}
-
-	/*Add contribution: */
-	for(i=0;i<NUMVERTICES;i++){
-		already=false;
-		for(j=0;j<i;j++){
-			if (partition[i]==partition[j]){
-				already=true;
-				break;
-			}
-		}
-		if(!already){
-			partition_contributions->SetValue(partition[i],mean*area,ADD_VAL);
-			partition_areas->SetValue(partition[i],area,ADD_VAL);
-		};
-	}
+void       Tria::ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
+	TriaRef::GetInputDerivativeValue(dvalue,values,xyz_list,gauss,P1Enum);
 }
 /*}}}*/
-IssmDouble Tria::IceVolume(void){/*{{{*/
-
-	/*The volume of a troncated prism is base * 1/3 sum(length of edges)*/
-	IssmDouble base,surface,bed;
-	IssmDouble xyz_list[NUMVERTICES][3];
-
-	if(!IsIceInElement())return 0;
-
-	/*First get back the area of the base*/
-	base=this->GetArea();
-
-	/*Now get the average height*/
-	Input* surface_input = inputs->GetInput(SurfaceEnum); _assert_(surface_input);
-	Input* base_input     = inputs->GetInput(BaseEnum);     _assert_(base_input);
-	surface_input->GetInputAverage(&surface);
-	base_input->GetInputAverage(&bed);
-
-	/*Return: */
-	int domaintype;
-	parameters->FindParam(&domaintype,DomainTypeEnum);
-	if(domaintype==Domain2DverticalEnum){
-	  return base;
-	}
-	else{
-	  return base*(surface-bed);
-	}
+void       Tria::ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss){/*{{{*/
+	TriaRef::GetInputValue(pvalue,values,gauss,P1Enum);
 }
 /*}}}*/
-IssmDouble Tria::IceVolumeAboveFloatation(void){/*{{{*/
-
-	/*The volume above floatation: H + rho_water/rho_ice * bathymetry */
-	IssmDouble rho_ice,rho_water;
-	IssmDouble base,surface,bed,bathymetry;
-	IssmDouble xyz_list[NUMVERTICES][3];
-
-	if(!IsIceInElement() || IsFloating())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/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 height and bathymetry*/
-	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.));
+int        Tria::VelocityInterpolation(void){/*{{{*/
+	return TriaRef::VelocityInterpolation(this->element_type);
 }
 /*}}}*/
-IssmDouble Tria::MassFlux( IssmDouble x1, IssmDouble y1, IssmDouble x2, IssmDouble y2,int segment_id){/*{{{*/
-
-	int        domaintype;
-	IssmDouble mass_flux=0.;
-	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble normal[2];
-	IssmDouble length,rho_ice;
-	IssmDouble h1,h2;
-	IssmDouble vx1,vx2,vy1,vy2;
-	GaussTria* gauss_1=NULL;
-	GaussTria* gauss_2=NULL;
-
-	/*Get material parameters :*/
-	rho_ice=matpar->GetRhoIce();
-
-	/*First off, check that this segment belongs to this element: */
-	if (segment_id!=this->id)_error_("error message: segment with id " << segment_id << " does not belong to element with id:" << this->id);
-
-	/*Get xyz list: */
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-
-	/*get area coordinates of 0 and 1 locations: */
-	gauss_1=new GaussTria();
-	gauss_1->GaussFromCoords(x1,y1,&xyz_list[0][0]);
-	gauss_2=new GaussTria();
-	gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]);
-
-	normal[0]=cos(atan2(x1-x2,y2-y1));
-	normal[1]=sin(atan2(x1-x2,y2-y1));
-
-	length=sqrt(pow(x2-x1,2)+pow(y2-y1,2));
-
-	Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
-	this->parameters->FindParam(&domaintype,DomainTypeEnum);
-	Input* vx_input=NULL;
-	Input* vy_input=NULL;
-	if(domaintype==Domain2DhorizontalEnum){
-		vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
-		vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
-	}
-	else{
-		vx_input=inputs->GetInput(VxAverageEnum); _assert_(vx_input);
-		vy_input=inputs->GetInput(VyAverageEnum); _assert_(vy_input);
-	}
-
-	thickness_input->GetInputValue(&h1, gauss_1);
-	thickness_input->GetInputValue(&h2, gauss_2);
-	vx_input->GetInputValue(&vx1,gauss_1);
-	vx_input->GetInputValue(&vx2,gauss_2);
-	vy_input->GetInputValue(&vy1,gauss_1);
-	vy_input->GetInputValue(&vy2,gauss_2);
-
-	mass_flux= rho_ice*length*(  
-				(ONETHIRD*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*normal[0]+
-				(ONETHIRD*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*normal[1]
-				);
-
-	/*clean up and return:*/
-	delete gauss_1;
-	delete gauss_2;
-	return mass_flux;
+int        Tria::VertexConnectivity(int vertexindex){/*{{{*/
+	_assert_(this->vertices);
+	return this->vertices[vertexindex]->Connectivity();
 }
 /*}}}*/
-IssmDouble Tria::MassFlux( IssmDouble* segment){/*{{{*/
+void       Tria::ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum){/*{{{*/
 
-	int        domaintype;
-	IssmDouble mass_flux=0.;
-	IssmDouble xyz_list[NUMVERTICES][3];
-	IssmDouble normal[2];
-	IssmDouble length,rho_ice;
-	IssmDouble x1,y1,x2,y2,h1,h2;
-	IssmDouble vx1,vx2,vy1,vy2;
-	GaussTria* gauss_1=NULL;
-	GaussTria* gauss_2=NULL;
+	int         normal_orientation=0;
+	IssmDouble  s1,s2;
+	IssmDouble  levelset[NUMVERTICES];
 
-	/*Get material parameters :*/
-	rho_ice=matpar->GetRhoIce();
+	/*Recover parameters and values*/
+	IssmDouble* xyz_zero = xNew<IssmDouble>(2*3);
+	GetInputListOnVertices(&levelset[0],levelsetenum);
 
-	/*First off, check that this segment belongs to this element: */
-	if (reCast<int>(*(segment+4))!=this->id)_error_("error message: segment with id " << reCast<int>(*(segment+4)) << " does not belong to element with id:" << this->id);
+	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]);
 
-	/*Recover segment node locations: */
-	x1=*(segment+0); y1=*(segment+1); x2=*(segment+2); y2=*(segment+3);
+		if(levelset[2]<0.) normal_orientation=1; //orientation of quadrangle 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 xyz list: */
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-
-	/*get area coordinates of 0 and 1 locations: */
-	gauss_1=new GaussTria();
-	gauss_1->GaussFromCoords(x1,y1,&xyz_list[0][0]);
-	gauss_2=new GaussTria();
-	gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]);
-
-	normal[0]=cos(atan2(x1-x2,y2-y1));
-	normal[1]=sin(atan2(x1-x2,y2-y1));
-
-	length=sqrt(pow(x2-x1,2)+pow(y2-y1,2));
-
-	Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
-	this->parameters->FindParam(&domaintype,DomainTypeEnum);
-	Input* vx_input=NULL;
-	Input* vy_input=NULL;
-	if(domaintype==Domain2DhorizontalEnum){
-		vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
-		vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
+		/*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]);
 	}
-	else{
-		vx_input=inputs->GetInput(VxAverageEnum); _assert_(vx_input);
-		vy_input=inputs->GetInput(VyAverageEnum); _assert_(vy_input);
-	}
+	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]);
 
-	thickness_input->GetInputValue(&h1, gauss_1);
-	thickness_input->GetInputValue(&h2, gauss_2);
-	vx_input->GetInputValue(&vx1,gauss_1);
-	vx_input->GetInputValue(&vx2,gauss_2);
-	vy_input->GetInputValue(&vy1,gauss_1);
-	vy_input->GetInputValue(&vy2,gauss_2);
+		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]);
 
-	mass_flux= rho_ice*length*(  
-				(ONETHIRD*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*normal[0]+
-				(ONETHIRD*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*normal[1]
-				);
-
-	/*clean up and return:*/
-	delete gauss_1;
-	delete gauss_2;
-	return mass_flux;
-}
-/*}}}*/
-void       Tria::ElementResponse(IssmDouble* presponse,int response_enum){/*{{{*/
-
-	switch(response_enum){
-		case MaterialsRheologyBbarEnum:
-			*presponse=this->material->GetBbar();
-			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 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]);
 	}
+	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]);
 
-}
-/*}}}*/
-IssmDouble Tria::TotalSmb(void){/*{{{*/
+		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]);
 
-	/*The smb[kg yr-1] of one element is area[m2] * smb [kg m^-2 yr^-1]*/
-	IssmDouble base,smb,rho_ice;
-	IssmDouble Total_Smb=0;
-	IssmDouble xyz_list[NUMVERTICES][3];
-
-	/*Get material parameters :*/
-	rho_ice=matpar->GetRhoIce();
-
-   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/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]));	// area of element in m2
-
-	/*Now get the average SMB over the element*/
-	Input* smb_input = inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(smb_input);
-	smb_input->GetInputAverage(&smb);																								// average smb on element in m ice s-1
-   Total_Smb=rho_ice*base*smb;																											// smb on element in kg s-1
-
-	/*Return: */
-	return Total_Smb;
-}
-/*}}}*/
-IssmDouble Tria::MisfitArea(int weightsenum){/*{{{*/
-
-	/*Intermediaries*/
-	IssmDouble weight;
-	IssmDouble Jdet;
-	IssmDouble Jelem = 0;
-	IssmDouble xyz_list[NUMVERTICES][3];
-	GaussTria *gauss = NULL;
-
-	/*If on water, return 0: */
-	if(!IsIceInElement())return 0;
-
-	/*Retrieve all inputs we will be needing: */
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	Input* weights_input     =inputs->GetInput(weightsenum);     _assert_(weights_input);
-
-	/* Start  looping on the number of gaussian points: */
-	gauss=new GaussTria(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
-
-		/* Get Jacobian determinant: */
-		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
-
-		/*Get parameters at gauss point*/
-		weights_input->GetInputValue(&weight,gauss);
-
-		/*compute misfit between model and observation */
-		Jelem+=Jdet*weight*gauss->weight;
+		/*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]);
 	}
+	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];
 
-	/* clean up and Return: */
-	delete gauss;
-	return Jelem;
-}
-/*}}}*/
-IssmDouble Tria::Misfit(int modelenum,int observationenum,int weightsenum){/*{{{*/
-
-	/*Intermediaries*/
-	IssmDouble model,observation,weight;
-	IssmDouble Jdet;
-	IssmDouble Jelem = 0;
-	IssmDouble xyz_list[NUMVERTICES][3];
-	GaussTria *gauss = NULL;
-
-	/*If on water, return 0: */
-	if(!IsIceInElement())return 0;
-
-	/*Retrieve all inputs we will be needing: */
-	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
-	Input* model_input=inputs->GetInput(modelenum);   _assert_(model_input);
-	Input* observation_input=inputs->GetInput(observationenum);_assert_(observation_input);
-	Input* weights_input     =inputs->GetInput(weightsenum);     _assert_(weights_input);
-
-	/* Start  looping on the number of gaussian points: */
-	gauss=new GaussTria(2);
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-
-		gauss->GaussPoint(ig);
-
-		/* Get Jacobian determinant: */
-		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss);
-
-		/*Get parameters at gauss point*/
-		model_input->GetInputValue(&model,gauss);
-		observation_input->GetInputValue(&observation,gauss);
-		weights_input->GetInputValue(&weight,gauss);
-
-		/*compute misfit between model and observation */
-		Jelem+=0.5*(model-observation)*(model-observation)*Jdet*weight*gauss->weight;
+		/*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];
 	}
+	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];
 
-	/* clean up and Return: */
-	delete gauss;
-	return Jelem;
-}
-/*}}}*/
-IssmDouble Tria::Masscon(IssmDouble* levelset){ /*{{{*/
-
-
-	/*intermediary: */
-	IssmDouble* values=NULL;
-	Input*      thickness_input=NULL;
-	IssmDouble  thickness;
-	IssmDouble  weight;
-	IssmDouble  Jdet;
-	IssmDouble  volume;
-	IssmDouble  rho_ice;
-	IssmDouble* xyz_list=NULL;
-	int         point1;
-	IssmDouble  fraction1,fraction2;
-	bool        mainlynegative=true;
-	
-	/*Output:*/
-	volume=0;
-
-	/* Get node coordinates and dof list: */
-	GetVerticesCoordinates(&xyz_list);
-
-	/*Retrieve inputs required:*/
-	thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input);
-	
-	/*Retrieve material parameters: */
-	rho_ice=matpar->GetRhoIce();
-
-	/*Retrieve values of the levelset defining the masscon: */
-	values = xNew<IssmDouble>(NUMVERTICES);
-	for(int i=0;i<NUMVERTICES;i++){
-		values[i]=levelset[this->vertices[i]->Sid()];
+		/*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];
 	}
-		
-	/*Ok, use the level set values to figure out where we put our gaussian points:*/
-	this->GetLevelsetPositivePart(&point1,&fraction1,&fraction2,&mainlynegative,values);
-	Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,mainlynegative,4);
+	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];
 
-	volume=0;
-
-	for(int ig=gauss->begin();ig<gauss->end();ig++){
-		gauss->GaussPoint(ig);
-
-		this->JacobianDeterminant(&Jdet,xyz_list,gauss);
-		thickness_input->GetInputValue(&thickness, gauss);
-
-		volume+=thickness*gauss->weight*Jdet;
+		/*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];
 	}
+	else _error_("Case not covered");
 
-	/* clean up and Return: */
-	xDelete<IssmDouble>(xyz_list);
-	xDelete<IssmDouble>(values);
-	delete gauss;
-	return rho_ice*volume;
+	/*Assign output pointer*/
+	*pxyz_zero= xyz_zero;
 }
 /*}}}*/
 
@@ -2958,179 +3175,45 @@
 /*}}}*/
 #endif
 
-void       Tria::InputControlUpdate(IssmDouble scalar,bool save_parameter){/*{{{*/
+#ifdef _HAVE_DAKOTA_
+void       Tria::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){/*{{{*/
 
-	/*Intermediary*/
-	int    num_controls;
-	int*   control_type=NULL;
-	Input* input=NULL;
+	int             i,t,row;
+	IssmDouble      time;
+	TransientInput *transientinput = NULL;
+	IssmDouble      values[3];
 
-	/*retrieve some parameters: */
-	this->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
-	this->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
+	/*Check that name is an element input*/
+	if (!IsInput(name)) return;
 
-	for(int i=0;i<num_controls;i++){
-		input=(Input*)this->inputs->GetInput(control_type[i]);   _assert_(input);
-		if (input->ObjectEnum()!=ControlInputEnum){
-			_error_("input " << EnumToStringx(control_type[i]) << " is not a ControlInput");
-		}
+	switch(type){
 
-		((ControlInput*)input)->UpdateValue(scalar);
-		((ControlInput*)input)->Constrain();
-		if (save_parameter) ((ControlInput*)input)->SaveValue();
+		case VertexEnum:
+			/*Create transient input: */
+			for(t=0;t<ncols;t++){ //ncols is the number of times
 
-	}
+				/*create input values: */
+				for(i=0;i<3;i++){
+					row=this->vertices[i]->Sid();
+					values[i]=matrix[ncols*row+t];
+				}
 
-	/*Clean up and return*/
-	xDelete<int>(control_type);
-}
-/*}}}*/
-void       Tria::ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){/*{{{*/
+				/*time:*/
+				time=matrix[(nrows-1)*ncols+t];
 
-	int    vertexpidlist[NUMVERTICES];
-	IssmDouble grad_list[NUMVERTICES];
-	Input* grad_input=NULL;
+				if(t==0) transientinput=new TransientInput(name);
+				transientinput->AddTimeInput(new TriaInput(name,values,P1Enum),time);
+				transientinput->Configure(parameters);
+			}
+			this->inputs->AddInput(transientinput);
+			break;
 
-	Input* 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 TriaInput(GradientEnum,grad_list,P1Enum);
-
-	((ControlInput*)input)->SetGradient(grad_input);
-
-}/*}}}*/
-void       Tria::ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum){/*{{{*/
-
-	Input* 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]);
-
-	GaussTria* gauss=new GaussTria();
-	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]);
+		default:
+			_error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
 	}
-	delete gauss;
 
-	vector_control->SetValues(NUMVERTICES,&sidlist[0],&values[0],ADD_VAL);
-	vector_gradient->SetValues(NUMVERTICES,&sidlist[0],&gradients[0],ADD_VAL);
-
-}/*}}}*/
-void       Tria::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,const char* data,bool onsid){/*{{{*/
-
-	int vertexidlist[NUMVERTICES];
-	Input *input=NULL;
-
-	/*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)*/
-	input=(Input*)this->inputs->GetInput(control_enum);   _assert_(input);
-
-	/*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       Tria::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){/*{{{*/
-
-	IssmDouble  values[NUMVERTICES];
-	int         vertexpidlist[NUMVERTICES],control_init;
-
-
-	/*Get Domain type*/
-	int domaintype;
-	parameters->FindParam(&domaintype,DomainTypeEnum);
-
-	/*Specific case for depth averaged quantities*/
-	control_init=control_enum;
-	if(domaintype==Domain2DverticalEnum){
-		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 TriaInput(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);
-}
-/*}}}*/
-void       Tria::GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution, int enum_type){/*{{{*/
-
-	int        *doflist = NULL;
-	IssmDouble  value;
-
-	/*Fetch number of nodes for this finite element*/
-	int numnodes = this->NumberofNodes(this->element_type);
-
-	/*Fetch dof list and allocate solution vector*/
-	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
-	IssmDouble* values = xNew<IssmDouble>(numnodes);
-
-	/*Get inputs*/
-	Input* enum_input=inputs->GetInput(enum_type); _assert_(enum_input);
-
-	/*Ok, we have the values, fill in the array: */
-	GaussTria* gauss=new GaussTria();
-	for(int i=0;i<numnodes;i++){
-		gauss->GaussNode(this->element_type,i);
-
-		enum_input->GetInputValue(&value,gauss);
-		values[i]=value;
-	}
-
-	solution->SetValues(numnodes,doflist,values,INS_VAL);
-
-	/*Free ressources:*/
-	xDelete<int>(doflist);
-	xDelete<IssmDouble>(values);
-	delete gauss;
-}
-/*}}}*/
-
-#ifdef _HAVE_DAKOTA_
 void       Tria::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){/*{{{*/
 
 	int i,j;
@@ -3222,89 +3305,4 @@
 
 }
 /*}}}*/
-void       Tria::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){/*{{{*/
-
-	int             i,t,row;
-	IssmDouble      time;
-	TransientInput *transientinput = NULL;
-	IssmDouble      values[3];
-
-	/*Check that name is an element input*/
-	if (!IsInput(name)) return;
-
-	switch(type){
-
-		case VertexEnum:
-			/*Create transient input: */
-			for(t=0;t<ncols;t++){ //ncols is the number of times
-
-				/*create input values: */
-				for(i=0;i<3;i++){
-					row=this->vertices[i]->Sid();
-					values[i]=matrix[ncols*row+t];
-				}
-
-				/*time:*/
-				time=matrix[(nrows-1)*ncols+t];
-
-				if(t==0) transientinput=new TransientInput(name);
-				transientinput->AddTimeInput(new TriaInput(name,values,P1Enum),time);
-				transientinput->Configure(parameters);
-			}
-			this->inputs->AddInput(transientinput);
-			break;
-
-		default:
-			_error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
-	}
-
-}
-/*}}}*/
 #endif
-
-void       Tria::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 grounded 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        Tria::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<3;i++){
-		if (reCast<bool>(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;
-}
-/*}}}*/
Index: ../trunk-jpl/src/c/classes/Elements/Tria.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Tria.h	(revision 18910)
+++ ../trunk-jpl/src/c/classes/Elements/Tria.h	(revision 18911)
@@ -41,105 +41,101 @@
 		Object *copy();
 		/*}}}*/
 		/*Update virtual functions resolution: {{{*/
-		void  InputUpdateFromVector(IssmDouble* vector, int name, int type);
 		#ifdef _HAVE_DAKOTA_
+		void  InputUpdateFromMatrixDakota(IssmDouble* matrix, int nows, int ncols, int name, int type);
 		void  InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type);
-		void  InputUpdateFromMatrixDakota(IssmDouble* matrix, int nows, int ncols, int name, int type);
 		#endif
 		void  InputUpdateFromIoModel(int index, IoModel* iomodel);
+		void  InputUpdateFromVector(IssmDouble* vector, int name, int type);
 		/*}}}*/
 		/*Element virtual functions definitions: {{{*/
+		void        AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part);
+		void			CalvingRateLevermann();
 		IssmDouble  CharacteristicLength(void);
 		void        ComputeBasalStress(Vector<IssmDouble>* sigma_b);
+		void        ComputeDeviatoricStressTensor();
 		void        ComputeSigmaNN();
 		void        ComputeStressTensor();
-		void        ComputeDeviatoricStressTensor();
-		void			StrainRateparallel();
-		void			StrainRateperpendicular();
 		void        ComputeSurfaceNormalVelocity();
-		void        StressIntensityFactor(void){_error_("not implemented yet");};
-		void			CalvingRateLevermann();
 		void        Configure(Elements* elements,Loads* loads,Nodes* nodesin,Vertices* verticesin,Materials* materials,Parameters* parameters);
-		void        SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);
-		void        ResetHooks();
+		void        ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index);
+		void        ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum);
 		void        Delta18oParameterization(void);
+		int         EdgeOnBaseIndex();
+		void        EdgeOnBaseIndices(int* pindex1,int* pindex);
+		int         EdgeOnSurfaceIndex();
+		void        EdgeOnSurfaceIndices(int* pindex1,int* pindex);
+		void        ElementResponse(IssmDouble* presponse,int response_enum);
 		void        ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);
+		int         FiniteElement(void);
 		void        FSContactMigration(Vector<IssmDouble>* vertexgrounded,Vector<IssmDouble>* vertexfloating);
-		int         FiniteElement(void);
-		Element*    GetUpperElement(void){_error_("not implemented yet");};
 		Element*    GetBasalElement(void){_error_("not implemented yet");};
 		void        GetLevelsetPositivePart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlynegative,IssmDouble* levelsetvalues);
 		void        GetGroundedPart(int* point1,IssmDouble* fraction1, IssmDouble* fraction2,bool* mainlyfloating);
 		IssmDouble  GetGroundedPortion(IssmDouble* xyz_list);
+		void	      GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum);
+		void	      GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level);
 		int         GetNodeIndex(Node* node);
 		int         GetNumberOfNodes(void);
 		int         GetNumberOfNodes(int enum_type);
 		int         GetNumberOfVertices(void);
-		bool        IsOnBase();
-		bool        IsOnSurface();
-		bool        HasEdgeOnBase();
-		bool        HasEdgeOnSurface();
-		void        EdgeOnSurfaceIndices(int* pindex1,int* pindex);
-		void        EdgeOnBaseIndices(int* pindex1,int* pindex);
-		int         EdgeOnBaseIndex();
-		int         EdgeOnSurfaceIndex();
-		bool        IsNodeOnShelfFromFlags(IssmDouble* flags);
-		int         NumberofNodesVelocity(void);
-		int         NumberofNodesPressure(void);
 		void        GetSolutionFromInputsOneDof(Vector<IssmDouble>* solution,int enum_type);
+		Element*    GetUpperElement(void){_error_("not implemented yet");};
+		void        GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data,bool onsid);
 		void        GetVerticesCoordinatesBase(IssmDouble** pxyz_list);
 		void        GetVerticesCoordinatesTop(IssmDouble** pxyz_list);
+		bool        HasEdgeOnBase();
+		bool        HasEdgeOnSurface();
+		IssmDouble  IceVolume(void);
+		IssmDouble  IceVolumeAboveFloatation(void);
+		void        InputControlUpdate(IssmDouble scalar,bool save_parameter);
 		void        InputDepthAverageAtBase(int enum_type,int average_enum_type);
 		void        InputExtrude(int enum_type,int start){_error_("not implemented"); /*For penta only*/};
 		void        InputScale(int enum_type,IssmDouble scale_factor);
+		bool	   	IsFaceOnBoundary(void);
+		bool	   	IsIcefront(void);
+		bool        IsNodeOnShelfFromFlags(IssmDouble* flags);
+		bool        IsOnBase();
+		bool        IsOnSurface();
+		bool        IsZeroLevelset(int levelset_enum);
+		IssmDouble  Masscon(IssmDouble* levelset);
+		IssmDouble  MassFlux(IssmDouble* segment);
+		IssmDouble  MassFlux(IssmDouble x1,IssmDouble y1, IssmDouble x2, IssmDouble y2,int segment_id);
 		void        MaterialUpdateFromTemperature(void){_error_("not implemented yet");};
+		IssmDouble  Misfit(int modelenum,int observationenum,int weightsenum);
+		IssmDouble  MisfitArea(int weightsenum);
 		int         NodalValue(IssmDouble* pvalue, int index, int natureofdataenum);
+		int         NumberofNodesPressure(void);
+		int         NumberofNodesVelocity(void);
 		void        PositiveDegreeDay(IssmDouble* pdds,IssmDouble* pds,IssmDouble signorm);
+		void        PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding);
+		int         PressureInterpolation();
 		void        ReduceMatrices(ElementMatrix* Ke,ElementVector* pe);
 		void        ResetFSBasalBoundaryCondition(void);
+		void        ResetHooks();
+		void        SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index);
+		void        SetCurrentConfiguration(Elements* elements,Loads* loads,Nodes* nodes,Materials* materials,Parameters* parameters);
 	   Element*    SpawnBasalElement(void);
 		Element*    SpawnTopElement(void);
-		int         VelocityInterpolation();
-		int         PressureInterpolation();
+		void			StrainRateparallel();
+		void			StrainRateperpendicular();
+		void        StressIntensityFactor(void){_error_("not implemented yet");};
+		IssmDouble  SurfaceArea(void);
 		int         TensorInterpolation();
-		IssmDouble  SurfaceArea(void);
+		IssmDouble  TimeAdapt();
+		IssmDouble  TotalSmb(void);
 		void        Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finitelement);
-		IssmDouble  TimeAdapt();
+		int         UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf);
+		void        ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss);
 		void        ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss);
-		void        ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss);
+		int         VelocityInterpolation();
 		int         VertexConnectivity(int vertexindex);
 		void        VerticalSegmentIndices(int** pindices,int* pnumseg){_error_("not implemented yet");};
 		void        ZeroLevelsetCoordinates(IssmDouble** pxyz_zero,IssmDouble* xyz_list,int levelsetenum);
-		void	      GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum);
-		void	      GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level);
-		bool        IsZeroLevelset(int levelset_enum);
-		bool	   	IsIcefront(void);
-		bool	   	IsFaceOnBoundary(void);
 
-		void       AverageOntoPartition(Vector<IssmDouble>* partition_contributions,Vector<IssmDouble>* partition_areas,IssmDouble* vertex_response,IssmDouble* qmu_part);
-		IssmDouble IceVolume(void);
-		IssmDouble IceVolumeAboveFloatation(void);
-		IssmDouble TotalSmb(void);
-		IssmDouble MassFlux(IssmDouble* segment);
-		IssmDouble MassFlux(IssmDouble x1,IssmDouble y1, IssmDouble x2, IssmDouble y2,int segment_id);
-		void       ElementResponse(IssmDouble* presponse,int response_enum);
-		IssmDouble Masscon(IssmDouble* levelset);
-		IssmDouble Misfit(int modelenum,int observationenum,int weightsenum);
-		IssmDouble MisfitArea(int weightsenum);
-
 		#ifdef _HAVE_GIA_
 		void   GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,IssmDouble* x,IssmDouble* y);
 		#endif
-
-		void       GetVectorFromControlInputs(Vector<IssmDouble>* gradient,int control_enum,int control_index,const char* data,bool onsid);
-		void       SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index);
-		void       ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index);
-		void       ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum);
-		void       InputControlUpdate(IssmDouble scalar,bool save_parameter);
-
-		void   PotentialUngrounding(Vector<IssmDouble>* potential_sheet_ungrounding);
-		int    UpdatePotentialUngrounding(IssmDouble* vertices_potentially_ungrounding,Vector<IssmDouble>* vec_nodes_on_iceshelf,IssmDouble* nodes_on_iceshelf);
-
 		/*}}}*/
 		/*Tria specific routines:{{{*/
 		void           AddBasalInput(int input_enum, IssmDouble* values, int interpolation_enum);
@@ -147,18 +143,15 @@
 		IssmDouble     GetArea(void);
 		void           GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints);
 		int            GetElementType(void);
-		void           NormalSection(IssmDouble* normal,IssmDouble* xyz_list);
-		void           NormalTop(IssmDouble* normal,IssmDouble* xyz_list);
-		void           NormalBase(IssmDouble* normal,IssmDouble* xyz_list);
 		void           GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);
 		void           GetMaterialInputValue(IssmDouble* pvalue,Node* node,int enumtype);
 		Node*          GetNode(int node_number);
 		void	         InputUpdateFromSolutionOneDof(IssmDouble* solution,int enum_type);
 		void	         InputUpdateFromSolutionOneDofCollapsed(IssmDouble* solution,int enum_type){_error_("not implemented yet");};
 		void           JacobianDeterminant(IssmDouble*  pJdet, IssmDouble* xyz_list,Gauss* gauss);
+		void           JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss);
 		void           JacobianDeterminantLine(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
 		void           JacobianDeterminantSurface(IssmDouble*  pJdet, IssmDouble* xyz_list,Gauss* gauss);
-		void           JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss);
 		void           JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss);
 		IssmDouble     MinEdgeLength(IssmDouble* xyz_list){_error_("not implemented yet");};
 		Gauss*         NewGauss(void);
@@ -170,23 +163,25 @@
 		Gauss*         NewGaussLine(int vertex1,int vertex2,int order){_error_("not implemented yet");};
 		Gauss*         NewGaussTop(int order);
 		void           NodalFunctions(IssmDouble* basis,Gauss* gauss);
-		void           NodalFunctionsP1(IssmDouble* basis,Gauss* gauss);
-		void           NodalFunctionsP2(IssmDouble* basis,Gauss* gauss);
 		void           NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
-		void           NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
+		void           NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
 		void           NodalFunctionsMINIDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){_error_("not implemented yet");};
-		void           NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
-		void           NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss);
 		void           NodalFunctionsPressure(IssmDouble* basis,Gauss* gauss);
+		void           NodalFunctionsP1(IssmDouble* basis,Gauss* gauss);
+		void           NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss);
+		void           NodalFunctionsP2(IssmDouble* basis,Gauss* gauss);
 		void           NodalFunctionsTensor(IssmDouble* basis,Gauss* gauss);
+		void           NodalFunctionsVelocity(IssmDouble* basis,Gauss* gauss);
+		void           NormalBase(IssmDouble* normal,IssmDouble* xyz_list);
+		void           NormalSection(IssmDouble* normal,IssmDouble* xyz_list);
+		void           NormalTop(IssmDouble* normal,IssmDouble* xyz_list);
 		void	         SetClone(int* minranks);
 		void           SetTemporaryElementType(int element_type_in){_error_("not implemented yet");};
 		Seg*	         SpawnSeg(int index1,int index2);
 		IssmDouble     StabilizationParameter(IssmDouble u, IssmDouble v, IssmDouble w, IssmDouble diameter, IssmDouble kappa){_error_("not implemented yet");};
+		void           UpdateConstraintsExtrudeFromBase(void);
+		void           UpdateConstraintsExtrudeFromTop(void);
 		void           ViscousHeating(IssmDouble* pphi,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not implemented yet");};
-
-		void UpdateConstraintsExtrudeFromBase(void);
-		void UpdateConstraintsExtrudeFromTop(void);
 		/*}}}*/
 
 };
