Index: /issm/trunk/src/c/DataSet/Inputs.cpp
===================================================================
--- /issm/trunk/src/c/DataSet/Inputs.cpp	(revision 3799)
+++ /issm/trunk/src/c/DataSet/Inputs.cpp	(revision 3800)
@@ -300,5 +300,5 @@
 	}
 
-	if (!xinput | !yinput){
+	if (!xinput || !yinput){
 		/*we could not find one input with the correct enum type. No defaults values were provided, 
 		 * error out: */
Index: /issm/trunk/src/c/objects/Elements/Tria.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 3799)
+++ /issm/trunk/src/c/objects/Elements/Tria.cpp	(revision 3800)
@@ -1211,4 +1211,5 @@
 		/*Compute thickness at gaussian point: */
 		inputs->GetParameterValue(&thickness, gauss_l1l2l3,ThicknessEnum);
+		printf("thickness = %g\n",thickness);
 
 		/*Get strain rate from velocity: */
Index: /issm/trunk/src/c/objects/Inputs/TriaVertexInput.cpp
===================================================================
--- /issm/trunk/src/c/objects/Inputs/TriaVertexInput.cpp	(revision 3799)
+++ /issm/trunk/src/c/objects/Inputs/TriaVertexInput.cpp	(revision 3800)
@@ -149,5 +149,17 @@
 /*}}}*/
 /*FUNCTION TriaVertexInput::GetParameterValue(double* pvalue,double* gauss){{{1*/
-void TriaVertexInput::GetParameterValue(double* pvalue,double* gauss){ISSMERROR(" not supported yet!");}
+void TriaVertexInput::GetParameterValue(double* pvalue,double* gauss){
+	/*P1 interpolation on Gauss point*/
+
+	/*intermediary*/
+	double l1l2l3[3];
+
+	/*nodal functions: */
+	GetNodalFunctions(l1l2l3,gauss);
+
+	/*Assign output pointers:*/
+	*pvalue=l1l2l3[0]*values[0]+l1l2l3[1]*values[1]+l1l2l3[2]*values[2];
+
+}
 /*}}}*/
 /*FUNCTION TriaVertexInput::GetParameterValue(double* pvalue,double* gauss,double defaultvalue){{{1*/
@@ -161,5 +173,7 @@
 /*}}}*/
 /*FUNCTION TriaVertexInput::GetStrainRate(double* epsilon,Input* yinput, double* xyz_list, double* gauss){{{1*/
-void TriaVertexInput::GetStrainRate(double* epsilon,Input* yinput,double* xyz_list, double* gauss){ISSMERROR(" not supported yet!");}
+void TriaVertexInput::GetStrainRate(double* epsilon,Input* yinput,double* xyz_list, double* gauss){
+	ISSMERROR("not implemented yet");
+}
 /*}}}*/
 /*FUNCTION TriaVertexInput::GetStrainRateStokes(double* epsilon,Input* yinput, Input* zinput, double* xyz_list, double* gauss){{{1*/
@@ -171,2 +185,149 @@
 }
 /*}}}*/
+
+/*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);
+
+}
+/*}}}*/
Index: /issm/trunk/src/c/objects/Inputs/TriaVertexInput.h
===================================================================
--- /issm/trunk/src/c/objects/Inputs/TriaVertexInput.h	(revision 3799)
+++ /issm/trunk/src/c/objects/Inputs/TriaVertexInput.h	(revision 3800)
@@ -63,4 +63,11 @@
 		void GetStrainRateStokes(double* epsilon,Input* yinput, Input* zinput, double* xyz_list, double* gauss);
 		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);
 		/*}}}*/
 
Index: /issm/trunk/src/c/parallel/prognostic_core.cpp
===================================================================
--- /issm/trunk/src/c/parallel/prognostic_core.cpp	(revision 3799)
+++ /issm/trunk/src/c/parallel/prognostic_core.cpp	(revision 3800)
@@ -47,5 +47,6 @@
 	
 	_printf_("extract result from extruded inputs: \n");
-	InputToResult(&result,fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,ThicknessEnum,results->Size()+1,0,1);
+	result=new Result(results->Size()+1,0,1,"h_g",h_g);
+	//InputToResult(&result,fem_p->elements,fem_p->nodes, fem_p->vertices,fem_p->loads, fem_p->materials,fem_p->parameters,ThicknessEnum,results->Size()+1,0,1);
 	results->AddObject(result);