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

#ifdef HAVE_CONFIG_H
	#include "config.h"
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#include "stdio.h"
#include "./Tria.h"
#include <string.h>
#include "../EnumDefinitions/EnumDefinitions.h"
#include "../shared/shared.h"
#include "../DataSet/DataSet.h"
#include "../include/typedefs.h"


/*For debugging purposes: */
#define ELID 141 //element for which to print debug statements
#define RANK 2 //rank of cpu for which to print debug statements.
//#define _DEBUGELEMENTS_
//#define _DEBUGGAUSS_

Tria::Tria(){
	return;
}

Tria::~Tria(){
	return;
}

Tria::Tria(int tria_id,int tria_mid,int tria_mparid,int tria_node_ids[3],double tria_h[3],double tria_s[3],double tria_b[3],double tria_k[3],
				int tria_friction_type,double tria_p,double tria_q,int tria_shelf,double tria_fit,double tria_meanvel,double tria_epsvel,int tria_acceleration){
	
	int i;
	
	id=tria_id;
	mid=tria_mid;
	mparid=tria_mparid;
	for(i=0;i<3;i++){
		node_ids[i]=tria_node_ids[i];
		node_offsets[i]=UNDEF;
		nodes[i]=NULL;
		h[i]=tria_h[i];
		s[i]=tria_s[i];
		b[i]=tria_b[i];
		k[i]=tria_k[i];
	}
	matice=NULL;
	matice_offset=UNDEF;
	matpar=NULL;
	matpar_offset=UNDEF;
	friction_type=tria_friction_type;
	p=tria_p;
	q=tria_q;
	shelf=tria_shelf;
	fit=tria_fit;
	meanvel=tria_meanvel;
	epsvel=tria_epsvel;
	acceleration=tria_acceleration;
	onbed=1;
	
	return;
}

#undef __FUNCT__
#define __FUNCT__ "Tria::Echo"

void Tria::Echo(void){

	printf("Tria:\n");
	printf("   id: %i\n",id);
	printf("   mid: %i\n",mid);
	printf("   mparid: %i\n",mparid);
	printf("   node_ids=[%i,%i,%i]\n",node_ids[0],node_ids[1],node_ids[2]);
	printf("   node_offsets=[%i,%i,%i]\n",node_offsets[0],node_offsets[1],node_offsets[2]);
	printf("   matice_offset=%i\n",matice_offset);
	printf("   matpar_offset=%i\n",matpar_offset);
	printf("   h=[%g,%g,%g]\n",h[0],h[1],h[2]);
	printf("   s=[%g,%g,%g]\n",s[0],s[1],s[2]);
	printf("   b=[%g,%g,%g]\n",b[0],b[1],b[2]);
	printf("   k=[%g,%g,%g]\n",k[0],k[1],k[2]);
	printf("   friction_type: %i\n",friction_type);
	printf("   p: %g\n",p);
	printf("   q: %g\n",q);
	printf("   shelf: %i\n",shelf);
	printf("   fit: %g\n",fit);
	printf("   meanvel: %g\n",meanvel);
	printf("   epsvel: %g\n",epsvel);
	printf("   acceleration: %i\n",acceleration);
	printf("   onbed: %i\n",onbed);
	printf("   nodes: \n");
	if(nodes[0])nodes[0]->Echo();
	if(nodes[1])nodes[1]->Echo();
	if(nodes[2])nodes[2]->Echo();
	if(matice)matice->Echo();
	if(matpar)matpar->Echo();

	return;
}

void  Tria::Marshall(char** pmarshalled_dataset){

	char* marshalled_dataset=NULL;
	int   enum_type=0;

	/*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,&mid,sizeof(mid));marshalled_dataset+=sizeof(mid);
	memcpy(marshalled_dataset,&mparid,sizeof(mparid));marshalled_dataset+=sizeof(mparid);
	memcpy(marshalled_dataset,&node_ids,sizeof(node_ids));marshalled_dataset+=sizeof(node_ids);
	memcpy(marshalled_dataset,&nodes,sizeof(nodes));marshalled_dataset+=sizeof(nodes);
	memcpy(marshalled_dataset,&node_offsets,sizeof(node_offsets));marshalled_dataset+=sizeof(node_offsets);
	memcpy(marshalled_dataset,&matice,sizeof(matice));marshalled_dataset+=sizeof(matice);
	memcpy(marshalled_dataset,&matice_offset,sizeof(matice_offset));marshalled_dataset+=sizeof(matice_offset);
	memcpy(marshalled_dataset,&matpar,sizeof(matpar));marshalled_dataset+=sizeof(matpar);
	memcpy(marshalled_dataset,&matpar_offset,sizeof(matpar_offset));marshalled_dataset+=sizeof(matpar_offset);
	memcpy(marshalled_dataset,&h,sizeof(h));marshalled_dataset+=sizeof(h);
	memcpy(marshalled_dataset,&s,sizeof(s));marshalled_dataset+=sizeof(s);
	memcpy(marshalled_dataset,&b,sizeof(b));marshalled_dataset+=sizeof(b);
	memcpy(marshalled_dataset,&k,sizeof(k));marshalled_dataset+=sizeof(k);
	memcpy(marshalled_dataset,&friction_type,sizeof(friction_type));marshalled_dataset+=sizeof(friction_type);
	memcpy(marshalled_dataset,&onbed,sizeof(onbed));marshalled_dataset+=sizeof(onbed);
	memcpy(marshalled_dataset,&p,sizeof(p));marshalled_dataset+=sizeof(p);
	memcpy(marshalled_dataset,&q,sizeof(q));marshalled_dataset+=sizeof(q);
	memcpy(marshalled_dataset,&shelf,sizeof(shelf));marshalled_dataset+=sizeof(shelf);
	memcpy(marshalled_dataset,&fit,sizeof(fit));marshalled_dataset+=sizeof(fit);
	memcpy(marshalled_dataset,&meanvel,sizeof(meanvel));marshalled_dataset+=sizeof(meanvel);
	memcpy(marshalled_dataset,&epsvel,sizeof(epsvel));marshalled_dataset+=sizeof(epsvel);
	memcpy(marshalled_dataset,&acceleration,sizeof(acceleration));marshalled_dataset+=sizeof(acceleration);
	
	*pmarshalled_dataset=marshalled_dataset;
	return;
}
		
int   Tria::MarshallSize(){
	return sizeof(id)
		+sizeof(mid)
		+sizeof(mparid)
		+sizeof(node_ids)
		+sizeof(nodes)
		+sizeof(node_offsets)
		+sizeof(matice)
		+sizeof(matice_offset)
		+sizeof(matpar)
		+sizeof(matpar_offset)
		+sizeof(h)
		+sizeof(s)
		+sizeof(b)
		+sizeof(k)
		+sizeof(friction_type)
		+sizeof(onbed)
		+sizeof(p)
		+sizeof(q)
		+sizeof(shelf)
		+sizeof(fit)
		+sizeof(meanvel)
		+sizeof(epsvel)
		+sizeof(acceleration)
		+sizeof(int); //sizeof(int) for enum type
}

char* Tria::GetName(void){
	return "tria";
}

void  Tria::Demarshall(char** pmarshalled_dataset){

	char* marshalled_dataset=NULL;
	int   i;

	/*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(&mid,marshalled_dataset,sizeof(mid));marshalled_dataset+=sizeof(mid);
	memcpy(&mparid,marshalled_dataset,sizeof(mparid));marshalled_dataset+=sizeof(mparid);
	memcpy(&node_ids,marshalled_dataset,sizeof(node_ids));marshalled_dataset+=sizeof(node_ids);
	memcpy(&nodes,marshalled_dataset,sizeof(nodes));marshalled_dataset+=sizeof(nodes);
	memcpy(&node_offsets,marshalled_dataset,sizeof(node_offsets));marshalled_dataset+=sizeof(node_offsets);
	memcpy(&matice,marshalled_dataset,sizeof(matice));marshalled_dataset+=sizeof(matice);
	memcpy(&matice_offset,marshalled_dataset,sizeof(matice_offset));marshalled_dataset+=sizeof(matice_offset);
	memcpy(&matpar,marshalled_dataset,sizeof(matpar));marshalled_dataset+=sizeof(matpar);
	memcpy(&matpar_offset,marshalled_dataset,sizeof(matpar_offset));marshalled_dataset+=sizeof(matpar_offset);
	memcpy(&h,marshalled_dataset,sizeof(h));marshalled_dataset+=sizeof(h);
	memcpy(&s,marshalled_dataset,sizeof(s));marshalled_dataset+=sizeof(s);
	memcpy(&b,marshalled_dataset,sizeof(b));marshalled_dataset+=sizeof(b);
	memcpy(&k,marshalled_dataset,sizeof(k));marshalled_dataset+=sizeof(k);
	memcpy(&friction_type,marshalled_dataset,sizeof(friction_type));marshalled_dataset+=sizeof(friction_type);
	memcpy(&onbed,marshalled_dataset,sizeof(onbed));marshalled_dataset+=sizeof(onbed);
	memcpy(&p,marshalled_dataset,sizeof(p));marshalled_dataset+=sizeof(p);
	memcpy(&q,marshalled_dataset,sizeof(q));marshalled_dataset+=sizeof(q);
	memcpy(&shelf,marshalled_dataset,sizeof(shelf));marshalled_dataset+=sizeof(shelf);
	memcpy(&fit,marshalled_dataset,sizeof(fit));marshalled_dataset+=sizeof(fit);
	memcpy(&meanvel,marshalled_dataset,sizeof(meanvel));marshalled_dataset+=sizeof(meanvel);
	memcpy(&epsvel,marshalled_dataset,sizeof(epsvel));marshalled_dataset+=sizeof(epsvel);
	memcpy(&acceleration,marshalled_dataset,sizeof(acceleration));marshalled_dataset+=sizeof(acceleration);

	/*nodes and materials are not pointing to correct objects anymore:*/
	for(i=0;i<3;i++)nodes[i]=NULL;
	matice=NULL;
	matpar=NULL;

	/*return: */
	*pmarshalled_dataset=marshalled_dataset;
	return;
}
int Tria::Enum(void){

	return TriaEnum();

}
int    Tria::GetId(void){ return id; }

int    Tria::MyRank(void){ 
	extern int my_rank;
	return my_rank; 
}


#undef __FUNCT__ 
#define __FUNCT__ "Tria::Configure"
void  Tria::Configure(void* ploadsin,void* pnodesin,void* pmaterialsin){

	int i;
	
	DataSet* loadsin=NULL;
	DataSet* nodesin=NULL;
	DataSet* materialsin=NULL;

	/*Recover pointers :*/
	loadsin=(DataSet*)ploadsin;
	nodesin=(DataSet*)pnodesin;
	materialsin=(DataSet*)pmaterialsin;

	/*Link this element with its nodes, ie find pointers to the nodes in the nodes dataset.: */
	ResolvePointers((Object**)nodes,node_ids,node_offsets,3,nodesin);
	
	/*Same for materials: */
	ResolvePointers((Object**)&matice,&mid,&matice_offset,1,materialsin);
	ResolvePointers((Object**)&matpar,&mparid,&matpar_offset,1,materialsin);

}

#undef __FUNCT__ 
#define __FUNCT__ "Tria::CreateKMatrix"

void  Tria::CreateKMatrix(Mat Kgg,ParameterInputs* inputs,int analysis_type){

	/*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
	if ((analysis_type==DiagnosticHorizAnalysisEnum()) || (analysis_type==ControlAnalysisEnum())){

		CreateKMatrixDiagnosticHoriz( Kgg,inputs,analysis_type);

	}
	else{
		throw ErrorException(__FUNCT__,exprintf("%s%i%s\n","analysis: ",analysis_type," not supported yet"));
	}

}


#undef __FUNCT__ 
#define __FUNCT__ "Tria::CreateKMatrixDiagnosticHoriz"

void  Tria::CreateKMatrixDiagnosticHoriz(Mat Kgg,ParameterInputs* inputs,int analysis_type){


	/* local declarations */
	int             i,j;

	/* node data: */
	const int    numgrids=3;
	const int    numdofs=2*numgrids;
	double       xyz_list[numgrids][3];
	int          doflist[numdofs];
	int          numberofdofspernode;
	
	/* gaussian points: */
	int     num_gauss,ig;
	double* first_gauss_area_coord  =  NULL;
	double* second_gauss_area_coord =  NULL;
	double* third_gauss_area_coord  =  NULL;
	double* gauss_weights           =  NULL;
	double  gauss_weight;
	double  gauss_l1l2l3[3];

	/* material data: */
	double viscosity; //viscosity
	
	/* strain rate: */
	double epsilon[3]; /* epsilon=[exx,eyy,exy];*/

	/* matrices: */
	double B[3][numdofs];
	double Bprime[3][numdofs];
	double L[2][numdofs];
	double D[3][3]={{ 0,0,0 },{0,0,0},{0,0,0}};              // material matrix, simple scalar matrix.
	double D_scalar;
	double DL[2][2]={{ 0,0 },{0,0}}; //for basal drag
	double DL_scalar;

	/* local element matrices: */
	double Ke_gg[numdofs][numdofs]; //local element stiffness matrix 
	double Ke_gg_gaussian[numdofs][numdofs]; //stiffness matrix evaluated at the gaussian point.
	double Ke_gg_drag_gaussian[numdofs][numdofs]; //stiffness matrix contribution from drag
	
	double Jdet;
	
	/*slope: */
	double  slope[2]={0.0,0.0};
	double  slope_magnitude;

	/*input parameters for structural analysis (diagnostic): */
	double* velocity_param=NULL;
	double  vxvy_list[numgrids][2];
	double  thickness;

	/*friction: */
	double alpha2_list[numgrids]={0.0,0.0,0.0};
	double alpha2;

	double MAXSLOPE=.06; // 6 %
	double MOUNTAINKEXPONENT=10;


	/*First check that acceleration is not on: */
	if (acceleration){
		/*do nothing: */
		return;
	}

	/*recover extra inputs from users, at current convergence iteration: */
	velocity_param=ParameterInputsRecover(inputs,"velocity"); 

	/* Get node coordinates and dof list: */
	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
	GetDofList(&doflist[0],&numberofdofspernode);

	/* Set Ke_gg to 0: */
	for(i=0;i<numdofs;i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]=0.0;


	/*Initialize vxvy_list: */
	for(i=0;i<numgrids;i++){
		if(velocity_param){
			vxvy_list[i][0]=velocity_param[doflist[i*numberofdofspernode+0]];
			vxvy_list[i][1]=velocity_param[doflist[i*numberofdofspernode+1]];
		}
		else{
			vxvy_list[i][0]=0;
			vxvy_list[i][1]=0;
		}
	}

	#ifdef _DEBUGELEMENTS_
	if(my_rank==RANK && id==ELID){ 
		printf("El id %i Rank %i TriaElemnet input list before gaussian loop: \n",ELID,RANK); 
		printf("   rho_ice: %g \n",matpar->GetRhoIce());
		printf("   gravity: %g \n",matpar->GetG())
		printf("   rho_water: %g \n",matpar->GetRhoWater());
		printf("   Velocity: \n");
		for (i=0;i<numgrids;i++){
			printf("      node %i  [%g,%g]\n",i,vxvy_list[i][0],vxvy_list[i][1]);
		}
		printf("   flow_law_parameter [%g ]\n",matice->GetB());
		printf("   drag [%g %g %g ]\n",k[0],k[1],k[2]);
		printf("   thickness [%g %g %g]\n",h[0],h[1],h[2]);
		printf("   surface [%g %g %g ]\n",s[0],s[1],s[2]);
		printf("   bed [%g %g %g]\n",b[0],b[1],b[2]);
	}
	#endif


	/*Build alpha2_list used by drag stiffness matrix*/
	if (!shelf && (friction_type==2)){
		
		/*Allocate friction object: */
		Friction* friction=NewFriction();
		
		/*Initialize all fields: */
		friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char));
		strcpy(friction->element_type,"2d");
		
		friction->gravity=matpar->GetG();
		friction->rho_ice=matpar->GetRhoIce();
		friction->rho_water=matpar->GetRhoWater();
		friction->K=&k[0];
		friction->bed=&b[0];
		friction->thickness=&h[0];
		friction->velocities=&vxvy_list[0][0];
		friction->p=p;
		friction->q=q;

		/*Compute alpha2_list: */
		FrictionGetAlpha2(&alpha2_list[0],friction);

		/*Erase friction object: */
		DeleteFriction(&friction);
	}

	#ifdef _DEBUGELEMENTS_
	if(my_rank==RANK && id==ELID){ 
		printf("   alpha2_list [%g %g %g ]\n",alpha2_list[0],alpha2_list[1],alpha2_list[2]);
	}
	#endif

	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
	GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);

	#ifdef _DEBUGELEMENTS_
	if(my_rank==RANK && id==ELID){ 
		printf("   gaussian points: \n");
		for(i=0;i<num_gauss;i++){
			printf("    %g %g %g : %g\n",first_gauss_area_coord[i],second_gauss_area_coord[i],third_gauss_area_coord[i],gauss_weights[i]);
		}
	}
	#endif

	/* Start  looping on the number of gaussian points: */
	for (ig=0; ig<num_gauss; ig++){
		/*Pick up the gaussian point: */
		gauss_weight=*(gauss_weights+ig);
		gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 
		gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
		gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);


		/*Compute thickness at gaussian point: */
		GetParameterValue(&thickness, &h[0],gauss_l1l2l3);

		// 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: */
		if(!shelf){
			GetParameterDerivativeValue(&slope[0], &s[0],&xyz_list[0][0], gauss_l1l2l3);
			slope_magnitude=sqrt(pow(slope[0],2)+pow(slope[1],2));

			if (slope_magnitude>MAXSLOPE){
				alpha2_list[0]=pow(10,MOUNTAINKEXPONENT);
				alpha2_list[1]=pow(10,MOUNTAINKEXPONENT);
				alpha2_list[2]=pow(10,MOUNTAINKEXPONENT);
			}
		}

		/*Get strain rate from velocity: */
		GetStrainRate(&epsilon[0],&vxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3);
		
		/*Get viscosity: */
		matice->GetViscosity2d(&viscosity, &epsilon[0]);
		
		/* Get Jacobian determinant: */
		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3);

		/* Build the D matrix: we plug the gaussian weight, the thickness, the viscosity, and the jacobian determinant 
		   onto this scalar matrix, so that we win some computational time: */
		D_scalar=viscosity*thickness*gauss_weight*Jdet;
		for (i=0;i<3;i++){
			D[i][i]=D_scalar;
		}
		
		#ifdef _DEBUGELEMENTS_
		if(my_rank==RANK && id==ELID){ 
			printf("   gaussian loop %i\n",ig);
			printf("      thickness %g\n",thickness);
			printf("      slope [%g,%g]\n",slope[0],slope[1]);
			printf("      alpha2_list [%g,%g,%g]\n",alpha2_list[0],alpha2_list[1],alpha2_list[2]);
			printf("      epsilon [%g,%g,%g]\n",epsilon[0],epsilon[1],epsilon[2]);
			printf("      Matice: \n");
			matice->Echo();
			printf("      Matpar: \n");
			matpar->Echo();
			printf("\n    viscosity: %g \n",viscosity);
			printf("      jacobian: %g \n",Jdet);
			printf("      gauss_weight: %g \n",gauss_weight);
		}
		#endif

		/*Get B and Bprime matrices: */
		GetB(&B[0][0], &xyz_list[0][0], gauss_l1l2l3);
		GetBPrime(&Bprime[0][0], &xyz_list[0][0], gauss_l1l2l3);

		/*Get L matrix if viscous basal drag present: */
		if((friction_type==2) && (!shelf)){
			GetL(&L[0][0], &xyz_list[0][0], gauss_l1l2l3,numberofdofspernode);
		}	

		/*  Do the triple product tB*D*Bprime: */
		TripleMultiply( &B[0][0],3,numdofs,1,
					  &D[0][0],3,3,0,
					  &Bprime[0][0],3,numdofs,0,
					  &Ke_gg_gaussian[0][0],0);

		/* Add the Ke_gg_gaussian, and optionally Ke_gg_drag_gaussian onto Ke_gg: */
		for( i=0; i<numdofs; i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
		
		/*Now, take care of the basal friction if there is any: */
		if((!shelf) && (friction_type==2)){

			GetParameterValue(&alpha2, &alpha2_list[0],gauss_l1l2l3);

			if (velocity_param){
				DL_scalar=alpha2*gauss_weight*Jdet;
			}
			else{
				DL_scalar=0;
			}
				
			for (i=0;i<2;i++){
				DL[i][i]=DL_scalar;
			}
			
			/*  Do the triple producte tL*D*L: */
			TripleMultiply( &L[0][0],2,numdofs,1,
						&DL[0][0],2,2,0,
						&L[0][0],2,numdofs,0,
						&Ke_gg_drag_gaussian[0][0],0);

			for( i=0; i<numdofs; i++) for(j=0;j<numdofs;j++) Ke_gg[i][j]+=Ke_gg_drag_gaussian[i][j];
		}


		#ifdef _DEBUGELEMENTS_
		if(my_rank==RANK && id==ELID){ 
			printf("      alpha2 %g\n",alpha2);
			printf("      B:\n");
			for(i=0;i<3;i++){
				for(j=0;j<numdofs;j++){
					printf("%g ",B[i][j]);
				}
				printf("\n");
			}
			printf("      Bprime:\n");
			for(i=0;i<3;i++){
				for(j=0;j<numdofs;j++){
					printf("%g ",Bprime[i][j]);
				}
				printf("\n");
			}
			printf("      L:\n");
			for(i=0;i<2;i++){
				for(j=0;j<numdofs;j++){
					printf("%g ",L[i][j]);
				}
				printf("\n");
			}

		}
		#endif
	} // for (ig=0; ig<num_gauss; ig++)

	/*Add Ke_gg to global matrix Kgg: */
	MatSetValues(Kgg,numdofs,doflist,numdofs,doflist,(const double*)Ke_gg,ADD_VALUES);

	#ifdef _DEBUGELEMENTS_
	if(my_rank==RANK && id==ELID){ 
		printf("      Ke_gg erms:\n");
		for( i=0; i<numdofs; i++){
			for (j=0;j<numdofs;j++){
				printf("%g ",Ke_gg[i][j]);
			}
			printf("\n");
		}
		printf("      Ke_gg row_indices:\n");
		for( i=0; i<numdofs; i++){
			printf("%i ",doflist[i]);
		}

	}
	#endif

	cleanup_and_return: 
	xfree((void**)&first_gauss_area_coord);
	xfree((void**)&second_gauss_area_coord);
	xfree((void**)&third_gauss_area_coord);
	xfree((void**)&gauss_weights);

}
	

#undef __FUNCT__ 
#define __FUNCT__ "Tria::CreatePVector"
void  Tria::CreatePVector(Vec pg,ParameterInputs* inputs,int analysis_type){
	
	/*Just branch to the correct load generator, according to the type of analysis we are carrying out: */
	if ((analysis_type==DiagnosticHorizAnalysisEnum()) || (analysis_type==ControlAnalysisEnum())){
		CreatePVectorDiagnosticHoriz( pg,inputs,analysis_type);
	}
	else{
		throw ErrorException(__FUNCT__,exprintf("%s%i%s"," analysis ",analysis_type," not supported yet"));
	}

}



#undef __FUNCT__ 
#define __FUNCT__ "Tria::CreatePVectorDiagnosticHoriz"

void Tria::CreatePVectorDiagnosticHoriz( Vec pg, ParameterInputs* inputs, int analysis_type){

	int             i,j;

	/* node data: */
	const int    numgrids=3;
	const int    numdofs=2*numgrids;
	const int    NDOF2=2;
	double       xyz_list[numgrids][3];
	int          doflist[numdofs];
	int          numberofdofspernode;
	
	/* parameters: */
	double  plastic_stress; 
	double  slope[NDOF2];
	double  driving_stress_baseline;

	/* gaussian points: */
	int     num_gauss,ig;
	double* first_gauss_area_coord  =  NULL;
	double* second_gauss_area_coord =  NULL;
	double* third_gauss_area_coord  =  NULL;
	double* gauss_weights           =  NULL;
	double  gauss_weight;
	double  gauss_l1l2l3[3];

	/* Jacobian: */
	double Jdet;

	/*nodal functions: */
	double l1l2l3[3];

	/*element vector at the gaussian points: */
	double  pe_g[numdofs];
	double  pe_g_gaussian[numdofs];

	/*input parameters for structural analysis (diagnostic): */
	double  thickness;

	/*First check that acceleration is not on: */
	if (acceleration){
		return;
	}

	/* Get node coordinates and dof list: */
	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
	GetDofList(&doflist[0],&numberofdofspernode);

	/* Set pe_g to 0: */
	for(i=0;i<numdofs;i++) pe_g[i]=0.0;


	#ifdef _DEBUGELEMENTS_
	if(my_rank==RANK && id==ELID){ 
		printf("gravity %g\n",matpar->GetG());
		printf("rho_ice %g\n",matpar->GetRhoIce());
		printf("thickness [%g,%g,%g]\n",h[0],h[1],h[2]);
		printf("surface[%g,%g,%g]\n",s[0],s[1],s[2]);
		printf("bed[%g,%g,%g]\n",b[0],b[1],b[2]);
		printf("drag [%g,%g,%g]\n",k[0],k[1],k[2]);
	}
	#endif


	/* Get gaussian points and weights: */
	GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2); /*We need higher order because our load is order 2*/

	#ifdef _DEBUGELEMENTS_
	if(my_rank==RANK && id==ELID){ 
		printf("   gaussian points: \n");
		for(i=0;i<num_gauss;i++){
			printf("    %g %g %g : %g\n",first_gauss_area_coord[i],second_gauss_area_coord[i],third_gauss_area_coord[i],gauss_weights[i]);
		}
	}
	#endif



	/* Start  looping on the number of gaussian points: */
	for (ig=0; ig<num_gauss; ig++){
		/*Pick up the gaussian point: */
		gauss_weight=*(gauss_weights+ig);
		gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 
		gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
		gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);

		/*Compute thickness at gaussian point: */
		GetParameterValue(&thickness, &h[0],gauss_l1l2l3);
	
		GetParameterDerivativeValue(&slope[0], &s[0],&xyz_list[0][0], gauss_l1l2l3);
		
		/*In case we have plastic basal drag, compute plastic stress at gaussian point from k1, k2 and k3 fields in the 
		 * element itself: */
		if(friction_type==1){
			GetParameterValue(&plastic_stress, &k[0],gauss_l1l2l3);
		}

		/* Get Jacobian determinant: */
		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
		
		 /*Get nodal functions: */
		GetNodalFunctions(l1l2l3, gauss_l1l2l3);

		/*Compute driving stress: */
		driving_stress_baseline=matpar->GetRhoIce()*matpar->GetG()*thickness;


	 	#ifdef _DEBUGELEMENTS_
		if(my_rank==RANK && id==ELID){ 
			printf("      gaussian %i\n",ig);
			printf("      thickness %g\n",thickness);
			printf("      slope(%g,%g)\n",slope[0],slope[1]);
			printf("      Jdet %g\n",Jdet);
			printf("      gaussweigth %g\n",gauss_weight);
			printf("      l1l2l3 (%g,%g,%g)\n",l1l2l3[0],l1l2l3[1],l1l2l3[2]);
			if(friction_type==1)printf("      plastic_stress(%g)\n",plastic_stress);
		}
		#endif

		/*Build pe_g_gaussian vector: */
		if(friction_type==1){
			for (i=0;i<numgrids;i++){
				for (j=0;j<NDOF2;j++){
					pe_g_gaussian[i*NDOF2+j]=(-driving_stress_baseline*slope[j]-plastic_stress)*Jdet*gauss_weight*l1l2l3[i]; 
				}
			}
		}
		else {
			for (i=0;i<numgrids;i++){
				for (j=0;j<NDOF2;j++){
					pe_g_gaussian[i*NDOF2+j]=-driving_stress_baseline*slope[j]*Jdet*gauss_weight*l1l2l3[i];
				}
			}
		}

		/*Add pe_g_gaussian vector to pe_g: */
		for( i=0; i<numdofs; i++)pe_g[i]+=pe_g_gaussian[i];

	} //for (ig=0; ig<num_gauss; ig++)

	#ifdef _DEBUGELEMENTS_
	if(my_rank==RANK && id==ELID){ 
		printf("      pe_g->terms\n",ig);
		for( i=0; i<pe_g->nrows; i++){
			printf("%g ",*(pe_g->terms+i));
		}
		printf("\n");
		printf("      pe_g->row_indices\n",ig);
		for( i=0; i<pe_g->nrows; i++){
			printf("%i ",*(pe_g->row_indices+i));
		}
	}
	#endif

	/*Add pe_g to global vector pg: */
	VecSetValues(pg,numdofs,doflist,(const double*)pe_g,ADD_VALUES);

	cleanup_and_return: 
	xfree((void**)&first_gauss_area_coord);
	xfree((void**)&second_gauss_area_coord);
	xfree((void**)&third_gauss_area_coord);
	xfree((void**)&gauss_weights);

}


#undef __FUNCT__ 
#define __FUNCT__ "Tria::UpdateFromInputs"
void  Tria::UpdateFromInputs(ParameterInputs* inputs){

	int i;
	int doflist[MAXDOFSPERNODE*3];
	int numberofdofspernode;

	double* thickness=NULL;
	double* bed=NULL;
	double* surface=NULL;
	double* drag=NULL;
	double* temperature_average_param=NULL;
	double* flow_law_param=NULL;
	double temperature_average;
	double B_average;



	/*Get dof list: */
	GetDofList(&doflist[0],&numberofdofspernode);

	/*Update internal data if inputs holds new values: */
	thickness=ParameterInputsRecover(inputs,"thickness");
	bed=ParameterInputsRecover(inputs,"bed");
	surface=ParameterInputsRecover(inputs,"surface");
	drag=ParameterInputsRecover(inputs,"drag");

	for(i=0;i<3;i++){
		if(thickness)h[i]=thickness[doflist[i*numberofdofspernode+0]];
		if(bed)b[i]=bed[doflist[i*numberofdofspernode+0]];
		if(surface)s[i]=surface[doflist[i*numberofdofspernode+0]];
		if(drag)k[i]=drag[doflist[i*numberofdofspernode+0]];
	}

	//Update material if necessary
	temperature_average_param=ParameterInputsRecover(inputs,"temperature_average");
	flow_law_param=ParameterInputsRecover(inputs,"B");
	if(temperature_average_param){
		for(i=0;i<3;i++) temperature_average+=temperature_average_param[doflist[i*numberofdofspernode+0]];
		temperature_average=temperature_average/3.0;
		B_average=Paterson(temperature_average);
		matice->SetB(B_average);
	}

	//Update material if B is provided.
	if(flow_law_param){
		for(i=0;i<3;i++) B_average+=flow_law_param[doflist[i*numberofdofspernode+0]];
		B_average=B_average/3.0;
		matice->SetB(B_average);
	}

}
		
void  Tria::GetDofList(int* doflist,int* pnumberofdofspernode){

	int i,j;
	int doflist_per_node[MAXDOFSPERNODE];
	int numberofdofspernode;
	
	for(i=0;i<3;i++){
		nodes[i]->GetDofList(&doflist_per_node[0],&numberofdofspernode);
		for(j=0;j<numberofdofspernode;j++){
			doflist[i*numberofdofspernode+j]=doflist_per_node[j];
		}
	}

	/*Assign output pointers:*/
	*pnumberofdofspernode=numberofdofspernode;

}


#undef __FUNCT__ 
#define __FUNCT__ "Tria::GetParameterValue"
void Tria::GetParameterValue(double* pp, double* plist, double* gauss_l1l2l3){
	
	/*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter value at gaussian 
	 * point specifie by gauss_l1l2l3: */
	
	/*nodal functions: */
	double l1l2l3[3];

	/*output: */
	double p;

	GetNodalFunctions(l1l2l3, gauss_l1l2l3);

	p=l1l2l3[0]*plist[0]+l1l2l3[1]*plist[1]+l1l2l3[2]*plist[2];

	/*Assign output pointers:*/
	*pp=p;
}


#undef __FUNCT__ 
#define __FUNCT__ "Tria::GetParameterDerivativeValue"
void Tria::GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, double* gauss_l1l2l3){
	 
	const int NDOF2=2;
	const int numgrids=3;
	/*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter derivative value at gaussian 
	 * point specified by gauss_l1l2l3:
	 *   dp/dx=plist[0]*dh1/dx+plist[1]*dh2/dx+plist[2]*dh3/dx
	 *   dp/dx=plist[0]*dh1/dx+plist[1]*dh2/dx+plist[2]*dh3/dx
	 *
	 * p is a vector of size 2x1 already allocated.
	 */
	
	double dh1dh2dh3_basic[NDOF2][numgrids]; //nodal derivative functions in basic coordinate system.

	/*Get dh1dh2dh3 in basic coordinate system: */
	GetNodalFunctionsDerivativesBasic(&dh1dh2dh3_basic[0][0],xyz_list, gauss_l1l2l3);

	*(p+0)=plist[0]*dh1dh2dh3_basic[0][0]+plist[1]*dh1dh2dh3_basic[0][1]+plist[2]*dh1dh2dh3_basic[0][2];
	*(p+1)=plist[0]*dh1dh2dh3_basic[1][0]+plist[1]*dh1dh2dh3_basic[1][1]+plist[2]*dh1dh2dh3_basic[1][2];

}


#undef __FUNCT__ 
#define __FUNCT__ "Tria::GetStrainRate"
void Tria::GetStrainRate(double* epsilon, double* velocity, double* xyz_list, double* gauss_l1l2l3){

	int i;
	const int NDOF2=2;
	const int numgrids=3;

	double B[3][NDOF2*numgrids];

	/*Get B matrix: */
	GetB(&B[0][0], xyz_list, gauss_l1l2l3);

	#ifdef _DEBUG_
	printf("B for grid1 : [ %lf   %lf  \n",B[0][0],B[0][1]);
	printf("              [ %lf   %lf  \n",B[1][0],B[1][1]);
	printf("              [ %lf   %lf  ]\n",B[2][0],B[2][1]);

	printf("B for grid2 : [ %lf   %lf  \n",B[0][2],B[0][3]);
	printf("              [ %lf   %lf  \n",B[1][2],B[1][3]);
	printf("              [ %lf   %lf  ]\n",B[2][2],B[2][3]);

	printf("B for grid3 : [ %lf   %lf  \n",B[0][4],B[0][5]);
	printf("              [ %lf   %lf  \n",B[1][4],B[1][5]);
	printf("              [ %lf   %lf  ]\n",B[2][4],B[2][5]);
		
	for (i=0;i<numgrids;i++){
		printf("Velocity for grid %i %lf %lf\n",i,*(vxvy_list+2*i+0),*(vxvy_list+2*i+1));
	}
	#endif

	/*Multiply B by velocity, to get strain rate: */
	MatrixMultiply( &B[0][0],3,NDOF2*numgrids,0,
			              velocity,NDOF2*numgrids,1,0,
						  epsilon,0);

}

#undef __FUNCT__ 
#define __FUNCT__ "Tria::GetJacobianDeterminant" 
void Tria::GetJacobianDeterminant(double*  Jdet, double* xyz_list,double* gauss_l1l2l3){

	/*The Jacobian determinant is constant over the element, discard the gaussian points. 
	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/

	double x1,x2,x3,y1,y2,y3;
	
	x1=*(xyz_list+3*0+0);
	y1=*(xyz_list+3*0+1);
	x2=*(xyz_list+3*1+0);
	y2=*(xyz_list+3*1+1);
	x3=*(xyz_list+3*2+0);
	y3=*(xyz_list+3*2+1);


	*Jdet=sqrt(3.0)/6.0*((x2-x1)*(y3-y1)-(x3-x1)*(y2-y1));

	if(Jdet<0){
		printf("%s%s\n",__FUNCT__," error message: negative jacobian determinant!");
	}
	
}

#undef __FUNCT__ 
#define __FUNCT__ "Tria::GetB"

void Tria::GetB(double* B, double* xyz_list, double* gauss_l1l2l3){

	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
	 * For grid i, Bi can be expressed in the basic coordinate system
	 * by: 
	 *       Bi_basic=[ dh/dx    0    ]
	 *                [   0    dh/dy  ]
	 *                [ 1/2*dh/dy  1/2*dh/dx  ]
	 * where h is the interpolation function for grid i.
	 *
	 * We assume B has been allocated already, of size: 3x(NDOF2*numgrids)
	 */
	
	int i;
	const int NDOF2=2;
	const int numgrids=3;

	double dh1dh2dh3_basic[NDOF2][numgrids];


	/*Get dh1dh2dh3 in basic coordinate system: */
	GetNodalFunctionsDerivativesBasic(&dh1dh2dh3_basic[0][0],xyz_list, gauss_l1l2l3);

	#ifdef _DEBUG_ 
	for (i=0;i<3;i++){
		printf("Node %i  dh/dx=%lf dh/dy=%lf \n",i,dh1dh2dh3_basic[0][i],dh1dh2dh3_basic[1][i]);
	}
	#endif

	/*Build B: */
	for (i=0;i<numgrids;i++){
		*(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh2dh3_basic[0][i]; //B[0][NDOF2*i]=dh1dh2dh3_basic[0][i];
		*(B+NDOF2*numgrids*0+NDOF2*i+1)=0;
		*(B+NDOF2*numgrids*1+NDOF2*i)=0;
		*(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh2dh3_basic[1][i];
		*(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh2dh3_basic[1][i]; 
		*(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh2dh3_basic[0][i]; 
	}
}

#undef __FUNCT__ 
#define __FUNCT__ "Tria::GetBPrime"

void Tria::GetBPrime(double* Bprime, double* xyz_list, double* gauss_l1l2l3){

	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
	 * For grid i, Bi' can be expressed in the basic coordinate system
	 * by: 
	 *       Bi_prime__basic=[ 2*dh/dx dh/dy ]
	 *                       [ dh/dx  2*dh/dy]
	 *                       [dh/dy dh/dx]
	 * where h is the interpolation function for grid i.
	 *
	 * We assume B' has been allocated already, of size: 3x(NDOF2*numgrids)
	 */
	
	int i;
	const int NDOF2=2;
	const int numgrids=3;

	/*Same thing in the basic coordinate system: */
	double dh1dh2dh3_basic[NDOF2][numgrids];


	/*Get dh1dh2dh3 in basic coordinates system : */
	GetNodalFunctionsDerivativesBasic(&dh1dh2dh3_basic[0][0],xyz_list,gauss_l1l2l3);

	/*Build B': */
	for (i=0;i<numgrids;i++){
		*(Bprime+NDOF2*numgrids*0+NDOF2*i)=2*dh1dh2dh3_basic[0][i]; 
		*(Bprime+NDOF2*numgrids*0+NDOF2*i+1)=dh1dh2dh3_basic[1][i]; 
		*(Bprime+NDOF2*numgrids*1+NDOF2*i)=dh1dh2dh3_basic[0][i]; 
		*(Bprime+NDOF2*numgrids*1+NDOF2*i+1)=2*dh1dh2dh3_basic[1][i]; 
		*(Bprime+NDOF2*numgrids*2+NDOF2*i)=dh1dh2dh3_basic[1][i]; 
		*(Bprime+NDOF2*numgrids*2+NDOF2*i+1)=dh1dh2dh3_basic[0][i]; 
	}
}

#undef __FUNCT__ 
#define __FUNCT__ "Tria::GetL"

void Tria::GetL(double* L, double* xyz_list, double* gauss_l1l2l3,int numdof){

	/*Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof. 
	 * For grid i, Li can be expressed in the basic coordinate system
	 * by: 
	 *       numdof=1: 
	 *       Li_basic=h;
	 *       numdof=2:
	 *       Li_basic=[ h    0    ]
	 *                [   0   h  ]
	 * where h is the interpolation function for grid i.
	 *
	 * We assume L has been allocated already, of size: numgrids (numdof=1), or numdofx(numdof*numgrids) (numdof=2)
	 */
	
	int i;
	const int NDOF2=2;
	const int numgrids=3;

	double l1l2l3[3];


	/*Get l1l2l3 in basic coordinate system: */
	GetNodalFunctions(l1l2l3, gauss_l1l2l3);

	#ifdef _DELUG_ 
	for (i=0;i<3;i++){
		printf("Node %i  h=%lf \n",i,l1l2l3[i]);
	}
	#endif

	/*Build L: */
	if(numdof==1){
		for (i=0;i<numgrids;i++){
			L[i]=l1l2l3[i]; 
		}
	}
	else{
		for (i=0;i<numgrids;i++){
			*(L+numdof*numgrids*0+numdof*i)=l1l2l3[i]; //L[0][NDOF2*i]=dh1dh2dh3_basic[0][i];
			*(L+numdof*numgrids*0+numdof*i+1)=0;
			*(L+numdof*numgrids*1+numdof*i)=0;
			*(L+numdof*numgrids*1+numdof*i+1)=l1l2l3[i];
		}
	}
}

#undef __FUNCT__ 
#define __FUNCT__ "Tria::GetNodalFunctions"
void Tria::GetNodalFunctions(double* l1l2l3, double* gauss_l1l2l3){
	
	/*This routine returns the values of the nodal functions  at the gaussian point.*/

	/*First nodal function: */
	l1l2l3[0]=gauss_l1l2l3[0];

	/*Second nodal function: */
	l1l2l3[1]=gauss_l1l2l3[1];

	/*Third nodal function: */
	l1l2l3[2]=gauss_l1l2l3[2];

}

#undef __FUNCT__ 
#define __FUNCT__ "Tria::GetNodalFunctionsDerivativesBasic"
void Tria::GetNodalFunctionsDerivativesBasic(double* dh1dh2dh3_basic,double* xyz_list, double* gauss_l1l2l3){
	
	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
	 * basic coordinate system: */

	int i;
	const int NDOF2=2;
	const int numgrids=3;

	double dh1dh2dh3_param[NDOF2][numgrids];
	double Jinv[NDOF2][NDOF2];


	/*Get derivative values with respect to parametric coordinate system: */
	GetNodalFunctionsDerivativesParams(&dh1dh2dh3_param[0][0], gauss_l1l2l3); 

	/*Get Jacobian invert: */
	GetJacobianInvert(&Jinv[0][0], xyz_list, gauss_l1l2l3);

	/*Build dh1dh2dh3_basic: 
	 *
	 * [dhi/dx]= Jinv*[dhi/dr]
	 * [dhi/dy]       [dhi/ds]
	 */

	for (i=0;i<numgrids;i++){
		*(dh1dh2dh3_basic+numgrids*0+i)=Jinv[0][0]*dh1dh2dh3_param[0][i]+Jinv[0][1]*dh1dh2dh3_param[1][i];
		*(dh1dh2dh3_basic+numgrids*1+i)=Jinv[1][0]*dh1dh2dh3_param[0][i]+Jinv[1][1]*dh1dh2dh3_param[1][i];
	}

}

#undef __FUNCT__ 
#define __FUNCT__ "Tria::GetNodalFunctionsDerivativesParams"
void Tria::GetNodalFunctionsDerivativesParams(double* dl1dl2dl3,double* gauss_l1l2l3){
	
	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
	 * natural coordinate system) at the gaussian point. */

	const int NDOF2=2;
	const int numgrids=3;

	double sqrt3=sqrt(3.0);

	/*First nodal function: */
	*(dl1dl2dl3+numgrids*0+0)=-1.0/2.0; 
	*(dl1dl2dl3+numgrids*1+0)=-1.0/(2.0*sqrt3);

	/*Second nodal function: */
	*(dl1dl2dl3+numgrids*0+1)=1.0/2.0;
	*(dl1dl2dl3+numgrids*1+1)=-1.0/(2.0*sqrt3);

	/*Third nodal function: */
	*(dl1dl2dl3+numgrids*0+2)=0;
	*(dl1dl2dl3+numgrids*1+2)=1.0/sqrt3;

}

#undef __FUNCT__ 
#define __FUNCT__ "Tria::GetJacobianInvert"
void Tria::GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss_l1l2l3){

	double Jdet;
	const int NDOF2=2;
	const int numgrids=3;
	
	/*Call Jacobian routine to get the jacobian:*/
	GetJacobian(Jinv, xyz_list, gauss_l1l2l3);

	/*Invert Jacobian matrix: */
	MatrixInverse(Jinv,NDOF2,NDOF2,NULL,0,&Jdet);
		
}

#undef __FUNCT__ 
#define __FUNCT__ "Tria::GetJacobian"
void Tria::GetJacobian(double* J, double* xyz_list,double* gauss_l1l2l3){

	/*The Jacobian is constant over the element, discard the gaussian points. 
	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/

	const int NDOF2=2;
	const int numgrids=3;
	double x1,y1,x2,y2,x3,y3;
	double sqrt3=sqrt(3.0);
	
	x1=*(xyz_list+numgrids*0+0);
	y1=*(xyz_list+numgrids*0+1);
	x2=*(xyz_list+numgrids*1+0);
	y2=*(xyz_list+numgrids*1+1);
	x3=*(xyz_list+numgrids*2+0);
	y3=*(xyz_list+numgrids*2+1);


	*(J+NDOF2*0+0)=1.0/2.0*(x2-x1);
	*(J+NDOF2*1+0)=sqrt3/6.0*(2*x3-x1-x2);
	*(J+NDOF2*0+1)=1.0/2.0*(y2-y1);
	*(J+NDOF2*1+1)=sqrt3/6.0*(2*y3-y1-y2);
}

Matpar* Tria::GetMatPar(){
	return matpar;
}

int   Tria::GetShelf(){
	return shelf;
}

		
void  Tria::GetNodes(Node** pnodes){
	int i;
	for(i=0;i<3;i++){
		pnodes[i]=nodes[i];
	}
}
		

int Tria::GetOnBed(){
	return onbed;
}

void          Tria::GetThicknessList(double* thickness_list){

	int i;
	for(i=0;i<3;i++)thickness_list[i]=h[i];
}
void          Tria::GetBedList(double* bed_list){
	
	int i;
	for(i=0;i<3;i++)bed_list[i]=b[i];

}
  

Object* Tria::copy() {
	return new Tria(*this); 
}


#undef __FUNCT__ 
#define __FUNCT__ "Tria::Du"
void Tria::Du(Vec du_g,double* velocity,double* obs_velocity,ParameterInputs* inputs,int analysis_type){

	int i;
	
	/* node data: */
	const int    numgrids=3;
	const int    numdofs=2*numgrids;
	const int    NDOF2=2;
	double       xyz_list[numgrids][3];
	int          doflist[numdofs];
	int          numberofdofspernode;

	/* grid data: */
	double vx_list[numgrids];
	double vy_list[numgrids];
	double obs_vx_list[numgrids];
	double obs_vy_list[numgrids];

	/* gaussian points: */
	int     num_gauss,ig;
	double* first_gauss_area_coord  =  NULL;
	double* second_gauss_area_coord =  NULL;
	double* third_gauss_area_coord  =  NULL;
	double* gauss_weights           =  NULL;
	double  gauss_weight;
	double  gauss_l1l2l3[3];

	/* parameters: */
	double  velocity_x,velocity_y,velocity_mag;
	double  obs_velocity_x,obs_velocity_y,obs_velocity_mag;

	/*element vector : */
	double  due_g[numdofs];
	double  due_g_gaussian[numdofs];

	/* Jacobian: */
	double Jdet;

	/*nodal functions: */
	double l1l2l3[3];

	/*relative and algorithmic fitting: */
	double meanvel=0;
	double epsvel=0;
	double scalex=1;
	double scaley=1;
	double scale=0;
	double* fit_param=NULL;
	double  fit=-1;

	/* Get node coordinates and dof list: */
	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
	GetDofList(&doflist[0],&numberofdofspernode);

	/* Set due_g to 0: */
	for(i=0;i<numdofs;i++) due_g[i]=0.0;

	/* Recover input data: */
	fit_param=ParameterInputsRecover(inputs,"fit");
	if(fit_param){
		fit=*fit_param;
	}

	/*Initialize velocities: */
	for(i=0;i<numgrids;i++){
		vx_list[i]=velocity[doflist[i*numberofdofspernode+0]];
		vy_list[i]=velocity[doflist[i*numberofdofspernode+1]];
		obs_vx_list[i]=obs_velocity[doflist[i*numberofdofspernode+0]];
		obs_vy_list[i]=obs_velocity[doflist[i*numberofdofspernode+1]];

	}

	
	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
	GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);

	#ifdef _DEBUGELEMENTS_
	if(my_rank==RANK && id==ELID){ 
		printf("   gaussian points: \n");
		for(i=0;i<num_gauss;i++){
			printf("    %g %g %g : %g\n",first_gauss_area_coord[i],second_gauss_area_coord[i],third_gauss_area_coord[i],gauss_weights[i]);
		}
	}
	#endif
	
	/* Start  looping on the number of gaussian points: */
	for (ig=0; ig<num_gauss; ig++){
		/*Pick up the gaussian point: */
		gauss_weight=*(gauss_weights+ig);
		gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 
		gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
		gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);

		
		/*Compute velocities at gaussian point: */
		GetParameterValue(&velocity_x, &vx_list[0],gauss_l1l2l3);
		GetParameterValue(&velocity_y, &vy_list[0],gauss_l1l2l3);
		#ifdef _DEBUG_ 
			printf("Velocity: %lf %lf\n", velocity_x,velocity_y);
		#endif
	
		/*Compute obs_velocities at gaussian point: */
		GetParameterValue(&obs_velocity_x, &obs_vx_list[0],gauss_l1l2l3);
		GetParameterValue(&obs_velocity_y, &obs_vy_list[0],gauss_l1l2l3);
		#ifdef _DEBUG_ 
			printf("Observed velocity: %g %g\n", obs_velocity_x,obs_velocity_y);
		#endif

		/* Get Jacobian determinant: */
		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
		#ifdef _DEBUG_ 
		printf("Element id %i Jacobian determinant: %lf\n",TriaElementGetID(this),Jdet);
		#endif
		
		/* Get nodal functions value at gaussian point:*/
		GetNodalFunctions(l1l2l3, gauss_l1l2l3);

		/*Build due_g_gaussian vector: we have three cases here, according to which type of misfit we are using. */
		if(fit==0){
			/*We are using an absolute misfit: */
			for (i=0;i<numgrids;i++){
				due_g_gaussian[i*NDOF2+0]=(velocity_x-obs_velocity_x)*Jdet*gauss_weight*l1l2l3[i]; 
				due_g_gaussian[i*NDOF2+1]=(velocity_y-obs_velocity_y)*Jdet*gauss_weight*l1l2l3[i]; 
			}
		}
		else if(fit==1){
			/*We are using a relative misfit: */
			scalex=pow(meanvel/(obs_velocity_x+epsvel),2);
			scaley=pow(meanvel/(obs_velocity_y+epsvel),2);
			if(obs_velocity_x==0)scalex=0;
			if(obs_velocity_y==0)scaley=0;
			for (i=0;i<numgrids;i++){
				due_g_gaussian[i*NDOF2+0]=scalex*(velocity_x-obs_velocity_x)*Jdet*gauss_weight*l1l2l3[i]; 
				due_g_gaussian[i*NDOF2+1]=scaley*(velocity_y-obs_velocity_y)*Jdet*gauss_weight*l1l2l3[i]; 
			}
		}
		else if(fit==2){

			/*We are using a logarithmic misfit: */
			velocity_mag=sqrt(pow(velocity_x,2)+pow(velocity_y,2))+epsvel; //epsvel to avoid velocity being nil.
			obs_velocity_mag=sqrt(pow(obs_velocity_x,2)+pow(obs_velocity_y,2))+epsvel; //epsvel to avoid observed velocity being nil.
			scale=8*pow(meanvel,2)/pow(velocity_mag,2)*log(velocity_mag/obs_velocity_mag);
			for (i=0;i<numgrids;i++){
				due_g_gaussian[i*NDOF2+0]=scale*velocity_x*Jdet*gauss_weight*l1l2l3[i]; 
				due_g_gaussian[i*NDOF2+1]=scale*velocity_y*Jdet*gauss_weight*l1l2l3[i]; 
			}
		}
		else{
			/*Not supported yet! : */
			throw ErrorException(__FUNCT__,exprintf("%s%g","unsupported type of fit: ",fit));
		}
		
		
		/*Add due_g_gaussian vector to due_g: */
		for( i=0; i<numdofs; i++){
			due_g[i]+=due_g_gaussian[i];
		}
	}
	
	
	/*Add due_g to global vector du_g: */
	VecSetValues(du_g,numdofs,doflist,(const double*)due_g,ADD_VALUES);
	
	cleanup_and_return: 
	xfree((void**)&first_gauss_area_coord);
	xfree((void**)&second_gauss_area_coord);
	xfree((void**)&third_gauss_area_coord);
	xfree((void**)&gauss_weights);

}
#undef __FUNCT__ 
#define __FUNCT__ "Tria::Gradj"
void  Tria::Gradj(Vec grad_g,double* velocity,double* adjoint,ParameterInputs* inputs,int analysis_type,char* control_type){

	if (strcmp(control_type,"drag")==0){
		GradjDrag( grad_g,velocity,adjoint,inputs,analysis_type);
	}
	else if (strcmp(control_type,"B")==0){
		GradjB( grad_g, velocity, adjoint, inputs,analysis_type);
	}
	else throw ErrorException(__FUNCT__,exprintf("%s%s","control type not supported yet: ",control_type));
}

#undef __FUNCT__ 
#define __FUNCT__ "Tria::GradjDrag"
void  Tria::GradjDrag(Vec grad_g,double* velocity,double* adjoint,ParameterInputs* inputs,int analysis_type){


	int i;

	/* node data: */
	const int    numgrids=3;
	const int    numdofs=2*numgrids;
	const int    NDOF2=2;
	double       xyz_list[numgrids][3];
	int          doflist[numdofs];
	int          numberofdofspernode;

	/* grid data: */
	double vx_list[numgrids];
	double vy_list[numgrids];
	double vxvy_list[numgrids][2];
	double adjx_list[numgrids];
	double adjy_list[numgrids];

	double* basal_drag_param=NULL;
	double K_list[numgrids];
	double drag;

	double* bed_param=NULL;
	double bed_list[numgrids];
	double* thickness_param=NULL;
	double thickness_list[numgrids];

	/* gaussian points: */
	int     num_gauss,ig;
	double* first_gauss_area_coord  =  NULL;
	double* second_gauss_area_coord =  NULL;
	double* third_gauss_area_coord  =  NULL;
	double* gauss_weights           =  NULL;
	double  gauss_weight;
	double  gauss_l1l2l3[3];

	/* parameters: */
	double  viscosity2;
	double  dvx[NDOF2];
	double  dvy[NDOF2]; 
	double  dadjx[NDOF2];
	double  dadjy[NDOF2];
	double  vx,vy;
	double  lambda,mu;
	double  bed,thickness,Neff;
	
	/*drag: */
	double  pcoeff,qcoeff;
	double alpha_complement_list[numgrids];
	double alpha_complement;


	/*element vector at the gaussian points: */
	double  grade_g[numgrids];
	double  grade_g_gaussian[numgrids];

	/* Jacobian: */
	double Jdet;

	/*nodal functions: */
	double l1l2l3[3];

	/* strain rate: */
	double epsilon[3]; /* epsilon=[exx,eyy,exy];*/

	/* Get node coordinates and dof list: */
	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
	GetDofList(&doflist[0],&numberofdofspernode);

	/* Set grade_g to 0: */
	for(i=0;i<numdofs;i++) grade_g[i]=0.0;

	/* recover input parameters: */
	basal_drag_param=ParameterInputsRecover(inputs,"drag");
	bed_param=ParameterInputsRecover(inputs,"bed");
	thickness_param=ParameterInputsRecover(inputs,"thickness");

	/*Initialize parameter lists: */
	for(i=0;i<numgrids;i++){
		vx_list[i]=velocity[doflist[i*numberofdofspernode+0]];
		vy_list[i]=velocity[doflist[i*numberofdofspernode+1]];
		vxvy_list[i][0]=vx_list[i];
		vxvy_list[i][1]=vy_list[i];
		adjx_list[i]=adjoint[doflist[i*numberofdofspernode+0]];
		adjy_list[i]=adjoint[doflist[i*numberofdofspernode+1]];
		if(basal_drag_param){
			K_list[i]=basal_drag_param[doflist[i*numberofdofspernode+0]];
		}
		else{
			K_list[i]=k[i];
		}
		if(bed_param){
			bed_list[i]=bed_param[doflist[i*numberofdofspernode+0]];
		}
		else{
			bed_list[i]=b[i];
		}
		if(thickness_param){
			thickness_list[i]=thickness_param[doflist[i*numberofdofspernode+0]];
		}
		else{
			thickness_list[i]=h[i];
		}
	}
	
	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
	GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 4);

	#ifdef _DEBUGELEMENTS_
	if(my_rank==RANK && id==ELID){ 
		printf("   gaussian points: \n");
		for(i=0;i<num_gauss;i++){
			printf("    %g %g %g : %g\n",first_gauss_area_coord[i],second_gauss_area_coord[i],third_gauss_area_coord[i],gauss_weights[i]);
		}
	}
	#endif
	
	/* Start  looping on the number of gaussian points: */
	for (ig=0; ig<num_gauss; ig++){
		/*Pick up the gaussian point: */
		gauss_weight=*(gauss_weights+ig);
		gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 
		gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
		gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);

		/*Build alpha_complement_list: */
		if (!shelf && (friction_type==2)){
		
			/*Allocate friction object: */
			Friction* friction=NewFriction();
			
			/*Initialize all fields: */
			friction->element_type=(char*)xmalloc((strlen("2d")+1)*sizeof(char));
			strcpy(friction->element_type,"2d");
			
			friction->gravity=matpar->GetG();
			friction->rho_ice=matpar->GetRhoIce();
			friction->rho_water=matpar->GetRhoWater();
			friction->K=&K_list[0];
			friction->bed=&bed_list[0];
			friction->thickness=&thickness_list[0];
			friction->velocities=&vxvy_list[0][0];
			friction->p=p;
			friction->q=q;

			if(friction->p!=1) throw ErrorException(__FUNCT__,"non-linear friction not supported yet in control methods!");

			/*Compute alpha complement: */
			FrictionGetAlphaComplement(&alpha_complement_list[0],friction);

			/*Erase friction object: */
			DeleteFriction(&friction);
		}
		else{
			alpha_complement_list[0]=0;
			alpha_complement_list[1]=0;
			alpha_complement_list[2]=0;
		}
	
		/*Recover alpha_complement and k: */
		GetParameterValue(&alpha_complement, &alpha_complement_list[0],gauss_l1l2l3);
		GetParameterValue(&drag, &K_list[0],gauss_l1l2l3);
		#ifdef _DEBUG_ 
			printf("Drag complement: %20.20lf Drag: %20.20lf\n",alpha_complement,drag);
		#endif

		/*recover lambda and mu: */
		GetParameterValue(&lambda, &adjx_list[0],gauss_l1l2l3);
		GetParameterValue(&mu, &adjy_list[0],gauss_l1l2l3);
		#ifdef _DEBUG_ 
			printf("Adjoint vector %20.20lf %20.20lf\n",lambda,mu);
		#endif
			
		/*recover vx and vy: */
		GetParameterValue(&vx, &vx_list[0],gauss_l1l2l3);
		GetParameterValue(&vy, &vy_list[0],gauss_l1l2l3);
		#ifdef _DEBUG_ 
			printf("Velocity vector %20.20lf %20.20lf\n",vx,vy);
		#endif

		/* Get Jacobian determinant: */
		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
		#ifdef _DEBUG_ 
		printf("Element id %i Jacobian determinant: %lf\n",TriaElementGetID(this),Jdet);
		#endif
		
		/* Get nodal functions value at gaussian point:*/
		GetNodalFunctions(l1l2l3, gauss_l1l2l3);

		/*Build gradje_g_gaussian vector: */
		for (i=0;i<numgrids;i++){
			grade_g_gaussian[i]=-2*drag*alpha_complement*( (lambda*vx+mu*vy))*Jdet*gauss_weight*l1l2l3[i]; 
		}
		
		/*Add gradje_g_gaussian vector to gradje_g: */
		for( i=0; i<numdofs; i++)grade_g[i]+=grade_g_gaussian[i];
	}

	/*Add grade_g to global vector grad_g: */
	VecSetValues(grad_g,numdofs,doflist,(const double*)grade_g,ADD_VALUES);
	
	cleanup_and_return: 
	xfree((void**)&first_gauss_area_coord);
	xfree((void**)&second_gauss_area_coord);
	xfree((void**)&third_gauss_area_coord);
	xfree((void**)&gauss_weights);

}

#undef __FUNCT__ 
#define __FUNCT__ "Tria::GradjB"
void  Tria::GradjB(Vec grad_g,double* velocity,double* adjoint,ParameterInputs* inputs,int analysis_type){

	int i;

	/* node data: */
	const int    numgrids=3;
	const int    numdofs=2*numgrids;
	const int    NDOF2=2;
	double       xyz_list[numgrids][3];
	int          doflist[numdofs];
	int          numberofdofspernode;

	/* grid data: */
	double vx_list[numgrids];
	double vy_list[numgrids];
	double vxvy_list[numgrids][2];
	double adjx_list[numgrids];
	double adjy_list[numgrids];

	double* thickness_param=NULL;
	double thickness_list[numgrids];

	/* gaussian points: */
	int     num_gauss,ig;
	double* first_gauss_area_coord  =  NULL;
	double* second_gauss_area_coord =  NULL;
	double* third_gauss_area_coord  =  NULL;
	double* gauss_weights           =  NULL;
	double  gauss_weight;
	double  gauss_l1l2l3[3];

	/*element vector at the gaussian points: */
	double  grade_g[numgrids];
	double  grade_g_gaussian[numgrids];

	/* Jacobian: */
	double Jdet;

	/*nodal functions: */
	double l1l2l3[3];

	/* strain rate: */
	double epsilon[3]; /* epsilon=[exx,eyy,exy];*/


	/* parameters: */
	double  viscosity2;
	double  dvx[NDOF2];
	double  dvy[NDOF2]; 
	double  dadjx[NDOF2];
	double  dadjy[NDOF2];
	double  vx,vy;
	double  lambda,mu;
	double  thickness;
	

	/* Get node coordinates and dof list: */
	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
	GetDofList(&doflist[0],&numberofdofspernode);

	/* Set grade_g to 0: */
	for(i=0;i<numdofs;i++) grade_g[i]=0.0;

	/* recover input parameters: */
	thickness_param=ParameterInputsRecover(inputs,"thickness");

	/*Initialize parameter lists: */
	for(i=0;i<numgrids;i++){
		vx_list[i]=velocity[doflist[i*numberofdofspernode+0]];
		vy_list[i]=velocity[doflist[i*numberofdofspernode+1]];
		vxvy_list[i][0]=vx_list[i];
		vxvy_list[i][1]=vy_list[i];
		adjx_list[i]=adjoint[doflist[i*numberofdofspernode+0]];
		adjy_list[i]=adjoint[doflist[i*numberofdofspernode+1]];
		if(thickness_param){
			thickness_list[i]=thickness_param[doflist[i*numberofdofspernode+0]];
		}
		else{
			thickness_list[i]=h[i];
		}
	}
	
	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
	GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 4);
	#ifdef _DEBUG_ 
	for (i=0;i<num_gauss;i++){
		printf("Gauss coord %i: %lf %lf %lf Weight: %lf\n",i,*(first_gauss_area_coord+i),*(second_gauss_area_coord+i),*(third_gauss_area_coord+i),*(gauss_weights+i));
	}
	#endif

	/* Start  looping on the number of gaussian points: */
	for (ig=0; ig<num_gauss; ig++){
		/*Pick up the gaussian point: */
		gauss_weight=*(gauss_weights+ig);
		gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 
		gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
		gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);

		/*Get thickness: */
		GetParameterValue(&thickness, &h[0],gauss_l1l2l3);
	
		/*Get strain rate, if velocity has been supplied: */
		GetStrainRate(&epsilon[0],&vxvy_list[0][0],&xyz_list[0][0],gauss_l1l2l3);
	
		/*Get viscosity2: */
		matice->GetViscosity2(&viscosity2, &epsilon[0]);
		
		/*Get dvx, dvy, dadjx and dadjx: */
		GetParameterDerivativeValue(&dvx[0], &vx_list[0],&xyz_list[0][0], gauss_l1l2l3);
		GetParameterDerivativeValue(&dvy[0], &vy_list[0],&xyz_list[0][0], gauss_l1l2l3);
		GetParameterDerivativeValue(&dadjx[0], &adjx_list[0],&xyz_list[0][0], gauss_l1l2l3);
		GetParameterDerivativeValue(&dadjy[0], &adjy_list[0],&xyz_list[0][0], gauss_l1l2l3);

		/* Get Jacobian determinant: */
		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
		
		/* Get nodal functions value at gaussian point:*/
		GetNodalFunctions(l1l2l3, gauss_l1l2l3);
		#ifdef _DEBUG_
			_printf_("viscosity2 %g thickness %g dvx [%g %g] dvy [%g %g]  dadjx [%g %g] dadjy[%g %g]\n",viscosity2,thickness,dvx[0],dvx[1],dvy[0],dvy[1],dadjx[0],dadjx[1],dadjy[0],dadjy[1]);
		#endif

		/*Build gradje_g_gaussian vector: */
		for (i=0;i<numgrids;i++){
			grade_g_gaussian[i]=viscosity2*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*l1l2l3[i]; 
		}

		/*Add grade_g_gaussian to grade_g: */
		for( i=0; i<numdofs;i++)grade_g[i]+=grade_g_gaussian[i];
	}
	
	/*Add grade_g to global vector grad_g: */
	VecSetValues(grad_g,numdofs,doflist,(const double*)grade_g,ADD_VALUES);
	
	cleanup_and_return: 
	xfree((void**)&first_gauss_area_coord);
	xfree((void**)&second_gauss_area_coord);
	xfree((void**)&third_gauss_area_coord);
	xfree((void**)&gauss_weights);
}


#undef __FUNCT__ 
#define __FUNCT__ "Tria::Misfit"
double Tria::Misfit(double* velocity,double* obs_velocity,ParameterInputs* inputs,int analysis_type){
	int i;
	
	/* output: */
	double Jelem=0;

	/* node data: */
	const int    numgrids=3;
	const int    numdofs=2*numgrids;
	const int    NDOF2=2;
	double       xyz_list[numgrids][3];
	int          doflist[numdofs];
	int          numberofdofspernode;

	/* grid data: */
	double vx_list[numgrids];
	double vy_list[numgrids];
	double obs_vx_list[numgrids];
	double obs_vy_list[numgrids];

	/* gaussian points: */
	int     num_gauss,ig;
	double* first_gauss_area_coord  =  NULL;
	double* second_gauss_area_coord =  NULL;
	double* third_gauss_area_coord  =  NULL;
	double* gauss_weights           =  NULL;
	double  gauss_weight;
	double  gauss_l1l2l3[3];

	/* parameters: */
	double  velocity_x,velocity_y,velocity_mag;
	double  obs_velocity_x,obs_velocity_y,obs_velocity_mag;

	/* Jacobian: */
	double Jdet;
	
	/*relative and logarithmic control method :*/
	double scalex=1;
	double scaley=1;
	double* fit_param=NULL;
	double  fit=-1;

	/* Get node coordinates and dof list: */
	GetElementNodeData( &xyz_list[0][0], nodes, numgrids);
	GetDofList(&doflist[0],&numberofdofspernode);
	
	/* Recover input data: */
	fit_param=ParameterInputsRecover(inputs,"fit");
	if(fit_param){
		fit=*fit_param;
	}

	/*Initialize velocities: */
	for(i=0;i<numgrids;i++){
		vx_list[i]=velocity[doflist[i*numberofdofspernode+0]];
		vy_list[i]=velocity[doflist[i*numberofdofspernode+1]];
		obs_vx_list[i]=obs_velocity[doflist[i*numberofdofspernode+0]];
		obs_vy_list[i]=obs_velocity[doflist[i*numberofdofspernode+1]];

	}

	
	/* Get gaussian points and weights (make this a statically initialized list of points? fstd): */
	GaussTria( &num_gauss, &first_gauss_area_coord, &second_gauss_area_coord, &third_gauss_area_coord, &gauss_weights, 2);
	
	#ifdef _DEBUG_ 
	for (i=0;i<num_gauss;i++){
		printf("Gauss coord %i: %lf %lf %lf Weight: %lf\n",i,*(first_gauss_area_coord+i),*(second_gauss_area_coord+i),*(third_gauss_area_coord+i),*(gauss_weights+i));
	}
	#endif

	/* Start  looping on the number of gaussian points: */
	for (ig=0; ig<num_gauss; ig++){
		/*Pick up the gaussian point: */
		gauss_weight=*(gauss_weights+ig);
		gauss_l1l2l3[0]=*(first_gauss_area_coord+ig); 
		gauss_l1l2l3[1]=*(second_gauss_area_coord+ig);
		gauss_l1l2l3[2]=*(third_gauss_area_coord+ig);

	
		/*Compute velocities at gaussian point: */
		GetParameterValue(&velocity_x, &vx_list[0],gauss_l1l2l3);
		GetParameterValue(&velocity_y, &vy_list[0],gauss_l1l2l3);
		#ifdef _DEBUG_ 
			printf("Velocity: %g %g\n", velocity_x,velocity_y);
		#endif
	
		/*Compute obs_velocities at gaussian point: */
		GetParameterValue(&obs_velocity_x, &obs_vx_list[0],gauss_l1l2l3);
		GetParameterValue(&obs_velocity_y, &obs_vy_list[0],gauss_l1l2l3);
		#ifdef _DEBUG_ 
			printf("Observed velocity: %g %g\n", obs_velocity_x,obs_velocity_y);
		#endif

		/* Get Jacobian determinant: */
		GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss_l1l2l3);
		#ifdef _DEBUG_ 
		printf("Element id %i Jacobian determinant: %lf\n",TriaElementGetID(this),Jdet);
		#endif
		
		/*Differents misfits are allowed: */
		if(fit==0){
			/*Absolute misfit: */
			Jelem+=.5*(pow((velocity_x-obs_velocity_x),2)+pow((velocity_y-obs_velocity_y),2))*Jdet*gauss_weight;
		}
		else if(fit==1){
			/*Relative misfit: */
			scalex=pow(meanvel/(obs_velocity_x+epsvel),2);
			scaley=pow(meanvel/(obs_velocity_y+epsvel),2);
			if(obs_velocity_x==0)scalex=0;
			if(obs_velocity_y==0)scaley=0;
			Jelem+=.5*(scalex*pow((velocity_x-obs_velocity_x),2)+scaley*pow((velocity_y-obs_velocity_y),2))*Jdet*gauss_weight;
		}	
		else if(fit==2){
			/*Logarithmic misfit: */
			velocity_mag=sqrt(pow(velocity_x,2)+pow(velocity_y,2))+epsvel; //epsvel to avoid velocity being nil.
			obs_velocity_mag=sqrt(pow(obs_velocity_x,2)+pow(obs_velocity_y,2))+epsvel; //epsvel to avoid observed velocity being nil.
			Jelem+=4*pow(meanvel,2)*pow(log(velocity_mag/obs_velocity_mag),2);
		}
		else throw ErrorException(__FUNCT__,exprintf("%s%i%s","fit type",fit," not supported yet!"));

	}
	cleanup_and_return: 
	xfree((void**)&first_gauss_area_coord);
	xfree((void**)&second_gauss_area_coord);
	xfree((void**)&third_gauss_area_coord);
	xfree((void**)&gauss_weights);

	/*Return: */
	return Jelem;
}
