Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 4768)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 4769)
@@ -4863,40 +4863,4 @@
 }
 /*}}}*/
-/*FUNCTION Tria::GetB {{{1*/
-
-void Tria::GetB(double* B, 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=[ dh/dx    0    ]
-	 *                [   0    dh/dy  ]
-	 *                [ 1/2*dh/dy  1/2*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;
-
-	double dh1dh3[NDOF2][numgrids];
-
-
-	/*Get dh1dh2dh3 in actual coordinate system: */
-	GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list, gauss_l1l2l3);
-
-	/*Build B: */
-	for (i=0;i<numgrids;i++){
-		*(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh3[0][i]; //B[0][NDOF2*i]=dh1dh3[0][i];
-		*(B+NDOF2*numgrids*0+NDOF2*i+1)=0;
-		*(B+NDOF2*numgrids*1+NDOF2*i)=0;
-		*(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh3[1][i];
-		*(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh3[1][i]; 
-		*(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh3[0][i]; 
-	}
-}
-/*}}}*/
 /*FUNCTION Tria::GetB_prog {{{1*/
 
@@ -5042,93 +5006,4 @@
 }
 /*}}}*/
-/*FUNCTION Tria::GetJacobian {{{1*/
-void Tria::GetJacobian(double* J, double* xyz_list,double* gauss_l1l2l3){
-
-	/*The Jacobian is constant over the element, discard the gaussian points. 
-	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
-
-	const int NDOF2=2;
-	const int numgrids=3;
-	double x1,y1,x2,y2,x3,y3;
-	
-	x1=*(xyz_list+numgrids*0+0);
-	y1=*(xyz_list+numgrids*0+1);
-	x2=*(xyz_list+numgrids*1+0);
-	y2=*(xyz_list+numgrids*1+1);
-	x3=*(xyz_list+numgrids*2+0);
-	y3=*(xyz_list+numgrids*2+1);
-
-
-	*(J+NDOF2*0+0)=0.5*(x2-x1);
-	*(J+NDOF2*1+0)=SQRT3/6.0*(2*x3-x1-x2);
-	*(J+NDOF2*0+1)=0.5*(y2-y1);
-	*(J+NDOF2*1+1)=SQRT3/6.0*(2*y3-y1-y2);
-}
-/*}}}*/
-/*FUNCTION Tria::GetJacobianDeterminant2d {{{1*/
-void Tria::GetJacobianDeterminant2d(double*  Jdet, double* xyz_list,double* gauss_l1l2l3){
-
-	/*The Jacobian determinant is constant over the element, discard the gaussian points. 
-	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
-
-	double x1,x2,x3,y1,y2,y3;
-	
-	x1=*(xyz_list+3*0+0);
-	y1=*(xyz_list+3*0+1);
-	x2=*(xyz_list+3*1+0);
-	y2=*(xyz_list+3*1+1);
-	x3=*(xyz_list+3*2+0);
-	y3=*(xyz_list+3*2+1);
-
-	*Jdet=SQRT3/6.0*((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1));
-
-	if(Jdet<0){
-		ISSMERROR("negative jacobian determinant!");
-	}
-	
-}
-/*}}}*/
-/*FUNCTION Tria::GetJacobianDeterminant3d {{{1*/
-void Tria::GetJacobianDeterminant3d(double*  Jdet, double* xyz_list,double* gauss_l1l2l3){
-
-	/*The Jacobian determinant is constant over the element, discard the gaussian points. 
-	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
-
-	double 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=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){
-		ISSMERROR("negative jacobian determinant!");
-	}
-	
-}
-/*}}}*/
-/*FUNCTION Tria::GetJacobianInvert {{{1*/
-void Tria::GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss_l1l2l3){
-
-	double Jdet;
-	const int NDOF2=2;
-	const int numgrids=3;
-
-	/*Call Jacobian routine to get the jacobian:*/
-	GetJacobian(Jinv, xyz_list, gauss_l1l2l3);
-
-	/*Invert Jacobian matrix: */
-	MatrixInverse(Jinv,NDOF2,NDOF2,NULL,0,&Jdet);
-
-}
-/*}}}*/
 /*FUNCTION Tria::GetL {{{1*/
 
@@ -5178,76 +5053,4 @@
 		}
 	}
-}
-/*}}}*/
-/*FUNCTION Tria::GetNodalFunctions {{{1*/
-void Tria::GetNodalFunctions(double* l1l2l3, double* gauss_l1l2l3){
-	
-	/*This routine returns the values of the nodal functions  at the gaussian point.*/
-
-	/*First nodal function: */
-	l1l2l3[0]=gauss_l1l2l3[0];
-
-	/*Second nodal function: */
-	l1l2l3[1]=gauss_l1l2l3[1];
-
-	/*Third nodal function: */
-	l1l2l3[2]=gauss_l1l2l3[2];
-
-}
-/*}}}*/
-/*FUNCTION Tria::GetNodalFunctionsDerivatives {{{1*/
-void Tria::GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, double* gauss_l1l2l3){
-	
-	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
-	 * actual coordinate system: */
-
-	int i;
-	const int NDOF2=2;
-	const int numgrids=3;
-
-	double dh1dh3_ref[NDOF2][numgrids];
-	double Jinv[NDOF2][NDOF2];
-
-
-	/*Get derivative values with respect to parametric coordinate system: */
-	GetNodalFunctionsDerivativesReference(&dh1dh3_ref[0][0], gauss_l1l2l3); 
-
-	/*Get Jacobian invert: */
-	GetJacobianInvert(&Jinv[0][0], xyz_list, gauss_l1l2l3);
-
-	/*Build dh1dh3: 
-	 *
-	 * [dhi/dx]= Jinv*[dhi/dr]
-	 * [dhi/dy]       [dhi/ds]
-	 */
-
-	for (i=0;i<numgrids;i++){
-		*(dh1dh3+numgrids*0+i)=Jinv[0][0]*dh1dh3_ref[0][i]+Jinv[0][1]*dh1dh3_ref[1][i];
-		*(dh1dh3+numgrids*1+i)=Jinv[1][0]*dh1dh3_ref[0][i]+Jinv[1][1]*dh1dh3_ref[1][i];
-	}
-
-}
-/*}}}*/
-/*FUNCTION Tria::GetNodalFunctionsDerivativesReference {{{1*/
-void Tria::GetNodalFunctionsDerivativesReference(double* dl1dl3,double* gauss_l1l2l3){
-	
-	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
-	 * natural coordinate system) at the gaussian point. */
-
-	const int NDOF2=2;
-	const int numgrids=3;
-
-	/*First nodal function: */
-	*(dl1dl3+numgrids*0+0)=-0.5; 
-	*(dl1dl3+numgrids*1+0)=-1.0/(2.0*SQRT3);
-
-	/*Second nodal function: */
-	*(dl1dl3+numgrids*0+1)=0.5;
-	*(dl1dl3+numgrids*1+1)=-1.0/(2.0*SQRT3);
-
-	/*Third nodal function: */
-	*(dl1dl3+numgrids*0+2)=0;
-	*(dl1dl3+numgrids*1+2)=1.0/SQRT3;
-
 }
 /*}}}*/
Index: /issm/trunk/src/c/objects/Elements/Tria.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.h	(revision 4768)
+++ /issm/trunk/src/c/objects/Elements/Tria.h	(revision 4769)
@@ -136,5 +136,4 @@
 		double  GetArea(void);
 		double  GetAreaCoordinate(double x, double y, int which_one);
-		void	  GetB(double* B, double* xyz_list, double* gauss_l1l2l3);
 		void	  GetBPrime(double* Bprime, double* xyz_list, double* gauss_l1l2l3);
 		void	  GetBPrime_prog(double* Bprime_prog, double* xyz_list, double* gauss_l1l2l3);
@@ -143,12 +142,5 @@
 		void	  GetDofList(int* doflist,int* pnumberofdofs);
 		void	  GetDofList1(int* doflist);
-		void	  GetJacobian(double* J, double* xyz_list,double* gauss_l1l2l3);
-		void	  GetJacobianDeterminant2d(double*  Jdet, double* xyz_list,double* gauss_l1l2l3);
-		void	  GetJacobianDeterminant3d(double*  Jdet, double* xyz_list,double* gauss_l1l2l3);
-		void	  GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss_l1l2l3);
 		void	  GetL(double* L, double* xyz_list, double* gauss_l1l2l3,int numdof);
-		void	  GetNodalFunctions(double* l1l2l3, double* gauss_l1l2l3);
-		void	  GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, double* gauss_l1l2l3);
-		void	  GetNodalFunctionsDerivativesReference(double* dl1dl3,double* gauss_l1l2l3);
 		void	  GetParameterDerivativeValue(double* p, double* plist,double* xyz_list, double* gauss_l1l2l3);
 		void	  GetParameterValue(double* pp, double* plist, double* gauss_l1l2l3);
Index: /issm/trunk/src/c/objects/Elements/TriaRef.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/TriaRef.cpp	(revision 4768)
+++ /issm/trunk/src/c/objects/Elements/TriaRef.cpp	(revision 4769)
@@ -50,2 +50,193 @@
 }
 /*}}}*/
+
+/*Reference Element numerics*/
+/*FUNCTION TriaRef::GetB {{{1*/
+void TriaRef::GetB(double* B, 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=[ dh/dx           0    ]
+	 *          [   0           dh/dy  ]
+	 *          [ 1/2*dh/dy  1/2*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;
+
+	double dh1dh3[NDOF2][numgrids];
+
+
+	/*Get dh1dh2dh3 in actual coordinate system: */
+	GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list,gauss);
+
+	/*Build B: */
+	for (i=0;i<numgrids;i++){
+		*(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh3[0][i]; //B[0][NDOF2*i]=dh1dh3[0][i];
+		*(B+NDOF2*numgrids*0+NDOF2*i+1)=0;
+		*(B+NDOF2*numgrids*1+NDOF2*i)=0;
+		*(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh3[1][i];
+		*(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::GetJacobian {{{1*/
+void TriaRef::GetJacobian(double* J, double* xyz_list,double* gauss){
+	/*The Jacobian is constant over the element, discard the gaussian points. 
+	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
+
+	const int NDOF2=2;
+	const int numgrids=3;
+	double x1,y1,x2,y2,x3,y3;
+
+	x1=*(xyz_list+numgrids*0+0);
+	y1=*(xyz_list+numgrids*0+1);
+	x2=*(xyz_list+numgrids*1+0);
+	y2=*(xyz_list+numgrids*1+1);
+	x3=*(xyz_list+numgrids*2+0);
+	y3=*(xyz_list+numgrids*2+1);
+
+
+	*(J+NDOF2*0+0)=0.5*(x2-x1);
+	*(J+NDOF2*1+0)=SQRT3/6.0*(2*x3-x1-x2);
+	*(J+NDOF2*0+1)=0.5*(y2-y1);
+	*(J+NDOF2*1+1)=SQRT3/6.0*(2*y3-y1-y2);
+}
+/*}}}*/
+/*FUNCTION TriaRef::GetJacobianDeterminant2d {{{1*/
+void TriaRef::GetJacobianDeterminant2d(double* Jdet, double* xyz_list,double* gauss){
+
+	/*The Jacobian determinant is constant over the element, discard the gaussian points. 
+	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
+
+	double x1,x2,x3,y1,y2,y3;
+
+	x1=*(xyz_list+3*0+0);
+	y1=*(xyz_list+3*0+1);
+	x2=*(xyz_list+3*1+0);
+	y2=*(xyz_list+3*1+1);
+	x3=*(xyz_list+3*2+0);
+	y3=*(xyz_list+3*2+1);
+
+	*Jdet=SQRT3/6.0*((x2-x1)*(y3-y1)-(y2-y1)*(x3-x1));
+
+	if(Jdet<0){
+		ISSMERROR("negative jacobian determinant!");
+	}
+
+}
+/*}}}*/
+/*FUNCTION TriaRef::GetJacobianDeterminant3d {{{1*/
+void TriaRef::GetJacobianDeterminant3d(double*  Jdet, double* xyz_list,double* gauss){
+	/*The Jacobian determinant is constant over the element, discard the gaussian points. 
+	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
+
+	double 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=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){
+		ISSMERROR("negative jacobian determinant!");
+	}
+
+}
+/*}}}*/
+/*FUNCTION TriaRef::GetJacobianInvert {{{1*/
+void TriaRef::GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss){
+
+	double Jdet;
+	const int NDOF2=2;
+	const int numgrids=3;
+
+	/*Call Jacobian routine to get the jacobian:*/
+	GetJacobian(Jinv, xyz_list, gauss);
+
+	/*Invert Jacobian matrix: */
+	MatrixInverse(Jinv,NDOF2,NDOF2,NULL,0,&Jdet);
+
+}
+/*}}}*/
+/*FUNCTION TriaRef::GetNodalFunctions {{{1*/
+void TriaRef::GetNodalFunctions(double* l1l2l3, double* gauss){
+
+	/*This routine returns the values of the nodal functions  at the gaussian point.*/
+
+	/*First nodal function: */
+	l1l2l3[0]=gauss[0];
+
+	/*Second nodal function: */
+	l1l2l3[1]=gauss[1];
+
+	/*Third nodal function: */
+	l1l2l3[2]=gauss[2];
+
+}
+/*}}}*/
+/*FUNCTION TriaRef::GetNodalFunctionsDerivatives {{{1*/
+void TriaRef::GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, double* gauss){
+
+	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
+	 * actual coordinate system: */
+	int i;
+	const int NDOF2=2;
+	const int numgrids=3;
+	double dh1dh3_ref[NDOF2][numgrids];
+	double Jinv[NDOF2][NDOF2];
+
+	/*Get derivative values with respect to parametric coordinate system: */
+	GetNodalFunctionsDerivativesReference(&dh1dh3_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]
+	 */
+	for (i=0;i<numgrids;i++){
+		dh1dh3[numgrids*0+i]=Jinv[0][0]*dh1dh3_ref[0][i]+Jinv[0][1]*dh1dh3_ref[1][i];
+		dh1dh3[numgrids*1+i]=Jinv[1][0]*dh1dh3_ref[0][i]+Jinv[1][1]*dh1dh3_ref[1][i];
+	}
+
+}
+/*}}}*/
+/*FUNCTION TriaRef::GetNodalFunctionsDerivativesReference {{{1*/
+void TriaRef::GetNodalFunctionsDerivativesReference(double* dl1dl3,double* gauss){
+
+	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
+	 * natural coordinate system) at the gaussian point. */
+
+	const int NDOF2=2;
+	const int numgrids=3;
+
+	/*First nodal function: */
+	*(dl1dl3+numgrids*0+0)=-0.5; 
+	*(dl1dl3+numgrids*1+0)=-1.0/(2.0*SQRT3);
+
+	/*Second nodal function: */
+	*(dl1dl3+numgrids*0+1)=0.5;
+	*(dl1dl3+numgrids*1+1)=-1.0/(2.0*SQRT3);
+
+	/*Third nodal function: */
+	*(dl1dl3+numgrids*0+2)=0;
+	*(dl1dl3+numgrids*1+2)=1.0/SQRT3;
+
+}
+/*}}}*/
Index: /issm/trunk/src/c/objects/Elements/TriaRef.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/TriaRef.h	(revision 4768)
+++ /issm/trunk/src/c/objects/Elements/TriaRef.h	(revision 4769)
@@ -19,5 +19,16 @@
 		~TriaRef();
 
+		/*Management*/
 		void SetElementType(int type,int type_counter);
+
+		/*Numerics*/
+		void GetB(double* B, double* xyz_list, double* gauss);
+		void GetJacobian(double* J, double* xyz_list,double* gauss);
+		void GetJacobianDeterminant2d(double* Jdet, double* xyz_list,double* gauss);
+		void GetJacobianDeterminant3d(double* Jdet, double* xyz_list,double* gauss);
+		void GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss);
+		void GetNodalFunctions(double* l1l2l3, double* gauss);
+		void GetNodalFunctionsDerivatives(double* l1l2l3,double* xyz_list, double* gauss);
+		void GetNodalFunctionsDerivativesReference(double* dl1dl3,double* gauss);
 
 };
Index: /issm/trunk/src/c/objects/Inputs/TriaVertexInput.cpp
===================================================================
--- /issm/trunk/src/c/objects/Inputs/TriaVertexInput.cpp	(revision 4768)
+++ /issm/trunk/src/c/objects/Inputs/TriaVertexInput.cpp	(revision 4769)
@@ -24,7 +24,16 @@
 /*}}}*/
 /*FUNCTION TriaVertexInput::TriaVertexInput(int in_enum_type,double* values){{{1*/
-TriaVertexInput::TriaVertexInput(int in_enum_type,double* in_values){
-
+TriaVertexInput::TriaVertexInput(int in_enum_type,double* in_values)
+	:TriaRef(1)
+{
+
+	/*Set TriaRef*/
+	this->SetElementType(P1Enum,0);
+	this->element_type=P1Enum;
+
+	/*Set Enum*/
 	enum_type=in_enum_type;
+
+	/*Set values*/
 	values[0]=in_values[0];
 	values[1]=in_values[1];
@@ -311,145 +320,4 @@
 
 /*Intermediary*/
-/*FUNCTION TriaVertexInput::GetNodalFunctions {{{1*/
-void TriaVertexInput::GetNodalFunctions(double* l1l2l3, double* gauss){
-	
-	/*This routine returns the values of the nodal functions  at the gaussian point.*/
-
-	/*First nodal function: */
-	l1l2l3[0]=gauss[0];
-
-	/*Second nodal function: */
-	l1l2l3[1]=gauss[1];
-
-	/*Third nodal function: */
-	l1l2l3[2]=gauss[2];
-
-}
-/*}}}*/
-/*FUNCTION TriaVertexInput::GetB {{{1*/
-void TriaVertexInput::GetB(double* B, 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=[ dh/dx           0    ]
-	 *          [   0           dh/dy  ]
-	 *          [ 1/2*dh/dy  1/2*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;
-
-	double dh1dh3[NDOF2][numgrids];
-
-
-	/*Get dh1dh2dh3 in actual coordinate system: */
-	GetNodalFunctionsDerivatives(&dh1dh3[0][0],xyz_list,gauss);
-
-	/*Build B: */
-	for (i=0;i<numgrids;i++){
-		*(B+NDOF2*numgrids*0+NDOF2*i)=dh1dh3[0][i]; //B[0][NDOF2*i]=dh1dh3[0][i];
-		*(B+NDOF2*numgrids*0+NDOF2*i+1)=0;
-		*(B+NDOF2*numgrids*1+NDOF2*i)=0;
-		*(B+NDOF2*numgrids*1+NDOF2*i+1)=dh1dh3[1][i];
-		*(B+NDOF2*numgrids*2+NDOF2*i)=(float).5*dh1dh3[1][i]; 
-		*(B+NDOF2*numgrids*2+NDOF2*i+1)=(float).5*dh1dh3[0][i]; 
-	}
-}
-/*}}}*/
-/*FUNCTION TriaVertexInput::GetNodalFunctionsDerivatives {{{1*/
-void TriaVertexInput::GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, double* gauss){
-
-	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
-	 * actual coordinate system: */
-	int i;
-	const int NDOF2=2;
-	const int numgrids=3;
-	double dh1dh3_ref[NDOF2][numgrids];
-	double Jinv[NDOF2][NDOF2];
-
-	/*Get derivative values with respect to parametric coordinate system: */
-	GetNodalFunctionsDerivativesReference(&dh1dh3_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]
-	 */
-	for (i=0;i<numgrids;i++){
-		dh1dh3[numgrids*0+i]=Jinv[0][0]*dh1dh3_ref[0][i]+Jinv[0][1]*dh1dh3_ref[1][i];
-		dh1dh3[numgrids*1+i]=Jinv[1][0]*dh1dh3_ref[0][i]+Jinv[1][1]*dh1dh3_ref[1][i];
-	}
-
-}
-/*}}}*/
-/*FUNCTION TriaVertexInput::GetNodalFunctionsDerivativesReference {{{1*/
-void TriaVertexInput::GetNodalFunctionsDerivativesReference(double* dl1dl3,double* gauss_l1l2l3){
-
-	/*This routine returns the values of the nodal functions derivatives  (with respect to the 
-	 * natural coordinate system) at the gaussian point. */
-
-	const int NDOF2=2;
-	const int numgrids=3;
-
-	/*First nodal function: */
-	*(dl1dl3+numgrids*0+0)=-0.5; 
-	*(dl1dl3+numgrids*1+0)=-1.0/(2.0*SQRT3);
-
-	/*Second nodal function: */
-	*(dl1dl3+numgrids*0+1)=0.5;
-	*(dl1dl3+numgrids*1+1)=-1.0/(2.0*SQRT3);
-
-	/*Third nodal function: */
-	*(dl1dl3+numgrids*0+2)=0;
-	*(dl1dl3+numgrids*1+2)=1.0/SQRT3;
-
-}
-/*}}}*/
-/*FUNCTION TriaVertexInput::GetJacobian {{{1*/
-void TriaVertexInput::GetJacobian(double* J, double* xyz_list,double* gauss){
-
-	/*The Jacobian is constant over the element, discard the gaussian points. 
-	 * J is assumed to have been allocated of size NDOF2xNDOF2.*/
-
-	const int NDOF2=2;
-	const int numgrids=3;
-	double x1,y1,x2,y2,x3,y3;
-
-	x1=*(xyz_list+numgrids*0+0);
-	y1=*(xyz_list+numgrids*0+1);
-	x2=*(xyz_list+numgrids*1+0);
-	y2=*(xyz_list+numgrids*1+1);
-	x3=*(xyz_list+numgrids*2+0);
-	y3=*(xyz_list+numgrids*2+1);
-
-
-	*(J+NDOF2*0+0)=0.5*(x2-x1);
-	*(J+NDOF2*1+0)=SQRT3/6.0*(2*x3-x1-x2);
-	*(J+NDOF2*0+1)=0.5*(y2-y1);
-	*(J+NDOF2*1+1)=SQRT3/6.0*(2*y3-y1-y2);
-}
-/*}}}*/
-/*FUNCTION TriaVertexInput::GetJacobianInvert {{{1*/
-void TriaVertexInput::GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss){
-
-	double Jdet;
-	const int NDOF2=2;
-	const int numgrids=3;
-
-	/*Call Jacobian routine to get the jacobian:*/
-	GetJacobian(Jinv, xyz_list, gauss);
-
-	/*Invert Jacobian matrix: */
-	MatrixInverse(Jinv,NDOF2,NDOF2,NULL,0,&Jdet);
-
-}
-/*}}}*/
 /*FUNCTION TriaVertexInput::SquareMin{{{1*/
 void TriaVertexInput::SquareMin(double* psquaremin, bool process_units,Parameters* parameters){
Index: /issm/trunk/src/c/objects/Inputs/TriaVertexInput.h
===================================================================
--- /issm/trunk/src/c/objects/Inputs/TriaVertexInput.h	(revision 4768)
+++ /issm/trunk/src/c/objects/Inputs/TriaVertexInput.h	(revision 4769)
@@ -65,10 +65,4 @@
 		void ChangeEnum(int newenumtype);
 
-		void GetNodalFunctions(double* l1l2l3, double* gauss);
-		void GetB(double* B, double* xyz_list, double* gauss);
-		void GetNodalFunctionsDerivatives(double* dh1dh3,double* xyz_list, double* gauss);
-		void GetNodalFunctionsDerivativesReference(double* dl1dl3,double* gauss_l1l2l3);
-		void GetJacobian(double* J, double* xyz_list,double* gauss);
-		void GetJacobianInvert(double*  Jinv, double* xyz_list,double* gauss);
 		void SquareMin(double* psquaremin, bool process_units,Parameters* parameters);
 		void Scale(double scale_factor);