/*!\file PentaRef.c
 * \brief: implementation of the PentaRef 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    6
#define NUMNODESP1_2d 3
#define NUMNODESP1b   7
#define NUMNODESP1xP2 9
#define NUMNODESP2xP1 12
#define NUMNODESP2    15

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

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

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

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

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

/*Reference Element numerics*/
/*FUNCTION PentaRef::GetBSSAHO {{{*/
void PentaRef::GetBSSAHO(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){
	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. 
	 * For node i, Bi can be expressed in the actual coordinate system
	 * by: 
	 *       Bi=[ dh/dx          0      ]
	 *          [   0           dh/dy   ]
	 *          [ 1/2*dh/dy  1/2*dh/dx  ]
	 * where h is the interpolation function for node i.
	 *
	 * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1)
	 */

	IssmDouble dbasis[3][NUMNODESP1];

	/*Get dbasis in actual coordinate system: */
	GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss);

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

		B[NDOF2*NUMNODESP1*1+NDOF2*i+0] = 0.;
		B[NDOF2*NUMNODESP1*1+NDOF2*i+1] = dbasis[1][i];

		B[NDOF2*NUMNODESP1*2+NDOF2*i+0] = .5*dbasis[1][i];
		B[NDOF2*NUMNODESP1*2+NDOF2*i+1] = .5*dbasis[0][i];
	}
}
/*}}}*/
/*FUNCTION PentaRef::GetBSSAFS{{{*/
void PentaRef::GetBSSAFS(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){
	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. 
	 * For node i, Bi can be expressed in the actual coordinate system
	 * by: 
	 *       Bi=[ dh/dx          0       0   0 ]
	 *          [   0           dh/dy    0   0 ]
	 *          [ 1/2*dh/dy  1/2*dh/dx   0   0 ]
	 *          [   0            0       0   h ]
	 * where h is the interpolation function for node i.
	 *
	 * We assume B has been allocated already, of size: 5x(NDOF2*NUMNODESP1)
	 */

	int i;
	IssmDouble dbasismini[3][NUMNODESP1b];
	IssmDouble basis[NUMNODESP1];

	/*Get dbasis in actual coordinate system: */
	GetNodalFunctionsMINIDerivatives(&dbasismini[0][0],xyz_list, gauss);
	GetNodalFunctionsP1(basis,gauss);

	/*Build B: */
	for(i=0;i<NUMNODESP1;i++){
		B[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*i+0] = dbasismini[0][i];
		B[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*i+1] = 0.;
		B[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*i+2] = 0.;
		B[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*i+0] = 0.;
		B[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*i+1] = dbasismini[1][i];
		B[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*i+2] = 0.;
		B[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*i+0] = 0.5*dbasismini[1][i];
		B[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*i+1] = 0.5*dbasismini[0][i];
		B[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*i+2] = 0.;
		B[(NDOF3*NUMNODESP1b+NUMNODESP1)*3+NDOF3*i+0] = 0.;
		B[(NDOF3*NUMNODESP1b+NUMNODESP1)*3+NDOF3*i+1] = 0.;
		B[(NDOF3*NUMNODESP1b+NUMNODESP1)*3+NDOF3*i+2] = 0.;
	}
	for(i=0;i<1;i++){
		B[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*(NUMNODESP1+i)+0] = 0.;
		B[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*(NUMNODESP1+i)+1] = 0.;
		B[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*(NUMNODESP1+i)+2] = 0.;
		B[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*(NUMNODESP1+i)+0] = 0.;
		B[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*(NUMNODESP1+i)+1] = 0.;
		B[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*(NUMNODESP1+i)+2] = 0.;
		B[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*(NUMNODESP1+i)+0] = 0.;
		B[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*(NUMNODESP1+i)+1] = 0.;
		B[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*(NUMNODESP1+i)+2] = 0.;
		B[(NDOF3*NUMNODESP1b+NUMNODESP1)*3+NDOF3*(NUMNODESP1+i)+0] = 0.;
		B[(NDOF3*NUMNODESP1b+NUMNODESP1)*3+NDOF3*(NUMNODESP1+i)+1] = 0.;
		B[(NDOF3*NUMNODESP1b+NUMNODESP1)*3+NDOF3*(NUMNODESP1+i)+2] = 0.;
	}

	for(i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
		B[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NUMNODESP1b*NDOF3+i] = 0;
		B[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NUMNODESP1b*NDOF3+i] = 0;
		B[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NUMNODESP1b*NDOF3+i] = 0;
		B[(NDOF3*NUMNODESP1b+NUMNODESP1)*3+NUMNODESP1b*NDOF3+i] = basis[i];
	}
}
/*}}}*/
/*FUNCTION PentaRef::GetBHO {{{*/
void PentaRef::GetBHO(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){
	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. 
	 * For node i, Bi can be expressed in the actual coordinate system
	 * by: 
	 *       Bi=[ dh/dx          0      ]
	 *          [   0           dh/dy   ]
	 *          [ 1/2*dh/dy  1/2*dh/dx  ]
	 *          [ 1/2*dh/dz      0      ]
	 *          [  0         1/2*dh/dz  ]
	 * where h is the interpolation function for node i.
	 *
	 * We assume B has been allocated already, of size: 5x(NDOF2*numnodes)
	 */

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

	/*Get nodal functions derivatives*/
	IssmDouble* dbasis=xNew<IssmDouble>(3*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];
		B[NDOF2*numnodes*3+NDOF2*i+0] = .5*dbasis[2*numnodes+i];
		B[NDOF2*numnodes*3+NDOF2*i+1] = 0.;
		B[NDOF2*numnodes*4+NDOF2*i+0] = 0.;
		B[NDOF2*numnodes*4+NDOF2*i+1] = .5*dbasis[2*numnodes+i];
	}

	/*Clean-up*/
	xDelete<IssmDouble>(dbasis);
}
/*}}}*/
/*FUNCTION PentaRef::GetBprimeHO {{{*/
void PentaRef::GetBprimeHO(IssmDouble* B,IssmDouble* xyz_list,GaussPenta* gauss){
	/*Compute B  prime matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF2. 
	 * For node i, Bi can be expressed in the actual coordinate system
	 * by: 
	 *       Bi=[ 2*dh/dx     dh/dy   ]
	 *          [   dh/dx    2*dh/dy  ]
	 *          [ dh/dy      dh/dx    ]
	 *          [ dh/dz         0     ]
	 *          [  0         dh/dz    ]
	 * where h is the interpolation function for node i.
	 *
	 * We assume B has been allocated already, of size: 5x(NDOF2*numnodes)
	 */

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

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

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

	/*Clean-up*/
	xDelete<IssmDouble>(dbasis);
}
/*}}}*/
/*FUNCTION PentaRef::GetBprimeSSAFS{{{*/
void PentaRef::GetBprimeSSAFS(IssmDouble* Bprime, IssmDouble* xyz_list, GaussPenta* gauss){
	/*Compute Bprime  matrix. Bprime=[Bprime1 Bprime2 Bprime3 Bprime4 Bprime5 Bprime6] where Bprimei is of size 5*NDOF2. 
	 * For node i, Bprimei can be expressed in the actual coordinate system
	 * by: 
	 *       Bprimei=[ 2*dh/dx    dh/dy   0   0 ]
	 *               [  dh/dx    2*dh/dy  0   0 ]
	 *               [  dh/dy     dh/dx   0   0 ]
	 * where h is the interpolation function for node i.
	 *
	 * We assume Bprime has been allocated already, of size: 5x(NDOF2*NUMNODESP1)
	 */

	int    i;
	IssmDouble dbasismini[3][NUMNODESP1b];

	/*Get dbasis in actual coordinate system: */
	GetNodalFunctionsMINIDerivatives(&dbasismini[0][0],xyz_list, gauss);

	/*Build Bprime: */
	for(i=0;i<NUMNODESP1;i++){
		Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*i+0] = 2.*dbasismini[0][i];
		Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*i+1] = dbasismini[1][i];
		Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*i+2] = 0.;
		Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*i+0] = dbasismini[0][i];
		Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*i+1] = 2.*dbasismini[1][i];
		Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*i+2] = 0.;
		Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*i+0] = dbasismini[1][i];
		Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*i+1] = dbasismini[0][i];
		Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*i+2] = 0.;
	}

	for(i=0;i<1;i++){ //Add zeros for the bubble function
		Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*(NUMNODESP1+i)+0] = 0.;
		Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*(NUMNODESP1+i)+1] = 0.;
		Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NDOF3*(NUMNODESP1+i)+2] = 0.;
		Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*(NUMNODESP1+i)+0] = 0.;
		Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*(NUMNODESP1+i)+1] = 0.;
		Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NDOF3*(NUMNODESP1+i)+2] = 0.;
		Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*(NUMNODESP1+i)+0] = 0.;
		Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*(NUMNODESP1+i)+1] = 0.;
		Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NDOF3*(NUMNODESP1+i)+2] = 0.;
	}

	for(i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
		Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*0+NUMNODESP1b*NDOF3+i] = 0.;
		Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*1+NUMNODESP1b*NDOF3+i] = 0.;
		Bprime[(NDOF3*NUMNODESP1b+NUMNODESP1)*2+NUMNODESP1b*NDOF3+i] = 0.;
	}

}
/*}}}*/
/*FUNCTION PentaRef::GetBFSstrainrate {{{*/
void PentaRef::GetBFSstrainrate(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){

	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*NDOF4. 
	 * For node i, Bi can be expressed in the actual coordinate system
	 * by: 		Bi=[ dh/dx          0              0     ]
	 *					[   0           dh/dy           0     ]
	 *					[   0             0           dh/dy   ]
	 *					[ 1/2*dh/dy    1/2*dh/dx        0     ]
	 *					[ 1/2*dh/dz       0         1/2*dh/dx ]
	 *					[   0          1/2*dh/dz    1/2*dh/dy ]
	 *	where h is the interpolation function for node i.
	 *	Same thing for Bb except the last column that does not exist.
	 */

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

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

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

		B[3*numnodes*1+3*i+0] = 0.;
		B[3*numnodes*1+3*i+1] = dbasis[1*numnodes+i+0];
		B[3*numnodes*1+3*i+2] = 0.;

		B[3*numnodes*2+3*i+0] = 0.;
		B[3*numnodes*2+3*i+1] = 0.;
		B[3*numnodes*2+3*i+2] = dbasis[2*numnodes+i+0];

		B[3*numnodes*3+3*i+0] = .5*dbasis[1*numnodes+i+0];
		B[3*numnodes*3+3*i+1] = .5*dbasis[0*numnodes+i+0];
		B[3*numnodes*3+3*i+2] = 0.;

		B[3*numnodes*4+3*i+0] = .5*dbasis[2*numnodes+i+0];
		B[3*numnodes*4+3*i+1] = 0.;
		B[3*numnodes*4+3*i+2] = .5*dbasis[0*numnodes+i+0];

		B[3*numnodes*5+3*i+0] = 0.;
		B[3*numnodes*5+3*i+1] = .5*dbasis[2*numnodes+i+0];
		B[3*numnodes*5+3*i+2] = .5*dbasis[1*numnodes+i+0];
	}

	/*Clean up*/
	xDelete<IssmDouble>(dbasis);
}
/*}}}*/
/*FUNCTION PentaRef::GetBFS {{{*/
void PentaRef::GetBFS(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){

	/*Compute B  matrix. B=[Bv1 Bv2 ... Bp1 Bp2 ...] where Bvi is of size 3*NDOF3. 
	 * For node i, Bvi can be expressed in the actual coordinate system
	 * by: 	   Bvi=[ dh/dx          0              0      ]
	 *					[   0           dh/dy           0      ]
	 *					[   0             0           dh/dy    ]
	 *					[ 1/2*dh/dy    1/2*dh/dx        0      ]
	 *					[ 1/2*dh/dz       0         1/2*dh/dx  ]
	 *					[   0          1/2*dh/dz    1/2*dh/dy  ]
	 *					[   0             0             0      ]
	 *					[ dh/dx         dh/dy         dh/dz    ]
	 *
	 * by: 	  Bpi=[ 0 ]
	 *					[ 0 ]
	 *					[ 0 ]
	 *					[ 0 ]
	 *					[ 0 ]
	 *					[ 0 ]
	 *					[ h ]
	 *					[ 0 ]
	 *	where h is the interpolation function for node i.
	 *	Same thing for Bb except the last column that does not exist.
	 */

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

	/*Get nodal functions derivatives*/
	IssmDouble* vdbasis=xNew<IssmDouble>(3*vnumnodes);
	IssmDouble* pbasis =xNew<IssmDouble>(pnumnodes);
	GetNodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);
	GetNodalFunctionsPressure(pbasis,gauss);

	/*Build B: */
	for(int i=0;i<vnumnodes;i++){
		B[(3*vnumnodes+pnumnodes)*0+3*i+0] = vdbasis[0*vnumnodes+i];
		B[(3*vnumnodes+pnumnodes)*0+3*i+1] = 0.;
		B[(3*vnumnodes+pnumnodes)*0+3*i+2] = 0.;
		B[(3*vnumnodes+pnumnodes)*1+3*i+0] = 0.;
		B[(3*vnumnodes+pnumnodes)*1+3*i+1] = vdbasis[1*vnumnodes+i];
		B[(3*vnumnodes+pnumnodes)*1+3*i+2] = 0.;
		B[(3*vnumnodes+pnumnodes)*2+3*i+0] = 0.;
		B[(3*vnumnodes+pnumnodes)*2+3*i+1] = 0.;
		B[(3*vnumnodes+pnumnodes)*2+3*i+2] = vdbasis[2*vnumnodes+i];
		B[(3*vnumnodes+pnumnodes)*3+3*i+0] = .5*vdbasis[1*vnumnodes+i];
		B[(3*vnumnodes+pnumnodes)*3+3*i+1] = .5*vdbasis[0*vnumnodes+i];
		B[(3*vnumnodes+pnumnodes)*3+3*i+2] = 0.;
		B[(3*vnumnodes+pnumnodes)*4+3*i+0] = .5*vdbasis[2*vnumnodes+i];
		B[(3*vnumnodes+pnumnodes)*4+3*i+1] = 0.;
		B[(3*vnumnodes+pnumnodes)*4+3*i+2] = .5*vdbasis[0*vnumnodes+i];
		B[(3*vnumnodes+pnumnodes)*5+3*i+0] = 0.;
		B[(3*vnumnodes+pnumnodes)*5+3*i+1] = .5*vdbasis[2*vnumnodes+i];
		B[(3*vnumnodes+pnumnodes)*5+3*i+2] = .5*vdbasis[1*vnumnodes+i];
		B[(3*vnumnodes+pnumnodes)*6+3*i+0] = 0.;
		B[(3*vnumnodes+pnumnodes)*6+3*i+1] = 0.;
		B[(3*vnumnodes+pnumnodes)*6+3*i+2] = 0.;
		B[(3*vnumnodes+pnumnodes)*7+3*i+0] = vdbasis[0*vnumnodes+i];
		B[(3*vnumnodes+pnumnodes)*7+3*i+1] = vdbasis[1*vnumnodes+i];
		B[(3*vnumnodes+pnumnodes)*7+3*i+2] = vdbasis[2*vnumnodes+i];
	}
	for(int i=0;i<pnumnodes;i++){
		B[(3*vnumnodes+pnumnodes)*0+(3*vnumnodes)+i] = 0.;
		B[(3*vnumnodes+pnumnodes)*1+(3*vnumnodes)+i] = 0.;
		B[(3*vnumnodes+pnumnodes)*2+(3*vnumnodes)+i] = 0.;
		B[(3*vnumnodes+pnumnodes)*3+(3*vnumnodes)+i] = 0.;
		B[(3*vnumnodes+pnumnodes)*4+(3*vnumnodes)+i] = 0.;
		B[(3*vnumnodes+pnumnodes)*5+(3*vnumnodes)+i] = 0.;
		B[(3*vnumnodes+pnumnodes)*6+(3*vnumnodes)+i] = pbasis[i];
		B[(3*vnumnodes+pnumnodes)*7+(3*vnumnodes)+i] = 0.;
	}

	/*Clean up*/
	xDelete<IssmDouble>(vdbasis);
	xDelete<IssmDouble>(pbasis);
}
/*}}}*/
/*FUNCTION PentaRef::GetBFSGLS {{{*/
void PentaRef::GetBFSGLS(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){

	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 3*NDOF4. 
	 * For node i, Bi can be expressed in the actual coordinate system
	 * by: 		Bi=[ dh/dx          0              0       0  ]
	 *					[   0           dh/dy           0       0  ]
	 *					[   0             0           dh/dy     0  ]
	 *					[ 1/2*dh/dy    1/2*dh/dx        0       0  ]
	 *					[ 1/2*dh/dz       0         1/2*dh/dx   0  ]
	 *					[   0          1/2*dh/dz    1/2*dh/dy   0  ]
	 *					[   0             0             0       h  ]
	 *					[ dh/dx         dh/dy         dh/dz     0  ]
	 *	where h is the interpolation function for node i.
	 *	Same thing for Bb except the last column that does not exist.
	 */

	int i;
	IssmDouble dbasis[3][NUMNODESP1];
	IssmDouble basis[NUMNODESP1];

	/*Get dbasismini in actual coordinate system: */
	GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss);
	GetNodalFunctionsP1(&basis[0], gauss);

	/*Build B: */
	for(i=0;i<NUMNODESP1;i++){
		B[(NDOF4*NUMNODESP1)*0+NDOF4*i+0] = dbasis[0][i];
		B[(NDOF4*NUMNODESP1)*0+NDOF4*i+1] = 0.;
		B[(NDOF4*NUMNODESP1)*0+NDOF4*i+2] = 0.;
		B[(NDOF4*NUMNODESP1)*1+NDOF4*i+0] = 0.;
		B[(NDOF4*NUMNODESP1)*1+NDOF4*i+1] = dbasis[1][i];
		B[(NDOF4*NUMNODESP1)*1+NDOF4*i+2] = 0.;
		B[(NDOF4*NUMNODESP1)*2+NDOF4*i+0] = 0.;
		B[(NDOF4*NUMNODESP1)*2+NDOF4*i+1] = 0.;
		B[(NDOF4*NUMNODESP1)*2+NDOF4*i+2] = dbasis[2][i];
		B[(NDOF4*NUMNODESP1)*3+NDOF4*i+0] = .5*dbasis[1][i];
		B[(NDOF4*NUMNODESP1)*3+NDOF4*i+1] = .5*dbasis[0][i];
		B[(NDOF4*NUMNODESP1)*3+NDOF4*i+2] = 0.;
		B[(NDOF4*NUMNODESP1)*4+NDOF4*i+0] = .5*dbasis[2][i];
		B[(NDOF4*NUMNODESP1)*4+NDOF4*i+1] = 0.;
		B[(NDOF4*NUMNODESP1)*4+NDOF4*i+2] = .5*dbasis[0][i];
		B[(NDOF4*NUMNODESP1)*5+NDOF4*i+0] = 0.;
		B[(NDOF4*NUMNODESP1)*5+NDOF4*i+1] = .5*dbasis[2][i];
		B[(NDOF4*NUMNODESP1)*5+NDOF4*i+2] = .5*dbasis[1][i];
		B[(NDOF4*NUMNODESP1)*6+NDOF4*i+0] = 0.;
		B[(NDOF4*NUMNODESP1)*6+NDOF4*i+1] = 0.;
		B[(NDOF4*NUMNODESP1)*6+NDOF4*i+2] = 0.;
		B[(NDOF4*NUMNODESP1)*7+NDOF4*i+0] = dbasis[0][i];
		B[(NDOF4*NUMNODESP1)*7+NDOF4*i+1] = dbasis[1][i];
		B[(NDOF4*NUMNODESP1)*7+NDOF4*i+2] = dbasis[2][i];
	}

	for(i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
		B[(NDOF4*NUMNODESP1)*0+NDOF4*i+3] = 0.;
		B[(NDOF4*NUMNODESP1)*1+NDOF4*i+3] = 0.;
		B[(NDOF4*NUMNODESP1)*2+NDOF4*i+3] = 0.;
		B[(NDOF4*NUMNODESP1)*3+NDOF4*i+3] = 0.;
		B[(NDOF4*NUMNODESP1)*4+NDOF4*i+3] = 0.;
		B[(NDOF4*NUMNODESP1)*5+NDOF4*i+3] = 0.;
		B[(NDOF4*NUMNODESP1)*6+NDOF4*i+3] = basis[i];
		B[(NDOF4*NUMNODESP1)*7+NDOF4*i+3] = 0.;
	}

}
/*}}}*/
/*FUNCTION PentaRef::GetBprimeFS {{{*/
void PentaRef::GetBprimeFS(IssmDouble* B_prime, IssmDouble* xyz_list, GaussPenta* gauss){
	/*	Compute B'  matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2. 
	 *	For node i, Bi' can be expressed in the actual coordinate system
	 *	by: 
	 *			Bvi' = [  dh/dx   0          0    ]
	 *					 [   0      dh/dy      0    ]
	 *					 [   0      0         dh/dz ]
	 *					 [  dh/dy   dh/dx      0    ]
	 *					 [  dh/dz   0        dh/dx  ]
	 *					 [   0      dh/dz    dh/dy  ]
	 *					 [  dh/dx   dh/dy    dh/dz  ]
	 *					 [   0      0          0    ]
	 *
	 * by: 	  Bpi=[ 0 ]
	 *					[ 0 ]
	 *					[ 0 ]
	 *					[ 0 ]
	 *					[ 0 ]
	 *					[ 0 ]
	 *					[ 0 ]
	 *					[ h ]
	 *	where h is the interpolation function for node i.
	 */

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

	/*Get nodal functions derivatives*/
	IssmDouble* vdbasis=xNew<IssmDouble>(3*vnumnodes);
	IssmDouble* pbasis =xNew<IssmDouble>(pnumnodes);
	GetNodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);
	GetNodalFunctionsPressure(pbasis,gauss);

	/*Build B_prime: */
	for(int i=0;i<vnumnodes;i++){
		B_prime[(3*vnumnodes+pnumnodes)*0+3*i+0] = vdbasis[0*vnumnodes+i];
		B_prime[(3*vnumnodes+pnumnodes)*0+3*i+1] = 0.;
		B_prime[(3*vnumnodes+pnumnodes)*0+3*i+2] = 0.;
		B_prime[(3*vnumnodes+pnumnodes)*1+3*i+0] = 0.;
		B_prime[(3*vnumnodes+pnumnodes)*1+3*i+1] = vdbasis[1*vnumnodes+i];
		B_prime[(3*vnumnodes+pnumnodes)*1+3*i+2] = 0.;
		B_prime[(3*vnumnodes+pnumnodes)*2+3*i+0] = 0.;
		B_prime[(3*vnumnodes+pnumnodes)*2+3*i+1] = 0.;
		B_prime[(3*vnumnodes+pnumnodes)*2+3*i+2] = vdbasis[2*vnumnodes+i];
		B_prime[(3*vnumnodes+pnumnodes)*3+3*i+0] = vdbasis[1*vnumnodes+i];
		B_prime[(3*vnumnodes+pnumnodes)*3+3*i+1] = vdbasis[0*vnumnodes+i];
		B_prime[(3*vnumnodes+pnumnodes)*3+3*i+2] = 0.;
		B_prime[(3*vnumnodes+pnumnodes)*4+3*i+0] = vdbasis[2*vnumnodes+i];
		B_prime[(3*vnumnodes+pnumnodes)*4+3*i+1] = 0.;
		B_prime[(3*vnumnodes+pnumnodes)*4+3*i+2] = vdbasis[0*vnumnodes+i];
		B_prime[(3*vnumnodes+pnumnodes)*5+3*i+0] = 0.;
		B_prime[(3*vnumnodes+pnumnodes)*5+3*i+1] = vdbasis[2*vnumnodes+i];
		B_prime[(3*vnumnodes+pnumnodes)*5+3*i+2] = vdbasis[1*vnumnodes+i];
		B_prime[(3*vnumnodes+pnumnodes)*6+3*i+0] = vdbasis[0*vnumnodes+i];
		B_prime[(3*vnumnodes+pnumnodes)*6+3*i+1] = vdbasis[1*vnumnodes+i];
		B_prime[(3*vnumnodes+pnumnodes)*6+3*i+2] = vdbasis[2*vnumnodes+i];
		B_prime[(3*vnumnodes+pnumnodes)*7+3*i+0] = 0.;
		B_prime[(3*vnumnodes+pnumnodes)*7+3*i+1] = 0.;
		B_prime[(3*vnumnodes+pnumnodes)*7+3*i+2] = 0.;
	}
	for(int i=0;i<pnumnodes;i++){
		B_prime[(3*vnumnodes+pnumnodes)*0+(3*vnumnodes)+i] = 0.;
		B_prime[(3*vnumnodes+pnumnodes)*1+(3*vnumnodes)+i] = 0.;
		B_prime[(3*vnumnodes+pnumnodes)*2+(3*vnumnodes)+i] = 0.;
		B_prime[(3*vnumnodes+pnumnodes)*3+(3*vnumnodes)+i] = 0.;
		B_prime[(3*vnumnodes+pnumnodes)*4+(3*vnumnodes)+i] = 0.;
		B_prime[(3*vnumnodes+pnumnodes)*5+(3*vnumnodes)+i] = 0.;
		B_prime[(3*vnumnodes+pnumnodes)*6+(3*vnumnodes)+i] = 0.;
		B_prime[(3*vnumnodes+pnumnodes)*7+(3*vnumnodes)+i] = pbasis[i];
	}

	/*Clean up*/
	xDelete<IssmDouble>(vdbasis);
	xDelete<IssmDouble>(pbasis);
}
/*}}}*/
/*FUNCTION PentaRef::GetBprimeFSGLS {{{*/
void PentaRef::GetBprimeFSGLS(IssmDouble* B_prime, IssmDouble* xyz_list, GaussPenta* gauss){
	/*	Compute B'  matrix. B'=[B1' B2' B3' B4' B5' B6' Bb'] where Bi' is of size 3*NDOF2. 
	 *	For node i, Bi' can be expressed in the actual coordinate system
	 *	by: 
	 *				Bi'=[  dh/dx   0          0       0]
	 *					 [   0      dh/dy      0       0]
	 *					 [   0      0         dh/dz    0]
	 *					 [  dh/dy   dh/dx      0       0]
	 *					 [  dh/dz   0        dh/dx     0]
	 *					 [   0      dh/dz    dh/dy     0]
	 *					 [  dh/dx   dh/dy    dh/dz     0]
	 *					 [   0      0          0       h]
	 *	where h is the interpolation function for node i.
	 *
	 * 	Same thing for the bubble fonction except that there is no fourth column
	 */

	int i;
	IssmDouble dbasis[3][NUMNODESP1];
	IssmDouble basis[NUMNODESP1];

	/*Get dbasismini in actual coordinate system: */
	GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss);
	GetNodalFunctionsP1(basis, gauss);

	/*B_primeuild B_prime: */
	for(i=0;i<NUMNODESP1;i++){
		B_prime[(NDOF4*NUMNODESP1)*0+NDOF4*i+0] = dbasis[0][i];
		B_prime[(NDOF4*NUMNODESP1)*0+NDOF4*i+1] = 0.;
		B_prime[(NDOF4*NUMNODESP1)*0+NDOF4*i+2] = 0.;
		B_prime[(NDOF4*NUMNODESP1)*1+NDOF4*i+0] = 0.;
		B_prime[(NDOF4*NUMNODESP1)*1+NDOF4*i+1] = dbasis[1][i];
		B_prime[(NDOF4*NUMNODESP1)*1+NDOF4*i+2] = 0.;
		B_prime[(NDOF4*NUMNODESP1)*2+NDOF4*i+0] = 0.;
		B_prime[(NDOF4*NUMNODESP1)*2+NDOF4*i+1] = 0.;
		B_prime[(NDOF4*NUMNODESP1)*2+NDOF4*i+2] = dbasis[2][i];
		B_prime[(NDOF4*NUMNODESP1)*3+NDOF4*i+0] = dbasis[1][i];
		B_prime[(NDOF4*NUMNODESP1)*3+NDOF4*i+1] = dbasis[0][i];
		B_prime[(NDOF4*NUMNODESP1)*3+NDOF4*i+2] = 0.;
		B_prime[(NDOF4*NUMNODESP1)*4+NDOF4*i+0] = dbasis[2][i];
		B_prime[(NDOF4*NUMNODESP1)*4+NDOF4*i+1] = 0.;
		B_prime[(NDOF4*NUMNODESP1)*4+NDOF4*i+2] = dbasis[0][i];
		B_prime[(NDOF4*NUMNODESP1)*5+NDOF4*i+0] = 0.;
		B_prime[(NDOF4*NUMNODESP1)*5+NDOF4*i+1] = dbasis[2][i];
		B_prime[(NDOF4*NUMNODESP1)*5+NDOF4*i+2] = dbasis[1][i];
		B_prime[(NDOF4*NUMNODESP1)*6+NDOF4*i+0] = dbasis[0][i];
		B_prime[(NDOF4*NUMNODESP1)*6+NDOF4*i+1] = dbasis[1][i];
		B_prime[(NDOF4*NUMNODESP1)*6+NDOF4*i+2] = dbasis[2][i];
		B_prime[(NDOF4*NUMNODESP1)*7+NDOF4*i+0] = 0.;
		B_prime[(NDOF4*NUMNODESP1)*7+NDOF4*i+1] = 0.;
		B_prime[(NDOF4*NUMNODESP1)*7+NDOF4*i+2] = 0.;
	}

	for(i=0;i<NUMNODESP1;i++){ //last column 
		B_prime[(NDOF4*NUMNODESP1)*0+NDOF4*i+3] = 0.;
		B_prime[(NDOF4*NUMNODESP1)*1+NDOF4*i+3] = 0.;
		B_prime[(NDOF4*NUMNODESP1)*2+NDOF4*i+3] = 0.;
		B_prime[(NDOF4*NUMNODESP1)*3+NDOF4*i+3] = 0.;
		B_prime[(NDOF4*NUMNODESP1)*4+NDOF4*i+3] = 0.;
		B_prime[(NDOF4*NUMNODESP1)*5+NDOF4*i+3] = 0.;
		B_prime[(NDOF4*NUMNODESP1)*6+NDOF4*i+3] = 0.;
		B_prime[(NDOF4*NUMNODESP1)*7+NDOF4*i+3] = - basis[i];
	}

}
/*}}}*/
/*FUNCTION PentaRef::GetBAdvec{{{*/
void PentaRef::GetBAdvec(IssmDouble* B_advec, IssmDouble* xyz_list, GaussPenta* gauss){
	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. 
	 * For node i, Bi' can be expressed in the actual coordinate system
	 * by: 
	 *       Bi_advec =[ h ]
	 *                 [ h ]
	 *                 [ h ]
	 * where h is the interpolation function for node i.
	 *
	 * We assume B has been allocated already, of size: 3x(NDOF1*NUMNODESP1)
	 */

	/*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 B: */
	for(int i=0;i<numnodes;i++){
		B_advec[numnodes*0+i] = basis[i];
		B_advec[numnodes*1+i] = basis[i];
		B_advec[numnodes*2+i] = basis[i];
	}

	/*Clean-up*/
	xDelete<IssmDouble>(basis);
}
/*}}}*/
/*FUNCTION PentaRef::GetBConduct{{{*/
void PentaRef::GetBConduct(IssmDouble* B_conduct, IssmDouble* xyz_list, GaussPenta* gauss){
	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. 
	 * For node i, Bi' can be expressed in the actual coordinate system
	 * by: 
	 *       Bi_conduct=[ dh/dx ]
	 *                  [ dh/dy ]
	 *                  [ dh/dz ]
	 * where h is the interpolation function for node i.
	 *
	 * We assume B has been allocated already, of size: 3x(NDOF1*numnodes)
	 */

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

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

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

	/*Clean-up*/
	xDelete<IssmDouble>(dbasis);
}
/*}}}*/
/*FUNCTION PentaRef::GetBVert{{{*/
void PentaRef::GetBVert(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){
	/*	Compute B  matrix. B=[dh1/dz dh2/dz dh3/dz dh4/dz dh5/dz dh6/dz];
		where hi is the interpolation function for node i.*/

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

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

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

	/*Clean-up*/
	xDelete<IssmDouble>(dbasis);
}
/*}}}*/
/*FUNCTION PentaRef::GetBprimeAdvec{{{*/
void PentaRef::GetBprimeAdvec(IssmDouble* Bprime_advec, IssmDouble* xyz_list, GaussPenta* gauss){
	/*Compute B  matrix. B=[B1 B2 B3 B4 B5 B6] where Bi is of size 5*NDOF1. 
	 * For node i, Bi' can be expressed in the actual coordinate system
	 * by: 
	 *       Biprime_advec=[ dh/dx ]
	 *                     [ dh/dy ]
	 *                     [ dh/dz ]
	 * where h is the interpolation function for node i.
	 *
	 * We assume B has been allocated already, of size: 3x(NDOF1*numnodes)
	 */

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

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

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

	/*Clean-up*/
	xDelete<IssmDouble>(dbasis);
}
/*}}}*/
/*FUNCTION PentaRef::GetBprimeVert{{{*/
void PentaRef::GetBprimeVert(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){

	GetNodalFunctions(B,gauss);

}
/*}}}*/
/*FUNCTION PentaRef::GetBHOFriction{{{*/
void PentaRef::GetBHOFriction(IssmDouble* B, GaussPenta* gauss){
	/*Compute B  matrix. B=[B1 B2 B3] where Bi is square and of size 2x2. 
	 ** 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 (2 x 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);

	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 PentaRef::GetLFS{{{*/
void PentaRef::GetLFS(IssmDouble* LFS, GaussPenta* gauss){
	/* Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof. 
	 * For node i, Li can be expressed in the actual coordinate system
	 * by: 
	 *       Li=[ h 0 0 0 ]
	 *	 	      [ 0 h 0 0 ]
	 * where h is the interpolation function for node i.
	 */

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

	/*Get nodal functions derivatives*/
	IssmDouble* vbasis=xNew<IssmDouble>(vnumnodes);
	GetNodalFunctionsVelocity(vbasis,gauss);

	/*Build LFS: */
	for(int i=0;i<vnumnodes;i++){
		LFS[(vnumdof+pnumdof)*0+3*i+0] = vbasis[i];
		LFS[(vnumdof+pnumdof)*0+3*i+1] = 0.;
		LFS[(vnumdof+pnumdof)*0+3*i+2] = 0.;

		LFS[(vnumdof+pnumdof)*1+3*i+0] = 0.;
		LFS[(vnumdof+pnumdof)*1+3*i+1] = vbasis[i];
		LFS[(vnumdof+pnumdof)*1+3*i+2] = 0.;
	}

	for(int i=0;i<pnumnodes;i++){
		LFS[(vnumdof+pnumdof)*0+i+vnumdof+0] = 0.;
		LFS[(vnumdof+pnumdof)*1+i+vnumdof+0] = 0.;
	}

	/*Clean-up*/
	xDelete<IssmDouble>(vbasis);
}
/*}}}*/
/*FUNCTION PentaRef::GetLprimeFS {{{*/
void PentaRef::GetLprimeFS(IssmDouble* LprimeFS, IssmDouble* xyz_list, GaussPenta* gauss){
	/* Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. 
	 * For node i, Lpi can be expressed in the actual coordinate system
	 * by: 
	 *       Lpi=[ h    0    0   0]1
	 *		       [ 0    h    0   0]2
	 *		       [ h    0    0   0]3
	 *		       [ 0    h    0   0]4
	 *		       [ 0    0    h   0]5
	 *		       [ 0    0    h   0]6
	 *		       [ 0    0  dh/dz 0]7
	 *		       [ 0    0  dh/dz 0]8
	 *		       [ 0    0  dh/dz 0]9
	 *		       [dh/dz 0  dh/dx 0]0
	 *		       [ 0 dh/dz dh/dy 0]1
	 *           [ 0    0    0   h]2
	 *           [ 0    0    0   h]3
	 *           [ 0    0    0   h]4
	 *
	 *       Li=[ h    0    0   0]1
	 *	 	      [ 0    h    0   0]2
	 *		      [ 0    0    h   0]3
	 *		      [ 0    0    h   0]4
	 *	 	      [ h    0    0   0]5
	 *	 	      [ 0    h    0   0]6
	 *	 	      [ h    0    0   0]7
	 *	 	      [ 0    h    0   0]8
	 *		      [ 0    0    h   0]9
	 *		      [ 0    0    h   0]0
	 *		      [ 0    0    h   0]1
	 *	 	      [ h    0    0   0]2
	 *	 	      [ 0    h    0   0]3
	 *		      [ 0    0    h   0]4
	 * where h is the interpolation function for node i.
	 */
	int i;
	int num_dof=4;

	IssmDouble L1L2l3[NUMNODESP1_2d];
	IssmDouble dbasis[3][NUMNODESP1];

	/*Get L1L2l3 in actual coordinate system: */
	L1L2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
	L1L2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
	L1L2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;

	GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);

	/*Build LprimeFS: */
	for(int i=0;i<3;i++){
		LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+0]  = L1L2l3[i];
		LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+1]  = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+2]  = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+3]  = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+0]  = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+1]  = L1L2l3[i];
		LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+2]  = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+3]  = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+0]  = L1L2l3[i];
		LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+1]  = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+2]  = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+3]  = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+0]  = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+1]  = L1L2l3[i];
		LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+2]  = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+3]  = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*4+num_dof*i+0]  = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*4+num_dof*i+1]  = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*4+num_dof*i+2]  = L1L2l3[i];
		LprimeFS[num_dof*NUMNODESP1_2d*4+num_dof*i+3]  = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*5+num_dof*i+0]  = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*5+num_dof*i+1]  = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*5+num_dof*i+2]  = L1L2l3[i];
		LprimeFS[num_dof*NUMNODESP1_2d*5+num_dof*i+3]  = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*6+num_dof*i+0]  = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*6+num_dof*i+1]  = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*6+num_dof*i+2]  = dbasis[2][i];
		LprimeFS[num_dof*NUMNODESP1_2d*6+num_dof*i+3]  = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*7+num_dof*i+0]  = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*7+num_dof*i+1]  = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*7+num_dof*i+2]  = dbasis[2][i];
		LprimeFS[num_dof*NUMNODESP1_2d*7+num_dof*i+3]  = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*8+num_dof*i+0]  = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*8+num_dof*i+1]  = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*8+num_dof*i+2]  = dbasis[2][i];
		LprimeFS[num_dof*NUMNODESP1_2d*8+num_dof*i+3]  = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*9+num_dof*i+0]  = dbasis[2][i];
		LprimeFS[num_dof*NUMNODESP1_2d*9+num_dof*i+1]  = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*9+num_dof*i+2]  = dbasis[0][i];
		LprimeFS[num_dof*NUMNODESP1_2d*9+num_dof*i+3]  = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*10+num_dof*i+0] = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*10+num_dof*i+1] = dbasis[2][i];
		LprimeFS[num_dof*NUMNODESP1_2d*10+num_dof*i+2] = dbasis[1][i];
		LprimeFS[num_dof*NUMNODESP1_2d*10+num_dof*i+3] = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*11+num_dof*i+0] = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*11+num_dof*i+1] = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*11+num_dof*i+2] = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*11+num_dof*i+3] = L1L2l3[i];
		LprimeFS[num_dof*NUMNODESP1_2d*12+num_dof*i+0] = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*12+num_dof*i+1] = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*12+num_dof*i+2] = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*12+num_dof*i+3] = L1L2l3[i];
		LprimeFS[num_dof*NUMNODESP1_2d*13+num_dof*i+0] = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*13+num_dof*i+1] = 0;
		LprimeFS[num_dof*NUMNODESP1_2d*13+num_dof*i+2] = 0;
		LprimeFS[num_dof*NUMNODESP1_2d*13+num_dof*i+3] = L1L2l3[i];
	}
}
/*}}}*/
/*FUNCTION PentaRef::GetLSSAFS {{{*/
void PentaRef::GetLSSAFS(IssmDouble* LFS, GaussPenta* gauss){
	/*
	 * Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof. 
	 * For node i, Li can be expressed in the actual coordinate system
	 * by: 
	 *       Li=[ h    0 ]
	 *	 	      [ 0    h ]
	 *	 	      [ h    0 ]
	 *	 	      [ 0    h ]
	 *	 	      [ h    0 ]
	 *	 	      [ 0    h ]
	 *	 	      [ h    0 ]
	 *	 	      [ 0    h ]
	 * where h is the interpolation function for node i.
	 */

	int num_dof=2;
	IssmDouble L1L2l3[NUMNODESP1_2d];

	/*Get L1L2l3 in actual coordinate system: */
	L1L2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
	L1L2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
	L1L2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;

	/*Build LFS: */
	for(int i=0;i<3;i++){
		LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = L1L2l3[i];
		LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0;
		LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0;
		LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = L1L2l3[i];
		LFS[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = L1L2l3[i];
		LFS[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0;
		LFS[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0;
		LFS[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = L1L2l3[i];
		LFS[num_dof*NUMNODESP1_2d*4+num_dof*i+0] = L1L2l3[i];
		LFS[num_dof*NUMNODESP1_2d*4+num_dof*i+1] = 0;
		LFS[num_dof*NUMNODESP1_2d*5+num_dof*i+0] = 0;
		LFS[num_dof*NUMNODESP1_2d*5+num_dof*i+1] = L1L2l3[i];
		LFS[num_dof*NUMNODESP1_2d*6+num_dof*i+0] = L1L2l3[i];
		LFS[num_dof*NUMNODESP1_2d*6+num_dof*i+1] = 0;
		LFS[num_dof*NUMNODESP1_2d*7+num_dof*i+0] = 0;
		LFS[num_dof*NUMNODESP1_2d*7+num_dof*i+1] = L1L2l3[i];
	}
}
/*}}}*/
/*FUNCTION PentaRef::GetLprimeSSAFS {{{*/
void PentaRef::GetLprimeSSAFS(IssmDouble* LprimeFS, IssmDouble* xyz_list, GaussPenta* gauss){
	/* Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. 
	 * For node i, Lpi can be expressed in the actual coordinate system
	 * by: 
	 *       Lpi=[ h    0    0   0]
	 *		       [ 0    h    0   0]
	 *		       [ 0    0    h   0]
	 *		       [ 0    0    h   0]
	 *		       [ 0    0  dh/dz 0]
	 *		       [ 0    0  dh/dz 0]
	 *           [ 0    0    0   h]
	 *           [ 0    0    0   h]
	 * where h is the interpolation function for node i.
	 */
	int num_dof=3;
	int num_dof_vel=3*NUMNODESP1b;
	int num_dof_total=3*NUMNODESP1b+1*NUMNODESP1;
	IssmDouble L1L2l3[NUMNODESP1_2d];
	IssmDouble dbasis[3][NUMNODESP1];

	/*Get L1L2l3 in actual coordinate system: */
	L1L2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
	L1L2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
	L1L2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;

	GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);

	/*Build LprimeFS: */
	for(int i=0;i<3;i++){
		LprimeFS[num_dof_total*0+num_dof*i+0] = L1L2l3[i];
		LprimeFS[num_dof_total*0+num_dof*i+1] = 0.;
		LprimeFS[num_dof_total*0+num_dof*i+2] = 0.;
		LprimeFS[num_dof_total*1+num_dof*i+0] = 0.;
		LprimeFS[num_dof_total*1+num_dof*i+1] = L1L2l3[i];
		LprimeFS[num_dof_total*1+num_dof*i+2] = 0.;
		LprimeFS[num_dof_total*2+num_dof*i+0] = 0.;
		LprimeFS[num_dof_total*2+num_dof*i+1] = 0.;
		LprimeFS[num_dof_total*2+num_dof*i+2] = L1L2l3[i];
		LprimeFS[num_dof_total*3+num_dof*i+0] = 0.;
		LprimeFS[num_dof_total*3+num_dof*i+1] = 0.;
		LprimeFS[num_dof_total*3+num_dof*i+2] = L1L2l3[i];
		LprimeFS[num_dof_total*4+num_dof*i+0] = 0.;
		LprimeFS[num_dof_total*4+num_dof*i+1] = 0.;
		LprimeFS[num_dof_total*4+num_dof*i+2] = dbasis[2][i];
		LprimeFS[num_dof_total*5+num_dof*i+0] = 0.;
		LprimeFS[num_dof_total*5+num_dof*i+1] = 0.;
		LprimeFS[num_dof_total*5+num_dof*i+2] = dbasis[2][i];
		LprimeFS[num_dof_total*6+num_dof*i+0] = 0.;
		LprimeFS[num_dof_total*6+num_dof*i+1] = 0.;
		LprimeFS[num_dof_total*6+num_dof*i+2] = 0.;
		LprimeFS[num_dof_total*7+num_dof*i+0] = 0.;
		LprimeFS[num_dof_total*7+num_dof*i+1] = 0.;
		LprimeFS[num_dof_total*7+num_dof*i+2] = 0.;
	}
	for(int i=3;i<7;i++){
		LprimeFS[num_dof_total*0+num_dof*i+0] = 0.;
		LprimeFS[num_dof_total*0+num_dof*i+1] = 0.;
		LprimeFS[num_dof_total*0+num_dof*i+2] = 0.;
		LprimeFS[num_dof_total*1+num_dof*i+0] = 0.;
		LprimeFS[num_dof_total*1+num_dof*i+1] = 0.;
		LprimeFS[num_dof_total*1+num_dof*i+2] = 0.;
		LprimeFS[num_dof_total*2+num_dof*i+0] = 0.;
		LprimeFS[num_dof_total*2+num_dof*i+1] = 0.;
		LprimeFS[num_dof_total*2+num_dof*i+2] = 0.;
		LprimeFS[num_dof_total*3+num_dof*i+0] = 0.;
		LprimeFS[num_dof_total*3+num_dof*i+1] = 0.;
		LprimeFS[num_dof_total*3+num_dof*i+2] = 0.;
		LprimeFS[num_dof_total*4+num_dof*i+0] = 0.;
		LprimeFS[num_dof_total*4+num_dof*i+1] = 0.;
		LprimeFS[num_dof_total*4+num_dof*i+2] = 0.;
		LprimeFS[num_dof_total*5+num_dof*i+0] = 0.;
		LprimeFS[num_dof_total*5+num_dof*i+1] = 0.;
		LprimeFS[num_dof_total*5+num_dof*i+2] = 0.;
		LprimeFS[num_dof_total*6+num_dof*i+0] = 0.;
		LprimeFS[num_dof_total*6+num_dof*i+1] = 0.;
		LprimeFS[num_dof_total*6+num_dof*i+2] = 0.;
		LprimeFS[num_dof_total*7+num_dof*i+0] = 0.;
		LprimeFS[num_dof_total*7+num_dof*i+1] = 0.;
		LprimeFS[num_dof_total*7+num_dof*i+2] = 0.;
	}
	for(int i=0;i<3;i++){
		LprimeFS[num_dof_total*0+num_dof_vel+i] = 0.;
		LprimeFS[num_dof_total*1+num_dof_vel+i] = 0.;
		LprimeFS[num_dof_total*2+num_dof_vel+i] = 0.;
		LprimeFS[num_dof_total*3+num_dof_vel+i] = 0.;
		LprimeFS[num_dof_total*4+num_dof_vel+i] = 0.;
		LprimeFS[num_dof_total*5+num_dof_vel+i] = 0.;
		LprimeFS[num_dof_total*6+num_dof_vel+i] = L1L2l3[i];
		LprimeFS[num_dof_total*7+num_dof_vel+i] = L1L2l3[i];
	}
	for(int i=3;i<6;i++){
		LprimeFS[num_dof_total*0+num_dof_vel+i] = 0.;
		LprimeFS[num_dof_total*1+num_dof_vel+i] = 0.;
		LprimeFS[num_dof_total*2+num_dof_vel+i] = 0.;
		LprimeFS[num_dof_total*3+num_dof_vel+i] = 0.;
		LprimeFS[num_dof_total*4+num_dof_vel+i] = 0.;
		LprimeFS[num_dof_total*5+num_dof_vel+i] = 0.;
		LprimeFS[num_dof_total*6+num_dof_vel+i] = 0.;
		LprimeFS[num_dof_total*7+num_dof_vel+i] = 0.;
	}
}
/*}}}*/
/*FUNCTION PentaRef::GetLFSSSA {{{*/
void PentaRef::GetLFSSSA(IssmDouble* LFS, GaussPenta* gauss){
	/* Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof. 
	 * For node i, Li can be expressed in the actual coordinate system
	 * by: 
	 *       Li=[ h    0    0 ]
	 *	 	      [ 0    h    0 ]
	 *		      [ 0    0    h ]
	 *		      [ 0    0    h ]
	 * where h is the interpolation function for node i.
	 */

	int num_dof=3;
	IssmDouble L1L2l3[NUMNODESP1_2d];

	/*Get L1L2l3 in actual coordinate system: */
	L1L2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
	L1L2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
	L1L2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;

	/*Build LFS: */
	for(int i=0;i<3;i++){
		LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = L1L2l3[i];
		LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;
		LFS[num_dof*NUMNODESP1_2d*0+num_dof*i+2] = 0.;
		LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;
		LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = L1L2l3[i];
		LFS[num_dof*NUMNODESP1_2d*1+num_dof*i+2] = 0.;
		LFS[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = 0.;
		LFS[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0.;
		LFS[num_dof*NUMNODESP1_2d*2+num_dof*i+2] = L1L2l3[i];
		LFS[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0.;
		LFS[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = 0.;
		LFS[num_dof*NUMNODESP1_2d*3+num_dof*i+2] = L1L2l3[i];
	}
}
/*}}}*/
/*FUNCTION PentaRef::GetLprimeFSSSA {{{*/
void PentaRef::GetLprimeFSSSA(IssmDouble* LprimeFS, IssmDouble* xyz_list, GaussPenta* gauss){
	/* Compute Lprime  matrix. Lprime=[Lp1 Lp2 Lp3] where Lpi is square and of size numdof. 
	 * For node i, Lpi can be expressed in the actual coordinate system
	 * by: 
	 *       Lpi=[ h    0 ]
	 *		       [ 0    h ]
	 *		       [ h    0 ]
	 *		       [ 0    h ]
	 * where h is the interpolation function for node i.
	 */
	int num_dof=2;
	IssmDouble L1L2l3[NUMNODESP1_2d];
	IssmDouble dbasis[3][NUMNODESP1];

	/*Get L1L2l3 in actual coordinate system: */
	L1L2l3[0]=gauss->coord1*(1-gauss->coord4)/2.0;
	L1L2l3[1]=gauss->coord2*(1-gauss->coord4)/2.0;
	L1L2l3[2]=gauss->coord3*(1-gauss->coord4)/2.0;
	GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);

	/*Build LprimeFS: */
	for(int i=0;i<3;i++){
		LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+0] = L1L2l3[i];
		LprimeFS[num_dof*NUMNODESP1_2d*0+num_dof*i+1] = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+0] = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*1+num_dof*i+1] = L1L2l3[i];
		LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+0] = L1L2l3[i];
		LprimeFS[num_dof*NUMNODESP1_2d*2+num_dof*i+1] = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+0] = 0.;
		LprimeFS[num_dof*NUMNODESP1_2d*3+num_dof*i+1] = L1L2l3[i];
	}
}
/*}}}*/
/*FUNCTION PentaRef::GetJacobian {{{*/
void PentaRef::GetJacobian(IssmDouble* J, IssmDouble* xyz_list,GaussPenta* gauss){
	/*The Jacobian is constant over the element, discard the gaussian points. 
	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/

	IssmDouble A1,A2,A3;  // area coordinates
	IssmDouble xi,eta,zi; // parametric coordinates
	IssmDouble x1,x2,x3,x4,x5,x6;
	IssmDouble y1,y2,y3,y4,y5,y6;
	IssmDouble z1,z2,z3,z4,z5,z6;

	/*Figure out xi,eta and zi (parametric coordinates), for this gaussian point: */
	A1  = gauss->coord1;
	A2  = gauss->coord2;
	A3  = gauss->coord3;
	xi  = A2-A1;
	eta = SQRT3*A3;
	zi  = gauss->coord4;

	x1=xyz_list[3*0+0];
	x2=xyz_list[3*1+0];
	x3=xyz_list[3*2+0];
	x4=xyz_list[3*3+0];
	x5=xyz_list[3*4+0];
	x6=xyz_list[3*5+0];

	y1=xyz_list[3*0+1];
	y2=xyz_list[3*1+1];
	y3=xyz_list[3*2+1];
	y4=xyz_list[3*3+1];
	y5=xyz_list[3*4+1];
	y6=xyz_list[3*5+1];

	z1=xyz_list[3*0+2];
	z2=xyz_list[3*1+2];
	z3=xyz_list[3*2+2];
	z4=xyz_list[3*3+2];
	z5=xyz_list[3*4+2];
	z6=xyz_list[3*5+2];

	J[NDOF3*0+0] = 0.25*(x1-x2-x4+x5)*zi+0.25*(-x1+x2-x4+x5);
	J[NDOF3*1+0] = SQRT3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*zi+SQRT3/12.0*(-x1-x2+2*x3-x4-x5+2*x6);
	J[NDOF3*2+0] = SQRT3/12.0*(x1+x2-2*x3-x4-x5+2*x6)*eta+1/4*(x1-x2-x4+x5)*xi +0.25*(-x1+x5-x2+x4);

	J[NDOF3*0+1] = 0.25*(y1-y2-y4+y5)*zi+0.25*(-y1+y2-y4+y5);
	J[NDOF3*1+1] = SQRT3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*zi+SQRT3/12.0*(-y1-y2+2*y3-y4-y5+2*y6);
	J[NDOF3*2+1] = SQRT3/12.0*(y1+y2-2*y3-y4-y5+2*y6)*eta+0.25*(y1-y2-y4+y5)*xi+0.25*(y4-y1+y5-y2);

	J[NDOF3*0+2] = 0.25*(z1-z2-z4+z5)*zi+0.25*(-z1+z2-z4+z5);
	J[NDOF3*1+2] = SQRT3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*zi+SQRT3/12.0*(-z1-z2+2*z3-z4-z5+2*z6);
	J[NDOF3*2+2] = SQRT3/12.0*(z1+z2-2*z3-z4-z5+2*z6)*eta+0.25*(z1-z2-z4+z5)*xi+0.25*(-z1+z5-z2+z4);
}
/*}}}*/
/*FUNCTION PentaRef::GetJacobianDeterminant {{{*/
void PentaRef::GetJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussPenta* gauss){
	/*On a penta, Jacobian varies according to coordinates. We need to get the Jacobian, and take 
	 * the determinant of it: */
	IssmDouble J[3][3];

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

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

}
/*}}}*/
/*FUNCTION PentaRef::GetTriaJacobianDeterminant{{{*/
void PentaRef::GetTriaJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussPenta* gauss){
	/*The Jacobian determinant 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 z1=xyz_list[3*0+2];
	IssmDouble x2=xyz_list[3*1+0];
	IssmDouble y2=xyz_list[3*1+1];
	IssmDouble z2=xyz_list[3*1+2];
	IssmDouble x3=xyz_list[3*2+0];
	IssmDouble y3=xyz_list[3*2+1];
	IssmDouble z3=xyz_list[3*2+2];

	/*Jdet = norm( AB ^ AC ) / (2 * area of the reference triangle), with areaRef=sqrt(3) */
	*Jdet=SQRT3/6.*pow(pow(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)),2)+pow(((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)),2)+pow(((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)),2),0.5);
	if(*Jdet<0) _error_("negative jacobian determinant!");
}
/*}}}*
/*FUNCTION PentaRef::GetSegmentJacobianDeterminant{{{*/
void PentaRef::GetSegmentJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,GaussPenta* gauss){
	/*The Jacobian determinant 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 z1=xyz_list[3*0+2];
	IssmDouble x2=xyz_list[3*1+0];
	IssmDouble y2=xyz_list[3*1+1];
	IssmDouble z2=xyz_list[3*1+2];

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

}
/*}}}*/
/*FUNCTION PentaRef::GetJacobianInvert {{{*/
void PentaRef::GetJacobianInvert(IssmDouble* Jinv, IssmDouble* xyz_list,GaussPenta* gauss){

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

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

	/*Invert Jacobian matrix: */
	Matrix3x3Invert(Jinv,&J[0][0]);
}
/*}}}*/
/*FUNCTION PentaRef::GetNodalFunctions{{{*/
void PentaRef::GetNodalFunctions(IssmDouble* basis,GaussPenta* gauss){
	/*This routine returns the values of the nodal functions  at the gaussian point.*/

	_assert_(basis);

	/*Get current coordinates in reference element*/
	IssmDouble zeta=gauss->coord4;

	switch(this->element_type){
		case P1Enum: case P1DGEnum:
			basis[0]=gauss->coord1*(1.-zeta)/2.;
			basis[1]=gauss->coord2*(1.-zeta)/2.;
			basis[2]=gauss->coord3*(1.-zeta)/2.;
			basis[3]=gauss->coord1*(1.+zeta)/2.;
			basis[4]=gauss->coord2*(1.+zeta)/2.;
			basis[5]=gauss->coord3*(1.+zeta)/2.;
			return;
		case P1bubbleEnum: case P1bubblecondensedEnum:
			basis[0]=gauss->coord1*(1.-zeta)/2.;
			basis[1]=gauss->coord2*(1.-zeta)/2.;
			basis[2]=gauss->coord3*(1.-zeta)/2.;
			basis[3]=gauss->coord1*(1.+zeta)/2.;
			basis[4]=gauss->coord2*(1.+zeta)/2.;
			basis[5]=gauss->coord3*(1.+zeta)/2.;
			basis[6]=27.*gauss->coord1*gauss->coord2*gauss->coord3*(1.+zeta)*(1.-zeta);
			return;
		case P2xP1Enum:
			/*Corner nodes*/
			basis[ 0]=gauss->coord1*(2.*gauss->coord1-1.)*(1.-zeta)/2.;
			basis[ 1]=gauss->coord2*(2.*gauss->coord2-1.)*(1.-zeta)/2.;
			basis[ 2]=gauss->coord3*(2.*gauss->coord3-1.)*(1.-zeta)/2.;
			basis[ 3]=gauss->coord1*(2.*gauss->coord1-1.)*(1.+zeta)/2.;
			basis[ 4]=gauss->coord2*(2.*gauss->coord2-1.)*(1.+zeta)/2.;
			basis[ 5]=gauss->coord3*(2.*gauss->coord3-1.)*(1.+zeta)/2.;
			/*mid-sides of triangles*/
			basis[ 6]=4.*gauss->coord3*gauss->coord2*(1.-zeta)/2.;
			basis[ 7]=4.*gauss->coord3*gauss->coord1*(1.-zeta)/2.;
			basis[ 8]=4.*gauss->coord1*gauss->coord2*(1.-zeta)/2.;
			basis[ 9]=4.*gauss->coord3*gauss->coord2*(1.+zeta)/2.;
			basis[10]=4.*gauss->coord3*gauss->coord1*(1.+zeta)/2.;
			basis[11]=4.*gauss->coord1*gauss->coord2*(1.+zeta)/2.;
			return;
		case P1xP2Enum:
			/*Corner nodes*/
			basis[ 0]=gauss->coord1*zeta*(zeta-1.)/2.;
			basis[ 1]=gauss->coord2*zeta*(zeta-1.)/2.;
			basis[ 2]=gauss->coord3*zeta*(zeta-1.)/2.;
			basis[ 3]=gauss->coord1*zeta*(zeta+1.)/2.;
			basis[ 4]=gauss->coord2*zeta*(zeta+1.)/2.;
			basis[ 5]=gauss->coord3*zeta*(zeta+1.)/2.;
			/*mid-sides of quads*/
			basis[ 6]=gauss->coord1*(1.-zeta*zeta);
			basis[ 7]=gauss->coord2*(1.-zeta*zeta);
			basis[ 8]=gauss->coord3*(1.-zeta*zeta);
			return;
		case P2Enum:
			/*Corner nodes*/
			basis[ 0]=gauss->coord1*(2.*gauss->coord1-1.)*zeta*(zeta-1.)/2.;
			basis[ 1]=gauss->coord2*(2.*gauss->coord2-1.)*zeta*(zeta-1.)/2.;
			basis[ 2]=gauss->coord3*(2.*gauss->coord3-1.)*zeta*(zeta-1.)/2.;
			basis[ 3]=gauss->coord1*(2.*gauss->coord1-1.)*zeta*(zeta+1.)/2.;
			basis[ 4]=gauss->coord2*(2.*gauss->coord2-1.)*zeta*(zeta+1.)/2.;
			basis[ 5]=gauss->coord3*(2.*gauss->coord3-1.)*zeta*(zeta+1.)/2.;
			/*mid-sides of quads*/
			basis[ 6]=gauss->coord1*(1.-zeta*zeta);
			basis[ 7]=gauss->coord2*(1.-zeta*zeta);
			basis[ 8]=gauss->coord3*(1.-zeta*zeta);
			/*mid-sides of triangles*/
			basis[ 9]=2.*gauss->coord3*gauss->coord2*zeta*(zeta-1.);
			basis[10]=2.*gauss->coord3*gauss->coord1*zeta*(zeta-1.);
			basis[11]=2.*gauss->coord1*gauss->coord2*zeta*(zeta-1.);
			basis[12]=2.*gauss->coord3*gauss->coord2*zeta*(zeta+1.);
			basis[13]=2.*gauss->coord3*gauss->coord1*zeta*(zeta+1.);
			basis[14]=2.*gauss->coord1*gauss->coord2*zeta*(zeta+1.);
			return;
		default:
			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
	}
}
/*}}}*/
/*FUNCTION PentaRef::GetNodalFunctionsVelocity{{{*/
void PentaRef::GetNodalFunctionsVelocity(IssmDouble* basis,GaussPenta* gauss){
	/*This routine returns the values of the nodal functions  at the gaussian point.*/

	switch(this->element_type){
		case P1P1Enum:
			this->element_type = P1Enum;
			this->GetNodalFunctions(basis,gauss);
			this->element_type = P1P1Enum;
			return;
		case P1P1GLSEnum:
			this->element_type = P1Enum;
			this->GetNodalFunctions(basis,gauss);
			this->element_type = P1P1GLSEnum;
			return;
		case MINIcondensedEnum:
			this->element_type = P1bubbleEnum;
			this->GetNodalFunctions(basis,gauss);
			this->element_type = MINIcondensedEnum;
			return;
		case MINIEnum:
			this->element_type = P1bubbleEnum;
			this->GetNodalFunctions(basis,gauss);
			this->element_type = MINIEnum;
			return;
		case TaylorHoodEnum:
			this->element_type = P2Enum;
			this->GetNodalFunctions(basis,gauss);
			this->element_type = TaylorHoodEnum;
			return;
		default:
			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
	}
}
/*}}}*/
/*FUNCTION PentaRef::GetNodalFunctionsPressure{{{*/
void PentaRef::GetNodalFunctionsPressure(IssmDouble* basis,GaussPenta* gauss){
	/*This routine returns the values of the nodal functions  at the gaussian point.*/

	switch(this->element_type){
		case P1P1Enum:
			this->element_type = P1Enum;
			this->GetNodalFunctions(basis,gauss);
			this->element_type = P1P1Enum;
			return;
		case P1P1GLSEnum:
			this->element_type = P1Enum;
			this->GetNodalFunctions(basis,gauss);
			this->element_type = P1P1GLSEnum;
			return;
		case MINIcondensedEnum:
			this->element_type = P1Enum;
			this->GetNodalFunctions(basis,gauss);
			this->element_type = MINIcondensedEnum;
			return;
		case MINIEnum:
			this->element_type = P1Enum;
			this->GetNodalFunctions(basis,gauss);
			this->element_type = MINIEnum;
			return;
		case TaylorHoodEnum:
			this->element_type = P1Enum;
			this->GetNodalFunctions(basis,gauss);
			this->element_type = TaylorHoodEnum;
			return;
		default:
			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
	}
}
/*}}}*/
/*FUNCTION PentaRef::GetNodalFunctionsDerivatives{{{*/
void PentaRef::GetNodalFunctionsDerivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussPenta* gauss){

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

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

	/*Get nodal functions derivatives in reference triangle*/
	IssmDouble* dbasis_ref=xNew<IssmDouble>(3*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]
	 * [dhi/dz]        [dhi/dzeta]
	 */

	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]+Jinv[0][2]*dbasis_ref[2*numnodes+i];
		dbasis[numnodes*1+i]=Jinv[1][0]*dbasis_ref[0*numnodes+i]+Jinv[1][1]*dbasis_ref[1*numnodes+i]+Jinv[1][2]*dbasis_ref[2*numnodes+i];
		dbasis[numnodes*2+i]=Jinv[2][0]*dbasis_ref[0*numnodes+i]+Jinv[2][1]*dbasis_ref[1*numnodes+i]+Jinv[2][2]*dbasis_ref[2*numnodes+i];
	}

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

}
/*}}}*/
/*FUNCTION PentaRef::GetNodalFunctionsDerivativesVelocity{{{*/
void PentaRef::GetNodalFunctionsDerivativesVelocity(IssmDouble* dbasis,IssmDouble* xyz_list, GaussPenta* gauss){
	switch(this->element_type){
		case P1P1Enum:
			this->element_type = P1Enum;
			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
			this->element_type = P1P1Enum;
			return;
		case P1P1GLSEnum:
			this->element_type = P1Enum;
			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
			this->element_type = P1P1GLSEnum;
			return;
		case MINIcondensedEnum:
			this->element_type = P1bubbleEnum;
			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
			this->element_type = MINIcondensedEnum;
			return;
		case MINIEnum:
			this->element_type = P1bubbleEnum;
			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
			this->element_type = MINIEnum;
			return;
		case TaylorHoodEnum:
			this->element_type = P2Enum;
			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
			this->element_type = TaylorHoodEnum;
			return;
		default:
			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
	}
}
/*}}}*/
/*FUNCTION PentaRef::GetNodalFunctionsDerivativesPressure{{{*/
void PentaRef::GetNodalFunctionsDerivativesPressure(IssmDouble* dbasis,IssmDouble* xyz_list, GaussPenta* gauss){
	switch(this->element_type){
		case P1P1Enum:
			this->element_type = P1Enum;
			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
			this->element_type = P1P1Enum;
			return;
		case P1P1GLSEnum:
			this->element_type = P1Enum;
			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
			this->element_type = P1P1GLSEnum;
			return;
		case MINIcondensedEnum:
			this->element_type = P1Enum;
			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
			this->element_type = MINIcondensedEnum;
			return;
		case MINIEnum:
			this->element_type = P1Enum;
			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
			this->element_type = MINIEnum;
			return;
		case TaylorHoodEnum:
			this->element_type = P1Enum;
			this->GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
			this->element_type = TaylorHoodEnum;
			return;
		default:
			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
	}
}
/*}}}*/
/*FUNCTION PentaRef::GetNodalFunctionsDerivativesReference{{{*/
void PentaRef::GetNodalFunctionsDerivativesReference(IssmDouble* dbasis,GaussPenta* 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);

	/*Get current coordinates in reference element*/
	IssmDouble zeta=gauss->coord4;

	switch(this->element_type){
		case P1Enum: case P1DGEnum:
			/*Nodal function 1*/
			dbasis[NUMNODESP1*0+0]   = (zeta-1.)/4.;
			dbasis[NUMNODESP1*1+0]   = SQRT3/12.*(zeta-1.);
			dbasis[NUMNODESP1*2+0]   = -.5*gauss->coord1;
			/*Nodal function 2*/
			dbasis[NUMNODESP1*0+1]   = (1.-zeta)/4.;
			dbasis[NUMNODESP1*1+1]   = SQRT3/12.*(zeta-1);
			dbasis[NUMNODESP1*2+1]   = -.5*gauss->coord2;
			/*Nodal function 3*/
			dbasis[NUMNODESP1*0+2]   = 0.;
			dbasis[NUMNODESP1*1+2]   = SQRT3/6.*(1.-zeta);
			dbasis[NUMNODESP1*2+2]   = -.5*gauss->coord3;
			/*Nodal function 4*/
			dbasis[NUMNODESP1*0+3]   = -(1.+zeta)/4.;
			dbasis[NUMNODESP1*1+3]   = -SQRT3/12.*(1.+zeta);
			dbasis[NUMNODESP1*2+3]   = .5*gauss->coord1;
			/*Nodal function 5*/
			dbasis[NUMNODESP1*0+4]   = (1.+zeta)/4.;
			dbasis[NUMNODESP1*1+4]   = -SQRT3/12.*(1.+zeta);
			dbasis[NUMNODESP1*2+4]   = .5*gauss->coord2;
			/*Nodal function 6*/
			dbasis[NUMNODESP1*0+5]   = 0.;
			dbasis[NUMNODESP1*1+5]   = SQRT3/6.*(1.+zeta);
			dbasis[NUMNODESP1*2+5]   = .5*gauss->coord3;
			return;
		case P1bubbleEnum: case P1bubblecondensedEnum:
			/*Nodal function 1*/
			dbasis[NUMNODESP1b*0+0]   = (zeta-1.)/4.;
			dbasis[NUMNODESP1b*1+0]   = SQRT3/12.*(zeta-1.);
			dbasis[NUMNODESP1b*2+0]   = -.5*gauss->coord1;
			/*Nodal function 2*/
			dbasis[NUMNODESP1b*0+1]   = (1.-zeta)/4.;
			dbasis[NUMNODESP1b*1+1]   = SQRT3/12.*(zeta-1);
			dbasis[NUMNODESP1b*2+1]   = -.5*gauss->coord2;
			/*Nodal function 3*/
			dbasis[NUMNODESP1b*0+2]   = 0.;
			dbasis[NUMNODESP1b*1+2]   = SQRT3/6.*(1.-zeta);
			dbasis[NUMNODESP1b*2+2]   = -.5*gauss->coord3;
			/*Nodal function 4*/
			dbasis[NUMNODESP1b*0+3]   = -(1.+zeta)/4.;
			dbasis[NUMNODESP1b*1+3]   = -SQRT3/12.*(1.+zeta);
			dbasis[NUMNODESP1b*2+3]   = .5*gauss->coord1;
			/*Nodal function 5*/
			dbasis[NUMNODESP1b*0+4]   = (1.+zeta)/4.;
			dbasis[NUMNODESP1b*1+4]   = -SQRT3/12.*(1.+zeta);
			dbasis[NUMNODESP1b*2+4]   = .5*gauss->coord2;
			/*Nodal function 6*/
			dbasis[NUMNODESP1b*0+5]   = 0.;
			dbasis[NUMNODESP1b*1+5]   = SQRT3/6.*(1.+zeta);
			dbasis[NUMNODESP1b*2+5]   = .5*gauss->coord3;
			/*Nodal function 7*/
			dbasis[NUMNODESP1b*0+6] = 27.*(1.+zeta)*(1.-zeta)*(-.5*gauss->coord2*gauss->coord3 + .5*gauss->coord1*gauss->coord3);
			dbasis[NUMNODESP1b*1+6] = 27.*(1.+zeta)*(1.-zeta)*SQRT3*(-1./6.*gauss->coord2*gauss->coord3 - 1./6.*gauss->coord1*gauss->coord3 +1./3.*gauss->coord1*gauss->coord2);
			dbasis[NUMNODESP1b*2+6] = -54*gauss->coord1*gauss->coord2*gauss->coord3*zeta;
			return;
		case P2xP1Enum:
			/*Nodal function 1*/
			dbasis[NUMNODESP2xP1*0+0 ] = .5*(1.-zeta)*(-2.*gauss->coord1 + 0.5);
			dbasis[NUMNODESP2xP1*1+0 ] = .5*(1.-zeta)*(-2.*SQRT3/3.*gauss->coord1 + SQRT3/6.);
			dbasis[NUMNODESP2xP1*2+0 ] = -.5*gauss->coord1*(2.*gauss->coord1-1.);
			/*Nodal function 2*/
			dbasis[NUMNODESP2xP1*0+1 ] = .5*(1.-zeta)*(+2.*gauss->coord2 - 0.5);
			dbasis[NUMNODESP2xP1*1+1 ] = .5*(1.-zeta)*(-2.*SQRT3/3.*gauss->coord2 + SQRT3/6.);
			dbasis[NUMNODESP2xP1*2+1 ] = -.5*gauss->coord2*(2.*gauss->coord2-1.);
			/*Nodal function 3*/
			dbasis[NUMNODESP2xP1*0+2 ] = 0.;
			dbasis[NUMNODESP2xP1*1+2 ] = .5*(1.-zeta)*(4.*SQRT3/3.*gauss->coord3 - SQRT3/3.);
			dbasis[NUMNODESP2xP1*2+2 ] = -.5*gauss->coord3*(2.*gauss->coord3-1.);
			/*Nodal function 4*/
			dbasis[NUMNODESP2xP1*0+3 ] = .5*(1.+zeta)*(-2.*gauss->coord1 + 0.5);
			dbasis[NUMNODESP2xP1*1+3 ] = .5*(1.+zeta)*(-2.*SQRT3/3.*gauss->coord1 + SQRT3/6.);
			dbasis[NUMNODESP2xP1*2+3 ] = .5*gauss->coord1*(2.*gauss->coord1-1.);
			/*Nodal function 5*/
			dbasis[NUMNODESP2xP1*0+4 ] = .5*(1.+zeta)*(+2.*gauss->coord2 - 0.5);
			dbasis[NUMNODESP2xP1*1+4 ] = .5*(1.+zeta)*(-2.*SQRT3/3.*gauss->coord2 + SQRT3/6.);
			dbasis[NUMNODESP2xP1*2+4 ] = .5*gauss->coord2*(2.*gauss->coord2-1.);
			/*Nodal function 6*/
			dbasis[NUMNODESP2xP1*0+5 ] = 0.;
			dbasis[NUMNODESP2xP1*1+5 ] = .5*(1.+zeta)*(4.*SQRT3/3.*gauss->coord3 - SQRT3/3.);
			dbasis[NUMNODESP2xP1*2+5 ] = .5*gauss->coord3*(2.*gauss->coord3-1.);

			/*Nodal function 7*/
			dbasis[NUMNODESP2xP1*0+6 ] = (1.-zeta)*gauss->coord3;
			dbasis[NUMNODESP2xP1*1+6 ] = .5*(1.-zeta)*(+4.*SQRT3/3.*gauss->coord2 - 2.*SQRT3/3.*gauss->coord3);
			dbasis[NUMNODESP2xP1*2+6 ] = -2.*gauss->coord3*gauss->coord2;
			/*Nodal function 8*/
			dbasis[NUMNODESP2xP1*0+7 ] = -(1.-zeta)*gauss->coord3;
			dbasis[NUMNODESP2xP1*1+7 ] = .5*(1.-zeta)*(+4.*SQRT3/3.*gauss->coord1 - 2.*SQRT3/3.*gauss->coord3);
			dbasis[NUMNODESP2xP1*2+7 ] = -2.*gauss->coord3*gauss->coord1;
			/*Nodal function 9*/
			dbasis[NUMNODESP2xP1*0+8 ] = (1.-zeta)*(gauss->coord1-gauss->coord2);
			dbasis[NUMNODESP2xP1*1+8 ] = .5*(1.-zeta)*(-2.*SQRT3/3.*(gauss->coord1+gauss->coord2));
			dbasis[NUMNODESP2xP1*2+8 ] = -2.*gauss->coord1*gauss->coord2;
			/*Nodal function 10*/
			dbasis[NUMNODESP2xP1*0+9 ] = (1.+zeta)*gauss->coord3;
			dbasis[NUMNODESP2xP1*1+9 ] = .5*(1.+zeta)*(+4.*SQRT3/3.*gauss->coord2 - 2.*SQRT3/3.*gauss->coord3);
			dbasis[NUMNODESP2xP1*2+9 ] = 2.*gauss->coord3*gauss->coord2;
			/*Nodal function 11*/
			dbasis[NUMNODESP2xP1*0+10] = -(1.+zeta)*gauss->coord3;
			dbasis[NUMNODESP2xP1*1+10] = .5*(1.+zeta)*(+4.*SQRT3/3.*gauss->coord1 - 2.*SQRT3/3.*gauss->coord3);
			dbasis[NUMNODESP2xP1*2+10] = 2.*gauss->coord3*gauss->coord1;
			/*Nodal function 12*/
			dbasis[NUMNODESP2xP1*0+11] = (1.+zeta)*(gauss->coord1-gauss->coord2);
			dbasis[NUMNODESP2xP1*1+11] = .5*(1.+zeta)*(-2.*SQRT3/3.*(gauss->coord1+gauss->coord2));
			dbasis[NUMNODESP2xP1*2+11] = 2.*gauss->coord1*gauss->coord2;
			return;
		case P1xP2Enum:
			/*Nodal function 1*/
			dbasis[NUMNODESP1xP2*0+0]   = -zeta*(zeta-1.)/4.;
			dbasis[NUMNODESP1xP2*1+0]   = -SQRT3/12.*zeta*(zeta-1.);
			dbasis[NUMNODESP1xP2*2+0]   = .5*(2.*zeta-1.)*gauss->coord1;
			/*Nodal function 2*/
			dbasis[NUMNODESP1xP2*0+1]   = zeta*(zeta-1.)/4.;
			dbasis[NUMNODESP1xP2*1+1]   = -SQRT3/12.*zeta*(zeta-1);
			dbasis[NUMNODESP1xP2*2+1]   = .5*(2.*zeta-1.)*gauss->coord2;
			/*Nodal function 3*/
			dbasis[NUMNODESP1xP2*0+2]   = 0.;
			dbasis[NUMNODESP1xP2*1+2]   = SQRT3/6.*zeta*(zeta-1.);
			dbasis[NUMNODESP1xP2*2+2]   = .5*(2.*zeta-1.)*gauss->coord3;
			/*Nodal function 4*/
			dbasis[NUMNODESP1xP2*0+3]   = -zeta*(zeta+1)/4.;
			dbasis[NUMNODESP1xP2*1+3]   = -SQRT3/12.*zeta*(zeta+1.);
			dbasis[NUMNODESP1xP2*2+3]   = .5*(2.*zeta+1.)*gauss->coord1;
			/*Nodal function 5*/
			dbasis[NUMNODESP1xP2*0+4]   = zeta*(zeta+1.)/4.;
			dbasis[NUMNODESP1xP2*1+4]   = -SQRT3/12.*zeta*(zeta+1.);
			dbasis[NUMNODESP1xP2*2+4]   = .5*(2.*zeta+1.)*gauss->coord2;
			/*Nodal function 6*/
			dbasis[NUMNODESP1xP2*0+5]   = 0.;
			dbasis[NUMNODESP1xP2*1+5]   = SQRT3/6.*zeta*(zeta+1.);
			dbasis[NUMNODESP1xP2*2+5]   = .5*(2.*zeta+1.)*gauss->coord3;

			/*Nodal function 7*/
			dbasis[NUMNODESP1xP2*0+6 ] = -0.5*(1.-zeta*zeta);
			dbasis[NUMNODESP1xP2*1+6 ] = -SQRT3/6.*(1.-zeta*zeta);
			dbasis[NUMNODESP1xP2*2+6 ] = -2.*zeta*gauss->coord1;
			/*Nodal function 8*/
			dbasis[NUMNODESP1xP2*0+7 ] = 0.5*(1.-zeta*zeta);
			dbasis[NUMNODESP1xP2*1+7 ] = -SQRT3/6.*(1.-zeta*zeta);
			dbasis[NUMNODESP1xP2*2+7 ] = -2.*zeta*gauss->coord2;
			/*Nodal function 9*/
			dbasis[NUMNODESP1xP2*0+8 ] = 0.;
			dbasis[NUMNODESP1xP2*1+8 ] = SQRT3/3.*(1.-zeta*zeta);
			dbasis[NUMNODESP1xP2*2+8 ] = -2.*zeta*gauss->coord3;
			return;
		case P2Enum:
			/*Nodal function 1*/
			dbasis[NUMNODESP2*0+0 ] = .5*zeta*(zeta-1.)*(-2.*gauss->coord1 + 0.5);
			dbasis[NUMNODESP2*1+0 ] = .5*zeta*(zeta-1.)*(-2.*SQRT3/3.*gauss->coord1 + SQRT3/6.);
			dbasis[NUMNODESP2*2+0 ] = .5*(2.*zeta-1.)*gauss->coord1*(2.*gauss->coord1-1.);
			/*Nodal function 2*/
			dbasis[NUMNODESP2*0+1 ] = .5*zeta*(zeta-1.)*(+2.*gauss->coord2 - 0.5);
			dbasis[NUMNODESP2*1+1 ] = .5*zeta*(zeta-1.)*(-2.*SQRT3/3.*gauss->coord2 + SQRT3/6.);
			dbasis[NUMNODESP2*2+1 ] = .5*(2.*zeta-1.)*gauss->coord2*(2.*gauss->coord2-1.);
			/*Nodal function 3*/
			dbasis[NUMNODESP2*0+2 ] = 0.;
			dbasis[NUMNODESP2*1+2 ] = .5*zeta*(zeta-1.)*(4.*SQRT3/3.*gauss->coord3 - SQRT3/3.);
			dbasis[NUMNODESP2*2+2 ] = .5*(2.*zeta-1.)*gauss->coord3*(2.*gauss->coord3-1.);
			/*Nodal function 4*/
			dbasis[NUMNODESP2*0+3 ] = .5*zeta*(zeta+1.)*(-2.*gauss->coord1 + 0.5);
			dbasis[NUMNODESP2*1+3 ] = .5*zeta*(zeta+1.)*(-2.*SQRT3/3.*gauss->coord1 + SQRT3/6.);
			dbasis[NUMNODESP2*2+3 ] = .5*(2.*zeta+1.)*gauss->coord1*(2.*gauss->coord1-1.);
			/*Nodal function 5*/
			dbasis[NUMNODESP2*0+4 ] = .5*zeta*(zeta+1.)*(+2.*gauss->coord2 - 0.5);
			dbasis[NUMNODESP2*1+4 ] = .5*zeta*(zeta+1.)*(-2.*SQRT3/3.*gauss->coord2 + SQRT3/6.);
			dbasis[NUMNODESP2*2+4 ] = .5*(2.*zeta+1.)*gauss->coord2*(2.*gauss->coord2-1.);
			/*Nodal function 6*/
			dbasis[NUMNODESP2*0+5 ] = 0.;
			dbasis[NUMNODESP2*1+5 ] = .5*zeta*(zeta+1.)*(4.*SQRT3/3.*gauss->coord3 - SQRT3/3.);
			dbasis[NUMNODESP2*2+5 ] = .5*(2.*zeta+1.)*gauss->coord3*(2.*gauss->coord3-1.);

			/*Nodal function 7*/
			dbasis[NUMNODESP2*0+6 ] = -0.5*(1.-zeta*zeta);
			dbasis[NUMNODESP2*1+6 ] = -SQRT3/6.*(1.-zeta*zeta);
			dbasis[NUMNODESP2*2+6 ] = -2.*zeta*gauss->coord1;
			/*Nodal function 8*/
			dbasis[NUMNODESP2*0+7 ] = 0.5*(1.-zeta*zeta);
			dbasis[NUMNODESP2*1+7 ] = -SQRT3/6.*(1.-zeta*zeta);
			dbasis[NUMNODESP2*2+7 ] = -2.*zeta*gauss->coord2;
			/*Nodal function 9*/
			dbasis[NUMNODESP2*0+8 ] = 0.;
			dbasis[NUMNODESP2*1+8 ] = SQRT3/3.*(1.-zeta*zeta);
			dbasis[NUMNODESP2*2+8 ] = -2.*zeta*gauss->coord3;

			/*Nodal function 10*/
			dbasis[NUMNODESP2*0+9 ] = zeta*(zeta-1.)*gauss->coord3;
			dbasis[NUMNODESP2*1+9 ] = .5*zeta*(zeta-1.)*(+4.*SQRT3/3.*gauss->coord2 - 2.*SQRT3/3.*gauss->coord3);
			dbasis[NUMNODESP2*2+9 ] = 2.*gauss->coord3*gauss->coord2*(2.*zeta-1.);
			/*Nodal function 11*/
			dbasis[NUMNODESP2*0+10] = -zeta*(zeta-1.)*gauss->coord3;
			dbasis[NUMNODESP2*1+10] = .5*zeta*(zeta-1.)*(+4.*SQRT3/3.*gauss->coord1 - 2.*SQRT3/3.*gauss->coord3);
			dbasis[NUMNODESP2*2+10] = 2.*gauss->coord3*gauss->coord1*(2.*zeta-1.);
			/*Nodal function 12*/
			dbasis[NUMNODESP2*0+11] = zeta*(zeta-1.)*(gauss->coord1-gauss->coord2);
			dbasis[NUMNODESP2*1+11] = .5*zeta*(zeta-1.)*(-2.*SQRT3/3.*(gauss->coord1+gauss->coord2));
			dbasis[NUMNODESP2*2+11] = 2.*gauss->coord1*gauss->coord2*(2.*zeta-1.);
			/*Nodal function 13*/
			dbasis[NUMNODESP2*0+12] = zeta*(zeta+1.)*gauss->coord3;
			dbasis[NUMNODESP2*1+12] = .5*zeta*(zeta+1.)*(+4.*SQRT3/3.*gauss->coord2 - 2.*SQRT3/3.*gauss->coord3);
			dbasis[NUMNODESP2*2+12] = 2.*gauss->coord3*gauss->coord2*(2.*zeta+1.);
			/*Nodal function 14*/
			dbasis[NUMNODESP2*0+13] = -zeta*(zeta+1.)*gauss->coord3;
			dbasis[NUMNODESP2*1+13] = .5*zeta*(zeta+1.)*(+4.*SQRT3/3.*gauss->coord1 - 2.*SQRT3/3.*gauss->coord3);
			dbasis[NUMNODESP2*2+13] = 2.*gauss->coord3*gauss->coord1*(2.*zeta+1.);
			/*Nodal function 15*/
			dbasis[NUMNODESP2*0+14] = zeta*(zeta+1.)*(gauss->coord1-gauss->coord2);
			dbasis[NUMNODESP2*1+14] = .5*zeta*(zeta+1.)*(-2.*SQRT3/3.*(gauss->coord1+gauss->coord2));
			dbasis[NUMNODESP2*2+14] = 2.*gauss->coord1*gauss->coord2*(2.*zeta+1.);
			return;
		default:
			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
	}

}
/*}}}*/
/*FUNCTION PentaRef::GetNodalFunctionsMINI{{{*/
void PentaRef::GetNodalFunctionsMINI(IssmDouble* l1l7, GaussPenta* gauss){
	/*This routine returns the values of the nodal functions  at the gaussian point.*/

	l1l7[0]=gauss->coord1*(1.0-gauss->coord4)/2.0;
	l1l7[1]=gauss->coord2*(1.0-gauss->coord4)/2.0;
	l1l7[2]=gauss->coord3*(1.0-gauss->coord4)/2.0;
	l1l7[3]=gauss->coord1*(1.0+gauss->coord4)/2.0;
	l1l7[4]=gauss->coord2*(1.0+gauss->coord4)/2.0;
	l1l7[5]=gauss->coord3*(1.0+gauss->coord4)/2.0;
	l1l7[6]=27*gauss->coord1*gauss->coord2*gauss->coord3*(1.0+gauss->coord4)*(1.0-gauss->coord4);

}
/*}}}*/
/*FUNCTION PentaRef::GetNodalFunctionsMINIDerivatives{{{*/
void PentaRef::GetNodalFunctionsMINIDerivatives(IssmDouble* dbasismini,IssmDouble* xyz_list, GaussPenta* gauss){

	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
	 * actual coordinate system): */

	IssmDouble    dbasismini_ref[3][NUMNODESP1b];
	IssmDouble    Jinv[3][3];

	/*Get derivative values with respect to parametric coordinate system: */
	GetNodalFunctionsMINIDerivativesReference(&dbasismini_ref[0][0], gauss); 

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

	/*Build dbasis: 
	 *
	 * [dhi/dx]= Jinv'*[dhi/dr]
	 * [dhi/dy]        [dhi/ds]
	 * [dhi/dz]        [dhi/dzeta]
	 */

	for(int i=0;i<NUMNODESP1b;i++){
		*(dbasismini+NUMNODESP1b*0+i)=Jinv[0][0]*dbasismini_ref[0][i]+Jinv[0][1]*dbasismini_ref[1][i]+Jinv[0][2]*dbasismini_ref[2][i];
		*(dbasismini+NUMNODESP1b*1+i)=Jinv[1][0]*dbasismini_ref[0][i]+Jinv[1][1]*dbasismini_ref[1][i]+Jinv[1][2]*dbasismini_ref[2][i];
		*(dbasismini+NUMNODESP1b*2+i)=Jinv[2][0]*dbasismini_ref[0][i]+Jinv[2][1]*dbasismini_ref[1][i]+Jinv[2][2]*dbasismini_ref[2][i];
	}

}
/*}}}*/
/*FUNCTION PentaRef::GetNodalFunctionsMINIDerivativesReference{{{*/
void PentaRef::GetNodalFunctionsMINIDerivativesReference(IssmDouble* dbasis,GaussPenta* gauss){

	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
	 * natural coordinate system) at the gaussian point. */
	IssmDouble zeta=gauss->coord4;

	/*Nodal function 1*/
	dbasis[NUMNODESP1b*0+0]=-0.5*(1.0-zeta)/2.0;
	dbasis[NUMNODESP1b*1+0]=-SQRT3/6.0*(1.0-zeta)/2.0;
	dbasis[NUMNODESP1b*2+0]=-0.5*gauss->coord1;
	/*Nodal function 2*/
	dbasis[NUMNODESP1b*0+1]=0.5*(1.0-zeta)/2.0;
	dbasis[NUMNODESP1b*1+1]=-SQRT3/6.0*(1.0-zeta)/2.0;
	dbasis[NUMNODESP1b*2+1]=-0.5*gauss->coord2;
	/*Nodal function 3*/
	dbasis[NUMNODESP1b*0+2]=0.;
	dbasis[NUMNODESP1b*1+2]=SQRT3/3.0*(1.0-zeta)/2.0;
	dbasis[NUMNODESP1b*2+2]=-0.5*gauss->coord3;
	/*Nodal function 4*/
	dbasis[NUMNODESP1b*0+3]=-0.5*(1.0+zeta)/2.0;
	dbasis[NUMNODESP1b*1+3]=-SQRT3/6.0*(1.0+zeta)/2.0;
	dbasis[NUMNODESP1b*2+3]=0.5*gauss->coord1;
	/*Nodal function 5*/
	dbasis[NUMNODESP1b*0+4]=0.5*(1.0+zeta)/2.0;
	dbasis[NUMNODESP1b*1+4]=-SQRT3/6.0*(1.0+zeta)/2.0;
	dbasis[NUMNODESP1b*2+4]=0.5*gauss->coord2;
	/*Nodal function 6*/
	dbasis[NUMNODESP1b*0+5]=0.;
	dbasis[NUMNODESP1b*1+5]=SQRT3/3.0*(1.0+zeta)/2.0;
	dbasis[NUMNODESP1b*2+5]=0.5*gauss->coord3;
	/*Nodal function 7*/
	dbasis[NUMNODESP1b*0+6]=27.*(1.+zeta)*(1.-zeta)*(-.5*gauss->coord2*gauss->coord3 + .5*gauss->coord1*gauss->coord3);
	dbasis[NUMNODESP1b*1+6]=27.*(1.+zeta)*(1.-zeta)*SQRT3*(-1./6.*gauss->coord2*gauss->coord3 - 1./6.*gauss->coord1*gauss->coord3 +1./3.*gauss->coord1*gauss->coord2);
	dbasis[NUMNODESP1b*2+6]=27*gauss->coord1*gauss->coord2*gauss->coord3*(-2.0*zeta);
}
/*}}}*/
/*FUNCTION PentaRef::GetNodalFunctionsP1 {{{*/
void PentaRef::GetNodalFunctionsP1(IssmDouble* basis, GaussPenta* gauss){
	/*This routine returns the values of the nodal functions  at the gaussian point.*/

	basis[0]=gauss->coord1*(1-gauss->coord4)/2.0;
	basis[1]=gauss->coord2*(1-gauss->coord4)/2.0;
	basis[2]=gauss->coord3*(1-gauss->coord4)/2.0;
	basis[3]=gauss->coord1*(1+gauss->coord4)/2.0;
	basis[4]=gauss->coord2*(1+gauss->coord4)/2.0;
	basis[5]=gauss->coord3*(1+gauss->coord4)/2.0;

}
/*}}}*/
/*FUNCTION PentaRef::GetNodalFunctionsP1Derivatives {{{*/
void PentaRef::GetNodalFunctionsP1Derivatives(IssmDouble* dbasis,IssmDouble* xyz_list, GaussPenta* gauss){

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

	/*Get derivative values with respect to parametric coordinate system: */
	GetNodalFunctionsP1DerivativesReference(&dbasis_ref[0][0], gauss); 

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

	/*Build dh1dh3: 
	 *
	 * [dhi/dx]= Jinv*[dhi/dr]
	 * [dhi/dy]       [dhi/ds]
	 * [dhi/dz]       [dhi/dn]
	 */

	for (int i=0;i<NUMNODESP1;i++){
		*(dbasis+NUMNODESP1*0+i)=Jinv[0][0]*dbasis_ref[0][i]+Jinv[0][1]*dbasis_ref[1][i]+Jinv[0][2]*dbasis_ref[2][i];
		*(dbasis+NUMNODESP1*1+i)=Jinv[1][0]*dbasis_ref[0][i]+Jinv[1][1]*dbasis_ref[1][i]+Jinv[1][2]*dbasis_ref[2][i];
		*(dbasis+NUMNODESP1*2+i)=Jinv[2][0]*dbasis_ref[0][i]+Jinv[2][1]*dbasis_ref[1][i]+Jinv[2][2]*dbasis_ref[2][i];
	}

}
/*}}}*/
/*FUNCTION PentaRef::GetNodalFunctionsP1DerivativesReference {{{*/
void PentaRef::GetNodalFunctionsP1DerivativesReference(IssmDouble* dbasis,GaussPenta* gauss){

	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
	 * natural coordinate system) at the gaussian point. Those values vary along xi,eta,z */

	IssmDouble zeta=gauss->coord4;

	/*Nodal function 1*/
	dbasis[NUMNODESP1*0+0]=-0.5*(1.0-zeta)/2.0;
	dbasis[NUMNODESP1*1+0]=-SQRT3/6.0*(1.0-zeta)/2.0;
	dbasis[NUMNODESP1*2+0]=-0.5*gauss->coord1;
	/*Nodal function 2*/
	dbasis[NUMNODESP1*0+1]=0.5*(1.0-zeta)/2.0;
	dbasis[NUMNODESP1*1+1]=-SQRT3/6.0*(1.0-zeta)/2.0;
	dbasis[NUMNODESP1*2+1]=-0.5*gauss->coord2;
	/*Nodal function 3*/
	dbasis[NUMNODESP1*0+2]=0.;
	dbasis[NUMNODESP1*1+2]=SQRT3/3.0*(1.0-zeta)/2.0;
	dbasis[NUMNODESP1*2+2]=-0.5*gauss->coord3;
	/*Nodal function 4*/
	dbasis[NUMNODESP1*0+3]=-0.5*(1.0+zeta)/2.0;
	dbasis[NUMNODESP1*1+3]=-SQRT3/6.0*(1.0+zeta)/2.0;
	dbasis[NUMNODESP1*2+3]=0.5*gauss->coord1;
	/*Nodal function 5*/
	dbasis[NUMNODESP1*0+4]=0.5*(1.0+zeta)/2.0;
	dbasis[NUMNODESP1*1+4]=-SQRT3/6.0*(1.0+zeta)/2.0;
	dbasis[NUMNODESP1*2+4]=0.5*gauss->coord2;
	/*Nodal function 6*/
	dbasis[NUMNODESP1*0+5]=0.;
	dbasis[NUMNODESP1*1+5]=SQRT3/3.0*(1.0+zeta)/2.0;
	dbasis[NUMNODESP1*2+5]=0.5*gauss->coord3;
}
/*}}}*/
/*FUNCTION PentaRef::GetQuadNodalFunctions {{{*/
void PentaRef::GetQuadNodalFunctions(IssmDouble* l1l4,GaussPenta* gauss,int index1,int index2,int index3,int index4){
	/*This routine returns the values of the nodal functions  at the gaussian point.*/

	IssmDouble BasisFunctions[6];

	GetNodalFunctionsP1(&BasisFunctions[0],gauss);

	_assert_(index1>=0 && index1<6);
	_assert_(index2>=0 && index2<6);
	_assert_(index3>=0 && index3<6);
	_assert_(index4>=0 && index4<6);

	l1l4[0]=BasisFunctions[index1];
	l1l4[1]=BasisFunctions[index2];
	l1l4[2]=BasisFunctions[index3];
	l1l4[3]=BasisFunctions[index4];

}
/*}}}*/
/*FUNCTION PentaRef::GetQuadJacobianDeterminant{{{*/
void PentaRef::GetQuadJacobianDeterminant(IssmDouble* Jdet,IssmDouble xyz_list[4][3],GaussPenta* gauss){
	/*This routine returns the values of the nodal functions  at the gaussian point.*/

	IssmDouble x1,x2,x3,x4,y1,y2,y3,y4,z1,z2,z3,z4;

	x1=xyz_list[0][0];
	y1=xyz_list[0][1];
	z1=xyz_list[0][2];
	x2=xyz_list[1][0];
	y2=xyz_list[1][1];
	z2=xyz_list[1][2];
	x3=xyz_list[2][0];
	y3=xyz_list[2][1];
	z3=xyz_list[2][2];
	x4=xyz_list[3][0];
	y4=xyz_list[3][1];
	z4=xyz_list[3][2];

	/*Jdet = (Area of the trapezoid)/(Area trapezoid ref) with AreaRef = 4*/
	/*Area of a trabezoid = altitude * (base1 + base2)/2 */
	*Jdet= sqrt((x2-x1)*(x2-x1) + (y2-y1)*(y2-y1)) * (z4-z1 + z3-z2)/8.;
	if(*Jdet<0.) _error_("negative jacobian determinant!");

}
/*}}}*/
/*FUNCTION PentaRef::GetInputValue{{{*/
void PentaRef::GetInputValue(IssmDouble* pvalue,IssmDouble* plist,GaussPenta* gauss){
	/*P1 interpolation on Gauss point*/

	/*intermediary*/
	IssmDouble basis[6];

	/*nodal functions: */
	GetNodalFunctionsP1(&basis[0],gauss);

	/*Assign output pointers:*/
	*pvalue=basis[0]*plist[0]+basis[1]*plist[1]+basis[2]*plist[2]+basis[3]*plist[3]+basis[4]*plist[4]+basis[5]*plist[5];

}
/*}}}*/
/*FUNCTION PentaRef::GetInputDerivativeValue{{{*/
void PentaRef::GetInputDerivativeValue(IssmDouble* p, IssmDouble* plist,IssmDouble* xyz_list, GaussPenta* gauss){
	/*From node values of parameter p (p_list[0], p_list[1], p_list[2],
	 * p_list[3], p_list[4] and p_list[4]), return parameter derivative value at
	 * gaussian point specified by gauss_coord:
	 *   dp/dx=p_list[0]*dh1/dx+p_list[1]*dh2/dx+p_list[2]*dh3/dx+p_list[3]*dh4/dx+p_list[4]*dh5/dx+p_list[5]*dh6/dx;
	 *   dp/dy=p_list[0]*dh1/dy+p_list[1]*dh2/dy+p_list[2]*dh3/dy+p_list[3]*dh4/dy+p_list[4]*dh5/dy+p_list[5]*dh6/dy;
	 *   dp/dz=p_list[0]*dh1/dz+p_list[1]*dh2/dz+p_list[2]*dh3/dz+p_list[3]*dh4/dz+p_list[4]*dh5/dz+p_list[5]*dh6/dz;
	 *
	 *   p is a vector of size 3x1 already allocated.
	 */
	IssmDouble dbasis[3][NUMNODESP1];

	/*Get nodal funnctions derivatives in actual coordinate system: */
	GetNodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list, gauss);

	/*Assign output*/
	p[0]=plist[0]*dbasis[0][0]+plist[1]*dbasis[0][1]+plist[2]*dbasis[0][2]+plist[3]*dbasis[0][3]+plist[4]*dbasis[0][4]+plist[5]*dbasis[0][5];
	p[1]=plist[0]*dbasis[1][0]+plist[1]*dbasis[1][1]+plist[2]*dbasis[1][2]+plist[3]*dbasis[1][3]+plist[4]*dbasis[1][4]+plist[5]*dbasis[1][5];
	p[2]=plist[0]*dbasis[2][0]+plist[1]*dbasis[2][1]+plist[2]*dbasis[2][2]+plist[3]*dbasis[2][3]+plist[4]*dbasis[2][4]+plist[5]*dbasis[2][5];

}
/*}}}*/
/*FUNCTION PentaRef::NumberofNodes{{{*/
int PentaRef::NumberofNodes(void){

	switch(this->element_type){
		case P1Enum:                return NUMNODESP1;
		case P1bubbleEnum:          return NUMNODESP1b;
		case P1bubblecondensedEnum: return NUMNODESP1b;
		case P2Enum:                return NUMNODESP2;
		case P2xP1Enum:             return NUMNODESP2xP1;
		case P1xP2Enum:             return NUMNODESP1xP2;
		case P1P1Enum:              return NUMNODESP1*2;
		case P1P1GLSEnum:           return NUMNODESP1*2;
		case MINIcondensedEnum:     return NUMNODESP1b+NUMNODESP1;
		case MINIEnum:              return NUMNODESP1b+NUMNODESP1;
		case TaylorHoodEnum:        return NUMNODESP2+NUMNODESP1;
		default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
	}

	return -1;
}
/*}}}*/
/*FUNCTION PentaRef::NumberofNodesPressure{{{*/
int PentaRef::NumberofNodesPressure(void){

	switch(this->element_type){
		case P1P1Enum:          return NUMNODESP1;
		case P1P1GLSEnum:       return NUMNODESP1;
		case MINIcondensedEnum: return NUMNODESP1;
		case MINIEnum:          return NUMNODESP1;
		case TaylorHoodEnum:    return NUMNODESP1;
		default: _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
	}

	return -1;
}
/*}}}*/
/*FUNCTION PentaRef::NumberofNodesVelocity{{{*/
int PentaRef::NumberofNodesVelocity(void){

	switch(this->element_type){
		case P1P1Enum:          return NUMNODESP1;
		case P1P1GLSEnum:       return NUMNODESP1;
		case MINIcondensedEnum: return NUMNODESP1b;
		case MINIEnum:          return NUMNODESP1b;
		case TaylorHoodEnum:    return NUMNODESP2;
		default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
	}

	return -1;
}
/*}}}*/
/*FUNCTION PentaRef::VelocityInterpolation{{{*/
int PentaRef::VelocityInterpolation(void){

	switch(this->element_type){
		case P1P1Enum:          return P1Enum;
		case P1P1GLSEnum:       return P1Enum;
		case MINIcondensedEnum: return P1bubbleEnum;
		case MINIEnum:          return P1bubbleEnum;
		case TaylorHoodEnum:    return P2Enum;
		default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
	}

	return -1;
}
/*}}}*/
/*FUNCTION PentaRef::PressureInterpolation{{{*/
int PentaRef::PressureInterpolation(void){

	switch(this->element_type){
		case P1P1Enum:          return P1Enum;
		case P1P1GLSEnum:       return P1Enum;
		case MINIcondensedEnum: return P1Enum;
		case MINIEnum:          return P1Enum;
		case TaylorHoodEnum:    return P1Enum;
		default:       _error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
	}

	return -1;
}
/*}}}*/
/*FUNCTION PentaRef::BasalNodeIndices{{{*/
void PentaRef::BasalNodeIndices(int* pnumindices,int** pindices,int finiteelement){

	/*Output*/
	int  numindices;
	int* indices = NULL;

	switch(finiteelement){
		case P1Enum: case P1DGEnum:
			numindices = 3;
			indices    = xNew<int>(numindices);
			indices[0] = 0;
			indices[1] = 1;
			indices[2] = 2;
			break;
		case P1bubbleEnum: case P1bubblecondensedEnum:
			numindices = 3;
			indices    = xNew<int>(numindices);
			indices[0] = 0;
			indices[1] = 1;
			indices[2] = 2;
			break;
		case P2xP1Enum:
			numindices = 6;
			indices    = xNew<int>(numindices);
			indices[0] = 0;
			indices[1] = 1;
			indices[2] = 2;
			indices[3] = 6;
			indices[4] = 7;
			indices[5] = 8;
			break;
		case P1xP2Enum:
			numindices = 3;
			indices    = xNew<int>(numindices);
			indices[0] = 0;
			indices[1] = 1;
			indices[2] = 2;
			return;
		case P2Enum:
			numindices = 6;
			indices    = xNew<int>(numindices);
			indices[0] = 0;
			indices[1] = 1;
			indices[2] = 2;
			indices[3] = 9;
			indices[4] = 10;
			indices[5] = 11;
			break;
		default:
			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
	}

	/*Assign output pointer*/
	*pnumindices = numindices;
	*pindices    = indices;
}
/*}}}*/
/*FUNCTION PentaRef::SurfaceNodeIndices{{{*/
void PentaRef::SurfaceNodeIndices(int* pnumindices,int** pindices,int finiteelement){

	/*Output*/
	int  numindices;
	int* indices = NULL;

	switch(finiteelement){
		case P1Enum: case P1DGEnum:
			numindices = 3;
			indices    = xNew<int>(numindices);
			indices[0] = 3;
			indices[1] = 4;
			indices[2] = 5;
			break;
		case P1bubbleEnum: case P1bubblecondensedEnum:
			numindices = 3;
			indices    = xNew<int>(numindices);
			indices[0] = 3;
			indices[1] = 4;
			indices[2] = 5;
			break;
		case P2xP1Enum:
			numindices = 6;
			indices    = xNew<int>(numindices);
			indices[0] = 3;
			indices[1] = 4;
			indices[2] = 5;
			indices[3] = 9;
			indices[4] = 10;
			indices[5] = 11;
			break;
		case P1xP2Enum:
			numindices = 3;
			indices    = xNew<int>(numindices);
			indices[0] = 3;
			indices[1] = 4;
			indices[2] = 5;
			return;
		case P2Enum:
			numindices = 6;
			indices    = xNew<int>(numindices);
			indices[0] = 3;
			indices[1] = 4;
			indices[2] = 5;
			indices[3] = 12;
			indices[4] = 13;
			indices[5] = 14;
			break;
		default:
			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
	}

	/*Assign output pointer*/
	*pnumindices = numindices;
	*pindices    = indices;
}
/*}}}*/
