/*!\file Tria.cpp
 * \brief: implementation of the Tria object
 */

/*Headers:*/
/*{{{1*/
#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 "../objects.h"
#include "../../shared/shared.h"
#include "../../Container/Container.h"
#include "../../include/include.h"
/*}}}*/

/*Element macros*/
#define NUMVERTICES 3

/*Constructors/destructor/copy*/
/*FUNCTION Tria::Tria(){{{1*/
Tria::Tria(){

	int i;

	this->nodes=NULL;
	this->matice=NULL;
	this->matpar=NULL;
	for(i=0;i<3;i++)this->horizontalneighborsids[i]=UNDEF;
	this->inputs=NULL;
	this->parameters=NULL;
	this->results=NULL;

}
/*}}}*/
/*FUNCTION Tria::Tria(int id, int sid,int index, IoModel* iomodel,int nummodels){{{1*/
Tria::Tria(int tria_id, int tria_sid, int index, IoModel* iomodel,int nummodels)
	:TriaRef(nummodels)
	,TriaHook(nummodels,index+1,iomodel){
		
		int i;
		/*id: */
		this->id=tria_id;
		this->sid=tria_sid;

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

		/*Build horizontalneighborsids list: */
		_assert_(iomodel->Data(MeshElementconnectivityEnum));
		//for (i=0;i<3;i++) this->horizontalneighborsids[i]=(int)iomodel->elementconnectivity[3*index+i]-1;

		/*intialize inputs and results: */
		this->inputs=new Inputs();
		this->results=new Results();

		/*initialize pointers:*/
		this->nodes=NULL;
		this->matice=NULL;
		this->matpar=NULL;

}
/*}}}*/
/*FUNCTION Tria::~Tria(){{{1*/
Tria::~Tria(){
	delete inputs;
	delete results;
	this->parameters=NULL;
}
/*}}}*/
/*FUNCTION Tria::copy {{{1*/
Object* Tria::copy() {

	int i;
	Tria* tria=NULL;

	tria=new Tria();

	//deal with TriaRef mother class
	tria->element_type_list=(int*)xmalloc(this->numanalyses*sizeof(int));
	for(i=0;i<this->numanalyses;i++) tria->element_type_list[i]=this->element_type_list[i];

	//deal with TriaHook mother class
	tria->numanalyses=this->numanalyses;
	tria->hnodes=new Hook*[tria->numanalyses];
	for(i=0;i<tria->numanalyses;i++)tria->hnodes[i]=(Hook*)this->hnodes[i]->copy();
	tria->hmatice=(Hook*)this->hmatice->copy();
	tria->hmatpar=(Hook*)this->hmatpar->copy();

	/*deal with Tria fields: */
	tria->id=this->id;
	tria->sid=this->sid;
	if(this->inputs){
		tria->inputs=(Inputs*)this->inputs->Copy();
	}
	else{
		tria->inputs=new Inputs();
	}
	if(this->results){
		tria->results=(Results*)this->results->Copy();
	}
	else{
		tria->results=new Results();
	}
	/*point parameters: */
	tria->parameters=this->parameters;

	/*recover objects: */
	tria->nodes=(Node**)xmalloc(3*sizeof(Node*)); //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<3;i++)tria->nodes[i]=this->nodes[i];
	tria->matice=(Matice*)tria->hmatice->delivers();
	tria->matpar=(Matpar*)tria->hmatpar->delivers();

	/*neighbors: */
	for(i=0;i<3;i++)tria->horizontalneighborsids[i]=this->horizontalneighborsids[i];

	return tria;
}
/*}}}*/

/*Marshall*/
#ifdef _SERIAL_
/*FUNCTION Tria::Marshall {{{1*/
void  Tria::Marshall(char** pmarshalled_dataset){

	int   i;
	char* marshalled_dataset=NULL;
	int   enum_type=0;
	char* marshalled_inputs=NULL;
	int   marshalled_inputs_size;
	char* marshalled_results=NULL;
	int   marshalled_results_size;
	int   flaghook; //to indicate if hook is NULL or exists

	/*recover marshalled_dataset: */
	marshalled_dataset=*pmarshalled_dataset;

	/*get enum type of Tria: */
	enum_type=TriaEnum;

	/*marshall enum: */
	memcpy(marshalled_dataset,&enum_type,sizeof(enum_type));marshalled_dataset+=sizeof(enum_type);

	/*marshall Tria data: */
	memcpy(marshalled_dataset,&id,sizeof(id));marshalled_dataset+=sizeof(id);
	memcpy(marshalled_dataset,&sid,sizeof(sid));marshalled_dataset+=sizeof(sid);
	memcpy(marshalled_dataset,&numanalyses,sizeof(numanalyses));marshalled_dataset+=sizeof(numanalyses);

	/*Mershall Ref: */
	for(i=0;i<numanalyses;i++){
		memcpy(marshalled_dataset,&element_type_list[i],sizeof(element_type_list[i]));marshalled_dataset+=sizeof(element_type_list[i]);
	}

	/*Marshall hooks: */
	for(i=0;i<numanalyses;i++){
		if(hnodes[i]){
			/*Set flag to 1 as there is a hook */
			flaghook=1;
			memcpy(marshalled_dataset,&flaghook,sizeof(flaghook));marshalled_dataset+=sizeof(flaghook);
			hnodes[i]->Marshall(&marshalled_dataset);
		}
		else{
			/*Set flag to 0 and do not marshall flag as there is no Hook */
			flaghook=0;
			memcpy(marshalled_dataset,&flaghook,sizeof(flaghook));marshalled_dataset+=sizeof(flaghook);
		}
	}
	hmatice->Marshall(&marshalled_dataset);
	hmatpar->Marshall(&marshalled_dataset);

	/*Marshall inputs: */
	marshalled_inputs_size=inputs->MarshallSize();
	marshalled_inputs=inputs->Marshall();
	memcpy(marshalled_dataset,marshalled_inputs,marshalled_inputs_size*sizeof(char));
	marshalled_dataset+=marshalled_inputs_size;

	/*Marshall results: */
	marshalled_results_size=results->MarshallSize();
	marshalled_results=results->Marshall();
	memcpy(marshalled_dataset,marshalled_results,marshalled_results_size*sizeof(char));
	marshalled_dataset+=marshalled_results_size;

	/*parameters: don't do anything about it. parameters are marshalled somewhere else!*/

	xfree((void**)&marshalled_inputs);
	xfree((void**)&marshalled_results);

	/*marshall horizontal neighbors: */
	memcpy(marshalled_dataset,horizontalneighborsids,3*sizeof(int));marshalled_dataset+=3*sizeof(int);

	*pmarshalled_dataset=marshalled_dataset;
	return;
}
/*}}}*/
/*FUNCTION Tria::MarshallSize {{{1*/
int   Tria::MarshallSize(){

	int i;
	int hnodes_size=0;;

	for(i=0;i<numanalyses;i++){
		hnodes_size+=sizeof(int); //Flag 0 or 1
		if (hnodes[i]) hnodes_size+=hnodes[i]->MarshallSize();
	}

	return sizeof(id)
	  +sizeof(sid)
	  +hnodes_size
	  +sizeof(numanalyses)
	  +numanalyses*sizeof(int) //element_type_lists
	  +hmatice->MarshallSize()
	  +hmatpar->MarshallSize()
	  +inputs->MarshallSize()
	  +results->MarshallSize()
	  +3*sizeof(int)
	  +sizeof(int); //sizeof(int) for enum type
}
/*}}}*/
/*FUNCTION Tria::Demarshall {{{1*/
void  Tria::Demarshall(char** pmarshalled_dataset){

	char* marshalled_dataset=NULL;
	int i;
	int flaghook;

	/*recover marshalled_dataset: */
	marshalled_dataset=*pmarshalled_dataset;

	/*this time, no need to get enum type, the pointer directly points to the beginning of the 
	 *object data (thanks to DataSet::Demarshall):*/
	memcpy(&id,marshalled_dataset,sizeof(id));marshalled_dataset+=sizeof(id);
	memcpy(&sid,marshalled_dataset,sizeof(sid));marshalled_dataset+=sizeof(sid);
	memcpy(&numanalyses,marshalled_dataset,sizeof(numanalyses));marshalled_dataset+=sizeof(numanalyses);

	/*demarshall Ref: */
	this->element_type_list=(int*)xmalloc(this->numanalyses*sizeof(int));
	for(i=0;i<numanalyses;i++){ memcpy(&element_type_list[i],marshalled_dataset,sizeof(int));marshalled_dataset+=sizeof(int);}

	/*allocate dynamic memory: */
	this->hnodes=new Hook*[this->numanalyses];
	/*demarshall hooks: */
	for(i=0;i<numanalyses;i++){
		memcpy(&flaghook,marshalled_dataset,sizeof(flaghook));marshalled_dataset+=sizeof(flaghook);
		if(flaghook){ // there is a hook so demarshall it
			hnodes[i]=new Hook();
			hnodes[i]->Demarshall(&marshalled_dataset);
		}
		else hnodes[i]=NULL; //There is no hook so it is NULL
	}
	hmatice=new Hook(); hmatice->Demarshall(&marshalled_dataset);
	hmatpar=new Hook(); hmatpar->Demarshall(&marshalled_dataset);

	/*pointers are garbabe, until configuration is carried out: */
	nodes=NULL;
	matice=NULL;
	matpar=NULL;
	
	/*demarshall inputs: */
	inputs=(Inputs*)DataSetDemarshallRaw(&marshalled_dataset); 
	results=(Results*)DataSetDemarshallRaw(&marshalled_dataset); 

	/*parameters: may not exist even yet, so let Configure handle it: */
	this->parameters=NULL;

	/*neighbors: */
	memcpy(&this->horizontalneighborsids,marshalled_dataset,3*sizeof(int));marshalled_dataset+=3*sizeof(int);

	/*return: */
	*pmarshalled_dataset=marshalled_dataset;
	return;
}
/*}}}*/
#endif

/*Other*/
/*FUNCTION Tria::AverageOntoPartition {{{1*/
void  Tria::AverageOntoPartition(Vec partition_contributions,Vec partition_areas,double* vertex_response,double* qmu_part){

	bool      already=false;
	int       i,j;
	int       partition[NUMVERTICES];
	int       offsetsid[NUMVERTICES];
	int       offsetdof[NUMVERTICES];
	double    area;
	double    mean;
	double    values[3];

	/*First, get the area: */
	area=this->GetArea();

	/*Figure out the average for this element: */
	this->GetSidList(&offsetsid[0]);
	this->GetDofList1(&offsetdof[0]);
	mean=0;
	for(i=0;i<NUMVERTICES;i++){
		partition[i]=(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){
			VecSetValue(partition_contributions,partition[i],mean*area,ADD_VALUES);
			VecSetValue(partition_areas,partition[i],area,ADD_VALUES);
		};
	}
}
/*}}}*/
/*FUNCTION Tria::CreateKMatrix {{{1*/
void  Tria::CreateKMatrix(Mat Kff, Mat Kfs,Vec df){

	/*retreive parameters: */
	ElementMatrix* Ke=NULL;
	int analysis_type;
	parameters->FindParam(&analysis_type,AnalysisTypeEnum);

	/*Checks in debugging mode{{{2*/
	_assert_(this->nodes && this->matice && this->matpar && this->parameters && this->inputs);
	/*}}}*/
	
	/*Skip if water element*/
	if(IsOnWater()) return;

	/*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
	switch(analysis_type){
		#ifdef _HAVE_DIAGNOSTIC_
		case DiagnosticHorizAnalysisEnum:
			Ke=CreateKMatrixDiagnosticMacAyeal();
			break;
		case AdjointHorizAnalysisEnum:
			Ke=CreateKMatrixAdjointMacAyeal();
			break;
		case DiagnosticHutterAnalysisEnum:
			Ke=CreateKMatrixDiagnosticHutter();
			break;
		 #endif
		case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum:
			Ke=CreateKMatrixSlope();
			break;
		case PrognosticAnalysisEnum:
			Ke=CreateKMatrixPrognostic();
			break;
		#ifdef _HAVE_HYDROLOGY_
		case HydrologyAnalysisEnum:
			Ke=CreateKMatrixHydrology();
			break;
		#endif
		#ifdef _HAVE_BALANCED_
		case BalancethicknessAnalysisEnum:
			Ke=CreateKMatrixBalancethickness();
			break;
		#endif
		#ifdef _HAVE_CONTROL_
		case AdjointBalancethicknessAnalysisEnum:
			Ke=CreateKMatrixAdjointBalancethickness();
			break;
		#endif
		default:
			_error_("analysis %i (%s) not supported yet",analysis_type,EnumToStringx(analysis_type));
	}

	/*Add to global matrix*/
	if(Ke){
		Ke->AddToGlobal(Kff,Kfs);
		delete Ke;
	}
}
/*}}}*/
/*FUNCTION Tria::CreateKMatrixMelting {{{1*/
ElementMatrix* Tria::CreateKMatrixMelting(void){

	/*Constants*/
	const int  numdof=NUMVERTICES*NDOF1;

	/*Intermediaries */
	int        i,j,ig;
	double     heatcapacity,latentheat;
	double     Jdet,D_scalar;
	double     xyz_list[NUMVERTICES][3];
	double     L[3];
	GaussTria *gauss=NULL;

	/*Initialize Element matrix*/
	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);

	/*Retrieve all inputs and parameters*/
	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
	latentheat=matpar->GetLatentHeat();
	heatcapacity=matpar->GetHeatCapacity();

	/* Start looping on the number of gauss  (nodes on the bedrock) */
	gauss=new GaussTria(2);
	for (ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0], gauss);

		D_scalar=latentheat/heatcapacity*gauss->weight*Jdet;

		TripleMultiply(&L[0],numdof,1,0,
					&D_scalar,1,1,0,
					&L[0],1,numdof,0,
					&Ke->values[0],1);
	}

	/*Clean up and return*/
	delete gauss;
	return Ke;
}
/*}}}*/
/*FUNCTION Tria::CreateKMatrixPrognostic {{{1*/
ElementMatrix* Tria::CreateKMatrixPrognostic(void){

	switch(GetElementType()){
		case P1Enum:
			return CreateKMatrixPrognostic_CG();
		case P1DGEnum:
			return CreateKMatrixPrognostic_DG();
		default:
			_error_("Element type %s not supported yet",EnumToStringx(GetElementType()));
	}

}
/*}}}*/
/*FUNCTION Tria::CreateKMatrixPrognostic_CG {{{1*/
ElementMatrix* Tria::CreateKMatrixPrognostic_CG(void){

	/*Constants*/
	const int    numdof=NDOF1*NUMVERTICES;

	/*Intermediaries */
	int        stabilization;
	int        i,j,ig,dim;
	double     Jdettria,DL_scalar,dt,h;
	double     vel,vx,vy,dvxdx,dvydy;
	double     dvx[2],dvy[2];
	double     v_gauss[2]={0.0};
	double     xyz_list[NUMVERTICES][3];
	double     L[NUMVERTICES];
	double     B[2][NUMVERTICES];
	double     Bprime[2][NUMVERTICES];
	double     K[2][2]                        ={0.0};
	double     KDL[2][2]                      ={0.0};
	double     DL[2][2]                        ={0.0};
	double     DLprime[2][2]                   ={0.0};
	GaussTria *gauss=NULL;

	/*Initialize Element matrix*/
	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);

	/*Retrieve all inputs and parameters*/
	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
	this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
	this->parameters->FindParam(&dim,MeshDimensionEnum);
	this->parameters->FindParam(&stabilization,PrognosticStabilizationEnum);
	Input* vxaverage_input=NULL;
	Input* vyaverage_input=NULL;
	if(dim==2){
		vxaverage_input=inputs->GetInput(VxEnum); _assert_(vxaverage_input);
		vyaverage_input=inputs->GetInput(VyEnum); _assert_(vyaverage_input);
	}
	else{
		vxaverage_input=inputs->GetInput(VxAverageEnum); _assert_(vxaverage_input);
		vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input);
	}
	h=sqrt(2*this->GetArea());

	/* Start  looping on the number of gaussian points: */
	gauss=new GaussTria(2);
	for (ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
		GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);

		vxaverage_input->GetInputValue(&vx,gauss);
		vyaverage_input->GetInputValue(&vy,gauss);
		vxaverage_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
		vyaverage_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);

		DL_scalar=gauss->weight*Jdettria;

		TripleMultiply( &L[0],1,numdof,1,
					&DL_scalar,1,1,0,
					&L[0],1,numdof,0,
					&Ke->values[0],1);

		GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss);
		GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);

		dvxdx=dvx[0];
		dvydy=dvy[1];
		DL_scalar=dt*gauss->weight*Jdettria;

		DL[0][0]=DL_scalar*dvxdx;
		DL[1][1]=DL_scalar*dvydy;
		DLprime[0][0]=DL_scalar*vx;
		DLprime[1][1]=DL_scalar*vy;

		TripleMultiply( &B[0][0],2,numdof,1,
					&DL[0][0],2,2,0,
					&B[0][0],2,numdof,0,
					&Ke->values[0],1);

		TripleMultiply( &B[0][0],2,numdof,1,
					&DLprime[0][0],2,2,0,
					&Bprime[0][0],2,numdof,0,
					&Ke->values[0],1);

		if(stabilization==2){
			/*Streamline upwinding*/
			vel=sqrt(pow(vx,2.)+pow(vy,2.));
			K[0][0]=h/(2*vel)*vx*vx;
			K[1][0]=h/(2*vel)*vy*vx;
			K[0][1]=h/(2*vel)*vx*vy;
			K[1][1]=h/(2*vel)*vy*vy;
		}
		else if(stabilization==1){
			/*MacAyeal*/
			vxaverage_input->GetInputAverage(&vx);
			vyaverage_input->GetInputAverage(&vy);
			K[0][0]=h/2.0*fabs(vx);
			K[0][1]=0.;
			K[1][0]=0.;
			K[1][1]=h/2.0*fabs(vy);
		}
		if(stabilization==1 || stabilization==2){
			KDL[0][0]=DL_scalar*K[0][0];
			KDL[1][0]=DL_scalar*K[1][0];
			KDL[0][1]=DL_scalar*K[0][1];
			KDL[1][1]=DL_scalar*K[1][1];
			TripleMultiply( &Bprime[0][0],2,numdof,1,
						&KDL[0][0],2,2,0,
						&Bprime[0][0],2,numdof,0,
						&Ke->values[0],1);
		}
	}

	/*Clean up and return*/
	delete gauss;
	return Ke;
}
/*}}}*/
/*FUNCTION Tria::CreateKMatrixPrognostic_DG {{{1*/
ElementMatrix* Tria::CreateKMatrixPrognostic_DG(void){

	/*Constants*/
	const int    numdof=NDOF1*NUMVERTICES;

	/*Intermediaries */
	int        i,j,ig,dim;
	double     xyz_list[NUMVERTICES][3];
	double     Jdettria,dt,vx,vy;
	double     L[NUMVERTICES];
	double     B[2][NUMVERTICES];
	double     Bprime[2][NUMVERTICES];
	double     DL[2][2]={0.0};
	double     DLprime[2][2]={0.0};
	double     DL_scalar;
	GaussTria  *gauss=NULL;

	/*Initialize Element matrix*/
	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);

	/*Retrieve all inputs and parameters*/
	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
	this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
	this->parameters->FindParam(&dim,MeshDimensionEnum);
	Input* vxaverage_input=NULL;
	Input* vyaverage_input=NULL;
	if(dim==2){
		vxaverage_input=inputs->GetInput(VxEnum); _assert_(vxaverage_input);
		vyaverage_input=inputs->GetInput(VyEnum); _assert_(vyaverage_input);
	}
	else{
		vxaverage_input=inputs->GetInput(VxAverageEnum); _assert_(vxaverage_input);
		vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input);
	}

	/* Start  looping on the number of gaussian points: */
	gauss=new GaussTria(2);
	for (ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		vxaverage_input->GetInputValue(&vx,gauss);
		vyaverage_input->GetInputValue(&vy,gauss);

		GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
		GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);

		DL_scalar=gauss->weight*Jdettria;

		TripleMultiply( &L[0],1,numdof,1,
					&DL_scalar,1,1,0,
					&L[0],1,numdof,0,
					&Ke->values[0],1);

		/*WARNING: B and Bprime are inverted compared to usual prognostic!!!!*/
		GetBPrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
		GetBprimePrognostic(&B[0][0], &xyz_list[0][0], gauss);

		DL_scalar=-dt*gauss->weight*Jdettria;

		DLprime[0][0]=DL_scalar*vx;
		DLprime[1][1]=DL_scalar*vy;

		TripleMultiply( &B[0][0],2,numdof,1,
					&DLprime[0][0],2,2,0,
					&Bprime[0][0],2,numdof,0,
					&Ke->values[0],1);
	}

	/*Clean up and return*/
	delete gauss;
	return Ke;
}
/*}}}*/
/*FUNCTION Tria::CreateKMatrixSlope {{{1*/
ElementMatrix* Tria::CreateKMatrixSlope(void){

	/*constants: */
	const int    numdof=NDOF1*NUMVERTICES;

	/* Intermediaries */
	int        i,j,ig;
	double     DL_scalar,Jdet;
	double     xyz_list[NUMVERTICES][3];
	double     L[1][3];
	GaussTria *gauss = NULL;

	/*Initialize Element matrix*/
	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);

	GetVerticesCoordinates(&xyz_list[0][0],nodes,NUMVERTICES);

	/* Start looping on the number of gaussian points: */
	gauss=new GaussTria(2);
	for (ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);
		
		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
		DL_scalar=gauss->weight*Jdet;

		GetL(&L[0][0], &xyz_list[0][0], gauss,NDOF1);

		TripleMultiply(&L[0][0],1,3,1,
					&DL_scalar,1,1,0,
					&L[0][0],1,3,0,
					&Ke->values[0],1);
	}

	/*Clean up and return*/
	delete gauss;
	return Ke;
}
/*}}}*/
/*FUNCTION Tria::CreatePVector {{{1*/
void  Tria::CreatePVector(Vec pf){

	/*retrive parameters: */
	ElementVector* pe=NULL;
	int analysis_type;
	parameters->FindParam(&analysis_type,AnalysisTypeEnum);

	/*asserts: {{{*/
	/*if debugging mode, check that all pointers exist*/
	_assert_(this->nodes && this->matice && this->matpar && this->parameters && this->inputs);
	/*}}}*/
	
	/*Skip if water element*/
	if(IsOnWater()) return;

	/*Just branch to the correct load generator, according to the type of analysis we are carrying out: */
	switch(analysis_type){
		#ifdef _HAVE_DIAGNOSTIC_
		case DiagnosticHorizAnalysisEnum:
			pe=CreatePVectorDiagnosticMacAyeal();
			break;
		case DiagnosticHutterAnalysisEnum:
			pe=CreatePVectorDiagnosticHutter();
			break;
		#endif
		case BedSlopeXAnalysisEnum: case SurfaceSlopeXAnalysisEnum: case BedSlopeYAnalysisEnum: case SurfaceSlopeYAnalysisEnum:
			pe=CreatePVectorSlope();
			break;
		case PrognosticAnalysisEnum:
			pe=CreatePVectorPrognostic();
			break;
		#ifdef _HAVE_HYDROLOGY_
		case HydrologyAnalysisEnum:
			pe=CreatePVectorHydrology();
			break;
		#endif
		#ifdef _HAVE_BALANCED_
		case BalancethicknessAnalysisEnum:
			pe=CreatePVectorBalancethickness();
			break;
		#endif
		#ifdef _HAVE_CONTROL_
		case AdjointBalancethicknessAnalysisEnum:
			pe=CreatePVectorAdjointBalancethickness();
			break;
		case AdjointHorizAnalysisEnum:
			pe=CreatePVectorAdjointHoriz();
			break;
		#endif
		default:
			_error_("analysis %i (%s) not supported yet",analysis_type,EnumToStringx(analysis_type));
	}

	/*Add to global Vector*/
	if(pe){
		pe->AddToGlobal(pf);
		delete pe;
	}
}
/*}}}*/
/*FUNCTION Tria::CreatePVectorPrognostic{{{1*/
ElementVector* Tria::CreatePVectorPrognostic(void){

	switch(GetElementType()){
		case P1Enum:
			return CreatePVectorPrognostic_CG();
		case P1DGEnum:
			return CreatePVectorPrognostic_DG();
		default:
			_error_("Element type %s not supported yet",EnumToStringx(GetElementType()));
	}
}
/*}}}*/
/*FUNCTION Tria::CreatePVectorPrognostic_CG {{{1*/
ElementVector* Tria::CreatePVectorPrognostic_CG(void){

	/*Constants*/
	const int    numdof=NDOF1*NUMVERTICES;

	/*Intermediaries */
	int        i,j,ig;
	double     Jdettria,dt;
	double     surface_mass_balance_g,basal_melting_g,basal_melting_correction_g,thickness_g;
	double     xyz_list[NUMVERTICES][3];
	double     L[NUMVERTICES];
	GaussTria* gauss=NULL;

	/*Initialize Element vector*/
	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);

	/*Retrieve all inputs and parameters*/
	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
	this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
	Input* surface_mass_balance_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(surface_mass_balance_input);
	Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum);          _assert_(basal_melting_input);
	Input* basal_melting_correction_input=inputs->GetInput(BasalforcingsMeltingRateCorrectionEnum);
	Input* thickness_input=inputs->GetInput(ThicknessEnum);                             _assert_(thickness_input);

	/*Initialize basal_melting_correction_g to 0, do not forget!:*/
	/* Start  looping on the number of gaussian points: */
	gauss=new GaussTria(2);
	for(ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
		GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);

		surface_mass_balance_input->GetInputValue(&surface_mass_balance_g,gauss);
		basal_melting_input->GetInputValue(&basal_melting_g,gauss);
		thickness_input->GetInputValue(&thickness_g,gauss);
		if(basal_melting_correction_input) basal_melting_correction_input->GetInputValue(&basal_melting_correction_g,gauss);

		for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(thickness_g+dt*(surface_mass_balance_g-basal_melting_g-basal_melting_correction_g))*L[i];
	}

	/*Clean up and return*/
	delete gauss;
	return pe;
}
/*}}}*/
/*FUNCTION Tria::CreatePVectorPrognostic_DG {{{1*/
ElementVector* Tria::CreatePVectorPrognostic_DG(void){

	/*Constants*/
	const int    numdof=NDOF1*NUMVERTICES;

	/*Intermediaries */
	int        i,j,ig;
	double     Jdettria,dt;
	double     surface_mass_balance_g,basal_melting_g,thickness_g;
	double     xyz_list[NUMVERTICES][3];
	double     L[NUMVERTICES];
	GaussTria* gauss=NULL;

	/*Initialize Element vector*/
	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);

	/*Retrieve all inputs and parameters*/
	this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
	Input* surface_mass_balance_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(surface_mass_balance_input);
	Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum);          _assert_(basal_melting_input);
	Input* thickness_input=inputs->GetInput(ThicknessEnum);                             _assert_(thickness_input);

	/* Start  looping on the number of gaussian points: */
	gauss=new GaussTria(2);
	for(ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
		GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);

		surface_mass_balance_input->GetInputValue(&surface_mass_balance_g,gauss);
		basal_melting_input->GetInputValue(&basal_melting_g,gauss);
		thickness_input->GetInputValue(&thickness_g,gauss);

		for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(thickness_g+dt*(surface_mass_balance_g-basal_melting_g))*L[i];
	}

	/*Clean up and return*/
	delete gauss;
	return pe;
}
/*}}}*/
/*FUNCTION Tria::CreatePVectorSlope {{{1*/
ElementVector* Tria::CreatePVectorSlope(void){

	/*Constants*/
	const int    numdof=NDOF1*NUMVERTICES;
	
	/*Intermediaries */
	int        i,j,ig;
	int        analysis_type;
	double     Jdet;
	double     xyz_list[NUMVERTICES][3];
	double     slope[2];
	double     basis[3];
	GaussTria* gauss=NULL;

	/*Initialize Element vector*/
	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);

	/*Retrieve all inputs and parameters*/
	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
	Input* slope_input=NULL;
	if ( (analysis_type==SurfaceSlopeXAnalysisEnum) || (analysis_type==SurfaceSlopeYAnalysisEnum)){
		slope_input=inputs->GetInput(SurfaceEnum); _assert_(slope_input);
	}
	if ( (analysis_type==BedSlopeXAnalysisEnum) || (analysis_type==BedSlopeYAnalysisEnum)){
		slope_input=inputs->GetInput(BedEnum);     _assert_(slope_input);
	}
		
	/* Start  looping on the number of gaussian points: */
	gauss=new GaussTria(2);
	for(ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
		GetNodalFunctions(basis, gauss);

		slope_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss);

		if ( (analysis_type==SurfaceSlopeXAnalysisEnum) || (analysis_type==BedSlopeXAnalysisEnum)){
			for(i=0;i<numdof;i++) pe->values[i]+=Jdet*gauss->weight*slope[0]*basis[i];
		}
		if ( (analysis_type==SurfaceSlopeYAnalysisEnum) || (analysis_type==BedSlopeYAnalysisEnum)){
			for(i=0;i<numdof;i++) pe->values[i]+=Jdet*gauss->weight*slope[1]*basis[i];
		}
	}

	/*Clean up and return*/
	delete gauss;
	return pe;
}
/*}}}*/
/*FUNCTION Tria::CreateJacobianMatrix{{{1*/
void  Tria::CreateJacobianMatrix(Mat Jff){

	/*retrieve parameters: */
	ElementMatrix* Ke=NULL;
	int analysis_type;
	parameters->FindParam(&analysis_type,AnalysisTypeEnum);

	/*Checks in debugging {{{2*/
	_assert_(this->nodes && this->matice && this->matpar && this->parameters && this->inputs);
	/*}}}*/

	/*Skip if water element*/
	if(IsOnWater()) return;

	/*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
	switch(analysis_type){
#ifdef _HAVE_DIAGNOSTIC_
		case DiagnosticHorizAnalysisEnum:
			Ke=CreateJacobianDiagnosticMacayeal();
			break;
#endif
		default:
			_error_("analysis %i (%s) not supported yet",analysis_type,EnumToStringx(analysis_type));
	}

	/*Add to global matrix*/
	if(Ke){
		Ke->AddToGlobal(Jff);
		delete Ke;
	}
}
/*}}}*/
/*FUNCTION Tria::ComputeBasalStress {{{1*/
void  Tria::ComputeBasalStress(Vec eps){
	_error_("Not Implemented yet");
}
/*}}}*/
/*FUNCTION Tria::ComputeStrainRate {{{1*/
void  Tria::ComputeStrainRate(Vec eps){
	_error_("Not Implemented yet");
}
/*}}}*/
/*FUNCTION Tria::ComputeStressTensor {{{1*/
void  Tria::ComputeStressTensor(){

	int         iv;
	double      xyz_list[NUMVERTICES][3];
	double      pressure,viscosity;
	double      epsilon[3]; /* epsilon=[exx,eyy,exy];*/
	double      sigma_xx[NUMVERTICES];
	double		sigma_yy[NUMVERTICES];
	double		sigma_zz[NUMVERTICES]={0,0,0};
	double      sigma_xy[NUMVERTICES];
	double		sigma_xz[NUMVERTICES]={0,0,0};
	double		sigma_yz[NUMVERTICES]={0,0,0};
	GaussTria* gauss=NULL;

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

	/*Retrieve all inputs we will be needing: */
	Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input);
	Input* vx_input=inputs->GetInput(VxEnum);             _assert_(vx_input);
	Input* vy_input=inputs->GetInput(VyEnum);             _assert_(vy_input);

	/* 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->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
		matice->GetViscosity2d(&viscosity,&epsilon[0]);
		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->inputs->AddInput(new TriaP1Input(StressTensorxxEnum,&sigma_xx[0]));
	this->inputs->AddInput(new TriaP1Input(StressTensorxyEnum,&sigma_xy[0]));
	this->inputs->AddInput(new TriaP1Input(StressTensorxzEnum,&sigma_xz[0]));
	this->inputs->AddInput(new TriaP1Input(StressTensoryyEnum,&sigma_yy[0]));
	this->inputs->AddInput(new TriaP1Input(StressTensoryzEnum,&sigma_yz[0]));
	this->inputs->AddInput(new TriaP1Input(StressTensorzzEnum,&sigma_zz[0]));

	/*Clean up and return*/
	delete gauss;
}
/*}}}*/
/*FUNCTION Tria::Configure {{{1*/
void  Tria::Configure(Elements* elementsin, Loads* loadsin, DataSet* nodesin, Materials* materialsin, Parameters* parametersin){
	
	/*go into parameters and get the analysis_counter: */
	int analysis_counter;
	parametersin->FindParam(&analysis_counter,AnalysisCounterEnum);

	/*Get Element type*/
	this->element_type=this->element_type_list[analysis_counter];

	/*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[analysis_counter]) this->hnodes[analysis_counter]->configure(nodesin);
	this->hmatice->configure(materialsin);
	this->hmatpar->configure(materialsin);

	/*Now, go pick up the objects inside the hooks: */
	if(this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp();
	else this->nodes=NULL;
	this->matice=(Matice*)this->hmatice->delivers();
	this->matpar=(Matpar*)this->hmatpar->delivers();

	/*point parameters to real dataset: */
	this->parameters=parametersin;

	/*get inputs configured too: */
	this->inputs->Configure(parameters);

}
/*}}}*/
/*FUNCTION Tria::DeepEcho{{{1*/
void Tria::DeepEcho(void){

	printf("Tria:\n");
	printf("   id: %i\n",id);
	if(nodes){
		nodes[0]->DeepEcho();
		nodes[1]->DeepEcho();
		nodes[2]->DeepEcho();
	}
	else printf("nodes = NULL\n");

	if (matice) matice->DeepEcho();
	else printf("matice = NULL\n");

	if (matpar) matpar->DeepEcho();
	else printf("matpar = NULL\n");

	printf("   parameters\n");
	if (parameters) parameters->DeepEcho();
	else printf("parameters = NULL\n");

	printf("   inputs\n");
	if (inputs) inputs->DeepEcho();
	else printf("inputs=NULL\n");

	if (results) results->DeepEcho();
	else printf("results=NULL\n");

	printf("neighboor sids: \n");
	printf(" %i %i %i\n",horizontalneighborsids[0],horizontalneighborsids[1],horizontalneighborsids[2]);
	
	return;
}
/*}}}*/
/*FUNCTION Tria::DeleteResults {{{1*/
void  Tria::DeleteResults(void){

	/*Delete and reinitialize results*/
	delete this->results;
	this->results=new Results();

}
/*}}}*/
/*FUNCTION Tria::Echo{{{1*/
void Tria::Echo(void){
	printf("Tria:\n");
	printf("   id: %i\n",id);
	if(nodes){
		nodes[0]->Echo();
		nodes[1]->Echo();
		nodes[2]->Echo();
	}
	else printf("nodes = NULL\n");

	if (matice) matice->Echo();
	else printf("matice = NULL\n");

	if (matpar) matpar->Echo();
	else printf("matpar = NULL\n");

	printf("   parameters\n");
	if (parameters) parameters->Echo();
	else printf("parameters = NULL\n");

	printf("   inputs\n");
	if (inputs) inputs->Echo();
	else printf("inputs=NULL\n");

	if (results) results->Echo();
	else printf("results=NULL\n");

	printf("neighboor sids: \n");
	printf(" %i %i %i\n",horizontalneighborsids[0],horizontalneighborsids[1],horizontalneighborsids[2]);
}
/*}}}*/
/*FUNCTION Tria::ObjectEnum{{{1*/
int Tria::ObjectEnum(void){

	return TriaEnum;

}
/*}}}*/
/*FUNCTION Tria::GetArea {{{1*/
double Tria::GetArea(void){

	double area=0;
	double xyz_list[NUMVERTICES][3];
	double x1,y1,x2,y2,x3,y3;

	/*Get xyz list: */
	GetVerticesCoordinates(&xyz_list[0][0], nodes, 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;
}
/*}}}*/
/*FUNCTION Tria::GetDofList {{{1*/
void  Tria::GetDofList(int** pdoflist, int approximation_enum,int setenum){

	int i,j;
	int count=0;
	int numberofdofs=0;
	int* doflist=NULL;

	/*First, figure out size of doflist and create it: */
	for(i=0;i<3;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
	doflist=(int*)xmalloc(numberofdofs*sizeof(int));

	/*Populate: */
	count=0;
	for(i=0;i<3;i++){
		nodes[i]->GetDofList(doflist+count,approximation_enum,setenum);
		count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
	}

	/*Assign output pointers:*/
	*pdoflist=doflist;
}
/*}}}*/
/*FUNCTION Tria::GetDofList1 {{{1*/
void  Tria::GetDofList1(int* doflist){

	int i;
	for(i=0;i<3;i++) doflist[i]=nodes[i]->GetDofList1();

}
/*}}}*/
/*FUNCTION Tria::GetElementType {{{1*/
int Tria::GetElementType(){

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

}
/*}}}*/
/*FUNCTION Tria::GetHorizontalNeighboorSids {{{1*/
int* Tria::GetHorizontalNeighboorSids(){

	/*return TriaRef field*/
	return &this->horizontalneighborsids[0];

}
/*}}}*/
/*FUNCTION Tria::GetNodeIndex {{{1*/
int Tria::GetNodeIndex(Node* node){

	_assert_(nodes);
	for(int i=0;i<NUMVERTICES;i++){
		if(node==nodes[i])
		 return i;
	}
	_error_("Node provided not found among element nodes");
}
/*}}}*/
/*FUNCTION Tria::GetInputListOnVertices(double* pvalue,int enumtype) {{{1*/
void Tria::GetInputListOnVertices(double* pvalue,int enumtype){

	/*Intermediaries*/
	double     value[NUMVERTICES];
	GaussTria *gauss              = NULL;

	/*Recover input*/
	Input* input=inputs->GetInput(enumtype);
	if (!input) _error_("Input %s not found in element",EnumToStringx(enumtype));

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

	/* Start looping on the number of vertices: */
	gauss=new GaussTria();
	for (int iv=0;iv<NUMVERTICES;iv++){
		gauss->GaussVertex(iv);
		input->GetInputValue(&pvalue[iv],gauss);
	}

	/*clean-up*/
	delete gauss;
}
/*}}}*/
/*FUNCTION Tria::GetInputListOnVertices(double* pvalue,int enumtype,double defaultvalue) {{{1*/
void Tria::GetInputListOnVertices(double* pvalue,int enumtype,double defaultvalue){

	double     value[NUMVERTICES];
	GaussTria *gauss = NULL;
	Input     *input = inputs->GetInput(enumtype);

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

	/* Start looping on the number of vertices: */
	if (input){
		gauss=new GaussTria();
		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]=defaultvalue;
	}

	/*clean-up*/
	delete gauss;
}
/*}}}*/
/*FUNCTION Tria::GetInputListOnVertices(double* pvalue,int enumtype,double defaultvalue,int index) TO BE REMOVED{{{1*/
void Tria::GetInputListOnVertices(double* pvalue,int enumtype,double defaultvalue,int index){

	double     value[NUMVERTICES];
	GaussTria *gauss = NULL;
	Input     *input = inputs->GetInput(enumtype);

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

	/* Start looping on the number of vertices: */
	if (input){
		gauss=new GaussTria();
		for (int iv=0;iv<NUMVERTICES;iv++){
			gauss->GaussVertex(iv);
			input->GetInputValue(&pvalue[iv],gauss,index);
		}
	}
	else{
		for (int iv=0;iv<NUMVERTICES;iv++) pvalue[iv]=defaultvalue;
	}

	/*clean-up*/
	delete gauss;
}
/*}}}*/
/*FUNCTION Tria::GetInputValue(double* pvalue,Node* node,int enumtype) {{{1*/
void Tria::GetInputValue(double* pvalue,Node* node,int enumtype){

	Input* input=inputs->GetInput(enumtype);
	if(!input) _error_("No input of type %s found in tria",EnumToStringx(enumtype));

	GaussTria* gauss=new GaussTria();
	gauss->GaussVertex(this->GetNodeIndex(node));

	input->GetInputValue(pvalue,gauss);
	delete gauss;
}
/*}}}*/
/*FUNCTION Tria::GetSidList {{{1*/
void  Tria::GetSidList(int* sidlist){
	for(int i=0;i<NUMVERTICES;i++) sidlist[i]=nodes[i]->GetSidList();
}
/*}}}*/
/*FUNCTION Tria::GetConnectivityList {{{1*/
void  Tria::GetConnectivityList(int* connectivity){
	for(int i=0;i<NUMVERTICES;i++) connectivity[i]=nodes[i]->GetConnectivity();
}
/*}}}*/
/*FUNCTION Tria::GetSolutionFromInputs{{{1*/
void  Tria::GetSolutionFromInputs(Vec solution){

	/*retrive parameters: */
	int analysis_type;
	parameters->FindParam(&analysis_type,AnalysisTypeEnum);

	/*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */
	switch(analysis_type){
	#ifdef _HAVE_DIAGNOSTIC_
	case DiagnosticHorizAnalysisEnum:
		GetSolutionFromInputsDiagnosticHoriz(solution);
		break;
	case DiagnosticHutterAnalysisEnum:
		GetSolutionFromInputsDiagnosticHutter(solution);
		break;
	#endif
	#ifdef _HAVE_HYDROLOGY_
	case HydrologyAnalysisEnum:
		GetSolutionFromInputsHydrology(solution);
		break;
	#endif
	default:
		_error_("analysis: %s not supported yet",EnumToStringx(analysis_type));
	}

}
/*}}}*/
/*FUNCTION Tria::GetStrainRate2d(double* epsilon,double* xyz_list, GaussTria* gauss, Input* vx_input, Input* vy_input){{{1*/
void Tria::GetStrainRate2d(double* epsilon,double* xyz_list, GaussTria* gauss, Input* vx_input, Input* vy_input){
	/*Compute the 2d Strain Rate (3 components):
	 * epsilon=[exx eyy exy] */

	int i;
	double epsilonvx[3];
	double epsilonvy[3];

	/*Check that both inputs have been found*/
	if (!vx_input || !vy_input){
		_error_("Input missing. Here are the input pointers we have for vx: %p, vy: %p\n",vx_input,vy_input);
	}

	/*Get strain rate assuming that epsilon has been allocated*/
	vx_input->GetVxStrainRate2d(epsilonvx,xyz_list,gauss);
	vy_input->GetVyStrainRate2d(epsilonvy,xyz_list,gauss);

	/*Sum all contributions*/
	for(i=0;i<3;i++) epsilon[i]=epsilonvx[i]+epsilonvy[i];
}
/*}}}*/
/*FUNCTION Tria::GetVectorFromInputs{{{1*/
void  Tria::GetVectorFromInputs(Vec vector,int input_enum){

	int doflist1[NUMVERTICES];

	/*Get out if this is not an element input*/
	if(!IsInput(input_enum)) return;

	/*Prepare index list*/
	this->GetDofList1(&doflist1[0]);

	/*Get input (either in element or material)*/
	Input* input=inputs->GetInput(input_enum);
	if(!input) _error_("Input %s not found in element",EnumToStringx(input_enum));

	/*We found the enum.  Use its values to fill into the vector, using the vertices ids: */
	input->GetVectorFromInputs(vector,&doflist1[0]);
}
/*}}}*/
/*FUNCTION Tria::GetVectorFromResults{{{1*/
void  Tria::GetVectorFromResults(Vec vector,int offset,int interp){

	/*Get result*/
	ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(offset);
	if(interp==P1Enum){
		int doflist1[NUMVERTICES];
		int connectivity[NUMVERTICES];
		this->GetSidList(&doflist1[0]);
		this->GetConnectivityList(&connectivity[0]);
		elementresult->GetVectorFromResults(vector,&doflist1[0],&connectivity[0],NUMVERTICES);
	}
	else if(interp==P0Enum){
		elementresult->GetElementVectorFromResults(vector,sid);
	}
	else{
		printf("Interpolation %s not supported\n",EnumToStringx(interp));
	}
}
/*}}}*/
/*FUNCTION Tria::Id {{{1*/
int    Tria::Id(){
	
	return id;

}
/*}}}*/
/*FUNCTION Tria::Sid {{{1*/
int    Tria::Sid(){
	
	return sid;

}
/*}}}*/
/*FUNCTION Tria::InputArtificialNoise{{{1*/
void  Tria::InputArtificialNoise(int enum_type,double min,double max){

	Input* input=NULL;

	/*Make a copy of the original input: */
	input=(Input*)this->inputs->GetInput(enum_type);
	if(!input)_error_(" could not find old input with enum: %s",EnumToStringx(enum_type));

	/*ArtificialNoise: */
	input->ArtificialNoise(min,max);
}
/*}}}*/
/*FUNCTION Tria::InputConvergence{{{1*/
bool Tria::InputConvergence(double* eps, int* enums,int num_enums,int* criterionenums,double* criterionvalues,int num_criterionenums){

	bool    converged=true;
	int     i;
	Input** new_inputs=NULL;
	Input** old_inputs=NULL;

	new_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the new inputs
	old_inputs=(Input**)xmalloc(num_enums/2*sizeof(Input*)); //half the enums are for the old inputs

	for(i=0;i<num_enums/2;i++){
		new_inputs[i]=(Input*)this->inputs->GetInput(enums[2*i+0]);
		old_inputs[i]=(Input*)this->inputs->GetInput(enums[2*i+1]);
		if(!new_inputs[i])_error_("%s%s"," could not find input with enum ",EnumToStringx(enums[2*i+0]));
		if(!old_inputs[i])_error_("%s%s"," could not find input with enum ",EnumToStringx(enums[2*i+0]));
	}

	/*ok, we've got the inputs (new and old), now loop throught the number of criterions and fill the eps array:*/
	for(i=0;i<num_criterionenums;i++){
		IsInputConverged(eps+i,new_inputs,old_inputs,num_enums/2,criterionenums[i]);
		if(eps[i]>criterionvalues[i]) converged=false; 
	}

	/*clean up and return*/
	xfree((void**)&new_inputs);
	xfree((void**)&old_inputs);
	return converged;
}
/*}}}*/
/*FUNCTION Tria::InputDepthAverageAtBase {{{1*/
void  Tria::InputDepthAverageAtBase(int enum_type,int average_enum_type,int object_enum){

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

	/*copy input of enum_type*/
	if (object_enum==MeshElementsEnum)
	 oldinput=(Input*)this->inputs->GetInput(enum_type);
	else if (object_enum==MaterialsEnum)
	 oldinput=(Input*)this->matice->inputs->GetInput(enum_type);
	else
	 _error_("object %s not supported yet",EnumToStringx(object_enum));
	if(!oldinput)_error_("%s%s"," could not find old input with enum: ",EnumToStringx(enum_type));
	newinput=(Input*)oldinput->copy();

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

	/*Add new input to current element*/
	if (object_enum==MeshElementsEnum)
	 this->inputs->AddInput((Input*)newinput);
	else if (object_enum==MaterialsEnum)
	 this->matice->inputs->AddInput((Input*)newinput);
	else
	 _error_("object %s not supported yet",EnumToStringx(object_enum));
}
/*}}}*/
/*FUNCTION Tria::InputDuplicate{{{1*/
void  Tria::InputDuplicate(int original_enum,int new_enum){

	/*Call inputs method*/
	if (IsInput(original_enum)) inputs->DuplicateInput(original_enum,new_enum);

}
/*}}}*/
/*FUNCTION Tria::InputScale{{{1*/
void  Tria::InputScale(int enum_type,double scale_factor){

	Input* input=NULL;

	/*Make a copy of the original input: */
	input=(Input*)this->inputs->GetInput(enum_type);
	if(!input)_error_(" could not find old input with enum: %s",EnumToStringx(enum_type));

	/*Scale: */
	input->Scale(scale_factor);
}
/*}}}*/
/*FUNCTION Tria::InputToResult{{{1*/
void  Tria::InputToResult(int enum_type,int step,double time){

	int    i;
	Input *input = NULL;

	/*Go through all the input objects, and find the one corresponding to enum_type, if it exists: */
	if (enum_type==MaterialsRheologyBbarEnum) input=this->matice->inputs->GetInput(enum_type);
	else input=this->inputs->GetInput(enum_type);
	//if (!input) _error_("Input %s not found in tria->inputs",EnumToStringx(enum_type));
	if(!input)return;

	/*If we don't find it, no big deal, just don't do the transfer. Otherwise, build a new Result 
	 * object out of the input, with the additional step and time information: */
	this->results->AddObject((Object*)input->SpawnResult(step,time));
	
	#ifdef _HAVE_CONTROL_
	if(input->ObjectEnum()==ControlInputEnum){
		if(((ControlInput*)input)->gradient!=NULL) this->results->AddObject((Object*)((ControlInput*)input)->SpawnGradient(step,time));
	}
	#endif
}
/*}}}*/
/*FUNCTION Tria::InputUpdateFromConstant(int value, int name);{{{1*/
void  Tria::InputUpdateFromConstant(int constant, int name){
	/*Check that name is an element input*/
	if (!IsInput(name)) return;

	/*update input*/
	this->inputs->AddInput(new IntInput(name,constant));
}
/*}}}*/
/*FUNCTION Tria::InputUpdateFromConstant(double value, int name);{{{1*/
void  Tria::InputUpdateFromConstant(double constant, int name){
	/*Check that name is an element input*/
	if (!IsInput(name)) return;

	/*update input*/
	this->inputs->AddInput(new DoubleInput(name,constant));
}
/*}}}*/
/*FUNCTION Tria::InputUpdateFromConstant(bool value, int name);{{{1*/
void  Tria::InputUpdateFromConstant(bool constant, int name){
	/*Check that name is an element input*/
	if (!IsInput(name)) return;

	/*update input*/
	this->inputs->AddInput(new BoolInput(name,constant));
}
/*}}}*/
/*FUNCTION Tria::InputUpdateFromIoModel{{{1*/
void Tria::InputUpdateFromIoModel(int index, IoModel* iomodel){ //i is the element index

	/*Intermediaries*/
	int    i,j;
	int    tria_vertex_ids[3];
	double nodeinputs[3];
	double cmmininputs[3];
	double cmmaxinputs[3];
	bool   control_analysis=false;
	int    num_control_type;
	double yts;
	int    num_cm_responses;
   
	/*Get parameters: */
	iomodel->Constant(&yts,ConstantsYtsEnum); 
	iomodel->Constant(&control_analysis,InversionIscontrolEnum);
	if(control_analysis) iomodel->Constant(&num_control_type,InversionNumControlParametersEnum);
	if(control_analysis) iomodel->Constant(&num_cm_responses,InversionNumCostFunctionsEnum);

	/*Recover vertices ids needed to initialize inputs*/
	for(i=0;i<3;i++){ 
		tria_vertex_ids[i]=(int)iomodel->Data(MeshElementsEnum)[3*index+i]; //ids for vertices are in the elements array from Matlab
	}

	/*Control Inputs*/
	#ifdef _HAVE_CONTROL_
	if (control_analysis && iomodel->Data(InversionControlParametersEnum)){
		for(i=0;i<num_control_type;i++){
			switch((int)iomodel->Data(InversionControlParametersEnum)[i]){
				case BalancethicknessThickeningRateEnum:
					if (iomodel->Data(BalancethicknessThickeningRateEnum)){
						for(j=0;j<3;j++)nodeinputs[j]=iomodel->Data(BalancethicknessThickeningRateEnum)[tria_vertex_ids[j]-1]/yts;
						for(j=0;j<3;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
						for(j=0;j<3;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
						this->inputs->AddInput(new ControlInput(BalancethicknessThickeningRateEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
					}
					break;
				case VxEnum:
					if (iomodel->Data(VxEnum)){
						for(j=0;j<3;j++)nodeinputs[j]=iomodel->Data(VxEnum)[tria_vertex_ids[j]-1]/yts;
						for(j=0;j<3;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
						for(j=0;j<3;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
						this->inputs->AddInput(new ControlInput(VxEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
					}
					break;
				case VyEnum:
					if (iomodel->Data(VyEnum)){
						for(j=0;j<3;j++)nodeinputs[j]=iomodel->Data(VyEnum)[tria_vertex_ids[j]-1]/yts;
						for(j=0;j<3;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
						for(j=0;j<3;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i]/yts;
						this->inputs->AddInput(new ControlInput(VyEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
					}
					break;
				case FrictionCoefficientEnum:
					if (iomodel->Data(FrictionCoefficientEnum)){
						for(j=0;j<3;j++)nodeinputs[j]=iomodel->Data(FrictionCoefficientEnum)[tria_vertex_ids[j]-1];
						for(j=0;j<3;j++)cmmininputs[j]=iomodel->Data(InversionMinParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i];
						for(j=0;j<3;j++)cmmaxinputs[j]=iomodel->Data(InversionMaxParametersEnum)[(tria_vertex_ids[j]-1)*num_control_type+i];
						this->inputs->AddInput(new ControlInput(FrictionCoefficientEnum,TriaP1InputEnum,nodeinputs,cmmininputs,cmmaxinputs,i+1));
					}
					break;
				case MaterialsRheologyBbarEnum:
					/*Matice will take care of it*/ break;
				default:
					_error_("Control %s not implemented yet",EnumToStringx((int)iomodel->Data(InversionControlParametersEnum)[i]));
			}
		}
	}
	#endif

	/*DatasetInputs*/
	if (control_analysis && iomodel->Data(InversionCostFunctionsCoefficientsEnum)){

		/*Create inputs and add to DataSetInput*/
		DatasetInput* datasetinput=new DatasetInput(InversionCostFunctionsCoefficientsEnum);
		for(i=0;i<num_cm_responses;i++){
			for(j=0;j<3;j++)nodeinputs[j]=iomodel->Data(InversionCostFunctionsCoefficientsEnum)[(tria_vertex_ids[j]-1)*num_cm_responses+i];
			datasetinput->inputs->AddObject(new TriaP1Input(InversionCostFunctionsCoefficientsEnum,nodeinputs));
		}

		/*Add datasetinput to element inputs*/
		this->inputs->AddInput(datasetinput);
	}
}
/*}}}*/
/*FUNCTION Tria::InputUpdateFromSolution {{{1*/
void  Tria::InputUpdateFromSolution(double* solution){

	/*retrive parameters: */
	int analysis_type;
	parameters->FindParam(&analysis_type,AnalysisTypeEnum);

	/*Just branch to the correct InputUpdateFromSolution generator, according to the type of analysis we are carrying out: */
	switch(analysis_type){
		#ifdef _HAVE_DIAGNOSTIC_
		case DiagnosticHorizAnalysisEnum:
			InputUpdateFromSolutionDiagnosticHoriz( solution);
			break;
		case DiagnosticHutterAnalysisEnum:
			InputUpdateFromSolutionDiagnosticHoriz( solution);
			break;
		#endif
		#ifdef _HAVE_CONTROL_
		case AdjointHorizAnalysisEnum:
			InputUpdateFromSolutionAdjointHoriz( solution);
			break;
		case AdjointBalancethicknessAnalysisEnum:
			InputUpdateFromSolutionAdjointBalancethickness( solution);
			break;
		#endif
		#ifdef _HAVE_HYDROLOGY_ 
		case HydrologyAnalysisEnum:
			InputUpdateFromSolutionHydrology(solution);
			break ;
		#endif
	 	#ifdef _HAVE_BALANCED_
		case BalancethicknessAnalysisEnum:
			InputUpdateFromSolutionOneDof(solution,ThicknessEnum);
			break;
		#endif
		case BedSlopeXAnalysisEnum:
			InputUpdateFromSolutionOneDof(solution,BedSlopeXEnum);
			break;
		case BedSlopeYAnalysisEnum:
			InputUpdateFromSolutionOneDof(solution,BedSlopeYEnum);
			break;
		case SurfaceSlopeXAnalysisEnum:
			InputUpdateFromSolutionOneDof(solution,SurfaceSlopeXEnum);
			break;
		case SurfaceSlopeYAnalysisEnum:
			InputUpdateFromSolutionOneDof(solution,SurfaceSlopeYEnum);
			break;
		case PrognosticAnalysisEnum:
			InputUpdateFromSolutionPrognostic(solution);
			break;
		default:
			_error_("analysis %i (%s) not supported yet",analysis_type,EnumToStringx(analysis_type));
	}
}
/*}}}*/
/*FUNCTION Tria::InputUpdateFromSolutionOneDof{{{1*/
void  Tria::InputUpdateFromSolutionOneDof(double* solution,int enum_type){

	const int numdof          = NDOF1*NUMVERTICES;

	int*      doflist=NULL;
	double    values[numdof];

	/*Get dof list: */
	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);

	/*Use the dof list to index into the solution vector: */
	for(int i=0;i<numdof;i++){
		values[i]=solution[doflist[i]];
		if(isnan(values[i])) _error_("NaN found in solution vector");
	}

	/*Add input to the element: */
	this->inputs->AddInput(new TriaP1Input(enum_type,values));

	/*Free ressources:*/
	xfree((void**)&doflist);
}
/*}}}*/
/*FUNCTION Tria::InputUpdateFromSolutionPrognostic{{{1*/
void  Tria::InputUpdateFromSolutionPrognostic(double* solution){

	/*Intermediaries*/
	const int numdof = NDOF1*NUMVERTICES;

	int       i,hydroadjustment;
	int*      doflist=NULL;
	double    rho_ice,rho_water,minthickness;
	double    newthickness[numdof];
	double    newbed[numdof];
	double    newsurface[numdof];
	double    oldbed[NUMVERTICES];
	double    oldsurface[NUMVERTICES];
	double    oldthickness[NUMVERTICES];

	/*Get dof list: */
	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);

	/*Use the dof list to index into the solution vector: */
	this->parameters->FindParam(&minthickness,PrognosticMinThicknessEnum);
	for(i=0;i<numdof;i++){
		newthickness[i]=solution[doflist[i]];
		if(isnan(newthickness[i])) _error_("NaN found in solution vector");
		/*Constrain thickness to be at least 1m*/
		if(newthickness[i]<minthickness) newthickness[i]=minthickness;
	}

	/*Get previous bed, thickness and surface*/
	GetInputListOnVertices(&oldbed[0],BedEnum);
	GetInputListOnVertices(&oldsurface[0],SurfaceEnum);
	GetInputListOnVertices(&oldthickness[0],ThicknessEnum);

	/*Fing PrognosticHydrostaticAdjustment to figure out how to update the geometry:*/
	this->parameters->FindParam(&hydroadjustment,PrognosticHydrostaticAdjustmentEnum);
	rho_ice=matpar->GetRhoIce();
	rho_water=matpar->GetRhoWater();

	for(i=0;i<numdof;i++) {
		/*If shelf: hydrostatic equilibrium*/
		if (this->nodes[i]->IsGrounded()){
			newsurface[i]=oldbed[i]+newthickness[i]; //surface = oldbed + newthickness
			newbed[i]=oldbed[i];               //same bed: do nothing
		}
		else{ //this is an ice shelf

			if(hydroadjustment==AbsoluteEnum){
				newsurface[i]=newthickness[i]*(1-rho_ice/rho_water);
				newbed[i]=newthickness[i]*(-rho_ice/rho_water);
			}
			else if(hydroadjustment==IncrementalEnum){
				newsurface[i]=oldsurface[i]+(1.0-rho_ice/rho_water)*(newthickness[i]-oldthickness[i]); //surface = oldsurface + (1-di) * dH 
				newbed[i]=oldbed[i]-rho_ice/rho_water*(newthickness[i]-oldthickness[i]); //bed = oldbed + di * dH
			}
			else _error_("Hydrostatic adjustment %i (%s) not supported yet",hydroadjustment,EnumToStringx(hydroadjustment));
		}
	}

	/*Add input to the element: */
	this->inputs->AddInput(new TriaP1Input(ThicknessEnum,newthickness));
	this->inputs->AddInput(new TriaP1Input(SurfaceEnum,newsurface));
	this->inputs->AddInput(new TriaP1Input(BedEnum,newbed));

	/*Free ressources:*/
	xfree((void**)&doflist);
}
/*}}}*/
/*FUNCTION Tria::InputUpdateFromVector(double* vector, int name, int type);{{{1*/
void  Tria::InputUpdateFromVector(double* vector, int name, int type){

	/*Check that name is an element input*/
	if (!IsInput(name)) return;

	switch(type){

		case VertexEnum:

			/*New TriaP1Input*/
			double values[3];

			/*Get values on the 3 vertices*/
			for (int i=0;i<3;i++){
				values[i]=vector[this->nodes[i]->GetVertexDof()];
			}

			/*update input*/
			if (name==MaterialsRheologyBbarEnum || name==MaterialsRheologyBEnum){
				matice->inputs->AddInput(new TriaP1Input(name,values));
			}
			else{
				this->inputs->AddInput(new TriaP1Input(name,values));
			}
			return;

		default:
			_error_("type %i (%s) not implemented yet",type,EnumToStringx(type));
	}
}
/*}}}*/
/*FUNCTION Tria::InputUpdateFromVector(int* vector, int name, int type);{{{1*/
void  Tria::InputUpdateFromVector(int* vector, int name, int type){
	_error_(" not supported yet!");
}
/*}}}*/
/*FUNCTION Tria::InputUpdateFromVector(bool* vector, int name, int type);{{{1*/
void  Tria::InputUpdateFromVector(bool* vector, int name, int type){
	_error_(" not supported yet!");
}
/*}}}*/
/*FUNCTION Tria::InputCreate(double scalar,int enum,int code);{{{1*/
void Tria::InputCreate(double scalar,int name,int code){

	/*Check that name is an element input*/
	if (!IsInput(name)) return;
	
	if ((code==5) || (code==1)){ //boolean
		this->inputs->AddInput(new BoolInput(name,(bool)scalar));
	}
	else if ((code==6) || (code==2)){ //integer
		this->inputs->AddInput(new IntInput(name,(int)scalar));
	}
	else if ((code==7) || (code==3)){ //double
		this->inputs->AddInput(new DoubleInput(name,(int)scalar));
	}
	else _error_("%s%i"," could not recognize nature of vector from code ",code);

}
/*}}}*/
/*FUNCTION Tria::InputCreate(double* vector,int index,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){{{1*/
void Tria::InputCreate(double* vector, int index,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){//index into elements

	/*Intermediaries*/
	int    i,j,t;
	int    tria_vertex_ids[3];
	int    row;
	double nodeinputs[3];
	double time;
	TransientInput* transientinput=NULL;
	int    numberofvertices;
	int    numberofelements;
	double yts;


	/*Fetch parameters: */
	iomodel->Constant(&numberofvertices,MeshNumberofverticesEnum);
	iomodel->Constant(&numberofelements,MeshNumberofelementsEnum);
	iomodel->Constant(&yts,ConstantsYtsEnum);

	/*Branch on type of vector: nodal or elementary: */
	if(vector_type==1){ //nodal vector

		/*Recover vertices ids needed to initialize inputs*/
		for(i=0;i<3;i++){ 
			tria_vertex_ids[i]=(int)iomodel->Data(MeshElementsEnum)[3*index+i]; //ids for vertices are in the elements array from Matlab
		}

		/*Are we in transient or static? */
		if(M==numberofvertices){

			/*create input values: */
			for(i=0;i<3;i++)nodeinputs[i]=(double)vector[tria_vertex_ids[i]-1];

			/*process units: */
			UnitConversion(&nodeinputs[0], 3 ,ExtToIuEnum, vector_enum);

			/*create static input: */
			this->inputs->AddInput(new TriaP1Input(vector_enum,nodeinputs));
		}
		else if(M==numberofvertices+1){
			/*create transient input: */
			for(t=0;t<N;t++){ //N is the number of times

				/*create input values: */
				for(i=0;i<3;i++){
					row=tria_vertex_ids[i]-1;
					nodeinputs[i]=(double)vector[N*row+t];
				}

				/*process units: */
				UnitConversion(&nodeinputs[0], 3 ,ExtToIuEnum, vector_enum);

				/*time? :*/
				time=(double)vector[(M-1)*N+t]*yts;

				if(t==0) transientinput=new TransientInput(vector_enum);
				transientinput->AddTimeInput(new TriaP1Input(vector_enum,nodeinputs),time);
			}
			this->inputs->AddInput(transientinput);
		}
		else _error_("nodal vector is either numberofnodes (%i), or numberofnodes+1 long. Field provided is %i long",numberofvertices,M);
	}
	else if(vector_type==2){ //element vector
		/*Are we in transient or static? */
		if(M==numberofelements){

			/*static mode: create an input out of the element value: */

			if (code==5){ //boolean
				this->inputs->AddInput(new BoolInput(vector_enum,(bool)vector[index]));
			}
			else if (code==6){ //integer
				this->inputs->AddInput(new IntInput(vector_enum,(int)vector[index]));
			}
			else if (code==7){ //double
				this->inputs->AddInput(new DoubleInput(vector_enum,(double)vector[index]));
			}
			else _error_("%s%i"," could not recognize nature of vector from code ",code);
		}
		else {
			_error_("transient elementary inputs not supported yet!");
		}
	}
	else{
		_error_("Cannot add input for vector type %i (not supported)",vector_type);
	}

}
/*}}}*/
/*FUNCTION Tria::IsInput{{{1*/
bool Tria::IsInput(int name){
	if (
				name==ThicknessEnum ||
				name==SurfaceEnum ||
				name==BedEnum ||
				name==SurfaceSlopeXEnum ||
				name==SurfaceSlopeYEnum ||
				name==BasalforcingsMeltingRateEnum ||
				name==WatercolumnEnum || 
				name==SurfaceforcingsMassBalanceEnum ||
				name==SurfaceAreaEnum||
				name==VxEnum ||
				name==VyEnum ||
				name==InversionVxObsEnum ||
				name==InversionVyObsEnum ||
				name==FrictionCoefficientEnum ||
				name==MaterialsRheologyBbarEnum ||
				name==GradientEnum ||
				name==OldGradientEnum ||
				name==QmuVxEnum ||
				name==QmuVyEnum ||
				name==QmuPressureEnum ||
				name==QmuBedEnum ||
				name==QmuThicknessEnum ||
				name==QmuSurfaceEnum ||
				name==QmuTemperatureEnum ||
				name==QmuMeltingEnum
		){
		return true;
	}
	else return false;
}
/*}}}*/
/*FUNCTION Tria::IsOnBed {{{1*/
bool Tria::IsOnBed(){
	
	bool onbed;
	inputs->GetInputValue(&onbed,MeshElementonbedEnum);
	return onbed;
}
/*}}}*/
/*FUNCTION Tria::IsFloating {{{1*/
bool   Tria::IsFloating(){

	bool shelf;
	inputs->GetInputValue(&shelf,MaskElementonfloatingiceEnum);
	return shelf;
}
/*}}}*/
/*FUNCTION Tria::IsNodeOnShelf {{{1*/
bool   Tria::IsNodeOnShelf(){

	int  i;
	bool shelf=false;

	for(i=0;i<3;i++){
		if (nodes[i]->IsFloating()){
			shelf=true;
			break;
		}
	}
	return shelf;
}
/*}}}*/
/*FUNCTION Tria::IsNodeOnShelfFromFlags {{{1*/
bool   Tria::IsNodeOnShelfFromFlags(double* flags){

	int  i;
	bool shelf=false;

	for(i=0;i<NUMVERTICES;i++){
		if (flags[nodes[i]->Sid()]){
			shelf=true;
			break;
		}
	}
	return shelf;
}
/*}}}*/
/*FUNCTION Tria::IsOnWater {{{1*/
bool   Tria::IsOnWater(){

	bool water;
	inputs->GetInputValue(&water,MaskElementonwaterEnum);
	return water;
}
/*}}}*/
/*FUNCTION Tria::ListResultsInfo{{{*/
void Tria::ListResultsInfo(int** in_resultsenums,int** in_resultssizes,double** in_resultstimes,int** in_resultssteps,int* in_num_results){

	/*Intermediaries*/
	int     i;
	int     numberofresults = 0;
	int     *resultsenums   = NULL;
	int     *resultssizes   = NULL;
	double  *resultstimes   = NULL;
	int     *resultssteps   = NULL;

	/*Checks*/
	_assert_(in_num_results);

	/*Count number of results*/
	for(i=0;i<this->results->Size();i++){
		ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
		numberofresults++;
	}

	if(numberofresults){

		/*Allocate output*/
		resultsenums=(int*)xmalloc(numberofresults*sizeof(int));
		resultssizes=(int*)xmalloc(numberofresults*sizeof(int));
		resultstimes=(double*)xmalloc(numberofresults*sizeof(double));
		resultssteps=(int*)xmalloc(numberofresults*sizeof(int));

		/*populate enums*/
		for(i=0;i<this->results->Size();i++){
			ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
			resultsenums[i]=elementresult->InstanceEnum();
			resultstimes[i]=elementresult->GetTime();
			resultssteps[i]=elementresult->GetStep();
			if(elementresult->ObjectEnum()==TriaP1ElementResultEnum){
				resultssizes[i]=P1Enum;
			}
			else{
				resultssizes[i]=P0Enum;
			}
		}
	}

	/*Assign output pointers:*/
	*in_num_results=numberofresults;
	*in_resultsenums=resultsenums;
	*in_resultssizes=resultssizes;
	*in_resultstimes=resultstimes;
	*in_resultssteps=resultssteps;

}/*}}}*/
/*FUNCTION Tria::MigrateGroundingLine{{{1*/
void  Tria::MigrateGroundingLine(double* old_floating_ice,double* sheet_ungrounding){

	int     i,migration_style,unground;
	bool    elementonshelf = false;
	double  bed_hydro,yts,gl_melting_rate;
	double  rho_water,rho_ice,density;
	double  melting[NUMVERTICES];
	double  h[NUMVERTICES],s[NUMVERTICES],b[NUMVERTICES],ba[NUMVERTICES];

	/*Recover info at the vertices: */
	parameters->FindParam(&migration_style,GroundinglineMigrationEnum);
	parameters->FindParam(&yts,ConstantsYtsEnum);
	GetInputListOnVertices(&h[0],ThicknessEnum);
	GetInputListOnVertices(&s[0],SurfaceEnum);
	GetInputListOnVertices(&b[0],BedEnum);
	GetInputListOnVertices(&ba[0],BathymetryEnum);
	rho_water=matpar->GetRhoWater();
	rho_ice=matpar->GetRhoIce();
	density=rho_ice/rho_water;
	
	/*go through vertices, and update inputs, considering them to be TriaVertex type: */
	for(i=0;i<NUMVERTICES;i++){
		/*Ice shelf: if bed below bathymetry, impose it at the bathymetry and update surface, elso do nothing */
		if(old_floating_ice[nodes[i]->Sid()]){
			if(b[i]<=ba[i]){ 
				b[i]=ba[i];
				s[i]=b[i]+h[i];
				nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,false));
				nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,true));
			}
		}
		/*Ice sheet: if hydrostatic bed above bathymetry, ice sheet starts to unground, elso do nothing */
		/*Change only if AgressiveMigration or if the ice sheet is in contact with the ocean*/
		else{
			bed_hydro=-density*h[i];
			if (bed_hydro>ba[i]){
				/*Unground only if the element is connected to the ice shelf*/
				if(migration_style==AgressiveMigrationEnum){
					s[i]=(1-density)*h[i];
					b[i]=-density*h[i];
					nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,true));
					nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false));
				}
				else if(migration_style==SoftMigrationEnum && sheet_ungrounding[nodes[i]->Sid()]){
					s[i]=(1-density)*h[i];
					b[i]=-density*h[i];
					nodes[i]->inputs->AddInput(new BoolInput(MaskVertexonfloatingiceEnum,true));
					nodes[i]->inputs->AddInput(new BoolInput(MaskVertexongroundediceEnum,false));
				}
			}
		}
	}

	/*If at least one vertex is now floating, the element is now floating*/
	for(i=0;i<NUMVERTICES;i++){
		if(nodes[i]->IsFloating()){
			elementonshelf=true;
			break;
		}
	}
	
   /*Add basal melting rate if element just ungrounded*/
	if(!this->IsFloating() && elementonshelf==true){
		for(i=0;i<NUMVERTICES;i++)melting[i]=gl_melting_rate/yts;
		this->inputs->AddInput(new TriaP1Input(BasalforcingsMeltingRateEnum,&melting[0]));
	} 

	/*Update inputs*/
   this->inputs->AddInput(new BoolInput(MaskElementonfloatingiceEnum,elementonshelf));

	/*Update inputs*/    
	this->inputs->AddInput(new TriaP1Input(SurfaceEnum,&s[0]));
	this->inputs->AddInput(new TriaP1Input(BedEnum,&b[0]));
}
/*}}}*/
/*FUNCTION Tria::MyRank {{{1*/
int    Tria::MyRank(void){ 
	extern int my_rank;
	return my_rank; 
}
/*}}}*/
/*FUNCTION Tria::PatchFill{{{1*/
void  Tria::PatchFill(int* prow, Patch* patch){

	int i,row;
	int vertices_ids[3];

	/*recover pointer: */
	row=*prow;
		
	for(i=0;i<3;i++) vertices_ids[i]=nodes[i]->GetVertexId(); //vertices id start at column 3 of the patch.

	for(i=0;i<this->results->Size();i++){
		ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);

		/*For this result,fill the information in the Patch object (element id + vertices ids), and then hand 
		 *it to the result object, to fill the rest: */
		patch->fillelementinfo(row,this->sid+1,vertices_ids,3);
		elementresult->PatchFill(row,patch);

		/*increment rower: */
		row++;
	}

	/*Assign output pointers:*/
	*prow=row;
}
/*}}}*/
/*FUNCTION Tria::PatchSize{{{1*/
void  Tria::PatchSize(int* pnumrows, int* pnumvertices,int* pnumnodes){

	int     i;
	int     numrows       = 0;
	int     numnodes      = 0;
	int     temp_numnodes = 0;

	/*Go through all the results objects, and update the counters: */
	for (i=0;i<this->results->Size();i++){
		ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
		/*first, we have one more result: */
		numrows++;
		/*now, how many vertices and how many nodal values for this result? :*/
		temp_numnodes=elementresult->NumberOfNodalValues(); //ask result object.
		if(temp_numnodes>numnodes)numnodes=temp_numnodes;
	}

	/*Assign output pointers:*/
	*pnumrows=numrows;
	*pnumvertices=NUMVERTICES;
	*pnumnodes=numnodes;
}
/*}}}*/
/*FUNCTION Tria::PotentialSheetUngrounding{{{1*/
void  Tria::PotentialSheetUngrounding(Vec potential_sheet_ungrounding){

	int     i;
	double  h[NUMVERTICES],ba[NUMVERTICES];
	double  bed_hydro;
	double  rho_water,rho_ice,density;
	bool    elementonshelf = false;

	/*material parameters: */
	rho_water=matpar->GetRhoWater();
	rho_ice=matpar->GetRhoIce();
	density=rho_ice/rho_water;
	GetInputListOnVertices(&h[0],ThicknessEnum);
	GetInputListOnVertices(&ba[0],BathymetryEnum);

	/*go through vertices, and figure out which ones are on the ice sheet, and want to unground: */
	for(i=0;i<NUMVERTICES;i++){
		/*Find if grounded vertices want to start floating*/
		if (!nodes[i]->IsFloating()){
			bed_hydro=-density*h[i];
			if (bed_hydro>ba[i]){
				/*Vertex that could potentially unground, flag it*/
				VecSetValue(potential_sheet_ungrounding,nodes[i]->Sid(),1,INSERT_VALUES);
			}
		}
	}
}
/*}}}*/
/*FUNCTION Tria::PositiveDegreeDay{{{1*/
void  Tria::PositiveDegreeDay(){

	_error_("Not implemented yet");
}
/*}}}*/
/*FUNCTION Tria::ProcessResultsUnits{{{1*/
void  Tria::ProcessResultsUnits(void){

	int i;

	for(i=0;i<this->results->Size();i++){
		ElementResult* elementresult=(ElementResult*)this->results->GetObjectByOffset(i);
		elementresult->ProcessUnits(this->parameters);
	}
}
/*}}}*/
/*FUNCTION Tria::RequestedOutput{{{1*/
void Tria::RequestedOutput(int output_enum,int step,double time){

	if(IsInput(output_enum)){
		/*just transfer this input to results, and we are done: */
		InputToResult(output_enum,step,time);
	}
	else{
		/*this input does not exist, compute it, and then transfer to results: */
		switch(output_enum){
			case StressTensorEnum: 
				this->ComputeStressTensor();
				InputToResult(StressTensorxxEnum,step,time);
				InputToResult(StressTensorxyEnum,step,time);
				InputToResult(StressTensorxzEnum,step,time);
				InputToResult(StressTensoryyEnum,step,time);
				InputToResult(StressTensoryzEnum,step,time);
				InputToResult(StressTensorzzEnum,step,time);
				break;

			default:
				/*do nothing, no need to derail the computation because one of the outputs requested cannot be found: */
				break;
		}
	}

}
/*}}}*/
/*FUNCTION Tria::SetClone {{{1*/
void  Tria::SetClone(int* minranks){

	_error_("not implemented yet");
}
/*}}}1*/
/*FUNCTION Tria::SmearFunction {{{1*/
void  Tria::SmearFunction(Vec smearedvector,double (*WeightFunction)(double distance,double radius),double radius){
	_error_("not implemented yet");

}
/*}}}1*/
/*FUNCTION Tria::SetCurrentConfiguration {{{1*/
void  Tria::SetCurrentConfiguration(Elements* elementsin, Loads* loadsin, DataSet* nodesin, Materials* materialsin, Parameters* parametersin){
	
	/*go into parameters and get the analysis_counter: */
	int analysis_counter;
	parametersin->FindParam(&analysis_counter,AnalysisCounterEnum);

	/*Get Element type*/
	this->element_type=this->element_type_list[analysis_counter];

	/*Pick up nodes*/
	if(this->hnodes[analysis_counter]) this->nodes=(Node**)this->hnodes[analysis_counter]->deliverp();
	else this->nodes=NULL;

}
/*}}}*/
/*FUNCTION Tria::SurfaceArea {{{1*/
double Tria::SurfaceArea(void){

	int    i;
	double S;
	double normal[3];
	double v13[3],v23[3];
	double xyz_list[NUMVERTICES][3];

	/*If on water, return 0: */
	if(IsOnWater())return 0;

	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);

	for (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(pow(normal[0],(double)2)+pow(normal[1],(double)2)+pow(normal[2],(double)2));

	/*Return: */
	return S;
}
/*}}}*/
/*FUNCTION Tria::SurfaceNormal{{{1*/
void Tria::SurfaceNormal(double* surface_normal, double xyz_list[3][3]){

	int i;
	double v13[3],v23[3];
	double normal[3];
	double normal_norm;

	for (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];

	normal_norm=sqrt( pow(normal[0],(double)2)+pow(normal[1],(double)2)+pow(normal[2],(double)2) );

	*(surface_normal)=normal[0]/normal_norm;
	*(surface_normal+1)=normal[1]/normal_norm;
	*(surface_normal+2)=normal[2]/normal_norm;
}
/*}}}*/
/*FUNCTION Tria::TimeAdapt{{{1*/
double  Tria::TimeAdapt(void){

	/*intermediary: */
	int    i;
	double C,dt;
	double dx,dy;
	double maxx,minx;
	double maxy,miny;
	double maxabsvx,maxabsvy;
	double xyz_list[NUMVERTICES][3];

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

	/*Get for Vx and Vy, the max of abs value: */
	#ifdef _HAVE_RESPONSES_
	this->MaxAbsVx(&maxabsvx,false);
	this->MaxAbsVy(&maxabsvy,false);
	#else
		_error_("ISSM was not compiled with responses compiled in, exiting!");
	#endif

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

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

	for(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];
	}
	dx=maxx-minx;
	dy=maxy-miny;

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

	return dt;
}
/*}}}*/
/*FUNCTION Tria::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type){{{1*/
void Tria::Update(int index, IoModel* iomodel,int analysis_counter,int analysis_type){ //i is the element index

	/*Intermediaries*/
	int    i,j;
	int    tria_node_ids[3];
	int    tria_vertex_ids[3];
	int    tria_type;
	double nodeinputs[3];
	double yts;
	int    progstabilization,balancestabilization;
	bool   dakota_analysis;

	/*Checks if debuging*/
	/*{{{2*/
	_assert_(iomodel->Data(MeshElementsEnum));
	/*}}}*/

	/*Fetch parameters: */
	iomodel->Constant(&yts,ConstantsYtsEnum);
	iomodel->Constant(&progstabilization,PrognosticStabilizationEnum);
	iomodel->Constant(&balancestabilization,BalancethicknessStabilizationEnum);
	iomodel->Constant(&dakota_analysis,QmuIsdakotaEnum);

	/*Recover element type*/
	if ((analysis_type==PrognosticAnalysisEnum && progstabilization==3) || (analysis_type==BalancethicknessAnalysisEnum && balancestabilization==3)){
		/*P1 Discontinuous Galerkin*/
		tria_type=P1DGEnum;
	}
	else{
		/*P1 Continuous Galerkin*/
		tria_type=P1Enum;
	}
	this->SetElementType(tria_type,analysis_counter);

	/*Recover vertices ids needed to initialize inputs*/
	for(i=0;i<3;i++){ 
		tria_vertex_ids[i]=(int)iomodel->Data(MeshElementsEnum)[3*index+i]; //ids for vertices are in the elements array from Matlab
	}

	/*Recover nodes ids needed to initialize the node hook.*/
	if (tria_type==P1DGEnum){
		/*Discontinuous Galerkin*/
		tria_node_ids[0]=iomodel->nodecounter+3*index+1;
		tria_node_ids[1]=iomodel->nodecounter+3*index+2;
		tria_node_ids[2]=iomodel->nodecounter+3*index+3;
	}
	else{
		/*Continuous Galerkin*/
		for(i=0;i<3;i++){ 
			tria_node_ids[i]=iomodel->nodecounter+(int)*(iomodel->Data(MeshElementsEnum)+3*index+i); //ids for vertices are in the elements array from Matlab
		}
	}

	/*hooks: */
	this->SetHookNodes(tria_node_ids,analysis_counter); this->nodes=NULL; //set hook to nodes, for this analysis type

	/*Fill with IoModel*/
	this->InputUpdateFromIoModel(index,iomodel);

	/*Defaults if not provided in iomodel*/
	switch(analysis_type){

		case DiagnosticHorizAnalysisEnum:

			/*default vx,vy and vz: either observation or 0 */
			if(!iomodel->Data(VxEnum)){
				for(i=0;i<3;i++)nodeinputs[i]=0;
				this->inputs->AddInput(new TriaP1Input(VxEnum,nodeinputs));
				if(dakota_analysis) this->inputs->AddInput(new TriaP1Input(QmuVxEnum,nodeinputs));
			}
			if(!iomodel->Data(VyEnum)){
				for(i=0;i<3;i++)nodeinputs[i]=0;
				this->inputs->AddInput(new TriaP1Input(VyEnum,nodeinputs));
				if(dakota_analysis) this->inputs->AddInput(new TriaP1Input(QmuVyEnum,nodeinputs));
			}
			if(!iomodel->Data(VzEnum)){
				for(i=0;i<3;i++)nodeinputs[i]=0;
				this->inputs->AddInput(new TriaP1Input(VzEnum,nodeinputs));
				if(dakota_analysis) this->inputs->AddInput(new TriaP1Input(QmuVzEnum,nodeinputs));
			}
			if(!iomodel->Data(PressureEnum)){
				for(i=0;i<3;i++)nodeinputs[i]=0;
				if(dakota_analysis){
					this->inputs->AddInput(new TriaP1Input(PressureEnum,nodeinputs));
					this->inputs->AddInput(new TriaP1Input(QmuPressureEnum,nodeinputs));
				}
			}
			break;

		default:
			/*No update for other solution types*/
			break;

	}

	//this->parameters: we still can't point to it, it may not even exist. Configure will handle this.
	this->parameters=NULL;
}
/*}}}*/
/*FUNCTION Tria::UpdatePotentialSheetUngrounding{{{1*/
int Tria::UpdatePotentialSheetUngrounding(double* vertices_potentially_ungrounding,Vec vec_nodes_on_iceshelf,double* nodes_on_iceshelf){

	int i;
	int nflipped=0;

	/*Go through nodes, and whoever is on the potential_sheet_ungrounding, ends up in nodes_on_iceshelf: */
	for(i=0;i<3;i++){
		if (vertices_potentially_ungrounding[nodes[i]->Sid()]){
			VecSetValue(vec_nodes_on_iceshelf,nodes[i]->Sid(),1,INSERT_VALUES);
		
			/*If node was not on ice shelf, we flipped*/
			if(nodes_on_iceshelf[nodes[i]->Sid()]==0){
				nflipped++;
			}
		}
	}
	return nflipped;
}
/*}}}*/

#ifdef _HAVE_RESPONSES_
/*FUNCTION Tria::IceVolume {{{1*/
double Tria::IceVolume(void){

	/*The volume of a troncated prism is base * 1/3 sum(length of edges)*/
	double base,surface,bed;
	double xyz_list[NUMVERTICES][3];

	if(IsOnWater())return 0;

	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);

	/*First calculate the area of the base (cross section triangle)
	 * http://en.wikipedia.org/wiki/Triangle
	 * base = 1/2 abs((xA-xC)(yB-yA)-(xA-xB)(yC-yA))*/
	base = 1./2. * fabs((xyz_list[0][0]-xyz_list[2][0])*(xyz_list[1][1]-xyz_list[0][1]) - (xyz_list[0][0]-xyz_list[1][0])*(xyz_list[2][1]-xyz_list[0][1]));

	/*Now get the average height*/
	Input* surface_input = inputs->GetInput(SurfaceEnum); _assert_(surface_input);
	Input* bed_input     = inputs->GetInput(BedEnum);     _assert_(bed_input);
	surface_input->GetInputAverage(&surface);
	bed_input->GetInputAverage(&bed);

	/*Return: */
	return base*(surface-bed);
}
/*}}}*/
/*FUNCTION Tria::MassFlux {{{1*/
double Tria::MassFlux( double* segment,bool process_units){

	const int    numdofs=2;

	int        i,dim;
	double     mass_flux=0;
	double     xyz_list[NUMVERTICES][3];
	double     normal[2];
	double     length,rho_ice;
	double     x1,y1,x2,y2,h1,h2;
	double     vx1,vx2,vy1,vy2;
	GaussTria* gauss_1=NULL;
	GaussTria* gauss_2=NULL;

	/*Get material parameters :*/
	rho_ice=matpar->GetRhoIce();

	/*First off, check that this segment belongs to this element: */
	if ((int)*(segment+4)!=this->id)_error_("%s%i%s%i","error message: segment with id ",(int)*(segment+4)," does not belong to element with id:",this->id);

	/*Recover segment node locations: */
	x1=*(segment+0); y1=*(segment+1); x2=*(segment+2); y2=*(segment+3);

	/*Get xyz list: */
	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);

	/*get area coordinates of 0 and 1 locations: */
	gauss_1=new GaussTria();
	gauss_1->GaussFromCoords(x1,y1,&xyz_list[0][0]);
	gauss_2=new GaussTria();
	gauss_2->GaussFromCoords(x2,y2,&xyz_list[0][0]);

	normal[0]=cos(atan2(x1-x2,y2-y1));
	normal[1]=sin(atan2(x1-x2,y2-y1));

	length=sqrt(pow(x2-x1,2.0)+pow(y2-y1,2));

	Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
	this->parameters->FindParam(&dim,MeshDimensionEnum);
	Input* vx_input=NULL;
	Input* vy_input=NULL;
	if(dim==2){
		vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
		vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);
	}
	else{
		vx_input=inputs->GetInput(VxAverageEnum); _assert_(vx_input);
		vy_input=inputs->GetInput(VyAverageEnum); _assert_(vy_input);
	}

	thickness_input->GetInputValue(&h1, gauss_1);
	thickness_input->GetInputValue(&h2, gauss_2);
	vx_input->GetInputValue(&vx1,gauss_1);
	vx_input->GetInputValue(&vx2,gauss_2);
	vy_input->GetInputValue(&vy1,gauss_1);
	vy_input->GetInputValue(&vy2,gauss_2);

	mass_flux= rho_ice*length*(  
				(ONETHIRD*(h1-h2)*(vx1-vx2)+0.5*h2*(vx1-vx2)+0.5*(h1-h2)*vx2+h2*vx2)*normal[0]+
				(ONETHIRD*(h1-h2)*(vy1-vy2)+0.5*h2*(vy1-vy2)+0.5*(h1-h2)*vy2+h2*vy2)*normal[1]
				);

	/*Process units: */
	mass_flux=UnitConversion(mass_flux,IuToExtEnum,MassFluxEnum);

	/*clean up and return:*/
	delete gauss_1;
	delete gauss_2;
	return mass_flux;
}
/*}}}*/
/*FUNCTION Tria::MaxAbsVx{{{1*/
void  Tria::MaxAbsVx(double* pmaxabsvx, bool process_units){

	/*Get maximum:*/
	double maxabsvx=this->inputs->MaxAbs(VxEnum);

	/*process units if requested: */
	if(process_units) maxabsvx=UnitConversion(maxabsvx,IuToExtEnum,VxEnum);

	/*Assign output pointers:*/
	*pmaxabsvx=maxabsvx;
}
/*}}}*/
/*FUNCTION Tria::MaxAbsVy{{{1*/
void  Tria::MaxAbsVy(double* pmaxabsvy, bool process_units){

	/*Get maximum:*/
	double maxabsvy=this->inputs->MaxAbs(VyEnum);

	/*process units if requested: */
	if(process_units) maxabsvy=UnitConversion(maxabsvy,IuToExtEnum,VyEnum);

	/*Assign output pointers:*/
	*pmaxabsvy=maxabsvy;
}
/*}}}*/
/*FUNCTION Tria::MaxAbsVz{{{1*/
void  Tria::MaxAbsVz(double* pmaxabsvz, bool process_units){

	/*Get maximum:*/
	double maxabsvz=this->inputs->MaxAbs(VzEnum);

	/*process units if requested: */
	if(process_units) maxabsvz=UnitConversion(maxabsvz,IuToExtEnum,VyEnum);

	/*Assign output pointers:*/
	*pmaxabsvz=maxabsvz;
}
/*}}}*/
/*FUNCTION Tria::MaxVel{{{1*/
void  Tria::MaxVel(double* pmaxvel, bool process_units){

	/*Get maximum:*/
	double maxvel=this->inputs->Max(VelEnum);

	/*process units if requested: */
	if(process_units) maxvel=UnitConversion(maxvel,IuToExtEnum,VelEnum);

	/*Assign output pointers:*/
	*pmaxvel=maxvel;
}
/*}}}*/
/*FUNCTION Tria::MaxVx{{{1*/
void  Tria::MaxVx(double* pmaxvx, bool process_units){

	/*Get maximum:*/
	double maxvx=this->inputs->Max(VxEnum);

	/*process units if requested: */
	if(process_units) maxvx=UnitConversion(maxvx,IuToExtEnum,VxEnum);

	/*Assign output pointers:*/
	*pmaxvx=maxvx;
}
/*}}}*/
/*FUNCTION Tria::MaxVy{{{1*/
void  Tria::MaxVy(double* pmaxvy, bool process_units){

	/*Get maximum:*/
	double maxvy=this->inputs->Max(VyEnum);

	/*process units if requested: */
	if(process_units) maxvy=UnitConversion(maxvy,IuToExtEnum,VyEnum);

	/*Assign output pointers:*/
	*pmaxvy=maxvy;

}
/*}}}*/
/*FUNCTION Tria::MaxVz{{{1*/
void  Tria::MaxVz(double* pmaxvz, bool process_units){

	/*Get maximum:*/
	double maxvz=this->inputs->Max(VzEnum);

	/*process units if requested: */
	if(process_units) maxvz=UnitConversion(maxvz,IuToExtEnum,VzEnum);

	/*Assign output pointers:*/
	*pmaxvz=maxvz;
}
/*}}}*/
/*FUNCTION Tria::MinVel{{{1*/
void  Tria::MinVel(double* pminvel, bool process_units){

	/*Get minimum:*/
	double minvel=this->inputs->Min(VelEnum);

	/*process units if requested: */
	if(process_units) minvel=UnitConversion(minvel,IuToExtEnum,VelEnum);

	/*Assign output pointers:*/
	*pminvel=minvel;
}
/*}}}*/
/*FUNCTION Tria::MinVx{{{1*/
void  Tria::MinVx(double* pminvx, bool process_units){

	/*Get minimum:*/
	double minvx=this->inputs->Min(VxEnum);

	/*process units if requested: */
	if(process_units) minvx=UnitConversion(minvx,IuToExtEnum,VxEnum);

	/*Assign output pointers:*/
	*pminvx=minvx;
}
/*}}}*/
/*FUNCTION Tria::MinVy{{{1*/
void  Tria::MinVy(double* pminvy, bool process_units){

	/*Get minimum:*/
	double minvy=this->inputs->Min(VyEnum);

	/*process units if requested: */
	if(process_units) minvy=UnitConversion(minvy,IuToExtEnum,VyEnum);

	/*Assign output pointers:*/
	*pminvy=minvy;
}
/*}}}*/
/*FUNCTION Tria::MinVz{{{1*/
void  Tria::MinVz(double* pminvz, bool process_units){

	/*Get minimum:*/
	double minvz=this->inputs->Min(VzEnum);

	/*process units if requested: */
	if(process_units) minvz=UnitConversion(minvz,IuToExtEnum,VzEnum);

	/*Assign output pointers:*/
	*pminvz=minvz;
}
/*}}}*/
/*FUNCTION Tria::NodalValue {{{1*/
int    Tria::NodalValue(double* pvalue, int index, int natureofdataenum,bool process_units){

	int i;
	int found=0;
	double value;
	Input* data=NULL;
	GaussTria *gauss                            = NULL;

	/*First, serarch the input: */
	data=inputs->GetInput(natureofdataenum); 

	/*figure out if we have the vertex id: */
	found=0;
	for(i=0;i<NUMVERTICES;i++){
		if(index==nodes[i]->GetVertexId()){
			/*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;
			}
		}
	}

	if(found)*pvalue=value;
	return found;
}
/*}}}*/
/*FUNCTION Tria::ElementResponse{{{1*/
void Tria::ElementResponse(double* presponse,int response_enum,bool process_units){

	switch(response_enum){
		case MaterialsRheologyBbarEnum:
			*presponse=this->matice->GetBbar();
			break;
		case VelEnum:

			/*Get input:*/
			double vel;
			Input* vel_input;

			vel_input=this->inputs->GetInput(VelEnum); _assert_(vel_input);
			vel_input->GetInputAverage(&vel);

			/*process units if requested: */
			if(process_units) vel=UnitConversion(vel,IuToExtEnum,VelEnum);

			/*Assign output pointers:*/
			*presponse=vel;
		default:  
			_error_("Response type %s not supported yet!",EnumToStringx(response_enum));
	}

}
/*}}}*/
#endif

#ifdef _HAVE_DIAGNOSTIC_
/*FUNCTION Tria::CreateKMatrixDiagnosticMacAyeal {{{1*/
ElementMatrix* Tria::CreateKMatrixDiagnosticMacAyeal(void){

	/*compute all stiffness matrices for this element*/
	ElementMatrix* Ke1=CreateKMatrixDiagnosticMacAyealViscous();
	ElementMatrix* Ke2=CreateKMatrixDiagnosticMacAyealFriction();
	ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
	
	/*clean-up and return*/
	delete Ke1;
	delete Ke2;
	return Ke;
}
/*}}}*/
/*FUNCTION Tria::CreateKMatrixDiagnosticMacAyealViscous{{{1*/
ElementMatrix* Tria::CreateKMatrixDiagnosticMacAyealViscous(void){

	/*Constants*/
	const int  numdof=NDOF2*NUMVERTICES;

	/*Intermediaries*/
	int        i,j,ig;
	double     xyz_list[NUMVERTICES][3];
	double     viscosity,newviscosity,oldviscosity;
	double     viscosity_overshoot,thickness,Jdet;
	double     epsilon[3],oldepsilon[3];    /* epsilon=[exx,eyy,exy];    */
	double     B[3][numdof];
	double     Bprime[3][numdof];
	double     D[3][3]   = {0.0};
	double     D_scalar;
	GaussTria *gauss = NULL;

	/*Initialize Element matrix*/
	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);

	/*Retrieve all inputs and parameters*/
	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
	Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);
	Input* vx_input=inputs->GetInput(VxEnum);               _assert_(vx_input);
	Input* vy_input=inputs->GetInput(VyEnum);               _assert_(vy_input);
	Input* vxold_input=inputs->GetInput(VxPicardEnum);      _assert_(vxold_input);
	Input* vyold_input=inputs->GetInput(VyPicardEnum);      _assert_(vyold_input);
	this->parameters->FindParam(&viscosity_overshoot,DiagnosticViscosityOvershootEnum);

	/* Start  looping on the number of gaussian points: */
	gauss=new GaussTria(2);
	for (ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
		GetBMacAyeal(&B[0][0], &xyz_list[0][0], gauss);
		GetBprimeMacAyeal(&Bprime[0][0], &xyz_list[0][0], gauss);

		this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
		this->GetStrainRate2d(&oldepsilon[0],&xyz_list[0][0],gauss,vxold_input,vyold_input);
		matice->GetViscosity2d(&viscosity, &epsilon[0]);
		matice->GetViscosity2d(&oldviscosity, &oldepsilon[0]);
		thickness_input->GetInputValue(&thickness, gauss);

		newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
		D_scalar=2*newviscosity*thickness*gauss->weight*Jdet;
		for (i=0;i<3;i++) D[i][i]=D_scalar;

		TripleMultiply(&B[0][0],3,numdof,1,
					&D[0][0],3,3,0,
					&Bprime[0][0],3,numdof,0,
					&Ke->values[0],1);
	}

	/*Transform Coordinate System*/
	TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);

	/*Clean up and return*/
	delete gauss;
	return Ke;
}
/*}}}*/
/*FUNCTION Tria::CreateKMatrixDiagnosticMacAyealFriction {{{1*/
ElementMatrix* Tria::CreateKMatrixDiagnosticMacAyealFriction(void){

	/*Constants*/
	const int  numdof=NDOF2*NUMVERTICES;

	/*Intermediaries*/
	int        i,j,ig;
	int        analysis_type;
	double     MAXSLOPE  = .06; // 6 %
	double     MOUNTAINKEXPONENT = 10;
	double     slope_magnitude,alpha2;
	double     Jdet;
	double     L[2][numdof];
	double     DL[2][2]  = {{ 0,0 },{0,0}};
	double     DL_scalar;
	double     slope[2]  = {0.0,0.0};
	double     xyz_list[NUMVERTICES][3];
	Friction  *friction = NULL;
	GaussTria *gauss    = NULL;

	/*Initialize Element matrix and return if necessary*/
	if(IsFloating()) return NULL;
	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);

	/*Retrieve all inputs and parameters*/
	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
	Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input);
	Input* vx_input=inputs->GetInput(VxEnum);           _assert_(vx_input);
	Input* vy_input=inputs->GetInput(VyEnum);           _assert_(vy_input);
	Input* vz_input=inputs->GetInput(VzEnum);           _assert_(vz_input);
	parameters->FindParam(&analysis_type,AnalysisTypeEnum);

	/*build friction object, used later on: */
	friction=new Friction("2d",inputs,matpar,analysis_type);

	/* Start  looping on the number of gaussian points: */
	gauss=new GaussTria(2);
	for (ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		// If we have a slope > 6% for this element,  it means  we are on a mountain. In this particular case, 
		//velocity should be = 0. To achieve this result, we set alpha2_list to a very high value: */
		surface_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
		slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2));
		if(slope_magnitude>MAXSLOPE) alpha2=pow((double)10,MOUNTAINKEXPONENT);
		else friction->GetAlpha2(&alpha2, gauss,VxEnum,VyEnum,VzEnum);

		GetL(&L[0][0], &xyz_list[0][0], gauss,NDOF2);
		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
		DL_scalar=alpha2*gauss->weight*Jdet;
		for (i=0;i<2;i++) DL[i][i]=DL_scalar;
		
		TripleMultiply( &L[0][0],2,numdof,1,
					&DL[0][0],2,2,0,
					&L[0][0],2,numdof,0,
					&Ke->values[0],1);
	}

	/*Transform Coordinate System*/
	TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);

	/*Clean up and return*/
	delete gauss;
	delete friction;
	return Ke;
}
/*}}}*/
/*FUNCTION Tria::CreateKMatrixDiagnosticHutter{{{1*/
ElementMatrix* Tria::CreateKMatrixDiagnosticHutter(void){

	/*Intermediaries*/
	const int numdof=NUMVERTICES*NDOF2;
	int    i,connectivity;

	/*Initialize Element matrix*/
	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);

	/*Create Element matrix*/
	for(i=0;i<NUMVERTICES;i++){
		connectivity=nodes[i]->GetConnectivity();
		Ke->values[(2*i)*numdof  +(2*i)  ]=1/(double)connectivity;
		Ke->values[(2*i+1)*numdof+(2*i+1)]=1/(double)connectivity;
	}

	/*Clean up and return*/
	return Ke;
}
/*}}}*/
/*FUNCTION Tria::CreatePVectorDiagnosticMacAyeal {{{1*/
ElementVector* Tria::CreatePVectorDiagnosticMacAyeal(){

	/*Constants*/
	const int    numdof=NDOF2*NUMVERTICES;

	/*Intermediaries */
	int            i,j,ig;
	double         driving_stress_baseline,thickness;
	double         Jdet;
	double         xyz_list[NUMVERTICES][3];
	double         slope[2];
	double         basis[3];
	double         pe_g_gaussian[numdof];
	GaussTria*     gauss=NULL;

	/*Initialize Element vector*/
	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,MacAyealApproximationEnum);

	/*Retrieve all inputs and parameters*/
	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
	Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input); 
	Input* surface_input=inputs->GetInput(SurfaceEnum);     _assert_(surface_input);
	Input* drag_input=inputs->GetInput(FrictionCoefficientEnum);_assert_(drag_input);

	/* Start  looping on the number of gaussian points: */
	gauss=new GaussTria(2);
	for(ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
		GetNodalFunctions(basis, gauss);

		thickness_input->GetInputValue(&thickness,gauss);
		surface_input->GetInputDerivativeValue(&slope[0],&xyz_list[0][0],gauss);
		driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG()*thickness;

		/*Build pe_g_gaussian vector: */
		for (i=0;i<NUMVERTICES;i++){
			for (j=0;j<NDOF2;j++){
				pe->values[i*NDOF2+j]+=-driving_stress_baseline*slope[j]*Jdet*gauss->weight*basis[i];
			}
		}
	}

	/*Transform coordinate system*/
	TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYEnum);

	/*Clean up and return*/
	delete gauss;
	return pe;
}
/*}}}*/
/*FUNCTION Tria::CreatePVectorDiagnosticHutter{{{1*/
ElementVector* Tria::CreatePVectorDiagnosticHutter(void){

	/*Intermediaries */
	int        i,connectivity;
	double     constant_part,ub,vb;
	double     rho_ice,gravity,n,B;
	double     slope2,thickness;
	double     slope[2];
	GaussTria* gauss=NULL;

	/*Initialize Element vector*/
	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);

	/*Retrieve all inputs and parameters*/
	rho_ice=matpar->GetRhoIce();
	gravity=matpar->GetG();
	n=matice->GetN();
	B=matice->GetBbar();
	Input* slopex_input=inputs->GetInput(SurfaceSlopeXEnum); _assert_(slopex_input);
	Input* slopey_input=inputs->GetInput(SurfaceSlopeYEnum); _assert_(slopey_input);
	Input* thickness_input=inputs->GetInput(ThicknessEnum);  _assert_(thickness_input);

	/*Spawn 3 sing elements: */
	gauss=new GaussTria();
	for(i=0;i<NUMVERTICES;i++){

		gauss->GaussVertex(i);

		connectivity=nodes[i]->GetConnectivity();

		thickness_input->GetInputValue(&thickness,gauss);
		slopex_input->GetInputValue(&slope[0],gauss);
		slopey_input->GetInputValue(&slope[1],gauss);
		slope2=pow(slope[0],2)+pow(slope[1],2);

		constant_part=-2*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2));

		ub=-1.58*pow((double)10.0,(double)-10.0)*rho_ice*gravity*thickness*slope[0];
		vb=-1.58*pow((double)10.0,(double)-10.0)*rho_ice*gravity*thickness*slope[1];

		pe->values[2*i]  =(ub-2.0*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2.0))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[0])/(double)connectivity;
		pe->values[2*i+1]=(vb-2.0*pow(rho_ice*gravity,n)*pow(slope2,((n-1)/2.0))*pow(thickness,n)/(pow(B,n)*(n+1))*slope[1])/(double)connectivity;
	}

	/*Clean up and return*/
	delete gauss;
	return pe;
}
/*}}}*/
/*FUNCTION Tria::CreateJacobianDiagnosticMacayeal{{{1*/
ElementMatrix* Tria::CreateJacobianDiagnosticMacayeal(void){

	/*Constants*/
	const int    numdof=NDOF2*NUMVERTICES;

	/*Intermediaries */
	int        i,j,ig;
	double     xyz_list[NUMVERTICES][3];
	double     Jdet,thickness;
	double     eps1dotdphii,eps1dotdphij;
	double     eps2dotdphii,eps2dotdphij;
	double     mu_prime;
	double     epsilon[3];/* epsilon=[exx,eyy,exy];*/
	double     eps1[2],eps2[2];
	double     phi[NUMVERTICES];
	double     dphi[2][NUMVERTICES];
	GaussTria *gauss=NULL;

	/*Initialize Jacobian with regular MacAyeal (first part of the Gateau derivative)*/
	ElementMatrix* Ke=CreateKMatrixDiagnosticMacAyeal();

	/*Retrieve all inputs and parameters*/
	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
	Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
	Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
	Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);

	/* Start  looping on the number of gaussian points: */
	gauss=new GaussTria(2);
	for (ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
		GetNodalFunctionsDerivatives(&dphi[0][0],&xyz_list[0][0],gauss);

		thickness_input->GetInputValue(&thickness, gauss);
		this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
		matice->GetViscosity2dDerivativeEpsSquare(&mu_prime,&epsilon[0]);
		eps1[0]=2*epsilon[0]+epsilon[1];   eps2[0]=epsilon[2];
		eps1[1]=epsilon[2];                eps2[1]=epsilon[0]+2*epsilon[1];

		for(i=0;i<3;i++){
			for(j=0;j<3;j++){
				eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i];
				eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j];
				eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i];
				eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j];

				Ke->values[6*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps1dotdphii;
				Ke->values[6*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps1dotdphii;
				Ke->values[6*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps2dotdphii;
				Ke->values[6*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps2dotdphii;
			}
		}
	}

	/*Transform Coordinate System*/
	TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);

	/*Clean up and return*/
	delete gauss;
	return Ke;
}
/*}}}*/
/*FUNCTION Tria::GetSolutionFromInputsDiagnosticHoriz{{{1*/
void  Tria::GetSolutionFromInputsDiagnosticHoriz(Vec solution){

	const int    numdof=NDOF2*NUMVERTICES;

	int          i;
	int*         doflist=NULL;
	double       vx,vy;
	double       values[numdof];
	GaussTria*   gauss=NULL;

	/*Get dof list: */
	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);

	/*Get inputs*/
	Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
	Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);

	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
	/*P1 element only for now*/
	gauss=new GaussTria();
	for(i=0;i<NUMVERTICES;i++){

		gauss->GaussVertex(i);

		/*Recover vx and vy*/
		vx_input->GetInputValue(&vx,gauss);
		vy_input->GetInputValue(&vy,gauss);
		values[i*NDOF2+0]=vx;
		values[i*NDOF2+1]=vy;
	}

	VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES);

	/*Free ressources:*/
	delete gauss;
	xfree((void**)&doflist);
}
/*}}}*/
/*FUNCTION Tria::GetSolutionFromInputsDiagnosticHutter{{{1*/
void  Tria::GetSolutionFromInputsDiagnosticHutter(Vec solution){

	const int    numdof=NDOF2*NUMVERTICES;

	int        i;
	double     vx,vy;
	double     values[numdof];
	int       *doflist = NULL;
	GaussTria *gauss   = NULL;

	/*Get dof list: */
	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);

	/*Get inputs*/
	Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
	Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);

	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
	/*P1 element only for now*/
	gauss=new GaussTria();
	for(i=0;i<NUMVERTICES;i++){

		gauss->GaussVertex(i);

		/*Recover vx and vy*/
		vx_input->GetInputValue(&vx,gauss);
		vy_input->GetInputValue(&vy,gauss);
		values[i*NDOF2+0]=vx;
		values[i*NDOF2+1]=vy;
	}

	VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES);

	/*Free ressources:*/
	delete gauss;
	xfree((void**)&doflist);
}
/*}}}*/
/*FUNCTION Tria::InputUpdateFromSolutionDiagnosticHoriz {{{1*/
void  Tria::InputUpdateFromSolutionDiagnosticHoriz(double* solution){
	
	const int numdof=NDOF2*NUMVERTICES;

	int       i;
	int*      doflist=NULL;
	double    rho_ice,g;
	double    values[numdof];
	double    vx[NUMVERTICES];
	double    vy[NUMVERTICES];
	double    vz[NUMVERTICES];
	double    vel[NUMVERTICES];
	double    pressure[NUMVERTICES];
	double    thickness[NUMVERTICES];
	
	/*Get dof list: */
	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);

	/*Use the dof list to index into the solution vector: */
	for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];

	/*Transform solution in Cartesian Space*/
	TransformSolutionCoord(&values[0],nodes,NUMVERTICES,XYEnum);

	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
	for(i=0;i<NUMVERTICES;i++){
		vx[i]=values[i*NDOF2+0];
		vy[i]=values[i*NDOF2+1];

		/*Check solution*/
		if(isnan(vx[i])) _error_("NaN found in solution vector");
		if(isnan(vy[i])) _error_("NaN found in solution vector");
	}

	/*Get Vz and compute vel*/
	GetInputListOnVertices(&vz[0],VzEnum,0);
	for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);

	/*For pressure: we have not computed pressure in this analysis, for this element. We are in 2D, 
	 *so the pressure is just the pressure at the bedrock: */
	rho_ice=matpar->GetRhoIce();
	g=matpar->GetG();
	GetInputListOnVertices(&thickness[0],ThicknessEnum);
	for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*thickness[i];

	/*Now, we have to move the previous Vx and Vy inputs  to old 
	 * status, otherwise, we'll wipe them off: */
	this->inputs->ChangeEnum(VxEnum,VxPicardEnum);
	this->inputs->ChangeEnum(VyEnum,VyPicardEnum);
	this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);

	/*Add vx and vy as inputs to the tria element: */
	this->inputs->AddInput(new TriaP1Input(VxEnum,vx));
	this->inputs->AddInput(new TriaP1Input(VyEnum,vy));
	this->inputs->AddInput(new TriaP1Input(VelEnum,vel));
	this->inputs->AddInput(new TriaP1Input(PressureEnum,pressure));

	/*Free ressources:*/
	xfree((void**)&doflist);

}
/*}}}*/
/*FUNCTION Tria::InputUpdateFromSolutionDiagnosticHutter {{{1*/
void  Tria::InputUpdateFromSolutionDiagnosticHutter(double* solution){
	
	const int numdof=NDOF2*NUMVERTICES;
	
	int       i;
	int*      doflist=NULL;
	double    rho_ice,g;
	double    values[numdof];
	double    vx[NUMVERTICES];
	double    vy[NUMVERTICES];
	double    vz[NUMVERTICES];
	double    vel[NUMVERTICES];
	double    pressure[NUMVERTICES];
	double    thickness[NUMVERTICES];
	
	/*Get dof list: */
	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);

	/*Use the dof list to index into the solution vector: */
	for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];

	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
	for(i=0;i<NUMVERTICES;i++){
		vx[i]=values[i*NDOF2+0];
		vy[i]=values[i*NDOF2+1];

		/*Check solution*/
		if(isnan(vx[i])) _error_("NaN found in solution vector");
		if(isnan(vy[i])) _error_("NaN found in solution vector");
	}

	/*Now Compute vel*/
	GetInputListOnVertices(&vz[0],VzEnum,0.0); //default is 0
	for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);

	/*For pressure: we have not computed pressure in this analysis, for this element. We are in 2D, 
	 *so the pressure is just the pressure at the bedrock: */
	rho_ice=matpar->GetRhoIce();
	g=matpar->GetG();
	GetInputListOnVertices(&thickness[0],ThicknessEnum);
	for(i=0;i<NUMVERTICES;i++) pressure[i]=rho_ice*g*thickness[i];

	/*Now, we have to move the previous Vx and Vy inputs  to old 
	 * status, otherwise, we'll wipe them off: */
	this->inputs->ChangeEnum(VxEnum,VxPicardEnum);
	this->inputs->ChangeEnum(VyEnum,VyPicardEnum);
	this->inputs->ChangeEnum(PressureEnum,PressurePicardEnum);

	/*Add vx and vy as inputs to the tria element: */
	this->inputs->AddInput(new TriaP1Input(VxEnum,vx));
	this->inputs->AddInput(new TriaP1Input(VyEnum,vy));
	this->inputs->AddInput(new TriaP1Input(VelEnum,vel));
	this->inputs->AddInput(new TriaP1Input(PressureEnum,pressure));

	/*Free ressources:*/
	xfree((void**)&doflist);
}
/*}}}*/
#endif

#ifdef _HAVE_CONTROL_
/*FUNCTION Tria::InputControlUpdate{{{1*/
void  Tria::InputControlUpdate(double scalar,bool save_parameter){

	/*Intermediary*/
	int    num_controls;
	int*   control_type=NULL;
	Input* input=NULL;

	/*retrieve some parameters: */
	this->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
	this->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);

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

		if(control_type[i]==MaterialsRheologyBbarEnum){
			input=(Input*)matice->inputs->GetInput(control_type[i]); _assert_(input);
		}
		else{
			input=(Input*)this->inputs->GetInput(control_type[i]);   _assert_(input);
		}

		if (input->ObjectEnum()!=ControlInputEnum){
			_error_("input %s is not a ControlInput",EnumToStringx(control_type[i]));
		}

		((ControlInput*)input)->UpdateValue(scalar);
		((ControlInput*)input)->Constrain();
		if (save_parameter) ((ControlInput*)input)->SaveValue();

	}

	/*Clean up and return*/
	xfree((void**)&control_type);
}
/*}}}*/
/*FUNCTION Tria::ControlInputGetGradient{{{1*/
void Tria::ControlInputGetGradient(Vec gradient,int enum_type,int control_index){

	int doflist1[NUMVERTICES];
	Input* input=NULL;

	if(enum_type==MaterialsRheologyBbarEnum){
		input=(Input*)matice->inputs->GetInput(enum_type);
	}
	else{
		input=inputs->GetInput(enum_type);
	}
	if (!input) _error_("Input %s not found",EnumToStringx(enum_type));
	if (input->ObjectEnum()!=ControlInputEnum) _error_("Input %s is not a ControlInput",EnumToStringx(enum_type));

	GradientIndexing(&doflist1[0],control_index);
	((ControlInput*)input)->GetGradient(gradient,&doflist1[0]);

}/*}}}*/
/*FUNCTION Tria::ControlInputScaleGradient{{{1*/
void Tria::ControlInputScaleGradient(int enum_type,double scale){

	Input* input=NULL;

	if(enum_type==MaterialsRheologyBbarEnum){
		input=(Input*)matice->inputs->GetInput(enum_type);
	}
	else{
		input=inputs->GetInput(enum_type);
	}
	if (!input) _error_("Input %s not found",EnumToStringx(enum_type));
	if (input->ObjectEnum()!=ControlInputEnum) _error_("Input %s is not a ControlInput",EnumToStringx(enum_type));

	((ControlInput*)input)->ScaleGradient(scale);
}/*}}}*/
/*FUNCTION Tria::ControlInputSetGradient{{{1*/
void Tria::ControlInputSetGradient(double* gradient,int enum_type,int control_index){

	int    doflist1[NUMVERTICES];
	double grad_list[NUMVERTICES];
	Input* grad_input=NULL;
	Input* input=NULL;

	if(enum_type==MaterialsRheologyBbarEnum){
		input=(Input*)matice->inputs->GetInput(enum_type);
	}
	else{
		input=inputs->GetInput(enum_type);
	}
	if (!input) _error_("Input %s not found",EnumToStringx(enum_type));
	if (input->ObjectEnum()!=ControlInputEnum) _error_("Input %s is not a ControlInput",EnumToStringx(enum_type));

	GradientIndexing(&doflist1[0],control_index);
	for(int i=0;i<NUMVERTICES;i++) grad_list[i]=gradient[doflist1[i]];
	grad_input=new TriaP1Input(GradientEnum,grad_list);

	((ControlInput*)input)->SetGradient(grad_input);

}/*}}}*/
/*FUNCTION Tria::Gradj {{{1*/
void  Tria::Gradj(Vec gradient,int control_type,int control_index){
	/*dJ/dalpha = ∂L/∂alpha = ∂J/∂alpha + ∂/∂alpha(KU-F)*/

	/*If on water, grad = 0: */
	if(IsOnWater()) return;

	/*First deal with ∂/∂alpha(KU-F)*/
	switch(control_type){
		case FrictionCoefficientEnum:
			GradjDragMacAyeal(gradient,control_index);
			break;
		case MaterialsRheologyBbarEnum:
			GradjBMacAyeal(gradient,control_index);
			break;
		case BalancethicknessThickeningRateEnum:
			GradjDhDtBalancedthickness(gradient,control_index);
			break;
		case VxEnum:
			GradjVxBalancedthickness(gradient,control_index);
			break;
		case VyEnum:
			GradjVyBalancedthickness(gradient,control_index);
			break;
		default:
			_error_("%s%i","control type not supported yet: ",control_type);
	}

	/*Now deal with ∂J/∂alpha*/
	int        *responses = NULL;
	int         num_responses,resp;
	this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
	this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum);

	for(resp=0;resp<num_responses;resp++) switch(responses[resp]){
		//FIXME: the control type should be checked somewhere (with respect to what variable are we taking the gradient!)

		case ThicknessAbsMisfitEnum:
		case ThicknessAbsGradientEnum:
		case SurfaceAbsVelMisfitEnum:
		case SurfaceRelVelMisfitEnum:
		case SurfaceLogVelMisfitEnum:
		case SurfaceLogVxVyMisfitEnum:
		case SurfaceAverageVelMisfitEnum:
			/*Nothing, J does not depends on the parameter being inverted for*/
			break;
		case DragCoefficientAbsGradientEnum:
			GradjDragGradient(gradient,resp,control_index);
			break;
		case RheologyBbarAbsGradientEnum:
			GradjBGradient(gradient,resp,control_index);
			break;
		default:
			_error_("response %s not supported yet",EnumToStringx(responses[resp]));
	}

	xfree((void**)&responses);
}
/*}}}*/
/*FUNCTION Tria::GradjBGradient{{{1*/
void  Tria::GradjBGradient(Vec gradient,int weight_index,int control_index){

	int        i,ig;
	int        doflist1[NUMVERTICES];
	double     Jdet,weight;
	double     xyz_list[NUMVERTICES][3];
	double     dbasis[NDOF2][NUMVERTICES];
	double     dk[NDOF2]; 
	double     grade_g[NUMVERTICES]={0.0};
	GaussTria  *gauss=NULL;

	/*Retrieve all inputs we will be needing: */
	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
	GradientIndexing(&doflist1[0],control_index);
	Input* rheologyb_input=matice->inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);
	Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum);                _assert_(weights_input);

	/* Start  looping on the number of gaussian points: */
	gauss=new GaussTria(2);
	for (ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
		GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
		weights_input->GetInputValue(&weight,gauss,weight_index);

		/*Build alpha_complement_list: */
		rheologyb_input->GetInputDerivativeValue(&dk[0],&xyz_list[0][0],gauss);

		/*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
		for (i=0;i<NUMVERTICES;i++) grade_g[i]+=-weight*Jdet*gauss->weight*(dbasis[0][i]*dk[0]+dbasis[1][i]*dk[1]);
	}
	VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES);

	/*Clean up and return*/
	delete gauss;
}
/*}}}*/
/*FUNCTION Tria::GradjBMacAyeal{{{1*/
void  Tria::GradjBMacAyeal(Vec gradient,int control_index){

	/*Intermediaries*/
	int        i,ig;
	int        doflist[NUMVERTICES];
	double     vx,vy,lambda,mu,thickness,Jdet;
	double     viscosity_complement;
	double     dvx[NDOF2],dvy[NDOF2],dadjx[NDOF2],dadjy[NDOF2],dB[NDOF2]; 
	double     xyz_list[NUMVERTICES][3];
	double     basis[3],epsilon[3];
	double     grad[NUMVERTICES]={0.0};
	GaussTria *gauss = NULL;

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

	/*Retrieve all inputs*/
	Input* thickness_input=inputs->GetInput(ThicknessEnum);                     _assert_(thickness_input);
	Input* vx_input=inputs->GetInput(VxEnum);                                   _assert_(vx_input);
	Input* vy_input=inputs->GetInput(VyEnum);                                   _assert_(vy_input);
	Input* adjointx_input=inputs->GetInput(AdjointxEnum);                       _assert_(adjointx_input);
	Input* adjointy_input=inputs->GetInput(AdjointyEnum);                       _assert_(adjointy_input);
	Input* rheologyb_input=matice->inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);

	/* Start  looping on the number of gaussian points: */
	gauss=new GaussTria(4);
	for (ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		thickness_input->GetInputValue(&thickness,gauss);
		rheologyb_input->GetInputDerivativeValue(&dB[0],&xyz_list[0][0],gauss);
		vx_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
		vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
		adjointx_input->GetInputDerivativeValue(&dadjx[0],&xyz_list[0][0],gauss);
		adjointy_input->GetInputDerivativeValue(&dadjy[0],&xyz_list[0][0],gauss);

		this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
		matice->GetViscosityComplement(&viscosity_complement,&epsilon[0]);

		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
		GetNodalFunctions(basis,gauss);

		/*standard gradient dJ/dki*/
		for (i=0;i<NUMVERTICES;i++) grad[i]+=-viscosity_complement*thickness*(
					(2*dvx[0]+dvy[1])*2*dadjx[0]+(dvx[1]+dvy[0])*(dadjx[1]+dadjy[0])+(2*dvy[1]+dvx[0])*2*dadjy[1]
					)*Jdet*gauss->weight*basis[i];
	}

	VecSetValues(gradient,NUMVERTICES,doflist,(const double*)grad,ADD_VALUES);

	/*clean-up*/
	delete gauss;
}
/*}}}*/
/*FUNCTION Tria::GradjDragMacAyeal {{{1*/
void  Tria::GradjDragMacAyeal(Vec gradient,int control_index){

	int        i,ig;
	int        analysis_type;
	int        doflist1[NUMVERTICES];
	double     vx,vy,lambda,mu,alpha_complement,Jdet;
	double     bed,thickness,Neff,drag;
	double     xyz_list[NUMVERTICES][3];
	double     dk[NDOF2]; 
	double     grade_g[NUMVERTICES]={0.0};
	double     grade_g_gaussian[NUMVERTICES];
	double     basis[3];
	double     epsilon[3]; /* epsilon=[exx,eyy,exy];*/
	Friction*  friction=NULL;
	GaussTria  *gauss=NULL;

	if(IsFloating())return;

	/*retrive parameters: */
	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
	GradientIndexing(&doflist1[0],control_index);

	/*Build frictoin element, needed later: */
	friction=new Friction("2d",inputs,matpar,analysis_type);

	/*Retrieve all inputs we will be needing: */
	Input* adjointx_input=inputs->GetInput(AdjointxEnum);                   _assert_(adjointx_input);
	Input* adjointy_input=inputs->GetInput(AdjointyEnum);                   _assert_(adjointy_input);
	Input* vx_input=inputs->GetInput(VxEnum);                               _assert_(vx_input);
	Input* vy_input=inputs->GetInput(VyEnum);                               _assert_(vy_input);
	Input* dragcoefficient_input=inputs->GetInput(FrictionCoefficientEnum); _assert_(dragcoefficient_input);

	/* Start  looping on the number of gaussian points: */
	gauss=new GaussTria(4);
	for (ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
		GetNodalFunctions(basis, gauss);

		/*Build alpha_complement_list: */
		friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum,VzEnum);
	
		dragcoefficient_input->GetInputValue(&drag, gauss);
		adjointx_input->GetInputValue(&lambda, gauss);
		adjointy_input->GetInputValue(&mu, gauss);
		vx_input->GetInputValue(&vx,gauss);
		vy_input->GetInputValue(&vy,gauss);
		dragcoefficient_input->GetInputDerivativeValue(&dk[0],&xyz_list[0][0],gauss);

		/*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
		for (i=0;i<NUMVERTICES;i++){
			grade_g_gaussian[i]=-2*drag*alpha_complement*((lambda*vx+mu*vy))*Jdet*gauss->weight*basis[i];
		}
		
		/*Add gradje_g_gaussian vector to gradje_g: */
		for(i=0;i<NUMVERTICES;i++){
			_assert_(!isnan(grade_g[i]));
			grade_g[i]+=grade_g_gaussian[i];
		}
	}
	/*Analytical gradient*/
	//delete gauss;
	//gauss=new GaussTria();
	//for (int iv=0;iv<NUMVERTICES;iv++){
	//	gauss->GaussVertex(iv);
	//	friction->GetAlphaComplement(&alpha_complement, gauss,VxEnum,VyEnum,VzEnum);
	//	dragcoefficient_input->GetInputValue(&drag, gauss);
	//	adjointx_input->GetInputValue(&lambda, gauss);
	//	adjointy_input->GetInputValue(&mu, gauss);
	//	vx_input->GetInputValue(&vx,gauss);
	//	vy_input->GetInputValue(&vy,gauss);
	//	grade_g[iv]=-2*drag*alpha_complement*((lambda*vx+mu*vy));
	//	VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,INSERT_VALUES);
	//}
	/*End Analytical gradient*/

	VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES);

	/*Clean up and return*/
	delete gauss;
	delete friction;
}
/*}}}*/
/*FUNCTION Tria::GradjDragGradient{{{1*/
void  Tria::GradjDragGradient(Vec gradient, int weight_index,int control_index){

	int        i,ig;
	int        doflist1[NUMVERTICES];
	double     Jdet,weight;
	double     xyz_list[NUMVERTICES][3];
	double     dbasis[NDOF2][NUMVERTICES];
	double     dk[NDOF2]; 
	double     grade_g[NUMVERTICES]={0.0};
	GaussTria  *gauss=NULL;

	/*Retrieve all inputs we will be needing: */
	if(IsFloating())return;
	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
	GradientIndexing(&doflist1[0],control_index);
	Input* dragcoefficient_input=inputs->GetInput(FrictionCoefficientEnum); _assert_(dragcoefficient_input);
	Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum);                 _assert_(weights_input);

	/* Start  looping on the number of gaussian points: */
	gauss=new GaussTria(2);
	for (ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
		GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);
		weights_input->GetInputValue(&weight,gauss,weight_index);

		/*Build alpha_complement_list: */
		dragcoefficient_input->GetInputDerivativeValue(&dk[0],&xyz_list[0][0],gauss);

		/*Build gradje_g_gaussian vector (actually -dJ/ddrag): */
		for (i=0;i<NUMVERTICES;i++){
			grade_g[i]+=-weight*Jdet*gauss->weight*(dbasis[0][i]*dk[0]+dbasis[1][i]*dk[1]);
			_assert_(!isnan(grade_g[i]));
		}
	}
	VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES);

	/*Clean up and return*/
	delete gauss;
}
/*}}}*/
/*FUNCTION Tria::GradjDhDtBalancedthickness{{{1*/
void  Tria::GradjDhDtBalancedthickness(Vec gradient,int control_index){

	/*Intermediaries*/
	int    doflist1[NUMVERTICES];
	double lambda[NUMVERTICES];
	double gradient_g[NUMVERTICES];

	/*Compute Gradient*/
	GradientIndexing(&doflist1[0],control_index);
	GetInputListOnVertices(&lambda[0],AdjointEnum);
	for(int i=0;i<NUMVERTICES;i++) gradient_g[i]=-lambda[i];

	VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)gradient_g,INSERT_VALUES);
}
/*}}}*/
/*FUNCTION Tria::GradjVxBalancedthickness{{{1*/
void  Tria::GradjVxBalancedthickness(Vec gradient,int control_index){

	/*Intermediaries*/
	int        i,ig;
	int        doflist1[NUMVERTICES];
	double     thickness,Jdet;
	double     basis[3];
	double     Dlambda[2],dp[2];
	double     xyz_list[NUMVERTICES][3];
	double     grade_g[NUMVERTICES] = {0.0};
	GaussTria *gauss                = NULL;

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

	/*Retrieve all inputs we will be needing: */
	Input* adjoint_input=inputs->GetInput(AdjointEnum);     _assert_(adjoint_input);
	Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);

	/* Start  looping on the number of gaussian points: */
	gauss=new GaussTria(2);
	for (ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
		GetNodalFunctions(basis, gauss);
		
		adjoint_input->GetInputDerivativeValue(&Dlambda[0],&xyz_list[0][0],gauss);
		thickness_input->GetInputValue(&thickness, gauss);
		thickness_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss);

		for(i=0;i<NUMVERTICES;i++) grade_g[i]+=thickness*Dlambda[0]*Jdet*gauss->weight*basis[i];
	}

	VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES);

	/*Clean up and return*/
	delete gauss;
}
/*}}}*/
/*FUNCTION Tria::GradjVyBalancedthickness{{{1*/
void  Tria::GradjVyBalancedthickness(Vec gradient,int control_index){

	/*Intermediaries*/
	int        i,ig;
	int        doflist1[NUMVERTICES];
	double     thickness,Jdet;
	double     basis[3];
	double     Dlambda[2],dp[2];
	double     xyz_list[NUMVERTICES][3];
	double     grade_g[NUMVERTICES] = {0.0};
	GaussTria *gauss                = NULL;

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

	/*Retrieve all inputs we will be needing: */
	Input* adjoint_input=inputs->GetInput(AdjointEnum);     _assert_(adjoint_input);
	Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);

	/* Start  looping on the number of gaussian points: */
	gauss=new GaussTria(2);
	for (ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
		GetNodalFunctions(basis, gauss);

		adjoint_input->GetInputDerivativeValue(&Dlambda[0],&xyz_list[0][0],gauss);
		thickness_input->GetInputValue(&thickness, gauss);
		thickness_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss);

		for(i=0;i<NUMVERTICES;i++) grade_g[i]+=thickness*Dlambda[1]*Jdet*gauss->weight*basis[i];
	}
	VecSetValues(gradient,NUMVERTICES,doflist1,(const double*)grade_g,ADD_VALUES);

	/*Clean up and return*/
	delete gauss;
}
/*}}}*/
/*FUNCTION Tria::GradientIndexing{{{1*/
void  Tria::GradientIndexing(int* indexing,int control_index){

	/*Get some parameters*/
	int num_controls;
	parameters->FindParam(&num_controls,InversionNumControlParametersEnum);

	/*get gradient indices*/
	for(int i=0;i<NUMVERTICES;i++){
		indexing[i]=num_controls*this->nodes[i]->GetVertexDof() + control_index;
	}

}
/*}}}*/
/*FUNCTION Tria::RheologyBbarAbsGradient{{{1*/
double Tria::RheologyBbarAbsGradient(bool process_units,int weight_index){

	/* Intermediaries */
	int        ig;
	double     Jelem = 0;
	double     weight;
	double     Jdet;
	double     xyz_list[NUMVERTICES][3];
	double     dp[NDOF2];
	GaussTria *gauss = NULL;

	/*retrieve parameters and inputs*/

	/*If on water, return 0: */
	if(IsOnWater()) return 0;

	/*Retrieve all inputs we will be needing: */
	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
	Input* weights_input  =inputs->GetInput(InversionCostFunctionsCoefficientsEnum);              _assert_(weights_input);
	Input* rheologyb_input=matice->inputs->GetInput(MaterialsRheologyBbarEnum); _assert_(rheologyb_input);

	/* Start looping on the number of gaussian points: */
	gauss=new GaussTria(2);
	for (ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		/* Get Jacobian determinant: */
		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);

		/*Get all parameters at gaussian point*/
		weights_input->GetInputValue(&weight,gauss,weight_index);
		rheologyb_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss);

		/*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ 
		Jelem+=weight*1/2*(pow(dp[0],2.)+pow(dp[1],2.))*Jdet*gauss->weight;
	}

	/*Clean up and return*/
	delete gauss;
	return Jelem;
}
/*}}}*/
/*FUNCTION Tria::SurfaceAverageVelMisfit {{{1*/
double Tria::SurfaceAverageVelMisfit(bool process_units,int weight_index){

	const int    numdof=2*NUMVERTICES;

	int        i,ig;
	double     Jelem=0,S,Jdet;
	double     misfit;
	double     vx,vy,vxobs,vyobs,weight;
	double     xyz_list[NUMVERTICES][3];
	GaussTria *gauss=NULL;

	/*If on water, return 0: */
	if(IsOnWater())return 0;

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

	/*Retrieve all inputs we will be needing: */
	inputs->GetInputValue(&S,SurfaceAreaEnum);
	Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum);   _assert_(weights_input);
	Input* vx_input     =inputs->GetInput(VxEnum);        _assert_(vx_input);
	Input* vy_input     =inputs->GetInput(VyEnum);        _assert_(vy_input);
	Input* vxobs_input  =inputs->GetInput(InversionVxObsEnum);     _assert_(vxobs_input);
	Input* vyobs_input  =inputs->GetInput(InversionVyObsEnum);     _assert_(vyobs_input);

	/* Start  looping on the number of gaussian points: */
	gauss=new GaussTria(3);
	for (ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		/* Get Jacobian determinant: */
		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);

		/*Get all parameters at gaussian point*/
		weights_input->GetInputValue(&weight,gauss,weight_index);
		vx_input->GetInputValue(&vx,gauss);
		vy_input->GetInputValue(&vy,gauss);
		vxobs_input->GetInputValue(&vxobs,gauss);
		vyobs_input->GetInputValue(&vyobs,gauss);

		/*Compute SurfaceAverageVelMisfitEnum:
		 *
		 *      1                    2              2
		 * J = ---  sqrt(  (u - u   )  +  (v - v   )  )
		 *      S                obs            obs
		 */
		misfit=1/S*pow( pow(vx-vxobs,2.) + pow(vy-vyobs,2.) ,0.5);

		if(process_units)UnitConversion(misfit,IuToExtEnum,SurfaceAverageVelMisfitEnum);

		/*Add to cost function*/
		Jelem+=misfit*weight*Jdet*gauss->weight;
	}

	/*clean-up and Return: */
	delete gauss;
	return Jelem;
}
/*}}}*/
/*FUNCTION Tria::SurfaceLogVelMisfit {{{1*/
double Tria::SurfaceLogVelMisfit(bool process_units,int weight_index){

	const int    numdof=NDOF2*NUMVERTICES;

	int        i,ig;
	double     Jelem=0;
	double     misfit,Jdet;
	double     epsvel=2.220446049250313e-16;
	double     meanvel=3.170979198376458e-05; /*1000 m/yr*/
	double     velocity_mag,obs_velocity_mag;
	double     xyz_list[NUMVERTICES][3];
	double     vx,vy,vxobs,vyobs,weight;
	GaussTria *gauss=NULL;

	/*If on water, return 0: */
	if(IsOnWater())return 0;

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

	/*Retrieve all inputs we will be needing: */
	Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum);   _assert_(weights_input);
	Input* vx_input     =inputs->GetInput(VxEnum);        _assert_(vx_input);
	Input* vy_input     =inputs->GetInput(VyEnum);        _assert_(vy_input);
	Input* vxobs_input  =inputs->GetInput(InversionVxObsEnum);     _assert_(vxobs_input);
	Input* vyobs_input  =inputs->GetInput(InversionVyObsEnum);     _assert_(vyobs_input);

	/* Start  looping on the number of gaussian points: */
	gauss=new GaussTria(4);
	for (ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		/* Get Jacobian determinant: */
		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);

		/*Get all parameters at gaussian point*/
		weights_input->GetInputValue(&weight,gauss,weight_index);
		vx_input->GetInputValue(&vx,gauss);
		vy_input->GetInputValue(&vy,gauss);
		vxobs_input->GetInputValue(&vxobs,gauss);
		vyobs_input->GetInputValue(&vyobs,gauss);

		/*Compute SurfaceLogVelMisfit:
		 *                 [        vel + eps     ] 2
		 * J = 4 \bar{v}^2 | log ( -----------  ) |  
		 *                 [       vel   + eps    ]
		 *                            obs
		 */
		velocity_mag    =sqrt(pow(vx,   2.)+pow(vy,   2.))+epsvel;
		obs_velocity_mag=sqrt(pow(vxobs,2.)+pow(vyobs,2.))+epsvel;
		misfit=4*pow(meanvel,2.)*pow(log(velocity_mag/obs_velocity_mag),2.);

		if(process_units)UnitConversion(misfit,IuToExtEnum,SurfaceLogVelMisfitEnum);

		/*Add to cost function*/
		Jelem+=misfit*weight*Jdet*gauss->weight;
	}

	/*clean-up and Return: */
	delete gauss;
	return Jelem;
}
/*}}}*/
/*FUNCTION Tria::SurfaceLogVxVyMisfit {{{1*/
double Tria::SurfaceLogVxVyMisfit(bool process_units,int weight_index){

	const int    numdof=NDOF2*NUMVERTICES;

	int        i,ig;
	int        fit=-1;
	double     Jelem=0, S=0;
	double     epsvel=2.220446049250313e-16;
	double     meanvel=3.170979198376458e-05; /*1000 m/yr*/
	double     misfit, Jdet;
	double     vx,vy,vxobs,vyobs,weight;
	double     xyz_list[NUMVERTICES][3];
	GaussTria *gauss=NULL;

	/*If on water, return 0: */
	if(IsOnWater())return 0;

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

	/*Retrieve all inputs we will be needing: */
	Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum);   _assert_(weights_input);
	Input* vx_input     =inputs->GetInput(VxEnum);        _assert_(vx_input);
	Input* vy_input     =inputs->GetInput(VyEnum);        _assert_(vy_input);
	Input* vxobs_input  =inputs->GetInput(InversionVxObsEnum);     _assert_(vxobs_input);
	Input* vyobs_input  =inputs->GetInput(InversionVyObsEnum);     _assert_(vyobs_input);
	
	/* Start  looping on the number of gaussian points: */
	gauss=new GaussTria(4);
	for (ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		/* Get Jacobian determinant: */
		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);

		/*Get all parameters at gaussian point*/
		weights_input->GetInputValue(&weight,gauss,weight_index);
		vx_input->GetInputValue(&vx,gauss);
		vy_input->GetInputValue(&vy,gauss);
		vxobs_input->GetInputValue(&vxobs,gauss);
		vyobs_input->GetInputValue(&vyobs,gauss);

		/*Compute SurfaceRelVelMisfit:
		 *
		 *      1            [        |u| + eps     2          |v| + eps     2  ]
		 * J = --- \bar{v}^2 | log ( -----------  )   +  log ( -----------  )   |  
		 *      2            [       |u    |+ eps              |v    |+ eps     ]
		 *                              obs                       obs
		 */
		misfit=0.5*pow(meanvel,2.)*(
					pow(log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)),2.) +
					pow(log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)),2.) );

		if(process_units)UnitConversion(misfit,IuToExtEnum,SurfaceLogVxVyMisfitEnum);

		/*Add to cost function*/
		Jelem+=misfit*weight*Jdet*gauss->weight;
	}

	/*clean-up and Return: */
	delete gauss;
	return Jelem;
}
/*}}}*/
/*FUNCTION Tria::SurfaceAbsVelMisfit {{{1*/
double Tria::SurfaceAbsVelMisfit(bool process_units,int weight_index){

	const int    numdof=NDOF2*NUMVERTICES;

	int        i,ig;
	double     Jelem=0;
	double     misfit,Jdet;
	double     vx,vy,vxobs,vyobs,weight;
	double     xyz_list[NUMVERTICES][3];
	GaussTria *gauss=NULL;

	/*If on water, return 0: */
	if(IsOnWater())return 0;

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

	/*Retrieve all inputs we will be needing: */
	Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum);   _assert_(weights_input);
	Input* vx_input     =inputs->GetInput(VxEnum);        _assert_(vx_input);
	Input* vy_input     =inputs->GetInput(VyEnum);        _assert_(vy_input);
	Input* vxobs_input  =inputs->GetInput(InversionVxObsEnum);     _assert_(vxobs_input);
	Input* vyobs_input  =inputs->GetInput(InversionVyObsEnum);     _assert_(vyobs_input);

	/* Start  looping on the number of gaussian points: */
	gauss=new GaussTria(2);
	for (ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		/* Get Jacobian determinant: */
		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);

		/*Get all parameters at gaussian point*/
		weights_input->GetInputValue(&weight,gauss,weight_index);
		vx_input->GetInputValue(&vx,gauss);
		vy_input->GetInputValue(&vy,gauss);
		vxobs_input->GetInputValue(&vxobs,gauss);
		vyobs_input->GetInputValue(&vyobs,gauss);

		/*Compute SurfaceAbsVelMisfitEnum:
		 *
		 *      1  [           2              2 ]
		 * J = --- | (u - u   )  +  (v - v   )  |
		 *      2  [       obs            obs   ]
		 *
		 */
		misfit=0.5*( pow(vx-vxobs,2.) + pow(vy-vyobs,2.) );

		if(process_units)UnitConversion(misfit,IuToExtEnum,SurfaceAverageVelMisfitEnum);

		/*Add to cost function*/
		Jelem+=misfit*weight*Jdet*gauss->weight;
	}

	/*clean up and Return: */
	delete gauss;
	return Jelem;
}
/*}}}*/
/*FUNCTION Tria::SurfaceRelVelMisfit {{{1*/
double Tria::SurfaceRelVelMisfit(bool process_units,int weight_index){
	const int  numdof=2*NUMVERTICES;

	int        i,ig;
	double     Jelem=0;
	double     scalex=1,scaley=1;
	double     misfit,Jdet;
	double     epsvel=2.220446049250313e-16;
	double     meanvel=3.170979198376458e-05; /*1000 m/yr*/
	double     vx,vy,vxobs,vyobs,weight;
	double     xyz_list[NUMVERTICES][3];
	GaussTria *gauss=NULL;

	/*If on water, return 0: */
	if(IsOnWater())return 0;

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

	/*Retrieve all inputs we will be needing: */
	Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum);   _assert_(weights_input);
	Input* vx_input     =inputs->GetInput(VxEnum);        _assert_(vx_input);
	Input* vy_input     =inputs->GetInput(VyEnum);        _assert_(vy_input);
	Input* vxobs_input  =inputs->GetInput(InversionVxObsEnum);     _assert_(vxobs_input);
	Input* vyobs_input  =inputs->GetInput(InversionVyObsEnum);     _assert_(vyobs_input);

	/* Start  looping on the number of gaussian points: */
	gauss=new GaussTria(4);
	for (ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		/* Get Jacobian determinant: */
		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);

		/*Get all parameters at gaussian point*/
		weights_input->GetInputValue(&weight,gauss,weight_index);
		vx_input->GetInputValue(&vx,gauss);
		vy_input->GetInputValue(&vy,gauss);
		vxobs_input->GetInputValue(&vxobs,gauss);
		vyobs_input->GetInputValue(&vyobs,gauss);

		/*Compute SurfaceRelVelMisfit:
		 *                        
		 *      1  [     \bar{v}^2             2   \bar{v}^2              2 ]
		 * J = --- | -------------  (u - u   ) + -------------  (v - v   )  |
		 *      2  [  (u   + eps)^2       obs    (v   + eps)^2       obs    ]
		 *              obs                        obs                      
		 */
		scalex=pow(meanvel/(vxobs+epsvel),2.); if(vxobs==0)scalex=0;
		scaley=pow(meanvel/(vyobs+epsvel),2.); if(vyobs==0)scaley=0;
		misfit=0.5*(scalex*pow((vx-vxobs),2.)+scaley*pow((vy-vyobs),2.));
		if(process_units)UnitConversion(misfit,IuToExtEnum,SurfaceRelVelMisfitEnum);

		/*Add to cost function*/
		Jelem+=misfit*weight*Jdet*gauss->weight;
	}

	/*clean up and Return: */
	delete gauss;
	return Jelem;
}
/*}}}*/
/*FUNCTION Tria::ThicknessAbsGradient{{{1*/
double Tria::ThicknessAbsGradient(bool process_units,int weight_index){

	/* Intermediaries */
	int        ig;
	double     Jelem = 0;
	double     weight;
	double     Jdet;
	double     xyz_list[NUMVERTICES][3];
	double     dp[NDOF2];
	GaussTria *gauss = NULL;

	/*retrieve parameters and inputs*/

	/*If on water, return 0: */
	if(IsOnWater()) return 0;

	/*Retrieve all inputs we will be needing: */
	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
	Input* weights_input  =inputs->GetInput(InversionCostFunctionsCoefficientsEnum);   _assert_(weights_input);
	Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);

	/* Start looping on the number of gaussian points: */
	gauss=new GaussTria(2);
	for (ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		/* Get Jacobian determinant: */
		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);

		/*Get all parameters at gaussian point*/
		weights_input->GetInputValue(&weight,gauss,weight_index);
		thickness_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss);

		/*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ 
		Jelem+=weight*1/2*(pow(dp[0],2.)+pow(dp[1],2.))*Jdet*gauss->weight;
	}

	/*Clean up and return*/
	delete gauss;
	return Jelem;
}
/*}}}*/
/*FUNCTION Tria::ThicknessAbsMisfit {{{1*/
double Tria::ThicknessAbsMisfit(bool process_units,int weight_index){

	/*Intermediaries*/
	int        i,ig;
	double     thickness,thicknessobs,weight;
	double     Jdet;
	double     Jelem = 0;
	double     xyz_list[NUMVERTICES][3];
	GaussTria *gauss = NULL;
	double     dH[2];

	/*If on water, return 0: */
	if(IsOnWater())return 0;

	/*Retrieve all inputs we will be needing: */
	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
	Input* thickness_input   =inputs->GetInput(ThicknessEnum);   _assert_(thickness_input);
	Input* thicknessobs_input=inputs->GetInput(InversionThicknessObsEnum);_assert_(thicknessobs_input);
	Input* weights_input     =inputs->GetInput(InversionCostFunctionsCoefficientsEnum);     _assert_(weights_input);

	/* Start  looping on the number of gaussian points: */
	gauss=new GaussTria(2);
	for (ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		/* Get Jacobian determinant: */
		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);

		/*Get parameters at gauss point*/
		thickness_input->GetInputValue(&thickness,gauss);
		thickness_input->GetInputDerivativeValue(&dH[0],&xyz_list[0][0],gauss);
		thicknessobs_input->GetInputValue(&thicknessobs,gauss);
		weights_input->GetInputValue(&weight,gauss,weight_index);

		/*compute ThicknessAbsMisfit*/
		Jelem+=0.5*pow(thickness-thicknessobs,2.0)*weight*Jdet*gauss->weight;
	}

	/* clean up and Return: */
	delete gauss;
	return Jelem;
}
/*}}}*/
/*FUNCTION Tria::CreatePVectorAdjointBalancethickness{{{1*/
ElementVector* Tria::CreatePVectorAdjointBalancethickness(void){

	/*Constants*/
	const int    numdof=1*NUMVERTICES;

	/*Intermediaries */
	int         i,ig,resp;
	double      Jdet;
	double      thickness,thicknessobs,weight;
	int        *responses = NULL;
	int         num_responses;
	double      xyz_list[NUMVERTICES][3];
	double      basis[3];
	double      dbasis[NDOF2][NUMVERTICES];
	double      dH[2];
	GaussTria*  gauss=NULL;

	/*Initialize Element vector*/
	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);

	/*Retrieve all inputs and parameters*/
	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
	this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
	this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum);
	Input* thickness_input    = inputs->GetInput(ThicknessEnum);   _assert_(thickness_input);
	Input* thicknessobs_input = inputs->GetInput(InversionThicknessObsEnum);_assert_(thicknessobs_input);
	Input* weights_input      = inputs->GetInput(InversionCostFunctionsCoefficientsEnum);     _assert_(weights_input);

	/* Start  looping on the number of gaussian points: */
	gauss=new GaussTria(2);
	for(ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
		GetNodalFunctions(basis, gauss);
		GetNodalFunctionsDerivatives(&dbasis[0][0],&xyz_list[0][0],gauss);

		thickness_input->GetInputValue(&thickness, gauss);
		thickness_input->GetInputDerivativeValue(&dH[0],&xyz_list[0][0],gauss);
		thicknessobs_input->GetInputValue(&thicknessobs, gauss);

		/*Loop over all requested responses*/
		for(resp=0;resp<num_responses;resp++) switch(responses[resp]){

			case ThicknessAbsMisfitEnum:
				weights_input->GetInputValue(&weight, gauss,resp);
				for(i=0;i<numdof;i++) pe->values[i]+=(thicknessobs-thickness)*weight*Jdet*gauss->weight*basis[i];
				break;
			case ThicknessAbsGradientEnum:
				weights_input->GetInputValue(&weight, gauss,resp);
				for(i=0;i<numdof;i++) pe->values[i]+= - weight*dH[0]*dbasis[0][i]*Jdet*gauss->weight;
				for(i=0;i<numdof;i++) pe->values[i]+= - weight*dH[1]*dbasis[1][i]*Jdet*gauss->weight;
				break;
			default:
				_error_("response %s not supported yet",EnumToStringx(responses[resp]));
		}
	}

	/*Clean up and return*/
	delete gauss;
	xfree((void**)&responses);
	return pe;
}
/*}}}*/
/*FUNCTION Tria::CreatePVectorAdjointHoriz{{{1*/
ElementVector* Tria::CreatePVectorAdjointHoriz(void){

	/*Constants*/
	const int    numdof=NDOF2*NUMVERTICES;

	/*Intermediaries */
	int        i,resp,ig;
	int       *responses=NULL;
	int        num_responses;
	double     Jdet;
	double     obs_velocity_mag,velocity_mag;
	double     dux,duy;
	double     epsvel=2.220446049250313e-16;
	double     meanvel=3.170979198376458e-05; /*1000 m/yr*/
	double     scalex=0,scaley=0,scale=0,S=0;
	double     vx,vy,vxobs,vyobs,weight;
	double     xyz_list[NUMVERTICES][3];
	double     basis[3];
	GaussTria* gauss=NULL;

	/*Initialize Element vector*/
	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);

	/*Retrieve all inputs and parameters*/
	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
	this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
	this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum);
	Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum);   _assert_(weights_input);
	Input* vx_input     =inputs->GetInput(VxEnum);        _assert_(vx_input);
	Input* vy_input     =inputs->GetInput(VyEnum);        _assert_(vy_input);
	Input* vxobs_input  =inputs->GetInput(InversionVxObsEnum);     _assert_(vxobs_input);
	Input* vyobs_input  =inputs->GetInput(InversionVyObsEnum);     _assert_(vyobs_input);

	/*Get Surface if required by one response*/
	for(resp=0;resp<num_responses;resp++){
		if(responses[resp]==SurfaceAverageVelMisfitEnum){
			inputs->GetInputValue(&S,SurfaceAreaEnum); break;
		}
	}

	/* Start  looping on the number of gaussian points: */
	gauss=new GaussTria(4);
	for (ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		/* Get Jacobian determinant: */
		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);

		/*Get all parameters at gaussian point*/
		vx_input->GetInputValue(&vx,gauss);
		vy_input->GetInputValue(&vy,gauss);
		vxobs_input->GetInputValue(&vxobs,gauss);
		vyobs_input->GetInputValue(&vyobs,gauss);
		GetNodalFunctions(basis, gauss);

		/*Loop over all requested responses*/
		for(resp=0;resp<num_responses;resp++){

			weights_input->GetInputValue(&weight,gauss,resp);

			switch(responses[resp]){
				case SurfaceAbsVelMisfitEnum:
					/*
					 *      1  [           2              2 ]
					 * J = --- | (u - u   )  +  (v - v   )  |
					 *      2  [       obs            obs   ]
					 *
					 *        dJ
					 * DU = - -- = (u   - u )
					 *        du     obs
					 */
					for (i=0;i<NUMVERTICES;i++){
						dux=vxobs-vx;
						duy=vyobs-vy;
						pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
						pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
					}
					break;
				case SurfaceRelVelMisfitEnum:
					/*
					 *      1  [     \bar{v}^2             2   \bar{v}^2              2 ]
					 * J = --- | -------------  (u - u   ) + -------------  (v - v   )  |
					 *      2  [  (u   + eps)^2       obs    (v   + eps)^2       obs    ]
					 *              obs                        obs                      
					 *
					 *        dJ     \bar{v}^2
					 * DU = - -- = ------------- (u   - u )
					 *        du   (u   + eps)^2    obs
					 *               obs
					 */
					for (i=0;i<NUMVERTICES;i++){
						scalex=pow(meanvel/(vxobs+epsvel),2.); if(vxobs==0)scalex=0;
						scaley=pow(meanvel/(vyobs+epsvel),2.); if(vyobs==0)scaley=0;
						dux=scalex*(vxobs-vx);
						duy=scaley*(vyobs-vy);
						pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
						pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
					}
					break;
				case SurfaceLogVelMisfitEnum:
					/*
					 *                 [        vel + eps     ] 2
					 * J = 4 \bar{v}^2 | log ( -----------  ) |  
					 *                 [       vel   + eps    ]
					 *                            obs
					 *
					 *        dJ                 2 * log(...)
					 * DU = - -- = - 4 \bar{v}^2 -------------  u
					 *        du                 vel^2 + eps
					 *            
					 */
					for (i=0;i<NUMVERTICES;i++){
						velocity_mag    =sqrt(pow(vx,   2.)+pow(vy,   2.))+epsvel;
						obs_velocity_mag=sqrt(pow(vxobs,2.)+pow(vyobs,2.))+epsvel;
						scale=-8*pow(meanvel,2.)/pow(velocity_mag,2.)*log(velocity_mag/obs_velocity_mag);
						dux=scale*vx;
						duy=scale*vy;
						pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
						pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
					}
					break;
				case SurfaceAverageVelMisfitEnum:
					/*
					 *      1                    2              2
					 * J = ---  sqrt(  (u - u   )  +  (v - v   )  )
					 *      S                obs            obs
					 *
					 *        dJ      1       1 
					 * DU = - -- = - --- ----------- * 2 (u - u   )
					 *        du      S  2 sqrt(...)           obs
					 */
					for (i=0;i<NUMVERTICES;i++){
						scale=1./(S*2*sqrt(pow(vx-vxobs,2.)+pow(vy-vyobs,2.))+epsvel);
						dux=scale*(vxobs-vx);
						duy=scale*(vyobs-vy);
						pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
						pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
					}
					break;
				case SurfaceLogVxVyMisfitEnum:
					/*
					 *      1            [        |u| + eps     2          |v| + eps     2  ]
					 * J = --- \bar{v}^2 | log ( -----------  )   +  log ( -----------  )   |  
					 *      2            [       |u    |+ eps              |v    |+ eps     ]
					 *                              obs                       obs
					 *        dJ                              1      u                             1
					 * DU = - -- = - \bar{v}^2 log(u...) --------- ----  ~ - \bar{v}^2 log(u...) ------
					 *        du                         |u| + eps  |u|                           u + eps
					 */
					for (i=0;i<NUMVERTICES;i++){
						dux = - pow(meanvel,2.) * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
						duy = - pow(meanvel,2.) * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
						pe->values[i*NDOF2+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
						pe->values[i*NDOF2+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
					}
					break;
				case DragCoefficientAbsGradientEnum:
					/*Nothing in P vector*/
					break;
				case ThicknessAbsGradientEnum:
					/*Nothing in P vector*/
					break;
				case RheologyBbarAbsGradientEnum:
					/*Nothing in P vector*/
					break;
				default:
					_error_("response %s not supported yet",EnumToStringx(responses[resp]));
			}
		}
	}

	/*Clean up and return*/
	delete gauss;
	xfree((void**)&responses);
	return pe;
}
/*}}}*/
/*FUNCTION Tria::CreatePVectorAdjointStokes{{{1*/
ElementVector* Tria::CreatePVectorAdjointStokes(void){

	/*Intermediaries */
	int        i,resp,ig;
	int       *responses=NULL;
	int        num_responses;
	double     Jdet;
	double     obs_velocity_mag,velocity_mag;
	double     dux,duy;
	double     epsvel=2.220446049250313e-16;
	double     meanvel=3.170979198376458e-05; /*1000 m/yr*/
	double     scalex=0,scaley=0,scale=0,S=0;
	double     vx,vy,vxobs,vyobs,weight;
	double     xyz_list[NUMVERTICES][3];
	double     basis[3];
	GaussTria* gauss=NULL;

	/*Initialize Element vector*/
	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,StokesApproximationEnum);

	/*Retrieve all inputs and parameters*/
	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
	this->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
	this->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum);
	Input* weights_input = inputs->GetInput(InversionCostFunctionsCoefficientsEnum);   _assert_(weights_input);
	Input* vx_input      = inputs->GetInput(VxEnum);        _assert_(vx_input);
	Input* vy_input      = inputs->GetInput(VyEnum);        _assert_(vy_input);
	Input* vxobs_input   = inputs->GetInput(InversionVxObsEnum);     _assert_(vxobs_input);
	Input* vyobs_input   = inputs->GetInput(InversionVyObsEnum);     _assert_(vyobs_input);

	/*Get Surface if required by one response*/
	for(resp=0;resp<num_responses;resp++){
		if(responses[resp]==SurfaceAverageVelMisfitEnum){
			inputs->GetInputValue(&S,SurfaceAreaEnum); break;
		}
	}

	/* Start  looping on the number of gaussian points: */
	gauss=new GaussTria(4);
	for (ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		/* Get Jacobian determinant: */
		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);

		/*Get all parameters at gaussian point*/
		vx_input->GetInputValue(&vx,gauss);
		vy_input->GetInputValue(&vy,gauss);
		vxobs_input->GetInputValue(&vxobs,gauss);
		vyobs_input->GetInputValue(&vyobs,gauss);
		GetNodalFunctions(basis, gauss);

		/*Loop over all requested responses*/
		for(resp=0;resp<num_responses;resp++){

			weights_input->GetInputValue(&weight,gauss,resp);

			switch(responses[resp]){

				case SurfaceAbsVelMisfitEnum:
					/*
					 *      1  [           2              2 ]
					 * J = --- | (u - u   )  +  (v - v   )  |
					 *      2  [       obs            obs   ]
					 *
					 *        dJ
					 * DU = - -- = (u   - u )
					 *        du     obs
					 */
					for (i=0;i<NUMVERTICES;i++){
						dux=vxobs-vx;
						duy=vyobs-vy;
						pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
						pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
					}
					break;
				case SurfaceRelVelMisfitEnum:
					/*
					 *      1  [     \bar{v}^2             2   \bar{v}^2              2 ]
					 * J = --- | -------------  (u - u   ) + -------------  (v - v   )  |
					 *      2  [  (u   + eps)^2       obs    (v   + eps)^2       obs    ]
					 *              obs                        obs                      
					 *
					 *        dJ     \bar{v}^2
					 * DU = - -- = ------------- (u   - u )
					 *        du   (u   + eps)^2    obs
					 *               obs
					 */
					for (i=0;i<NUMVERTICES;i++){
						scalex=pow(meanvel/(vxobs+epsvel),2.); if(vxobs==0)scalex=0;
						scaley=pow(meanvel/(vyobs+epsvel),2.); if(vyobs==0)scaley=0;
						dux=scalex*(vxobs-vx);
						duy=scaley*(vyobs-vy);
						pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
						pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
					}
					break;
				case SurfaceLogVelMisfitEnum:
					/*
					 *                 [        vel + eps     ] 2
					 * J = 4 \bar{v}^2 | log ( -----------  ) |  
					 *                 [       vel   + eps    ]
					 *                            obs
					 *
					 *        dJ                 2 * log(...)
					 * DU = - -- = - 4 \bar{v}^2 -------------  u
					 *        du                 vel^2 + eps
					 *            
					 */
					for (i=0;i<NUMVERTICES;i++){
						velocity_mag    =sqrt(pow(vx,   2.)+pow(vy,   2.))+epsvel;
						obs_velocity_mag=sqrt(pow(vxobs,2.)+pow(vyobs,2.))+epsvel;
						scale=-8*pow(meanvel,2.)/pow(velocity_mag,2.)*log(velocity_mag/obs_velocity_mag);
						dux=scale*vx;
						duy=scale*vy;
						pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
						pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
					}
					break;
				case SurfaceAverageVelMisfitEnum:
					/*
					 *      1                    2              2
					 * J = ---  sqrt(  (u - u   )  +  (v - v   )  )
					 *      S                obs            obs
					 *
					 *        dJ      1       1 
					 * DU = - -- = - --- ----------- * 2 (u - u   )
					 *        du      S  2 sqrt(...)           obs
					 */
					for (i=0;i<NUMVERTICES;i++){
						scale=1./(S*2*sqrt(pow(vx-vxobs,2.)+pow(vy-vyobs,2.))+epsvel);
						dux=scale*(vxobs-vx);
						duy=scale*(vyobs-vy);
						pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
						pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
					}
					break;
				case SurfaceLogVxVyMisfitEnum:
					/*
					 *      1            [        |u| + eps     2          |v| + eps     2  ]
					 * J = --- \bar{v}^2 | log ( -----------  )   +  log ( -----------  )   |  
					 *      2            [       |u    |+ eps              |v    |+ eps     ]
					 *                              obs                       obs
					 *        dJ                              1      u                             1
					 * DU = - -- = - \bar{v}^2 log(u...) --------- ----  ~ - \bar{v}^2 log(u...) ------
					 *        du                         |u| + eps  |u|                           u + eps
					 */
					for (i=0;i<NUMVERTICES;i++){
						dux = - pow(meanvel,2.) * log((fabs(vx)+epsvel)/(fabs(vxobs)+epsvel)) / (vx+epsvel);
						duy = - pow(meanvel,2.) * log((fabs(vy)+epsvel)/(fabs(vyobs)+epsvel)) / (vy+epsvel);
						pe->values[i*NDOF4+0]+=dux*weight*Jdet*gauss->weight*basis[i]; 
						pe->values[i*NDOF4+1]+=duy*weight*Jdet*gauss->weight*basis[i]; 
					}
					break;
				case DragCoefficientAbsGradientEnum:
					/*Nothing in P vector*/
					break;
				case ThicknessAbsGradientEnum:
					/*Nothing in P vector*/
					break;
				case RheologyBbarAbsGradientEnum:
					/*Nothing in P vector*/
					break;
				default:
					_error_("response %s not supported yet",EnumToStringx(responses[resp]));
			}
		}
	}

	/*Clean up and return*/
	delete gauss;
	xfree((void**)&responses);
	return pe;
}
/*}}}*/
/*FUNCTION Tria::DragCoefficientAbsGradient{{{1*/
double Tria::DragCoefficientAbsGradient(bool process_units,int weight_index){

	/* Intermediaries */
	int        ig;
	double     Jelem = 0;
	double     weight;
	double     Jdet;
	double     xyz_list[NUMVERTICES][3];
	double     dp[NDOF2];
	GaussTria *gauss = NULL;

	/*retrieve parameters and inputs*/

	/*If on water, return 0: */
	if(IsOnWater()) return 0;

	/*Retrieve all inputs we will be needing: */
	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
	Input* weights_input=inputs->GetInput(InversionCostFunctionsCoefficientsEnum);         _assert_(weights_input);
	Input* drag_input   =inputs->GetInput(FrictionCoefficientEnum); _assert_(drag_input);

	/* Start looping on the number of gaussian points: */
	gauss=new GaussTria(2);
	for (ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		/* Get Jacobian determinant: */
		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);

		/*Get all parameters at gaussian point*/
		weights_input->GetInputValue(&weight,gauss,weight_index);
		drag_input->GetInputDerivativeValue(&dp[0],&xyz_list[0][0],gauss);

		/*Tikhonov regularization: J = 1/2 ((dp/dx)^2 + (dp/dy)^2) */ 
		Jelem+=weight*1/2*(pow(dp[0],2.)+pow(dp[1],2.))*Jdet*gauss->weight;
	}

	/*Clean up and return*/
	delete gauss;
	return Jelem;
}
/*}}}*/
/*FUNCTION Tria::CreateKMatrixAdjointBalancethickness {{{1*/
ElementMatrix* Tria::CreateKMatrixAdjointBalancethickness(void){

	ElementMatrix* Ke=NULL;

	/*Get Element Matrix of the forward model*/
	switch(GetElementType()){
		case P1Enum:
			Ke=CreateKMatrixBalancethickness_CG();
			break;
		case P1DGEnum:
			Ke=CreateKMatrixBalancethickness_DG();
			break;
		default:
			_error_("Element type %s not supported yet",EnumToStringx(GetElementType()));
	}

	/*Transpose and return Ke*/
	Ke->Transpose();
	return Ke;
}
/*}}}*/
/*FUNCTION Tria::CreateKMatrixAdjointMacAyeal{{{1*/
ElementMatrix* Tria::CreateKMatrixAdjointMacAyeal(void){

	/*Constants*/
	const int    numdof=NDOF2*NUMVERTICES;

	/*Intermediaries */
	int        i,j,ig;
	bool       incomplete_adjoint;
	double     xyz_list[NUMVERTICES][3];
	double     Jdet,thickness;
	double     eps1dotdphii,eps1dotdphij;
	double     eps2dotdphii,eps2dotdphij;
	double     mu_prime;
	double     epsilon[3];/* epsilon=[exx,eyy,exy];*/
	double     eps1[2],eps2[2];
	double     phi[NUMVERTICES];
	double     dphi[2][NUMVERTICES];
	GaussTria *gauss=NULL;

	/*Initialize Jacobian with regular MacAyeal (first part of the Gateau derivative)*/
	parameters->FindParam(&incomplete_adjoint,InversionIncompleteAdjointEnum);
	ElementMatrix* Ke=CreateKMatrixDiagnosticMacAyeal();
	if(incomplete_adjoint) return Ke;

	/*Retrieve all inputs and parameters*/
	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
	Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
	Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
	Input* thickness_input=inputs->GetInput(ThicknessEnum); _assert_(thickness_input);

	/* Start  looping on the number of gaussian points: */
	gauss=new GaussTria(2);
	for (ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		GetJacobianDeterminant2d(&Jdet, &xyz_list[0][0],gauss);
		GetNodalFunctionsDerivatives(&dphi[0][0],&xyz_list[0][0],gauss);

		thickness_input->GetInputValue(&thickness, gauss);
		this->GetStrainRate2d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
		matice->GetViscosity2dDerivativeEpsSquare(&mu_prime,&epsilon[0]);
		eps1[0]=2*epsilon[0]+epsilon[1];   eps2[0]=epsilon[2];
		eps1[1]=epsilon[2];                eps2[1]=epsilon[0]+2*epsilon[1];

		for(i=0;i<3;i++){
			for(j=0;j<3;j++){
				eps1dotdphii=eps1[0]*dphi[0][i]+eps1[1]*dphi[1][i];
				eps1dotdphij=eps1[0]*dphi[0][j]+eps1[1]*dphi[1][j];
				eps2dotdphii=eps2[0]*dphi[0][i]+eps2[1]*dphi[1][i];
				eps2dotdphij=eps2[0]*dphi[0][j]+eps2[1]*dphi[1][j];

				Ke->values[6*(2*i+0)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps1dotdphii;
				Ke->values[6*(2*i+0)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps1dotdphii;
				Ke->values[6*(2*i+1)+2*j+0]+=gauss->weight*Jdet*2*mu_prime*thickness*eps1dotdphij*eps2dotdphii;
				Ke->values[6*(2*i+1)+2*j+1]+=gauss->weight*Jdet*2*mu_prime*thickness*eps2dotdphij*eps2dotdphii;
			}
		}
	}

	/*Transform Coordinate System*/
	TransformStiffnessMatrixCoord(Ke,nodes,NUMVERTICES,XYEnum);

	/*Clean up and return*/
	delete gauss;
	//Ke->Transpose();
	return Ke;
}
/*}}}*/
/*FUNCTION Tria::InputUpdateFromSolutionAdjointHoriz {{{1*/
void  Tria::InputUpdateFromSolutionAdjointHoriz(double* solution){

	const int numdof=NDOF2*NUMVERTICES;

	int       i;
	int*      doflist=NULL;
	double    values[numdof];
	double    lambdax[NUMVERTICES];
	double    lambday[NUMVERTICES];

	/*Get dof list: */
	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);

	/*Use the dof list to index into the solution vector: */
	for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];

	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
	for(i=0;i<NUMVERTICES;i++){
		lambdax[i]=values[i*NDOF2+0];
		lambday[i]=values[i*NDOF2+1];

		/*Check solution*/
		if(isnan(lambdax[i])) _error_("NaN found in solution vector");
		if(isnan(lambday[i])) _error_("NaN found in solution vector");
	}

	/*Add vx and vy as inputs to the tria element: */
	this->inputs->AddInput(new TriaP1Input(AdjointxEnum,lambdax));
	this->inputs->AddInput(new TriaP1Input(AdjointyEnum,lambday));

	/*Free ressources:*/
	xfree((void**)&doflist);
}
/*}}}*/
/*FUNCTION Tria::InputUpdateFromSolutionAdjointBalancethickness {{{1*/
void  Tria::InputUpdateFromSolutionAdjointBalancethickness(double* solution){

	const int numdof=NDOF1*NUMVERTICES;

	int       i;
	int*      doflist=NULL;
	double    values[numdof];
	double    lambda[NUMVERTICES];

	/*Get dof list: */
	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);

	/*Use the dof list to index into the solution vector: */
	for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];

	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
	for(i=0;i<numdof;i++){
		lambda[i]=values[i];
		if(isnan(lambda[i])) _error_("NaN found in solution vector");
	}

	/*Add vx and vy as inputs to the tria element: */
	this->inputs->AddInput(new TriaP1Input(AdjointEnum,lambda));

	/*Free ressources:*/
	xfree((void**)&doflist);
}
/*}}}*/
/*FUNCTION Tria::GetVectorFromControlInputs{{{1*/
void  Tria::GetVectorFromControlInputs(Vec vector,int control_enum,int control_index,const char* data){

	int doflist1[NUMVERTICES];
	Input *input=NULL;

	/*Get out if this is not an element input*/
	if(!IsInput(control_enum)) return;

	/*Prepare index list*/
	GradientIndexing(&doflist1[0],control_index);

	/*Get input (either in element or material)*/
	if(control_enum==MaterialsRheologyBbarEnum){
		input=(Input*)matice->inputs->GetInput(control_enum); _assert_(input);
	}
	else{
		input=(Input*)this->inputs->GetInput(control_enum);   _assert_(input);
	}

	/*Check that it is a ControlInput*/
	if (input->ObjectEnum()!=ControlInputEnum){
		_error_("input %s is not a ControlInput",EnumToStringx(control_enum));
	}

	((ControlInput*)input)->GetVectorFromInputs(vector,&doflist1[0],data);
}
/*}}}*/
/*FUNCTION Tria::SetControlInputsFromVector{{{1*/
void  Tria::SetControlInputsFromVector(double* vector,int control_enum,int control_index){

	double  values[NUMVERTICES];
	int     doflist1[NUMVERTICES];
	Input  *input     = NULL;
	Input  *new_input = NULL;

	/*Get out if this is not an element input*/
	if(!IsInput(control_enum)) return;

	/*Prepare index list*/
	GradientIndexing(&doflist1[0],control_index);

	/*Get values on vertices*/
	for (int i=0;i<NUMVERTICES;i++){
		values[i]=vector[doflist1[i]];
	}
	new_input = new TriaP1Input(control_enum,values);

	if(control_enum==MaterialsRheologyBbarEnum){
		input=(Input*)matice->inputs->GetInput(control_enum); _assert_(input);
	}
	else{
		input=(Input*)this->inputs->GetInput(control_enum);   _assert_(input);
	}

	if (input->ObjectEnum()!=ControlInputEnum){
		_error_("input %s is not a ControlInput",EnumToStringx(control_enum));
	}

	((ControlInput*)input)->SetInput(new_input);
}
/*}}}*/
#endif

#ifdef _HAVE_HYDROLOGY_
/*FUNCTION Tria::CreateHydrologyWaterVelocityInput {{{1*/
void Tria::CreateHydrologyWaterVelocityInput(void){

	/*material parameters: */
	double mu_water;
	double VelocityFactor;  // This factor represents the number 12 in laminar flow velocity which can vary by differnt hydrology.CR
	double n_man,CR;
	double w;
	double rho_ice, rho_water, g;
	double dsdx,dsdy,dbdx,dbdy;
	double vx[NUMVERTICES];
	double vy[NUMVERTICES];
	GaussTria *gauss = NULL;

	/*Retrieve all inputs and parameters*/
	rho_ice=matpar->GetRhoIce();
	rho_water=matpar->GetRhoWater();
	g=matpar->GetG();
	CR=matpar->GetHydrologyCR(); // To have Lebrocq equavalent equation: CR=0.01,n_man=0.02
	n_man=matpar->GetHydrologyN(); 
	mu_water=matpar->GetMuWater();
	Input* surfaceslopex_input=inputs->GetInput(SurfaceSlopeXEnum); _assert_(surfaceslopex_input);
	Input* surfaceslopey_input=inputs->GetInput(SurfaceSlopeYEnum); _assert_(surfaceslopey_input);
	Input* bedslopex_input=inputs->GetInput(BedSlopeXEnum);         _assert_(bedslopex_input);
	Input* bedslopey_input=inputs->GetInput(BedSlopeYEnum);         _assert_(bedslopey_input);
	Input* watercolumn_input=inputs->GetInput(WatercolumnEnum);     _assert_(watercolumn_input);

	/* compute VelocityFactor */
	VelocityFactor= n_man*pow(CR,2)*rho_water*g/mu_water;
	
	gauss=new GaussTria();
	for (int iv=0;iv<NUMVERTICES;iv++){
		gauss->GaussVertex(iv);
		surfaceslopex_input->GetInputValue(&dsdx,gauss);
		surfaceslopey_input->GetInputValue(&dsdy,gauss);
		bedslopex_input->GetInputValue(&dbdx,gauss);
		bedslopey_input->GetInputValue(&dbdy,gauss);
		watercolumn_input->GetInputValue(&w,gauss);

		/* Water velocity x and y components */
	//	vx[iv]= - pow(w,2)/(12 * mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx);
	//	vy[iv]= - pow(w,2)/(12 * mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy);
	
		vx[iv]= - pow(w,2)/(VelocityFactor* mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx);
		vy[iv]= - pow(w,2)/(VelocityFactor* mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy);
	}

	/*clean-up*/
	delete gauss;

	/*Add to inputs*/
	this->inputs->AddInput(new TriaP1Input(HydrologyWaterVxEnum,vx));
	this->inputs->AddInput(new TriaP1Input(HydrologyWaterVyEnum,vy));
}
/*}}}*/
/*FUNCTION Tria::CreateKMatrixHydrology{{{1*/
ElementMatrix* Tria::CreateKMatrixHydrology(void){

	/*Constants*/
	const int    numdof=NDOF1*NUMVERTICES;

	/*Intermediaries */
	double     diffusivity;
	int        i,j,ig;
	double     Jdettria,DL_scalar,dt,h;
	double     vx,vy,vel,dvxdx,dvydy;
	double     dvx[2],dvy[2];
	double     v_gauss[2]={0.0};
	double     xyz_list[NUMVERTICES][3];
	double     L[NUMVERTICES];
	double     B[2][NUMVERTICES];
	double     Bprime[2][NUMVERTICES];
	double     K[2][2]                        ={0.0};
	double     KDL[2][2]                      ={0.0};
	double     DL[2][2]                        ={0.0};
	double     DLprime[2][2]                   ={0.0};
	GaussTria *gauss=NULL;

	/*Skip if water or ice shelf element*/
	if(IsOnWater() | IsFloating()) return NULL;

	/*Initialize Element matrix*/
	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);

	/*Create water velocity vx and vy from current inputs*/
	CreateHydrologyWaterVelocityInput();

	/*Retrieve all inputs and parameters*/
	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
	this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
	this->parameters->FindParam(&diffusivity,HydrologyStabilizationEnum);
	Input* vx_input=inputs->GetInput(HydrologyWaterVxEnum); _assert_(vx_input);
	Input* vy_input=inputs->GetInput(HydrologyWaterVyEnum); _assert_(vy_input);
	h=sqrt(2*this->GetArea());

	/* Start  looping on the number of gaussian points: */
	gauss=new GaussTria(2);
	for (ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
		GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);

		vx_input->GetInputValue(&vx,gauss);
		vy_input->GetInputValue(&vy,gauss);
		vx_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
		vy_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);

		DL_scalar=gauss->weight*Jdettria;

		TripleMultiply( &L[0],1,numdof,1,
					&DL_scalar,1,1,0,
					&L[0],1,numdof,0,
					&Ke->values[0],1);

		GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss);
		GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);

		dvxdx=dvx[0];
		dvydy=dvy[1];
		DL_scalar=dt*gauss->weight*Jdettria;

		DL[0][0]=DL_scalar*dvxdx;
		DL[1][1]=DL_scalar*dvydy;
		DLprime[0][0]=DL_scalar*vx;
		DLprime[1][1]=DL_scalar*vy;

		TripleMultiply( &B[0][0],2,numdof,1,
					&DL[0][0],2,2,0,
					&B[0][0],2,numdof,0,
					&Ke->values[0],1);

		TripleMultiply( &B[0][0],2,numdof,1,
					&DLprime[0][0],2,2,0,
					&Bprime[0][0],2,numdof,0,
					&Ke->values[0],1);

		/*Artificial diffusivity*/
		vel=sqrt(pow(vx,2.)+pow(vy,2.));
		K[0][0]=diffusivity*h/(2*vel)*vx*vx;
		K[1][0]=diffusivity*h/(2*vel)*vy*vx;
		K[0][1]=diffusivity*h/(2*vel)*vx*vy;
		K[1][1]=diffusivity*h/(2*vel)*vy*vy;
		KDL[0][0]=DL_scalar*K[0][0];
		KDL[1][0]=DL_scalar*K[1][0];
		KDL[0][1]=DL_scalar*K[0][1];
		KDL[1][1]=DL_scalar*K[1][1];

		TripleMultiply( &Bprime[0][0],2,numdof,1,
					&KDL[0][0],2,2,0,
					&Bprime[0][0],2,numdof,0,
					&Ke->values[0],1);
	}

	/*Clean up and return*/
	delete gauss;
	return Ke;
}
/*}}}*/
/*FUNCTION Tria::CreatePVectorHydrology {{{1*/
ElementVector* Tria::CreatePVectorHydrology(void){

	/*Constants*/
	const int    numdof=NDOF1*NUMVERTICES;

	/*Intermediaries */
	int        i,j,ig;
	double     Jdettria,dt;
	double     basal_melting_g;
	double     old_watercolumn_g;
	double     xyz_list[NUMVERTICES][3];
	double     basis[numdof];
	GaussTria* gauss=NULL;

	/*Skip if water or ice shelf element*/
	if(IsOnWater() | IsFloating()) return NULL;

	/*Initialize Element vector*/
	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);

	/*Retrieve all inputs and parameters*/
	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
	this->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
	Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum); _assert_(basal_melting_input);
	Input* old_watercolumn_input=inputs->GetInput(WaterColumnOldEnum);         _assert_(old_watercolumn_input);

	/*Initialize basal_melting_correction_g to 0, do not forget!:*/
	/* Start  looping on the number of gaussian points: */
	gauss=new GaussTria(2);
	for(ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
		GetNodalFunctions(basis, gauss);

		basal_melting_input->GetInputValue(&basal_melting_g,gauss);
		old_watercolumn_input->GetInputValue(&old_watercolumn_g,gauss);

		if(dt)for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(old_watercolumn_g+dt*basal_melting_g)*basis[i];
		else  for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*basal_melting_g*basis[i];
	}
		
	/*Clean up and return*/
	delete gauss;
	return pe;
}
/*}}}*/
/*FUNCTION Tria::GetSolutionFromInputsHydrology{{{1*/
void  Tria::GetSolutionFromInputsHydrology(Vec solution){

	const int    numdof=NDOF1*NUMVERTICES;

	int i;
	int*         doflist=NULL;
	double       watercolumn;
	double       values[numdof];
	GaussTria*   gauss=NULL;

	/*Get dof list: */
	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);

	/*Get inputs*/
	Input* watercolumn_input=inputs->GetInput(WatercolumnEnum); _assert_(watercolumn_input);

	/*Ok, we have watercolumn values, fill in watercolumn array: */
	/*P1 element only for now*/
	gauss=new GaussTria();
	for(i=0;i<NUMVERTICES;i++){

		gauss->GaussVertex(i);

		/*Recover watercolumn*/
		watercolumn_input->GetInputValue(&watercolumn,gauss);
		values[i]=watercolumn;
	}

	VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES);

	/*Free ressources:*/
	delete gauss;
	xfree((void**)&doflist);
}
/*}}}*/
/*FUNCTION Tria::InputUpdateFromSolutionHydrology{{{1*/
void  Tria::InputUpdateFromSolutionHydrology(double* solution){

	/*Intermediaries*/
	const int numdof = NDOF1*NUMVERTICES;

	int       i;
	int*      doflist=NULL;
	double    values[numdof];

	/*Get dof list: */
	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);

	/*Use the dof list to index into the solution vector: */
	for(i=0;i<numdof;i++){
		values[i]=solution[doflist[i]];
		if(isnan(values[i])) _error_("NaN found in solution vector");
		if (values[i]<pow((double)10,(double)-10))values[i]=pow((double)10,(double)-10); //correcting the water column to positive values
 
	}

	/*Add input to the element: */
	this->inputs->AddInput(new TriaP1Input(WatercolumnEnum,values));

	/*Free ressources:*/
	xfree((void**)&doflist);
}
/*}}}*/
#endif

#ifdef _HAVE_DAKOTA_
/*FUNCTION Tria::InputUpdateFromVectorDakota(double* vector, int name, int type);{{{1*/
void  Tria::InputUpdateFromVectorDakota(double* vector, int name, int type){
	
	int i,j;

	/*Check that name is an element input*/
	if (!IsInput(name)) return;

	switch(type){

		case VertexEnum:

			/*New TriaP1Input*/
			double values[3];

			/*Get values on the 3 vertices*/
			for (i=0;i<3;i++){
				values[i]=vector[this->nodes[i]->GetSidList()]; //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:
					/*Update thickness + surface: assume bed is constant. On ice shelves, takes hydrostatic equilibrium {{{2*/
					double  thickness[3];
					double  thickness_init[3];
					double  hydrostatic_ratio[3];
					double  surface[3];
					double  bed[3];
					
					/*retrieve inputs: */
					GetInputListOnVertices(&thickness_init[0],ThicknessEnum);
					GetInputListOnVertices(&hydrostatic_ratio[0],GeometryHydrostaticRatioEnum);
					GetInputListOnVertices(&bed[0],BedEnum);
					GetInputListOnVertices(&surface[0],SurfaceEnum);

					/*build new thickness: */
//					for(j=0;j<3;j++)thickness[j]=values[j];

					/*build new bed and surface: */
					if (this->IsFloating()){
						/*hydrostatic equilibrium: */
						double rho_ice,rho_water,di;
						rho_ice=this->matpar->GetRhoIce();
						rho_water=this->matpar->GetRhoWater();

						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];
						}

//						for(j=0;j<3;j++){
//							surface[j]=(1-di)*thickness[j];
//							bed[j]=-di*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++)surface[j]=bed[j]+thickness[j];
						for(j=0;j<3;j++)bed[j]=surface[j]-thickness[j];
					}

					/*Add new inputs: */
					this->inputs->AddInput(new TriaP1Input(ThicknessEnum,thickness));
					this->inputs->AddInput(new TriaP1Input(BedEnum,bed));
					this->inputs->AddInput(new TriaP1Input(SurfaceEnum,surface));

					/*}}}*/
					break;
				default:
					this->inputs->AddInput(new TriaP1Input(name,values));
			}
			break;

		default:
			_error_("type %i (%s) not implemented yet",type,EnumToStringx(type));
	}

}
/*}}}*/
/*FUNCTION Tria::InputUpdateFromVectorDakota(int* vector, int name, int type);{{{1*/
void  Tria::InputUpdateFromVectorDakota(int* vector, int name, int type){
	_error_(" not supported yet!");
}
/*}}}*/
/*FUNCTION Tria::InputUpdateFromVectorDakota(bool* vector, int name, int type);{{{1*/
void  Tria::InputUpdateFromVectorDakota(bool* vector, int name, int type){
	_error_(" not supported yet!");
}
/*}}}*/
/*FUNCTION Tria::InputUpdateFromMatrixDakota(double* matrix, int nrows, int ncols, int name, int type);{{{1*/
void  Tria::InputUpdateFromMatrixDakota(double* matrix, int nrows, int ncols, int name, int type){
	
	int i,j,t;
	TransientInput* transientinput=NULL;
	double values[3];
	double time;
	int row;
	double yts;

	/*Check that name is an element input*/
	if (!IsInput(name)) return;

	switch(type){

		case VertexEnum:
			
			/*Create transient input: */
						
			parameters->FindParam(&yts,ConstantsYtsEnum);
			for(t=0;t<ncols;t++){ //ncols is the number of times

				/*create input values: */
				for(i=0;i<3;i++){
					row=this->nodes[i]->GetSidList();
					values[i]=(double)matrix[ncols*row+t];
				}

				/*time? :*/
				time=(double)matrix[(nrows-1)*ncols+t]*yts;

				if(t==0) transientinput=new TransientInput(name);
				transientinput->AddTimeInput(new TriaP1Input(name,values),time);
				transientinput->Configure(parameters);
			}
			this->inputs->AddInput(transientinput);
			break;

		default:
			_error_("type %i (%s) not implemented yet",type,EnumToStringx(type));
	}

}
/*}}}*/
#endif

#ifdef _HAVE_BALANCED_
/*FUNCTION Tria::CreateKMatrixBalancethickness {{{1*/
ElementMatrix* Tria::CreateKMatrixBalancethickness(void){

	switch(GetElementType()){
		case P1Enum:
			return CreateKMatrixBalancethickness_CG();
		case P1DGEnum:
			return CreateKMatrixBalancethickness_DG();
		default:
			_error_("Element type %s not supported yet",EnumToStringx(GetElementType()));
	}

}
/*}}}*/
/*FUNCTION Tria::CreateKMatrixBalancethickness_CG {{{1*/
ElementMatrix* Tria::CreateKMatrixBalancethickness_CG(void){

	/*Constants*/
	const int    numdof=NDOF1*NUMVERTICES;

	/*Intermediaries */
	int        stabilization;
	int        i,j,ig,dim;
	double     Jdettria,vx,vy,dvxdx,dvydy,vel,h;
	double     dvx[2],dvy[2];
	double     xyz_list[NUMVERTICES][3];
	double     L[NUMVERTICES];
	double     B[2][NUMVERTICES];
	double     Bprime[2][NUMVERTICES];
	double     K[2][2]                          = {0.0};
	double     KDL[2][2]                        = {0.0};
	double     DL[2][2]                         = {0.0};
	double     DLprime[2][2]                    = {0.0};
	double     DL_scalar;
	GaussTria *gauss                            = NULL;

	/*Initialize Element matrix*/
	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);

	/*Retrieve all Inputs and parameters: */
	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
	this->parameters->FindParam(&stabilization,BalancethicknessStabilizationEnum);
	this->parameters->FindParam(&dim,MeshDimensionEnum);
	Input* vxaverage_input=NULL;
	Input* vyaverage_input=NULL;
	if(dim==2){
		vxaverage_input=inputs->GetInput(VxEnum); _assert_(vxaverage_input);
		vyaverage_input=inputs->GetInput(VyEnum); _assert_(vyaverage_input);
	}
	else{
		vxaverage_input=inputs->GetInput(VxAverageEnum); _assert_(vxaverage_input);
		vyaverage_input=inputs->GetInput(VyAverageEnum); _assert_(vyaverage_input);
	}
	h=sqrt(2*this->GetArea());

	/*Start looping on the number of gaussian points:*/
	gauss=new GaussTria(2);
	for (ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
		GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss);
		GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);

		vxaverage_input->GetInputValue(&vx,gauss);
		vyaverage_input->GetInputValue(&vy,gauss);
		vxaverage_input->GetInputDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
		vyaverage_input->GetInputDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);

		dvxdx=dvx[0];
		dvydy=dvy[1];
		DL_scalar=gauss->weight*Jdettria;

		DL[0][0]=DL_scalar*dvxdx;
		DL[1][1]=DL_scalar*dvydy;

		DLprime[0][0]=DL_scalar*vx;
		DLprime[1][1]=DL_scalar*vy;

		TripleMultiply( &B[0][0],2,numdof,1,
					&DL[0][0],2,2,0,
					&B[0][0],2,numdof,0,
					&Ke->values[0],1);

		TripleMultiply( &B[0][0],2,numdof,1,
					&DLprime[0][0],2,2,0,
					&Bprime[0][0],2,numdof,0,
					&Ke->values[0],1);

		if(stabilization==1){
			/*Streamline upwinding*/
			vel=sqrt(pow(vx,2.)+pow(vy,2.));
			K[0][0]=h/(2*vel)*vx*vx;
			K[1][0]=h/(2*vel)*vy*vx;
			K[0][1]=h/(2*vel)*vx*vy;
			K[1][1]=h/(2*vel)*vy*vy;
		}
		else if(stabilization==2){
			/*MacAyeal*/
			vxaverage_input->GetInputAverage(&vx);
			vyaverage_input->GetInputAverage(&vy);
			K[0][0]=h/2.0*fabs(vx);
			K[0][1]=0.;
			K[1][0]=0.;
			K[1][1]=h/2.0*fabs(vy);
		}
		if(stabilization==1 || stabilization==2){
			KDL[0][0]=DL_scalar*K[0][0];
			KDL[1][0]=DL_scalar*K[1][0];
			KDL[0][1]=DL_scalar*K[0][1];
			KDL[1][1]=DL_scalar*K[1][1];
			TripleMultiply( &Bprime[0][0],2,numdof,1,
						&KDL[0][0],2,2,0,
						&Bprime[0][0],2,numdof,0,
						&Ke->values[0],1);
		}
	}

	/*Clean up and return*/
	delete gauss;
	return Ke;
}
/*}}}*/
/*FUNCTION Tria::CreateKMatrixBalancethickness_DG {{{1*/
ElementMatrix* Tria::CreateKMatrixBalancethickness_DG(void){

	/*Constants*/
	const int  numdof=NDOF1*NUMVERTICES;

	/*Intermediaries*/
	int        i,j,ig,dim;
	double     vx,vy,Jdettria;
	double     xyz_list[NUMVERTICES][3];
	double     B[2][NUMVERTICES];
	double     Bprime[2][NUMVERTICES];
	double     DL[2][2]={0.0};
	double     DL_scalar;
	GaussTria  *gauss=NULL;

	/*Initialize Element matrix*/
	ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);

	/*Retrieve all inputs and parameters*/
	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
	this->parameters->FindParam(&dim,MeshDimensionEnum);
	Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input);
	Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input);

	/*Start looping on the number of gaussian points:*/
	gauss=new GaussTria(2);
	for (ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
		/*WARNING: B and Bprime are inverted compared to usual prognostic!!!!*/
		GetBPrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
		GetBprimePrognostic(&B[0][0], &xyz_list[0][0], gauss);

		vx_input->GetInputValue(&vx,gauss);
		vy_input->GetInputValue(&vy,gauss);

		DL_scalar=-gauss->weight*Jdettria;
		DL[0][0]=DL_scalar*vx;
		DL[1][1]=DL_scalar*vy;

		TripleMultiply( &B[0][0],2,numdof,1,
					&DL[0][0],2,2,0,
					&Bprime[0][0],2,numdof,0,
					&Ke->values[0],1);
	}

	/*Clean up and return*/
	delete gauss;
	return Ke;
}
/*}}}*/
/*FUNCTION Tria::CreatePVectorBalancethickness{{{1*/
ElementVector* Tria::CreatePVectorBalancethickness(void){

	switch(GetElementType()){
		case P1Enum:
			return CreatePVectorBalancethickness_CG();
			break;
		case P1DGEnum:
			return CreatePVectorBalancethickness_DG();
		default:
			_error_("Element type %s not supported yet",EnumToStringx(GetElementType()));
	}
}
/*}}}*/
/*FUNCTION Tria::CreatePVectorBalancethickness_CG{{{1*/
ElementVector* Tria::CreatePVectorBalancethickness_CG(void){

	/*Constants*/
	const int    numdof=NDOF1*NUMVERTICES;
	
	/*Intermediaries */
	int        i,j,ig;
	double     xyz_list[NUMVERTICES][3];
	double     dhdt_g,basal_melting_g,surface_mass_balance_g,Jdettria;
	double     L[NUMVERTICES];
	GaussTria* gauss=NULL;

	/*Initialize Element vector*/
	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);

	/*Retrieve all inputs and parameters*/
	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
	Input* surface_mass_balance_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(surface_mass_balance_input);
	Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum);          _assert_(basal_melting_input);
	Input* dhdt_input=inputs->GetInput(BalancethicknessThickeningRateEnum);             _assert_(dhdt_input);
	
	/* Start  looping on the number of gaussian points: */
	gauss=new GaussTria(2);
	for(ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		surface_mass_balance_input->GetInputValue(&surface_mass_balance_g,gauss);
		basal_melting_input->GetInputValue(&basal_melting_g,gauss);
		dhdt_input->GetInputValue(&dhdt_g,gauss);

		GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
		GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);

		for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*L[i];
	}

	/*Clean up and return*/
	delete gauss;
	return pe;
}
/*}}}*/
/*FUNCTION Tria::CreatePVectorBalancethickness_DG {{{1*/
ElementVector* Tria::CreatePVectorBalancethickness_DG(void){

	/*Constants*/
	const int    numdof=NDOF1*NUMVERTICES;

	/*Intermediaries */
	int        i,j,ig;
	double     xyz_list[NUMVERTICES][3];
	double     basal_melting_g,surface_mass_balance_g,dhdt_g,Jdettria;
	double     L[NUMVERTICES];
	GaussTria* gauss=NULL;

	/*Initialize Element vector*/
	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);

	/*Retrieve all inputs and parameters*/
	GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
	Input* surface_mass_balance_input=inputs->GetInput(SurfaceforcingsMassBalanceEnum); _assert_(surface_mass_balance_input);
	Input* basal_melting_input=inputs->GetInput(BasalforcingsMeltingRateEnum);          _assert_(basal_melting_input);
	Input* dhdt_input=inputs->GetInput(BalancethicknessThickeningRateEnum);                                       _assert_(dhdt_input);

	/* Start  looping on the number of gaussian points: */
	gauss=new GaussTria(2);
	for(ig=gauss->begin();ig<gauss->end();ig++){

		gauss->GaussPoint(ig);

		surface_mass_balance_input->GetInputValue(&surface_mass_balance_g,gauss);
		basal_melting_input->GetInputValue(&basal_melting_g,gauss);
		dhdt_input->GetInputValue(&dhdt_g,gauss);

		GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
		GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);

		for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(surface_mass_balance_g-basal_melting_g-dhdt_g)*L[i];
	}

	/*Clean up and return*/
	delete gauss;
	return pe;
}
/*}}}*/
#endif
