Index: /issm/trunk/src/c/objects/Elements/BeamRef.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/BeamRef.cpp	(revision 4941)
+++ /issm/trunk/src/c/objects/Elements/BeamRef.cpp	(revision 4942)
@@ -60,4 +60,13 @@
 }
 /*}}}*/
+/*FUNCTION BeamRef::GetJacobianDeterminant2d {{{1*/
+void BeamRef::GetJacobianDeterminant2d(double*  Jdet, double* xyz_list,double gauss){
+
+	/*Jdet = 0.5 * length*/
+	*Jdet=1.0/2.0*(sqrt(pow(xyz_list[3*1+0]-xyz_list[3*0+0],2)+pow(xyz_list[3*1+1]-xyz_list[3*0+1],2)));
+	if(*Jdet<0) ISSMERROR(" negative jacobian determinant!");
+
+}
+/*}}}*/
 /*FUNCTION BeamRef::GetNodalFunctions {{{1*/
 void BeamRef::GetNodalFunctions(double* l1l2, double gauss){
Index: /issm/trunk/src/c/objects/Elements/BeamRef.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/BeamRef.h	(revision 4941)
+++ /issm/trunk/src/c/objects/Elements/BeamRef.h	(revision 4942)
@@ -24,4 +24,5 @@
 		/*Numerics*/
 		void GetJacobianDeterminant(double* Jdet, double* z_list,double gauss);
+		void GetJacobianDeterminant2d(double* Jdet, double* xyz_list,double gauss);
 		void GetNodalFunctions(double* l1l2, double gauss);
 		void GetParameterValue(double* pvalue, double* value_list,double gauss_coord);
Index: /issm/trunk/src/c/objects/Elements/Element.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Element.h	(revision 4941)
+++ /issm/trunk/src/c/objects/Elements/Element.h	(revision 4942)
@@ -36,4 +36,6 @@
 		virtual bool   GetOnBed()=0;
 		virtual void   GetThicknessList(double* thickness_list)=0;
+		virtual void   GetParameterValue(double* pvalue,Node* node,int enumtype)=0;
+		virtual void   GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,Input* input_in)=0;
 		virtual void   GetBedList(double* bed_list)=0;
 		virtual void   Gradj(Vec gradient,int control_type)=0;
Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 4941)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 4942)
@@ -4148,5 +4148,4 @@
 		delete tria;
 	}
-	extern int my_rank;
 
 	xfree((void**)&first_gauss_area_coord);
@@ -4236,5 +4235,5 @@
 }
 /*}}}*/
-/*FUNCTION Penta::GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,int enumtype) {{{1*/
+/*FUNCTION Penta::GetParameterValue(double* pvalue,Node* node,int enumtype) {{{1*/
 void Penta::GetParameterValue(double* pvalue,Node* node,int enumtype){
 
@@ -4264,4 +4263,100 @@
 	return;
 
+}
+/*}}}*/
+/*FUNCTION Penta::GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,Input* input_in) {{{1*/
+void Penta::GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,Input* input_in){
+
+	/*Output*/
+	double value;
+
+	/*Intermediaries*/
+	const int numnodes=6;
+	int       grid1=-1,grid2=-1;
+	int       grid3;
+	int       i;
+	double    gauss_penta[numnodes];
+
+	/*go through 6 nodes (all nodes for penta) and identify 1st and 2nd nodes: */
+	ISSMASSERT(nodes);
+	for(i=0;i<numnodes;i++){
+		if (node1==nodes[i]) grid1=i;
+		if (node2==nodes[i]) grid2=i;
+	}
+
+	/*Reverse grid1 and 2 if necessary*/
+	if (grid1>grid2){
+
+		/*Reverse grid1 and grid2*/
+		grid3=grid1; grid1=grid2; grid2=grid3;
+
+		/*Change segment gauss point (between -1 and +1)*/
+		gauss_seg=-gauss_seg;
+	}
+
+	/*Build Penta Gauss point*/
+	if (grid2<3){
+
+		/*We are on the basal triangle*/
+		gauss_penta[3]=-1.;
+
+		if (grid1==0 && grid2==1){
+			/*gauss_seg is between 0 and 1*/
+			gauss_penta[0]=0.5*(1.-gauss_seg);
+			gauss_penta[1]=1.-0.5*(1.-gauss_seg);
+			gauss_penta[2]=0.;
+		}
+		else if (grid1==0 && grid2==2){
+			/*gauss_seg is between 0 and 2*/
+			gauss_penta[0]=0.5*(1.-gauss_seg);
+			gauss_penta[1]=0.;
+			gauss_penta[2]=1.-0.5*(1.-gauss_seg);
+		}
+		else if (grid1==1 && grid2==2){
+			/*gauss_seg is between 1 and 2*/
+			gauss_penta[0]=0.;
+			gauss_penta[1]=0.5*(1.-gauss_seg);
+			gauss_penta[2]=1.-0.5*(1.-gauss_seg);
+		}
+		else{
+			ISSMERROR("The 2 nodes provided are either the same or did not match current Penta nodes");
+		}
+	}
+	else if(grid1>2){
+
+		/*We are on the surface triangle*/
+		gauss_penta[3]=+1.;
+
+		if (grid1==3 && grid2==4){
+			/*gauss_seg is between 0 and 1*/
+			gauss_penta[0]=0.5*(1.-gauss_seg);
+			gauss_penta[1]=1.-0.5*(1.-gauss_seg);
+			gauss_penta[2]=0.;
+		}
+		else if (grid1==3 && grid2==5){
+			/*gauss_seg is between 0 and 2*/
+			gauss_penta[0]=0.5*(1.-gauss_seg);
+			gauss_penta[1]=0.;
+			gauss_penta[2]=1.-0.5*(1.-gauss_seg);
+		}
+		else if (grid1==4 && grid2==5){
+			/*gauss_seg is between 1 and 2*/
+			gauss_penta[0]=0.;
+			gauss_penta[1]=0.5*(1.-gauss_seg);
+			gauss_penta[2]=1.-0.5*(1.-gauss_seg);
+		}
+		else{
+			ISSMERROR("The 2 nodes provided are either the same or did not match current Penta nodes");
+		}
+	}
+	else{
+		ISSMERROR("vertical segments not implemented yet");
+	}
+
+	/*OK, now we can call input method*/
+	input_in->GetParameterValue(&value, &gauss_penta[0]);
+
+	/*Assign output pointers:*/
+	*pvalue=value;
 }
 /*}}}*/
Index: /issm/trunk/src/c/objects/Elements/Penta.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.h	(revision 4941)
+++ /issm/trunk/src/c/objects/Elements/Penta.h	(revision 4942)
@@ -147,4 +147,5 @@
 		void	  GetNodalFunctionsStokes(double* l1l7, double* gauss_coord);
 		void    GetParameterValue(double* pvalue,Node* node,int enumtype);
+		void    GetParameterValue(double* pvalue,Node* node1,Node* node2,double gauss_seg,Input* input_in);
 		void	  GetPhi(double* phi, double*  epsilon, double viscosity);
 		void	  GetSolutionFromInputsDiagnosticHoriz(Vec solutiong);
Index: /issm/trunk/src/c/objects/Elements/Tria.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.h	(revision 4941)
+++ /issm/trunk/src/c/objects/Elements/Tria.h	(revision 4942)
@@ -142,5 +142,4 @@
 		void	  GetDofList(int* doflist,int* pnumberofdofs);
 		void	  GetDofList1(int* doflist);
-		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);
Index: /issm/trunk/src/c/objects/Loads/Numericalflux.cpp
===================================================================
--- /issm/trunk/src/c/objects/Loads/Numericalflux.cpp	(revision 4941)
+++ /issm/trunk/src/c/objects/Loads/Numericalflux.cpp	(revision 4942)
@@ -880,5 +880,5 @@
 }
 /*}}}*/
-/*FUNCTION Numericalflux::GetParameterValue {{{1*/
+/*FUNCTION Numericalflux::GetParameterValue(double* pvalue, double gauss_coord,int enumtype) {{{1*/
 void Numericalflux::GetParameterValue(double* pvalue, double gauss_coord,int enumtype){
 
@@ -901,5 +901,5 @@
 }
 /*}}}*/
-/*FUNCTION Numericalflux::GetParameterValue {{{1*/
+/*FUNCTION Numericalflux::GetParameterValue(double* pvalue, double gauss_coord,Input* input_in) {{{1*/
 void Numericalflux::GetParameterValue(double* pvalue, double gauss_coord,Input* input_in){