Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 4769)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 4770)
@@ -4863,104 +4863,4 @@
 }
 /*}}}*/
-/*FUNCTION Tria::GetB_prog {{{1*/
-
-void Tria::GetB_prog(double* B_prog, double* xyz_list, double* gauss_l1l2l3){
-
-	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
-	 * For grid i, Bi can be expressed in the actual coordinate system
-	 * by: 
-	 *       Bi=[ h ]
-	 *                [ h ]
-	 * where h is the interpolation function for grid i.
-	 *
-	 * We assume B_prog has been allocated already, of size: 2x(NDOF1*numgrids)
-	 */
-
-	int i;
-	const int NDOF1=1;
-	const int numgrids=3;
-
-	double l1l2l3[numgrids];
-
-
-	/*Get dh1dh2dh3 in actual coordinate system: */
-	GetNodalFunctions(&l1l2l3[0],gauss_l1l2l3);
-
-	/*Build B_prog: */
-	for (i=0;i<numgrids;i++){
-		*(B_prog+NDOF1*numgrids*0+NDOF1*i)=l1l2l3[i];
-		*(B_prog+NDOF1*numgrids*1+NDOF1*i)=l1l2l3[i];
-	}
-}
-/*}}}*/
-/*FUNCTION Tria::GetBPrime {{{1*/
-
-void Tria::GetBPrime(double* Bprime, double* xyz_list, double* gauss_l1l2l3){
-
-	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
-	 * For grid i, Bi' can be expressed in the actual coordinate system
-	 * by: 
-	 *       Bi_prime=[ 2*dh/dx dh/dy ]
-	 *                       [ dh/dx  2*dh/dy]
-	 *                       [dh/dy dh/dx]
-	 * where h is the interpolation function for grid i.
-	 *
-	 * We assume B' has been allocated already, of size: 3x(NDOF2*numgrids)
-	 */
-	
-	int i;
-	const int NDOF2=2;
-	const int numgrids=3;
-
-	/*Same thing in the actual coordinate system: */
-	double dh1dh3[NDOF2][numgrids];
-
-
-	/*Get dh1dh2dh3 in actual coordinates system : */
-	GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list,gauss_l1l2l3);
-
-	/*Build B': */
-	for (i=0;i<numgrids;i++){
-		*(Bprime+NDOF2*numgrids*0+NDOF2*i)=2*dh1dh3[0][i]; 
-		*(Bprime+NDOF2*numgrids*0+NDOF2*i+1)=dh1dh3[1][i]; 
-		*(Bprime+NDOF2*numgrids*1+NDOF2*i)=dh1dh3[0][i]; 
-		*(Bprime+NDOF2*numgrids*1+NDOF2*i+1)=2*dh1dh3[1][i]; 
-		*(Bprime+NDOF2*numgrids*2+NDOF2*i)=dh1dh3[1][i]; 
-		*(Bprime+NDOF2*numgrids*2+NDOF2*i+1)=dh1dh3[0][i]; 
-	}
-}
-/*}}}*/
-/*FUNCTION Tria::GetBPrime_prog {{{1*/
-
-void Tria::GetBPrime_prog(double* Bprime_prog, double* xyz_list, double* gauss_l1l2l3){
-
-	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
-	 * For grid i, Bi' can be expressed in the actual coordinate system
-	 * by: 
-	 *       Bi_prime=[ dh/dx ]
-	 *                       [ dh/dy ]
-	 * where h is the interpolation function for grid i.
-	 *
-	 * We assume B' has been allocated already, of size: 3x(NDOF2*numgrids)
-	 */
-
-	int i;
-	const int NDOF1=1;
-	const int NDOF2=2;
-	const int numgrids=3;
-
-	/*Same thing in the actual coordinate system: */
-	double dh1dh3[NDOF2][numgrids];
-
-	/*Get dh1dh2dh3 in actual coordinates system : */
-	GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list,gauss_l1l2l3);
-
-	/*Build B': */
-	for (i=0;i<numgrids;i++){
-		*(Bprime_prog+NDOF1*numgrids*0+NDOF1*i)=dh1dh3[0][i]; 
-		*(Bprime_prog+NDOF1*numgrids*1+NDOF1*i)=dh1dh3[1][i]; 
-	}
-}
-/*}}}*/
 /*FUNCTION Tria::GetElementType {{{1*/
 int Tria::GetElementType(){
@@ -5004,53 +4904,4 @@
 	}
 
-}
-/*}}}*/
-/*FUNCTION Tria::GetL {{{1*/
-
-void Tria::GetL(double* L, double* xyz_list, double* gauss_l1l2l3,int numdof){
-
-	/*Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof. 
-	 * For grid 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 grid i.
-	 *
-	 * We assume L has been allocated already, of size: numgrids (numdof=1), or numdofx(numdof*numgrids) (numdof=2)
-	 */
-
-	int i;
-	const int NDOF2=2;
-	const int numgrids=3;
-
-	double l1l2l3[3];
-
-
-	/*Get l1l2l3 in actual coordinate system: */
-	GetNodalFunctions(l1l2l3, gauss_l1l2l3);
-
-#ifdef _DELUG_ 
-	for (i=0;i<3;i++){
-		printf("Node %i  h=%lf \n",i,l1l2l3[i]);
-	}
-#endif
-
-	/*Build L: */
-	if(numdof==1){
-		for (i=0;i<numgrids;i++){
-			L[i]=l1l2l3[i]; 
-		}
-	}
-	else{
-		for (i=0;i<numgrids;i++){
-			*(L+numdof*numgrids*0+numdof*i)=l1l2l3[i]; //L[0][NDOF2*i]=dh1dh3[0][i];
-			*(L+numdof*numgrids*0+numdof*i+1)=0;
-			*(L+numdof*numgrids*1+numdof*i)=0;
-			*(L+numdof*numgrids*1+numdof*i+1)=l1l2l3[i];
-		}
-	}
 }
 /*}}}*/
Index: /issm/trunk/src/c/objects/Elements/Tria.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.h	(revision 4769)
+++ /issm/trunk/src/c/objects/Elements/Tria.h	(revision 4770)
@@ -136,19 +136,15 @@
 		double  GetArea(void);
 		double  GetAreaCoordinate(double x, double y, int which_one);
-		void	  GetBPrime(double* Bprime, double* xyz_list, double* gauss_l1l2l3);
-		void	  GetBPrime_prog(double* Bprime_prog, double* xyz_list, double* gauss_l1l2l3);
-		void	  GetB_prog(double* B_prog, double* xyz_list, double* gauss_l1l2l3);
 		int     GetElementType(void);
 		void	  GetDofList(int* doflist,int* pnumberofdofs);
 		void	  GetDofList1(int* doflist);
-		void	  GetL(double* L, double* xyz_list, double* gauss_l1l2l3,int numdof);
 		void	  GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, double* gauss_l1l2l3);
 		void	  GetParameterValue(double* pp, double* plist, double* gauss_l1l2l3);
-		void      GetParameterValue(double* pvalue,Node* node,int enumtype);
-		void      GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,int enumtype);
+		void    GetParameterValue(double* pvalue,Node* node,int enumtype);
+		void    GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,int enumtype);
 		void	  GetSolutionFromInputsDiagnosticHoriz(Vec solution);
 		void	  GetSolutionFromInputsAdjointHoriz(Vec solution);
 		void	  GetSolutionFromInputsDiagnosticHutter(Vec solution);
-		void      GetStrainRate2d(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input);
+		void    GetStrainRate2d(double* epsilon,double* xyz_list, double* gauss, Input* vx_input, Input* vy_input);
 		void	  GradjDragStokes(Vec gradient);
 		void	  InputUpdateFromSolutionAdjointHoriz( double* solution);
Index: /issm/trunk/src/c/objects/Elements/TriaRef.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/TriaRef.cpp	(revision 4769)
+++ /issm/trunk/src/c/objects/Elements/TriaRef.cpp	(revision 4770)
@@ -83,4 +83,139 @@
 		*(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh3[1][i]; 
 		*(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh3[0][i]; 
+	}
+}
+/*}}}*/
+/*FUNCTION TriaRef::GetB_prog {{{1*/
+void TriaRef::GetB_prog(double* B_prog, double* xyz_list, double* gauss){
+	/*Compute B  matrix. B=[B1 B2 B3] where Bi is of size 3*NDOF2. 
+	 * For grid i, Bi can be expressed in the actual coordinate system
+	 * by: 
+	 *       Bi=[ h ]
+	 *          [ h ]
+	 * where h is the interpolation function for grid i.
+	 *
+	 * We assume B_prog has been allocated already, of size: 2x(NDOF1*numgrids)
+	 */
+
+	int i;
+	const int NDOF1=1;
+	const int numgrids=3;
+
+	double l1l2l3[numgrids];
+
+
+	/*Get dh1dh2dh3 in actual coordinate system: */
+	GetNodalFunctions(&l1l2l3[0],gauss);
+
+	/*Build B_prog: */
+	for (i=0;i<numgrids;i++){
+		*(B_prog+NDOF1*numgrids*0+NDOF1*i)=l1l2l3[i];
+		*(B_prog+NDOF1*numgrids*1+NDOF1*i)=l1l2l3[i];
+	}
+}
+/*}}}*/
+/*FUNCTION TriaRef::GetBPrime {{{1*/
+void TriaRef::GetBPrime(double* Bprime, double* xyz_list, double* gauss){
+
+	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
+	 * For grid i, Bi' can be expressed in the actual coordinate system
+	 * by: 
+	 *       Bi_prime=[ 2*dh/dx    dh/dy ]
+	 *                [   dh/dx  2*dh/dy ]
+	 *                [   dh/dy    dh/dx ]
+	 * where h is the interpolation function for grid i.
+	 *
+	 * We assume B' has been allocated already, of size: 3x(NDOF2*numgrids)
+	 */
+
+	int i;
+	const int NDOF2=2;
+	const int numgrids=3;
+
+	/*Same thing in the actual coordinate system: */
+	double dh1dh3[NDOF2][numgrids];
+
+	/*Get dh1dh2dh3 in actual coordinates system : */
+	GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list,gauss);
+
+	/*Build B': */
+	for (i=0;i<numgrids;i++){
+		*(Bprime+NDOF2*numgrids*0+NDOF2*i)=2*dh1dh3[0][i]; 
+		*(Bprime+NDOF2*numgrids*0+NDOF2*i+1)=dh1dh3[1][i]; 
+		*(Bprime+NDOF2*numgrids*1+NDOF2*i)=dh1dh3[0][i]; 
+		*(Bprime+NDOF2*numgrids*1+NDOF2*i+1)=2*dh1dh3[1][i]; 
+		*(Bprime+NDOF2*numgrids*2+NDOF2*i)=dh1dh3[1][i]; 
+		*(Bprime+NDOF2*numgrids*2+NDOF2*i+1)=dh1dh3[0][i]; 
+	}
+}
+/*}}}*/
+/*FUNCTION TriaRef::GetBPrime_prog {{{1*/
+void TriaRef::GetBPrime_prog(double* Bprime_prog, double* xyz_list, double* gauss){
+	/*Compute B'  matrix. B'=[B1' B2' B3'] where Bi' is of size 3*NDOF2. 
+	 * For grid i, Bi' can be expressed in the actual coordinate system
+	 * by: 
+	 *       Bi_prime=[ dh/dx ]
+	 *                       [ dh/dy ]
+	 * where h is the interpolation function for grid i.
+	 *
+	 * We assume B' has been allocated already, of size: 3x(NDOF2*numgrids)
+	 */
+
+	int i;
+	const int NDOF1=1;
+	const int NDOF2=2;
+	const int numgrids=3;
+
+	/*Same thing in the actual coordinate system: */
+	double dh1dh3[NDOF2][numgrids];
+
+	/*Get dh1dh2dh3 in actual coordinates system : */
+	GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list,gauss);
+
+	/*Build B': */
+	for (i=0;i<numgrids;i++){
+		*(Bprime_prog+NDOF1*numgrids*0+NDOF1*i)=dh1dh3[0][i]; 
+		*(Bprime_prog+NDOF1*numgrids*1+NDOF1*i)=dh1dh3[1][i]; 
+	}
+}
+/*}}}*/
+/*FUNCTION TriaRef::GetL {{{1*/
+void TriaRef::GetL(double* L, double* xyz_list, double* gauss,int numdof){
+	/*Compute L  matrix. L=[L1 L2 L3] where Li is square and of size numdof. 
+	 * For grid 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 grid i.
+	 *
+	 * We assume L has been allocated already, of size: numgrids (numdof=1), or numdofx(numdof*numgrids) (numdof=2)
+	 */
+
+	int i;
+	const int NDOF2=2;
+	const int numgrids=3;
+
+	double l1l2l3[3];
+
+
+	/*Get l1l2l3 in actual coordinate system: */
+	GetNodalFunctions(l1l2l3, gauss);
+
+	/*Build L: */
+	if(numdof==1){
+		for (i=0;i<numgrids;i++){
+			L[i]=l1l2l3[i]; 
+		}
+	}
+	else{
+		for (i=0;i<numgrids;i++){
+			*(L+numdof*numgrids*0+numdof*i)=l1l2l3[i]; //L[0][NDOF2*i]=dh1dh3[0][i];
+			*(L+numdof*numgrids*0+numdof*i+1)=0;
+			*(L+numdof*numgrids*1+numdof*i)=0;
+			*(L+numdof*numgrids*1+numdof*i+1)=l1l2l3[i];
+		}
 	}
 }
Index: /issm/trunk/src/c/objects/Elements/TriaRef.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/TriaRef.h	(revision 4769)
+++ /issm/trunk/src/c/objects/Elements/TriaRef.h	(revision 4770)
@@ -24,4 +24,8 @@
 		/*Numerics*/
 		void GetB(double* B, double* xyz_list, double* gauss);
+		void GetBPrime(double* Bprime, double* xyz_list, double* gauss);
+		void GetBPrime_prog(double* Bprime_prog, double* xyz_list, double* gauss);
+		void GetB_prog(double* B_prog, double* xyz_list, double* gauss);
+		void GetL(double* L, double* xyz_list, double* gauss,int numdof);
 		void GetJacobian(double* J, double* xyz_list,double* gauss);
 		void GetJacobianDeterminant2d(double* Jdet, double* xyz_list,double* gauss);