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

/*Headers:*/
/*{{{*/
#ifdef HAVE_CONFIG_H
	#include <config.h>
#else
#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
#endif

#include "../classes.h"
#include "../../shared/shared.h"
/*}}}*/

/*Element macros*/
#define NUMNODESP1  3
#define NUMNODESP1b 4
#define NUMNODESP2  6

/*Object constructors and destructor*/
/*FUNCTION TriaRef::TriaRef(){{{*/
TriaRef::TriaRef(){
	this->element_type_list=NULL;
}
/*}}}*/
/*FUNCTION TriaRef::TriaRef(int* types,int nummodels){{{*/
TriaRef::TriaRef(const int nummodels){

	/*Only allocate pointer*/
	element_type_list=xNew<int>(nummodels);

}
/*}}}*/
/*FUNCTION TriaRef::~TriaRef(){{{*/
TriaRef::~TriaRef(){
	xDelete<int>(element_type_list);
}
/*}}}*/

/*Management*/
/*FUNCTION TriaRef::SetElementType{{{*/
void TriaRef::SetElementType(int type,int type_counter){

	_assert_(type==P1Enum || type==P1DGEnum || type==P1bubbleEnum || type==P1bubblecondensedEnum || type==P2Enum);

	/*initialize element type*/
	this->element_type_list[type_counter]=type;
}
/*}}}*/

/*Reference Element numerics*/
/*FUNCTION TriaRef::GetBHydro {{{*/
void TriaRef::GetBHydro(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss){
	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
	 * For node i, Bi can be expressed in the actual coordinate system
	 * by: 
	 *       Bi=[ dN/dx ]
	 *          [ dN/dy ]
	 * where N is the interpolation function for node i.
	 *
	 * We assume B has been allocated already, of size: 3x(NDOF2*numnodes)
	 */

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

	/*Get nodal functions derivatives*/
	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);

	/*Build B: */
	for(int i=0;i<numnodes;i++){
		B[numnodes*0+i] = dbasis[0*numnodes+i];
		B[numnodes*1+i] = dbasis[1*numnodes+i];
	}

	/*Clean-up*/
	xDelete<IssmDouble>(dbasis);
}
/*}}}*/
/*FUNCTION TriaRef::GetBSSA {{{*/
void TriaRef::GetBSSA(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss){
	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
	 * For node i, Bi can be expressed in the actual coordinate system
	 * by: 
	 *       Bi=[ dN/dx           0    ]
	 *          [   0           dN/dy  ]
	 *          [ 1/2*dN/dy  1/2*dN/dx ]
	 * where N is the interpolation function for node i.
	 *
	 * We assume B has been allocated already, of size: 3x(NDOF2*numnodes)
	 */

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

	/*Get nodal functions derivatives*/
	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);

	/*Build B: */
	for(int i=0;i<numnodes;i++){
		B[NDOF2*numnodes*0+NDOF2*i+0] = dbasis[0*numnodes+i];
		B[NDOF2*numnodes*0+NDOF2*i+1] = 0.;
		B[NDOF2*numnodes*1+NDOF2*i+0] = 0.;
		B[NDOF2*numnodes*1+NDOF2*i+1] = dbasis[1*numnodes+i];
		B[NDOF2*numnodes*2+NDOF2*i+0] = .5*dbasis[1*numnodes+i];
		B[NDOF2*numnodes*2+NDOF2*i+1] = .5*dbasis[0*numnodes+i];
	}

	/*Clean-up*/
	xDelete<IssmDouble>(dbasis);
}
/*}}}*/
/*FUNCTION TriaRef::GetBSSAFS {{{*/
void TriaRef::GetBSSAFS(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss){

	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
	 * For node i, Bi can be expressed in the actual coordinate system
	 * by: 
	 *       Bi=[   dN/dx         0     ]
	 *          [       0       dN/dy   ]
	 *          [  1/2*dN/dy  1/2*dN/dx ]
	 * where N is the interpolation function for node i.
	 *
	 * We assume B has been allocated already, of size: 3x(NDOF2*numnodes)
	 */

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

	/*Get nodal functions derivatives*/
	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);

	/*Build B': */
	for(int i=0;i<numnodes;i++){
		B[NDOF2*numnodes*0+NDOF2*i+0] = dbasis[0*numnodes+i];
		B[NDOF2*numnodes*0+NDOF2*i+1] = 0.;
		B[NDOF2*numnodes*1+NDOF2*i+0] = 0.;
		B[NDOF2*numnodes*1+NDOF2*i+1] = dbasis[1*numnodes+i];
		B[NDOF2*numnodes*2+NDOF2*i+0] = 0.5*dbasis[1*numnodes+i];
		B[NDOF2*numnodes*2+NDOF2*i+1] = 0.5*dbasis[0*numnodes+i];
	}

	/*Clean-up*/
	xDelete<IssmDouble>(dbasis);
}
/*}}}*/
/*FUNCTION TriaRef::GetSegmentBFlux{{{*/
void TriaRef::GetSegmentBFlux(IssmDouble* B,GaussTria* gauss, int index1,int index2){
	/*Compute B  matrix. B=[phi1 phi2 -phi3 -phi4]
	 *
	 * and phi1=phi3 phi2=phi4
	 *
	 * We assume B has been allocated already, of size: 1x4
	 */

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

	/*Get nodal functions*/
	IssmDouble* basis=xNew<IssmDouble>(numnodes);
	GetNodalFunctions(basis,gauss);

	/*Build B for this segment*/
	B[0] = +basis[index1];
	B[1] = +basis[index2];
	B[2] = -basis[index1];
	B[3] = -basis[index2];

	/*Clean-up*/
	xDelete<IssmDouble>(basis);
}
/*}}}*/
/*FUNCTION TriaRef::GetSegmentBprimeFlux{{{*/
void TriaRef::GetSegmentBprimeFlux(IssmDouble* Bprime,GaussTria* gauss, int index1,int index2){
	/*Compute Bprime  matrix. Bprime=[phi1 phi2 phi3 phi4]
	 *
	 * and phi1=phi3 phi2=phi4
	 *
	 * We assume Bprime has been allocated already, of size: 1x4
	 */

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

	/*Get nodal functions*/
	IssmDouble* basis=xNew<IssmDouble>(numnodes);
	GetNodalFunctions(basis,gauss);

	/*Build B'*/
	Bprime[0] = basis[index1];
	Bprime[1] = basis[index2];
	Bprime[2] = basis[index1];
	Bprime[3] = basis[index2];

	/*Clean-up*/
	xDelete<IssmDouble>(basis);
}
/*}}}*/
/*FUNCTION TriaRef::GetBMasstransport{{{*/
void TriaRef::GetBMasstransport(IssmDouble* B, IssmDouble* xyz_list, GaussTria* gauss){
	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
	 * For node i, Bi can be expressed in the actual coordinate system
	 * by: 
	 *       Bi=[ N ]
	 *          [ N ]
	 * where N is the interpolation function for node i.
	 *
	 * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes)
	 */

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

	/*Get nodal functions*/
	IssmDouble* basis=xNew<IssmDouble>(numnodes);
	GetNodalFunctions(basis,gauss);

	/*Build B: */
	for(int i=0;i<numnodes;i++){
		B[numnodes*0+i] = basis[i];
		B[numnodes*1+i] = basis[i];
	}

	/*Clean-up*/
	xDelete<IssmDouble>(basis);
}
/*}}}*/
/*FUNCTION TriaRef::GetBprimeSSA {{{*/
void TriaRef::GetBprimeSSA(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss){

	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
	 * For node i, Bi' can be expressed in the actual coordinate system
	 * by: 
	 *       Bi_prime=[ 2*dN/dx    dN/dy ]
	 *                [   dN/dx  2*dN/dy ]
	 *                [   dN/dy    dN/dx ]
	 * where hNis the interpolation function for node i.
	 *
	 * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
	 */

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

	/*Get nodal functions derivatives*/
	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);

	/*Build B': */
	for(int i=0;i<numnodes;i++){
		Bprime[NDOF2*numnodes*0+NDOF2*i+0] = 2.*dbasis[0*numnodes+i];
		Bprime[NDOF2*numnodes*0+NDOF2*i+1] =    dbasis[1*numnodes+i];
		Bprime[NDOF2*numnodes*1+NDOF2*i+0] =    dbasis[0*numnodes+i];
		Bprime[NDOF2*numnodes*1+NDOF2*i+1] = 2.*dbasis[1*numnodes+i];
		Bprime[NDOF2*numnodes*2+NDOF2*i+0] =    dbasis[1*numnodes+i];
		Bprime[NDOF2*numnodes*2+NDOF2*i+1] =    dbasis[0*numnodes+i];
	}

	/*Clean-up*/
	xDelete<IssmDouble>(dbasis);
}
/*}}}*/
/*FUNCTION TriaRef::GetBprimeSSAFS {{{*/
void TriaRef::GetBprimeSSAFS(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss){
	/*Compute Bprime  matrix. Bprime=[Bprime1 Bprime2 Bprime3] where Bprimei is of size 3*NDOF2. 
	 * For node i, Bprimei can be expressed in the actual coordinate system
	 * by: 
	 *       Bprimei=[  dN/dx    0   ]
	 *               [    0    dN/dy ]
	 *               [  dN/dy  dN/dx ]
	 N               [  dN/dx  dN/dy ]
	 * where N is the interpolation function for node i.
	 *
	 * We assume Bprime has been allocated already, of size: 3x(NDOF2*numnodes)
	 */

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

	/*Get nodal functions*/
	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);

	/*Build Bprime: */
	for(int i=0;i<numnodes;i++){
		Bprime[NDOF2*numnodes*0+NDOF2*i+0] = dbasis[0*numnodes+i];
		Bprime[NDOF2*numnodes*0+NDOF2*i+1] = 0.;
		Bprime[NDOF2*numnodes*1+NDOF2*i+0] = 0.;
		Bprime[NDOF2*numnodes*1+NDOF2*i+1] = dbasis[1*numnodes+i];
		Bprime[NDOF2*numnodes*2+NDOF2*i+0] = dbasis[1*numnodes+i];
		Bprime[NDOF2*numnodes*2+NDOF2*i+1] = dbasis[0*numnodes+i];
		Bprime[NDOF2*numnodes*3+NDOF2*i+0] = dbasis[0*numnodes+i];
		Bprime[NDOF2*numnodes*3+NDOF2*i+1] = dbasis[1*numnodes+i];
	}

	/*Clean-up*/
	xDelete<IssmDouble>(dbasis);
}
/*}}}*/
/*FUNCTION TriaRef::GetBprimeMasstransport{{{*/
void TriaRef::GetBprimeMasstransport(IssmDouble* Bprime, IssmDouble* xyz_list, GaussTria* gauss){
	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
	 * For node i, Bi' can be expressed in the actual coordinate system
	 * by: 
	 *       Bi_prime=[ dN/dx ]
	 *                [ dN/dy ]
	 * where N is the interpolation function for node i.
	 *
	 * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
	 */

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

	/*Get nodal functions derivatives*/
	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);

	/*Build B': */
	for(int i=0;i<numnodes;i++){
		Bprime[numnodes*0+i] = dbasis[0*numnodes+i];
		Bprime[numnodes*1+i] = dbasis[1*numnodes+i];
	}

	/*Clean-up*/
	xDelete<IssmDouble>(dbasis);
}
/*}}}*/
/*FUNCTION TriaRef::GetBSSAFriction{{{*/
void TriaRef::GetBSSAFriction(IssmDouble* B, IssmDouble* xyz_list,GaussTria* gauss){
	/*Compute B  matrix. B=[B1 B2 B3] where Bi is square and of size 2. 
	 * For node i, Bi can be expressed in the actual coordinate system
	 * by: 
	 *                 Bi=[ N   0 ]
	 *                    [ 0   N ]
	 * where N is the interpolation function for node i.
	 *
	 * We assume B has been allocated already, of size: 2 x (numdof*numnodes)
	 */

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

	/*Get nodal functions derivatives*/
	IssmDouble* basis=xNew<IssmDouble>(numnodes);
	GetNodalFunctions(basis,gauss);

	/*Build L: */
	for(int i=0;i<numnodes;i++){
		B[2*numnodes*0+2*i+0] = basis[i];
		B[2*numnodes*0+2*i+1] = 0.;
		B[2*numnodes*1+2*i+0] = 0.;
		B[2*numnodes*1+2*i+1] = basis[i];
	}

	/*Clean-up*/
	xDelete<IssmDouble>(basis);
}
/*}}}*/
/*FUNCTION TriaRef::GetJacobian{{{*/
void TriaRef::GetJacobian(IssmDouble* J, IssmDouble* xyz_list,GaussTria* gauss){
	/*The Jacobian is constant over the element, discard the gaussian points. 
	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/

	IssmDouble x1 = xyz_list[3*0+0];
	IssmDouble y1 = xyz_list[3*0+1];
	IssmDouble x2 = xyz_list[3*1+0];
	IssmDouble y2 = xyz_list[3*1+1];
	IssmDouble x3 = xyz_list[3*2+0];
	IssmDouble y3 = xyz_list[3*2+1];

	J[2*0+0] = 0.5*(x2-x1);
	J[2*1+0] = SQRT3/6.0*(2*x3-x1-x2);
	J[2*0+1] = 0.5*(y2-y1);
	J[2*1+1] = SQRT3/6.0*(2*y3-y1-y2);
}
/*}}}*/
/*FUNCTION TriaRef::GetSegmentJacobianDeterminant{{{*/
void TriaRef::GetSegmentJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,GaussTria* gauss){
	/*The Jacobian determinant is constant over the element, discard the gaussian points. 
	 * J is assumed to have been allocated*/

	IssmDouble x1 = xyz_list[3*0+0];
	IssmDouble y1 = xyz_list[3*0+1];
	IssmDouble x2 = xyz_list[3*1+0];
	IssmDouble y2 = xyz_list[3*1+1];

	*Jdet = .5*sqrt(pow(x2-x1,2) + pow(y2-y1,2));
	if(*Jdet<0) _error_("negative jacobian determinant!");

}
/*}}}*/
/*FUNCTION TriaRef::GetJacobianDeterminant{{{*/
void TriaRef::GetJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,GaussTria* gauss){
	/*The Jacobian determinant is constant over the element, discard the gaussian points. 
	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
	IssmDouble J[2][2];

	/*Get Jacobian*/
	GetJacobian(&J[0][0],xyz_list,gauss);

	/*Get Determinant*/
	Matrix2x2Determinant(Jdet,&J[0][0]);
	if(*Jdet<0) _error_("negative jacobian determinant!");

}
/*}}}*/
/*FUNCTION TriaRef::GetJacobianInvert{{{*/
void TriaRef::GetJacobianInvert(IssmDouble*  Jinv, IssmDouble* xyz_list,GaussTria* gauss){

	/*Jacobian*/
	IssmDouble J[2][2];

	/*Call Jacobian routine to get the jacobian:*/
	GetJacobian(&J[0][0], xyz_list, gauss);

	/*Invert Jacobian matrix: */
	Matrix2x2Invert(Jinv,&J[0][0]);

}
/*}}}*/
/*FUNCTION TriaRef::GetNodalFunctions{{{*/
void TriaRef::GetNodalFunctions(IssmDouble* basis,GaussTria* gauss){
	/*This routine returns the values of the nodal functions  at the gaussian point.*/

	_assert_(basis);

	switch(this->element_type){
		case P1Enum: case P1DGEnum:
			basis[0]=gauss->coord1;
			basis[1]=gauss->coord2;
			basis[2]=gauss->coord3;
			return;
		case P1bubbleEnum: case P1bubblecondensedEnum:
			/*Corner nodes*/
			basis[0]=gauss->coord1;
			basis[1]=gauss->coord2;
			basis[2]=gauss->coord3;
			/*bubble*/
			basis[3]=27.*gauss->coord1*gauss->coord2*gauss->coord3;
			return;
		case P2Enum:
			/*Corner nodes*/
			basis[0]=gauss->coord1*(2.*gauss->coord1-1.);
			basis[1]=gauss->coord2*(2.*gauss->coord2-1.);
			basis[2]=gauss->coord3*(2.*gauss->coord3-1.);
			/*Mid-sides*/
			basis[3]=4.*gauss->coord3*gauss->coord2;
			basis[4]=4.*gauss->coord3*gauss->coord1;
			basis[5]=4.*gauss->coord1*gauss->coord2;
			return;
		default:
			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
	}
}
/*}}}*/
/*FUNCTION TriaRef::GetSegmentNodalFunctions{{{*/
void TriaRef::GetSegmentNodalFunctions(IssmDouble* basis,GaussTria* gauss,int index1,int index2){
	/*This routine returns the values of the nodal functions  at the gaussian point.*/

	_assert_(index1>=0 && index1<3);
	_assert_(index2>=0 && index2<3);

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

	/*Get nodal functions*/
	IssmDouble* triabasis=xNew<IssmDouble>(numnodes);
	GetNodalFunctions(triabasis,gauss);

	switch(this->element_type){
		case P1Enum: case P1DGEnum:
			basis[0]=triabasis[index1];
			basis[1]=triabasis[index2];
			xDelete<IssmDouble>(triabasis);
			return;
		case P1bubbleEnum: case P1bubblecondensedEnum:
			basis[0]=triabasis[index1];
			basis[1]=triabasis[index2];
			xDelete<IssmDouble>(triabasis);
			return;
		case P2Enum:
			_assert_(index2<index1);
			basis[0]=triabasis[index1];
			basis[1]=triabasis[index2];
			basis[2]=triabasis[3+index2-1];
			xDelete<IssmDouble>(triabasis);
			return;
		default:
			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
	}

	/*Clean up*/
	xDelete<IssmDouble>(triabasis);
}
/*}}}*/
/*FUNCTION TriaRef::GetNodalFunctionsDerivatives{{{*/
void TriaRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussTria* gauss){

	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
	 * actual coordinate system): */
	IssmDouble    Jinv[2][2];

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

	/*Get nodal functions derivatives in reference triangle*/
	IssmDouble* dbasis_ref=xNew<IssmDouble>(2*numnodes);
	GetNodalFunctionsDerivativesReference(dbasis_ref,gauss); 

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

	/*Build dbasis: 
	 * [dhi/dx]= Jinv*[dhi/dr]
	 * [dhi/dy]       [dhi/ds]
	 */
	for(int i=0;i<numnodes;i++){
		dbasis[numnodes*0+i] = Jinv[0][0]*dbasis_ref[0*numnodes+i]+Jinv[0][1]*dbasis_ref[1*numnodes+i];
		dbasis[numnodes*1+i] = Jinv[1][0]*dbasis_ref[0*numnodes+i]+Jinv[1][1]*dbasis_ref[1*numnodes+i];
	}

	/*Clean up*/
	xDelete<IssmDouble>(dbasis_ref);

}
/*}}}*/
/*FUNCTION TriaRef::GetNodalFunctionsDerivativesReference{{{*/
void TriaRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussTria* gauss){
	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
	 * natural coordinate system) at the gaussian point. */

	_assert_(dbasis && gauss);

	switch(this->element_type){
		case P1Enum: case P1DGEnum:
			/*Nodal function 1*/
			dbasis[NUMNODESP1*0+0] = -0.5;
			dbasis[NUMNODESP1*1+0] = -SQRT3/6.;
			/*Nodal function 2*/
			dbasis[NUMNODESP1*0+1] = 0.5;
			dbasis[NUMNODESP1*1+1] = -SQRT3/6.;
			/*Nodal function 3*/
			dbasis[NUMNODESP1*0+2] = 0;
			dbasis[NUMNODESP1*1+2] = SQRT3/3.;
			return;
		case P1bubbleEnum: case P1bubblecondensedEnum:
			/*Nodal function 1*/
			dbasis[NUMNODESP1b*0+0] = -0.5;
			dbasis[NUMNODESP1b*1+0] = -SQRT3/6.;
			/*Nodal function 2*/
			dbasis[NUMNODESP1b*0+1] = 0.5;
			dbasis[NUMNODESP1b*1+1] = -SQRT3/6.;
			/*Nodal function 3*/
			dbasis[NUMNODESP1b*0+2] = 0;
			dbasis[NUMNODESP1b*1+2] = SQRT3/3.;
			/*Nodal function 4*/
			dbasis[NUMNODESP1b*0+3] = 27.*(-.5*gauss->coord2*gauss->coord3 + .5*gauss->coord1*gauss->coord3);
			dbasis[NUMNODESP1b*1+3] = 27.*SQRT3*(-1./6.*gauss->coord2*gauss->coord3 - 1./6.*gauss->coord1*gauss->coord3 +1./3.*gauss->coord1*gauss->coord2);
			return;
		case P2Enum:
			/*Nodal function 1*/
			dbasis[NUMNODESP2*0+0] = -2.*gauss->coord1 + 0.5;
			dbasis[NUMNODESP2*1+0] = -2.*SQRT3/3.*gauss->coord1 + SQRT3/6.;
			/*Nodal function 2*/
			dbasis[NUMNODESP2*0+1] = +2.*gauss->coord2 - 0.5;
			dbasis[NUMNODESP2*1+1] = -2.*SQRT3/3.*gauss->coord2 + SQRT3/6.;
			/*Nodal function 3*/
			dbasis[NUMNODESP2*0+2] = 0.;
			dbasis[NUMNODESP2*1+2] = +4.*SQRT3/3.*gauss->coord3 - SQRT3/3.;
			/*Nodal function 4*/
			dbasis[NUMNODESP2*0+3] = +2.*gauss->coord3;
			dbasis[NUMNODESP2*1+3] = +4.*SQRT3/3.*gauss->coord2 - 2.*SQRT3/3.*gauss->coord3;
			/*Nodal function 5*/
			dbasis[NUMNODESP2*0+4] = -2.*gauss->coord3;
			dbasis[NUMNODESP2*1+4] = +4.*SQRT3/3.*gauss->coord1 - 2.*SQRT3/3.*gauss->coord3;
			/*Nodal function 6*/
			dbasis[NUMNODESP2*0+5] = 2.*(gauss->coord1-gauss->coord2);
			dbasis[NUMNODESP2*1+5] = -2.*SQRT3/3.*(gauss->coord1+gauss->coord2);
			return;
		default:
			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
	}

}
/*}}}*/
/*FUNCTION TriaRef::GetInputDerivativeValue{{{*/
void TriaRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussTria* gauss){

	/*From node values of parameter p (plist[0],plist[1],plist[2]), return parameter derivative value at gaussian 
	 * point specified by gauss_basis:
	 *   dp/dx=plist[0]*dh1/dx+plist[1]*dh2/dx+plist[2]*dh3/dx
	 *   dp/dx=plist[0]*dh1/dx+plist[1]*dh2/dx+plist[2]*dh3/dx
	 *
	 * p is a vector already allocated.
	 */

	/*Output*/
	IssmDouble dpx=0.;
	IssmDouble dpy=0.;

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

	/*Get nodal functions derivatives*/
	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);

	/*Calculate parameter for this Gauss point*/
	for(int i=0;i<numnodes;i++) dpx += dbasis[0*numnodes+i]*plist[i];
	for(int i=0;i<numnodes;i++) dpy += dbasis[1*numnodes+i]*plist[i];

	/*Assign values*/
	xDelete<IssmDouble>(dbasis);
	p[0]=dpx;
	p[1]=dpy;

}
/*}}}*/
/*FUNCTION TriaRef::GetInputValue{{{*/
void TriaRef::GetInputValue(IssmDouble* p, IssmDouble* plist, GaussTria* gauss){

	/*Output*/
	IssmDouble value =0.;

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

	/*Get nodal functions*/
	IssmDouble* basis=xNew<IssmDouble>(numnodes);
	GetNodalFunctions(basis, gauss);

	/*Calculate parameter for this Gauss point*/
	for(int i=0;i<numnodes;i++) value += basis[i]*plist[i];

	/*Assign output pointer*/
	xDelete<IssmDouble>(basis);
	*p = value;
}
/*}}}*/
/*FUNCTION TriaRef::NumberofNodes{{{*/
int TriaRef::NumberofNodes(void){

	switch(this->element_type){
		case P1Enum:                return NUMNODESP1;
		case P1DGEnum:              return NUMNODESP1;
		case P1bubbleEnum:          return NUMNODESP1b;
		case P1bubblecondensedEnum: return NUMNODESP1b;
		case P2Enum:                return NUMNODESP2;
		default: _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
	}

	return -1;
}
/*}}}*/
