/*!\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;
	this->element_type_list=NULL;
}/*}}}*/
Element::~Element(){/*{{{*/
	xDelete<int>(element_type_list);
	delete inputs;
}
/*}}}*/

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

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

	IssmDouble *xyz_list=NULL;
	IssmDouble  eps_xx,eps_xy,eps_yy,eps_xz,eps_yz,eps_zz,eps_eff;
	IssmDouble  epsmin=1.e-27;
	IssmDouble  eps_0,eps_f,sigma_0,B,D,n;
	int         dim,counter=0;
	IssmDouble  k1,k2,threshold=1.e-12;

	/* Retrieve parameters */
	this->GetVerticesCoordinates(&xyz_list);
	this->ComputeStrainRate();
	parameters->FindParam(&dim,DomainDimensionEnum);
	parameters->FindParam(&eps_f,DamageC1Enum);
	parameters->FindParam(&sigma_0,DamageStressThresholdEnum);

	/* Retrieve inputs */
	Input* eps_xx_input=this->GetInput(StrainRatexxEnum); _assert_(eps_xx_input);
	Input* eps_yy_input=this->GetInput(StrainRateyyEnum); _assert_(eps_yy_input);
	Input* eps_xy_input=this->GetInput(StrainRatexyEnum); _assert_(eps_xy_input);
	Input* eps_xz_input=NULL;
	Input* eps_yz_input=NULL;
	Input* eps_zz_input=NULL;
	if(dim==3){
		eps_xz_input=this->GetInput(StrainRatexzEnum); _assert_(eps_xz_input);
		eps_yz_input=this->GetInput(StrainRateyzEnum); _assert_(eps_yz_input);
		eps_zz_input=this->GetInput(StrainRatezzEnum); _assert_(eps_zz_input);
	}

	/* Fetch number of nodes and allocate output*/
   int numnodes = this->GetNumberOfNodes();
   IssmDouble* newD = xNew<IssmDouble>(numnodes);

	/* Retrieve domain-dependent inputs */
	Input* n_input=this->GetInput(MaterialsRheologyNEnum); _assert_(n_input);
   Input* damage_input = NULL;
   Input* B_input = NULL;
	int domaintype;
   parameters->FindParam(&domaintype,DomainTypeEnum);
	if(domaintype==Domain2DhorizontalEnum){
      damage_input = this->GetInput(DamageDbarEnum);  _assert_(damage_input);
      B_input=this->GetInput(MaterialsRheologyBbarEnum); _assert_(B_input);
   }
   else{
      damage_input = this->GetInput(DamageDEnum);   _assert_(damage_input);
      B_input=this->GetInput(MaterialsRheologyBEnum); _assert_(B_input);
   }
	
	/* Start looping on the number of nodes: */
	Gauss* gauss=this->NewGauss();
	for (int i=0;i<numnodes;i++){
		gauss->GaussNode(this->GetElementType(),i);

		eps_xx_input->GetInputValue(&eps_xx,gauss);
		eps_yy_input->GetInputValue(&eps_yy,gauss);
		eps_xy_input->GetInputValue(&eps_xy,gauss);
		if(dim==3){
			eps_xz_input->GetInputValue(&eps_xz,gauss);
			eps_yz_input->GetInputValue(&eps_yz,gauss);
			eps_zz_input->GetInputValue(&eps_zz,gauss);
		}
		else{eps_xz=0; eps_yz=0; eps_zz=0;}

		/* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */
		eps_eff=sqrt(eps_xx*eps_xx+eps_yy*eps_yy+eps_xy*eps_xy+eps_xz*eps_xz+eps_yz*eps_yz+eps_xx*eps_yy+epsmin*epsmin);

		B_input->GetInputValue(&B,gauss);
      n_input->GetInputValue(&n,gauss);
      damage_input->GetInputValue(&D,gauss);

		/* Compute threshold strain rate from threshold stress */
		eps_0=pow(sigma_0/B,n); 
		_assert_(eps_f>eps_0);

		/* Compute kappa (k) from pre-existing level of damage, using Newton-Raphson */
		/* provide a reasonable initial guess */
		if(D==0){
			k1=eps_0;
		}
		else{
			k1=exp(n*eps_0/(eps_f-eps_0)-n*log(1.-D)+log(eps_0)); /* initial guess */
		}
	
		counter=0;
		while(true){
      	/*Newton step k2=k1-f(k1)/f'(k1) */
			k2=k1-(k1+(eps_f-eps_0)/n*log(k1)-eps_0+(eps_f-eps_0)*(log(1.-D)-1./n*log(eps_0)))/(1.+(eps_f-eps_0)/n/k1);

      	if( fabs(k2-k1)/(fabs(k2))<threshold ){
         	break;
      	}
      	else{
				k1=k2;
      	   counter++;
      	}

      	if(counter>50) break;
   	}

		if(eps_eff>k2){
			newD[i]=1.-pow(eps_0/eps_eff,1./n)*exp(-(eps_eff-eps_0)/(eps_f-eps_0));
		}
		else newD[i]=D;
	}

	/* Add new damage input to DamageEnum and NewDamageEnum */
	this->AddInput(NewDamageEnum,newD,this->GetElementType());
	if(domaintype==Domain2DhorizontalEnum){
		this->AddInput(DamageDbarEnum,newD,this->GetElementType());
	}
	else{
		this->AddInput(DamageDEnum,newD,this->GetElementType());
	}

	/*Clean up and return*/
	xDelete<IssmDouble>(xyz_list);
	xDelete<IssmDouble>(newD);
	delete gauss;

}/*}}}*/
void       Element::ComputeStrainRate(){/*{{{*/

	int         dim;
	IssmDouble *xyz_list = NULL;
	IssmDouble  epsilon[6];

	/*Retrieve all inputs we will be needing: */
	this->GetVerticesCoordinates(&xyz_list);
	parameters->FindParam(&dim,DomainDimensionEnum);
	Input* vx_input=this->GetInput(VxEnum); _assert_(vx_input);
	Input* vy_input=this->GetInput(VyEnum); _assert_(vy_input);
	Input* vz_input=NULL;
	if(dim==3){vz_input=this->GetInput(VzEnum); _assert_(vz_input);}

	/*Allocate arrays*/
	int numvertices = this->GetNumberOfVertices();
	IssmDouble* eps_xx = xNew<IssmDouble>(numvertices);
	IssmDouble* eps_yy = xNew<IssmDouble>(numvertices);
	IssmDouble* eps_zz = xNew<IssmDouble>(numvertices);
	IssmDouble* eps_xy = xNew<IssmDouble>(numvertices);
	IssmDouble* eps_xz = xNew<IssmDouble>(numvertices);
	IssmDouble* eps_yz = xNew<IssmDouble>(numvertices);

	/* Start looping on the number of vertices: */
	Gauss* gauss=this->NewGauss();
	for (int iv=0;iv<numvertices;iv++){
		gauss->GaussVertex(iv);

		/*Compute strain rate viscosity and pressure: */
		if(dim==2)
		 this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
		else
		 this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);

		if(dim==2){
			 /* epsilon=[exx,eyy,exy];*/
			eps_xx[iv]=epsilon[0]; 
			eps_yy[iv]=epsilon[1];
			eps_xy[iv]=epsilon[2];
		}
		else{
			/*epsilon=[exx eyy ezz exy exz eyz]*/
			eps_xx[iv]=epsilon[0]; 
			eps_yy[iv]=epsilon[1];
			eps_zz[iv]=epsilon[2];
			eps_xy[iv]=epsilon[3]; 
			eps_xz[iv]=epsilon[4];
			eps_yz[iv]=epsilon[5];
		}
	}

	/*Add Stress tensor components into inputs*/
	this->AddInput(StrainRatexxEnum,eps_xx,P1Enum);
	this->AddInput(StrainRatexyEnum,eps_xy,P1Enum);
	this->AddInput(StrainRatexzEnum,eps_xz,P1Enum);
	this->AddInput(StrainRateyyEnum,eps_yy,P1Enum);
	this->AddInput(StrainRateyzEnum,eps_yz,P1Enum);
	this->AddInput(StrainRatezzEnum,eps_zz,P1Enum);

	/*Clean up and return*/
	delete gauss;
	xDelete<IssmDouble>(xyz_list);
	xDelete<IssmDouble>(eps_xx);
	xDelete<IssmDouble>(eps_yy);
	xDelete<IssmDouble>(eps_zz);
	xDelete<IssmDouble>(eps_xy);
	xDelete<IssmDouble>(eps_xz);
	xDelete<IssmDouble>(eps_yz);

}
/*}}}*/
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++) {
			_printf_("nodes[" << i << "] = " << nodes[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*/
	int        dim;
	IssmDouble Jdet;
	IssmDouble divergence=0.;
	IssmDouble dvx[3],dvy[3],dvz[3];
	IssmDouble *xyz_list = NULL;

	/*Get inputs and parameters*/
	this->FindParam(&dim,DomainDimensionEnum);
	Input* vx_input = this->GetInput(VxEnum); _assert_(vx_input);
	Input* vy_input = this->GetInput(VyEnum); _assert_(vy_input);
	Input* vz_input = NULL;
	if(dim==3){
		vz_input = this->GetInput(VzEnum); _assert_(vz_input);
	}
	this->GetVerticesCoordinates(&xyz_list);

	Gauss* gauss=this->NewGauss(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);
		if(dim==2){
			divergence += (dvx[0]+dvy[1])*gauss->weight*Jdet;
		}
		else{
			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::dViscositydBSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/

	/*Intermediaries*/
	IssmDouble dmudB;
	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 = 1/2 exx^2*/
		 this->StrainRateSSA1d(&epsilon1d,xyz_list,gauss,vx_input);
		 eps_eff = sqrt(epsilon1d*epsilon1d/2.);
	 }

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

	/*Assign output pointer*/
	*pdmudB=dmudB;

}
/*}}}*/
void       Element::dViscositydDSSA(IssmDouble* pdmudB,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/

	/*Intermediaries*/
	IssmDouble dmudB;
	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 = 1/2 exx^2*/
		this->StrainRateSSA1d(&epsilon1d,xyz_list,gauss,vx_input);
		eps_eff = sqrt(epsilon1d*epsilon1d/2.);
	}

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

	/*Assign output pointer*/
	*pdmudB=dmudB;

}
/*}}}*/
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::GetVectorFromInputs(Vector<IssmDouble>* vector,int input_enum){/*{{{*/

	/*Fetch number vertices for this element and allocate arrays*/
	int numvertices = this->GetNumberOfVertices();
	int*        vertexpidlist = xNew<int>(numvertices);
	IssmDouble* values        = xNew<IssmDouble>(numvertices);

	/*Fill in values*/
	this->GetVertexPidList(vertexpidlist);
	this->GetInputListOnVertices(values,input_enum);
	vector->SetValues(numvertices,vertexpidlist,values,INS_VAL);

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

}
/*}}}*/
void       Element::GetVertexPidList(int* pidlist){/*{{{*/

	int numvertices = this->GetNumberOfVertices();
	for(int i=0;i<numvertices;i++) pidlist[i]=vertices[i]->Pid();

}
/*}}}*/
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();
}
/*}}}*/
IssmDouble Element::GetXcoord(IssmDouble* xyz_list,Gauss* gauss){/*{{{*/

	/*output*/
	IssmDouble x;

	/*Create list of x*/
	int         numvertices = this->GetNumberOfVertices();
	IssmDouble* x_list      = xNew<IssmDouble>(numvertices);

	for(int i=0;i<numvertices;i++) x_list[i]=xyz_list[i*3+0];
	ValueP1OnGauss(&x,x_list,gauss);

	xDelete<IssmDouble>(x_list);
	return x;
}/*}}}*/
IssmDouble Element::GetYcoord(IssmDouble* xyz_list,Gauss* gauss){/*{{{*/

	/*output*/
	IssmDouble y;

	/*Create list of y*/
	int         numvertices = this->GetNumberOfVertices();
	IssmDouble* y_list      = xNew<IssmDouble>(numvertices);

	for(int i=0;i<numvertices;i++) y_list[i]=xyz_list[i*3+1];
	ValueP1OnGauss(&y,y_list,gauss);

	xDelete<IssmDouble>(y_list);
	return y;
}/*}}}*/
IssmDouble Element::GetZcoord(IssmDouble* xyz_list,Gauss* gauss){/*{{{*/

	/*output*/
	IssmDouble z;

	/*Create list of z*/
	int         numvertices = this->GetNumberOfVertices();
	IssmDouble* z_list      = xNew<IssmDouble>(numvertices);

	for(int i=0;i<numvertices;i++) z_list[i]=xyz_list[i*3+2];
	ValueP1OnGauss(&z,z_list,gauss);

	xDelete<IssmDouble>(z_list);
	return z;
}/*}}}*/
void       Element::GradientIndexing(int* indexing,int control_index){/*{{{*/

	/*Get number of controls*/
	int num_controls;
	parameters->FindParam(&num_controls,InversionNumControlParametersEnum);

	/*Get number of vertices*/
	int numvertices = this->GetNumberOfVertices();

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

}
/*}}}*/
bool       Element::HasNodeOnBase(){/*{{{*/
	return (this->inputs->Max(MeshVertexonbaseEnum)>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;
	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::InputDuplicate(int original_enum,int new_enum){/*{{{*/

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

}
/*}}}*/
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 || migration_style==ContactEnum){ //Floating if all nodes are floating
		if(this->inputs->Max(MaskGroundediceLevelsetEnum) <= 0.) shelf=true;
		else shelf=false;
	}
	else if(migration_style==NoneEnum || migration_style==AggressiveMigrationEnum || 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==BalancethicknessThickeningRateEnum ||
				name==BalancethicknessApparentMassbalanceEnum ||
				name==SigmaNNEnum ||
				name==SurfaceSlopeXEnum ||
				name==SurfaceSlopeYEnum ||
				name==SurfaceforcingsMassBalanceEnum ||
				name==BasalforcingsGroundediceMeltingRateEnum ||
				name==BasalforcingsFloatingiceMeltingRateEnum ||
				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==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==MeshVertexonbaseEnum 

				) {
					return true;
				}
	else return false;
}
/*}}}*/
void       Element::LinearFloatingiceMeltingRate(){/*{{{*/

	int numvertices      = this->GetNumberOfVertices();
	IssmDouble  deepwaterel,upperwaterel,deepwatermelt;
	IssmDouble* base     = xNew<IssmDouble>(numvertices);
	IssmDouble* values   = xNew<IssmDouble>(numvertices);

	parameters->FindParam(&deepwaterel,BasalforcingsDeepwaterElevationEnum);
	parameters->FindParam(&upperwaterel,BasalforcingsUpperwaterElevationEnum);
	parameters->FindParam(&deepwatermelt,BasalforcingsDeepwaterMeltingRateEnum);

	this->GetInputListOnVertices(base,BaseEnum);
	for(int i=0;i<numvertices;i++){
		if(base[i]>upperwaterel)      values[i]=0;
		else if (base[i]<deepwaterel) values[i]=deepwatermelt;
		else values[i]=deepwatermelt*(base[i]-upperwaterel)/(deepwaterel-upperwaterel);
	}

	this->AddInput(BasalforcingsFloatingiceMeltingRateEnum,values,P1Enum);

}/*}}}*/
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: 
			case StressTensorxyEnum: 
			case StressTensorxzEnum: 
			case StressTensoryyEnum: 
			case StressTensoryzEnum: 
			case StressTensorzzEnum: 
				this->ComputeStressTensor();
				input=this->inputs->GetInput(output_enum);
				break;
			case StrainRatexxEnum:
			case StrainRatexyEnum:
			case StrainRatexzEnum:
			case StrainRateyyEnum:
			case StrainRateyzEnum:
			case StrainRatezzEnum:
				this->ComputeStrainRate();
				input=this->inputs->GetInput(output_enum);
				break;
			case DeviatoricStressxxEnum: 
			case DeviatoricStressxyEnum: 
			case DeviatoricStressxzEnum: 
			case DeviatoricStressyyEnum: 
			case DeviatoricStressyzEnum: 
			case DeviatoricStresszzEnum: 
				this->ComputeDeviatoricStressTensor();
				input=this->inputs->GetInput(output_enum);
				break;
			case SigmaNNEnum: 
				this->ComputeSigmaNN();
				input=this->inputs->GetInput(output_enum);
				break;
			case NewDamageEnum:
				this->ComputeNewDamage();
				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:
					Gauss* gauss = this->NewGauss();
					input->GetInputValue(&value,gauss);
					delete gauss;
			}
			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){/*{{{*/
	/*The effective strain rate is defined in Paterson 3d Ed p 91 eq 9,
	 * and Cuffey p 303 eq 8.18:
	 *
	 *  2 eps_eff^2 = eps_xx^2 + eps_yy^2 + eps_zz^2 + 2(eps_xy^2 + eps_xz^2 + eps_yz^2)
	 *
	 *  or
	 *
	 *  eps_eff = 1/sqrt(2) sqrt( \sum_ij eps_ij^2 )
	 *
	 *          = 1/sqrt(2) ||eps||_F
	 *
	 *  where ||.||_F is the Frobenius norm */

	/*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 = 1/2 ( exx^2 + eyy^2 + 2*exy^2 )*/
		this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
		eps_eff = 1./sqrt(2.)*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=this->GetZcoord(xyz_list,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 = 1/2 (exx^2 + 2*exy^2 )*/
		this->StrainRateHO2dvertical(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
		eps_eff = 1./sqrt(2.)*sqrt(epsilon2d[0]*epsilon2d[0] + 2.*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 = 1/2 exx^2*/
		 this->StrainRateSSA1d(&epsilon1d,xyz_list,gauss,vx_input);
		 eps_eff = sqrt(epsilon1d*epsilon1d/2.);
	 }

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