/*!\file Tria.cpp
 * \brief: implementation of the Tria object
 */
/*Headers:*/
/*{{{*/
#ifdef HAVE_CONFIG_H
	#include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#include <stdio.h>
#include <string.h>
#include <math.h>
#include "../classes.h"
#include "../Inputs2/TriaInput2.h"
#include "../Inputs2/PentaInput2.h"
#include "../Inputs2/ControlInput2.h"
#include "../Inputs2/DatasetInput2.h"
#include "../Inputs2/TransientInput2.h"
#include "../../shared/shared.h"
#ifdef _HAVE_GIAIVINS_
#include "../../modules/GiaDeflectionCorex/GiaDeflectionCorex.h"
#endif
/*}}}*/

/*Element macros*/
#define NUMVERTICES   3
#define NUMVERTICES1D 2

/*Constructors/destructor/copy*/
Tria::Tria(int tria_id,int tria_sid,int tria_lid,IoModel* iomodel,int nummodels)/*{{{*/
	:ElementHook(nummodels,tria_id,NUMVERTICES,iomodel){

		this->iscollapsed = 0;

		/*id: */
		this->id  = tria_id;
		this->sid = tria_sid;
		this->lid = tria_lid;

		/*this->parameters: we still can't point to it, it may not even exist. Configure will handle this.*/
		this->parameters = NULL;

		/*initialize pointers:*/
		this->nodes    = NULL;
		this->vertices = NULL;
		this->material = NULL;
		if(nummodels>0){
			this->element_type_list=xNew<int>(nummodels);
			for(int i=0;i<nummodels;i++) this->element_type_list[i] = 0;
		}
		else this->element_type_list = NULL;

		/*surface and base*/
		IssmDouble sum;
		this->isonsurface = false;
		this->isonbase    = false;
		switch(iomodel->domaintype){
			case Domain2DverticalEnum:
				_assert_(iomodel->Data("md.mesh.vertexonsurface"));
				_assert_(iomodel->Data("md.mesh.vertexonbase"));
				sum = 0.;
				for(int i=0;i<NUMVERTICES;i++) sum += iomodel->Data("md.mesh.vertexonsurface")[reCast<int>(iomodel->elements[(tria_id-1)*NUMVERTICES+i])-1];
				_assert_(sum>=0 && sum<3);
				if(sum>1.) this->isonsurface = true;
				sum = 0.;
				for(int i=0;i<NUMVERTICES;i++) sum += iomodel->Data("md.mesh.vertexonbase")[reCast<int>(iomodel->elements[(tria_id-1)*NUMVERTICES+i])-1];
				_assert_(sum>=0 && sum<3);
				if(sum>1.) this->isonbase = true;
				break;
			case Domain2DhorizontalEnum:
			case Domain3DsurfaceEnum:
				this->isonsurface = true;
				this->isonbase    = true;
				break;
			default: _error_("mesh "<<EnumToStringx(iomodel->domaintype)<<" not supported yet");
		}

}/*}}}*/
Tria::~Tria(){/*{{{*/
	this->parameters=NULL;
}
/*}}}*/
Object* Tria::copy() {/*{{{*/

	int i;
	Tria* tria=NULL;

	tria=new Tria();

	tria->iscollapsed=this->iscollapsed;

	//deal with TriaRef mother class
	int nanalyses = this->numanalyses;
	if(nanalyses > 0){
		tria->element_type_list=xNew<int>(nanalyses);
		for(i=0;i<nanalyses;i++){
			if (this->element_type_list[i]) tria->element_type_list[i]=this->element_type_list[i];
			else tria->element_type_list[i] = 0;
		}
	}
	else tria->element_type_list = NULL;
	tria->element_type=this->element_type;
	tria->numanalyses=nanalyses;

	//deal with ElementHook mother class
	if (this->hnodes){
		tria->hnodes=xNew<Hook*>(tria->numanalyses);
		for(i=0;i<tria->numanalyses;i++){
			if (this->hnodes[i]) tria->hnodes[i] = (Hook*)(this->hnodes[i]->copy());
			else tria->hnodes[i] = NULL;
		}
	}
	else tria->hnodes = NULL;

	tria->hvertices = (Hook*)this->hvertices->copy();
	tria->hmaterial = (Hook*)this->hmaterial->copy();
	tria->hneighbors = NULL;

	/*deal with Tria fields: */
	tria->id  = this->id;
	tria->sid = this->sid;
	tria->lid = this->lid;
	tria->isonbase  = this->isonbase;
	tria->isonsurface  = this->isonsurface;

	/*point parameters: */
	tria->parameters=this->parameters;

	/*recover objects: */
	if (this->nodes){
		unsigned int num_nodes = 3;
		tria->nodes = xNew<Node*>(num_nodes); //we cannot rely on an analysis_counter to tell us which analysis_type we are running, so we just copy the nodes.
		for(i=0;i<num_nodes;i++) if(this->nodes[i]) tria->nodes[i]=this->nodes[i]; else tria->nodes[i] = NULL;
	}
	else tria->nodes = NULL;

	tria->vertices = (Vertex**)this->hvertices->deliverp();
	tria->material = (Material*)this->hmaterial->delivers();

	return tria;
}
/*}}}*/
void Tria::Marshall(char** pmarshalled_data,int* pmarshalled_data_size, int marshall_direction){ /*{{{*/

	MARSHALLING_ENUM(TriaEnum);
	MARSHALLING(this->iscollapsed);
	MARSHALLING(this->isonsurface);
	MARSHALLING(this->isonbase);

	/*Call parent classes: */
	ElementHook::Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);
	Element::MarshallElement(pmarshalled_data,pmarshalled_data_size,marshall_direction,this->numanalyses);
	TriaRef::Marshall(pmarshalled_data,pmarshalled_data_size,marshall_direction);

	vertices = (Vertex**)this->hvertices->deliverp();
	material = (Material*)this->hmaterial->delivers();

}
/*}}}*/

/*Other*/
void       Tria::AddBasalInput2(int input_enum,IssmDouble* values, int interpolation_enum){/*{{{*/

	/*Call inputs method*/
	_assert_(this->inputs2);

	int domaintype;
	parameters->FindParam(&domaintype,DomainTypeEnum);
	switch(domaintype){
		case Domain2DhorizontalEnum:
			this->AddInput2(input_enum,values,interpolation_enum);
			break;
		case Domain2DverticalEnum:{
			_error_("not implemented yet");
										  }
			break;
		default: _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");
	}

}
/*}}}*/
void       Tria::AddInput2(int input_enum,IssmDouble* values, int interpolation_enum){/*{{{*/

	/*Intermediaries*/
	int vertexlids[NUMVERTICES];

	/*Call inputs method*/
	if(!this->inputs2){
		int* temp = xNew<int>(3);
		_error_("inputs2 not set");
	}
	_assert_(this->inputs2);
	switch(interpolation_enum){
		case P1Enum:
			for(int i=0;i<NUMVERTICES;i++) vertexlids[i]=this->vertices[i]->lid;
			inputs2->SetTriaInput(input_enum,interpolation_enum,NUMVERTICES,vertexlids,values);
			break;
		case P1DGEnum:
			for(int i=0;i<NUMVERTICES;i++) vertexlids[i]=this->vertices[i]->lid;
			inputs2->SetTriaInput(input_enum,interpolation_enum,this->lid,NUMVERTICES,values);
			break;
		default:
			inputs2->SetTriaInput(input_enum,interpolation_enum,this->lid,this->GetNumberOfNodes(interpolation_enum),values);
	}

}
/*}}}*/
void       Tria::AddControlInput(int input_enum,Inputs2* inputs2,IoModel* iomodel,IssmDouble* values,IssmDouble* values_min,IssmDouble* values_max, int interpolation_enum,int id){/*{{{*/

	/*Intermediaries*/
	int vertexlids[NUMVERTICES];

	_assert_(iomodel->elements);
	for(int i=0;i<NUMVERTICES;i++){
		int vertexid =reCast<int>(iomodel->elements[NUMVERTICES*this->Sid()+i]); //ids for vertices are in the elements array from Matlab
		vertexlids[i]=iomodel->my_vertices_lids[vertexid-1];
	}

	/*Call inputs method*/
	switch(interpolation_enum){
		case P1Enum:
			inputs2->SetTriaControlInput(input_enum,TriaInput2Enum,interpolation_enum,id,NUMVERTICES,vertexlids,values,values_min,values_max);
			break;
		default:
			_error_("Cannot add \""<<EnumToStringx(input_enum)<<"\" interpolation "<<EnumToStringx(interpolation_enum)<<" not supported");
	}

}
/*}}}*/
void       Tria::DatasetInputCreate(IssmDouble* array,int M,int N,int* individual_enums,int num_inputs,Inputs2* inputs2,IoModel* iomodel,int input_enum){/*{{{*/

	/*Intermediaries*/
	int        vertexsids[NUMVERTICES];
	int        vertexlids[NUMVERTICES];
	IssmDouble nodeinputs[NUMVERTICES];

	/*Some sanity checks*/
	if(num_inputs<1)                 _error_("Cannot create a DatasetInput of size <1");
	if(M!=iomodel->numberofvertices) _error_("Input size not supported yet");
	if(N!=num_inputs)                _error_("Sizes are not consistent");

	/*Get indices*/
	_assert_(iomodel->elements);
	for(int i=0;i<NUMVERTICES;i++){
		vertexsids[i] = reCast<int>(iomodel->elements[NUMVERTICES*this->Sid()+i])-1;
		vertexlids[i] = iomodel->my_vertices_lids[vertexsids[i]];
	}

	/*Create inputs and add to DataSetInput*/
	for(int i=0;i<num_inputs;i++){
		for(int j=0;j<NUMVERTICES;j++) nodeinputs[j]=array[vertexsids[j]*N+i];
		inputs2->SetTriaDatasetInput(input_enum,individual_enums[i],P1Enum,NUMVERTICES,vertexlids,nodeinputs);
	}
}
/*}}}*/
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->GetVerticesPidList(&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::CalvingRateVonmises(){/*{{{*/

	IssmDouble  xyz_list[NUMVERTICES][3];
	IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
	IssmDouble  calvingratex[NUMVERTICES];
	IssmDouble  calvingratey[NUMVERTICES];
	IssmDouble  calvingrate[NUMVERTICES];
	IssmDouble  lambda1,lambda2,ex,ey,vx,vy,vel;
	IssmDouble  sigma_vm[NUMVERTICES];
	IssmDouble  B,sigma_max,sigma_max_floating,sigma_max_grounded,n;
	IssmDouble  epse_2,groundedice,bed,sealevel;		// added sealevel


	/* Get node coordinates and dof list: */
	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);

	/*Retrieve all inputs and parameters we will need*/
	Input2* vx_input = this->GetInput2(VxEnum); _assert_(vx_input);
	Input2* vy_input = this->GetInput2(VyEnum); _assert_(vy_input);
	Input2* B_input  = this->GetInput2(MaterialsRheologyBbarEnum);   _assert_(B_input);
	Input2* gr_input = this->GetInput2(MaskGroundediceLevelsetEnum); _assert_(gr_input);
	Input2* bs_input = this->GetInput2(BaseEnum);                    _assert_(bs_input);
	Input2* smax_fl_input = this->GetInput2(CalvingStressThresholdFloatingiceEnum); _assert_(smax_fl_input);
	Input2* smax_gr_input = this->GetInput2(CalvingStressThresholdGroundediceEnum); _assert_(smax_gr_input);
	Input2* n_input  = this->GetInput2(MaterialsRheologyNEnum); _assert_(n_input);
	Input2* sl_input  = this->GetInput2(SealevelEnum); _assert_(sl_input);



	/* Start looping on the number of vertices: */
	GaussTria* gauss=new GaussTria();
	for(int iv=0;iv<NUMVERTICES;iv++){
		gauss->GaussVertex(iv);

		/*Get velocity components and thickness*/
		B_input->GetInputValue(&B,gauss);
		n_input->GetInputValue(&n,gauss);
		vx_input->GetInputValue(&vx,gauss);
		vy_input->GetInputValue(&vy,gauss);
		gr_input->GetInputValue(&groundedice,gauss);
		bs_input->GetInputValue(&bed,gauss);
		smax_fl_input->GetInputValue(&sigma_max_floating,gauss);
		smax_gr_input->GetInputValue(&sigma_max_grounded,gauss);
		vel=sqrt(vx*vx+vy*vy)+1.e-14;
		sl_input->GetInputValue(&sealevel,gauss);

		/*Compute strain rate and viscosity: */
		this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);

		/*Get Eigen values*/
		Matrix2x2Eigen(&lambda1,&lambda2,&ex,&ey,epsilon[0],epsilon[2],epsilon[1]);
		_assert_(!xIsNan<IssmDouble>(lambda1));
		_assert_(!xIsNan<IssmDouble>(lambda2));

		/*Process Eigen values (only account for extension)*/
		lambda1 = max(lambda1,0.);
		lambda2 = max(lambda2,0.);

		/*Calculate sigma_vm*/
		epse_2    = 1./2. *(lambda1*lambda1 + lambda2*lambda2);
		sigma_vm[iv]  = sqrt(3.) * B * pow(epse_2,1./(2.*n));

		/*OLD (keep for a little bit)*/
		//sigma_max = 800.e+3; //IUGG previous test
		//sigma_max = 1000.e+3; //GRL
		//if(groundedice<0) sigma_max=150.e+3;

		/*Tensile stress threshold*/
		if(groundedice<0)
		 sigma_max = sigma_max_floating;
		else
		 sigma_max = sigma_max_grounded;
		/*Assign values*/
		if(bed>sealevel){		// Changed 0. to sealevel
			calvingratex[iv]=0.;
			calvingratey[iv]=0.;
		}
		else{
			calvingratex[iv]=vx*sigma_vm[iv]/sigma_max;
			calvingratey[iv]=vy*sigma_vm[iv]/sigma_max;
		}
		calvingrate[iv] =sqrt(calvingratex[iv]*calvingratex[iv] + calvingratey[iv]*calvingratey[iv]);

	}

	/*Add input*/
	this->AddInput2(CalvingratexEnum,&calvingratex[0],P1DGEnum);
	this->AddInput2(CalvingrateyEnum,&calvingratey[0],P1DGEnum);
	this->AddInput2(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum);
	this->AddInput2(SigmaVMEnum,&sigma_vm[0],P1DGEnum);

	/*Clean up and return*/
	delete gauss;
}
/*}}}*/
void       Tria::CalvingCrevasseDepth(){/*{{{*/

	IssmDouble  xyz_list[NUMVERTICES][3];
	IssmDouble  calvingrate[NUMVERTICES];
	IssmDouble  vx,vy,vel;
	IssmDouble  water_height, bed,Ho,thickness,surface;
	IssmDouble  surface_crevasse[NUMVERTICES], basal_crevasse[NUMVERTICES], crevasse_depth[NUMVERTICES], H_surf, H_surfbasal;
	IssmDouble  B, strainparallel, straineffective,n;
	IssmDouble  s_xx,s_xy,s_yy,s1,s2,stmp;
	int crevasse_opening_stress;

	/* Get node coordinates and dof list: */
	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);

	/*retrieve the type of crevasse_opening_stress*/
	this->parameters->FindParam(&crevasse_opening_stress,CalvingCrevasseDepthEnum);

	IssmDouble rho_ice        = this->FindParam(MaterialsRhoIceEnum);
	IssmDouble rho_seawater   = this->FindParam(MaterialsRhoSeawaterEnum);
	IssmDouble rho_freshwater = this->FindParam(MaterialsRhoFreshwaterEnum);
	IssmDouble constant_g     = this->FindParam(ConstantsGEnum);

	Input2*   H_input                 = this->GetInput2(ThicknessEnum); _assert_(H_input);
	Input2*   bed_input               = this->GetInput2(BedEnum); _assert_(bed_input);
	Input2*   surface_input           = this->GetInput2(SurfaceEnum); _assert_(surface_input);
	Input2*	strainrateparallel_input  = this->GetInput2(StrainRateparallelEnum);  _assert_(strainrateparallel_input);
	Input2*	strainrateeffective_input = this->GetInput2(StrainRateeffectiveEnum); _assert_(strainrateeffective_input);
	Input2*	vx_input                  = this->GetInput2(VxEnum); _assert_(vx_input);
	Input2*	vy_input                  = this->GetInput2(VxEnum); _assert_(vy_input);
	Input2*   waterheight_input       = this->GetInput2(WaterheightEnum); _assert_(waterheight_input);
	Input2*   s_xx_input              = this->GetInput2(DeviatoricStressxxEnum);     _assert_(s_xx_input);
	Input2*   s_xy_input              = this->GetInput2(DeviatoricStressxyEnum);     _assert_(s_xy_input);
	Input2*   s_yy_input              = this->GetInput2(DeviatoricStressyyEnum);     _assert_(s_yy_input);
	Input2*	B_input  = this->GetInput2(MaterialsRheologyBbarEnum);   _assert_(B_input);
	Input2*	n_input  = this->GetInput2(MaterialsRheologyNEnum);   _assert_(n_input);

	/*Loop over all elements of this partition*/
	GaussTria* gauss=new GaussTria();
	for (int iv=0;iv<NUMVERTICES;iv++){
		gauss->GaussVertex(iv);

		H_input->GetInputValue(&thickness,gauss);
		bed_input->GetInputValue(&bed,gauss);
		surface_input->GetInputValue(&surface,gauss);
		strainrateparallel_input->GetInputValue(&strainparallel,gauss);
		strainrateeffective_input->GetInputValue(&straineffective,gauss);
		vx_input->GetInputValue(&vx,gauss);
		vy_input->GetInputValue(&vy,gauss);
		waterheight_input->GetInputValue(&water_height,gauss);
		s_xx_input->GetInputValue(&s_xx,gauss);
		s_xy_input->GetInputValue(&s_xy,gauss);
		s_yy_input->GetInputValue(&s_yy,gauss);
		B_input->GetInputValue(&B,gauss);
		n_input->GetInputValue(&n,gauss);

		vel=sqrt(vx*vx+vy*vy)+1.e-14;

		s1=(s_xx+s_yy)/2.+sqrt(pow((s_xx-s_yy)/2.,2)+pow(s_xy,2));
		s2=(s_xx+s_yy)/2.-sqrt(pow((s_xx-s_yy)/2.,2)+pow(s_xy,2));
		if(fabs(s2)>fabs(s1)){stmp=s2; s2=s1; s1=stmp;}

		Ho = thickness - (rho_seawater/rho_ice) * (-bed);
		if(Ho<0.)  Ho=0.;

		if(crevasse_opening_stress==0){		/*Otero2010: balance between the tensile deviatoric stress and ice overburden pressure*/
			surface_crevasse[iv] = B * strainparallel * pow(straineffective, ((1 / n)-1)) / (rho_ice * constant_g);
			basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice)) * (B * strainparallel * pow(straineffective,((1/n)-1)) / (rho_ice*constant_g) - Ho);
		}
		else if(crevasse_opening_stress==1){	 /*Benn2017,Todd2018: maximum principal stress */
			surface_crevasse[iv] = s1 / (rho_ice*constant_g);
			basal_crevasse[iv] = (rho_ice/(rho_seawater-rho_ice))* (s1/ (rho_ice*constant_g)-Ho);
		}

		/* some constraints */
		if (surface_crevasse[iv]<0.) {
			surface_crevasse[iv]=0.;
			water_height = 0.;
		}
		if (basal_crevasse[iv]<0.) basal_crevasse[iv]=0.;
		if (bed>0.) basal_crevasse[iv] = 0.;

		//if (surface_crevasse[iv]<water_height){
		//	water_height = surface_crevasse[iv];
		//}

		/* add water in surface crevasse */
		surface_crevasse[iv] = surface_crevasse[iv] + (rho_freshwater/rho_ice)*water_height; /* surface crevasse + water */
		crevasse_depth[iv] = surface_crevasse[iv] + (rho_freshwater/rho_ice)*water_height + basal_crevasse[iv]; /* surface crevasse + basal crevasse + water */

	}

	this->AddInput2(SurfaceCrevasseEnum,&surface_crevasse[0],P1DGEnum);
	this->AddInput2(BasalCrevasseEnum,&basal_crevasse[0],P1DGEnum);
	this->AddInput2(CrevasseDepthEnum,&crevasse_depth[0],P1DGEnum);

	delete gauss;
}
/*}}}*/
void       Tria::CalvingRateLevermann(){/*{{{*/

	IssmDouble  xyz_list[NUMVERTICES][3];
	IssmDouble  vx,vy,vel;
	IssmDouble  strainparallel;
	IssmDouble  propcoeff,bed;
	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*/
	Input2* vx_input=this->GetInput2(VxEnum);													_assert_(vx_input);
	Input2* vy_input=this->GetInput2(VyEnum);													_assert_(vy_input);
	Input2* bs_input = this->GetInput2(BaseEnum);                                 _assert_(bs_input);
	Input2* strainparallel_input=this->GetInput2(StrainRateparallelEnum);			_assert_(strainparallel_input);
	Input2* strainperpendicular_input=this->GetInput2(StrainRateperpendicularEnum);_assert_(strainperpendicular_input);
	Input2* levermanncoeff_input=this->GetInput2(CalvinglevermannCoeffEnum);      _assert_(levermanncoeff_input);

	/* Start looping on the number of vertices: */
	GaussTria* 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);
		bs_input->GetInputValue(&bed,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 */
		if(strainparallel>0. && strainperpendicular>0. && bed<=0.){
			calvingrate[iv]=propcoeff*strainparallel*strainperpendicular;
		}
		else
			calvingrate[iv]=0.;

		calvingratex[iv]=calvingrate[iv]*vx/(sqrt(vel)+1.e-14);
		calvingratey[iv]=calvingrate[iv]*vy/(sqrt(vel)+1.e-14);
	}

	/*Add input*/
	this->AddInput2(CalvingratexEnum,&calvingratex[0],P1DGEnum);
	this->AddInput2(CalvingrateyEnum,&calvingratey[0],P1DGEnum);
	this->AddInput2(CalvingCalvingrateEnum,&calvingrate[0],P1DGEnum);

	/*Clean up and return*/
	delete gauss;

}
/*}}}*/
void       Tria::CalvingFluxLevelset(){/*{{{*/

	/*Make sure there is an ice front here*/
	if(!IsIceInElement() || !IsZeroLevelset(MaskIceLevelsetEnum)){
		IssmDouble flux_per_area=0;
		this->AddInput2(CalvingFluxLevelsetEnum,&flux_per_area,P0Enum);
	}
	else{
		int               domaintype,index1,index2;
		const IssmPDouble epsilon = 1.e-15;
		IssmDouble        s1,s2;
		IssmDouble        gl[NUMVERTICES];
		IssmDouble        xyz_front[2][3];

		IssmDouble *xyz_list = NULL;
		this->GetVerticesCoordinates(&xyz_list);

		/*Recover parameters and values*/
		parameters->FindParam(&domaintype,DomainTypeEnum);
		Element::GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);

		/*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;

		if(domaintype==Domain2DverticalEnum){
			_error_("not implemented");
		}
		else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum || domaintype==Domain3DsurfaceEnum){
			int pt1 = 0;
			int pt2 = 1;
			if(gl[0]*gl[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=gl[2]/(gl[2]-gl[1]);
				s2=gl[2]/(gl[2]-gl[0]);
				if(gl[2]<0.){
					pt1 = 1; pt2 = 0;
				}
				xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
				xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
				xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
				xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
				xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
				xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
			}
			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

				/*Portion of the segments*/
				s1=gl[0]/(gl[0]-gl[1]);
				s2=gl[0]/(gl[0]-gl[2]);
				if(gl[0]<0.){
					pt1 = 1; pt2 = 0;
				}

				xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
				xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
				xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
				xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
				xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
				xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+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

				/*Portion of the segments*/
				s1=gl[1]/(gl[1]-gl[0]);
				s2=gl[1]/(gl[1]-gl[2]);
				if(gl[1]<0.){
					pt1 = 1; pt2 = 0;
				}

				xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
				xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
				xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
				xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
				xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
				xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
			}
			else{
				_error_("case not possible");
			}

		}
		else _error_("mesh type "<<EnumToStringx(domaintype)<<"not supported yet ");

		/*Some checks in debugging mode*/
		_assert_(s1>=0 && s1<=1.);
		_assert_(s2>=0 && s2<=1.);

		/*Get normal vector*/
		IssmDouble normal[3];
		this->NormalSection(&normal[0],&xyz_front[0][0]);
		normal[0] = -normal[0];
		normal[1] = -normal[1];

		/*Get inputs*/
		IssmDouble flux = 0.;
		IssmDouble area = 0.;
		IssmDouble calvingratex,calvingratey,thickness,Jdet,flux_per_area;
		IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
		Input2* thickness_input=this->GetInput2(ThicknessEnum); _assert_(thickness_input);
		Input2* calvingratex_input=NULL;
		Input2* calvingratey_input=NULL;
		if(domaintype==Domain2DhorizontalEnum){
			calvingratex_input=this->GetInput2(CalvingratexEnum); _assert_(calvingratex_input);
			calvingratey_input=this->GetInput2(CalvingrateyEnum); _assert_(calvingratey_input);
		}
		else{
			calvingratex_input=this->GetInput2(CalvingratexAverageEnum); _assert_(calvingratex_input);
			calvingratey_input=this->GetInput2(CalvingrateyAverageEnum); _assert_(calvingratey_input);
		}

		/*Start looping on Gaussian points*/
		Gauss* gauss=this->NewGauss(xyz_list,&xyz_front[0][0],3);
		for(int ig=gauss->begin();ig<gauss->end();ig++){

			gauss->GaussPoint(ig);
			thickness_input->GetInputValue(&thickness,gauss);
			calvingratex_input->GetInputValue(&calvingratex,gauss);
			calvingratey_input->GetInputValue(&calvingratey,gauss);
			this->JacobianDeterminantSurface(&Jdet,&xyz_front[0][0],gauss);

			flux += rho_ice*Jdet*gauss->weight*thickness*(calvingratex*normal[0] + calvingratey*normal[1]);
			area += Jdet*gauss->weight*thickness;

			flux_per_area=flux/area;
		}

		this->AddInput2(CalvingFluxLevelsetEnum,&flux_per_area,P0Enum);

		/*Clean up and return*/
		delete gauss;
	}
}
/*}}}*/
void       Tria::CalvingMeltingFluxLevelset(){/*{{{*/

	/*Make sure there is an ice front here*/
	if(!IsIceInElement() || !IsZeroLevelset(MaskIceLevelsetEnum)){
		IssmDouble flux_per_area=0;
		this->AddInput2(CalvingMeltingFluxLevelsetEnum,&flux_per_area,P0Enum);
	}
	else{
		int               domaintype,index1,index2;
		const IssmPDouble epsilon = 1.e-15;
		IssmDouble        s1,s2;
		IssmDouble        gl[NUMVERTICES];
		IssmDouble        xyz_front[2][3];


		IssmDouble *xyz_list = NULL;
		this->GetVerticesCoordinates(&xyz_list);

		/*Recover parameters and values*/
		parameters->FindParam(&domaintype,DomainTypeEnum);
		Element::GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);

		/*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;

		if(domaintype==Domain2DverticalEnum){
			_error_("not implemented");
		}
		else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum || domaintype==Domain3DsurfaceEnum){
			int pt1 = 0;
			int pt2 = 1;
			if(gl[0]*gl[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=gl[2]/(gl[2]-gl[1]);
				s2=gl[2]/(gl[2]-gl[0]);
				if(gl[2]<0.){
					pt1 = 1; pt2 = 0;
				}
				xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
				xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
				xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
				xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
				xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
				xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
			}
			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

				/*Portion of the segments*/
				s1=gl[0]/(gl[0]-gl[1]);
				s2=gl[0]/(gl[0]-gl[2]);
				if(gl[0]<0.){
					pt1 = 1; pt2 = 0;
				}

				xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
				xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
				xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
				xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
				xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
				xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+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

				/*Portion of the segments*/
				s1=gl[1]/(gl[1]-gl[0]);
				s2=gl[1]/(gl[1]-gl[2]);
				if(gl[1]<0.){
					pt1 = 1; pt2 = 0;
				}

				xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
				xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
				xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
				xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
				xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
				xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
			}
			else{
				_error_("case not possible");
			}

		}
		else _error_("mesh type "<<EnumToStringx(domaintype)<<"not supported yet ");

		/*Some checks in debugging mode*/
		_assert_(s1>=0 && s1<=1.);
		_assert_(s2>=0 && s2<=1.);

		/*Get normal vector*/
		IssmDouble normal[3];
		this->NormalSection(&normal[0],&xyz_front[0][0]);
		normal[0] = -normal[0];
		normal[1] = -normal[1];

		/*Get inputs*/
		IssmDouble flux = 0.;
		IssmDouble area = 0.;
		IssmDouble calvingratex,calvingratey,vx,vy,vel,meltingrate,meltingratex,meltingratey,thickness,Jdet,flux_per_area;
		IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
		Input2* thickness_input=this->GetInput2(ThicknessEnum); _assert_(thickness_input);
		Input2* calvingratex_input=NULL;
		Input2* calvingratey_input=NULL;
		Input2* vx_input=NULL;
		Input2* vy_input=NULL;
		Input2* meltingrate_input=NULL;
		if(domaintype==Domain2DhorizontalEnum){
			calvingratex_input=this->GetInput2(CalvingratexEnum); _assert_(calvingratex_input);
			calvingratey_input=this->GetInput2(CalvingrateyEnum); _assert_(calvingratey_input);
			vx_input=this->GetInput2(VxEnum); _assert_(vx_input);
			vy_input=this->GetInput2(VyEnum); _assert_(vy_input);
			meltingrate_input=this->GetInput2(CalvingMeltingrateEnum); _assert_(meltingrate_input);
		}
		else{
			calvingratex_input=this->GetInput2(CalvingratexAverageEnum); _assert_(calvingratex_input);
			calvingratey_input=this->GetInput2(CalvingrateyAverageEnum); _assert_(calvingratey_input);
		}

		/*Start looping on Gaussian points*/
		Gauss* gauss=this->NewGauss(xyz_list,&xyz_front[0][0],3);
		for(int ig=gauss->begin();ig<gauss->end();ig++){

			gauss->GaussPoint(ig);
			thickness_input->GetInputValue(&thickness,gauss);
			calvingratex_input->GetInputValue(&calvingratex,gauss);
			calvingratey_input->GetInputValue(&calvingratey,gauss);
			vx_input->GetInputValue(&vx,gauss);
			vy_input->GetInputValue(&vy,gauss);
			vel=vx*vx+vy*vy;
			meltingrate_input->GetInputValue(&meltingrate,gauss);
			meltingratex=meltingrate*vx/(sqrt(vel)+1.e-14);
			meltingratey=meltingrate*vy/(sqrt(vel)+1.e-14);
			this->JacobianDeterminantSurface(&Jdet,&xyz_front[0][0],gauss);

			flux += rho_ice*Jdet*gauss->weight*thickness*((calvingratex+meltingratex)*normal[0] + (calvingratey+meltingratey)*normal[1]);
			area += Jdet*gauss->weight*thickness;

			flux_per_area=flux/area;
		}

		this->AddInput2(CalvingMeltingFluxLevelsetEnum,&flux_per_area,P0Enum);

		/*Clean up and return*/
		delete gauss;
	}
}
/*}}}*/
IssmDouble Tria::CharacteristicLength(void){/*{{{*/

	return sqrt(2*this->GetArea());
}
/*}}}*/
void       Tria::ComputeBasalStress(void){/*{{{*/
	_error_("Not Implemented yet");
}
/*}}}*/
void       Tria::ComputeDeviatoricStressTensor(){/*{{{*/

	IssmDouble  xyz_list[NUMVERTICES][3];
	IssmDouble  viscosity,lambda1,lambda2;
	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};
	IssmDouble  tau_e[NUMVERTICES];
	IssmDouble  tau_1[NUMVERTICES];
	IssmDouble  tau_2[NUMVERTICES];
	GaussTria*  gauss=NULL;
	int domaintype,dim=2;

	/*Get approximation*/
	int approximation;
	this->GetInput2Value(&approximation,ApproximationEnum);

	/* Get node coordinates and dof list: */
	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);

	/*Retrieve all inputs we will be needing: */
	this->FindParam(&domaintype,DomainTypeEnum);
	Input2* vx_input=this->GetInput2(VxEnum); _assert_(vx_input);
	Input2* vy_input=this->GetInput2(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);
		switch(approximation){
			case SSAApproximationEnum:
				this->material->ViscositySSA(&viscosity,dim,&xyz_list[0][0],gauss,vx_input,vy_input);
				break;
			case HOApproximationEnum:
				this->material->ViscosityHO(&viscosity,dim,&xyz_list[0][0],gauss,vx_input,vy_input);
				break;
			case FSApproximationEnum:
				this->material->ViscosityFS(&viscosity,dim,&xyz_list[0][0],gauss,vx_input,vy_input,NULL);
				break;
			default:
				_error_("not supported yet");
		}

		/*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];
		tau_e[iv]=1/sqrt(2)*sqrt(pow(tau_xx[iv],2)+pow(tau_yy[iv],2)+2*pow(tau_xy[iv],2));

		/*Get Eigen values*/
		Matrix2x2Eigen(&tau_2[iv],&tau_1[iv],NULL,NULL,tau_xx[iv],tau_xy[iv],tau_yy[iv]);
	}

	/*Add Stress tensor components into inputs*/
	this->AddInput2(DeviatoricStressxxEnum,&tau_xx[0],P1DGEnum);
	this->AddInput2(DeviatoricStressxyEnum,&tau_xy[0],P1DGEnum);
	this->AddInput2(DeviatoricStressxzEnum,&tau_xz[0],P1DGEnum);
	this->AddInput2(DeviatoricStressyyEnum,&tau_yy[0],P1DGEnum);
	this->AddInput2(DeviatoricStressyzEnum,&tau_yz[0],P1DGEnum);
	this->AddInput2(DeviatoricStresszzEnum,&tau_zz[0],P1DGEnum);
	this->AddInput2(DeviatoricStresseffectiveEnum,&tau_e[0],P1DGEnum);
	this->AddInput2(DeviatoricStress1Enum,&tau_1[0],P1DGEnum);
	this->AddInput2(DeviatoricStress2Enum,&tau_2[0],P1DGEnum);

	/*Clean up and return*/
	delete gauss;
}
/*}}}*/
void			    Tria::ComputeEsaStrainAndVorticity(){ /*{{{*/

	IssmDouble  xyz_list[NUMVERTICES][3];
	IssmDouble  epsilon[4]; /* epsilon=[exx,eyy,exy+ (shear),exy- (rotation)];*/
	IssmDouble  strain_xx[NUMVERTICES];
	IssmDouble  strain_yy[NUMVERTICES];
	IssmDouble  strain_xy[NUMVERTICES];
	IssmDouble  vorticity_xy[NUMVERTICES];
	GaussTria*  gauss=NULL;

	/* Get node coordinates and dof list: */
	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);

	/*Retrieve all inputs we will be needing: */
	Input2* vx_input=this->GetInput2(EsaXmotionEnum); _assert_(vx_input);
	Input2* vy_input=this->GetInput2(EsaYmotionEnum); _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 vorticity rate: */
		this->StrainRateESA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);

		/*Compute Stress*/
		strain_xx[iv]=epsilon[0];
		strain_yy[iv]=epsilon[1];
		strain_xy[iv]=epsilon[2];
		vorticity_xy[iv]=epsilon[3];
	}

	/*Add Stress tensor components into inputs*/
	this->AddInput2(EsaStrainratexxEnum,&strain_xx[0],P1DGEnum);
	this->AddInput2(EsaStrainrateyyEnum,&strain_yy[0],P1DGEnum);
	this->AddInput2(EsaStrainratexyEnum,&strain_xy[0],P1DGEnum);
	this->AddInput2(EsaRotationrateEnum,&vorticity_xy[0],P1DGEnum);

	/*Clean up and return*/
	delete gauss;
}
/*}}}*/
void       Tria::ComputeSigmaNN(){/*{{{*/

	if(!IsOnBase()){
		IssmDouble sigma_nn[3]={0.};
		this->AddInput2(SigmaNNEnum,&sigma_nn[0],P1Enum);
		return;
	}
	else{
		IssmDouble* xyz_list=NULL;
		IssmDouble *xyz_list_base=NULL;
		IssmDouble  pressure,viscosity;
		IssmDouble  sigma_nn[3];
		IssmDouble  sigma_xx,sigma_xy,sigma_yy;
		IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
		IssmDouble  base_normal[2];
		int domaintype,dim=2;

		/* Get node coordinates and dof list: */
		GetVerticesCoordinates(&xyz_list);
	   GetVerticesCoordinatesBase(&xyz_list_base);

		/*Retrieve all inputs we will be needing: */
		this->FindParam(&domaintype,DomainTypeEnum);
		if(domaintype==Domain2DhorizontalEnum) _error_("stress tensor calculation not supported for mesh of type " <<EnumToStringx(domaintype)<<", extrude mesh or call ComputeDeviatoricStressTensor");
		Input2* pressure_input=this->GetInput2(PressureEnum); _assert_(pressure_input);
		Input2* vx_input=this->GetInput2(VxEnum);             _assert_(vx_input);
		Input2* vy_input=this->GetInput2(VyEnum);             _assert_(vy_input);

		/* Start looping on the number of vertices: */
		Gauss* gauss = this->NewGauss();
		for(int i=0;i<NUMVERTICES;i++){
			gauss->GaussNode(P1Enum,i);

			/*Compute strain rate viscosity and pressure: */
			this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
			this->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,NULL);
			pressure_input->GetInputValue(&pressure,gauss);

			/*Compute Stress*/
			sigma_xx=2*viscosity*epsilon[0]-pressure; // sigma = nu eps - pressure
			sigma_yy=2*viscosity*epsilon[1]-pressure;
			sigma_xy=2*viscosity*epsilon[2];

			/*Get normal vector to the bed */
			NormalBase(&base_normal[0],xyz_list_base);

			/*Compute sigma_nn*/
			sigma_nn[i]=sigma_xx*base_normal[0]*base_normal[0] + 2*sigma_xy*base_normal[0]*base_normal[1] + sigma_yy*base_normal[1]*base_normal[1];
		}

		/*Add Stress tensor components into inputs*/
		this->AddInput2(SigmaNNEnum,&sigma_nn[0],P1Enum);

		/*Clean up and return*/
		xDelete<IssmDouble>(xyz_list);
		xDelete<IssmDouble>(xyz_list_base);
		delete gauss;
	}
}
/*}}}*/
void       Tria::ComputeStressTensor(){/*{{{*/

	IssmDouble  xyz_list[NUMVERTICES][3];
	IssmDouble  pressure,viscosity;
	IssmDouble  epsilon[3]; /* epsilon=[exx,eyy,exy];*/
	IssmDouble  sigma_xx[NUMVERTICES];
	IssmDouble	sigma_yy[NUMVERTICES];
	IssmDouble	sigma_zz[NUMVERTICES]={0,0,0};
	IssmDouble  sigma_xy[NUMVERTICES];
	IssmDouble	sigma_xz[NUMVERTICES]={0,0,0};
	IssmDouble	sigma_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_("stress tensor calculation not supported for mesh of type " <<EnumToStringx(domaintype)<<", extrude mesh or call ComputeDeviatoricStressTensor");
	Input2* pressure_input=this->GetInput2(PressureEnum); _assert_(pressure_input);
	Input2* vx_input=this->GetInput2(VxEnum);             _assert_(vx_input);
	Input2* vy_input=this->GetInput2(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 viscosity and pressure: */
		this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
		this->material->ViscositySSA(&viscosity,dim,&xyz_list[0][0],gauss,vx_input,vy_input);
		pressure_input->GetInputValue(&pressure,gauss);

		/*Compute Stress*/
		sigma_xx[iv]=2*viscosity*epsilon[0]-pressure; // sigma = nu eps - pressure
		sigma_yy[iv]=2*viscosity*epsilon[1]-pressure;
		sigma_xy[iv]=2*viscosity*epsilon[2];
	}

	/*Add Stress tensor components into inputs*/
	this->AddInput2(StressTensorxxEnum,&sigma_xx[0],P1DGEnum);
	this->AddInput2(StressTensorxyEnum,&sigma_xy[0],P1DGEnum);
	this->AddInput2(StressTensorxzEnum,&sigma_xz[0],P1DGEnum);
	this->AddInput2(StressTensoryyEnum,&sigma_yy[0],P1DGEnum);
	this->AddInput2(StressTensoryzEnum,&sigma_yz[0],P1DGEnum);
	this->AddInput2(StressTensorzzEnum,&sigma_zz[0],P1DGEnum);

	/*Clean up and return*/
	delete gauss;
}
/*}}}*/
void       Tria::Configure(Elements* elementsin, Loads* loadsin,Nodes* nodesin,Vertices *verticesin,Materials* materialsin, Parameters* parametersin,Inputs2* inputs2in){/*{{{*/

	/*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];

	/*Take care of hooking up all objects for this element, ie links the objects in the hooks to their respective
	 * datasets, using internal ids and offsets hidden in hooks: */
	if(this->hnodes){
		if (this->hnodes[analysis_counter]) this->hnodes[analysis_counter]->configure(nodesin);
		else this->hnodes[analysis_counter] = NULL;
	}
	else this->hnodes = NULL;
	this->hvertices->configure(verticesin);
	this->hmaterial->configure(materialsin);

	/*Now, go pick up the objects inside the hooks: */
	if(this->hnodes && this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp();
	else this->nodes=NULL;
	this->vertices = (Vertex**)this->hvertices->deliverp();
	this->material = (Material*)this->hmaterial->delivers();

	/*point parameters to real dataset: */
	this->parameters=parametersin;
	this->inputs2=inputs2in;
}/*}}}*/
void       Tria::ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index,int offset,int N, int M){/*{{{*/

	int         idlist[NUMVERTICES];
	int         vertexlids[NUMVERTICES];
	int         gradidlist[NUMVERTICES];
	IssmDouble  grad_list[NUMVERTICES];

	GradientIndexing(&gradidlist[0],control_index);
	for(int i=0;i<NUMVERTICES;i++) vertexlids[i]=this->vertices[i]->lid;

	if(N==1){
		this->inputs2->SetTriaControlInputGradient(enum_type,P1Enum,NUMVERTICES,&vertexlids[0],&grad_list[0]);
	}
	else{
		for(int n=0;n<N;n++){
			for(int i=0;i<NUMVERTICES;i++){
				idlist[i] = offset + this->vertices[i]->Sid()+n*M;
				grad_list[i]=gradient[idlist[i]];
			}
			this->inputs2->SetTriaControlInputGradient(enum_type,P1Enum,NUMVERTICES,&vertexlids[0],&grad_list[0],n);
		}
	}

}/*}}}*/
void       Tria::ControlInputSetGradient(IssmDouble* gradient,int enum_type,int control_index){/*{{{*/

	int        idlist[NUMVERTICES];
	int        vertexlids[NUMVERTICES];
	IssmDouble grad_list[NUMVERTICES];

	GradientIndexing(&idlist[0],control_index);
	for(int i=0;i<NUMVERTICES;i++) grad_list[i]=gradient[idlist[i]];
	for(int i=0;i<NUMVERTICES;i++) vertexlids[i]=this->vertices[i]->lid;

	this->inputs2->SetTriaControlInputGradient(enum_type,P1Enum,NUMVERTICES,&vertexlids[0],&grad_list[0]);

}/*}}}*/
void       Tria::ControlToVectors(Vector<IssmPDouble>* vector_control, Vector<IssmPDouble>* vector_gradient,int control_enum){/*{{{*/

	int         sidlist[NUMVERTICES];
	int         lidlist[NUMVERTICES];
	int         connectivity[NUMVERTICES];
	IssmPDouble values[NUMVERTICES];
	IssmPDouble gradients[NUMVERTICES];
	IssmDouble  value,gradient;

	this->GetVerticesConnectivityList(&connectivity[0]);
	this->GetVerticesSidList(&sidlist[0]);
	this->GetVerticesLidList(&lidlist[0]);

	ElementInput2* control_value    = this->inputs2->GetControlInput2Data(control_enum,"value");    _assert_(control_value);
	ElementInput2* control_gradient = this->inputs2->GetControlInput2Data(control_enum,"gradient"); _assert_(control_gradient);
	control_value->Serve(NUMVERTICES,&lidlist[0]);
	control_gradient->Serve(NUMVERTICES,&lidlist[0]);

	GaussTria* gauss=new GaussTria();
	for (int iv=0;iv<NUMVERTICES;iv++){
		gauss->GaussVertex(iv);

		control_value->GetInputValue(&value,gauss);
		control_gradient->GetInputValue(&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::CreateDistanceInputFromSegmentlist(IssmDouble* distances,int distanceenum){/*{{{*/

	/*Get current field and vertex coordinates*/
	IssmDouble ls[NUMVERTICES],distance;
	Element::GetInputListOnVertices(&ls[0],distanceenum);

	/*Get distance from list of segments and reset ls*/
	for(int j=0;j<NUMVERTICES;j++){
		distance=distances[this->vertices[j]->Lid()];
		if(xIsNan<IssmDouble>(distance)) _error_("NaN found in vector");
		if(xIsInf<IssmDouble>(distance)) _error_("Inf found in vector");

		/*FIXME: do we really need this?*/
		if(distanceenum==MaskIceLevelsetEnum) if(distance>10000) distance=10000;
		if(ls[j]>0){
			ls[j] = distance;
		}
		else{
			ls[j] = - distance;
		}
	}

	/*Update Levelset*/
	this->AddInput2(distanceenum,&ls[0],P1Enum);
}
/*}}}*/
int        Tria::EdgeOnBaseIndex(void){/*{{{*/

	IssmDouble values[NUMVERTICES];
	int        indices[3][2] = {{1,2},{2,0},{0,1}};

	/*Retrieve all inputs and parameters*/
	Element::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*/
	Element::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*/
	Element::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*/
	Element::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(NULL);
			break;

		case VelEnum:{

			/*Get input:*/
			IssmDouble vel;
			Input2* vel_input=this->GetInput2(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];
	IssmDouble xmin,ymin;
	IssmDouble xmax,ymax;

	/*Get xyz list: */
	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
	xmin=xyz_list[0][0]; xmax=xyz_list[0][0];
	ymin=xyz_list[0][1]; ymax=xyz_list[0][1];

	for(int i=1;i<NUMVERTICES;i++){
		if(xyz_list[i][0]<xmin) xmin=xyz_list[i][0];
		if(xyz_list[i][0]>xmax) xmax=xyz_list[i][0];
		if(xyz_list[i][1]<ymin) ymin=xyz_list[i][1];
		if(xyz_list[i][1]>ymax) ymax=xyz_list[i][1];
	}

	*hx=xmax-xmin;
	*hy=ymax-ymin;
	*hz=0.;
}
/*}}}*/
int        Tria::FiniteElement(void){/*{{{*/
	return this->element_type;
}
/*}}}*/
IssmDouble Tria::FloatingArea(bool scaled){/*{{{*/

	/*Intermediaries*/
	int         domaintype;
	IssmDouble  phi,scalefactor,floatingarea;
	IssmDouble *xyz_list  = NULL;

	if(!IsIceInElement())return 0.;

	/*Get problem dimension*/
	this->FindParam(&domaintype,DomainTypeEnum);
	if(domaintype!=Domain2DhorizontalEnum && domaintype!=Domain3DEnum) _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");

	this->GetVerticesCoordinates(&xyz_list);
	phi=this->GetGroundedPortion(xyz_list);
	floatingarea=(1-phi)*this->GetArea();
	if(scaled==true){
		Input2* scalefactor_input = this->GetInput2(MeshScaleFactorEnum); _assert_(scalefactor_input);
		scalefactor_input->GetInputAverage(&scalefactor);
		floatingarea=floatingarea*scalefactor;
	}

	/*Clean up and return*/
	xDelete<IssmDouble>(xyz_list);
	return floatingarea;
}
/*}}}*/
void       Tria::FSContactMigration(Vector<IssmDouble>* vertex_sigmann,Vector<IssmDouble>* vertex_waterpressure){/*{{{*/

	if(!IsOnBase()) return;

	int approximation;
	this->GetInput2Value(&approximation,ApproximationEnum);

	if(approximation==HOApproximationEnum || approximation==SSAApproximationEnum || approximation==SSAHOApproximationEnum){
		_error_(" contact contiditon only works for FS elements");
	}
	/*Intermediaries*/
	IssmDouble* xyz_list = NULL;
	IssmDouble  bed_normal[2],base[NUMVERTICES],bed[NUMVERTICES],surface[NUMVERTICES],phi[NUMVERTICES];
	IssmDouble  water_pressure[NUMVERTICES],pressureice[NUMVERTICES],pressure[NUMVERTICES];
	IssmDouble  sigmaxx[NUMVERTICES],sigmayy[NUMVERTICES],sigmaxy[NUMVERTICES],sigma_nn[NUMVERTICES];
	IssmDouble  viscosity,epsilon[NUMVERTICES];
	Element::GetInputListOnVertices(&base[0],BaseEnum);
	Element::GetInputListOnVertices(&bed[0],BedEnum);
	Element::GetInputListOnVertices(&surface[0],SurfaceEnum);
	Element::GetInputListOnVertices(&pressure[0],PressureEnum);
	Element::GetInputListOnVertices(&phi[0],MaskGroundediceLevelsetEnum);
	IssmDouble rho_ice   = FindParam(MaterialsRhoIceEnum);
	IssmDouble rho_water = FindParam(MaterialsRhoSeawaterEnum);
	IssmDouble gravity   = FindParam(ConstantsGEnum);

	/* Get node coordinates and dof list: */
	GetVerticesCoordinates(&xyz_list);
	/*Retrieve all inputs we will be needing: */
	Input2* vx_input       = this->GetInput2(VxEnum);       _assert_(vx_input);
	Input2* vy_input       = this->GetInput2(VyEnum);       _assert_(vy_input);

	/*1. Recover stresses at the base*/
	GaussTria* gauss=new GaussTria();
	for (int iv=0;iv<NUMVERTICES;iv++){
		gauss->GaussVertex(iv);

		/*Compute strain rate viscosity and pressure: */
		this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
		this->material->ViscosityFS(&viscosity,2,xyz_list,gauss,vx_input,vy_input,NULL);
		/*FIXME: this is for Hongju only*/
	//	pressureice[iv]=gravity*rho_ice*(surface[iv]-base[iv]);
	//	if (pressure[iv]/pressureice[iv]>1) pressure[iv]=pressureice[iv];

		/*Compute Stress*/
		sigmaxx[iv]=2*viscosity*epsilon[0]-pressure[iv];
		sigmayy[iv]=2*viscosity*epsilon[1]-pressure[iv];
		sigmaxy[iv]=2*viscosity*epsilon[2];
	}

	/*2. compute contact condition*/
	for(int i=0;i<NUMVERTICES;i++){
		/*If was grounded*/
		if (phi[i]>=0.){
			NormalBase(&bed_normal[0],xyz_list);
			sigma_nn[i]=-1*(sigmaxx[i]*bed_normal[0]*bed_normal[0] + sigmayy[i]*bed_normal[1]*bed_normal[1]+2*sigmaxy[i]*bed_normal[0]*bed_normal[1]);
			water_pressure[i]=-gravity*rho_water*base[i];
			vertex_sigmann->SetValue(vertices[i]->Pid(),sigma_nn[i],ADD_VAL);
			vertex_waterpressure->SetValue(vertices[i]->Pid(),water_pressure[i],ADD_VAL);
		}
		/*If was floating*/
		else{
			/*Tricky part:
			 * 1. if base is now touching, we put 1 for sigma_nn and leave water pressure at 0 so that the rest of the module will reground this vertex
			 * 2. if base is still above bed, water pressure is set as 1, sigma_nn is left as 0, so the GL module will keep it afloat*/
			if(base[i]<bed[i]) vertex_sigmann->SetValue(vertices[i]->Pid(),+1.,ADD_VAL);
			vertex_waterpressure->SetValue(vertices[i]->Pid(),+1.,ADD_VAL);
		}
	}

	/*clean up*/
	delete gauss;
	xDelete<IssmDouble>(xyz_list);
}
/*}}}*/
IssmDouble Tria::GetArea(void){/*{{{*/

	IssmDouble xyz_list[NUMVERTICES][3];
	IssmDouble x1,y1,x2,y2,x3,y3;

	/*Get xyz list: */
	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
	x1=xyz_list[0][0]; y1=xyz_list[0][1];
	x2=xyz_list[1][0]; y2=xyz_list[1][1];
	x3=xyz_list[2][0]; y3=xyz_list[2][1];

	_assert_(x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1>0);
	return (x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1)/2;
}
/*}}}*/
IssmDouble Tria::GetHorizontalSurfaceArea(void){/*{{{*/

	return this->GetArea();
}
/*}}}*/
IssmDouble Tria::GetArea3D(void){/*{{{*/

	IssmDouble xyz_list[NUMVERTICES][3];
	IssmDouble x1,y1,z1,x2,y2,z2,x3,y3,z3;
	IssmDouble detm1,detm2,detm3;

	/*Get xyz list: */
	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
	x1=xyz_list[0][0]; y1=xyz_list[0][1]; z1=xyz_list[0][2];
	x2=xyz_list[1][0]; y2=xyz_list[1][1]; z2=xyz_list[1][2];
	x3=xyz_list[2][0]; y3=xyz_list[2][1]; z3=xyz_list[2][2];

	detm1=x1*y2 - x2*y1 - x1*y3 + x3*y1 + x2*y3 - x3*y2;
	detm2=y1*z2 - y2*z1 - y1*z3 + y3*z1 + y2*z3 - y3*z2;
	detm3=x2*z1 - x1*z2 + x1*z3 - x3*z1 - x2*z3 + x3*z2;

	return sqrt(pow(detm1,2) + pow(detm2,2) + pow(detm3,2))/2;
}
/*}}}*/
IssmDouble Tria::GetAreaIce(void){/*{{{*/

	/*return area of element covered by ice*/
	/*Intermediaries*/
	int numiceverts;
	IssmDouble area_fraction;
	IssmDouble s[2]; // s:fraction of intersected triangle edges that lie inside ice
	int* indices=NULL;

	this->GetLevelsetIntersection(&indices, &numiceverts, s, MaskIceLevelsetEnum, 0.);

	switch (numiceverts){
		case 0: // no vertex has ice: element is ice free
			area_fraction=0.;
			break;
		case 1: // one vertex has ice: get area of triangle
			area_fraction=s[0]*s[1];
			break;
		case 2: // two vertices have ice: get area of quadrangle
			area_fraction=s[0]+s[1]-s[0]*s[1];
			break;
		case NUMVERTICES: // all vertices have ice: return triangle area
			area_fraction=1.;
			break;
		default:
			_error_("Wrong number of ice vertices in Tria::GetAreaIce!");
			break;
	}
	_assert_((area_fraction>=0.) && (area_fraction<=1.));

	xDelete<int>(indices);
	return area_fraction*this->GetArea();
}/*}}}*/
IssmDouble Tria::GetAreaSpherical(void){/*{{{*/

	bool spherical=true;
	IssmDouble llr_list[NUMVERTICES][3];
	IssmDouble x1,y1,z1,x2,y2,z2,x3,y3,z3;
	IssmDouble arc12,arc23,arc31,semi_peri,excess;

	/*retrieve coordinates: lat,long,radius */
	::GetVerticesCoordinates(&llr_list[0][0],vertices,NUMVERTICES,spherical);
	x1=llr_list[0][0]/180*PI; y1=llr_list[0][1]/180*PI; z1=llr_list[0][2];
	x2=llr_list[1][0]/180*PI; y2=llr_list[1][1]/180*PI; z2=llr_list[1][2];
	x3=llr_list[2][0]/180*PI; y3=llr_list[2][1]/180*PI; z3=llr_list[2][2];

	/*compute great circle distance between vertices */
	arc12=2.*asin(sqrt(pow(sin((x2-x1)/2),2.0)+cos(x1)*cos(x2)*pow(sin((y2-y1)/2),2)));
	arc23=2.*asin(sqrt(pow(sin((x3-x2)/2),2.0)+cos(x2)*cos(x3)*pow(sin((y3-y2)/2),2)));
	arc31=2.*asin(sqrt(pow(sin((x1-x3)/2),2.0)+cos(x3)*cos(x1)*pow(sin((y1-y3)/2),2)));

	/*semi parameter */
	semi_peri=(arc12+arc23+arc31)/2;

	/*spherical excess */
	excess=4.*atan(sqrt(tan(semi_peri/2)*tan((semi_peri-arc12)/2)*tan((semi_peri-arc23)/2)*tan((semi_peri-arc31)/2)));

	/*area = excess*radius^2 */
	return excess*pow((z1+z2+z3)/3,2);
}
/*}}}*/
void       Tria::GetAreaCoordinates(IssmDouble* area_coordinates,IssmDouble* xyz_zero,IssmDouble* xyz_list,int numpoints){/*{{{*/
	/*Computeportion of the element that is grounded*/

	int         i,j,k;
	IssmDouble  area_init,area_portion;
	IssmDouble  xyz_bis[NUMVERTICES][3];

	area_init=GetArea();

	/*Initialize xyz_list with original xyz_list of triangle coordinates*/
	for(j=0;j<3;j++){
		for(k=0;k<3;k++){
			xyz_bis[j][k]=xyz_list[j*3+k];
		}
	}
	for(i=0;i<numpoints;i++){
		for(j=0;j<3;j++){
			for(k=0;k<3;k++){
				/*Change appropriate line*/
				xyz_bis[j][k]=xyz_zero[i*3+k];
			}

			/*Compute area fraction*/
			area_portion=fabs(xyz_bis[1][0]*xyz_bis[2][1] - xyz_bis[1][1]*xyz_bis[2][0] + xyz_bis[0][0]*xyz_bis[1][1] - xyz_bis[0][1]*xyz_bis[1][0] + xyz_bis[2][0]*xyz_bis[0][1] - xyz_bis[2][1]*xyz_bis[0][0])/2.;
			*(area_coordinates+3*i+j)=area_portion/area_init;

			/*Reinitialize xyz_list*/
			for(k=0;k<3;k++){
				/*Reinitialize xyz_list with original coordinates*/
				xyz_bis[j][k]=xyz_list[j*3+k];
			}
		}
	}
}
/*}}}*/
int        Tria::GetElementType(){/*{{{*/

	/*return TriaRef field*/
	return this->element_type;

}
/*}}}*/
void       Tria::GetGroundedPart(int* point1,IssmDouble* fraction1,IssmDouble* fraction2, bool* mainlyfloating){/*{{{*/
	/*Computeportion of the element that is grounded*/

	bool               floating=true;
	int                point;
	const IssmPDouble  epsilon= 1.e-15;
	IssmDouble         gl[NUMVERTICES];
	IssmDouble         f1,f2;

	/*Recover parameters and values*/
	Element::GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum);

	/*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 grounded or floating*/
	if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All grounded
		point=0;
		f1=1.;
		f2=1.;
	}
	else if(gl[0]<0 && gl[1]<0 && gl[2]<0){ //All floating
		point=0;
		f1=0.;
		f2=0.;
	}
	else{
		if(gl[0]*gl[1]*gl[2]<0) floating=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]);
		}
		else _error_("case not possible");
	}
	*point1=point;
	*fraction1=f1;
	*fraction2=f2;
	*mainlyfloating=floating;
}
/*}}}*/
IssmDouble Tria::GetGroundedPortion(IssmDouble* xyz_list){/*{{{*/
	/*Computeportion of the element that is grounded*/

	bool              mainlyfloating = true;
	int               domaintype,index1,index2;
	const IssmPDouble epsilon        = 1.e-15;
	IssmDouble        phi,s1,s2,area_init,area_grounded;
	IssmDouble        gl[NUMVERTICES];
	IssmDouble        xyz_bis[3][3];

	/*Recover parameters and values*/
	parameters->FindParam(&domaintype,DomainTypeEnum);
	Element::GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum);

	/*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;

	if(domaintype==Domain2DverticalEnum){
		this->EdgeOnBaseIndices(&index1,&index2);
		if(gl[index1]>0 && gl[index2]>0) phi=1; // All grounded
		else if(gl[index1]<0 && gl[index2]<0) phi=0; // All floating
		else if(gl[index1]<0 && gl[index2]>0){ //index2 grounded
			phi=1./(1.-gl[index1]/gl[index2]);
		}
		else if(gl[index2]<0 && gl[index1]>0){ //index1 grounded
			phi=1./(1.-gl[index2]/gl[index1]);
		}

	}
	else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum || domaintype==Domain3DsurfaceEnum){
		/*Check that not all nodes are grounded or floating*/
		if(gl[0]>0 && gl[1]>0 && gl[2]>0){ // All grounded
			phi=1;
		}
		else if(gl[0]<0 && gl[1]<0 && gl[2]<0){ //All floating
			phi=0;
		}
		else{
			/*Figure out if two nodes are floating or grounded*/
			if(gl[0]*gl[1]*gl[2]>0) mainlyfloating=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
				/*Coordinates of point 2: same as initial point 2*/
				xyz_bis[2][0]=xyz_list[3*2+0];
				xyz_bis[2][1]=xyz_list[3*2+1];
				xyz_bis[2][2]=xyz_list[3*2+2];

				/*Portion of the segments*/
				s1=gl[2]/(gl[2]-gl[1]);
				s2=gl[2]/(gl[2]-gl[0]);

				/*New point 1*/
				xyz_bis[1][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
				xyz_bis[1][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
				xyz_bis[1][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);

				/*New point 0*/
				xyz_bis[0][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
				xyz_bis[0][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
				xyz_bis[0][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
			}
			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
				/*Coordinates of point 0: same as initial point 2*/
				xyz_bis[0][0]=xyz_list[3*0+0];
				xyz_bis[0][1]=xyz_list[3*0+1];
				xyz_bis[0][2]=xyz_list[3*0+2];

				/*Portion of the segments*/
				s1=gl[0]/(gl[0]-gl[1]);
				s2=gl[0]/(gl[0]-gl[2]);

				/*New point 1*/
				xyz_bis[1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
				xyz_bis[1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
				xyz_bis[1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);

				/*New point 2*/
				xyz_bis[2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
				xyz_bis[2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
				xyz_bis[2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+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
				/*Coordinates of point 1: same as initial point 2*/
				xyz_bis[1][0]=xyz_list[3*1+0];
				xyz_bis[1][1]=xyz_list[3*1+1];
				xyz_bis[1][2]=xyz_list[3*1+2];

				/*Portion of the segments*/
				s1=gl[1]/(gl[1]-gl[0]);
				s2=gl[1]/(gl[1]-gl[2]);

				/*New point 0*/
				xyz_bis[0][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
				xyz_bis[0][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
				xyz_bis[0][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);

				/*New point 2*/
				xyz_bis[2][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
				xyz_bis[2][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
				xyz_bis[2][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
			}
			else _error_("case not possible");

			/*Compute fraction of grounded element*/
			if(domaintype==Domain3DsurfaceEnum){ //hack, need to be implemented in a Tria 3D
				GetJacobianDeterminant3D(&area_init, xyz_list,NULL);
				GetJacobianDeterminant3D(&area_grounded, &xyz_bis[0][0],NULL);
			}
			else{
				GetJacobianDeterminant(&area_init, xyz_list,NULL);
				GetJacobianDeterminant(&area_grounded, &xyz_bis[0][0],NULL);
			}
			if(mainlyfloating==true) area_grounded=area_init-area_grounded;
			phi=area_grounded/area_init;
		}
	}
	else _error_("mesh type "<<EnumToStringx(domaintype)<<"not supported yet ");

	if(phi>1 || phi<0) _error_("Error. Problem with portion of grounded element: value should be between 0 and 1");

	return phi;
}
/*}}}*/
IssmDouble Tria::GetIcefrontArea(){/*{{{*/

	IssmDouble  bed[NUMVERTICES]; //basinId[NUMVERTICES];
	IssmDouble	Haverage,frontarea;
	IssmDouble  x1,y1,x2,y2,distance;
	IssmDouble lsf[NUMVERTICES], Haux[NUMVERTICES], surfaces[NUMVERTICES], bases[NUMVERTICES];
	int* indices=NULL;
	IssmDouble* H=NULL;;
	int nrfrontbed,numiceverts;

	if(!IsZeroLevelset(MaskIceLevelsetEnum)) return 0;

	/*Retrieve all inputs and parameters*/
	Element::GetInputListOnVertices(&bed[0],BedEnum);
	Element::GetInputListOnVertices(&surfaces[0],SurfaceEnum);
	Element::GetInputListOnVertices(&bases[0],BaseEnum);
	Element::GetInputListOnVertices(&lsf[0],MaskIceLevelsetEnum);

	nrfrontbed=0;
	for(int i=0;i<NUMVERTICES;i++){
		/*Find if bed<0*/
		if(bed[i]<0.) nrfrontbed++;
	}

	if(nrfrontbed==3){
		/*2. Find coordinates of where levelset crosses 0*/
		int         numiceverts;
		IssmDouble  s[2],x[2],y[2];
		this->GetLevelsetIntersection(&indices, &numiceverts,&s[0],MaskIceLevelsetEnum,0.);
		_assert_(numiceverts);

		/*3 Write coordinates*/
		IssmDouble  xyz_list[NUMVERTICES][3];
		::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
		int counter = 0;
		if((numiceverts>0) && (numiceverts<NUMVERTICES)){
			for(int i=0;i<numiceverts;i++){
				for(int n=numiceverts;n<NUMVERTICES;n++){ // iterate over no-ice vertices
					x[counter] = xyz_list[indices[i]][0]+s[counter]*(xyz_list[indices[n]][0]-xyz_list[indices[i]][0]);
					y[counter] = xyz_list[indices[i]][1]+s[counter]*(xyz_list[indices[n]][1]-xyz_list[indices[i]][1]);
					counter++;
				}
			}
		}
		else if(numiceverts==NUMVERTICES){ //NUMVERTICES ice vertices: calving front lies on element edge

			for(int i=0;i<NUMVERTICES;i++){
				if(lsf[indices[i]]==0.){
					x[counter]=xyz_list[indices[i]][0];
					y[counter]=xyz_list[indices[i]][1];
					counter++;
				}
				if(counter==2) break;
			}
			if(counter==1){
				/*We actually have only 1 vertex on levelset, write a single point as a segment*/
				x[counter]=x[0];
				y[counter]=y[0];
				counter++;
			}
		}
		else{
			_error_("not sure what's going on here...");
		}
		x1=x[0]; y1=y[0]; x2=x[1]; y2=y[1];
		distance=sqrt(pow((x1-x2),2)+pow((y1-y2),2));

		int numthk=numiceverts+2;
		H=xNew<IssmDouble>(numthk);
		for(int iv=0;iv<NUMVERTICES;iv++) Haux[iv]=-bed[indices[iv]]; //sort bed in ice/noice

		switch(numiceverts){
			case 1: // average over triangle
				H[0]=Haux[0];
				H[1]=Haux[0]+s[0]*(Haux[1]-Haux[0]);
				H[2]=Haux[0]+s[1]*(Haux[2]-Haux[0]);
				Haverage=(H[1]+H[2])/2;
				break;
			case 2: // average over quadrangle
				H[0]=Haux[0];
				H[1]=Haux[1];
				H[2]=Haux[0]+s[0]*(Haux[2]-Haux[0]);
				H[3]=Haux[1]+s[1]*(Haux[2]-Haux[1]);
				Haverage=(H[2]+H[3])/2;
				break;
			default:
				_error_("Number of ice covered vertices wrong in Tria::GetIceFrontArea(void)");
				break;
		}
		frontarea=distance*Haverage;
	}
	else return 0;

	xDelete<int>(indices);
	xDelete<IssmDouble>(H);

	_assert_(frontarea>0);
	return frontarea;
}
/*}}}*/
void       Tria::GetIcefrontCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum){/*{{{*/

	/* Intermediaries */
	IssmDouble  levelset[NUMVERTICES];
	int         indicesfront[NUMVERTICES];

	/*Recover parameters and values*/
	Element::GetInputListOnVertices(&levelset[0],levelsetenum);

	/* Get nodes where there is no ice */
	int num_frontnodes=0;
	for(int i=0;i<NUMVERTICES;i++){
		if(levelset[i]>=0.){
			indicesfront[num_frontnodes]=i;
			num_frontnodes++;
		}
	}
	_assert_(num_frontnodes==2);

	/* arrange order of frontnodes such that they are oriented counterclockwise */
	if((NUMVERTICES+indicesfront[0]-indicesfront[1])%NUMVERTICES!=NUMVERTICES-1){
		int index=indicesfront[0];
		indicesfront[0]=indicesfront[1];
		indicesfront[1]=index;
	}

	IssmDouble* xyz_front = xNew<IssmDouble>(3*2);
	/* Return nodes */
	for(int i=0;i<2;i++){
		for(int j=0;j<3;j++){
			xyz_front[3*i+j]=xyz_list[3*indicesfront[i]+j];
		}
	}

	*pxyz_front=xyz_front;
}/*}}}*/
Input2*    Tria::GetInput2(int inputenum){/*{{{*/

	/*Get Input from dataset*/
	if(this->iscollapsed){
		PentaInput2* input = this->inputs2->GetPentaInput(inputenum);
		if(!input) return input;

		this->InputServe(input);
		return input;
	}
	else{
		TriaInput2* input = this->inputs2->GetTriaInput(inputenum);
		if(!input) return input;

		this->InputServe(input);
		return input;
	}
}/*}}}*/
Input2*    Tria::GetInput2(int inputenum,IssmDouble time){/*{{{*/

	/*Get Input from dataset*/
	if(this->iscollapsed){
		PentaInput2* input = this->inputs2->GetPentaInput(inputenum,time);
		if(!input) return input;

		this->InputServe(input);
		return input;
	}
	else{
		TriaInput2* input = this->inputs2->GetTriaInput(inputenum,time);
		if(!input) return input;

		this->InputServe(input);
		return input;
	}
}/*}}}*/
Input2*    Tria::GetInput2(int inputenum,IssmDouble start_time, IssmDouble end_time){/*{{{*/

	/*Get Input from dataset*/
	if(this->iscollapsed){
		_error_("Get Average input not implemented in Penta yet");
	}
	else{
		TriaInput2* input = this->inputs2->GetTriaInput(inputenum,start_time, end_time);
		if(!input) return input;

		this->InputServe(input);
		return input;
	}
}/*}}}*/
void       Tria::GetInputListOnVertices(IssmDouble* pvalue,Input2* input,IssmDouble default_value){/*{{{*/

	/*Checks in debugging mode*/
	_assert_(pvalue);

	/* Start looping on the number of vertices: */
	if(input){
		GaussTria gauss;
		for(int iv=0;iv<NUMVERTICES;iv++){
			gauss.GaussVertex(iv);
			input->GetInputValue(&pvalue[iv],&gauss);
		}
	}
	else{
		for(int iv=0;iv<NUMVERTICES;iv++) pvalue[iv] = default_value;
	}
}
/*}}}*/
void       Tria::GetInputListOnNodes(IssmDouble* pvalue,Input2* input,IssmDouble default_value){/*{{{*/

	/*Checks in debugging mode*/
	_assert_(pvalue);

	/*What type of finite element are we dealing with?*/
	int fe       = this->FiniteElement();
	int numnodes = this->GetNumberOfNodes();

	/* Start looping on the number of vertices: */
	if(input){
		GaussTria gauss;
		for(int iv=0;iv<numnodes;iv++){
			gauss.GaussNode(fe,iv);
			input->GetInputValue(&pvalue[iv],&gauss);
		}
	}
	else{
		for(int iv=0;iv<numnodes;iv++) pvalue[iv] = default_value;
	}
}
/*}}}*/
void       Tria::InputServe(Input2* input_in){/*{{{*/

	/*Return NULL pointer if input is NULL*/
	if(!input_in) return;

	/*Get Input from dataset*/
	if(this->iscollapsed){
		_assert_(input_in->ObjectEnum()==PentaInput2Enum);
		PentaInput2* input = xDynamicCast<PentaInput2*>(input_in);

		/*Intermediaries*/
		int numindices;
		int indices[3];

		/*Check interpolation*/
		int interpolation = input->GetInterpolation();
		switch(interpolation){
			case P0Enum:
				numindices = 1;
				indices[0] = this->lid;
				input->Serve(numindices,&indices[0]);
				break;
			case P1Enum:
				numindices = 3;
				for(int i=0;i<3;i++) indices[i] = vertices[i]->lid;
				input->Serve(numindices,&indices[0]);
				break;
			case P1DGEnum:
			case P1bubbleEnum:
				input->ServeCollapsed(this->lid,this->iscollapsed);
				break;
			default: _error_("interpolation "<<EnumToStringx(interpolation)<<" not supported");
		}

		/*Flag as collapsed for later use*/
		input->SetServeCollapsed(true);
		return;
	}
	else{
		_assert_(input_in->ObjectEnum()==TriaInput2Enum);
		TriaInput2* input = xDynamicCast<TriaInput2*>(input_in);

		/*Intermediaries*/
		int numindices;
		int indices[7];

		/*Check interpolation*/
		int interpolation = input->GetInterpolation();
		switch(interpolation){
			case P0Enum:
				numindices = 1;
				indices[0] = this->lid;
				input->Serve(numindices,&indices[0]);
				break;
			case P1Enum:
				numindices = 3;
				for(int i=0;i<3;i++) indices[i] = vertices[i]->lid;
				input->Serve(numindices,&indices[0]);
				break;
			case P1DGEnum:
				numindices = 3;
				input->Serve(this->lid,numindices);
				break;
			default:
				input->Serve(this->lid,this->GetNumberOfNodes(interpolation));
				break;
		}
		return;
	}
}/*}}}*/
DatasetInput2* Tria::GetDatasetInput2(int inputenum){/*{{{*/

	DatasetInput2* datasetinput = this->inputs2->GetDatasetInput2(inputenum);
	if(!datasetinput) return NULL;

	for(int i=0;i<datasetinput->GetNumIds();i++){

		/*Get Input from dataset*/
		if(this->iscollapsed){

			PentaInput2* input = datasetinput->GetPentaInputByOffset(i); _assert_(input);

			/*Intermediaries*/
			int numindices;
			int indices[3];

			/*Check interpolation*/
			int interpolation = input->GetInterpolation();
			switch(interpolation){
				case P0Enum:
					numindices = 1;
					indices[0] = this->lid;
					input->Serve(numindices,&indices[0]);
					break;
				case P1Enum:
					numindices = 3;
					for(int i=0;i<3;i++) indices[i] = vertices[i]->lid;
					input->Serve(numindices,&indices[0]);
					break;
				case P1DGEnum:
				case P1bubbleEnum:
					input->ServeCollapsed(this->lid,this->iscollapsed);
					break;
				default: _error_("interpolation "<<EnumToStringx(interpolation)<<" not supported");
			}

			/*Flag as collapsed for later use*/
			input->SetServeCollapsed(true);
		}
		else{

			TriaInput2* input = datasetinput->GetTriaInputByOffset(i); _assert_(input);

			/*Intermediaries*/
			int numindices;
			int indices[7];

			/*Check interpolation*/
			int interpolation = input->GetInterpolation();
			switch(interpolation){
				case P0Enum:
					numindices = 1;
					indices[0] = this->lid;
					input->Serve(numindices,&indices[0]);
					break;
				case P1Enum:
					numindices = 3;
					for(int i=0;i<3;i++) indices[i] = vertices[i]->lid;
					input->Serve(numindices,&indices[0]);
					break;
				case P1DGEnum:
					numindices = 3;
					input->Serve(this->lid,numindices);
					break;
				default: _error_("interpolation "<<EnumToStringx(interpolation)<<" not supported");
			}

		}
	}

	return datasetinput;
}/*}}}*/
void       Tria::CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble start_time,IssmDouble end_time){/*{{{*/

	_assert_(end_time>start_time);

	/*Intermediaries*/
	IssmDouble averaged_values[NUMVERTICES];
	IssmDouble current_values[NUMVERTICES];
	IssmDouble dt;
	int        found,start_offset,end_offset;
	int        averaging_method=0;


	/*Get transient input time steps*/
	int         numtimesteps;
	IssmDouble *timesteps    = NULL;
	TransientInput2* transient_input  = this->inputs2->GetTransientInput(transientinput_enum);
	transient_input->GetAllTimes(&timesteps,&numtimesteps);

	/*go through the timesteps, and grab offset for start and end*/
	found=binary_search(&start_offset,start_time,timesteps,numtimesteps);
	if(!found) _error_("Input not found (is TransientInput sorted ?)");
	found=binary_search(&end_offset,end_time,timesteps,numtimesteps);
	if(!found) _error_("Input not found (is TransientInput sorted ?)");

	Gauss* gauss=this->NewGauss();

	/*stack the input for each timestep in the slice*/
	int offset = start_offset;
	while(offset < end_offset ){
		if(offset==-1){
			/*get values for the first time: */
			_assert_(start_time<timesteps[0]);
			TriaInput2* input = transient_input->GetTriaInput(0);
			_assert_(input->GetInterpolation()==P1Enum);
			this->InputServe(input);
			for(int iv=0;iv<NUMVERTICES;iv++){
				gauss->GaussVertex(iv);
				input->GetInputValue(&current_values[iv],gauss);
			}
			dt = timesteps[0] - start_time; _assert_(dt>0.);
		}
		else{
			TriaInput2* input = transient_input->GetTriaInput(offset+1);
			_assert_(input->GetInterpolation()==P1Enum);
			this->InputServe(input);
			for(int iv=0;iv<NUMVERTICES;iv++){
				gauss->GaussVertex(iv);
				input->GetInputValue(&current_values[iv],gauss);
			}
			if(offset == numtimesteps-1){
				dt = end_time - timesteps[offset]; _assert_(dt>0.);
			}
			else{
				dt = timesteps[offset+1] - timesteps[offset]; _assert_(dt>0.);
			}
		}

		switch(averaging_method){
			case 0: /*Arithmetic mean*/
				if(offset==start_offset){
					for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv]  = dt*current_values[iv];
				}
				else{
					for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] += dt*current_values[iv];
				}
				break;
			case 1: /*Geometric mean*/
				if(offset==start_offset){
					for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv]  = dt*current_values[iv];
				}
				else{
					for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] *= dt*current_values[iv];
				}
				break;
			case 2: /*Harmonic mean*/
				if(offset==start_offset){
					for(int iv=0;iv<NUMVERTICES;iv++){
						_assert_(current_values[iv]>1.e-50);
						averaged_values[iv]  = dt*1./current_values[iv];
					}
				}
				else{
					for(int iv=0;iv<NUMVERTICES;iv++){
						_assert_(current_values[iv]>1.e-50);
						averaged_values[iv]  += dt*1./current_values[iv];
					}
				}
				break;
			default:
				_error_("averaging method is not recognised");
		}

		offset+=1;
	}

	/*Integration done, now normalize*/
	switch(averaging_method){
		case 0: //Arithmetic mean
			for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] =  averaged_values[iv]/(end_time - start_time);
			break;
		case 1: /*Geometric mean*/
			for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] = pow(averaged_values[iv], 1./(end_time - start_time));
			break;
		case 2: /*Harmonic mean*/
			for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] = 1./(averaged_values[iv]/(end_time - start_time));
			break;
		default:
			_error_("averaging method is not recognised");
	}

	this->AddInput2(averagedinput_enum,&averaged_values[0],P1Enum);

	/*Cleanup*/
	delete gauss;
	xDelete<IssmDouble>(timesteps);
}
/*}}}*/
void       Tria::GetInputAveragesUpToCurrentTime(int input_enum,IssmDouble** pvalues, IssmDouble** ptimes, int* pnumtimes, IssmDouble currenttime){/*{{{*/

	/*Get transient input time steps*/
	int         numtimesteps;
	IssmDouble *timesteps    = NULL;
	TransientInput2* transient_input  = this->inputs2->GetTransientInput(input_enum);

	transient_input->GetAllTimes(&timesteps,&numtimesteps);

	/*Figure out how many time steps we are going to return: */
	int  numsteps               = 0;
	bool iscurrenttime_included = false;
	for(int i=0;i<numtimesteps;i++){
		if(timesteps[i]==currenttime) iscurrenttime_included=true;
		if(timesteps[i]>currenttime)  break;
		else numsteps++;
	}
	if(iscurrenttime_included==false)numsteps++;

	/*allocate: */
	IssmDouble* times=xNew<IssmDouble>(numsteps);
	IssmDouble* values=xNew<IssmDouble>(numsteps);

	for(int i=0;i<numsteps;i++){
		if((iscurrenttime_included==false) && (i==(numsteps-1))){
			Input2* input = this->GetInput2(input_enum,currenttime);
			input->GetInputAverage(&values[i]);
			times[i]=currenttime;
		}
		else{
			TriaInput2* input = transient_input->GetTriaInput(i);
			this->InputServe(input);
			input->GetInputAverage(&values[i]);
			times[i]=timesteps[i];
		}
	}

	/*Assign output pointers*/
	*pvalues=values;
	*ptimes=times;
	*pnumtimes=numtimesteps;
}
/*}}}*/
void       Tria::GetInputValue(IssmDouble* pvalue,Node* node,int enumtype){/*{{{*/

	Input2* input=this->GetInput2(enumtype);
	if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria");

	int index = this->GetNodeIndex(node);

	GaussTria* gauss=new GaussTria();
	gauss->GaussNode(this->element_type,index);

	input->GetInputValue(pvalue,gauss);
	delete gauss;
}
/*}}}*/
void       Tria::GetInputValue(IssmDouble* pvalue,Vertex* vertex,int enumtype){/*{{{*/

	Input2* input=this->GetInput2(enumtype);
	if(!input) _error_("No input of type " << EnumToStringx(enumtype) << " found in tria");

	int index = this->GetVertexIndex(vertex);

	GaussTria* gauss=new GaussTria();
	gauss->GaussVertex(index);

	input->GetInputValue(pvalue,gauss);
	delete gauss;
}
/*}}}*/
void       Tria::GetLevelCoordinates(IssmDouble** pxyz_front,IssmDouble* xyz_list,int levelsetenum,IssmDouble level){/*{{{*/

	/* Intermediaries */
	int i, dir,nrfrontnodes;
	IssmDouble  levelset[NUMVERTICES];
	int indicesfront[NUMVERTICES];

	/*Recover parameters and values*/
	Element::GetInputListOnVertices(&levelset[0],levelsetenum);

	/* Get nodes where there is no ice */
	nrfrontnodes=0;
	for(i=0;i<NUMVERTICES;i++){
		if(levelset[i]==level){
			indicesfront[nrfrontnodes]=i;
			nrfrontnodes++;
		}
	}

	_assert_(nrfrontnodes==2);

	/* arrange order of frontnodes such that they are oriented counterclockwise */
	if((NUMVERTICES+indicesfront[0]-indicesfront[1])%NUMVERTICES!=NUMVERTICES-1){
		int index=indicesfront[0];
		indicesfront[0]=indicesfront[1];
		indicesfront[1]=index;
	}

	IssmDouble* xyz_front = xNew<IssmDouble>(3*nrfrontnodes);
	/* Return nodes */
	for(i=0;i<nrfrontnodes;i++){
		for(dir=0;dir<3;dir++){
			xyz_front[3*i+dir]=xyz_list[3*indicesfront[i]+dir];
		}
	}

	*pxyz_front=xyz_front;

}/*}}}*/
void       Tria::GetLevelsetIntersection(int** pindices, int* pnumiceverts, IssmDouble* fraction, int levelset_enum, IssmDouble level){/*{{{*/

	/* GetLevelsetIntersection computes:
	 * 1. indices of element, sorted in [iceverts, noiceverts] in counterclockwise fashion,
	 * 2. fraction of intersected triangle edges intersected by levelset, lying below level*/

	/*Intermediaries*/
	int i, numiceverts, numnoiceverts;
	int ind0, ind1, lastindex;
	int indices_ice[NUMVERTICES],indices_noice[NUMVERTICES];
	IssmDouble lsf[NUMVERTICES];
	int* indices = xNew<int>(NUMVERTICES);

	/*Retrieve all inputs and parameters*/
	Element::GetInputListOnVertices(&lsf[0],levelset_enum);

	/* Determine distribution of ice over element.
	 * Exploit: ice/no-ice parts are connected, so find starting vertex of segment*/
	lastindex=0;
	for(i=0;i<NUMVERTICES;i++){ // go backwards along vertices, and check for sign change
		ind0=(NUMVERTICES-i)%NUMVERTICES;
		ind1=(ind0-1+NUMVERTICES)%NUMVERTICES;
		if((lsf[ind0]-level)*(lsf[ind1]-level)<=0.){ // levelset has been crossed, find last index belonging to segment
			if(lsf[ind1]==level) //if levelset intersects 2nd vertex, choose this vertex as last
				lastindex=ind1;
			else
				lastindex=ind0;
			break;
		}
	}

	numiceverts=0;
	numnoiceverts=0;
	for(i=0;i<NUMVERTICES;i++){
		ind0=(lastindex+i)%NUMVERTICES;
		if(lsf[i]<=level){
			indices_ice[numiceverts]=i;
			numiceverts++;
		}
		else{
			indices_noice[numnoiceverts]=i;
			numnoiceverts++;
		}
	}
	//merge indices
	for(i=0;i<numiceverts;i++){indices[i]=indices_ice[i];}
	for(i=0;i<numnoiceverts;i++){indices[numiceverts+i]=indices_noice[i];}

	switch (numiceverts){
		case 0: // no vertex has ice: element is ice free, no intersection
			for(i=0;i<2;i++)
				fraction[i]=0.;
			break;
		case 1: // one vertex has ice:
			for(i=0;i<2;i++){
				fraction[i]=(level-lsf[indices[0]])/(lsf[indices[numiceverts+i]]-lsf[indices[0]]);
			}
			break;
		case 2: // two vertices have ice: fraction is computed from first ice vertex to last in CCW fashion
			for(i=0;i<2;i++){
				fraction[i]=(level-lsf[indices[i]])/(lsf[indices[numiceverts]]-lsf[indices[i]]);
			}
			break;
		case NUMVERTICES: // all vertices have ice: return triangle area
			for(i=0;i<2;i++)
				fraction[i]=1.;
			break;
		default:
			_error_("Wrong number of ice vertices in Tria::GetLevelsetIntersection!");
			break;
	}

	*pindices=indices;
	*pnumiceverts=numiceverts;
}
/*}}}*/
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=0;
	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]);
		}
		else{
			_error_("This case should NOT be happening");
		}
	}
	*point1=point;
	*fraction1=f1;
	*fraction2=f2;
	*mainlynegative=negative;
}
/*}}}*/
int        Tria::GetVertexIndex(Vertex* vertex){/*{{{*/
	_assert_(vertices);
	for(int i=0;i<NUMVERTICES;i++){
		if(vertex==vertices[i])
		 return i;
	}
	_error_("Vertex provided not found among element vertices");
}
/*}}}*/
int        Tria::GetNumberOfNodes(void){/*{{{*/
	if (this->nodes) return this->NumberofNodes(this->element_type);
	else return 0;
}
/*}}}*/
int        Tria::GetNumberOfNodes(int enum_type){/*{{{*/
	return this->NumberofNodes(enum_type);
}
/*}}}*/
int        Tria::GetNumberOfVertices(void){/*{{{*/
	return NUMVERTICES;
}
/*}}}*/
void       Tria::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,const char* data){/*{{{*/

	/*Get out if this is not an element input*/
	if(!IsInputEnum(control_enum)) _error_("Enum "<<EnumToStringx(control_enum)<<" is not in IsInput");

	/*Prepare index list*/
	int idlist[NUMVERTICES];
	GradientIndexing(&idlist[0],control_index);

	/*Get input (either in element or material)*/
	ElementInput2* input=this->inputs2->GetControlInput2Data(control_enum,data);   _assert_(input);

	/*Intermediaries*/
	int numindices;
	int indices[NUMVERTICES];

	/*Check interpolation*/
	int interpolation = input->GetInterpolation();
	switch(interpolation){
		case P1Enum:
			numindices = NUMVERTICES;
			for(int i=0;i<NUMVERTICES;i++) indices[i] = vertices[i]->lid;
			input->Serve(numindices,&indices[0]);
			break;
		default: _error_("interpolation "<<EnumToStringx(interpolation)<<" not supported");
	}

	/*Flag as collapsed for later use*/
	if(this->iscollapsed){
		xDynamicCast<PentaInput2*>(input)->SetServeCollapsed(true);
	}

	/* Start looping on the number of vertices: */
	IssmDouble values[NUMVERTICES];
	Gauss*gauss=this->NewGauss();
	for(int iv=0;iv<NUMVERTICES;iv++){
		gauss->GaussVertex(iv);
		input->GetInputValue(&values[iv],gauss);
	}
	delete gauss;

	vector->SetValues(NUMVERTICES,idlist,&values[0],INS_VAL);
}
/*}}}*/
void       Tria::GetVectorFromControlInputs(Vector<IssmDouble>* vector,int control_enum,int control_index,const char* data,int offset){/*{{{*/

	/*Get input*/
	ElementInput2* input=this->inputs2->GetControlInput2Data(control_enum,data);   _assert_(input);

	/*Lid list once for all*/
	int lidlist[NUMVERTICES];
	for(int i=0;i<NUMVERTICES;i++) lidlist[i] = vertices[i]->lid;

	/*Check what input we are dealing with*/
	switch(input->ObjectEnum()){
		case TriaInput2Enum:
			  {
				IssmDouble values[NUMVERTICES];
				int        idlist[NUMVERTICES];

				TriaInput2* triainput = xDynamicCast<TriaInput2*>(input);
				if(triainput->GetInputInterpolationType()!=P1Enum) _error_("not supported yet");
				input->Serve(NUMVERTICES,&lidlist[0]);

				/*Create list of indices and values for global vector*/
				GradientIndexing(&idlist[0],control_index);
				for(int i=0;i<NUMVERTICES;i++) values[i] = triainput->element_values[i];
				vector->SetValues(NUMVERTICES,idlist,values,INS_VAL);
				break;
			  }

		case TransientInputEnum:
				{
					TransientInput2* transientinput = xDynamicCast<TransientInput2*>(input);
					int N = transientinput->numtimesteps;
					int* M=NULL;
					parameters->FindParam(&M,NULL,ControlInputSizeMEnum);
					int* idlist = xNew<int>(NUMVERTICES*N);
					IssmDouble* values = xNew<IssmDouble>(NUMVERTICES*N);
					for(int t=0;t<transientinput->numtimesteps;t++) {
						IssmDouble time = transientinput->GetTimeByOffset(t);
						_error_("not implemented");
						//TriaInput* timeinput = xDynamicCast<TriaInput*>(transientinput->GetTimeInput(time));
						//if(timeinput->interpolation_type!=P1Enum) _error_("not supported yet");
						//input->Serve(NUMVERTICES,&lidlist[0]);
						///*Create list of indices and values for global vector*/
						//for(int i=0;i<NUMVERTICES;i++){
						//		idlist[N*i+t] = offset + this->vertices[i]->Sid()+t*M[control_index];
						//		values[N*i+t] = timeinput->values[i];
						//}
					}
					vector->SetValues(NUMVERTICES*transientinput->numtimesteps,idlist,values,INS_VAL);
					xDelete<int>(M);
					xDelete<int>(idlist);
					xDelete<IssmDouble>(values);
					break;
				}
		default: _error_("input "<<EnumToStringx(input->ObjectEnum())<<" not supported yet");
	}
}
/*}}}*/
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;

}/*}}}*/
IssmDouble Tria::GroundedArea(bool scaled){/*{{{*/

	/*Intermediaries*/
	int         domaintype;
	IssmDouble  phi,scalefactor,groundedarea;
	IssmDouble *xyz_list  = NULL;

	if(!IsIceInElement())return 0.;

	/*Get problem dimension*/
	this->FindParam(&domaintype,DomainTypeEnum);
	if(domaintype!=Domain2DhorizontalEnum && domaintype!=Domain3DEnum) _error_("mesh "<<EnumToStringx(domaintype)<<" not supported yet");

	this->GetVerticesCoordinates(&xyz_list);
	phi=this->GetGroundedPortion(xyz_list);
	groundedarea=phi*this->GetArea();
	if(scaled==true){
		Input2* scalefactor_input = this->GetInput2(MeshScaleFactorEnum); _assert_(scalefactor_input);
		scalefactor_input->GetInputAverage(&scalefactor);
		groundedarea=groundedarea*scalefactor;
	}

	/*Clean up and return*/
	xDelete<IssmDouble>(xyz_list);
	return groundedarea;
}
/*}}}*/
bool       Tria::HasEdgeOnBase(){/*{{{*/

	IssmDouble values[NUMVERTICES];
	IssmDouble sum;

	/*Retrieve all inputs and parameters*/
	Element::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*/
	Element::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::IcefrontMassFlux(bool scaled){/*{{{*/

	/*Make sure there is an ice front here*/
	if(!IsIceInElement() || !IsIcefront()) return 0;

	/*Scaled not implemented yet...*/
	_assert_(!scaled);

	/*Get domain type*/
	int domaintype;
	parameters->FindParam(&domaintype,DomainTypeEnum);

	/*Get ice front coordinates*/
	IssmDouble *xyz_list = NULL;
	IssmDouble* xyz_front = NULL;
	this->GetVerticesCoordinates(&xyz_list);
	this->GetIcefrontCoordinates(&xyz_front,xyz_list,MaskIceLevelsetEnum);

	/*Get normal vector*/
	IssmDouble normal[3];
	this->NormalSection(&normal[0],xyz_front);
	//normal[0] = -normal[0];
	//normal[1] = -normal[1];

	/*Get inputs*/
	IssmDouble flux = 0.;
	IssmDouble vx,vy,thickness,Jdet;
	IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
	Input2* thickness_input=this->GetInput2(ThicknessEnum); _assert_(thickness_input);
	Input2* vx_input=NULL;
	Input2* vy_input=NULL;
	if(domaintype==Domain2DhorizontalEnum){
		vx_input=this->GetInput2(VxEnum); _assert_(vx_input);
		vy_input=this->GetInput2(VyEnum); _assert_(vy_input);
	}
	else{
		vx_input=this->GetInput2(VxAverageEnum); _assert_(vx_input);
		vy_input=this->GetInput2(VyAverageEnum); _assert_(vy_input);
	}

	/*Start looping on Gaussian points*/
	Gauss* gauss=this->NewGauss(xyz_list,xyz_front,3);
	for(int ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);
		thickness_input->GetInputValue(&thickness,gauss);
		vx_input->GetInputValue(&vx,gauss);
		vy_input->GetInputValue(&vy,gauss);
		this->JacobianDeterminantSurface(&Jdet,xyz_front,gauss);

		flux += rho_ice*Jdet*gauss->weight*thickness*(vx*normal[0] + vy*normal[1]);
	}
	delete gauss;
	return flux;
}
/*}}}*/
IssmDouble Tria::IcefrontMassFluxLevelset(bool scaled){/*{{{*/

	/*Make sure there is an ice front here*/
	if(!IsIceInElement() || !IsZeroLevelset(MaskIceLevelsetEnum)) return 0;

	/*Scaled not implemented yet...*/
	_assert_(!scaled);

	int               domaintype,index1,index2;
	const IssmPDouble epsilon = 1.e-15;
	IssmDouble        s1,s2;
	IssmDouble        gl[NUMVERTICES];
	IssmDouble        xyz_front[2][3];


	IssmDouble *xyz_list = NULL;
	this->GetVerticesCoordinates(&xyz_list);

	/*Recover parameters and values*/
	parameters->FindParam(&domaintype,DomainTypeEnum);
	Element::GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);

	/*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;

	if(domaintype==Domain2DverticalEnum){
		_error_("not implemented");
	}
	else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum || domaintype==Domain3DsurfaceEnum){
		int pt1 = 0;
		int pt2 = 1;
		if(gl[0]*gl[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=gl[2]/(gl[2]-gl[1]);
			s2=gl[2]/(gl[2]-gl[0]);
			if(gl[2]<0.){
				pt1 = 1; pt2 = 0;
			}
			xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
			xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
			xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
			xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
			xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
			xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
		}
		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

			/*Portion of the segments*/
			s1=gl[0]/(gl[0]-gl[1]);
			s2=gl[0]/(gl[0]-gl[2]);
			if(gl[0]<0.){
				pt1 = 1; pt2 = 0;
			}

			xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
			xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
			xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
			xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
			xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
			xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+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

			/*Portion of the segments*/
			s1=gl[1]/(gl[1]-gl[0]);
			s2=gl[1]/(gl[1]-gl[2]);
			if(gl[1]<0.){
				pt1 = 1; pt2 = 0;
			}

			xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
			xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
			xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
			xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
			xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
			xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
		}
		else{
			_error_("case not possible");
		}

	}
	else _error_("mesh type "<<EnumToStringx(domaintype)<<"not supported yet ");

	/*Some checks in debugging mode*/
	_assert_(s1>=0 && s1<=1.);
	_assert_(s2>=0 && s2<=1.);

	/*Get normal vector*/
	IssmDouble normal[3];
	this->NormalSection(&normal[0],&xyz_front[0][0]);
	normal[0] = -normal[0];
	normal[1] = -normal[1];

	/*Get inputs*/
	IssmDouble flux = 0.;
	IssmDouble vx,vy,thickness,Jdet;
	IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
	Input2* thickness_input=this->GetInput2(ThicknessEnum); _assert_(thickness_input);
	Input2* vx_input=NULL;
	Input2* vy_input=NULL;
	if(domaintype==Domain2DhorizontalEnum){
		vx_input=this->GetInput2(VxEnum); _assert_(vx_input);
		vy_input=this->GetInput2(VyEnum); _assert_(vy_input);
	}
	else{
		vx_input=this->GetInput2(VxAverageEnum); _assert_(vx_input);
		vy_input=this->GetInput2(VyAverageEnum); _assert_(vy_input);
	}

	/*Start looping on Gaussian points*/
	Gauss* gauss=this->NewGauss(xyz_list,&xyz_front[0][0],3);
	for(int ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);
		thickness_input->GetInputValue(&thickness,gauss);
		vx_input->GetInputValue(&vx,gauss);
		vy_input->GetInputValue(&vy,gauss);
		this->JacobianDeterminantSurface(&Jdet,&xyz_front[0][0],gauss);

		flux += rho_ice*Jdet*gauss->weight*thickness*(vx*normal[0] + vy*normal[1]);
	}
	delete gauss;
	return flux;
}
/*}}}*/
IssmDouble Tria::GroundinglineMassFlux(bool scaled){/*{{{*/

	/*Make sure there is a grounding line here*/
	if(!IsIceInElement()) return 0;
	if(!IsZeroLevelset(MaskGroundediceLevelsetEnum)) return 0;

	/*Scaled not implemented yet...*/
	_assert_(!scaled);

	int               domaintype,index1,index2;
	const IssmPDouble epsilon = 1.e-15;
	IssmDouble        s1,s2;
	IssmDouble        gl[NUMVERTICES];
	IssmDouble        xyz_front[2][3];

	IssmDouble *xyz_list = NULL;
	this->GetVerticesCoordinates(&xyz_list);

	/*Recover parameters and values*/
	parameters->FindParam(&domaintype,DomainTypeEnum);
	Element::GetInputListOnVertices(&gl[0],MaskGroundediceLevelsetEnum);

	/*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;

	if(domaintype==Domain2DverticalEnum){
		_error_("not implemented");
	}
	else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum || domaintype==Domain3DsurfaceEnum){
		int pt1 = 0;
		int pt2 = 1;
		if(gl[0]*gl[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=gl[2]/(gl[2]-gl[1]);
			s2=gl[2]/(gl[2]-gl[0]);
			if(gl[2]<0.){
				pt1 = 1; pt2 = 0;
			}
			xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
			xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
			xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
			xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
			xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
			xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
		}
		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

			/*Portion of the segments*/
			s1=gl[0]/(gl[0]-gl[1]);
			s2=gl[0]/(gl[0]-gl[2]);
			if(gl[0]<0.){
				pt1 = 1; pt2 = 0;
			}

			xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
			xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
			xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
			xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
			xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
			xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+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

			/*Portion of the segments*/
			s1=gl[1]/(gl[1]-gl[0]);
			s2=gl[1]/(gl[1]-gl[2]);
			if(gl[1]<0.){
				pt1 = 1; pt2 = 0;
			}

			xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
			xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
			xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
			xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
			xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
			xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
		}
		else{
			_error_("case not possible");
		}

	}
	else _error_("mesh type "<<EnumToStringx(domaintype)<<"not supported yet ");

	/*Some checks in debugging mode*/
	_assert_(s1>=0 && s1<=1.);
	_assert_(s2>=0 && s2<=1.);

	/*Get normal vector*/
	IssmDouble normal[3];
	this->NormalSection(&normal[0],&xyz_front[0][0]);

	/*Get inputs*/
	IssmDouble flux = 0.;
	IssmDouble vx,vy,thickness,Jdet;
	IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
	Input2* thickness_input=this->GetInput2(ThicknessEnum); _assert_(thickness_input);
	Input2* vx_input=NULL;
	Input2* vy_input=NULL;
	if(domaintype==Domain2DhorizontalEnum){
		vx_input=this->GetInput2(VxEnum); _assert_(vx_input);
		vy_input=this->GetInput2(VyEnum); _assert_(vy_input);
	}
	else{
		vx_input=this->GetInput2(VxAverageEnum); _assert_(vx_input);
		vy_input=this->GetInput2(VyAverageEnum); _assert_(vy_input);
	}

	/*Start looping on Gaussian points*/
	Gauss* gauss=this->NewGauss(xyz_list,&xyz_front[0][0],3);
	for(int ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);
		thickness_input->GetInputValue(&thickness,gauss);
		vx_input->GetInputValue(&vx,gauss);
		vy_input->GetInputValue(&vy,gauss);
		this->JacobianDeterminantSurface(&Jdet,&xyz_front[0][0],gauss);

		flux += rho_ice*Jdet*gauss->weight*thickness*(vx*normal[0] + vy*normal[1]);
	}

	return flux;
}
/*}}}*/
IssmDouble Tria::IceVolume(bool scaled){/*{{{*/

	/*The volume of a truncated prism is area_base * 1/numedges sum(length of edges)*/

	/*Intermediaries*/
	int i, numiceverts;
	IssmDouble area_base,surface,base,Haverage,scalefactor;
	IssmDouble Haux[NUMVERTICES], surfaces[NUMVERTICES], bases[NUMVERTICES];
	IssmDouble SFaux[NUMVERTICES], scalefactors[NUMVERTICES];
	IssmDouble s[2]; // s:fraction of intersected triangle edges, that lies inside ice
	int* indices=NULL;
	IssmDouble* H=NULL;
	IssmDouble* SF=NULL;

	if(!IsIceInElement())return 0.;
	//if(!IsIceOnlyInElement()) return 0;

	int domaintype;
	parameters->FindParam(&domaintype,DomainTypeEnum);

	if(false && IsIcefront()){
		//Assumption: linear ice thickness profile on element.
		//Hence ice thickness at intersection of levelset function with triangle edge is linear interpolation of ice thickness at vertices.
		this->GetLevelsetIntersection(&indices, &numiceverts, s, MaskIceLevelsetEnum, 0.);
		int numthk=numiceverts+2;
		H=xNew<IssmDouble>(numthk);
		//Correct area distortion caused by projection if requestion
		area_base=this->GetAreaIce();
		if(scaled==true){
			Element::GetInputListOnVertices(&scalefactors[0],MeshScaleFactorEnum);
			for(i=0;i<NUMVERTICES;i++) SFaux[i]= scalefactors[indices[i]]; //sort thicknesses in ice/noice
			switch(numiceverts){
				case 1: // average over triangle
					SF[0]=SFaux[0];
					SF[1]=SFaux[0]+s[0]*(SFaux[1]-SFaux[0]);
					SF[2]=SFaux[0]+s[1]*(SFaux[2]-SFaux[0]);
					break;
				case 2: // average over quadrangle
					SF[0]=SFaux[0];
					SF[1]=SFaux[1];
					SF[2]=SFaux[0]+s[0]*(SFaux[2]-SFaux[0]);
					SF[3]=SFaux[1]+s[1]*(SFaux[2]-SFaux[1]);
					break;
				default:
					_error_("Number of ice covered vertices wrong in Tria::IceVolume()");
					break;
			}
			scalefactor=0.;
			for(i=0;i<numthk;i++)	scalefactor+=SF[i];
			scalefactor/=IssmDouble(numthk);
			area_base=area_base*scalefactor;
		}
		Element::GetInputListOnVertices(&surfaces[0],SurfaceEnum);
		Element::GetInputListOnVertices(&bases[0],BaseEnum);
		for(i=0;i<NUMVERTICES;i++) Haux[i]= surfaces[indices[i]]-bases[indices[i]]; //sort thicknesses in ice/noice
		switch(numiceverts){
			case 1: // average over triangle
				H[0]=Haux[0];
				H[1]=Haux[0]+s[0]*(Haux[1]-Haux[0]);
				H[2]=Haux[0]+s[1]*(Haux[2]-Haux[0]);
				break;
			case 2: // average over quadrangle
				H[0]=Haux[0];
				H[1]=Haux[1];
				H[2]=Haux[0]+s[0]*(Haux[2]-Haux[0]);
				H[3]=Haux[1]+s[1]*(Haux[2]-Haux[1]);
				break;
			default:
				_error_("Number of ice covered vertices wrong in Tria::IceVolume()");
				break;
		}
		Haverage=0.;
		for(i=0;i<numthk;i++)	Haverage+=H[i];
		Haverage/=IssmDouble(numthk);
	}
	else{
		/*First get back the area of the base*/
		area_base=this->GetArea();
		if(scaled==true){
			Input2* scalefactor_input = this->GetInput2(MeshScaleFactorEnum); _assert_(scalefactor_input);
			scalefactor_input->GetInputAverage(&scalefactor);
			area_base=area_base*scalefactor;
		}

		/*Now get the average height*/
		Input2* surface_input = this->GetInput2(SurfaceEnum); _assert_(surface_input);
		Input2* base_input    = this->GetInput2(BaseEnum);    _assert_(base_input);
		surface_input->GetInputAverage(&surface);
		base_input->GetInputAverage(&base);
		Haverage=surface-base;
	}

	/*Cleanup & return: */
	xDelete<int>(indices);
	xDelete<IssmDouble>(H);
	xDelete<IssmDouble>(SF);

	if(domaintype==Domain2DverticalEnum){
	  return area_base;
	}
	else{
	  return area_base*Haverage;
	}
}
/*}}}*/
IssmDouble Tria::IceVolumeAboveFloatation(bool scaled){/*{{{*/

	/*The volume above floatation: H + rho_water/rho_ice * bathymetry */
	IssmDouble rho_ice,rho_water;
	IssmDouble base,surface,bed,bathymetry,scalefactor;
	IssmDouble xyz_list[NUMVERTICES][3];

	if(!IsIceInElement() || IsFloating())return 0;

	rho_ice=FindParam(MaterialsRhoIceEnum);
	rho_water=FindParam(MaterialsRhoSeawaterEnum);
	::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]));
	if(scaled==true){
		Input2* scalefactor_input = this->GetInput2(MeshScaleFactorEnum); _assert_(scalefactor_input);
		scalefactor_input->GetInputAverage(&scalefactor);
		base=base*scalefactor;
	}

	/*Now get the average height and bathymetry*/
	Input2* surface_input = this->GetInput2(SurfaceEnum); _assert_(surface_input);
	Input2* base_input    = this->GetInput2(BaseEnum);    _assert_(base_input);
	Input2* bed_input     = this->GetInput2(BedEnum);     _assert_(bed_input);
	if(!bed_input) _error_("Could not find bed");
	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::InputDepthAverageAtBase(int enum_type,int average_enum_type){/*{{{*/

	/*New input*/
	Input2* oldinput=NULL;
	Input2* newinput=NULL;

	/*copy input of enum_type*/
	oldinput=this->GetInput2(enum_type);
	if(!oldinput)_error_("could not find old input with enum: " << EnumToStringx(enum_type));
	newinput=oldinput->copy();

	/*Assign new name (average)*/
	newinput->ChangeEnum(average_enum_type);

	/*Add new input to current element*/
	_error_("not implemented");
}
/*}}}*/
void       Tria::InputUpdateFromIoModel(int index, IoModel* iomodel){ //i is the element index/*{{{*/

	/*Intermediaries*/
	int        i,j;
	int        tria_vertex_ids[3];
	IssmDouble nodeinputs[3];
	IssmDouble cmmininputs[3];
	IssmDouble cmmaxinputs[3];
	bool       control_analysis,ad_analysis   = false;
	int        num_control_type,num_responses;
	char**     controls = NULL;
	IssmDouble yts;

	/*Get parameters: */
	iomodel->FindConstant(&yts,"md.constants.yts");
	iomodel->FindConstant(&control_analysis,"md.inversion.iscontrol");
	iomodel->FindConstant(&ad_analysis, "md.autodiff.isautodiff");
	if(control_analysis && !ad_analysis) iomodel->FindConstant(&num_control_type,"md.inversion.num_control_parameters");
	if(control_analysis && !ad_analysis) iomodel->FindConstant(&num_responses,"md.inversion.num_cost_functions");
	if(control_analysis && ad_analysis) iomodel->FindConstant(&num_control_type,"md.autodiff.num_independent_objects");
	if(control_analysis && ad_analysis) iomodel->FindConstant(&num_responses,"md.autodiff.num_dependent_objects");

	/*Recover vertices ids needed to initialize inputs*/
	for(i=0;i<3;i++){
		tria_vertex_ids[i]=reCast<int>(iomodel->elements[3*index+i]); //ids for vertices are in the elements array from Matlab
	}
}
/*}}}*/
void       Tria::InputUpdateFromSolutionOneDof(IssmDouble* solution,int enum_type){/*{{{*/

	/*Intermediary*/
	int* doflist = NULL;

	/*Fetch number of nodes for this finite element*/
	int numnodes = this->NumberofNodes(this->element_type);

	/*Fetch dof list and allocate solution vector*/
	GetDofListLocal(&doflist,NoneApproximationEnum,GsetEnum);
	IssmDouble* values    = xNew<IssmDouble>(numnodes);

	/*Use the dof list to index into the solution vector: */
	for(int i=0;i<numnodes;i++){
		values[i]=solution[doflist[i]];
		if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in solution vector");
		if(xIsInf<IssmDouble>(values[i])) _error_("Inf found in solution vector, SID = " << this->Sid());
	}

	/*Add input to the element: */
	this->AddInput2(enum_type,values,this->element_type);

	/*Free ressources:*/
	xDelete<IssmDouble>(values);
	xDelete<int>(doflist);
}
/*}}}*/
void       Tria::InputUpdateFromVector(IssmDouble* vector, int name, int type){/*{{{*/

	/*Check that name is an element input*/
	if(!IsInputEnum(name)) _error_("Enum "<<EnumToStringx(name)<<" is not in IsInput");

	int         numnodes;
	IssmDouble  value;
	int         lidlist[NUMVERTICES];
	int        *doflist = NULL;
	IssmDouble *values  = NULL;

	GetVerticesLidList(&lidlist[0]);

	switch(type){
		case VertexLIdEnum:
			values = xNew<IssmDouble>(NUMVERTICES);
			for(int i=0;i<NUMVERTICES;i++){
				values[i]=vector[this->vertices[i]->Lid()];
				if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in vector");
				if(xIsInf<IssmDouble>(values[i])) _error_("Inf found in vector");
			}
			/*update input*/
			inputs2->SetTriaInput(name,P1Enum,NUMVERTICES,lidlist,values);
			break;

		case VertexPIdEnum:
			values = xNew<IssmDouble>(NUMVERTICES);
			for(int i=0;i<NUMVERTICES;i++){
				values[i]=vector[this->vertices[i]->Pid()];
				if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in vector");
				if(xIsInf<IssmDouble>(values[i])) _error_("Inf found in vector");
			}
			/*update input*/
			inputs2->SetTriaInput(name,P1Enum,NUMVERTICES,lidlist,values);
			break;

		case VertexSIdEnum:
			values = xNew<IssmDouble>(NUMVERTICES);
			for(int i=0;i<NUMVERTICES;i++){
				values[i]=vector[this->vertices[i]->Sid()];
				if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in vector");
				if(xIsInf<IssmDouble>(values[i])) _error_("Inf found in vector");
			}
			/*update input*/
			inputs2->SetTriaInput(name,P1Enum,NUMVERTICES,lidlist,values);
			break;

		case NodesEnum:
			/*Get number of nodes and dof list: */
			numnodes = this->NumberofNodes(this->element_type);
			values   = xNew<IssmDouble>(numnodes);
			GetDofList(&doflist,NoneApproximationEnum,GsetEnum);

			for(int i=0;i<numnodes;i++){
				values[i]=vector[doflist[i]];
				if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in vector");
				if(xIsInf<IssmDouble>(values[i])) _error_("Inf found in vector");
			}
			//this->inputs->AddInput(new TriaInput(name,values,this->element_type));
			_error_("not implemented");
			break;

		case NodeSIdEnum:
			/*Get number of nodes and dof list: */
			numnodes = this->NumberofNodes(this->element_type);
			values   = xNew<IssmDouble>(numnodes);

			for(int i=0;i<numnodes;i++){
				values[i]=vector[nodes[i]->Sid()];
				if(xIsNan<IssmDouble>(values[i])) _error_("NaN found in vector");
				if(xIsInf<IssmDouble>(values[i])) _error_("Inf found in vector");
			}
			if(this->element_type==P1Enum){
				inputs2->SetTriaInput(name,P1Enum,NUMVERTICES,lidlist,values);
			}
			else{
				inputs2->SetTriaInput(name,this->element_type,this->lid,numnodes,values);
			}
			break;

		case ElementEnum:
			value=vector[this->Sid()];
			if(xIsNan<IssmDouble>(value)) _error_("NaN found in vector");
			if(xIsInf<IssmDouble>(value)) _error_("Inf found in vector");
			/*update input*/
			//this->inputs->AddInput(new DoubleInput(name,value));
			//inputs2->SetTriaInput(name,P1Enum,NUMVERTICES,lidlist,values);
			_error_("not implemented");
			break;

		default:
			_error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
	}

	/*Clean-up*/
	xDelete<int>(doflist);
	xDelete<IssmDouble>(values);

}
/*}}}*/
bool       Tria::IsFaceOnBoundary(void){/*{{{*/

	IssmDouble values[NUMVERTICES];
	IssmDouble sum;

	/*Retrieve all inputs and parameters*/
	Element::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*/
	Element::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::IsZeroLevelset(int levelset_enum){/*{{{*/

	bool iszerols;
	IssmDouble ls[NUMVERTICES];

	/*Retrieve all inputs and parameters*/
	Element::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);
	this->GetJacobianDeterminant(pJdet,xyz_list,(GaussTria*)gauss);

}
/*}}}*/
void       Tria::JacobianDeterminantBase(IssmDouble* pJdet,IssmDouble* xyz_list_base,Gauss* gauss){/*{{{*/

	_assert_(gauss->Enum()==GaussTriaEnum);
	this->GetSegmentJacobianDeterminant(pJdet,xyz_list_base,(GaussTria*)gauss);

}
/*}}}*/
void       Tria::JacobianDeterminantSurface(IssmDouble* pJdet,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/

	_assert_(gauss->Enum()==GaussTriaEnum);
	this->GetSegmentJacobianDeterminant(pJdet,xyz_list,(GaussTria*)gauss);

}
/*}}}*/
void       Tria::JacobianDeterminantTop(IssmDouble* pJdet,IssmDouble* xyz_list_top,Gauss* gauss){/*{{{*/

	_assert_(gauss->Enum()==GaussTriaEnum);
	this->GetSegmentJacobianDeterminant(pJdet,xyz_list_top,(GaussTria*)gauss);

}
/*}}}*/
IssmDouble Tria::Masscon(IssmDouble* levelset){ /*{{{*/

	/*intermediary: */
	IssmDouble* values=NULL;
	Input2*     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->GetInput2(ThicknessEnum); _assert_(thickness_input);

	/*Retrieve material parameters: */
	rho_ice=FindParam(MaterialsRhoIceEnum);

	/*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()];
	}

	/*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;
}
/*}}}*/
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 vx1,vx2,vy1,vy2,h1,h2;

	/*Get material parameters :*/
	IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);

	/*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: */
	GaussTria* gauss_1=new GaussTria();
	gauss_1->GaussFromCoords(x1,y1,&xyz_list[0][0]);
	GaussTria* gauss_2=new GaussTria();
	gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]);

	/*Get segment length and normal (needs to be properly oriented)*/
	IssmDouble nx=cos(atan2(x1-x2,y2-y1));
	IssmDouble ny=sin(atan2(x1-x2,y2-y1));
	IssmDouble length=sqrt(pow(x2-x1,2)+pow(y2-y1,2));

	/*Get velocity and thickness*/
	this->parameters->FindParam(&domaintype,DomainTypeEnum);
	Input2* thickness_input=this->GetInput2(ThicknessEnum); _assert_(thickness_input);
	Input2* vx_input=NULL;
	Input2* vy_input=NULL;
	if(domaintype==Domain2DhorizontalEnum){
		vx_input=this->GetInput2(VxEnum); _assert_(vx_input);
		vy_input=this->GetInput2(VyEnum); _assert_(vy_input);
	}
	else{
		vx_input=this->GetInput2(VxAverageEnum); _assert_(vx_input);
		vy_input=this->GetInput2(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*(
				(1./3.*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*nx+
				(1./3.*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*ny
				);

	/*clean up and return:*/
	delete gauss_1;
	delete gauss_2;
	return mass_flux;
}
/*}}}*/
IssmDouble Tria::MassFlux(IssmDouble* segment){/*{{{*/
	return this->MassFlux(segment[0],segment[1],segment[2],segment[3],reCast<int>(segment[4]));
}
/*}}}*/
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);
	Input2* model_input=this->GetInput2(modelenum);   _assert_(model_input);
	Input2* observation_input=this->GetInput2(observationenum);_assert_(observation_input);
	Input2* weights_input     =this->GetInput2(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;
	}

	/* clean up and Return: */
	delete gauss;
	return Jelem;
}
/*}}}*/
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);
	Input2* weights_input     =this->GetInput2(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;
	}

	/* clean up and Return: */
	delete gauss;
	return Jelem;
}
/*}}}*/
Gauss*     Tria::NewGauss(void){/*{{{*/
	return new GaussTria();
}
/*}}}*/
Gauss*     Tria::NewGauss(int order){/*{{{*/
	return new GaussTria(order);
}
/*}}}*/
Gauss*     Tria::NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order){/*{{{*/

	IssmDouble  area_coordinates[2][3];
	GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,2);
	return new GaussTria(area_coordinates,order);
}
/*}}}*/
Gauss*     Tria::NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,bool mainlyfloating,int order){/*{{{*/

	return new GaussTria(point1,fraction1,fraction2,mainlyfloating,order);
}
/*}}}*/
Gauss*     Tria::NewGauss(int point1,IssmDouble fraction1,IssmDouble fraction2,int order){/*{{{*/

	return new GaussTria(point1,fraction1,fraction2,order);
}
/*}}}*/
Gauss*     Tria::NewGauss(IssmDouble* xyz_list, IssmDouble* xyz_list_front,int order_horiz,int order_vert){/*{{{*/

	IssmDouble  area_coordinates[2][3];
	GetAreaCoordinates(&area_coordinates[0][0],xyz_list_front,xyz_list,2);
	return new GaussTria(area_coordinates,order_vert);
}
/*}}}*/
Gauss*     Tria::NewGaussBase(int order){/*{{{*/

	int indices[2];
	this->EdgeOnBaseIndices(&indices[0],&indices[1]);
	return new GaussTria(indices[0],indices[1],order);
}
/*}}}*/
Gauss*     Tria::NewGaussTop(int order){/*{{{*/

	int indices[2];
	this->EdgeOnSurfaceIndices(&indices[0],&indices[1]);
	return new GaussTria(indices[0],indices[1],order);
}
/*}}}*/
void       Tria::NodalFunctions(IssmDouble* basis, Gauss* gauss){/*{{{*/

	_assert_(gauss->Enum()==GaussTriaEnum);
	this->GetNodalFunctions(basis,(GaussTria*)gauss,this->element_type);

}
/*}}}*/
void       Tria::NodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/

	_assert_(gauss->Enum()==GaussTriaEnum);
	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss,this->element_type);

}
/*}}}*/
void       Tria::NodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/

	_assert_(gauss->Enum()==GaussTriaEnum);
	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss,this->VelocityInterpolation());

}
/*}}}*/
void       Tria::NodalFunctionsPressure(IssmDouble* basis, Gauss* gauss){/*{{{*/

	_assert_(gauss->Enum()==GaussTriaEnum);
	this->GetNodalFunctions(basis,(GaussTria*)gauss,this->PressureInterpolation());

}
/*}}}*/
void       Tria::NodalFunctionsP1(IssmDouble* basis, Gauss* gauss){/*{{{*/

	_assert_(gauss->Enum()==GaussTriaEnum);
	this->GetNodalFunctions(basis,(GaussTria*)gauss,P1Enum);

}
/*}}}*/
void       Tria::NodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/

	_assert_(gauss->Enum()==GaussTriaEnum);
	this->GetNodalFunctionsDerivatives(dbasis,xyz_list,(GaussTria*)gauss,P1Enum);

}
/*}}}*/
void       Tria::NodalFunctionsP2(IssmDouble* basis, Gauss* gauss){/*{{{*/

	_assert_(gauss->Enum()==GaussTriaEnum);
	this->GetNodalFunctions(basis,(GaussTria*)gauss,P2Enum);

}
/*}}}*/
void       Tria::NodalFunctionsTensor(IssmDouble* basis, Gauss* gauss){/*{{{*/

	_assert_(gauss->Enum()==GaussTriaEnum);
	this->GetNodalFunctions(basis,(GaussTria*)gauss,this->TensorInterpolation());

}
/*}}}*/
void       Tria::NodalFunctionsVelocity(IssmDouble* basis, Gauss* gauss){/*{{{*/

	_assert_(gauss->Enum()==GaussTriaEnum);
	this->GetNodalFunctions(basis,(GaussTria*)gauss,this->VelocityInterpolation());

}
/*}}}*/
int        Tria::NodalValue(IssmDouble* pvalue, int index, int natureofdataenum){/*{{{*/

	int         found = 0;
	IssmDouble  value;
	GaussTria  *gauss = NULL;

	/*First, serarch the input: */
	Input2* data=this->GetInput2(natureofdataenum);

	/*figure out if we have the vertex id: */
	found=0;
	for(int i=0;i<NUMVERTICES;i++){
		if(index==vertices[i]->Sid()){
			/*Do we have natureofdataenum in our inputs? :*/
			if(data){
				/*ok, we are good. retrieve value of input at vertex :*/
				gauss=new GaussTria(); gauss->GaussVertex(i);
				data->GetInputValue(&value,gauss);
				found=1;
				break;
			}
		}
	}

	/*clean-up*/
	delete gauss;

	if(found)*pvalue=value;
	return found;
}
/*}}}*/
void       Tria::NormalBase(IssmDouble* bed_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]);

	bed_normal[0]= + vector[1]/norm;
	bed_normal[1]= - vector[0]/norm;
	_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*/
	int index1,index2;
	IssmDouble vector[2];
	IssmDouble norm;

	this->EdgeOnSurfaceIndices(&index1,&index2);
	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]);

	top_normal[0]= + vector[1]/norm;
	top_normal[1]= - vector[0]/norm;
	_assert_(top_normal[1]>0);
}
/*}}}*/
int        Tria::ObjectEnum(void){/*{{{*/

	return TriaEnum;

}
/*}}}*/
int        Tria::NumberofNodesPressure(void){/*{{{*/
	return TriaRef::NumberofNodes(this->PressureInterpolation());
}
/*}}}*/
int        Tria::NumberofNodesVelocity(void){/*{{{*/
	return TriaRef::NumberofNodes(this->VelocityInterpolation());
}
/*}}}*/
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=FindParam(MaterialsRhoSeawaterEnum);
	rho_ice=FindParam(MaterialsRhoIceEnum);
	density=rho_ice/rho_water;
	Element::GetInputListOnVertices(&h[0],ThicknessEnum);
	Element::GetInputListOnVertices(&r[0],BedEnum);
	Element::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::PressureInterpolation(void){/*{{{*/
	return TriaRef::PressureInterpolation(this->element_type);
}
/*}}}*/
void       Tria::ReduceMatrices(ElementMatrix* Ke,ElementVector* pe){/*{{{*/

	/*Static condensation if requested*/
	if(pe){
		if(this->element_type==MINIcondensedEnum){
			int indices[2]={6,7};
			pe->StaticCondensation(Ke,2,&indices[0]);
		}
		else if(this->element_type==P1bubblecondensedEnum){
			int size   = nodes[3]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
			int offset = 0;
			for(int i=0;i<3;i++) offset+=nodes[i]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
			int* indices=xNew<int>(size);
			for(int i=0;i<size;i++) indices[i] = offset+i;
			pe->StaticCondensation(Ke,size,indices);
			xDelete<int>(indices);
		}
	}

	if(Ke){
		if(this->element_type==MINIcondensedEnum){
			int indices[2]={6,7};
			Ke->StaticCondensation(2,&indices[0]);
		}
		else if(this->element_type==P1bubblecondensedEnum){
			int size   = nodes[3]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
			int offset = 0;
			for(int i=0;i<3;i++) offset+=nodes[i]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
			int* indices=xNew<int>(size);
			for(int i=0;i<size;i++) indices[i] = offset+i;
			Ke->StaticCondensation(size,indices);
			xDelete<int>(indices);
		}
	}

}
/*}}}*/
void       Tria::ResetFSBasalBoundaryCondition(void){/*{{{*/

	int numnodes = this->NumberofNodesVelocity();

	int          approximation;
	IssmDouble*  vertexonbase= NULL;
	IssmDouble   slope,groundedice;
	IssmDouble   xz_plane[6];

	/*For FS only: we want the CS to be tangential to the bedrock*/
	this->GetInput2Value(&approximation,ApproximationEnum);
	if(!HasNodeOnBase() ||  approximation!=FSApproximationEnum) return;

	/*Get inputs*/
	Input2* slope_input=this->GetInput2(BedSlopeXEnum);                             _assert_(slope_input);
	Input2* groundedicelevelset_input=this->GetInput2(MaskGroundediceLevelsetEnum); _assert_(groundedicelevelset_input);
	vertexonbase = xNew<IssmDouble>(numnodes);
	this->GetInputListOnNodesVelocity(&vertexonbase[0],MeshVertexonbaseEnum);

	/*Loop over basal nodes and update their CS*/
	GaussTria* gauss = new GaussTria();
	for(int i=0;i<this->NumberofNodesVelocity();i++){

		if(vertexonbase[i]==1){
			gauss->GaussNode(this->VelocityInterpolation(),i);
			slope_input->GetInputValue(&slope,gauss);
			groundedicelevelset_input->GetInputValue(&groundedice,gauss);
			IssmDouble theta = atan(slope);

			/*New X axis                  New Z axis*/
			xz_plane[0]=cos(theta);       xz_plane[3]=0.;
			xz_plane[1]=sin(theta);       xz_plane[4]=0.;
			xz_plane[2]=0.;               xz_plane[5]=1.;

			if(groundedice>=0){
				this->nodes[i]->DofInSSet(1); //vy
			}
			else{
				this->nodes[i]->DofInFSet(1); //vy
			}

			XZvectorsToCoordinateSystem(&this->nodes[i]->coord_system[0][0],&xz_plane[0]);
		}
	}

	/*cleanup*/
	xDelete<IssmDouble>(vertexonbase);
	delete gauss;
}
/*}}}*/
void       Tria::ResetHooks(){/*{{{*/

	if(this->nodes) xDelete<Node*>(this->nodes);
	this->nodes=NULL;
	this->vertices=NULL;
	this->material=NULL;
	this->parameters=NULL;

	//deal with ElementHook mother class
	for(int i=0;i<this->numanalyses;i++) if(this->hnodes[i]) this->hnodes[i]->reset();
	this->hvertices->reset();
	this->hmaterial->reset();
	if(this->hneighbors) this->hneighbors->reset();

}
/*}}}*/
void       Tria::RignotMeltParameterization(){/*{{{*/

   IssmDouble A, B, alpha, beta;
	IssmDouble bed,qsg,qsg_basin,TF,yts;
	int numbasins;
	IssmDouble basinid[NUMVERTICES];
	IssmDouble* basin_icefront_area=NULL;

	/* Coefficients */
	A    = 3e-4;
	B    = 0.15;
	alpha = 0.39;
	beta = 1.18;

	/*Get inputs*/
	Input2* bed_input = this->GetInput2(BedEnum);                     _assert_(bed_input);
	Input2* qsg_input = this->GetInput2(FrontalForcingsSubglacialDischargeEnum);		 _assert_(qsg_input);
	Input2* TF_input  = this->GetInput2(FrontalForcingsThermalForcingEnum);          _assert_(TF_input);
	Element::GetInputListOnVertices(&basinid[0],FrontalForcingsBasinIdEnum);

	this->FindParam(&yts, ConstantsYtsEnum);
	this->parameters->FindParam(&numbasins,FrontalForcingsNumberofBasinsEnum);
	this->parameters->FindParam(&basin_icefront_area,&numbasins,FrontalForcingsBasinIcefrontAreaEnum);

	IssmDouble meltrates[NUMVERTICES];  //frontal melt-rate

	/* Start looping on the number of vertices: */
	GaussTria* gauss=new GaussTria();
	for(int iv=0;iv<NUMVERTICES;iv++){
		gauss->GaussVertex(iv);

		/* Get variables */
		bed_input->GetInputValue(&bed,gauss);
		qsg_input->GetInputValue(&qsg,gauss);
		TF_input->GetInputValue(&TF,gauss);

		if(basin_icefront_area[reCast<int>(basinid[iv])-1]==0.) meltrates[iv]=0.;
		else{
			/* change the unit of qsg (m^3/d -> m/d) with ice front area */
			qsg_basin=qsg/basin_icefront_area[reCast<int>(basinid[iv])-1];

			/* calculate melt rates */
			meltrates[iv]=((A*max(-bed,0.)*pow(max(qsg_basin,0.),alpha)+B)*pow(max(TF,0.),beta))/86400; //[m/s]
		}

		if(xIsNan<IssmDouble>(meltrates[iv])) _error_("NaN found in vector");
		if(xIsInf<IssmDouble>(meltrates[iv])) _error_("Inf found in vector");
	}

	/*Add input*/
	this->AddInput2(CalvingMeltingrateEnum,&meltrates[0],P1Enum);

	/*Cleanup and return*/
	xDelete<IssmDouble>(basin_icefront_area);
	delete gauss;
}
/*}}}*/
void       Tria::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index,int offset,int N, int M){/*{{{*/

	IssmDouble  values[NUMVERTICES];
	int         lidlist[NUMVERTICES];
	int         idlist[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(!IsInputEnum(control_enum)) return;

	this->GetVerticesLidList(&lidlist[0]);
	ElementInput2* input=this->inputs2->GetControlInput2Data(control_enum,"value");   _assert_(input);

	/*Get values on vertices*/
	for(int n=0;n<N;n++){
		for(int i=0;i<NUMVERTICES;i++){
			idlist[i]=offset + this->vertices[i]->Sid()+n*M;
			values[i]=vector[idlist[i]];
		}
		if(input->ObjectEnum()==TriaInput2Enum){
			input->SetInput(P1Enum,NUMVERTICES,&lidlist[0],&values[0]);
		}
		else if(input->ObjectEnum()==TransientInput2Enum){
			_error_("not implemented");
			//Input* new_input = new TriaInput(control_enum,values,P1Enum);
			//controlinput->SetInput(new_input,n);
			//controlinput->Configure(parameters);
		}
		else _error_("Type not supported");
	}
}
/*}}}*/
void       Tria::SetControlInputsFromVector(IssmDouble* vector,int control_enum,int control_index){/*{{{*/

	IssmDouble  values[NUMVERTICES];
	int         idlist[NUMVERTICES];
	int         lidlist[NUMVERTICES];

	/*Get Domain type*/
	int domaintype;
	parameters->FindParam(&domaintype,DomainTypeEnum);

	/*Specific case for depth averaged quantities*/
	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(!IsInputEnum(control_enum)) return;

	/*prepare index list*/
	this->GetVerticesLidList(&lidlist[0]);
	GradientIndexing(&idlist[0],control_index);

	/*Get values on vertices*/
	for(int i=0;i<NUMVERTICES;i++){
		values[i]=vector[idlist[i]];
	}

	/*Set Input*/
	ElementInput2* input=this->inputs2->GetControlInput2Data(control_enum,"value");   _assert_(input);
	input->SetInput(P1Enum,NUMVERTICES,&lidlist[0],&values[0]);
}
/*}}}*/
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();
	}

}
/*}}}*/
void       Tria::SetElementInput(int enum_in,IssmDouble value){/*{{{*/

	this->SetElementInput(this->inputs2,enum_in,value);

}
/*}}}*/
void       Tria::SetElementInput(Inputs2* inputs2,int enum_in,IssmDouble value){/*{{{*/

	_assert_(inputs2);
	inputs2->SetTriaInput(enum_in,P0Enum,this->lid,value);

}
/*}}}*/
void       Tria::SetElementInput(Inputs2* inputs2,int numindices,int* indices,IssmDouble* values,int enum_in){/*{{{*/

	_assert_(inputs2);
	inputs2->SetTriaInput(enum_in,P1Enum,numindices,indices,values);

}
/*}}}*/
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;

	/*go into parameters and get the analysis_counter: */
	this->parameters->FindParam(&analysis_counter,AnalysisCounterEnum);

	/*Create Seg*/
	Seg* seg=new Seg();
	seg->id=this->id;
	seg->sid=this->sid;
	seg->lid=this->lid;
	seg->inputs2=this->inputs2;
	seg->parameters=this->parameters;
	seg->element_type=P1Enum; //Only P1 CG for now (TO BE CHANGED)
	this->SpawnSegHook(xDynamicCast<ElementHook*>(seg),index1,index2);

	seg->iscollapsed = 1;
	seg->collapsed_ids[0] = index1;
	seg->collapsed_ids[1] = index2;

	/*Spawn material*/
	seg->material=(Material*)this->material->copy2(seg);

	/*recover nodes, material*/
	seg->nodes    = (Node**)seg->hnodes[analysis_counter]->deliverp();
	seg->vertices = (Vertex**)seg->hvertices->deliverp();

	/*Return new Seg*/
	return seg;
}
/*}}}*/
Element*   Tria::SpawnTopElement(void){/*{{{*/

	int index1,index2;
	int domaintype;

	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");
	}
}
/*}}}*/
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*/
	Input2* vx_input=this->GetInput2(VxEnum);                                  _assert_(vx_input);
	Input2* vy_input=this->GetInput2(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-14);
	}

	/*Add input*/
	this->AddInput2(StrainRateparallelEnum,&strainparallel[0],P1DGEnum);

	/*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*/
	Input2* vx_input=this->GetInput2(VxEnum);                                  _assert_(vx_input);
	Input2* vy_input=this->GetInput2(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-14);
	}

	/*Add input*/
	this->AddInput2(StrainRateperpendicularEnum,&strainperpendicular[0],P1DGEnum);

	/*Clean up and return*/
	delete gauss;
	xDelete<IssmDouble>(xyz_list);
}
/*}}}*/
IssmDouble Tria::SurfaceArea(void){/*{{{*/

	IssmDouble S;
	IssmDouble normal[3];
	IssmDouble v13[3],v23[3];
	IssmDouble xyz_list[NUMVERTICES][3];

	/*If on water, return 0: */
	if(!IsIceInElement()) return 0.;

	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);

	for(int i=0;i<3;i++){
		v13[i]=xyz_list[0][i]-xyz_list[2][i];
		v23[i]=xyz_list[1][i]-xyz_list[2][i];
	}

	normal[0]=v13[1]*v23[2]-v13[2]*v23[1];
	normal[1]=v13[2]*v23[0]-v13[0]*v23[2];
	normal[2]=v13[0]*v23[1]-v13[1]*v23[0];

	S = 0.5 * sqrt(normal[0]*normal[0] + normal[1]*normal[1] + normal[2]*normal[2]);

	/*Return: */
	return S;
}
/*}}}*/
int        Tria::TensorInterpolation(void){/*{{{*/
	return TriaRef::TensorInterpolation(this->element_type);
}
/*}}}*/
IssmDouble Tria::TimeAdapt(void){/*{{{*/

	/*intermediary: */
	IssmDouble C;
	IssmDouble xyz_list[NUMVERTICES][3];

	/*get CFL coefficient:*/
	this->parameters->FindParam(&C,TimesteppingCflCoefficientEnum);

	/*Get for Vx and Vy, the max of abs value: */
	Input2* vx_input = this->GetInput2(VxEnum); _assert_(vx_input);
	Input2* vy_input = this->GetInput2(VyEnum); _assert_(vy_input);
	IssmDouble maxabsvx = vx_input->GetInputMaxAbs();
	IssmDouble maxabsvy = vy_input->GetInputMaxAbs();

	/* Get node coordinates and dof list: */
	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);

	IssmDouble minx=xyz_list[0][0];
	IssmDouble maxx=xyz_list[0][0];
	IssmDouble miny=xyz_list[0][1];
	IssmDouble maxy=xyz_list[0][1];

	for(int i=1;i<NUMVERTICES;i++){
		if(xyz_list[i][0]<minx) minx=xyz_list[i][0];
		if(xyz_list[i][0]>maxx) maxx=xyz_list[i][0];
		if(xyz_list[i][1]<miny) miny=xyz_list[i][1];
		if(xyz_list[i][1]>maxy) maxy=xyz_list[i][1];
	}
	IssmDouble dx=maxx-minx;
	IssmDouble dy=maxy-miny;

	/*CFL criterion: */
	IssmDouble dt = C/(maxabsvx/dx+maxabsvy/dy);

	return dt;
}
/*}}}*/
IssmDouble Tria::TotalCalvingFluxLevelset(bool scaled){/*{{{*/

	/*Make sure there is an ice front here*/
	if(!IsIceInElement() || !IsZeroLevelset(MaskIceLevelsetEnum)) return 0;

	/*Scaled not implemented yet...*/
	_assert_(!scaled);

	int               domaintype,index1,index2;
	const IssmPDouble epsilon = 1.e-15;
	IssmDouble        s1,s2;
	IssmDouble        gl[NUMVERTICES];
	IssmDouble        xyz_front[2][3];


	IssmDouble *xyz_list = NULL;
	this->GetVerticesCoordinates(&xyz_list);

	/*Recover parameters and values*/
	parameters->FindParam(&domaintype,DomainTypeEnum);
	Element::GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);

	/*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;

	if(domaintype==Domain2DverticalEnum){
		_error_("not implemented");
	}
	else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum || domaintype==Domain3DsurfaceEnum){
		int pt1 = 0;
		int pt2 = 1;
		if(gl[0]*gl[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=gl[2]/(gl[2]-gl[1]);
			s2=gl[2]/(gl[2]-gl[0]);
			if(gl[2]<0.){
				pt1 = 1; pt2 = 0;
			}
			xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
			xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
			xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
			xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
			xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
			xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
		}
		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

			/*Portion of the segments*/
			s1=gl[0]/(gl[0]-gl[1]);
			s2=gl[0]/(gl[0]-gl[2]);
			if(gl[0]<0.){
				pt1 = 1; pt2 = 0;
			}

			xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
			xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
			xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
			xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
			xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
			xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+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

			/*Portion of the segments*/
			s1=gl[1]/(gl[1]-gl[0]);
			s2=gl[1]/(gl[1]-gl[2]);
			if(gl[1]<0.){
				pt1 = 1; pt2 = 0;
			}

			xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
			xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
			xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
			xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
			xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
			xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
		}
		else{
			_error_("case not possible");
		}

	}
	else _error_("mesh type "<<EnumToStringx(domaintype)<<"not supported yet ");

	/*Some checks in debugging mode*/
	_assert_(s1>=0 && s1<=1.);
	_assert_(s2>=0 && s2<=1.);

	/*Get normal vector*/
	IssmDouble normal[3];
	this->NormalSection(&normal[0],&xyz_front[0][0]);
	normal[0] = -normal[0];
	normal[1] = -normal[1];

	/*Get inputs*/
	IssmDouble flux = 0.;
	IssmDouble calvingratex,calvingratey,thickness,Jdet;
	IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
	Input2* thickness_input=this->GetInput2(ThicknessEnum); _assert_(thickness_input);
	Input2* calvingratex_input=NULL;
	Input2* calvingratey_input=NULL;
	if(domaintype==Domain2DhorizontalEnum){
		calvingratex_input=this->GetInput2(CalvingratexEnum); _assert_(calvingratex_input);
		calvingratey_input=this->GetInput2(CalvingrateyEnum); _assert_(calvingratey_input);
	}
	else{
		calvingratex_input=this->GetInput2(CalvingratexAverageEnum); _assert_(calvingratex_input);
		calvingratey_input=this->GetInput2(CalvingrateyAverageEnum); _assert_(calvingratey_input);
	}

	/*Start looping on Gaussian points*/
	Gauss* gauss=this->NewGauss(xyz_list,&xyz_front[0][0],3);
	for(int ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);
		thickness_input->GetInputValue(&thickness,gauss);
		calvingratex_input->GetInputValue(&calvingratex,gauss);
		calvingratey_input->GetInputValue(&calvingratey,gauss);
		this->JacobianDeterminantSurface(&Jdet,&xyz_front[0][0],gauss);

		flux += rho_ice*Jdet*gauss->weight*thickness*(calvingratex*normal[0] + calvingratey*normal[1]);
	}

	return flux;
}
/*}}}*/
IssmDouble Tria::TotalCalvingMeltingFluxLevelset(bool scaled){/*{{{*/

	/*Make sure there is an ice front here*/
	if(!IsIceInElement() || !IsZeroLevelset(MaskIceLevelsetEnum)) return 0;

	/*Scaled not implemented yet...*/
	_assert_(!scaled);

	int               domaintype,index1,index2;
	const IssmPDouble epsilon = 1.e-15;
	IssmDouble        s1,s2;
	IssmDouble        gl[NUMVERTICES];
	IssmDouble        xyz_front[2][3];


	IssmDouble *xyz_list = NULL;
	this->GetVerticesCoordinates(&xyz_list);

	/*Recover parameters and values*/
	parameters->FindParam(&domaintype,DomainTypeEnum);
	Element::GetInputListOnVertices(&gl[0],MaskIceLevelsetEnum);

	/*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;

	if(domaintype==Domain2DverticalEnum){
		_error_("not implemented");
	}
	else if(domaintype==Domain2DhorizontalEnum || domaintype==Domain3DEnum || domaintype==Domain3DsurfaceEnum){
		int pt1 = 0;
		int pt2 = 1;
		if(gl[0]*gl[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=gl[2]/(gl[2]-gl[1]);
			s2=gl[2]/(gl[2]-gl[0]);
			if(gl[2]<0.){
				pt1 = 1; pt2 = 0;
			}
			xyz_front[pt2][0]=xyz_list[3*2+0]+s1*(xyz_list[3*1+0]-xyz_list[3*2+0]);
			xyz_front[pt2][1]=xyz_list[3*2+1]+s1*(xyz_list[3*1+1]-xyz_list[3*2+1]);
			xyz_front[pt2][2]=xyz_list[3*2+2]+s1*(xyz_list[3*1+2]-xyz_list[3*2+2]);
			xyz_front[pt1][0]=xyz_list[3*2+0]+s2*(xyz_list[3*0+0]-xyz_list[3*2+0]);
			xyz_front[pt1][1]=xyz_list[3*2+1]+s2*(xyz_list[3*0+1]-xyz_list[3*2+1]);
			xyz_front[pt1][2]=xyz_list[3*2+2]+s2*(xyz_list[3*0+2]-xyz_list[3*2+2]);
		}
		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

			/*Portion of the segments*/
			s1=gl[0]/(gl[0]-gl[1]);
			s2=gl[0]/(gl[0]-gl[2]);
			if(gl[0]<0.){
				pt1 = 1; pt2 = 0;
			}

			xyz_front[pt1][0]=xyz_list[3*0+0]+s1*(xyz_list[3*1+0]-xyz_list[3*0+0]);
			xyz_front[pt1][1]=xyz_list[3*0+1]+s1*(xyz_list[3*1+1]-xyz_list[3*0+1]);
			xyz_front[pt1][2]=xyz_list[3*0+2]+s1*(xyz_list[3*1+2]-xyz_list[3*0+2]);
			xyz_front[pt2][0]=xyz_list[3*0+0]+s2*(xyz_list[3*2+0]-xyz_list[3*0+0]);
			xyz_front[pt2][1]=xyz_list[3*0+1]+s2*(xyz_list[3*2+1]-xyz_list[3*0+1]);
			xyz_front[pt2][2]=xyz_list[3*0+2]+s2*(xyz_list[3*2+2]-xyz_list[3*0+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

			/*Portion of the segments*/
			s1=gl[1]/(gl[1]-gl[0]);
			s2=gl[1]/(gl[1]-gl[2]);
			if(gl[1]<0.){
				pt1 = 1; pt2 = 0;
			}

			xyz_front[pt2][0]=xyz_list[3*1+0]+s1*(xyz_list[3*0+0]-xyz_list[3*1+0]);
			xyz_front[pt2][1]=xyz_list[3*1+1]+s1*(xyz_list[3*0+1]-xyz_list[3*1+1]);
			xyz_front[pt2][2]=xyz_list[3*1+2]+s1*(xyz_list[3*0+2]-xyz_list[3*1+2]);
			xyz_front[pt1][0]=xyz_list[3*1+0]+s2*(xyz_list[3*2+0]-xyz_list[3*1+0]);
			xyz_front[pt1][1]=xyz_list[3*1+1]+s2*(xyz_list[3*2+1]-xyz_list[3*1+1]);
			xyz_front[pt1][2]=xyz_list[3*1+2]+s2*(xyz_list[3*2+2]-xyz_list[3*1+2]);
		}
		else{
			_error_("case not possible");
		}

	}
	else _error_("mesh type "<<EnumToStringx(domaintype)<<"not supported yet ");

	/*Some checks in debugging mode*/
	_assert_(s1>=0 && s1<=1.);
	_assert_(s2>=0 && s2<=1.);

	/*Get normal vector*/
	IssmDouble normal[3];
	this->NormalSection(&normal[0],&xyz_front[0][0]);
	normal[0] = -normal[0];
	normal[1] = -normal[1];

	/*Get inputs*/
	IssmDouble flux = 0.;
	IssmDouble calvingratex,calvingratey,vx,vy,vel,meltingrate,meltingratex,meltingratey,thickness,Jdet;
	IssmDouble rho_ice=FindParam(MaterialsRhoIceEnum);
	Input2* thickness_input=this->GetInput2(ThicknessEnum); _assert_(thickness_input);
	Input2* calvingratex_input=NULL;
	Input2* calvingratey_input=NULL;
	Input2* vx_input=NULL;
	Input2* vy_input=NULL;
	Input2* meltingrate_input=NULL;
	if(domaintype==Domain2DhorizontalEnum){
		calvingratex_input=this->GetInput2(CalvingratexEnum); _assert_(calvingratex_input);
		calvingratey_input=this->GetInput2(CalvingrateyEnum); _assert_(calvingratey_input);
		vx_input=this->GetInput2(VxEnum); _assert_(vx_input);
		vy_input=this->GetInput2(VyEnum); _assert_(vy_input);
		meltingrate_input=this->GetInput2(CalvingMeltingrateEnum); _assert_(meltingrate_input);
	}
	else{
		calvingratex_input=this->GetInput2(CalvingratexAverageEnum); _assert_(calvingratex_input);
		calvingratey_input=this->GetInput2(CalvingrateyAverageEnum); _assert_(calvingratey_input);
	}

	/*Start looping on Gaussian points*/
	Gauss* gauss=this->NewGauss(xyz_list,&xyz_front[0][0],3);
	for(int ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);
		thickness_input->GetInputValue(&thickness,gauss);
		calvingratex_input->GetInputValue(&calvingratex,gauss);
		calvingratey_input->GetInputValue(&calvingratey,gauss);
		vx_input->GetInputValue(&vx,gauss);
		vy_input->GetInputValue(&vy,gauss);
		vel=vx*vx+vy*vy;
		meltingrate_input->GetInputValue(&meltingrate,gauss);
		meltingratex=meltingrate*vx/(sqrt(vel)+1.e-14);
		meltingratey=meltingrate*vy/(sqrt(vel)+1.e-14);
		this->JacobianDeterminantSurface(&Jdet,&xyz_front[0][0],gauss);

		flux += rho_ice*Jdet*gauss->weight*thickness*((calvingratex+meltingratex)*normal[0] + (calvingratey+meltingratey)*normal[1]);
	}

	return flux;
}
/*}}}*/
IssmDouble Tria::TotalFloatingBmb(bool scaled){/*{{{*/

	/*The fbmb[kg yr-1] of one element is area[m2] * melting_rate [kg m^-2 yr^-1]*/
	int        point1;
	bool       mainlyfloating;
	IssmDouble fbmb=0;
	IssmDouble rho_ice,fraction1,fraction2,floatingmelt,Jdet,scalefactor;
	IssmDouble Total_Fbmb=0;
	IssmDouble xyz_list[NUMVERTICES][3];
	Gauss*     gauss     = NULL;

   if(!IsIceInElement())return 0;

	/*Get material parameters :*/
	rho_ice=FindParam(MaterialsRhoIceEnum);
	Input2* floatingmelt_input = this->GetInput2(BasalforcingsFloatingiceMeltingRateEnum); _assert_(floatingmelt_input);
	Input2* gllevelset_input   = this->GetInput2(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
	Input2* scalefactor_input  = NULL;
	if(scaled==true){
		scalefactor_input = this->GetInput2(MeshScaleFactorEnum); _assert_(scalefactor_input);
	}
	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);

	this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
	/* Start  looping on the number of gaussian points: */
	gauss = this->NewGauss(point1,fraction1,fraction2,1-mainlyfloating,3);
	for(int ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);
		this->JacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
		floatingmelt_input->GetInputValue(&floatingmelt,gauss);
		if(scaled==true){
			scalefactor_input->GetInputValue(&scalefactor,gauss);
		}
		else scalefactor=1;
		fbmb+=floatingmelt*Jdet*gauss->weight*scalefactor;
	}

   Total_Fbmb=rho_ice*fbmb;	        // from volume to mass

	/*Return: */
	delete gauss;
	return Total_Fbmb;
}
/*}}}*/
IssmDouble Tria::TotalGroundedBmb(bool scaled){/*{{{*/

	/*The gbmb[kg yr-1] of one element is area[m2] * gounded melting rate [kg m^-2 yr^-1]*/
	int        point1;
	bool       mainlyfloating;
	IssmDouble gbmb=0;
	IssmDouble rho_ice,fraction1,fraction2,groundedmelt,Jdet,scalefactor;
	IssmDouble Total_Gbmb=0;
	IssmDouble xyz_list[NUMVERTICES][3];
	Gauss*     gauss     = NULL;

   if(!IsIceInElement())return 0;

	/*Get material parameters :*/
	rho_ice=FindParam(MaterialsRhoIceEnum);
	Input2* groundedmelt_input = this->GetInput2(BasalforcingsGroundediceMeltingRateEnum); _assert_(groundedmelt_input);
	Input2* gllevelset_input = this->GetInput2(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
	Input2* scalefactor_input = NULL;
	if(scaled==true){
		scalefactor_input = this->GetInput2(MeshScaleFactorEnum); _assert_(scalefactor_input);
	}
	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);

	this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
	/* Start  looping on the number of gaussian points: */
	gauss = this->NewGauss(point1,fraction1,fraction2,mainlyfloating,2);
	for(int ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);
		this->JacobianDeterminant(&Jdet,&xyz_list[0][0],gauss);
		groundedmelt_input->GetInputValue(&groundedmelt,gauss);
		if(scaled==true){
			scalefactor_input->GetInputValue(&scalefactor,gauss);
		}
		else scalefactor=1;
		gbmb+=groundedmelt*Jdet*gauss->weight*scalefactor;
	}

   Total_Gbmb=rho_ice*gbmb;	        // from volume to mass

	/*Return: */
	delete gauss;
	return Total_Gbmb;
}
/*}}}*/
IssmDouble Tria::TotalSmb(bool scaled){/*{{{*/

	/*The smb[kg yr-1] of one element is area[m2] * smb [kg m^-2 yr^-1]*/
	IssmDouble base,smb,rho_ice,scalefactor;
	IssmDouble Total_Smb=0;
	IssmDouble xyz_list[NUMVERTICES][3];

	/*Get material parameters :*/
	rho_ice=FindParam(MaterialsRhoIceEnum);

   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*/
	Input2* smb_input = this->GetInput2(SmbMassBalanceEnum); _assert_(smb_input);
	smb_input->GetInputAverage(&smb);	// average smb on element in m ice s-1
	if(scaled==true){
		Input2* scalefactor_input = this->GetInput2(MeshScaleFactorEnum); _assert_(scalefactor_input);
		scalefactor_input->GetInputAverage(&scalefactor);// average scalefactor on element
	}
	else{
		scalefactor=1.;
	}
   Total_Smb=rho_ice*base*smb*scalefactor;	// smb on element in kg s-1

	/*Return: */
	return Total_Smb;
}
/*}}}*/
void       Tria::Update(Inputs2* inputs2,int index, IoModel* iomodel,int analysis_counter,int analysis_type,int finiteelement_type){/*{{{*/

	/*Intermediaries*/
	int  numnodes;
	int* tria_node_ids = NULL;

	/*Checks if debuging*/
	_assert_(iomodel->elements);
	_assert_(index==this->sid);

	/*Recover element type*/
	this->element_type_list[analysis_counter]=finiteelement_type;

	/*Recover nodes ids needed to initialize the node hook.*/
	switch(finiteelement_type){
		case P0DGEnum:
			numnodes        = 1;
			tria_node_ids   = xNew<int>(numnodes);
			tria_node_ids[0]= index + 1;
			break;
		case P1Enum:
			numnodes        = 3;
			tria_node_ids   = xNew<int>(numnodes);
			tria_node_ids[0]=iomodel->elements[3*index+0];
			tria_node_ids[1]=iomodel->elements[3*index+1];
			tria_node_ids[2]=iomodel->elements[3*index+2];
			break;
		case P1DGEnum:
			numnodes        = 3;
			tria_node_ids   = xNew<int>(numnodes);
			tria_node_ids[0]=3*index+1;
			tria_node_ids[1]=3*index+2;
			tria_node_ids[2]=3*index+3;
			break;
		case P1bubbleEnum: case P1bubblecondensedEnum:
			numnodes        = 4;
			tria_node_ids   = xNew<int>(numnodes);
			tria_node_ids[0]=iomodel->elements[3*index+0];
			tria_node_ids[1]=iomodel->elements[3*index+1];
			tria_node_ids[2]=iomodel->elements[3*index+2];
			tria_node_ids[3]=iomodel->numberofvertices+index+1;
			break;
		case P2Enum:
			numnodes        = 6;
			tria_node_ids   = xNew<int>(numnodes);
			tria_node_ids[0]=iomodel->elements[3*index+0];
			tria_node_ids[1]=iomodel->elements[3*index+1];
			tria_node_ids[2]=iomodel->elements[3*index+2];
			tria_node_ids[3]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+0]+1;
			tria_node_ids[4]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+1]+1;
			tria_node_ids[5]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+2]+1;
			break;
		case P2bubbleEnum: case P2bubblecondensedEnum:
			numnodes        = 7;
			tria_node_ids   = xNew<int>(numnodes);
			tria_node_ids[0]=iomodel->elements[3*index+0];
			tria_node_ids[1]=iomodel->elements[3*index+1];
			tria_node_ids[2]=iomodel->elements[3*index+2];
			tria_node_ids[3]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+0]+1;
			tria_node_ids[4]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+1]+1;
			tria_node_ids[5]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+2]+1;
			tria_node_ids[6]=iomodel->numberofvertices+iomodel->numberofedges+index+1;
			break;
		case P1P1Enum: case P1P1GLSEnum:
			numnodes        = 6;
			tria_node_ids   = xNew<int>(numnodes);
			tria_node_ids[0]=iomodel->elements[3*index+0];
			tria_node_ids[1]=iomodel->elements[3*index+1];
			tria_node_ids[2]=iomodel->elements[3*index+2];

			tria_node_ids[3]=iomodel->numberofvertices+iomodel->elements[3*index+0];
			tria_node_ids[4]=iomodel->numberofvertices+iomodel->elements[3*index+1];
			tria_node_ids[5]=iomodel->numberofvertices+iomodel->elements[3*index+2];
			break;
		case MINIEnum: case MINIcondensedEnum:
			numnodes       = 7;
			tria_node_ids  = xNew<int>(numnodes);
			tria_node_ids[0]=iomodel->elements[3*index+0];
			tria_node_ids[1]=iomodel->elements[3*index+1];
			tria_node_ids[2]=iomodel->elements[3*index+2];
			tria_node_ids[3]=iomodel->numberofvertices+index+1;

			tria_node_ids[4]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[3*index+0];
			tria_node_ids[5]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[3*index+1];
			tria_node_ids[6]=iomodel->numberofvertices+iomodel->numberofelements+iomodel->elements[3*index+2];
			break;
		case TaylorHoodEnum:
		case XTaylorHoodEnum:
			numnodes        = 9;
			tria_node_ids   = xNew<int>(numnodes);
			tria_node_ids[0]=iomodel->elements[3*index+0];
			tria_node_ids[1]=iomodel->elements[3*index+1];
			tria_node_ids[2]=iomodel->elements[3*index+2];
			tria_node_ids[3]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+0]+1;
			tria_node_ids[4]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+1]+1;
			tria_node_ids[5]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+2]+1;

			tria_node_ids[6]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[3*index+0];
			tria_node_ids[7]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[3*index+1];
			tria_node_ids[8]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->elements[3*index+2];
			break;
		case LATaylorHoodEnum:
			numnodes        = 6;
			tria_node_ids   = xNew<int>(numnodes);
			tria_node_ids[0]=iomodel->elements[3*index+0];
			tria_node_ids[1]=iomodel->elements[3*index+1];
			tria_node_ids[2]=iomodel->elements[3*index+2];
			tria_node_ids[3]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+0]+1;
			tria_node_ids[4]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+1]+1;
			tria_node_ids[5]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+2]+1;
			break;
		case CrouzeixRaviartEnum:
			numnodes        = 10;
			tria_node_ids   = xNew<int>(numnodes);
			tria_node_ids[0]=iomodel->elements[3*index+0];
			tria_node_ids[1]=iomodel->elements[3*index+1];
			tria_node_ids[2]=iomodel->elements[3*index+2];
			tria_node_ids[3]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+0]+1;
			tria_node_ids[4]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+1]+1;
			tria_node_ids[5]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+2]+1;
			tria_node_ids[6]=iomodel->numberofvertices+iomodel->numberofedges+index+1;

			tria_node_ids[7]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofelements+3*index+1;
			tria_node_ids[8]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofelements+3*index+2;
			tria_node_ids[9]=iomodel->numberofvertices+iomodel->numberofedges+iomodel->numberofelements+3*index+3;
			break;
		case LACrouzeixRaviartEnum:
			numnodes        = 7;
			tria_node_ids   = xNew<int>(numnodes);
			tria_node_ids[0]=iomodel->elements[3*index+0];
			tria_node_ids[1]=iomodel->elements[3*index+1];
			tria_node_ids[2]=iomodel->elements[3*index+2];
			tria_node_ids[3]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+0]+1;
			tria_node_ids[4]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+1]+1;
			tria_node_ids[5]=iomodel->numberofvertices+iomodel->elementtoedgeconnectivity[3*index+2]+1;
			tria_node_ids[6]=iomodel->numberofvertices+iomodel->numberofedges+index+1;
			break;
		default:
			_error_("Finite element "<<EnumToStringx(finiteelement_type)<<" not supported yet");
	}

	/*hooks: */
	this->SetHookNodes(tria_node_ids,numnodes,analysis_counter); this->nodes=NULL;
	xDelete<int>(tria_node_ids);
}
/*}}}*/
void       Tria::UpdateConstraintsExtrudeFromBase(void){/*{{{*/

	if(!HasNodeOnBase()) return;

	int        extrusioninput;
	IssmDouble value,isonbase;

	this->parameters->FindParam(&extrusioninput,InputToExtrudeEnum);
	Input2* input = this->GetInput2(extrusioninput);      _assert_(input);
	Input2* onbase = this->GetInput2(MeshVertexonbaseEnum); _assert_(onbase);

	GaussTria* gauss=new GaussTria();
	for(int iv=0;iv<this->NumberofNodes(this->element_type);iv++){
		gauss->GaussNode(this->element_type,iv);
		onbase->GetInputValue(&isonbase,gauss);
		if(isonbase==1.){
			input->GetInputValue(&value,gauss);
			this->nodes[iv]->ApplyConstraint(0,value);
		}
	}
	delete gauss;

}
/*}}}*/
void       Tria::UpdateConstraintsExtrudeFromTop(void){/*{{{*/

	if(!HasNodeOnSurface()) return;

	int        extrusioninput;
	IssmDouble value,isonsurface;

	this->parameters->FindParam(&extrusioninput,InputToExtrudeEnum);
	Input2* input = this->GetInput2(extrusioninput); _assert_(input);
	Input2* onsurf = this->GetInput2(MeshVertexonsurfaceEnum); _assert_(onsurf);

	GaussTria* gauss=new GaussTria();
	for(int iv=0;iv<this->NumberofNodes(this->element_type);iv++){
		gauss->GaussNode(this->element_type,iv);
		onsurf->GetInputValue(&isonsurface,gauss);
		if(isonsurface==1.){
			input->GetInputValue(&value,gauss);
			this->nodes[iv]->ApplyConstraint(0,value);
		}
	}
	delete gauss;
}
/*}}}*/
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;
}
/*}}}*/
void       Tria::ValueP1DerivativesOnGauss(IssmDouble* dvalue,IssmDouble* values,IssmDouble* xyz_list,Gauss* gauss){/*{{{*/
	TriaRef::GetInputDerivativeValue(dvalue,values,xyz_list,gauss,P1Enum);
}
/*}}}*/
void       Tria::ValueP1OnGauss(IssmDouble* pvalue,IssmDouble* values,Gauss* gauss){/*{{{*/
	TriaRef::GetInputValue(pvalue,values,gauss,P1Enum);
}
/*}}}*/
int        Tria::VelocityInterpolation(void){/*{{{*/
	return TriaRef::VelocityInterpolation(this->element_type);
}
/*}}}*/
int        Tria::VertexConnectivity(int vertexindex){/*{{{*/
	_assert_(this->vertices);
	return this->vertices[vertexindex]->Connectivity();
}
/*}}}*/
void       Tria::WriteFieldIsovalueSegment(DataSet* segments,int fieldenum,IssmDouble fieldvalue){/*{{{*/

	_assert_(fieldvalue==0.); //field value != 0 not implemented yet

	/*Get field on vertices (we do not allow for higher order elements!!)*/
	IssmDouble lsf[NUMVERTICES];
	Element::GetInputListOnVertices(&lsf[0],fieldenum);

	/*1. check that we do cross fieldvalue in this element*/
	IssmDouble minvalue = lsf[0];
	IssmDouble maxvalue = lsf[0];
	for(int i=1;i<NUMVERTICES;i++){
		if(lsf[i]>maxvalue) maxvalue = lsf[i];
		if(lsf[i]<minvalue) minvalue = lsf[i];
	}
	if(minvalue>fieldvalue) return;
	if(maxvalue<fieldvalue) return;

	/*2. Find coordinates of where levelset crosses 0*/
	int         numiceverts;
	IssmDouble  s[2],x[2],y[2];
	int        *indices = NULL;
	this->GetLevelsetIntersection(&indices, &numiceverts,&s[0],fieldenum,fieldvalue);
	_assert_(numiceverts);

	/*3 Write coordinates*/
	IssmDouble  xyz_list[NUMVERTICES][3];
	::GetVerticesCoordinates(&xyz_list[0][0],this->vertices,NUMVERTICES);
	int counter = 0;
	if((numiceverts>0) && (numiceverts<NUMVERTICES)){
		for(int i=0;i<numiceverts;i++){
			for(int n=numiceverts;n<NUMVERTICES;n++){ // iterate over no-ice vertices
				x[counter] = xyz_list[indices[i]][0]+s[counter]*(xyz_list[indices[n]][0]-xyz_list[indices[i]][0]);
				y[counter] = xyz_list[indices[i]][1]+s[counter]*(xyz_list[indices[n]][1]-xyz_list[indices[i]][1]);
				counter++;
			}
		}
	}
	else if(numiceverts==NUMVERTICES){ //NUMVERTICES ice vertices: calving front lies on element edge

		for(int i=0;i<NUMVERTICES;i++){
			if(lsf[indices[i]]==0.){
				x[counter]=xyz_list[indices[i]][0];
				y[counter]=xyz_list[indices[i]][1];
				counter++;
			}
			if(counter==2) break;
		}
		if(counter==1){
			/*We actually have only 1 vertex on levelset, write a single point as a segment*/
			x[counter]=x[0];
			y[counter]=y[0];
			counter++;
		}
	}
	else{
		_error_("not sure what's going on here...");
	}

	/*4. Write segment*/
	_assert_(counter==2);
	segments->AddObject(new Contour<IssmDouble>(segments->Size()+1,2,&x[0],&y[0],false));

	/*Cleanup and return*/
	xDelete<int>(indices);
}
/*}}}*/

#ifdef _HAVE_GIAIVINS_
void       Tria::GiaDeflection(Vector<IssmDouble>* wg,Vector<IssmDouble>* dwgdt,IssmDouble* x, IssmDouble* y){/*{{{*/

	IssmDouble xyz_list[NUMVERTICES][3];

	/*gia solution parameters:*/
	IssmDouble lithosphere_thickness,mantle_viscosity;

	/*output: */
	IssmDouble  wi;
	IssmDouble  dwidt;

	/*arguments to GiaDeflectionCorex: */
	GiaDeflectionCoreArgs arguments;

	/*how many dofs are we working with here? */
	int gsize;
	IssmDouble yts;
	this->parameters->FindParam(&gsize,MeshNumberofverticesEnum);
	this->parameters->FindParam(&yts,ConstantsYtsEnum);

	/*recover gia solution parameters: */
	int cross_section_shape;
	this->parameters->FindParam(&cross_section_shape,GiaCrossSectionShapeEnum);

	/*what time is it? :*/
	IssmDouble currenttime;
	this->parameters->FindParam(&currenttime,TimeEnum);

	/*recover material parameters: */
	IssmDouble lithosphere_shear_modulus = FindParam(MaterialsLithosphereShearModulusEnum);
	IssmDouble lithosphere_density       = FindParam(MaterialsLithosphereDensityEnum);
	IssmDouble mantle_shear_modulus      = FindParam(MaterialsMantleShearModulusEnum);
	IssmDouble mantle_density            = FindParam(MaterialsMantleDensityEnum);
	IssmDouble rho_ice                   = FindParam(MaterialsRhoIceEnum);

	/*pull thickness averages! */
	IssmDouble *hes      = NULL;
	IssmDouble *times    = NULL;
	int         numtimes;
	this->GetInputAveragesUpToCurrentTime(ThicknessEnum,&hes,&times,&numtimes,currenttime);

	/*recover mantle viscosity: */
	Input2* mantle_viscosity_input=this->GetInput2(GiaMantleViscosityEnum);
	if (!mantle_viscosity_input)_error_("mantle viscosity input needed to compute gia deflection!");
	mantle_viscosity_input->GetInputAverage(&mantle_viscosity);

	/*recover lithosphere thickness: */
	Input2* lithosphere_thickness_input=this->GetInput2(GiaLithosphereThicknessEnum);
	if (!lithosphere_thickness_input)_error_("lithosphere thickness input needed to compute gia deflection!");
	lithosphere_thickness_input->GetInputAverage(&lithosphere_thickness);

	/*pull area of this Tria: */
	IssmDouble area=this->GetArea();

	/*element radius: */
	IssmDouble re=sqrt(area/PI);

	/*figure out gravity center of our element: */
	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
	IssmDouble x0=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0;
	IssmDouble y0=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0;

	/*start loading GiaDeflectionCore arguments: */
	arguments.re=re;
	arguments.hes=hes;
	arguments.times=times;
	arguments.numtimes=numtimes;
	arguments.currenttime=currenttime;
	arguments.lithosphere_shear_modulus=lithosphere_shear_modulus;
	arguments.lithosphere_density=lithosphere_density;
	arguments.mantle_shear_modulus=mantle_shear_modulus;
	arguments.mantle_viscosity=mantle_viscosity;
	arguments.mantle_density=mantle_density;
	arguments.lithosphere_thickness=lithosphere_thickness;
	arguments.rho_ice=rho_ice;
	arguments.idisk=this->id;
	arguments.iedge=cross_section_shape;
	arguments.yts=yts;

	for(int i=0;i<gsize;i++){
		/*compute distance from the center of the tria to the vertex i: */
		IssmDouble xi=x[i];
		IssmDouble yi=y[i];
		IssmDouble ri=sqrt(pow(xi-x0,2)+pow(yi-y0,2));

		/*load ri onto arguments for this vertex i: */
		arguments.ri=ri;

		/*for this Tria, compute contribution to rebound at vertex i: */
		GiaDeflectionCorex(&wi,&dwidt,&arguments);

		/*plug value into solution vector: */
		wg->SetValue(i,wi,ADD_VAL);
		dwgdt->SetValue(i,dwidt,ADD_VAL);
	}

	/*Free ressources: */
	xDelete<IssmDouble>(hes);
	xDelete<IssmDouble>(times);

	return;
}
/*}}}*/
#endif
#ifdef _HAVE_ESA_
void    Tria::EsaGeodetic2D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,Vector<IssmDouble>* pX,Vector<IssmDouble>* pY,IssmDouble* xx,IssmDouble* yy){ /*{{{*/

	/*diverse:*/
	int gsize;
	IssmDouble xyz_list[NUMVERTICES][3];
	IssmDouble area;
	IssmDouble earth_radius = 6371012.0;	// Earth's radius [m]
	IssmDouble I;		//ice/water loading
	IssmDouble rho_ice, rho_earth;

	/*precomputed elastic green functions:*/
	IssmDouble* U_elastic_precomputed = NULL;
	IssmDouble* H_elastic_precomputed = NULL;
	int         M, hemi;

	/*computation of Green functions:*/
	IssmDouble* U_elastic= NULL;
	IssmDouble* N_elastic= NULL;
	IssmDouble* E_elastic= NULL;
	IssmDouble* X_elastic= NULL;
	IssmDouble* Y_elastic= NULL;

	/*optimization:*/
	bool store_green_functions=false;

	/*Compute ice thickness change: */
	Input2* deltathickness_input=this->GetInput2(EsaDeltathicknessEnum);
	if (!deltathickness_input)_error_("delta thickness input needed to compute elastic adjustment!");
	deltathickness_input->GetInputAverage(&I);

	/*early return if we are not on the (ice) loading point: */
	if(I==0) return;

	/*recover material parameters: */
	rho_ice=FindParam(MaterialsRhoIceEnum);
	rho_earth=FindParam(MaterialsEarthDensityEnum);

	/*how many dofs are we working with here? */
	this->parameters->FindParam(&gsize,MeshNumberofverticesEnum);

	/*which hemisphere? for north-south, east-west components*/
	this->parameters->FindParam(&hemi,EsaHemisphereEnum);

	/*compute area of element:*/
	area=GetArea();

	/*figure out gravity center of our element (Cartesian): */
	IssmDouble x_element, y_element;
	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
	x_element=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0;
	y_element=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0;

	/*recover elastic Green's functions for displacement:*/
	DoubleVecParam* U_parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(EsaUElasticEnum)); _assert_(U_parameter);
	DoubleVecParam* H_parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(EsaHElasticEnum)); _assert_(H_parameter);
	U_parameter->GetParameterValueByPointer(&U_elastic_precomputed,&M);
	H_parameter->GetParameterValueByPointer(&H_elastic_precomputed,&M);

	/*initialize: */
	U_elastic=xNewZeroInit<IssmDouble>(gsize);
	N_elastic=xNewZeroInit<IssmDouble>(gsize);
	E_elastic=xNewZeroInit<IssmDouble>(gsize);
	X_elastic=xNewZeroInit<IssmDouble>(gsize);
	Y_elastic=xNewZeroInit<IssmDouble>(gsize);

	int* indices=xNew<int>(gsize);
	IssmDouble* U_values=xNewZeroInit<IssmDouble>(gsize);
	IssmDouble* N_values=xNewZeroInit<IssmDouble>(gsize);
	IssmDouble* E_values=xNewZeroInit<IssmDouble>(gsize);
	IssmDouble* X_values=xNewZeroInit<IssmDouble>(gsize);
	IssmDouble* Y_values=xNewZeroInit<IssmDouble>(gsize);
	IssmDouble dx, dy, dist, alpha, ang, ang2;
	IssmDouble N_azim, E_azim, X_azim, Y_azim;

	for(int i=0;i<gsize;i++){

		indices[i]=i;

		IssmDouble N_azim=0;
		IssmDouble E_azim=0;

		/*Compute alpha angle between centroid and current vertex: */
		dx = x_element - xx[i];		dy = y_element - yy[i];
		dist = sqrt(pow(dx,2)+pow(dy,2));						// distance between vertex and elemental centroid [m]
		alpha = dist*360.0/(2*PI*earth_radius) * PI/180.0;	// [in radians] 360 degree = 2*pi*earth_radius

		/*Compute azimuths, both north and east components: */
		ang = PI/2 - atan2(dy,dx);		// this is bearing angle!
		Y_azim = cos(ang);
		X_azim = sin(ang);

		/*Elastic component  (from Eq 17 in Adhikari et al, GMD 2015): */
		int index=reCast<int,IssmDouble>(alpha/PI*(M-1));
		U_elastic[i] += U_elastic_precomputed[index];
		Y_elastic[i] += H_elastic_precomputed[index]*Y_azim;
		X_elastic[i] += H_elastic_precomputed[index]*X_azim;

		/*Add all components to the pUp solution vectors:*/
		U_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*U_elastic[i];
		Y_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*Y_elastic[i];
		X_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*X_elastic[i];

		/*North-south, East-west components */
		if (hemi == -1) {
			ang2 = PI/2 - atan2(yy[i],xx[i]);
		}
		else if (hemi == 1) {
			ang2 = PI/2 - atan2(-yy[i],-xx[i]);
		}
		if (hemi != 0){
			N_azim = Y_azim*cos(ang2) + X_azim*sin(ang2);
			E_azim = X_azim*cos(ang2) - Y_azim*sin(ang2);
			N_elastic[i] += H_elastic_precomputed[index]*N_azim;
			E_elastic[i] += H_elastic_precomputed[index]*E_azim;
			N_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*N_elastic[i];
			E_values[i]+=3*rho_ice/rho_earth*area/(4*PI*pow(earth_radius,2))*I*E_elastic[i];
		}
	}

	pUp->SetValues(gsize,indices,U_values,ADD_VAL);
	pNorth->SetValues(gsize,indices,N_values,ADD_VAL);
	pEast->SetValues(gsize,indices,E_values,ADD_VAL);
	pX->SetValues(gsize,indices,X_values,ADD_VAL);
	pY->SetValues(gsize,indices,Y_values,ADD_VAL);

	/*free ressources:*/
	xDelete<int>(indices);
	xDelete<IssmDouble>(U_values); xDelete<IssmDouble>(N_values); xDelete<IssmDouble>(E_values);
	xDelete<IssmDouble>(U_elastic); xDelete<IssmDouble>(N_elastic); xDelete<IssmDouble>(E_elastic);
	xDelete<IssmDouble>(X_values); xDelete<IssmDouble>(Y_values);
	xDelete<IssmDouble>(X_elastic); xDelete<IssmDouble>(Y_elastic);

	return;
}
/*}}}*/
void    Tria::EsaGeodetic3D(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea){ /*{{{*/

	/*diverse:*/
	int gsize;
	bool spherical=true;
	IssmDouble llr_list[NUMVERTICES][3];
	IssmDouble xyz_list[NUMVERTICES][3];
	IssmDouble area;
	IssmDouble I;		//ice/water loading
	IssmDouble late,longe,re;
	IssmDouble lati,longi,ri;
	IssmDouble rho_ice,rho_earth;
	IssmDouble minlong=400;
	IssmDouble maxlong=-20;

	/*precomputed elastic green functions:*/
	IssmDouble* U_elastic_precomputed = NULL;
	IssmDouble* H_elastic_precomputed = NULL;
	int         M;

	/*computation of Green functions:*/
	IssmDouble* U_elastic= NULL;
	IssmDouble* N_elastic= NULL;
	IssmDouble* E_elastic= NULL;

	/*optimization:*/
	bool store_green_functions=false;

	/*Compute ice thickness change: */
	Input2* deltathickness_input=this->GetInput2(EsaDeltathicknessEnum);
	if (!deltathickness_input)_error_("delta thickness input needed to compute elastic adjustment!");
	deltathickness_input->GetInputAverage(&I);

	/*early return if we are not on the (ice) loading point: */
	if(I==0) return;

	/*recover material parameters: */
	rho_ice=FindParam(MaterialsRhoIceEnum);
	rho_earth=FindParam(MaterialsEarthDensityEnum);

	/*how many dofs are we working with here? */
	this->parameters->FindParam(&gsize,MeshNumberofverticesEnum);

	/*compute area of element:*/
	area=GetAreaSpherical();

	/*element centroid (spherical): */
	/* Where is the centroid of this element?:{{{*/
	::GetVerticesCoordinates(&llr_list[0][0],this->vertices,NUMVERTICES,spherical);

	minlong=400; maxlong=-20;
	for (int i=0;i<NUMVERTICES;i++){
		llr_list[i][0]=(90-llr_list[i][0]);
		if(llr_list[i][1]<0)llr_list[i][1]=180+(180+llr_list[i][1]);
		if(llr_list[i][1]>maxlong)maxlong=llr_list[i][1];
		if(llr_list[i][1]<minlong)minlong=llr_list[i][1];
	}
	if(minlong==0 && maxlong>180){
		if (llr_list[0][1]==0)llr_list[0][1]=360;
		if (llr_list[1][1]==0)llr_list[1][1]=360;
		if (llr_list[2][1]==0)llr_list[2][1]=360;
	}

	// correction at the north pole
	if(llr_list[0][0]==0)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
	if(llr_list[1][0]==0)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
	if(llr_list[2][0]==0)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;

	//correction at the south pole
	if(llr_list[0][0]==180)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
	if(llr_list[1][0]==180)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
	if(llr_list[2][0]==180)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;

	late=(llr_list[0][0]+llr_list[1][0]+llr_list[2][0])/3.0;
	longe=(llr_list[0][1]+llr_list[1][1]+llr_list[2][1])/3.0;

	late=90-late;
	if(longe>180)longe=(longe-180)-180;

	late=late/180*PI;
	longe=longe/180*PI;
	/*}}}*/

	/*figure out gravity center of our element (Cartesian): */
	IssmDouble x_element, y_element, z_element;
	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
	x_element=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0;
	y_element=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0;
	z_element=(xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2])/3.0;

	/*recover elastic Green's functions for displacement:*/
	DoubleVecParam* U_parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(EsaUElasticEnum)); _assert_(U_parameter);
	DoubleVecParam* H_parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(EsaHElasticEnum)); _assert_(H_parameter);
	U_parameter->GetParameterValueByPointer(&U_elastic_precomputed,&M);
	H_parameter->GetParameterValueByPointer(&H_elastic_precomputed,&M);

	/*initialize: */
	U_elastic=xNewZeroInit<IssmDouble>(gsize);
	N_elastic=xNewZeroInit<IssmDouble>(gsize);
	E_elastic=xNewZeroInit<IssmDouble>(gsize);

	int* indices=xNew<int>(gsize);
	IssmDouble* U_values=xNewZeroInit<IssmDouble>(gsize);
	IssmDouble* N_values=xNewZeroInit<IssmDouble>(gsize);
	IssmDouble* E_values=xNewZeroInit<IssmDouble>(gsize);
	IssmDouble alpha;
	IssmDouble delPhi,delLambda;
	IssmDouble dx, dy, dz, x, y, z;
	IssmDouble N_azim, E_azim;

	for(int i=0;i<gsize;i++){

		indices[i]=i;

		/*Compute alpha angle between centroid and current vertex: */
		lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;

		delPhi=fabs(lati-late); delLambda=fabs(longi-longe);
		alpha=2.*asin(sqrt(pow(sin(delPhi/2),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));

		/*Compute azimuths, both north and east components: */
		x = xx[i]; y = yy[i]; z = zz[i];
		if(latitude[i]==90){
			x=1e-12; y=1e-12;
		}
		if(latitude[i]==-90){
			x=1e-12; y=1e-12;
		}
		dx = x_element-x; dy = y_element-y; dz = z_element-z;
		N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
		E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);

		/*Elastic component  (from Eq 17 in Adhikari et al, GMD 2015): */
		int index=reCast<int,IssmDouble>(alpha/PI*(M-1));
		U_elastic[i] += U_elastic_precomputed[index];
		N_elastic[i] += H_elastic_precomputed[index]*N_azim;
		E_elastic[i] += H_elastic_precomputed[index]*E_azim;

		/*Add all components to the pUp solution vectors:*/
		U_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*U_elastic[i];
		N_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*N_elastic[i];
		E_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*E_elastic[i];
	}
	pUp->SetValues(gsize,indices,U_values,ADD_VAL);
	pNorth->SetValues(gsize,indices,N_values,ADD_VAL);
	pEast->SetValues(gsize,indices,E_values,ADD_VAL);

	/*free ressources:*/
	xDelete<int>(indices);
	xDelete<IssmDouble>(U_values); xDelete<IssmDouble>(N_values); xDelete<IssmDouble>(E_values);
	xDelete<IssmDouble>(U_elastic); xDelete<IssmDouble>(N_elastic); xDelete<IssmDouble>(E_elastic);

	return;
}
/*}}}*/
#endif
#ifdef _HAVE_SEALEVELRISE_
IssmDouble    Tria::OceanArea(void){ /*{{{*/

	if(IsWaterInElement()) return GetAreaSpherical();
	else return 0;

}
/*}}}*/
IssmDouble Tria::OceanAverage(IssmDouble* Sg){ /*{{{*/

	if(IsWaterInElement()){

		IssmDouble area;

		/*Compute area of element:*/
		area=GetAreaSpherical();

		/*Average Sg over vertices:*/
		IssmDouble Sg_avg=0; for(int i=0;i<NUMVERTICES;i++) Sg_avg+=Sg[this->vertices[i]->Sid()]/NUMVERTICES;

		/*return: */
		return area*Sg_avg;
	}
	else return 0;

}
/*}}}*/
void	Tria::SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,IssmDouble eartharea){/*{{{*/
	/*early return if we are not on an ice cap OR ocean:*/
	if(!IsIceOnlyInElement() && !IsWaterInElement()){
		dI_list[0] = 0.0; // this is important!!!
		dI_list[1] = 0.0; // this is important!!!
		dI_list[2] = 0.0; // this is important!!!
		return;
	}

	/*Compute area of element:*/
	IssmDouble area;
	area=GetAreaSpherical();

	/*Compute lat,long,radius of elemental centroid: */
	bool spherical=true;
	IssmDouble llr_list[NUMVERTICES][3];
	IssmDouble late,longe,re;
	/* Where is the centroid of this element?:{{{*/
	::GetVerticesCoordinates(&llr_list[0][0],this->vertices,NUMVERTICES,spherical);

	IssmDouble minlong=400;
	IssmDouble maxlong=-20;
	for (int i=0;i<NUMVERTICES;i++){
		llr_list[i][0]=(90-llr_list[i][0]);
		if(llr_list[i][1]<0)llr_list[i][1]=180+(180+llr_list[i][1]);
		if(llr_list[i][1]>maxlong)maxlong=llr_list[i][1];
		if(llr_list[i][1]<minlong)minlong=llr_list[i][1];
	}
	if(minlong==0 && maxlong>180){
		if (llr_list[0][1]==0)llr_list[0][1]=360;
		if (llr_list[1][1]==0)llr_list[1][1]=360;
		if (llr_list[2][1]==0)llr_list[2][1]=360;
	}

	// correction at the north pole
	if(llr_list[0][0]==0)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
	if(llr_list[1][0]==0)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
	if(llr_list[2][0]==0)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;

	//correction at the south pole
	if(llr_list[0][0]==180)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
	if(llr_list[1][0]==180)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
	if(llr_list[2][0]==180)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;

	late=(llr_list[0][0]+llr_list[1][0]+llr_list[2][0])/3.0;
	longe=(llr_list[0][1]+llr_list[1][1]+llr_list[2][1])/3.0;

	late=90-late;
	if(longe>180)longe=(longe-180)-180;

	late=late/180*PI;
	longe=longe/180*PI;
	/*}}}*/
	re=(llr_list[0][2]+llr_list[1][2]+llr_list[2][2])/3.0;

	if(IsWaterInElement()){
		IssmDouble rho_water, S;

		/*recover material parameters: */
		rho_water=FindParam(MaterialsRhoFreshwaterEnum);

		/*From Sg_old, recover water sea level rise:*/
		S=0; for(int i=0;i<NUMVERTICES;i++) S+=Sg_old[this->vertices[i]->Sid()]/NUMVERTICES;

		/* Perturbation terms for moment of inertia (moi_list):
		 * computed analytically (see Wu & Peltier, eqs 10 & 32)
		 * also consistent with my GMD formulation!
		 * ALL in geographic coordinates
		 * */
		dI_list[0] = -4*PI*(rho_water*S*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/eartharea;
		dI_list[1] = -4*PI*(rho_water*S*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/eartharea;
		dI_list[2] = +4*PI*(rho_water*S*area)*pow(re,4)*(1-pow(sin(late),2))/eartharea;
	}
	else if(IsIceOnlyInElement()){
		IssmDouble rho_ice, I;

		/*recover material parameters: */
		rho_ice=FindParam(MaterialsRhoIceEnum);

		/*Compute ice thickness change: */
		Input2* deltathickness_input=this->GetInput2(SealevelriseDeltathicknessEnum);
		if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!");
		deltathickness_input->GetInputAverage(&I);

		dI_list[0] = -4*PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*cos(longe))/eartharea;
		dI_list[1] = -4*PI*(rho_ice*I*area)*pow(re,4)*(sin(late)*cos(late)*sin(longe))/eartharea;
		dI_list[2] = +4*PI*(rho_ice*I*area)*pow(re,4)*(1-pow(sin(late),2))/eartharea;
	}

	return;
}/*}}}*/
void    Tria::SealevelriseEustatic(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/

	/*Computational flags:*/
	int bp_compute_fingerprints= 0;

	/*some paramters first: */
	this->parameters->FindParam(&bp_compute_fingerprints,DslComputeFingerprintsEnum);

	if(!IsOceanInElement()){
		/*ok, there is ocean in this element, we should compute eustatic loads for the ocean if we have requested
		 *bottom pressure fingerprints:*/
		if(bp_compute_fingerprints)this->SealevelriseEustaticBottomPressure(pSgi,peustatic,latitude,longitude,radius,oceanarea,eartharea);
	}
	//if(!IsIceInElement()){
		/*there is ice in this eleemnt, let's compute the eustatic response for ice changes:*/
		this->SealevelriseEustaticIce(pSgi,peustatic,latitude,longitude,radius,oceanarea,eartharea);
	//}

}
/*}}}*/
void    Tria::SealevelriseEustaticIce(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/

	/*diverse:*/
	int gsize;
	bool spherical=true;
	IssmDouble llr_list[NUMVERTICES][3];
	IssmDouble area;
	IssmDouble I;  //change in ice thickness or water level(Farrel and Clarke, Equ. 4)
	IssmDouble rho;
	IssmDouble late,longe,re;
	IssmDouble lati,longi,ri;
	bool notfullygrounded=false;

	/*elastic green function:*/
	IssmDouble* G_elastic_precomputed=NULL;
	int         M;

	/*ice properties: */
	IssmDouble rho_ice,rho_water,rho_earth;

	/*constants:*/
	IssmDouble constant=0;

	/*Initialize eustatic component: do not skip this step :):*/
	IssmDouble eustatic = 0.;

	/*Computational flags:*/
	bool computerigid = true;
	bool computeelastic= true;
	bool scaleoceanarea= false;

	/*early return if we are not on an ice cap:*/
	if(!IsIceOnlyInElement()){
		constant=0; this->AddInput2(SealevelEustaticMaskEnum,&constant,P0Enum);
		*peustatic=0; //do not forget to assign this pointer, otherwise, global eustatic will be garbage!
		return;
	}

	/*early return if we are fully floating:*/
	Input2* gr_input=this->GetInput2(MaskGroundediceLevelsetEnum); _assert_(gr_input);
	if (gr_input->GetInputMax()<=0){
		constant=0; this->AddInput2(SealevelEustaticMaskEnum,&constant,P0Enum);
		*peustatic=0; //do not forget to assign this pointer, otherwise, global eustatic will be garbage!
		return;
	}

	/*If we are here, we are on ice that is fully grounded or half-way to floating:*/
	if ((gr_input->GetInputMin())<0){
		notfullygrounded=true; //used later on.
	}

	/*Inform mask: */
	constant=1; this->AddInput2(SealevelEustaticMaskEnum,&constant,P0Enum);

	/*recover material parameters: */
	rho_ice=FindParam(MaterialsRhoIceEnum);
	rho_water=FindParam(MaterialsRhoFreshwaterEnum);
	rho_earth=FindParam(MaterialsEarthDensityEnum);

	/*recover love numbers and computational flags: */
	this->parameters->FindParam(&computerigid,SealevelriseRigidEnum);
	this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum);
	this->parameters->FindParam(&scaleoceanarea,SealevelriseOceanAreaScalingEnum);

	/*recover elastic green function:*/
	if(computeelastic){
		DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGElasticEnum));
		_assert_(parameter);
		parameter->GetParameterValueByPointer(&G_elastic_precomputed,&M);
	}

	/*how many dofs are we working with here? */
	this->parameters->FindParam(&gsize,MeshNumberofverticesEnum);

	/* Where is the centroid of this element?:{{{*/

	/*retrieve coordinates: */
	::GetVerticesCoordinates(&llr_list[0][0],this->vertices,NUMVERTICES,spherical);

	IssmDouble minlong=400;
	IssmDouble maxlong=-20;
	for (int i=0;i<NUMVERTICES;i++){
		llr_list[i][0]=(90-llr_list[i][0]);
		if(llr_list[i][1]<0)llr_list[i][1]=180+(180+llr_list[i][1]);
		if(llr_list[i][1]>maxlong)maxlong=llr_list[i][1];
		if(llr_list[i][1]<minlong)minlong=llr_list[i][1];
	}
	if(minlong==0 && maxlong>180){
		if (llr_list[0][1]==0)llr_list[0][1]=360;
		if (llr_list[1][1]==0)llr_list[1][1]=360;
		if (llr_list[2][1]==0)llr_list[2][1]=360;
	}

	// correction at the north pole
	if(llr_list[0][0]==0)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
	if(llr_list[1][0]==0)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
	if(llr_list[2][0]==0)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;

	//correction at the south pole
	if(llr_list[0][0]==180)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
	if(llr_list[1][0]==180)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
	if(llr_list[2][0]==180)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;

	late=(llr_list[0][0]+llr_list[1][0]+llr_list[2][0])/3.0;
	longe=(llr_list[0][1]+llr_list[1][1]+llr_list[2][1])/3.0;

	late=90-late;
	if(longe>180)longe=(longe-180)-180;

	late=late/180*PI;
	longe=longe/180*PI;
	/*}}}*/

	/*Compute area of element. Scale it by grounded fraction if not fully grounded: */
	area=GetAreaSpherical();
	if(notfullygrounded){
		IssmDouble phi=0;
		IssmDouble xyz_list[NUMVERTICES][3];
		::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);

		phi=this->GetGroundedPortion(&xyz_list[0][0]); //watch out, this only works because of the Thales theorem! We are in 3D, but this routine is inherently for 2D trias
		area*=phi;
	}

	/*Compute ice thickness: */
	Input2* deltathickness_input=this->GetInput2(SealevelriseDeltathicknessEnum);
	if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!");

	/*If we are fully grounded, take the average over the element: */
	if(!notfullygrounded)deltathickness_input->GetInputAverage(&I);
	else{
		IssmDouble total_weight=0;
		bool mainlyfloating = true;
		int         point1;
		IssmDouble  fraction1,fraction2;

		/*Recover portion of element that is grounded*/
		this->GetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating);
		Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,mainlyfloating,2);

		/* Start  looping on the number of gaussian points and average over these gaussian points: */
		total_weight=0;
		I=0;
		for(int ig=gauss->begin();ig<gauss->end();ig++){
			IssmDouble Ig=0;
			gauss->GaussPoint(ig);
			deltathickness_input->GetInputValue(&Ig,gauss);
			I+=Ig*gauss->weight;
			total_weight+=gauss->weight;
		}
		I=I/total_weight;
		delete gauss;
	}

	/*Compute eustatic compoent:*/
	_assert_(oceanarea>0.);
	if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2
	eustatic += rho_ice*area*I/(oceanarea*rho_water);

	if(computeelastic | computerigid){
		int* indices=xNew<int>(gsize);
		IssmDouble* values=xNew<IssmDouble>(gsize);
		IssmDouble alpha;
		IssmDouble delPhi,delLambda;
		for(int i=0;i<gsize;i++){
			indices[i]=i;

			IssmDouble G_rigid=0;  //do not remove =0!
			IssmDouble G_elastic=0;  //do not remove =0!

			/*Compute alpha angle between centroid and current vertex : */
			lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;

		   delPhi=fabs(lati-late); delLambda=fabs(longi-longe);
			alpha=2.*asin(sqrt(pow(sin(delPhi/2),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));

			//Rigid earth gravitational perturbation:
			if(computerigid)G_rigid=1.0/2.0/sin(alpha/2.0);

			//Elastic component  (from Eq 17 in Adhikari et al, GMD 2015)
			if(computeelastic){
				int index=reCast<int,IssmDouble>(alpha/PI*reCast<IssmDouble,int>(M-1));
				G_elastic += G_elastic_precomputed[index];
			}

			/*Add all components to the pSgi or pSgo solution vectors:*/
			values[i]=3*rho_ice/rho_earth*area/eartharea*I*(G_rigid+G_elastic);
		}
		pSgi->SetValues(gsize,indices,values,ADD_VAL);

		/*free ressources:*/
		xDelete<IssmDouble>(values);
		xDelete<int>(indices);
	}

	/*Assign output pointer:*/
	_assert_(!xIsNan<IssmDouble>(eustatic));
	_assert_(!xIsInf<IssmDouble>(eustatic));
	*peustatic=eustatic;
	return;
}
/*}}}*/
void    Tria::SealevelriseEustaticBottomPressure(Vector<IssmDouble>* pSgi,IssmDouble* peustatic,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble oceanarea,IssmDouble eartharea){ /*{{{*/

	/*diverse:*/
	int gsize;
	bool spherical=true;
	IssmDouble llr_list[NUMVERTICES][3];
	IssmDouble area;
	IssmDouble I;  //change in ice thickness or water level(Farrel and Clarke, Equ. 4)
	IssmDouble rho;
	IssmDouble late,longe,re;
	IssmDouble lati,longi,ri;

	/*elastic green function:*/
	IssmDouble* G_elastic_precomputed=NULL;
	int         M;

	/*ice properties: */
	IssmDouble rho_water,rho_earth;

	/*constants:*/
	IssmDouble constant=0;

	/*Initialize eustatic component: do not skip this step :):*/
	IssmDouble eustatic = 0.;

	/*Computational flags:*/
	bool computerigid = true;
	bool computeelastic= true;
	bool scaleoceanarea= false;
	bool bp_compute_fingerprints= false;

	/*we are here to compute fingerprints originating fromn bottom pressure loads:*/

	/*Inform mask: */
	constant=1; this->AddInput2(SealevelEustaticMaskEnum,&constant,P0Enum);

	/*recover material parameters: */
	rho_water=FindParam(MaterialsRhoFreshwaterEnum);
	rho_earth=FindParam(MaterialsEarthDensityEnum);

	/*recover love numbers and computational flags: */
	this->parameters->FindParam(&computerigid,SealevelriseRigidEnum);
	this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum);
	this->parameters->FindParam(&scaleoceanarea,SealevelriseOceanAreaScalingEnum);

	/*recover elastic green function:*/
	if(computeelastic){
		DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGElasticEnum));
		_assert_(parameter);
		parameter->GetParameterValueByPointer(&G_elastic_precomputed,&M);
	}

	/*how many dofs are we working with here? */
	this->parameters->FindParam(&gsize,MeshNumberofverticesEnum);

	/* Where is the centroid of this element?:{{{*/

	/*retrieve coordinates: */
	::GetVerticesCoordinates(&llr_list[0][0],this->vertices,NUMVERTICES,spherical);

	IssmDouble minlong=400;
	IssmDouble maxlong=-20;
	for (int i=0;i<NUMVERTICES;i++){
		llr_list[i][0]=(90-llr_list[i][0]);
		if(llr_list[i][1]<0)llr_list[i][1]=180+(180+llr_list[i][1]);
		if(llr_list[i][1]>maxlong)maxlong=llr_list[i][1];
		if(llr_list[i][1]<minlong)minlong=llr_list[i][1];
	}
	if(minlong==0 && maxlong>180){
		if (llr_list[0][1]==0)llr_list[0][1]=360;
		if (llr_list[1][1]==0)llr_list[1][1]=360;
		if (llr_list[2][1]==0)llr_list[2][1]=360;
	}

	// correction at the north pole
	if(llr_list[0][0]==0)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
	if(llr_list[1][0]==0)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
	if(llr_list[2][0]==0)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;

	//correction at the south pole
	if(llr_list[0][0]==180)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
	if(llr_list[1][0]==180)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
	if(llr_list[2][0]==180)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;

	late=(llr_list[0][0]+llr_list[1][0]+llr_list[2][0])/3.0;
	longe=(llr_list[0][1]+llr_list[1][1]+llr_list[2][1])/3.0;

	late=90-late;
	if(longe>180)longe=(longe-180)-180;

	late=late/180*PI;
	longe=longe/180*PI;
	/*}}}*/

	/*Compute area of element. For now, we dont do partially grounded elements:*/
	area=GetAreaSpherical();

	/*Compute bottom pressure change: */
	Input2* bottompressure_change_input=this->GetInput2(DslSeaWaterPressureChangeAtSeaFloor);
	if (!bottompressure_change_input)_error_("bottom pressure input needed to compute sea level rise fingerprint!");

	/*If we are fully grounded, take the average over the element: */
	bottompressure_change_input->GetInputAverage(&I);

	/*Compute eustatic compoent:*/
	_assert_(oceanarea>0.);
	if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2

	/*We do not need to add the bottom pressure component to the eustatic value: */
	eustatic += 0;

	if(computeelastic | computerigid){
		int* indices=xNew<int>(gsize);
		IssmDouble* values=xNew<IssmDouble>(gsize);
		IssmDouble alpha;
		IssmDouble delPhi,delLambda;
		for(int i=0;i<gsize;i++){
			indices[i]=i;

			IssmDouble G_rigid=0;  //do not remove =0!
			IssmDouble G_elastic=0;  //do not remove =0!

			/*Compute alpha angle between centroid and current vertex : */
			lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;

		   delPhi=fabs(lati-late); delLambda=fabs(longi-longe);
			alpha=2.*asin(sqrt(pow(sin(delPhi/2),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));

			//Rigid earth gravitational perturbation:
			if(computerigid)G_rigid=1.0/2.0/sin(alpha/2.0);

			//Elastic component  (from Eq 17 in Adhikari et al, GMD 2015)
			if(computeelastic){
				int index=reCast<int,IssmDouble>(alpha/PI*reCast<IssmDouble,int>(M-1));
				G_elastic += G_elastic_precomputed[index];
			}

			/*Add all components to the pSgi or pSgo solution vectors:*/
			values[i]=3*rho_water/rho_earth*area/eartharea*I*(G_rigid+G_elastic);
		}
		pSgi->SetValues(gsize,indices,values,ADD_VAL);

		/*free ressources:*/
		xDelete<IssmDouble>(values);
		xDelete<int>(indices);
	}

	/*Assign output pointer:*/
	_assert_(!xIsNan<IssmDouble>(eustatic));
	_assert_(!xIsInf<IssmDouble>(eustatic));
	*peustatic=eustatic;
	return;
}
/*}}}*/
void    Tria::SealevelriseNonEustatic(Vector<IssmDouble>* pSgo,IssmDouble* Sg_old,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble eartharea){ /*{{{*/

	/*diverse:*/
	int gsize;
	bool spherical=true;
	IssmDouble llr_list[NUMVERTICES][3];
	IssmDouble area;
	IssmDouble S;  //change in water water level(Farrel and Clarke, Equ. 4)
	IssmDouble late,longe;
	IssmDouble lati,longi,ri;
	IssmDouble minlong=400;
	IssmDouble maxlong=-20;
	IssmDouble constant=0;

	/*precomputed elastic green functions:*/
	IssmDouble* G_elastic_precomputed = NULL;
	int         M;

	/*computation of Green functions:*/
	IssmDouble* G_elastic= NULL;
	IssmDouble* G_rigid= NULL;

	/*optimization:*/
	bool store_green_functions=false;

	/*ice properties: */
	IssmDouble rho_ice,rho_water,rho_earth;

	/*Computational flags:*/
	bool computerigid = true;
	bool computeelastic= true;

	/*early return if we are not on the ocean:*/
	if (!IsWaterInElement()){
		constant=0; this->AddInput2(SealevelEustaticOceanMaskEnum,&constant,P0Enum);
		return;
	}
	constant=1; this->AddInput2(SealevelEustaticOceanMaskEnum,&constant,P0Enum);

	/*recover computational flags: */
	this->parameters->FindParam(&computerigid,SealevelriseRigidEnum);
	this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum);

	/*early return if rigid or elastic not requested:*/
	if(!computerigid && !computeelastic) return;

	/*recover material parameters: */
	rho_water=FindParam(MaterialsRhoFreshwaterEnum);
	rho_earth=FindParam(MaterialsEarthDensityEnum);

	/*how many dofs are we working with here? */
	this->parameters->FindParam(&gsize,MeshNumberofverticesEnum);

	/*From Sg_old, recover water sea level rise:*/
	S=0; for(int i=0;i<NUMVERTICES;i++) S+=Sg_old[this->vertices[i]->Sid()]/NUMVERTICES;

	/*Compute area of element:*/
	area=GetAreaSpherical();

	/* Where is the centroid of this element?:{{{*/
	::GetVerticesCoordinates(&llr_list[0][0],this->vertices,NUMVERTICES,spherical);

	minlong=400; maxlong=-20;
	for (int i=0;i<NUMVERTICES;i++){
		llr_list[i][0]=(90-llr_list[i][0]);
		if(llr_list[i][1]<0)llr_list[i][1]=180+(180+llr_list[i][1]);
		if(llr_list[i][1]>maxlong)maxlong=llr_list[i][1];
		if(llr_list[i][1]<minlong)minlong=llr_list[i][1];
	}
	if(minlong==0 && maxlong>180){
		if (llr_list[0][1]==0)llr_list[0][1]=360;
		if (llr_list[1][1]==0)llr_list[1][1]=360;
		if (llr_list[2][1]==0)llr_list[2][1]=360;
	}

	// correction at the north pole
	if(llr_list[0][0]==0)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
	if(llr_list[1][0]==0)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
	if(llr_list[2][0]==0)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;

	//correction at the south pole
	if(llr_list[0][0]==180)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
	if(llr_list[1][0]==180)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
	if(llr_list[2][0]==180)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;

	late=(llr_list[0][0]+llr_list[1][0]+llr_list[2][0])/3.0;
	longe=(llr_list[0][1]+llr_list[1][1]+llr_list[2][1])/3.0;

	late=90-late;
	if(longe>180)longe=(longe-180)-180;

	late=late/180*PI;
	longe=longe/180*PI;
	/*}}}*/

	if(computeelastic){

		/*recover elastic green function:*/
		DoubleVecParam* parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseGElasticEnum)); _assert_(parameter);
		parameter->GetParameterValueByPointer(&G_elastic_precomputed,&M);

		/*initialize G_elastic:*/
		G_elastic=xNewZeroInit<IssmDouble>(gsize);
	}
	if(computerigid) G_rigid=xNewZeroInit<IssmDouble>(gsize);

	int* indices=xNew<int>(gsize);
	IssmDouble* values=xNewZeroInit<IssmDouble>(gsize);
	IssmDouble alpha;
	IssmDouble delPhi,delLambda;

	for(int i=0;i<gsize;i++){

		indices[i]=i;

		/*Compute alpha angle between centroid and current vertex : */
		lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;

		delPhi=fabs(lati-late); delLambda=fabs(longi-longe);
		alpha=2.*asin(sqrt(pow(sin(delPhi/2),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));

		/*Rigid earth gravitational perturbation: */
		if(computerigid){
			G_rigid[i]=1.0/2.0/sin(alpha/2.0);
			values[i]+=3*rho_water/rho_earth*area/eartharea*S*G_rigid[i];
		}

		/*Elastic component  (from Eq 17 in Adhikari et al, GMD 2015): */
		if(computeelastic){
			int index=reCast<int,IssmDouble>(alpha/PI*(M-1));
			G_elastic[i] += G_elastic_precomputed[index];
			values[i]+=3*rho_water/rho_earth*area/eartharea*S*G_elastic[i];
		}
	}

	pSgo->SetValues(gsize,indices,values,ADD_VAL);

	/*free ressources:*/
	xDelete<IssmDouble>(values);
	xDelete<int>(indices);

	/*Free ressources:*/
	if(computeelastic) xDelete<IssmDouble>(G_elastic);
	if(computerigid) xDelete<IssmDouble>(G_rigid);

	return;
}
/*}}}*/
void    Tria::SealevelriseGeodetic(Vector<IssmDouble>* pUp,Vector<IssmDouble>* pNorth,Vector<IssmDouble>* pEast,IssmDouble* Sg,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz,IssmDouble eartharea,int horiz){ /*{{{*/

	/*diverse:*/
	int gsize;
	bool spherical=true;
	IssmDouble llr_list[NUMVERTICES][3];
	IssmDouble xyz_list[NUMVERTICES][3];
	IssmDouble area;
	IssmDouble I, S;		//change in relative ice thickness and sea level
	IssmDouble late,longe,re;
	IssmDouble lati,longi,ri;
	IssmDouble rho_ice,rho_water,rho_earth;
	IssmDouble minlong=400;
	IssmDouble maxlong=-20;

	/*precomputed elastic green functions:*/
	IssmDouble* U_elastic_precomputed = NULL;
	IssmDouble* H_elastic_precomputed = NULL;
	int         M;

	/*computation of Green functions:*/
	IssmDouble* U_elastic= NULL;
	IssmDouble* N_elastic= NULL;
	IssmDouble* E_elastic= NULL;
	DoubleVecParam* U_parameter = NULL;
	DoubleVecParam* H_parameter = NULL;
	IssmDouble* U_values=NULL;
	IssmDouble* N_values=NULL;
	IssmDouble* E_values=NULL;

	/*optimization:*/
	bool store_green_functions=false;

	/*computational flags:*/
	bool computerigid = true;
	bool computeelastic= true;

	/*early return if we are not on the ocean or on an ice cap:*/
	if(!IsIceOnlyInElement() && !IsWaterInElement()) return;

	/*early return if we are fully floating: */
	Input2* gr_input=this->GetInput2(MaskGroundediceLevelsetEnum); _assert_(gr_input);
	if(gr_input->GetInputMax()<=0)return;

	/*recover computational flags: */
	this->parameters->FindParam(&computerigid,SealevelriseRigidEnum);
	this->parameters->FindParam(&computeelastic,SealevelriseElasticEnum);

	/*early return if elastic not requested:*/
	if(!computeelastic) return;

	/*recover material parameters: */
	rho_ice=FindParam(MaterialsRhoIceEnum);
	rho_water=FindParam(MaterialsRhoFreshwaterEnum);
	rho_earth=FindParam(MaterialsEarthDensityEnum);

	/*how many dofs are we working with here? */
	this->parameters->FindParam(&gsize,MeshNumberofverticesEnum);

	/*compute area of element:*/
	area=GetAreaSpherical();

	/*element centroid (spherical): */
	/* Where is the centroid of this element?:{{{*/
	::GetVerticesCoordinates(&llr_list[0][0],this->vertices,NUMVERTICES,spherical);

	minlong=400; maxlong=-20;
	for (int i=0;i<NUMVERTICES;i++){
		llr_list[i][0]=(90-llr_list[i][0]);
		if(llr_list[i][1]<0)llr_list[i][1]=180+(180+llr_list[i][1]);
		if(llr_list[i][1]>maxlong)maxlong=llr_list[i][1];
		if(llr_list[i][1]<minlong)minlong=llr_list[i][1];
	}
	if(minlong==0 && maxlong>180){
		if (llr_list[0][1]==0)llr_list[0][1]=360;
		if (llr_list[1][1]==0)llr_list[1][1]=360;
		if (llr_list[2][1]==0)llr_list[2][1]=360;
	}

	// correction at the north pole
	if(llr_list[0][0]==0)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
	if(llr_list[1][0]==0)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
	if(llr_list[2][0]==0)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;

	//correction at the south pole
	if(llr_list[0][0]==180)llr_list[0][1]=(llr_list[1][1]+llr_list[2][1])/2.0;
	if(llr_list[1][0]==180)llr_list[1][1]=(llr_list[0][1]+llr_list[2][1])/2.0;
	if(llr_list[2][0]==180)llr_list[2][1]=(llr_list[0][1]+llr_list[1][1])/2.0;

	late=(llr_list[0][0]+llr_list[1][0]+llr_list[2][0])/3.0;
	longe=(llr_list[0][1]+llr_list[1][1]+llr_list[2][1])/3.0;

	late=90-late;
	if(longe>180)longe=(longe-180)-180;

	late=late/180*PI;
	longe=longe/180*PI;
	/*}}}*/

	/*figure out gravity center of our element (Cartesian): */
	IssmDouble x_element, y_element, z_element;
	::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
	x_element=(xyz_list[0][0]+xyz_list[1][0]+xyz_list[2][0])/3.0;
	y_element=(xyz_list[0][1]+xyz_list[1][1]+xyz_list[2][1])/3.0;
	z_element=(xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2])/3.0;

	/*recover elastic Green's functions for displacement:*/
	U_parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(U_parameter);
	U_parameter->GetParameterValueByPointer(&U_elastic_precomputed,&M);
	if(horiz){
		H_parameter = static_cast<DoubleVecParam*>(this->parameters->FindParamObject(SealevelriseHElasticEnum)); _assert_(H_parameter);
		H_parameter->GetParameterValueByPointer(&H_elastic_precomputed,&M);
	}

	/*From Sg, recover water sea level rise:*/
	S=0; for(int i=0;i<NUMVERTICES;i++) S+=Sg[this->vertices[i]->Sid()]/NUMVERTICES;

	/*Compute ice thickness change: */
	Input2* deltathickness_input=this->GetInput2(SealevelriseDeltathicknessEnum);
	if (!deltathickness_input)_error_("delta thickness input needed to compute sea level rise!");
	deltathickness_input->GetInputAverage(&I);

	/*initialize: */
	U_elastic=xNewZeroInit<IssmDouble>(gsize);
	if(horiz){
		N_elastic=xNewZeroInit<IssmDouble>(gsize);
		E_elastic=xNewZeroInit<IssmDouble>(gsize);
	}

	int* indices=xNew<int>(gsize);
	U_values=xNewZeroInit<IssmDouble>(gsize);
	if(horiz){
		N_values=xNewZeroInit<IssmDouble>(gsize);
		E_values=xNewZeroInit<IssmDouble>(gsize);
	}
	IssmDouble alpha;
	IssmDouble delPhi,delLambda;
	IssmDouble dx, dy, dz, x, y, z;
	IssmDouble N_azim, E_azim;

	for(int i=0;i<gsize;i++){

		indices[i]=i;

		/*Compute alpha angle between centroid and current vertex: */
		lati=latitude[i]/180*PI; longi=longitude[i]/180*PI;

		delPhi=fabs(lati-late); delLambda=fabs(longi-longe);
		alpha=2.*asin(sqrt(pow(sin(delPhi/2),2.0)+cos(lati)*cos(late)*pow(sin(delLambda/2),2)));

		/*Compute azimuths, both north and east components: */
		x = xx[i]; y = yy[i]; z = zz[i];
		if(latitude[i]==90){
			x=1e-12; y=1e-12;
		}
		if(latitude[i]==-90){
			x=1e-12; y=1e-12;
		}
		dx = x_element-x; dy = y_element-y; dz = z_element-z;
		if(horiz){
			N_azim = (-z*x*dx-z*y*dy+(pow(x,2)+pow(y,2))*dz) /pow((pow(x,2)+pow(y,2))*(pow(x,2)+pow(y,2)+pow(z,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
			E_azim = (-y*dx+x*dy) /pow((pow(x,2)+pow(y,2))*(pow(dx,2)+pow(dy,2)+pow(dz,2)),0.5);
		}

		/*Elastic component  (from Eq 17 in Adhikari et al, GMD 2015): */
		int index=reCast<int,IssmDouble>(alpha/PI*(M-1));
		U_elastic[i] += U_elastic_precomputed[index];
		if(horiz){
			N_elastic[i] += H_elastic_precomputed[index]*N_azim;
			E_elastic[i] += H_elastic_precomputed[index]*E_azim;
		}

		/*Add all components to the pUp solution vectors:*/
		if(IsIceOnlyInElement()){
			U_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*U_elastic[i];
			if(horiz){
				N_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*N_elastic[i];
				E_values[i]+=3*rho_ice/rho_earth*area/eartharea*I*E_elastic[i];
			}
		}
		else if(IsWaterInElement()) {
			U_values[i]+=3*rho_water/rho_earth*area/eartharea*S*U_elastic[i];
			if(horiz){
				N_values[i]+=3*rho_water/rho_earth*area/eartharea*S*N_elastic[i];
				E_values[i]+=3*rho_water/rho_earth*area/eartharea*S*E_elastic[i];
			}
		}
	}
	pUp->SetValues(gsize,indices,U_values,ADD_VAL);
	if(horiz){
		pNorth->SetValues(gsize,indices,N_values,ADD_VAL);
		pEast->SetValues(gsize,indices,E_values,ADD_VAL);
	}

	/*free ressources:*/
	xDelete<int>(indices);
	xDelete<IssmDouble>(U_values);
	xDelete<IssmDouble>(U_elastic);
	if(horiz){
		xDelete<IssmDouble>(N_values); xDelete<IssmDouble>(E_values);
		xDelete<IssmDouble>(N_elastic); xDelete<IssmDouble>(E_elastic);
	}

	return;
}
/*}}}*/
#endif

#ifdef _HAVE_DAKOTA_
void       Tria::InputUpdateFromMatrixDakota(IssmDouble* matrix, int nrows, int ncols, int name, int type){/*{{{*/

	/*Check that name is an element input*/
	if(!IsInputEnum(name)) _error_("Enum "<<EnumToStringx(name)<<" is not in IsInput");
	TransientInput2* transientinput = inputs2->GetTransientInput(name);

	switch(type){

		case VertexEnum:

			/*Get LID lists once for all*/
			IssmDouble  values[NUMVERTICES];
			int         lidlist[NUMVERTICES];
			this->GetVerticesLidList(&lidlist[0]);

			/*Create transient input: */
			for(int t=0;t<ncols;t++){ //ncols is the number of times
				for(int i=0;i<3;i++){
					int row=this->vertices[i]->Sid();
					values[i]=matrix[ncols*row+t];
				}

				/*time:*/
				IssmDouble time=matrix[(nrows-1)*ncols+t];

				transientinput->AddTriaTimeInput(t,NUMVERTICES,&lidlist[0],&values[0],P1Enum);
			}
			break;

		case ElementEnum:
			/*Get value for the element: */
			for(int t=0;t<ncols;t++){ //ncols is the number of times
				IssmDouble value=matrix[ncols*(this->Sid())+t];
				IssmDouble time=matrix[(nrows-1)*ncols+t];
				transientinput->AddTriaTimeInput(t,1,&(this->lid),&value,P0Enum);
			}
			break;

		default:
			_error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
	}
}
/*}}}*/
void       Tria::InputUpdateFromVectorDakota(IssmDouble* vector, int name, int type){/*{{{*/

	int i,j;

	/*Check that name is an element input*/
	if(!IsInputEnum(name)) _error_("Enum "<<EnumToStringx(name)<<" is not in IsInput");

	switch(type){

		case VertexEnum:

			/*New TriaInput*/
			IssmDouble values[3];

			/*Get values on the 3 vertices*/
			for (i=0;i<3;i++){
				values[i]=vector[this->vertices[i]->Sid()]; //careful, vector of values here is not parallel distributed, but serial distributed (from a serial Dakota core!)
			}

			/*Branch on the specified type of update: */
			switch(name){
				case ThicknessEnum:
					IssmDouble  thickness[3];
					IssmDouble  thickness_init[3];
					IssmDouble  hydrostatic_ratio[3];
					IssmDouble  surface[3];
					IssmDouble  bed[3];

					/*retrieve inputs: */
					Element::GetInputListOnVertices(&thickness_init[0],ThicknessEnum);
					Element::GetInputListOnVertices(&hydrostatic_ratio[0],GeometryHydrostaticRatioEnum);
					Element::GetInputListOnVertices(&bed[0],BaseEnum);
					Element::GetInputListOnVertices(&surface[0],SurfaceEnum);

					/*build new bed and surface: */
					if (this->IsFloating()){
						/*hydrostatic equilibrium: */
						IssmDouble rho_ice,rho_water,di;
						rho_ice   = this->FindParam(MaterialsRhoIceEnum);
						rho_water = this->FindParam(MaterialsRhoSeawaterEnum);
						di        = rho_ice/rho_water;

						/*build new thickness: */
						for (j=0; j<3; j++) {
							/*  for observed/interpolated/hydrostatic thickness, remove scaling from any hydrostatic thickness  */
							if (hydrostatic_ratio[j] >= 0.)
								thickness[j]=values[j]-(values[j]/thickness_init[j]-1.)*hydrostatic_ratio[j]*surface[j]/(1.-di);
							/*  for minimum thickness, don't scale  */
							else
								thickness[j]=thickness_init[j];

							/*  check the computed thickness and update bed*/
							if (thickness[j] < 0.) thickness[j]=1./(1.-di);
							bed[j]=surface[j]-thickness[j];
						}
					}
					else{
						/*build new thickness: */
						for (j=0; j<3; j++) {
							/*  for observed thickness, use scaled value  */
							if (hydrostatic_ratio[j] >= 0.)
								thickness[j]=values[j];
							/*  for minimum thickness, don't scale  */
							else
								thickness[j]=thickness_init[j];
						}

						/*update bed on grounded ice: */
						for(j=0;j<3;j++)bed[j]=surface[j]-thickness[j];
					}

					/*Add new inputs: */
					this->AddInput2(ThicknessEnum,thickness,P1Enum);
					this->AddInput2(BaseEnum,bed,P1Enum);
					this->AddInput2(SurfaceEnum,surface,P1Enum);

					break;
				case MaterialsRheologyBEnum:
					this->AddInput2(MaterialsRheologyBbarEnum,values,P1Enum);
					break;
				default:
					this->AddInput2(name,values,P1Enum);
			}
			break;

		case ElementEnum:
			IssmDouble value;
			/*Get value for the element: */
			value=vector[this->Sid()]; //careful, vector of values here is not parallel distributed, but serial distributed (from a serial Dakota core!)
			this->AddInput2(name,&value,P0Enum);
			break;
		default:
			_error_("type " << type << " (" << EnumToStringx(type) << ") not implemented yet");
	}

}
/*}}}*/
#endif
