Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 17094)
+++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 17095)
@@ -56,1110 +56,4 @@
 
 /*Reference Element numerics*/
-/*FUNCTION PentaRef::GetBSSAHO {{{*/
-void PentaRef::GetBSSAHO(IssmDouble* B, IssmDouble* xyz_list, Gauss* 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, Gauss* 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, Gauss* 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,Gauss* 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, Gauss* 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, Gauss* 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, Gauss* 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/dz    ]
-	 *					[ 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, Gauss* 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, Gauss* 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, Gauss* 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, Gauss* 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, Gauss* 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, Gauss* 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, Gauss* 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, Gauss* gauss){
-
-	GetNodalFunctions(B,gauss);
-
-}
-/*}}}*/
-/*FUNCTION PentaRef::GetBHOFriction{{{*/
-void PentaRef::GetBHOFriction(IssmDouble* B, Gauss* 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, Gauss* 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, Gauss* gauss_in){
-	/* 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        num_dof=4;
-	IssmDouble L1L2l3[NUMNODESP1_2d];
-	IssmDouble dbasis[3][NUMNODESP1];
-
-	/*Cast gauss to GaussPenta*/
-	_assert_(gauss_in->Enum()==GaussPentaEnum);
-	GaussPenta* gauss = dynamic_cast<GaussPenta*>(gauss_in);
-
-	/*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, Gauss* gauss_in){
-	/*
-	 * 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];
-
-	/*Cast gauss to GaussPenta*/
-	_assert_(gauss_in->Enum()==GaussPentaEnum);
-	GaussPenta* gauss = dynamic_cast<GaussPenta*>(gauss_in);
-
-	/*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, Gauss* gauss_in){
-	/* 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];
-
-	/*Cast gauss to GaussPenta*/
-	_assert_(gauss_in->Enum()==GaussPentaEnum);
-	GaussPenta* gauss = dynamic_cast<GaussPenta*>(gauss_in);
-
-	/*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, Gauss* gauss_in){
-	/* 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];
-
-	/*Cast gauss to GaussPenta*/
-	_assert_(gauss_in->Enum()==GaussPentaEnum);
-	GaussPenta* gauss = dynamic_cast<GaussPenta*>(gauss_in);
-
-	/*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, Gauss* gauss_in){
-	/* 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];
-
-	/*Cast gauss to GaussPenta*/
-	_assert_(gauss_in->Enum()==GaussPentaEnum);
-	GaussPenta* gauss = dynamic_cast<GaussPenta*>(gauss_in);
-
-	/*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,Gauss* gauss_in){
Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h	(revision 17094)
+++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h	(revision 17095)
@@ -41,25 +41,4 @@
 		void GetSegmentJacobianDeterminant(IssmDouble*  Jdet, IssmDouble* xyz_list,Gauss* gauss);
 		void GetJacobianInvert(IssmDouble*  Jinv, IssmDouble* xyz_list,Gauss* gauss);
-		void GetBSSAHO(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss);
-		void GetBSSAFS(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss);
-		void GetBHO(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss);
-		void GetBFSstrainrate(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss);
-		void GetBFS(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss);
-		void GetBFSGLS(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss);
-		void GetBprimeSSAFS(IssmDouble* Bprime, IssmDouble* xyz_list, Gauss* gauss);
-		void GetBprimeHO(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss);
-		void GetBprimeFS(IssmDouble* B_prime, IssmDouble* xyz_list, Gauss* gauss);
-		void GetBprimeFSGLS(IssmDouble* B_prime, IssmDouble* xyz_list, Gauss* gauss);
-		void GetBprimeVert(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss);
-		void GetBAdvec(IssmDouble* B_advec, IssmDouble* xyz_list, Gauss* gauss);
-		void GetBConduct(IssmDouble* B_conduct, IssmDouble* xyz_list, Gauss* gauss);
-		void GetBVert(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss);
-		void GetBprimeAdvec(IssmDouble* Bprime_advec, IssmDouble* xyz_list, Gauss* gauss);
-		void GetBHOFriction(IssmDouble* L, Gauss* gauss);
-		void GetLFS(IssmDouble* LFS, Gauss* gauss);
-		void GetLprimeFS(IssmDouble* LprimeFS, IssmDouble* xyz_list, Gauss* gauss);
-		void GetLSSAFS(IssmDouble* LSSAFS, Gauss* gauss);
-		void GetLprimeSSAFS(IssmDouble* LprimeSSAFS, IssmDouble* xyz_list, Gauss* gauss);
-		void GetLFSSSA(IssmDouble* LFSSSA, Gauss* gauss);
 		void GetLprimeFSSSA(IssmDouble* LprimeFSSSA, IssmDouble* xyz_list, Gauss* gauss);
 		void GetInputValue(IssmDouble* pvalue,IssmDouble* plist, Gauss* gauss);
Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 17094)
+++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.cpp	(revision 17095)
@@ -51,102 +51,4 @@
 
 /*Reference Element numerics*/
-/*FUNCTION TriaRef::GetBHydro {{{*/
-void TriaRef::GetBHydro(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss){
-	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
-	 * For node i, Bi can be expressed in the actual coordinate system
-	 * by: 
-	 *       Bi=[ dN/dx ]
-	 *          [ dN/dy ]
-	 * where N is the finiteelement function for node i.
-	 *
-	 * We assume B has been allocated already, of size: 3x(NDOF2*numnodes)
-	 */
-
-	/*Fetch number of nodes for this finite element*/
-	int numnodes = this->NumberofNodes();
-
-	/*Get nodal functions derivatives*/
-	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
-	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-
-	/*Build B: */
-	for(int i=0;i<numnodes;i++){
-		B[numnodes*0+i] = dbasis[0*numnodes+i];
-		B[numnodes*1+i] = dbasis[1*numnodes+i];
-	}
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(dbasis);
-}
-/*}}}*/
-/*FUNCTION TriaRef::GetBSSA {{{*/
-void TriaRef::GetBSSA(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss){
-	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
-	 * For node i, Bi can be expressed in the actual coordinate system
-	 * by: 
-	 *       Bi=[ dN/dx           0    ]
-	 *          [   0           dN/dy  ]
-	 *          [ 1/2*dN/dy  1/2*dN/dx ]
-	 * where N is the finiteelement function for node i.
-	 *
-	 * We assume B has been allocated already, of size: 3x(NDOF2*numnodes)
-	 */
-
-	/*Fetch number of nodes for this finite element*/
-	int numnodes = this->NumberofNodes();
-
-	/*Get nodal functions derivatives*/
-	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
-	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-
-	/*Build B: */
-	for(int i=0;i<numnodes;i++){
-		B[NDOF2*numnodes*0+NDOF2*i+0] = dbasis[0*numnodes+i];
-		B[NDOF2*numnodes*0+NDOF2*i+1] = 0.;
-		B[NDOF2*numnodes*1+NDOF2*i+0] = 0.;
-		B[NDOF2*numnodes*1+NDOF2*i+1] = dbasis[1*numnodes+i];
-		B[NDOF2*numnodes*2+NDOF2*i+0] = .5*dbasis[1*numnodes+i];
-		B[NDOF2*numnodes*2+NDOF2*i+1] = .5*dbasis[0*numnodes+i];
-	}
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(dbasis);
-}
-/*}}}*/
-/*FUNCTION TriaRef::GetBSSAFS {{{*/
-void TriaRef::GetBSSAFS(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss){
-
-	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
-	 * For node i, Bi can be expressed in the actual coordinate system
-	 * by: 
-	 *       Bi=[   dN/dx         0     ]
-	 *          [       0       dN/dy   ]
-	 *          [  1/2*dN/dy  1/2*dN/dx ]
-	 * where N is the finiteelement function for node i.
-	 *
-	 * We assume B has been allocated already, of size: 3x(NDOF2*numnodes)
-	 */
-
-	/*Fetch number of nodes for this finite element*/
-	int numnodes = this->NumberofNodes();
-
-	/*Get nodal functions derivatives*/
-	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
-	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-
-	/*Build B': */
-	for(int i=0;i<numnodes;i++){
-		B[NDOF2*numnodes*0+NDOF2*i+0] = dbasis[0*numnodes+i];
-		B[NDOF2*numnodes*0+NDOF2*i+1] = 0.;
-		B[NDOF2*numnodes*1+NDOF2*i+0] = 0.;
-		B[NDOF2*numnodes*1+NDOF2*i+1] = dbasis[1*numnodes+i];
-		B[NDOF2*numnodes*2+NDOF2*i+0] = 0.5*dbasis[1*numnodes+i];
-		B[NDOF2*numnodes*2+NDOF2*i+1] = 0.5*dbasis[0*numnodes+i];
-	}
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(dbasis);
-}
-/*}}}*/
 /*FUNCTION TriaRef::GetSegmentBFlux{{{*/
 void TriaRef::GetSegmentBFlux(IssmDouble* B,Gauss* gauss, int index1,int index2){
@@ -199,329 +101,4 @@
 	/*Clean-up*/
 	xDelete<IssmDouble>(basis);
-}
-/*}}}*/
-/*FUNCTION TriaRef::GetBExtrusion{{{*/
-void TriaRef::GetBExtrusion(IssmDouble* B, IssmDouble* xyz_list, Gauss* 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>(2*numnodes);
-	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-
-	/*Build B: */
-	for(int i=0;i<numnodes;i++){
-		B[i] = dbasis[1*numnodes+i];  
-	}
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(dbasis);
-}
-/*}}}*/
-/*FUNCTION TriaRef::GetBFS {{{*/
-void TriaRef::GetBFS(IssmDouble* B, IssmDouble* xyz_list, Gauss* 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=[ dphi/dx          0        ]
-	 *					 [   0           dphi/dy     ]
-	 *					 [ 1/2*dphi/dy    1/2*dphi/dx]
-	 *					 [   0             0         ]
-	 *					 [ dphi/dx         dphi/dy   ]
-	 *
-	 * by: 	  Bpi=[  0    ]
-	 *					[  0    ]
-	 *					[  0    ]
-	 *					[ phi_p ]
-	 *					[  0    ]
-	 *	where phi is the finiteelement 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>(2*vnumnodes);
-	IssmDouble* pbasis =xNew<IssmDouble>(pnumnodes);
-	GetNodalFunctionsDerivativesVelocity(vdbasis,xyz_list,gauss);
-	GetNodalFunctionsPressure(pbasis,gauss);
-
-	/*Build B: */
-	for(int i=0;i<vnumnodes;i++){
-		B[(2*vnumnodes+pnumnodes)*0+2*i+0] = vdbasis[0*vnumnodes+i];
-		B[(2*vnumnodes+pnumnodes)*0+2*i+1] = 0.;
-		B[(2*vnumnodes+pnumnodes)*1+2*i+0] = 0.;
-		B[(2*vnumnodes+pnumnodes)*1+2*i+1] = vdbasis[1*vnumnodes+i];
-		B[(2*vnumnodes+pnumnodes)*2+2*i+0] = .5*vdbasis[1*vnumnodes+i];
-		B[(2*vnumnodes+pnumnodes)*2+2*i+1] = .5*vdbasis[0*vnumnodes+i];
-		B[(2*vnumnodes+pnumnodes)*3+2*i+0] = 0.;
-		B[(2*vnumnodes+pnumnodes)*3+2*i+1] = 0.;
-		B[(2*vnumnodes+pnumnodes)*4+2*i+0] = vdbasis[0*vnumnodes+i];
-		B[(2*vnumnodes+pnumnodes)*4+2*i+1] = vdbasis[1*vnumnodes+i];
-	}
-	for(int i=0;i<pnumnodes;i++){
-		B[(2*vnumnodes+pnumnodes)*0+(2*vnumnodes)+i] = 0.;
-		B[(2*vnumnodes+pnumnodes)*1+(2*vnumnodes)+i] = 0.;
-		B[(2*vnumnodes+pnumnodes)*2+(2*vnumnodes)+i] = 0.;
-		B[(2*vnumnodes+pnumnodes)*3+(2*vnumnodes)+i] = pbasis[i];
-		B[(2*vnumnodes+pnumnodes)*4+(2*vnumnodes)+i] = 0.;
-	}
-
-	/*Clean up*/
-	xDelete<IssmDouble>(vdbasis);
-	xDelete<IssmDouble>(pbasis);
-}
-/*}}}*/
-/*FUNCTION TriaRef::GetBprimeFS {{{*/
-void TriaRef::GetBprimeFS(IssmDouble* B_prime, IssmDouble* xyz_list, Gauss* 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' = [  dphi/dx     0     ]
-	 *					 [     0      dphi/dy ]
-	 *					 [  dphi/dy   dphi/dx ]
-	 *					 [  dphi/dx   dphi/dy ]
-	 *					 [     0      0       ]
-	 *
-	 * by: 	  Bpi=[  0  ]
-	 *					[  0  ]
-	 *					[  0  ]
-	 *					[  0  ]
-	 *					[ phi ]
-	 *	where phi is the finiteelement 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>(2*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[(2*vnumnodes+pnumnodes)*0+2*i+0] = vdbasis[0*vnumnodes+i];
-		B_prime[(2*vnumnodes+pnumnodes)*0+2*i+1] = 0.;
-		B_prime[(2*vnumnodes+pnumnodes)*1+2*i+0] = 0.;
-		B_prime[(2*vnumnodes+pnumnodes)*1+2*i+1] = vdbasis[1*vnumnodes+i];
-		B_prime[(2*vnumnodes+pnumnodes)*2+2*i+0] = vdbasis[1*vnumnodes+i];
-		B_prime[(2*vnumnodes+pnumnodes)*2+2*i+1] = vdbasis[0*vnumnodes+i];
-		B_prime[(2*vnumnodes+pnumnodes)*3+2*i+0] = vdbasis[0*vnumnodes+i];
-		B_prime[(2*vnumnodes+pnumnodes)*3+2*i+1] = vdbasis[1*vnumnodes+i];
-		B_prime[(2*vnumnodes+pnumnodes)*4+2*i+0] = 0.;
-		B_prime[(2*vnumnodes+pnumnodes)*4+2*i+1] = 0.;
-	}
-	for(int i=0;i<pnumnodes;i++){
-		B_prime[(2*vnumnodes+pnumnodes)*0+(2*vnumnodes)+i] = 0.;
-		B_prime[(2*vnumnodes+pnumnodes)*1+(2*vnumnodes)+i] = 0.;
-		B_prime[(2*vnumnodes+pnumnodes)*2+(2*vnumnodes)+i] = 0.;
-		B_prime[(2*vnumnodes+pnumnodes)*3+(2*vnumnodes)+i] = 0.;
-		B_prime[(2*vnumnodes+pnumnodes)*4+(2*vnumnodes)+i] = pbasis[i];
-	}
-
-	/*Clean up*/
-	xDelete<IssmDouble>(vdbasis);
-	xDelete<IssmDouble>(pbasis);
-}
-/*}}}*/
-/*FUNCTION TriaRef::GetBMasstransport{{{*/
-void TriaRef::GetBMasstransport(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss){
-	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
-	 * For node i, Bi can be expressed in the actual coordinate system
-	 * by: 
-	 *       Bi=[ N ]
-	 *          [ N ]
-	 * where N is the finiteelement function for node i.
-	 *
-	 * We assume B_prog has been allocated already, of size: 2x(NDOF1*numnodes)
-	 */
-
-	/*Fetch number of nodes for this finite element*/
-	int numnodes = this->NumberofNodes();
-
-	/*Get nodal functions*/
-	IssmDouble* basis=xNew<IssmDouble>(numnodes);
-	GetNodalFunctions(basis,gauss);
-
-	/*Build B: */
-	for(int i=0;i<numnodes;i++){
-		B[numnodes*0+i] = basis[i];
-		B[numnodes*1+i] = basis[i];
-	}
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(basis);
-}
-/*}}}*/
-/*FUNCTION TriaRef::GetBprimeSSA {{{*/
-void TriaRef::GetBprimeSSA(IssmDouble* Bprime, IssmDouble* xyz_list, Gauss* gauss){
-
-	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
-	 * For node i, Bi' can be expressed in the actual coordinate system
-	 * by: 
-	 *       Bi_prime=[ 2*dN/dx    dN/dy ]
-	 *                [   dN/dx  2*dN/dy ]
-	 *                [   dN/dy    dN/dx ]
-	 * where hNis the finiteelement function for node i.
-	 *
-	 * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
-	 */
-
-	/*Fetch number of nodes for this finite element*/
-	int numnodes = this->NumberofNodes();
-
-	/*Get nodal functions derivatives*/
-	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
-	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-
-	/*Build B': */
-	for(int i=0;i<numnodes;i++){
-		Bprime[NDOF2*numnodes*0+NDOF2*i+0] = 2.*dbasis[0*numnodes+i];
-		Bprime[NDOF2*numnodes*0+NDOF2*i+1] =    dbasis[1*numnodes+i];
-		Bprime[NDOF2*numnodes*1+NDOF2*i+0] =    dbasis[0*numnodes+i];
-		Bprime[NDOF2*numnodes*1+NDOF2*i+1] = 2.*dbasis[1*numnodes+i];
-		Bprime[NDOF2*numnodes*2+NDOF2*i+0] =    dbasis[1*numnodes+i];
-		Bprime[NDOF2*numnodes*2+NDOF2*i+1] =    dbasis[0*numnodes+i];
-	}
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(dbasis);
-}
-/*}}}*/
-/*FUNCTION TriaRef::GetBprimeSSAFS {{{*/
-void TriaRef::GetBprimeSSAFS(IssmDouble* Bprime, IssmDouble* xyz_list, Gauss* gauss){
-	/*Compute Bprime  matrix. Bprime=[Bprime1 Bprime2 Bprime3] where Bprimei is of size 3*NDOF2. 
-	 * For node i, Bprimei can be expressed in the actual coordinate system
-	 * by: 
-	 *       Bprimei=[  dN/dx    0   ]
-	 *               [    0    dN/dy ]
-	 *               [  dN/dy  dN/dx ]
-	 N               [  dN/dx  dN/dy ]
-	 * where N is the finiteelement function for node i.
-	 *
-	 * We assume Bprime has been allocated already, of size: 3x(NDOF2*numnodes)
-	 */
-
-	/*Fetch number of nodes for this finite element*/
-	int numnodes = this->NumberofNodes();
-
-	/*Get nodal functions*/
-	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
-	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-
-	/*Build Bprime: */
-	for(int i=0;i<numnodes;i++){
-		Bprime[NDOF2*numnodes*0+NDOF2*i+0] = dbasis[0*numnodes+i];
-		Bprime[NDOF2*numnodes*0+NDOF2*i+1] = 0.;
-		Bprime[NDOF2*numnodes*1+NDOF2*i+0] = 0.;
-		Bprime[NDOF2*numnodes*1+NDOF2*i+1] = dbasis[1*numnodes+i];
-		Bprime[NDOF2*numnodes*2+NDOF2*i+0] = dbasis[1*numnodes+i];
-		Bprime[NDOF2*numnodes*2+NDOF2*i+1] = dbasis[0*numnodes+i];
-		Bprime[NDOF2*numnodes*3+NDOF2*i+0] = dbasis[0*numnodes+i];
-		Bprime[NDOF2*numnodes*3+NDOF2*i+1] = dbasis[1*numnodes+i];
-	}
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(dbasis);
-}
-/*}}}*/
-/*FUNCTION TriaRef::GetBprimeMasstransport{{{*/
-void TriaRef::GetBprimeMasstransport(IssmDouble* Bprime, IssmDouble* xyz_list, Gauss* gauss){
-	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
-	 * For node i, Bi' can be expressed in the actual coordinate system
-	 * by: 
-	 *       Bi_prime=[ dN/dx ]
-	 *                [ dN/dy ]
-	 * where N is the finiteelement function for node i.
-	 *
-	 * We assume B' has been allocated already, of size: 3x(NDOF2*numnodes)
-	 */
-
-	/*Fetch number of nodes for this finite element*/
-	int numnodes = this->NumberofNodes();
-
-	/*Get nodal functions derivatives*/
-	IssmDouble* dbasis=xNew<IssmDouble>(2*numnodes);
-	GetNodalFunctionsDerivatives(dbasis,xyz_list,gauss);
-
-	/*Build B': */
-	for(int i=0;i<numnodes;i++){
-		Bprime[numnodes*0+i] = dbasis[0*numnodes+i];
-		Bprime[numnodes*1+i] = dbasis[1*numnodes+i];
-	}
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(dbasis);
-}
-/*}}}*/
-/*FUNCTION TriaRef::GetBSSAFriction{{{*/
-void TriaRef::GetBSSAFriction(IssmDouble* B, IssmDouble* xyz_list,Gauss* gauss){
-	/*Compute B  matrix. B=[B1 B2 B3] where Bi is square and of size 2. 
-	 * For node i, Bi can be expressed in the actual coordinate system
-	 * by: 
-	 *                 Bi=[ N   0 ]
-	 *                    [ 0   N ]
-	 * where N is the finiteelement function for node i.
-	 *
-	 * We assume B has been allocated already, of size: 2 x (numdof*numnodes)
-	 */
-
-	/*Fetch number of nodes for this finite element*/
-	int numnodes = this->NumberofNodes();
-
-	/*Get nodal functions derivatives*/
-	IssmDouble* basis=xNew<IssmDouble>(numnodes);
-	GetNodalFunctions(basis,gauss);
-
-	/*Build L: */
-	for(int i=0;i<numnodes;i++){
-		B[2*numnodes*0+2*i+0] = basis[i];
-		B[2*numnodes*0+2*i+1] = 0.;
-		B[2*numnodes*1+2*i+0] = 0.;
-		B[2*numnodes*1+2*i+1] = basis[i];
-	}
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(basis);
-}
-/*}}}*/
-/*FUNCTION TriaRef::GetLFS{{{*/
-void TriaRef::GetLFS(IssmDouble* LFS, Gauss* 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 ]
-	 * 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*NDOF2;
-
-	/*Get nodal functions derivatives*/
-	IssmDouble* vbasis=xNew<IssmDouble>(vnumnodes);
-	GetNodalFunctionsVelocity(vbasis,gauss);
-
-	/*Build LFS: */
-	for(int i=0;i<vnumnodes;i++){
-		LFS[2*i+0] = vbasis[i];
-		LFS[2*i+1] = 0.;
-	}
-
-	for(int i=0;i<pnumnodes;i++){
-		LFS[i+vnumdof+0] = 0.;
-	}
-
-	/*Clean-up*/
-	xDelete<IssmDouble>(vbasis);
 }
 /*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h	(revision 17094)
+++ /issm/trunk-jpl/src/c/classes/Elements/TriaRef.h	(revision 17095)
@@ -23,20 +23,8 @@
 
 		/*Numerics*/
-		void GetBExtrusion(IssmDouble* B_prime, IssmDouble* xyz_list, Gauss* gauss);
-		void GetBFS(IssmDouble* B_prime, IssmDouble* xyz_list, Gauss* gauss);
-		void GetBSSA(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss);
-		void GetBSSAFS(IssmDouble* B , IssmDouble* xyz_list, Gauss* gauss);
-		void GetBprimeFS(IssmDouble* B_prime, IssmDouble* xyz_list, Gauss* gauss);
-		void GetBprimeSSA(IssmDouble* Bprime, IssmDouble* xyz_list, Gauss* gauss);
-		void GetBprimeSSAFS(IssmDouble* Bprime, IssmDouble* xyz_list, Gauss* gauss);
-		void GetBprimeMasstransport(IssmDouble* Bprime, IssmDouble* xyz_list, Gauss* gauss);
-		void GetBMasstransport(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss);
-		void GetBHydro(IssmDouble* B, IssmDouble* xyz_list, Gauss* gauss);
-		void GetBSSAFriction(IssmDouble* L, IssmDouble* xyz_list,Gauss* gauss);
 		void GetJacobian(IssmDouble* J, IssmDouble* xyz_list,Gauss* gauss);
 		void GetSegmentJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss);
 		void GetJacobianDeterminant(IssmDouble* Jdet, IssmDouble* xyz_list,Gauss* gauss);
 		void GetJacobianInvert(IssmDouble*  Jinv, IssmDouble* xyz_list,Gauss* gauss);
-		void GetLFS(IssmDouble* LFS, Gauss* gauss);
 		void GetNodalFunctions(IssmDouble* basis,Gauss* gauss);
 		void GetNodalFunctions(IssmDouble* basis,Gauss* gauss,int finiteelement);