/*!\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 <stdio.h>
#include <string.h>
#include "../objects.h"
#include "../../EnumDefinitions/EnumDefinitions.h"
#include "../../shared/shared.h"
#include "../../Container/Container.h"
#include "../../include/include.h"
/*}}}*/

/*Element macros*/
#define NUMNODESP1    6
#define NUMNODESP1_2d 3
#define NUMNODESMINI  7

/*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){

	_assert_(type==P1Enum || type==P1DGEnum);

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

/*Reference Element numerics*/
/*FUNCTION PentaRef::GetBMacAyealPattyn {{{*/
void PentaRef::GetBMacAyealPattyn(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)=dbasis[0][i]; 
		*(B+NDOF2*NUMNODESP1*0+NDOF2*i+1)=0.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)=(float).5*dbasis[1][i]; 
		*(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=(float).5*dbasis[0][i]; 
	}
}
/*}}}*/
/*FUNCTION PentaRef::GetBMacAyealStokes{{{*/
void PentaRef::GetBMacAyealStokes(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 dh1dh7[3][NUMNODESMINI];
	IssmDouble l1l6[NUMNODESP1];

	/*Get dh1dh6 in actual coordinate system: */
	GetNodalFunctionsMINIDerivatives(&dh1dh7[0][0],xyz_list, gauss);
	GetNodalFunctionsP1(l1l6, gauss);

	/*Build B: */
	for (i=0;i<NUMNODESMINI;i++){
		*(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i)=dh1dh7[0][i]; //B[0][NDOF4*i]=dh1dh6[0][i];
		*(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1)=0;
		*(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2)=0;
		*(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i)=0;
		*(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1)=dh1dh7[1][i];
		*(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2)=0;
		*(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i)=0.5*dh1dh7[1][i];
		*(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1)=0.5*dh1dh7[0][i];
		*(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2)=0;
		*(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i)=0;
		*(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+1)=0;
		*(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+2)=0;
	}

	for (i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
		*(B+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3)=0;
		*(B+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3)=0;
		*(B+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3)=0;
		*(B+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3)=l1l6[i];
	}
}
/*}}}*/
/*FUNCTION PentaRef::GetBPattyn {{{*/
void PentaRef::GetBPattyn(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*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)=dbasis[0][i]; 
		*(B+NDOF2*NUMNODESP1*0+NDOF2*i+1)=0.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)=(float).5*dbasis[1][i]; 
		*(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=(float).5*dbasis[0][i]; 

		*(B+NDOF2*NUMNODESP1*3+NDOF2*i)=(float).5*dbasis[2][i]; 
		*(B+NDOF2*NUMNODESP1*3+NDOF2*i+1)=0.0;

		*(B+NDOF2*NUMNODESP1*4+NDOF2*i)=0.0;
		*(B+NDOF2*NUMNODESP1*4+NDOF2*i+1)=(float).5*dbasis[2][i]; 
	}

}
/*}}}*/
/*FUNCTION PentaRef::GetBprimePattyn {{{*/
void PentaRef::GetBprimePattyn(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss_coord){
	/*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*NUMNODESP1)
	 */
	IssmDouble dbasis[3][NUMNODESP1];

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

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

		*(B+NDOF2*NUMNODESP1*1+NDOF2*i)=dbasis[0][i];
		*(B+NDOF2*NUMNODESP1*1+NDOF2*i+1)=2.0*dbasis[1][i];

		*(B+NDOF2*NUMNODESP1*2+NDOF2*i)=dbasis[1][i]; 
		*(B+NDOF2*NUMNODESP1*2+NDOF2*i+1)=dbasis[0][i]; 

		*(B+NDOF2*NUMNODESP1*3+NDOF2*i)=dbasis[2][i]; 
		*(B+NDOF2*NUMNODESP1*3+NDOF2*i+1)=0.0;

		*(B+NDOF2*NUMNODESP1*4+NDOF2*i)=0.0;
		*(B+NDOF2*NUMNODESP1*4+NDOF2*i+1)=dbasis[2][i]; 
	}
}
/*}}}*/
/*FUNCTION PentaRef::GetBprimeMacAyealStokes{{{*/
void PentaRef::GetBprimeMacAyealStokes(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 dh1dh7[3][NUMNODESMINI];

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

	/*Build Bprime: */
	for (i=0;i<NUMNODESMINI;i++){
		*(Bprime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i)=2*dh1dh7[0][i]; //Bprime[0][NDOF4*i]=dh1dh6[0][i];
		*(Bprime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+1)=dh1dh7[1][i];
		*(Bprime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+2)=0;
		*(Bprime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i)=dh1dh7[0][i];
		*(Bprime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+1)=2*dh1dh7[1][i];
		*(Bprime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+2)=0;
		*(Bprime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i)=dh1dh7[1][i];
		*(Bprime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+1)=dh1dh7[0][i];
		*(Bprime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+2)=0;
	}

	for (i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
		*(Bprime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3)=0;
		*(Bprime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3)=0;
		*(Bprime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3)=0;
	}

}
/*}}}*/
/*FUNCTION PentaRef::GetBStokes {{{*/
void PentaRef::GetBStokes(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 dh1dh7[3][NUMNODESMINI];
	IssmDouble l1l6[NUMNODESP1];

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

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

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

}
/*}}}*/
/*FUNCTION PentaRef::GetBprimeStokes {{{*/
void PentaRef::GetBprimeStokes(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 dh1dh7[3][NUMNODESMINI];
	IssmDouble l1l6[NUMNODESP1];

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

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

	for (i=0;i<NUMNODESP1;i++){ //last column not for the bubble function
		*(B_prime+(NDOF4*NUMNODESP1+3)*0+NDOF4*i+3)=0;
		*(B_prime+(NDOF4*NUMNODESP1+3)*1+NDOF4*i+3)=0;
		*(B_prime+(NDOF4*NUMNODESP1+3)*2+NDOF4*i+3)=0;
		*(B_prime+(NDOF4*NUMNODESP1+3)*3+NDOF4*i+3)=0;
		*(B_prime+(NDOF4*NUMNODESP1+3)*4+NDOF4*i+3)=0;
		*(B_prime+(NDOF4*NUMNODESP1+3)*5+NDOF4*i+3)=0;
		*(B_prime+(NDOF4*NUMNODESP1+3)*6+NDOF4*i+3)=0;
		*(B_prime+(NDOF4*NUMNODESP1+3)*7+NDOF4*i+3)=l1l6[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)
	 */

	/*Same thing in the actual coordinate system: */
	IssmDouble l1l6[6];

	/*Get dh1dh2dh3 in actual coordinates system : */
	GetNodalFunctionsP1(l1l6, gauss);

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

	/*Same thing in the actual coordinate system: */
	IssmDouble dh1dh6[3][NUMNODESP1];

	/*Get dh1dh2dh3 in actual coordinates system : */
	GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);

	/*Build B': */
	for (int i=0;i<NUMNODESP1;i++){
		*(B_conduct+NDOF1*NUMNODESP1*0+NDOF1*i)=dh1dh6[0][i]; 
		*(B_conduct+NDOF1*NUMNODESP1*1+NDOF1*i)=dh1dh6[1][i]; 
		*(B_conduct+NDOF1*NUMNODESP1*2+NDOF1*i)=dh1dh6[2][i]; 
	}
}
/*}}}*/
/*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.*/

	int i;
	IssmDouble dh1dh6[3][NUMNODESP1];

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

	/*Build B: */
	for (i=0;i<NUMNODESP1;i++){
		B[i]=dh1dh6[2][i];  
	}

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

	/*Same thing in the actual coordinate system: */
	IssmDouble dh1dh6[3][NUMNODESP1];

	/*Get dh1dh2dh3 in actual coordinates system : */
	GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],xyz_list,gauss);

	/*Build B': */
	for (int i=0;i<NUMNODESP1;i++){
		*(Bprime_advec+NDOF1*NUMNODESP1*0+NDOF1*i)=dh1dh6[0][i]; 
		*(Bprime_advec+NDOF1*NUMNODESP1*1+NDOF1*i)=dh1dh6[1][i]; 
		*(Bprime_advec+NDOF1*NUMNODESP1*2+NDOF1*i)=dh1dh6[2][i]; 
	}
}
/*}}}*/
/*FUNCTION PentaRef::GetBprimeVert{{{*/
void PentaRef::GetBprimeVert(IssmDouble* B, IssmDouble* xyz_list, GaussPenta* gauss){
	/* Compute Bprime  matrix. Bprime=[L1 L2 L3 L4 L5 L6] where Li is the nodal function for node i*/

	GetNodalFunctionsP1(B, gauss);

}
/*}}}*/
/*FUNCTION PentaRef::GetL{{{*/
void PentaRef::GetL(IssmDouble* L, GaussPenta* gauss, int numdof){
	/*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: 
	 **       numdof=1: 
	 **                 Li=h;
	 **       numdof=2:
	 **                 Li=[ h   0 ]
	 **                    [ 0   h ]
	 ** where h is the interpolation function for node i.
	 **
	 ** We assume L has been allocated already, of size: NUMNODESP1 (numdof=1), or numdofx(numdof*NUMNODESP1) (numdof=2)
	 **/

	int i;
	IssmDouble l1l6[6];

	/*Get l1l6 in actual coordinate system: */
	GetNodalFunctionsP1(l1l6,gauss);

	/*Build L: */
	if(numdof==1){
		for (i=0;i<NUMNODESP1;i++){
			L[i]=l1l6[i]; 
		}
	}
	else{
		for (i=0;i<NUMNODESP1;i++){
			*(L+numdof*NUMNODESP1*0+numdof*i)=l1l6[i]; 
			*(L+numdof*NUMNODESP1*0+numdof*i+1)=0;
			*(L+numdof*NUMNODESP1*1+numdof*i)=0;
			*(L+numdof*NUMNODESP1*1+numdof*i+1)=l1l6[i];
		}
	}
} 
/*}}}*/
/*FUNCTION PentaRef::GetLStokes{{{*/
void PentaRef::GetLStokes(IssmDouble* LStokes, 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 ]
	 *		      [ 0 0 ]
	 *		      [ 0 0 ]
	 * where h is the interpolation function for node i.
	 */

	const int num_dof=4;
	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 LStokes: */
	for (int i=0;i<3;i++){
		*(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+0)=l1l2l3[i];
		*(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0.;
		*(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+2)=0.;
		*(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+3)=0.;

		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+0)=0.;
		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+2)=0.;
		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+3)=0.;
	}
}
/*}}}*/
/*FUNCTION PentaRef::GetLprimeStokes {{{*/
void PentaRef::GetLprimeStokes(IssmDouble* LprimeStokes, 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 dh1dh6[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(&dh1dh6[0][0],xyz_list,gauss);

	/*Build LprimeStokes: */
	for (i=0;i<3;i++){
		*(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i)=l1l2l3[i]; //LprimeStokes[0][NDOF2*i]=dh1dh3[0][i];
		*(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+2)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+3)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
		*(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+2)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+3)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i)=l1l2l3[i];
		*(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+2)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+3)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=l1l2l3[i];
		*(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+2)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+3)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+1)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+2)=l1l2l3[i];
		*(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+3)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+1)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+2)=l1l2l3[i];
		*(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+3)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+1)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+2)=dh1dh6[2][i];
		*(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+3)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+1)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+2)=dh1dh6[2][i];
		*(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+3)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*8+num_dof*i)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*8+num_dof*i+1)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*8+num_dof*i+2)=dh1dh6[2][i];
		*(LprimeStokes+num_dof*NUMNODESP1_2d*8+num_dof*i+3)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*9+num_dof*i)=dh1dh6[2][i];
		*(LprimeStokes+num_dof*NUMNODESP1_2d*9+num_dof*i+1)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*9+num_dof*i+2)=dh1dh6[0][i];
		*(LprimeStokes+num_dof*NUMNODESP1_2d*9+num_dof*i+3)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*10+num_dof*i)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*10+num_dof*i+1)=dh1dh6[2][i];
		*(LprimeStokes+num_dof*NUMNODESP1_2d*10+num_dof*i+2)=dh1dh6[1][i];
		*(LprimeStokes+num_dof*NUMNODESP1_2d*10+num_dof*i+3)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*11+num_dof*i)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*11+num_dof*i+1)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*11+num_dof*i+2)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*11+num_dof*i+3)=l1l2l3[i];
		*(LprimeStokes+num_dof*NUMNODESP1_2d*12+num_dof*i)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*12+num_dof*i+1)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*12+num_dof*i+2)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*12+num_dof*i+3)=l1l2l3[i];
		*(LprimeStokes+num_dof*NUMNODESP1_2d*13+num_dof*i)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*13+num_dof*i+1)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*13+num_dof*i+2)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*13+num_dof*i+3)=l1l2l3[i];
	}
}
/*}}}*/
/*FUNCTION PentaRef::GetLMacAyealStokes {{{*/
void PentaRef::GetLMacAyealStokes(IssmDouble* LStokes, 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 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 LStokes: */
	for (i=0;i<3;i++){
		*(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh3[0][i];
		*(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0;
		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i)=0;
		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
		*(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i)=l1l2l3[i];
		*(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0;
		*(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i)=0;
		*(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=l1l2l3[i];
		*(LStokes+num_dof*NUMNODESP1_2d*4+num_dof*i)=l1l2l3[i];
		*(LStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+1)=0;
		*(LStokes+num_dof*NUMNODESP1_2d*5+num_dof*i)=0;
		*(LStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+1)=l1l2l3[i];
		*(LStokes+num_dof*NUMNODESP1_2d*6+num_dof*i)=l1l2l3[i];
		*(LStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+1)=0;
		*(LStokes+num_dof*NUMNODESP1_2d*7+num_dof*i)=0;
		*(LStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+1)=l1l2l3[i];

	}
}
/*}}}*/
/*FUNCTION PentaRef::GetLprimeMacAyealStokes {{{*/
void PentaRef::GetLprimeMacAyealStokes(IssmDouble* LprimeStokes, 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 i;
	int num_dof=4;

	IssmDouble l1l2l3[NUMNODESP1_2d];
	IssmDouble dh1dh6[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(&dh1dh6[0][0],xyz_list,gauss);

	/*Build LprimeStokes: */
	for (i=0;i<3;i++){
		*(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i)=l1l2l3[i]; //LprimeStokes[0][NDOF2*i]=dh1dh3[0][i];
		*(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+2)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+3)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
		*(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+2)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+3)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+2)=l1l2l3[i];
		*(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+3)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+2)=l1l2l3[i];
		*(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+3)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+1)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+2)=dh1dh6[2][i];
		*(LprimeStokes+num_dof*NUMNODESP1_2d*4+num_dof*i+3)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+1)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+2)=dh1dh6[2][i];
		*(LprimeStokes+num_dof*NUMNODESP1_2d*5+num_dof*i+3)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+1)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+2)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*6+num_dof*i+3)=l1l2l3[i];
		*(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+1)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+2)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*7+num_dof*i+3)=l1l2l3[i];
	}
}
/*}}}*/
/*FUNCTION PentaRef::GetLStokesMacAyeal {{{*/
void PentaRef::GetLStokesMacAyeal(IssmDouble* LStokes, 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]
	 *		      [ 0    0    h   0]
	 *		      [ 0    0    h   0]
	 * where h is the interpolation function for node i.
	 */

	int i;
	int num_dof=4;

	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 LStokes: */
	for (i=0;i<3;i++){
		*(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i)=l1l2l3[i]; //LStokes[0][NDOF2*i]=dh1dh3[0][i];
		*(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0;
		*(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+2)=0;
		*(LStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+3)=0;
		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i)=0;
		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+2)=0;
		*(LStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+3)=0;
		*(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i)=0;
		*(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0;
		*(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+2)=l1l2l3[i];
		*(LStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+3)=0;
		*(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i)=0;
		*(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=0;
		*(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+2)=l1l2l3[i];
		*(LStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+3)=0;

	}
}
/*}}}*/
/*FUNCTION PentaRef::GetLprimeStokesMacAyeal {{{*/
void PentaRef::GetLprimeStokesMacAyeal(IssmDouble* LprimeStokes, 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 i;
	int num_dof=2;

	IssmDouble l1l2l3[NUMNODESP1_2d];
	IssmDouble dh1dh6[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(&dh1dh6[0][0],xyz_list,gauss);

	/*Build LprimeStokes: */
	for (i=0;i<3;i++){
		*(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i)=l1l2l3[i]; //LprimeStokes[0][NDOF2*i]=dh1dh3[0][i];
		*(LprimeStokes+num_dof*NUMNODESP1_2d*0+num_dof*i+1)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*1+num_dof*i+1)=l1l2l3[i];
		*(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i)=l1l2l3[i];
		*(LprimeStokes+num_dof*NUMNODESP1_2d*2+num_dof*i+1)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i)=0;
		*(LprimeStokes+num_dof*NUMNODESP1_2d*3+num_dof*i+1)=l1l2l3[i];
	}
}
/*}}}*/
/*FUNCTION PentaRef::GetJacobian {{{*/
void PentaRef::GetJacobian(IssmDouble* J, IssmDouble* xyz_list,GaussPenta* gauss){

	int i,j;

	/*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) _error2_("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,x2,x3,y1,y2,y3,z1,z2,z3;

	x1=*(xyz_list+3*0+0);
	y1=*(xyz_list+3*0+1);
	z1=*(xyz_list+3*0+2);
	x2=*(xyz_list+3*1+0);
	y2=*(xyz_list+3*1+1);
	z2=*(xyz_list+3*1+2);
	x3=*(xyz_list+3*2+0);
	y3=*(xyz_list+3*2+1);
	z3=*(xyz_list+3*2+2);

	/*Jdet = norm( AB ^ AC ) / (2 * area of the reference triangle), with areaRef=sqrt(3) */
	*Jdet=SQRT3/6.0*pow(pow(((y2-y1)*(z3-z1)-(z2-z1)*(y3-y1)),2.0)+pow(((z2-z1)*(x3-x1)-(x2-x1)*(z3-z1)),2.0)+pow(((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1)),2.0),0.5);
	if(*Jdet<0) _error2_("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,x2,y1,y2,z1,z2;

	x1=*(xyz_list+3*0+0);
	y1=*(xyz_list+3*0+1);
	z1=*(xyz_list+3*0+2);
	x2=*(xyz_list+3*1+0);
	y2=*(xyz_list+3*1+1);
	z2=*(xyz_list+3*1+2);

	*Jdet=1.0/2.0*sqrt(pow(x2-x1,2.) + pow(y2-y1,2.) + pow(z2-z1,2.));
	if(*Jdet<0) _error2_("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::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* dh1dh7,IssmDouble* xyz_list, GaussPenta* gauss){

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

	int       i;
	IssmDouble    dh1dh7_ref[3][NUMNODESMINI];
	IssmDouble    Jinv[3][3];

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

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

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

	for (i=0;i<NUMNODESMINI;i++){
		*(dh1dh7+NUMNODESMINI*0+i)=Jinv[0][0]*dh1dh7_ref[0][i]+Jinv[0][1]*dh1dh7_ref[1][i]+Jinv[0][2]*dh1dh7_ref[2][i];
		*(dh1dh7+NUMNODESMINI*1+i)=Jinv[1][0]*dh1dh7_ref[0][i]+Jinv[1][1]*dh1dh7_ref[1][i]+Jinv[1][2]*dh1dh7_ref[2][i];
		*(dh1dh7+NUMNODESMINI*2+i)=Jinv[2][0]*dh1dh7_ref[0][i]+Jinv[2][1]*dh1dh7_ref[1][i]+Jinv[2][2]*dh1dh7_ref[2][i];
	}

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

	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
	 * natural coordinate system) at the gaussian point. */
	IssmDouble r=gauss->coord2-gauss->coord1;
	IssmDouble s=-3.0/SQRT3*(gauss->coord1+gauss->coord2-2.0/3.0);
	IssmDouble zeta=gauss->coord4;

	/*First nodal function: */
	*(dl1dl7+NUMNODESMINI*0+0)=-0.5*(1.0-zeta)/2.0;
	*(dl1dl7+NUMNODESMINI*1+0)=-SQRT3/6.0*(1.0-zeta)/2.0;
	*(dl1dl7+NUMNODESMINI*2+0)=-0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD);

	/*Second nodal function: */
	*(dl1dl7+NUMNODESMINI*0+1)=0.5*(1.0-zeta)/2.0;
	*(dl1dl7+NUMNODESMINI*1+1)=-SQRT3/6.0*(1.0-zeta)/2.0;
	*(dl1dl7+NUMNODESMINI*2+1)=-0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD);

	/*Third nodal function: */
	*(dl1dl7+NUMNODESMINI*0+2)=0;
	*(dl1dl7+NUMNODESMINI*1+2)=SQRT3/3.0*(1.0-zeta)/2.0;
	*(dl1dl7+NUMNODESMINI*2+2)=-0.5*(SQRT3/3.0*s+ONETHIRD);

	/*Fourth nodal function: */
	*(dl1dl7+NUMNODESMINI*0+3)=-0.5*(1.0+zeta)/2.0;
	*(dl1dl7+NUMNODESMINI*1+3)=-SQRT3/6.0*(1.0+zeta)/2.0;
	*(dl1dl7+NUMNODESMINI*2+3)=0.5*(-0.5*r-SQRT3/6.0*s+ONETHIRD);

	/*Fith nodal function: */
	*(dl1dl7+NUMNODESMINI*0+4)=0.5*(1.0+zeta)/2.0;
	*(dl1dl7+NUMNODESMINI*1+4)=-SQRT3/6.0*(1.0+zeta)/2.0;
	*(dl1dl7+NUMNODESMINI*2+4)=0.5*(0.5*r-SQRT3/6.0*s+ONETHIRD);

	/*Sixth nodal function: */
	*(dl1dl7+NUMNODESMINI*0+5)=0;
	*(dl1dl7+NUMNODESMINI*1+5)=SQRT3/3.0*(1.0+zeta)/2.0;
	*(dl1dl7+NUMNODESMINI*2+5)=0.5*(SQRT3/3.0*s+ONETHIRD);

	/*Seventh nodal function: */
	*(dl1dl7+NUMNODESMINI*0+6)=9.0/2.0*r*(1.0+zeta)*(zeta-1.0)*(SQRT3*s+1.0);
	*(dl1dl7+NUMNODESMINI*1+6)=9.0/4.0*(1+zeta)*(1-zeta)*(SQRT3*pow(s,2.0)-2.0*s-SQRT3*pow(r,2.0));
	*(dl1dl7+NUMNODESMINI*2+6)=27*gauss->coord1*gauss->coord2*gauss->coord3*(-2.0*zeta);

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

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

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

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

	/*Get derivative values with respect to parametric coordinate system: */
	GetNodalFunctionsP1DerivativesReference(&dh1dh6_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++){
		*(dh1dh6+NUMNODESP1*0+i)=Jinv[0][0]*dh1dh6_ref[0][i]+Jinv[0][1]*dh1dh6_ref[1][i]+Jinv[0][2]*dh1dh6_ref[2][i];
		*(dh1dh6+NUMNODESP1*1+i)=Jinv[1][0]*dh1dh6_ref[0][i]+Jinv[1][1]*dh1dh6_ref[1][i]+Jinv[1][2]*dh1dh6_ref[2][i];
		*(dh1dh6+NUMNODESP1*2+i)=Jinv[2][0]*dh1dh6_ref[0][i]+Jinv[2][1]*dh1dh6_ref[1][i]+Jinv[2][2]*dh1dh6_ref[2][i];
	}

}
/*}}}*/
/*FUNCTION PentaRef::GetNodalFunctionsP1DerivativesReference {{{*/
void PentaRef::GetNodalFunctionsP1DerivativesReference(IssmDouble* dl1dl6,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 A1,A2,A3,z;

	A1=gauss->coord1; _assert_(A1>=0 && A1<=1);//first area coordinate value. In term of xi and eta: A1=(1-xi)/2-eta/(2*SQRT3);
	A2=gauss->coord2; _assert_(A2>=0 && A2<=1);//second area coordinate value In term of xi and eta: A2=(1+xi)/2-eta/(2*SQRT3);
	A3=gauss->coord3; _assert_(A3>=0 && A3<=1);//third area coordinate value  In term of xi and eta: A3=y/SQRT3;
	z =gauss->coord4; _assert_(z>=-1 &&  z<=1);//fourth vertical coordinate value. Corresponding nodal function: (1-z)/2 and (1+z)/2

	/*First nodal function derivatives. The corresponding nodal function is N=A1*(1-z)/2. Its derivatives follow*/
	*(dl1dl6+NUMNODESP1*0+0)=-0.5*(1.0-z)/2.0;
	*(dl1dl6+NUMNODESP1*1+0)=-0.5/SQRT3*(1.0-z)/2.0;
	*(dl1dl6+NUMNODESP1*2+0)=-0.5*A1;

	/*Second nodal function: The corresponding nodal function is N=A2*(1-z)/2. Its derivatives follow*/
	*(dl1dl6+NUMNODESP1*0+1)=0.5*(1.0-z)/2.0;
	*(dl1dl6+NUMNODESP1*1+1)=-0.5/SQRT3*(1.0-z)/2.0;
	*(dl1dl6+NUMNODESP1*2+1)=-0.5*A2;

	/*Third nodal function: The corresponding nodal function is N=A3*(1-z)/2. Its derivatives follow*/
	*(dl1dl6+NUMNODESP1*0+2)=0.0;
	*(dl1dl6+NUMNODESP1*1+2)=1.0/SQRT3*(1.0-z)/2.0;
	*(dl1dl6+NUMNODESP1*2+2)=-0.5*A3;

	/*Fourth nodal function: The corresponding nodal function is N=A1*(1+z)/2. Its derivatives follow*/
	*(dl1dl6+NUMNODESP1*0+3)=-0.5*(1.0+z)/2.0;
	*(dl1dl6+NUMNODESP1*1+3)=-0.5/SQRT3*(1.0+z)/2.0;
	*(dl1dl6+NUMNODESP1*2+3)=0.5*A1;

	/*Fifth nodal function: The corresponding nodal function is N=A2*(1+z)/2. Its derivatives follow*/
	*(dl1dl6+NUMNODESP1*0+4)=0.5*(1.0+z)/2.0;
	*(dl1dl6+NUMNODESP1*1+4)=-0.5/SQRT3*(1.0+z)/2.0;
	*(dl1dl6+NUMNODESP1*2+4)=0.5*A2;

	/*Sixth nodal function: The corresponding nodal function is N=A3*(1+z)/2. Its derivatives follow*/
	*(dl1dl6+NUMNODESP1*0+5)=0.0;
	*(dl1dl6+NUMNODESP1*1+5)=1.0/SQRT3*(1.0+z)/2.0;
	*(dl1dl6+NUMNODESP1*2+5)=0.5*A3;
}
/*}}}*/
/*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= pow(pow(x2-x1,2.) + pow(y2-y1,2.),0.5) * (z4-z1 + z3-z2)/8;
	if(*Jdet<0) _error2_("negative jacobian determinant!");

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

	/*intermediary*/
	IssmDouble l1l6[6];

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

	/*Assign output pointers:*/
	*pvalue=l1l6[0]*plist[0]+l1l6[1]*plist[1]+l1l6[2]*plist[2]+l1l6[3]*plist[3]+l1l6[4]*plist[4]+l1l6[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 dh1dh6[3][NUMNODESP1];

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

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

}
/*}}}*/
