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

/*Constructors/destructor/copy*/
Element::Element(){/*{{{*/
	this->id  = -1;
	this->sid = -1;
	this->inputs     = NULL;
	this->nodes      = NULL;
	this->vertices   = NULL;
	this->material   = NULL;
	this->matpar     = NULL;
	this->inputs     = NULL;
	this->parameters = NULL;
}/*}}}*/
Element::~Element(){/*{{{*/
	delete inputs;
}
/*}}}*/

/*Other*/
void Element::AddInput(Input* input_in){/*{{{*/

	/*Call inputs method*/
	_assert_(this->inputs);
	this->inputs->AddInput(input_in);
}/*}}}*/
/* bool Element::AllActive(void){/\*{{{*\/ */

/* 	int numnodes = this->GetNumberOfNodes(); */
/* 	for(int i=0;i<numnodes;i++){ */
/* 		if(!this->nodes[i]->IsActive()) return false; */
/* 	} */
/* 	return true; */
/* }/\*}}}*\/ */
/* bool Element::AnyActive(void){/\*{{{*\/ */

/* 	int numnodes = this->GetNumberOfNodes(); */
/* 	for(int i=0;i<numnodes;i++){ */
/* 		if(this->nodes[i]->IsActive()) return true; */
/* 	} */
/* 	return false; */
/* }/\*}}}*\/ */
void Element::CoordinateSystemTransform(IssmDouble** ptransform,Node** nodes_list,int numnodes,int* cs_array){/*{{{*/

	int         i,counter;
	int         numdofs   = 0;
	IssmDouble  norm;
	IssmDouble *transform = NULL;
	IssmDouble  coord_system[3][3];

	/*Some checks in debugging mode*/
	_assert_(numnodes && nodes_list);

	/*Get total number of dofs*/
	for(i=0;i<numnodes;i++){
		switch(cs_array[i]){
			case PressureEnum: numdofs+=1; break;
			case XYEnum:       numdofs+=2; break;
			case XYZEnum:      numdofs+=3; break;
			default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
		}
	}

	/*Allocate and initialize transform matrix*/
	transform=xNew<IssmDouble>(numdofs*numdofs);
	for(i=0;i<numdofs*numdofs;i++) transform[i]=0.0;

	/*Create transform matrix for all nodes (x,y for 2d and x,y,z for 3d). It is a block matrix
	 *for 3 nodes:

	 *     | T1 0  0 |
	 * Q = | 0  T2 0 |
	 *     | 0  0  T3|
	 *
	 * Where T1 is the transform matrix for node 1. It is a simple copy of the coordinate system
	 * associated to this node*/
	counter=0;
	for(i=0;i<numnodes;i++){
		nodes_list[i]->GetCoordinateSystem(&coord_system[0][0]);
		switch(cs_array[i]){
			case PressureEnum:
				/*DO NOT change anything*/
				transform[(numdofs)*(counter) + counter] = 1.;
				counter+=1;
				break;
			case XYEnum:
				/*We remove the z component, we need to renormalize x and y: x=[x1 x2 0] y=[-x2 x1 0]*/
				norm = sqrt( coord_system[0][0]*coord_system[0][0] + coord_system[1][0]*coord_system[1][0]); _assert_(norm>1.e-4);
				transform[(numdofs)*(counter+0) + counter+0] =   coord_system[0][0]/norm;
				transform[(numdofs)*(counter+0) + counter+1] = - coord_system[1][0]/norm;
				transform[(numdofs)*(counter+1) + counter+0] =   coord_system[1][0]/norm;
				transform[(numdofs)*(counter+1) + counter+1] =   coord_system[0][0]/norm;
				counter+=2;
				break;
			case XYZEnum:
				/*The 3 coordinates are changed (x,y,z)*/
				transform[(numdofs)*(counter+0) + counter+0] = coord_system[0][0];
				transform[(numdofs)*(counter+0) + counter+1] = coord_system[0][1];
				transform[(numdofs)*(counter+0) + counter+2] = coord_system[0][2];
				transform[(numdofs)*(counter+1) + counter+0] = coord_system[1][0];
				transform[(numdofs)*(counter+1) + counter+1] = coord_system[1][1];
				transform[(numdofs)*(counter+1) + counter+2] = coord_system[1][2];
				transform[(numdofs)*(counter+2) + counter+0] = coord_system[2][0];
				transform[(numdofs)*(counter+2) + counter+1] = coord_system[2][1];
				transform[(numdofs)*(counter+2) + counter+2] = coord_system[2][2];
				counter+=3;
				break;
			default:
				_error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
		}
	}

	/*Assign output pointer*/
	*ptransform=transform;
}
/*}}}*/
void Element::DeleteMaterials(void){/*{{{*/
	delete this->material;
}/*}}}*/
void Element::DeepEcho(void){/*{{{*/

	_printf_(EnumToStringx(this->ObjectEnum())<<" element:\n");
	_printf_("   id : "<<this->id <<"\n");
	_printf_("   sid: "<<this->sid<<"\n");
	if(vertices){
		int numvertices = this->GetNumberOfVertices();
		for(int i=0;i<numvertices;i++) vertices[i]->Echo();
	}
	else _printf_("vertices = NULL\n");

	if(nodes){
		int numnodes = this->GetNumberOfNodes();
		for(int i=0;i<numnodes;i++) nodes[i]->DeepEcho();
	}
	else _printf_("nodes = NULL\n");

	if (material) material->DeepEcho();
	else _printf_("material = 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");

	return;
}
/*}}}*/
void Element::Echo(void){/*{{{*/
	_printf_(EnumToStringx(this->ObjectEnum())<<" element:\n");
	_printf_("   id : "<<this->id <<"\n");
	_printf_("   sid: "<<this->sid<<"\n");
	if(vertices){
		int numvertices = this->GetNumberOfVertices();
		for(int i=0;i<numvertices;i++) vertices[i]->Echo();
	}
	else _printf_("vertices = NULL\n");

	if(nodes){
		int numnodes = this->GetNumberOfNodes();
		for(int i=0;i<numnodes;i++) nodes[i]->Echo();
	}
	else _printf_("nodes = NULL\n");

	if (material) material->Echo();
	else _printf_("material = 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");
}
/*}}}*/
IssmDouble Element::Divergence(void){/*{{{*/
	/*Compute element divergence*/

	/*Intermediaries*/
	IssmDouble Jdet;
	IssmDouble divergence=0.;
	IssmDouble dvx[3],dvy[3],dvz[3];
	IssmDouble *xyz_list = NULL;

	/*Get inputs and parameters*/
	Input* vx_input        = this->GetInput(VxEnum); _assert_(vx_input);
	Input* vy_input        = this->GetInput(VyEnum); _assert_(vy_input);
	Input* vz_input        = this->GetInput(VzEnum); _assert_(vz_input);
	this->GetVerticesCoordinates(&xyz_list);

	Gauss* gauss=this->NewGaussBase(5);
	for(int ig=gauss->begin();ig<gauss->end();ig++){
		gauss->GaussPoint(ig);
		this->JacobianDeterminant(&Jdet,xyz_list,gauss);

		/*Get strain rate assuming that epsilon has been allocated*/
		vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
		vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
		vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss);

		divergence += (dvx[0]+dvy[1]+dvz[2])*gauss->weight*Jdet;
	}

	/*Clean up and return*/
	delete gauss;
	return divergence;
}/*}}}*/
void Element::ThermalToEnthalpy(IssmDouble* penthalpy,IssmDouble temperature,IssmDouble waterfraction,IssmDouble pressure){/*{{{*/
	matpar->ThermalToEnthalpy(penthalpy,temperature,waterfraction,pressure);
}/*}}}*/
void Element::EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){/*{{{*/
	matpar->EnthalpyToThermal(ptemperature,pwaterfraction,enthalpy,pressure);
}/*}}}*/
IssmDouble Element::EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){/*{{{*/
	return matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure);
}/*}}}*/
IssmDouble Element::EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){/*{{{*/
	return matpar->GetEnthalpyDiffusionParameterVolume(numvertices,enthalpy,pressure);
}/*}}}*/
void Element::FindParam(bool* pvalue,int paramenum){/*{{{*/
	this->parameters->FindParam(pvalue,paramenum);
}/*}}}*/
void Element::FindParam(int* pvalue,int paramenum){/*{{{*/
	this->parameters->FindParam(pvalue,paramenum);
}/*}}}*/
void Element::FindParam(IssmDouble* pvalue,int paramenum){/*{{{*/
	this->parameters->FindParam(pvalue,paramenum);
}/*}}}*/
void Element::FindParam(int** pvalues,int* psize,int paramenum){/*{{{*/
	this->parameters->FindParam(pvalues,psize,paramenum);
}/*}}}*/
void Element::GetDofList(int** pdoflist,int approximation_enum,int setenum){/*{{{*/

	/*Fetch number of nodes and dof for this finite element*/
	int numnodes = this->GetNumberOfNodes();

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

	/*Allocate output*/
	int* doflist=xNew<int>(numberofdofs);

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

	/*Assign output pointers:*/
	*pdoflist=doflist;
}
/*}}}*/
void Element::GetDofListVelocity(int** pdoflist,int setenum){/*{{{*/

	/*Fetch number of nodes and dof for this finite element*/
	int numnodes = this->NumberofNodesVelocity();

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

	/*Allocate output*/
	int* doflist=xNew<int>(numberofdofs);

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

	/*Assign output pointers:*/
	*pdoflist=doflist;
}
/*}}}*/
void Element::GetDofListPressure(int** pdoflist,int setenum){/*{{{*/

	/*Fetch number of nodes and dof for this finite element*/
	int vnumnodes = this->NumberofNodesVelocity();
	int pnumnodes = this->NumberofNodesPressure();

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

	/*Allocate output*/
	int* doflist=xNew<int>(numberofdofs);

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

	/*Assign output pointers:*/
	*pdoflist=doflist;
}
/*}}}*/
IssmDouble Element::GetMaterialParameter(int enum_in){/*{{{*/

	_assert_(this->matpar);
	switch(enum_in){ // FIXME: change this to material
		case MaterialsRheologyNEnum:
			return this->material->GetN();
		case MaterialsRheologyBEnum:
			return this->material->GetB();
		case MaterialsRheologyBbarEnum:
			return this->material->GetBbar();
		default:
			return this->matpar->GetMaterialParameter(enum_in);
	}
}/*}}}*/
void Element::GetPhi(IssmDouble* phi, IssmDouble*  epsilon, IssmDouble viscosity){/*{{{*/
	/*Compute deformational heating from epsilon and viscosity */

	IssmDouble epsilon_matrix[3][3];
	IssmDouble epsilon_eff;
	IssmDouble epsilon_sqr[3][3];

	/* Build epsilon matrix */
	epsilon_matrix[0][0]=epsilon[0];
	epsilon_matrix[1][0]=epsilon[3];
	epsilon_matrix[2][0]=epsilon[4];
	epsilon_matrix[0][1]=epsilon[3];
	epsilon_matrix[1][1]=epsilon[1];
	epsilon_matrix[2][1]=epsilon[5];
	epsilon_matrix[0][2]=epsilon[4];
	epsilon_matrix[1][2]=epsilon[5];
	epsilon_matrix[2][2]=epsilon[2];

	/* Effective value of epsilon_matrix */
	epsilon_sqr[0][0]=epsilon_matrix[0][0]*epsilon_matrix[0][0];
	epsilon_sqr[1][0]=epsilon_matrix[1][0]*epsilon_matrix[1][0];
	epsilon_sqr[2][0]=epsilon_matrix[2][0]*epsilon_matrix[2][0];
	epsilon_sqr[0][1]=epsilon_matrix[0][1]*epsilon_matrix[0][1];
	epsilon_sqr[1][1]=epsilon_matrix[1][1]*epsilon_matrix[1][1];
	epsilon_sqr[2][1]=epsilon_matrix[2][1]*epsilon_matrix[2][1];
	epsilon_sqr[0][2]=epsilon_matrix[0][2]*epsilon_matrix[0][2];
	epsilon_sqr[1][2]=epsilon_matrix[1][2]*epsilon_matrix[1][2];
	epsilon_sqr[2][2]=epsilon_matrix[2][2]*epsilon_matrix[2][2];
	epsilon_eff=1/sqrt(2.)*sqrt(epsilon_sqr[0][0]+epsilon_sqr[0][1]+ epsilon_sqr[0][2]+ epsilon_sqr[1][0]+ epsilon_sqr[1][1]+ epsilon_sqr[1][2]+ epsilon_sqr[2][0]+ epsilon_sqr[2][1]+ epsilon_sqr[2][2]);

	/*Phi = Tr(sigma * eps) 
	 *    = Tr(sigma'* eps)
	 *    = 2 * eps_eff * sigma'_eff
	 *    = 4 * mu * eps_eff ^2*/
	*phi=4.*epsilon_eff*epsilon_eff*viscosity;
}
/*}}}*/
Input* Element::GetInput(int inputenum){/*{{{*/
	return inputs->GetInput(inputenum);
}/*}}}*/
void Element::GetInputListOnVertices(IssmDouble* pvalue,int enumtype){/*{{{*/

	/*Recover input*/
	Input* input=this->GetInput(enumtype);
	if (!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element");

	/*Fetch number vertices for this element*/
	int numvertices = this->GetNumberOfVertices();

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

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

	/*clean-up*/
	delete gauss;
}
/*}}}*/
void Element::GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue){/*{{{*/

	/*Recover input*/
	Input* input=this->GetInput(enumtype);

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

	/*Fetch number vertices for this element*/
	int numvertices = this->GetNumberOfVertices();

	/* Start looping on the number of vertices: */
	if (input){
		Gauss* gauss=this->NewGauss();
		for (int iv=0;iv<numvertices;iv++){
			gauss->GaussVertex(iv);
			input->GetInputValue(&pvalue[iv],gauss);
		}
		delete gauss;
	}
	else{
		for(int iv=0;iv<numvertices;iv++) pvalue[iv]=defaultvalue;
	}
}
/*}}}*/
void Element::GetInputListOnNodes(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue){/*{{{*/

	_assert_(pvalue);

	Input *input    = this->GetInput(enumtype);
	int    numnodes = this->GetNumberOfNodes();

	/* Start looping on the number of vertices: */
	if(input){
		Gauss* gauss=this->NewGauss();
		for(int iv=0;iv<numnodes;iv++){
			gauss->GaussNode(this->FiniteElement(),iv);
			input->GetInputValue(&pvalue[iv],gauss);
		}
		delete gauss;
	}
	else{
		for(int iv=0;iv<numnodes;iv++) pvalue[iv]=defaultvalue;
	}
}
/*}}}*/
void Element::GetInputListOnNodes(IssmDouble* pvalue,int enumtype){/*{{{*/

	_assert_(pvalue);

	int    numnodes = this->GetNumberOfNodes();
	Input *input    = this->GetInput(enumtype);
	if(!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element");

	/* Start looping on the number of vertices: */
	Gauss* gauss=this->NewGauss();
	for(int iv=0;iv<numnodes;iv++){
		gauss->GaussNode(this->FiniteElement(),iv);
		input->GetInputValue(&pvalue[iv],gauss);
	}
	delete gauss;
}
/*}}}*/
void Element::GetInputListOnNodesVelocity(IssmDouble* pvalue,int enumtype){/*{{{*/

	_assert_(pvalue);

	int    numnodes = this->NumberofNodesVelocity();
	Input *input    = this->GetInput(enumtype);
	if(!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element");

	/* Start looping on the number of vertices: */
	Gauss* gauss=this->NewGauss();
	for(int iv=0;iv<numnodes;iv++){
		gauss->GaussNode(this->VelocityInterpolation(),iv);
		input->GetInputValue(&pvalue[iv],gauss);
	}
	delete gauss;
}
/*}}}*/
void Element::GetInputValue(bool* pvalue,int inputenum){/*{{{*/

	Input* input=inputs->GetInput(inputenum);
	if(!input) _error_("Input " << EnumToStringx(inputenum) << " not found in element");
	input->GetInputValue(pvalue);

}/*}}}*/
void Element::GetInputValue(int* pvalue,int inputenum){/*{{{*/

	Input* input=inputs->GetInput(inputenum);
	if(!input) _error_("Input " << EnumToStringx(inputenum) << " not found in element");
	input->GetInputValue(pvalue);

}/*}}}*/
void Element::GetInputValue(IssmDouble* pvalue,int inputenum){/*{{{*/

	Input* input=inputs->GetInput(inputenum);
	if(!input) _error_("Input " << EnumToStringx(inputenum) << " not found in element");
	input->GetInputValue(pvalue);

}/*}}}*/
void Element::GetInputValue(IssmDouble* pvalue,Gauss* gauss,int inputenum){/*{{{*/

	Input* input=inputs->GetInput(inputenum);
	if(!input) _error_("Input " << EnumToStringx(inputenum) << " not found in element");
	input->GetInputValue(pvalue,gauss);

}/*}}}*/
void Element::GetNodesSidList(int* sidlist){/*{{{*/

	_assert_(sidlist);
	_assert_(nodes);
	int numnodes = this->GetNumberOfNodes();
	for(int i=0;i<numnodes;i++){
		sidlist[i]=nodes[i]->Sid();
	}
}
/*}}}*/
void Element::GetNodesLidList(int* lidlist){/*{{{*/

	_assert_(lidlist);
	_assert_(nodes);
	int numnodes = this->GetNumberOfNodes();
	for(int i=0;i<numnodes;i++){
		lidlist[i]=nodes[i]->Lid();
	}
}
/*}}}*/
void Element::GetVerticesCoordinates(IssmDouble** pxyz_list){/*{{{*/

	int         numvertices = this->GetNumberOfVertices();
	IssmDouble* xyz_list    = xNew<IssmDouble>(numvertices*3);
	::GetVerticesCoordinates(xyz_list,this->vertices,numvertices);

	*pxyz_list = xyz_list;

}/*}}}*/
void Element::GetVerticesSidList(int* sidlist){/*{{{*/

	int numvertices = this->GetNumberOfVertices();
	for(int i=0;i<numvertices;i++) sidlist[i]=this->vertices[i]->Sid();
}
/*}}}*/
void Element::GetVerticesConnectivityList(int* connectivity){/*{{{*/

	int numvertices = this->GetNumberOfVertices();
	for(int i=0;i<numvertices;i++) connectivity[i]=this->vertices[i]->Connectivity();
}
/*}}}*/
bool Element::HasNodeOnBed(){/*{{{*/
	return (this->inputs->Max(MeshVertexonbedEnum)>0.);
}/*}}}*/
bool Element::HasNodeOnSurface(){/*{{{*/
	return (this->inputs->Max(MeshVertexonsurfaceEnum)>0.);
}/*}}}*/
int  Element::Id(){/*{{{*/

	return this->id;

}
/*}}}*/
void Element::InputChangeName(int original_enum,int new_enum){/*{{{*/
	this->inputs->ChangeEnum(original_enum,new_enum);
}
/*}}}*/
void Element::InputCreate(IssmDouble* vector,IoModel* iomodel,int M,int N,int vector_type,int vector_enum,int code){/*{{{*/

	/*Intermediaries*/
	int        i,t,row;
	IssmDouble time;

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

		int         numvertices = this->GetNumberOfVertices();
		int        *vertexids   = xNew<int>(numvertices);
		IssmDouble *values      = xNew<IssmDouble>(numvertices);

		/*Recover vertices ids needed to initialize inputs*/
		_assert_(iomodel->elements);
		for(i=0;i<numvertices;i++){ 
			vertexids[i]=reCast<int>(iomodel->elements[numvertices*this->Sid()+i]); //ids for vertices are in the elements array from Matlab
		}

		/*Are we in transient or static? */
		if(M==iomodel->numberofvertices){
			for(i=0;i<numvertices;i++) values[i]=vector[vertexids[i]-1];
			this->AddInput(vector_enum,values,P1Enum);
		}
		else if(M==iomodel->numberofvertices+1){
			/*create transient input: */
			IssmDouble* times = xNew<IssmDouble>(N);
			for(t=0;t<N;t++) times[t] = vector[(M-1)*N+t];
			TransientInput* transientinput=new TransientInput(vector_enum,times,N);
			for(t=0;t<N;t++){
				for(i=0;i<numvertices;i++) values[i]=vector[N*(vertexids[i]-1)+t];
				switch(this->ObjectEnum()){
					case TriaEnum:  transientinput->AddTimeInput(new TriaInput( vector_enum,values,P1Enum)); break;
					case PentaEnum: transientinput->AddTimeInput(new PentaInput(vector_enum,values,P1Enum)); break;
					case TetraEnum: transientinput->AddTimeInput(new TetraInput(vector_enum,values,P1Enum)); break;
					default: _error_("Not implemented yet");
				}
			}
			this->inputs->AddInput(transientinput);
			xDelete<IssmDouble>(times);
		}
		else _error_("nodal vector is either numberofvertices or numberofvertices+1 long. Field provided (" << EnumToStringx(vector_enum) << ") is " << M << " long");

		xDelete<IssmDouble>(values);
		xDelete<int>(vertexids);
	}
	else if(vector_type==2){ //element vector
		/*Are we in transient or static? */
		if(M==iomodel->numberofelements){
			if (code==5){ //boolean
				this->inputs->AddInput(new BoolInput(vector_enum,reCast<bool>(vector[this->Sid()])));
			}
			else if (code==6){ //integer
				this->inputs->AddInput(new IntInput(vector_enum,reCast<int>(vector[this->Sid()])));
			}
			else if (code==7){ //IssmDouble
				this->inputs->AddInput(new DoubleInput(vector_enum,vector[this->Sid()]));
			}
			else _error_("could not recognize nature of vector from code " << code);
		}
		else {
			_error_("transient element inputs not supported yet!");
		}
	}
	else{
		_error_("Cannot add input for vector type " << vector_type << " (not supported)");
	}
}/*}}}*/
void Element::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));
}
/*}}}*/
void Element::InputUpdateFromConstant(IssmDouble constant, int name){/*{{{*/

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

	/*update input*/
	this->inputs->AddInput(new DoubleInput(name,constant));
}
/*}}}*/
void Element::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));
}
/*}}}*/
bool Element::IsFloating(){/*{{{*/

	bool shelf;
	int  migration_style;
	parameters->FindParam(&migration_style,GroundinglineMigrationEnum);

	if(migration_style==SubelementMigrationEnum || migration_style==SubelementMigration2Enum){ //Floating if all nodes are floating
		if(this->inputs->Max(MaskGroundediceLevelsetEnum) <= 0.) shelf=true;
		else shelf=false;
	}
	else if(migration_style==NoneEnum || migration_style==AgressiveMigrationEnum || migration_style==SoftMigrationEnum){ //Floating if all nodes are floating
		if(this->inputs->Min(MaskGroundediceLevelsetEnum) > 0.) shelf=false;
		else shelf=true;
	}
	else _error_("migration_style not implemented yet");

	return shelf;
}/*}}}*/
bool Element::IsIceInElement(){/*{{{*/
	return (this->inputs->Min(MaskIceLevelsetEnum)<0.);
}
/*}}}*/
bool Element::IsInput(int name){/*{{{*/
	if (
				name==ThicknessEnum ||
				name==SurfaceEnum ||
				name==BaseEnum ||
				name==BedEnum ||
				name==SurfaceSlopeXEnum ||
				name==SurfaceSlopeYEnum ||
				name==SurfaceforcingsMassBalanceEnum ||
				name==BasalforcingsMeltingRateEnum ||
				name==BasalforcingsGeothermalfluxEnum ||
				name==SurfaceAreaEnum||
				name==DamageDEnum ||
				name==DamageDbarEnum ||
				name==PressureEnum ||
				name==VxEnum ||
				name==VyEnum ||
				name==VzEnum ||
				name==VxMeshEnum ||
				name==VyMeshEnum ||
				name==VzMeshEnum ||
				name==InversionVxObsEnum ||
				name==InversionVyObsEnum ||
				name==InversionVzObsEnum ||
				name==TemperatureEnum ||
				name==EnthalpyEnum ||
				name==EnthalpyPicardEnum ||
				name==WaterfractionEnum||
				name==WatercolumnEnum || 
				name==FrictionCoefficientEnum ||
				name==MaskGroundediceLevelsetEnum ||
				name==MaskIceLevelsetEnum ||
				name==IceMaskNodeActivationEnum ||
				name==LevelsetfunctionSlopeXEnum ||
				name==LevelsetfunctionSlopeYEnum ||
				name==LevelsetfunctionPicardEnum ||
				name==MasstransportCalvingrateEnum ||
				name==GradientEnum ||
				name==OldGradientEnum  ||
				name==ConvergedEnum || 
				name==QmuVxEnum ||
				name==QmuVyEnum ||
				name==QmuVzEnum ||
				name==QmuVxMeshEnum ||
				name==QmuVyMeshEnum ||
				name==QmuVzMeshEnum ||
				name==QmuPressureEnum ||
				name==QmuBaseEnum ||
				name==QmuThicknessEnum ||
				name==QmuSurfaceEnum ||
				name==QmuTemperatureEnum ||
				name==QmuMeltingEnum ||
				name==QmuMaskGroundediceLevelsetEnum ||
				name==QmuMaskIceLevelsetEnum ||
				name==QmuMaterialsRheologyBEnum ||
				name==MaterialsRheologyBEnum ||
				name==MaterialsRheologyBbarEnum ||
				name==MaterialsRheologyNEnum ||
				name==GiaWEnum || 
				name==GiadWdtEnum ||
				name==SedimentHeadEnum ||
				name==EplHeadEnum ||
				name==SedimentHeadOldEnum ||
				name==EplHeadOldEnum ||
				name==HydrologydcEplThicknessOldEnum ||
				name==HydrologydcEplInitialThicknessEnum ||
				name==HydrologydcEplThicknessEnum ||
				name==HydrologydcMaskEplactiveNodeEnum ||
				name==MeshVertexonbedEnum 

				) {
					return true;
				}
	else return false;
}
/*}}}*/
ElementVector* Element::NewElementVector(int approximation_enum){/*{{{*/
	return new ElementVector(nodes,this->GetNumberOfNodes(),this->parameters,approximation_enum);
}
/*}}}*/
ElementMatrix* Element::NewElementMatrix(int approximation_enum){/*{{{*/
	return new ElementMatrix(nodes,this->GetNumberOfNodes(),this->parameters,approximation_enum);
}
/*}}}*/
ElementMatrix* Element::NewElementMatrixCoupling(int number_nodes,int approximation_enum){/*{{{*/
	return new ElementMatrix(nodes,number_nodes,this->parameters,approximation_enum);
}
/*}}}*/
IssmDouble Element::PureIceEnthalpy(IssmDouble pressure){/*{{{*/
	return this->matpar->PureIceEnthalpy(pressure);
}/*}}}*/
void Element::ResultInterpolation(int* pinterpolation,int* pnodesperelement,int output_enum){/*{{{*/

	Input* input=this->inputs->GetInput(output_enum);

	/*If this input is not already in Inputs, maybe it needs to be computed?*/
	if(!input){
		switch(output_enum){
			case ViscousHeatingEnum:
				this->ViscousHeatingCreateInput();
				input=this->inputs->GetInput(output_enum);
				break;
			case StressTensorxxEnum: 
				this->ComputeStressTensor();
				input=this->inputs->GetInput(output_enum);
				break;
			case StressTensorxyEnum: 
				this->ComputeStressTensor();
				input=this->inputs->GetInput(output_enum);
				break;
			case StressTensorxzEnum: 
				this->ComputeStressTensor();
				input=this->inputs->GetInput(output_enum);
				break;
			case StressTensoryyEnum: 
				this->ComputeStressTensor();
				input=this->inputs->GetInput(output_enum);
				break;
			case StressTensoryzEnum: 
				this->ComputeStressTensor();
				input=this->inputs->GetInput(output_enum);
				break;
			case StressTensorzzEnum: 
				this->ComputeStressTensor();
				input=this->inputs->GetInput(output_enum);
				break;
			default:
				_error_("input "<<EnumToStringx(output_enum)<<" not found in element");
		}
	}

	/*Assign output pointer*/
	*pinterpolation   = input->GetResultInterpolation();
	*pnodesperelement = input->GetResultNumberOfNodes();
}/*}}}*/
void Element::ResultToVector(Vector<IssmDouble>* vector,int output_enum){/*{{{*/

	Input* input=this->inputs->GetInput(output_enum);
	if(!input) _error_("input "<<EnumToStringx(output_enum)<<" not found in element");

	switch(input->GetResultInterpolation()){
		case P0Enum:{
			IssmDouble  value;
			bool        bvalue;
			Input*      input = this->GetInput(output_enum); _assert_(input);
			switch(input->ObjectEnum()){
				case DoubleInputEnum:
					input->GetInputValue(&value);
					break;
				case BoolInputEnum:
					input->GetInputValue(&bvalue);
					value=reCast<IssmDouble>(bvalue);
					break;
				default:
					_error_("Input of type "<<EnumToStringx(input->GetResultInterpolation())<<" not supported");
			}
			vector->SetValue(this->Sid(),value,INS_VAL);
			break;
		}
		case P1Enum:{
			int         numvertices = this->GetNumberOfVertices();
			IssmDouble *values      = xNew<IssmDouble>(numvertices);
			int        *connectivity= xNew<int>(numvertices);
			int        *sidlist     = xNew<int>(numvertices);

			this->GetVerticesSidList(sidlist);
			this->GetVerticesConnectivityList(connectivity);
			this->GetInputListOnVertices(values,output_enum);
			for(int i=0;i<numvertices;i++) values[i] = values[i]/reCast<IssmDouble>(connectivity[i]);

			vector->SetValues(numvertices,sidlist,values,ADD_VAL);

			xDelete<IssmDouble>(values);
			xDelete<int>(connectivity);
			xDelete<int>(sidlist);
			break;
		}
		default:
					 _error_("interpolation "<<EnumToStringx(input->GetResultInterpolation())<<" not supported yet");
	}
} /*}}}*/
void Element::ResultToPatch(IssmDouble* values,int nodesperelement,int output_enum){/*{{{*/

	Input* input=this->inputs->GetInput(output_enum);
	if(!input) _error_("input "<<EnumToStringx(output_enum)<<" not found in element");

	input->ResultToPatch(values,nodesperelement,this->Sid());

} /*}}}*/
void Element::SetwiseNodeConnectivity(int* pd_nz,int* po_nz,Node* node,bool* flags,int* flagsindices,int set1_enum,int set2_enum){/*{{{*/

	/*Intermediaries*/
	const int numnodes = this->GetNumberOfNodes();

	/*Output */
	int d_nz = 0;
	int o_nz = 0;

	/*Loop over all nodes*/
	for(int i=0;i<numnodes;i++){

		if(!flags[this->nodes[i]->Lid()]){

			/*flag current node so that no other element processes it*/
			flags[this->nodes[i]->Lid()]=true;

			int counter=0;
			while(flagsindices[counter]>=0) counter++;
			flagsindices[counter]=this->nodes[i]->Lid();

			/*if node is clone, we have an off-diagonal non-zero, else it is a diagonal non-zero*/
			switch(set2_enum){
				case FsetEnum:
					if(nodes[i]->indexing.fsize){
						if(this->nodes[i]->IsClone())
						 o_nz += 1;
						else
						 d_nz += 1;
					}
					break;
				case GsetEnum:
					if(nodes[i]->indexing.gsize){
						if(this->nodes[i]->IsClone())
						 o_nz += 1;
						else
						 d_nz += 1;
					}
					break;
				case SsetEnum:
					if(nodes[i]->indexing.ssize){
						if(this->nodes[i]->IsClone())
						 o_nz += 1;
						else
						 d_nz += 1;
					}
					break;
				default: _error_("not supported");
			}
		}
	}

	/*Special case: 2d/3d coupling, the node of this element might be connected
	 *to the basal element*/
	int analysis_type,approximation,numlayers;
	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
	if(analysis_type==StressbalanceAnalysisEnum){
		inputs->GetInputValue(&approximation,ApproximationEnum);
		if(approximation==SSAHOApproximationEnum || approximation==SSAFSApproximationEnum){
			parameters->FindParam(&numlayers,MeshNumberoflayersEnum);
			o_nz += numlayers*3;
			d_nz += numlayers*3;
		}
	}

	/*Assign output pointers: */
	*pd_nz=d_nz;
	*po_nz=o_nz;
}
/*}}}*/
int  Element::Sid(){/*{{{*/

	return this->sid;

}
/*}}}*/
IssmDouble Element::TMeltingPoint(IssmDouble pressure){/*{{{*/
	_assert_(matpar);
	return this->matpar->TMeltingPoint(pressure);
}/*}}}*/
void Element::ViscousHeatingCreateInput(void){/*{{{*/

	/*Intermediaries*/
	IssmDouble phi;
	IssmDouble viscosity;
	IssmDouble epsilon[6];
	IssmDouble thickness;
	IssmDouble *xyz_list = NULL;

	/*Fetch number vertices and allocate memory*/
	int         numvertices    = this->GetNumberOfVertices();
	IssmDouble* viscousheating = xNew<IssmDouble>(numvertices);

	/*Retrieve all inputs and parameters*/
	this->GetVerticesCoordinatesBase(&xyz_list);
	Input* vx_input        = this->GetInput(VxEnum); _assert_(vx_input);
	Input* vy_input        = this->GetInput(VyEnum); _assert_(vy_input);
	Input* vz_input        = this->GetInput(VzEnum); _assert_(vz_input);
	Input* thickness_input = this->GetInput(ThicknessEnum); _assert_(thickness_input);

	/*loop over vertices: */
	Gauss* gauss=this->NewGauss();
	for (int iv=0;iv<numvertices;iv++){
		gauss->GaussVertex(iv);

		thickness_input->GetInputValue(&thickness,gauss);

		this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
		this->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input);
		this->GetPhi(&phi,&epsilon[0],viscosity);

		viscousheating[iv]=phi*thickness;
	}

	/*Create PentaVertex input, which will hold the basal friction:*/
	this->AddInput(ViscousHeatingEnum,viscousheating,P1Enum);

	/*Clean up and return*/
	xDelete<IssmDouble>(viscousheating);
	delete gauss;
}
/*}}}*/
void Element::ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/

	/*Intermediaries*/
	IssmDouble viscosity;
	IssmDouble epsilon3d[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
	IssmDouble epsilon2d[3]; /* epsilon=[exx,eyy,exy];            */
	IssmDouble eps_eff;
	IssmDouble eps0=1.e-27;

	if(dim==3){
		/* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */
		this->StrainRateFS(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input,vz_input);
		eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[3]*epsilon3d[3] +  epsilon3d[4]*epsilon3d[4] + epsilon3d[5]*epsilon3d[5] + epsilon3d[0]*epsilon3d[1]+eps0*eps0);
	}
	else{
		/* eps_eff^2 = exx^2 + eyy^2 + 2*exy^2 */
		this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
		eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + 2.*epsilon2d[2]*epsilon2d[2]);
	}

	/*Get viscosity*/
	material->GetViscosity(&viscosity,eps_eff);

	/*Assign output pointer*/
	*pviscosity=viscosity;
}
/*}}}*/
void Element::ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surface_input){/*{{{*/
	/*Compute the L1L2 viscosity
	 *
	 *      1
	 * mu = - A^-1 (sigma'_e)^(1-n)
	 *      2
	 *
	 * sigma'_e^2 = |sigma'_//|^2 + |sigma'_perp|^2 (see Perego 2012 eq. 17,18)
	 *
	 * L1L2 assumptions:
	 *
	 * (1) |eps_b|_// = A (|sigma'_//|^2 + |sigma'_perp|^2)^((n-1)/2) |sigma'_//|
	 * (2) |sigma'_perp|^2 = |rho g (s-z) grad(s)|^2
	 *
	 * Assuming that n = 3, we have a polynom of degree 3 to solve (the only unkown is X=|sigma'_//|)
	 *
	 * A X^3 + A |rho g (s-z) grad(s)|^2 X - |eps_b|_// = 0     */

	IssmDouble z,s,viscosity,p,q,delta;
	IssmDouble tau_perp,tau_par,eps_b,A;
	IssmDouble epsilon[5];   /*exx eyy exy exz eyz*/
	IssmDouble slope[3];

	/*Check that both inputs have been found*/
	if (!vx_input || !vy_input || !surface_input) _error_("Input missing");

	/*Get tau_perp*/
	surface_input->GetInputValue(&s,gauss);
	surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
	z=GetZcoord(gauss);
	tau_perp = matpar->GetRhoIce() * matpar->GetG() * fabs(s-z)*sqrt(slope[0]*slope[0]+slope[1]*slope[1]);

	/* Get eps_b*/
	this->StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
	eps_b = sqrt(epsilon[0]*epsilon[0] + epsilon[1]*epsilon[1] + epsilon[0]*epsilon[1] + epsilon[2]*epsilon[2]);
	if(eps_b==0.){
		*pviscosity = 2.5e+17;
		return;
	}

	/*Get A*/
	_assert_(material->GetN()==3.0);
	A=material->GetA();

	/*Solve for tau_perp (http://fr.wikipedia.org/wiki/Méthode_de_Cardan)*/
	p     = tau_perp *tau_perp;
	q     = - eps_b/A;
	delta = q *q + p*p*p*4./27.;
	_assert_(delta>0);
	tau_par = pow(0.5*(-q+sqrt(delta)),1./3.) - pow(0.5*(q+sqrt(delta)),1./3.);

	/*Viscosity*/
	viscosity = 1./(2.*A*(tau_par*tau_par + tau_perp*tau_perp));
	_assert_(!xIsNan(viscosity));
	_assert_(viscosity > 0.);

	/*Assign output pointer*/
	*pviscosity = viscosity;
}/*}}}*/
void Element::ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/

	/*Intermediaries*/
	IssmDouble viscosity;
	IssmDouble epsilon3d[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/
	IssmDouble epsilon2d[2]; /* epsilon=[exx,exy];           */
	IssmDouble eps_eff;

	if(dim==3){
		/* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */
		this->StrainRateHO(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input);
		eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[2]*epsilon3d[2] +  epsilon3d[3]*epsilon3d[3] + epsilon3d[4]*epsilon3d[4] + epsilon3d[0]*epsilon3d[1]);
	}
	else{
		/* eps_eff^2 = exx^2 + exy^2 */
		this->StrainRateHO2dvertical(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
		eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1]);
	}

	/*Get viscosity*/
	material->GetViscosity(&viscosity,eps_eff);

	/*Assign output pointer*/
	*pviscosity=viscosity;
}/*}}}*/
void Element::ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/

	/*Intermediaries*/
	IssmDouble viscosity;
	IssmDouble epsilon2d[3];/* epsilon=[exx,eyy,exy];    */
	IssmDouble epsilon1d;   /* epsilon=[exx];    */
	IssmDouble eps_eff;

	 if(dim==2){
		 /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy+*/
		 this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
		 eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + epsilon2d[2]*epsilon2d[2] + epsilon2d[0]*epsilon2d[1]);
	 }
	 else{
		 /* eps_eff^2 = exx^2*/
		 this->StrainRateSSA1d(&epsilon1d,xyz_list,gauss,vx_input);
		 eps_eff = sqrt(epsilon1d*epsilon1d);
	 }

	/*Get viscosity*/
	material->GetViscosityBar(&viscosity,eps_eff);

	/*Assign output pointer*/
	*pviscosity=viscosity;
}/*}}}*/
void Element::ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
	this->material->GetViscosity2dDerivativeEpsSquare(pmu_prime,epsilon);
}/*}}}*/
void Element::ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
	this->material->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon);
}/*}}}*/
void Element::ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
	this->material->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon);
}/*}}}*/
void Element::StrainRateFS(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/
	/*Compute the 3d Strain Rate (6 components):
	 *
	 * epsilon=[exx eyy ezz exy exz eyz]
	 */

	/*Intermediaries*/
	IssmDouble dvx[3];
	IssmDouble dvy[3];
	IssmDouble dvz[3];

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

	/*Get strain rate assuming that epsilon has been allocated*/
	vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
	vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
	vz_input->GetInputDerivativeValue(&dvz[0],xyz_list,gauss);
	epsilon[0] = dvx[0];
	epsilon[1] = dvy[1];
	epsilon[2] = dvz[2];
	epsilon[3] = 0.5*(dvx[1] + dvy[0]);
	epsilon[4] = 0.5*(dvx[2] + dvz[0]);
	epsilon[5] = 0.5*(dvy[2] + dvz[1]);

}/*}}}*/
void Element::StrainRateHO(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
	/*Compute the 3d Blatter/HOStrain Rate (5 components):
	 *
	 * epsilon=[exx eyy exy exz eyz]
	 *
	 * with exz=1/2 du/dz
	 *      eyz=1/2 dv/dz
	 *
	 * the contribution of vz is neglected
	 */

	/*Intermediaries*/
	IssmDouble dvx[3];
	IssmDouble dvy[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: " << vx_input << ", vy: " << vy_input << "\n");
	}

	/*Get strain rate assuming that epsilon has been allocated*/
	vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
	vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
	epsilon[0] = dvx[0];
	epsilon[1] = dvy[1];
	epsilon[2] = 0.5*(dvx[1] + dvy[0]);
	epsilon[3] = 0.5*dvx[2];
	epsilon[4] = 0.5*dvy[2];

}/*}}}*/
void Element::StrainRateHO2dvertical(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
	/*Compute the 2d Blatter/HOStrain Rate (2 components):
	 *
	 * epsilon=[exx exz]
	 *
	 * with exz=1/2 du/dz
	 *
	 * the contribution of vz is neglected
	 */

	/*Intermediaries*/
	IssmDouble dvx[3];

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

	/*Get strain rate assuming that epsilon has been allocated*/
	vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
	epsilon[0] = dvx[0];
	epsilon[1] = 0.5*dvx[1];

}/*}}}*/
void Element::StrainRateSSA(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/

	/*Intermediaries*/
	IssmDouble dvx[3];
	IssmDouble dvy[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: " << vx_input << ", vy: " << vy_input << "\n");
	}

	/*Get strain rate assuming that epsilon has been allocated*/
	vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
	vy_input->GetInputDerivativeValue(&dvy[0],xyz_list,gauss);
	epsilon[0] = dvx[0];
	epsilon[1] = dvy[1];
	epsilon[2] = 0.5*(dvx[1] + dvy[0]);

}/*}}}*/
void Element::StrainRateSSA1d(IssmDouble* epsilon,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input){/*{{{*/

	/*Intermediaries*/
	IssmDouble dvx[3];

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

	/*Get strain rate assuming that epsilon has been allocated*/
	vx_input->GetInputDerivativeValue(&dvx[0],xyz_list,gauss);
	*epsilon = dvx[0];

}/*}}}*/
void Element::TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,int transformenum){/*{{{*/

	/*All nodes have the same Coordinate System*/
	int numnodes  = this->GetNumberOfNodes();
	int* cs_array = xNew<int>(numnodes);
	for(int i=0;i<numnodes;i++) cs_array[i]=transformenum;

	/*Call core*/
	TransformInvStiffnessMatrixCoord(Ke,this->nodes,numnodes,cs_array);

	/*Clean-up*/
	xDelete<int>(cs_array);
}/*}}}*/
void Element::TransformInvStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes_list,int numnodes,int* cs_array){/*{{{*/

	int         i,j;
	int         numdofs   = 0;
	IssmDouble *transform = NULL;
	IssmDouble *values    = NULL;

	/*Get total number of dofs*/
	for(i=0;i<numnodes;i++){
		switch(cs_array[i]){
			case PressureEnum: numdofs+=1; break;
			case XYEnum:       numdofs+=2; break;
			case XYZEnum:      numdofs+=3; break;
			default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
		}
	}

	/*Copy current stiffness matrix*/
	values=xNew<IssmDouble>(Ke->nrows*Ke->ncols);
	for(i=0;i<Ke->nrows;i++) for(j=0;j<Ke->ncols;j++) values[i*Ke->ncols+j]=Ke->values[i*Ke->ncols+j];

	/*Get Coordinate Systems transform matrix*/
	CoordinateSystemTransform(&transform,nodes_list,numnodes,cs_array);

	/*Transform matrix: R*Ke*R^T */
	TripleMultiply(transform,numdofs,numdofs,0,
				values,Ke->nrows,Ke->ncols,0,
				transform,numdofs,numdofs,1,
				&Ke->values[0],0);

	/*Free Matrix*/
	xDelete<IssmDouble>(transform);
	xDelete<IssmDouble>(values);
}/*}}}*/
void Element::TransformLoadVectorCoord(ElementVector* pe,int transformenum){/*{{{*/

	/*All nodes have the same Coordinate System*/
	int  numnodes = this->GetNumberOfNodes();
	int* cs_array = xNew<int>(numnodes);
	for(int i=0;i<numnodes;i++) cs_array[i]=transformenum;

	/*Call core*/
	this->TransformLoadVectorCoord(pe,this->nodes,numnodes,cs_array);

	/*Clean-up*/
	xDelete<int>(cs_array);
}/*}}}*/
void Element::TransformLoadVectorCoord(ElementVector* pe,int* cs_array){/*{{{*/

	this->TransformLoadVectorCoord(pe,this->nodes,this->GetNumberOfNodes(),cs_array);

}/*}}}*/
void Element::TransformLoadVectorCoord(ElementVector* pe,Node** nodes_list,int numnodes,int* cs_array){/*{{{*/

	int         i;
	int         numdofs   = 0;
	IssmDouble *transform = NULL;
	IssmDouble *values    = NULL;

	/*Get total number of dofs*/
	for(i=0;i<numnodes;i++){
		switch(cs_array[i]){
			case PressureEnum: numdofs+=1; break;
			case XYEnum:       numdofs+=2; break;
			case XYZEnum:      numdofs+=3; break;
			default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
		}
	}

	/*Copy current load vector*/
	values=xNew<IssmDouble>(pe->nrows);
	for(i=0;i<pe->nrows;i++) values[i]=pe->values[i];

	/*Get Coordinate Systems transform matrix*/
	CoordinateSystemTransform(&transform,nodes_list,numnodes,cs_array);

	/*Transform matrix: R^T*pe */
	MatrixMultiply(transform,numdofs,numdofs,1,
				values,pe->nrows,1,0,
				&pe->values[0],0);

	/*Free Matrices*/
	xDelete<IssmDouble>(transform);
	xDelete<IssmDouble>(values);
}/*}}}*/
void Element::TransformSolutionCoord(IssmDouble* values,int transformenum){/*{{{*/

	/*All nodes have the same Coordinate System*/
	int  numnodes = this->GetNumberOfNodes();
	int* cs_array = xNew<int>(numnodes);
	for(int i=0;i<numnodes;i++) cs_array[i]=transformenum;

	/*Call core*/
	this->TransformSolutionCoord(values,this->nodes,numnodes,cs_array);

	/*Clean-up*/
	xDelete<int>(cs_array);
}/*}}}*/
void Element::TransformSolutionCoord(IssmDouble* values,int* transformenum_list){/*{{{*/
	this->TransformSolutionCoord(values,this->nodes,this->GetNumberOfNodes(),transformenum_list);
}/*}}}*/
void Element::TransformSolutionCoord(IssmDouble* values,int numnodes,int transformenum){/*{{{*/

	/*All nodes have the same Coordinate System*/
	int* cs_array = xNew<int>(numnodes);
	for(int i=0;i<numnodes;i++) cs_array[i]=transformenum;

	/*Call core*/
	this->TransformSolutionCoord(values,this->nodes,numnodes,cs_array);

	/*Clean-up*/
	xDelete<int>(cs_array);
}/*}}}*/
void Element::TransformSolutionCoord(IssmDouble* solution,int numnodes,int* cs_array){/*{{{*/
	this->TransformSolutionCoord(solution,this->nodes,numnodes,cs_array);
}/*}}}*/
void Element::TransformSolutionCoord(IssmDouble* values,Node** nodes_list,int numnodes,int transformenum){/*{{{*/
	/*NOT NEEDED*/
	/*All nodes have the same Coordinate System*/
	int* cs_array = xNew<int>(numnodes);
	for(int i=0;i<numnodes;i++) cs_array[i]=transformenum;

	/*Call core*/
	this->TransformSolutionCoord(values,nodes_list,numnodes,cs_array);

	/*Clean-up*/
	xDelete<int>(cs_array);
}/*}}}*/
void Element::TransformSolutionCoord(IssmDouble* solution,Node** nodes_list,int numnodes,int* cs_array){/*{{{*/

	int         i;
	int         numdofs   = 0;
	IssmDouble *transform = NULL;
	IssmDouble *values    = NULL;

	/*Get total number of dofs*/
	for(i=0;i<numnodes;i++){
		switch(cs_array[i]){
			case PressureEnum: numdofs+=1; break;
			case XYEnum:       numdofs+=2; break;
			case XYZEnum:      numdofs+=3; break;
			default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
		}
	}

	/*Copy current solution vector*/
	values=xNew<IssmDouble>(numdofs);
	for(i=0;i<numdofs;i++) values[i]=solution[i];

	/*Get Coordinate Systems transform matrix*/
	CoordinateSystemTransform(&transform,nodes_list,numnodes,cs_array);

	/*Transform matrix: R*U */
	MatrixMultiply(transform,numdofs,numdofs,0,
				values,numdofs,1,0,
				&solution[0],0);

	/*Free Matrices*/
	xDelete<IssmDouble>(transform);
	xDelete<IssmDouble>(values);
}/*}}}*/
void Element::TransformStiffnessMatrixCoord(ElementMatrix* Ke,int transformenum){/*{{{*/

	/*All nodes have the same Coordinate System*/
	int  numnodes = this->GetNumberOfNodes();
	int* cs_array = xNew<int>(numnodes);
	for(int i=0;i<numnodes;i++) cs_array[i]=transformenum;

	/*Call core*/
	this->TransformStiffnessMatrixCoord(Ke,this->nodes,numnodes,cs_array);

	/*Clean-up*/
	xDelete<int>(cs_array);
}/*}}}*/
void Element::TransformStiffnessMatrixCoord(ElementMatrix* Ke,int* transformenum_list){/*{{{*/
	this->TransformStiffnessMatrixCoord(Ke,this->nodes,this->GetNumberOfNodes(),transformenum_list);
}/*}}}*/
void Element::TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes_list,int numnodes,int* cs_array){/*{{{*/

	int         numdofs = 0;
	IssmDouble *transform = NULL;
	IssmDouble *values    = NULL;

	/*Get total number of dofs*/
	for(int i=0;i<numnodes;i++){
		switch(cs_array[i]){
			case PressureEnum: numdofs+=1; break;
			case XYEnum:       numdofs+=2; break;
			case XYZEnum:      numdofs+=3; break;
			default: _error_("Coordinate system " << EnumToStringx(cs_array[i]) << " not supported yet");
		}
	}

	/*Copy current stiffness matrix*/
	values=xNew<IssmDouble>(Ke->nrows*Ke->ncols);
	for(int i=0;i<Ke->nrows*Ke->ncols;i++) values[i]=Ke->values[i];

	/*Get Coordinate Systems transform matrix*/
	CoordinateSystemTransform(&transform,nodes_list,numnodes,cs_array);

	/*Transform matrix: R^T*Ke*R */
	TripleMultiply(transform,numdofs,numdofs,1,
				values,Ke->nrows,Ke->ncols,0,
				transform,numdofs,numdofs,0,
				&Ke->values[0],0);

	/*Free Matrix*/
	xDelete<IssmDouble>(transform);
	xDelete<IssmDouble>(values);
}/*}}}*/
