Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15624)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15625)
@@ -845,17 +845,17 @@
 void  Penta::GetDofList(int** pdoflist,int approximation_enum,int setenum){
 
-	int  i,count=0;
-	int  numberofdofs=0;
-	int* doflist=NULL;
-
-	/*First, figure out size of doflist: */
-	for(i=0;i<6;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
-
-	/*Allocate: */
-	doflist=xNew<int>(numberofdofs);
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = this->NumberofNodes();
+
+	/*First, figure out size of doflist and create it: */
+	int numberofdofs=0;
+	for(int i=0;i<numnodes;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
+
+	/*Allocate output*/
+	int* doflist=xNew<int>(numberofdofs);
 
 	/*Populate: */
-	count=0;
-	for(i=0;i<6;i++){
+	int count=0;
+	for(int i=0;i<numnodes;i++){
 		nodes[i]->GetDofList(doflist+count,approximation_enum,setenum);
 		count+=nodes[i]->GetNumberOfDofs(approximation_enum,setenum);
@@ -1102,4 +1102,45 @@
 		for (int iv=0;iv<NUMVERTICES;iv++) pvalue[iv]=defaultvalue;
 	}
+}
+/*}}}*/
+/*FUNCTION Penta::GetInputListOnNodes(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue) {{{*/
+void Penta::GetInputListOnNodes(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue){
+
+	_assert_(pvalue);
+
+	Input *input    = inputs->GetInput(enumtype);
+	int    numnodes = this->NumberofNodes();
+
+	/* Start looping on the number of vertices: */
+	if(input){
+		GaussPenta* gauss=new GaussPenta();
+		for(int iv=0;iv<this->NumberofNodes();iv++){
+			gauss->GaussNode(numnodes,iv);
+			input->GetInputValue(&pvalue[iv],gauss);
+		}
+		delete gauss;
+	}
+	else{
+		for(int iv=0;iv<numnodes;iv++) pvalue[iv]=defaultvalue;
+	}
+}
+/*}}}*/
+/*FUNCTION Penta::GetInputListOnNodes(IssmDouble* pvalue,int enumtype) {{{*/
+void Penta::GetInputListOnNodes(IssmDouble* pvalue,int enumtype){
+
+	_assert_(pvalue);
+
+	/*Recover input*/
+	Input* input=inputs->GetInput(enumtype);
+	if (!input) _error_("Input " << EnumToStringx(enumtype) << " not found in element");
+	int    numnodes = this->NumberofNodes();
+
+	/* Start looping on the number of vertices: */
+	GaussPenta* gauss=new GaussPenta();
+	for (int iv=0;iv<this->NumberofNodes();iv++){
+		gauss->GaussNode(numnodes,iv);
+		input->GetInputValue(&pvalue[iv],gauss);
+	}
+	delete gauss;
 }
 /*}}}*/
@@ -9407,23 +9448,22 @@
 void  Penta::InputUpdateFromSolutionDiagnosticHO(IssmDouble* solution){
 
-	const int    numdof=NDOF2*NUMVERTICES;
-
-	int    i;
-	IssmDouble rho_ice,g;
-	IssmDouble values[numdof];
-	IssmDouble vx[NUMVERTICES];
-	IssmDouble vy[NUMVERTICES];
-	IssmDouble vz[NUMVERTICES];
-	IssmDouble vel[NUMVERTICES];
-	IssmDouble pressure[NUMVERTICES];
-	IssmDouble surface[NUMVERTICES];
-	IssmDouble xyz_list[NUMVERTICES][3];
-	int*   doflist = NULL;
-
-	/*Get dof list: */
+	int         i;
+	IssmDouble  rho_ice,g;
+	IssmDouble  xyz_list[NUMVERTICES][3];
+	int        *doflist = NULL;
+
+	/*Fetch number of nodes and dof for this finite element*/
+	int numnodes = this->NumberofNodes();
+	int numdof   = numnodes*NDOF2;
+
+	/*Fetch dof list and allocate solution vectors*/
 	GetDofList(&doflist,HOApproximationEnum,GsetEnum);
-
-	/*Get node data: */
-	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
+	IssmDouble* values    = xNew<IssmDouble>(numdof);
+	IssmDouble* vx        = xNew<IssmDouble>(numnodes);
+	IssmDouble* vy        = xNew<IssmDouble>(numnodes);
+	IssmDouble* vz        = xNew<IssmDouble>(numnodes);
+	IssmDouble* vel       = xNew<IssmDouble>(numnodes);
+	IssmDouble* pressure  = xNew<IssmDouble>(numnodes);
+	IssmDouble* surface   = xNew<IssmDouble>(numnodes);
 
 	/*Use the dof list to index into the solution vector: */
@@ -9432,7 +9472,8 @@
 	/*Transform solution in Cartesian Space*/
 	TransformSolutionCoord(&values[0],nodes,NUMVERTICES,XYEnum);
+	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
-	for(i=0;i<NUMVERTICES;i++){
+	for(i=0;i<numnodes;i++){
 		vx[i]=values[i*NDOF2+0];
 		vy[i]=values[i*NDOF2+1];
@@ -9443,15 +9484,7 @@
 	}
 
-	/*Get Vz*/
-	Input* vz_input=inputs->GetInput(VzEnum);
-	if (vz_input){
-		GetInputListOnVertices(&vz[0],VzEnum);
-	}
-	else{
-		for(i=0;i<NUMVERTICES;i++) vz[i]=0.0;
-	}
-
-	/*Now Compute vel*/
-	for(i=0;i<NUMVERTICES;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
+	/*Get Vz and compute vel*/
+	GetInputListOnNodes(&vz[0],VzEnum,0.);
+	for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
 
 	/*For pressure: we have not computed pressure in this analysis, for this element. We are in 3D, 
@@ -9475,4 +9508,11 @@
 
 	/*Free ressources:*/
+	xDelete<IssmDouble>(surface);
+	xDelete<IssmDouble>(pressure);
+	xDelete<IssmDouble>(vel);
+	xDelete<IssmDouble>(vz);
+	xDelete<IssmDouble>(vy);
+	xDelete<IssmDouble>(vx);
+	xDelete<IssmDouble>(values);
 	xDelete<int>(doflist);
 }
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 15624)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h	(revision 15625)
@@ -192,4 +192,6 @@
 		void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype);
 		void           GetInputListOnVertices(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue);
+		void           GetInputListOnNodes(IssmDouble* pvalue,int enumtype);
+		void           GetInputListOnNodes(IssmDouble* pvalue,int enumtype,IssmDouble defaultvalue);
 		void           GetInputValue(IssmDouble* pvalue,Node* node,int enumtype);
 		void	         GetPhi(IssmDouble* phi, IssmDouble*  epsilon, IssmDouble viscosity);
Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15624)
+++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp	(revision 15625)
@@ -3425,10 +3425,10 @@
 	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
 	IssmDouble* values    = xNew<IssmDouble>(numdof);
-	IssmDouble* vx        = xNew<IssmDouble>(numdof);
-	IssmDouble* vy        = xNew<IssmDouble>(numdof);
-	IssmDouble* vz        = xNew<IssmDouble>(numdof);
-	IssmDouble* vel       = xNew<IssmDouble>(numdof);
-	IssmDouble* pressure  = xNew<IssmDouble>(numdof);
-	IssmDouble* thickness = xNew<IssmDouble>(numdof);
+	IssmDouble* vx        = xNew<IssmDouble>(numnodes);
+	IssmDouble* vy        = xNew<IssmDouble>(numnodes);
+	IssmDouble* vz        = xNew<IssmDouble>(numnodes);
+	IssmDouble* vel       = xNew<IssmDouble>(numnodes);
+	IssmDouble* pressure  = xNew<IssmDouble>(numnodes);
+	IssmDouble* thickness = xNew<IssmDouble>(numnodes);
 
 	/*Use the dof list to index into the solution vector: */
Index: /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.cpp	(revision 15624)
+++ /issm/trunk-jpl/src/c/classes/Inputs/PentaInput.cpp	(revision 15625)
@@ -141,5 +141,4 @@
 	IssmDouble B_reduced[6][DOFVELOCITY*numnodes];
 	IssmDouble velocity[numnodes][DOFVELOCITY];
-
 	_assert_(this->NumberofNodes()==6); //Check Tria too
 
@@ -313,5 +312,4 @@
 	/*Get B matrix: */
 	GetBHO(&B[0][0], xyz_list, gauss);
-
 	_assert_(this->NumberofNodes()==6); //Check Tria too
 
Index: /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp	(revision 15624)
+++ /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.cpp	(revision 15625)
@@ -366,24 +366,11 @@
 	/*update static arrays*/
 	switch(iv){
-		case 0:
-			coord1=1; coord2=0; coord3=0; coord4= -1;
-			break;
-		case 1:
-			coord1=0; coord2=1; coord3=0; coord4= -1;
-			break;
-		case 2:
-			coord1=0; coord2=0; coord3=1; coord4= -1;
-			break;
-		case 3:
-			coord1=1; coord2=0; coord3=0; coord4= +1;
-			break;
-		case 4:
-			coord1=0; coord2=1; coord3=0; coord4= +1;
-			break;
-		case 5:
-			coord1=0; coord2=0; coord3=1; coord4= +1;
-			break;
-		default:
-			_error_("vertex index should be in [0 5]");
+		case 0: coord1=1.; coord2=0.; coord3=0.; coord4= -1.; break;
+		case 1: coord1=0.; coord2=1.; coord3=0.; coord4= -1.; break;
+		case 2: coord1=0.; coord2=0.; coord3=1.; coord4= -1.; break;
+		case 3: coord1=1.; coord2=0.; coord3=0.; coord4= +1.; break;
+		case 4: coord1=0.; coord2=1.; coord3=0.; coord4= +1.; break;
+		case 5: coord1=0.; coord2=0.; coord3=1.; coord4= +1.; break;
+		default: _error_("vertex index should be in [0 5]");
 
 	}
@@ -405,4 +392,50 @@
 	else{
 		_error_("Tria not supported yet");
+	}
+
+}
+/*}}}*/
+/*FUNCTION GaussPenta::GaussNode{{{*/
+void GaussPenta::GaussNode(int numnodes,int iv){
+
+	/*in debugging mode: check that the default constructor has been called*/
+	_assert_(numgauss==-1);
+
+	/*update static arrays*/
+	switch(numnodes){
+		case 6: //P1 Lagrange element
+			switch(iv){
+				case 0: coord1=1.; coord2=0.; coord3=0.; coord4=-1.; break;
+				case 1: coord1=0.; coord2=1.; coord3=0.; coord4=-1.; break;
+				case 2: coord1=0.; coord2=0.; coord3=1.; coord4=-1.; break;
+				case 3: coord1=1.; coord2=0.; coord3=0.; coord4=+1.; break;
+				case 4: coord1=0.; coord2=1.; coord3=0.; coord4=+1.; break;
+				case 5: coord1=0.; coord2=0.; coord3=1.; coord4=+1.; break;
+				default: _error_("node index should be in [0 5]");
+			}
+			break;
+		case 15: //P2 Lagrange element
+			switch(iv){
+				case 0: coord1=1.; coord2=0.; coord3=0.; coord4=-1.; break;
+				case 1: coord1=0.; coord2=1.; coord3=0.; coord4=-1.; break;
+				case 2: coord1=0.; coord2=0.; coord3=1.; coord4=-1.; break;
+				case 3: coord1=1.; coord2=0.; coord3=0.; coord4=+1.; break;
+				case 4: coord1=0.; coord2=1.; coord3=0.; coord4=+1.; break;
+				case 5: coord1=0.; coord2=0.; coord3=1.; coord4=+1.; break;
+
+				case 6: coord1=1.; coord2=0.; coord3=0.; coord4=0.; break;
+				case 7: coord1=0.; coord2=1.; coord3=0.; coord4=0.; break;
+				case 8: coord1=0.; coord2=0.; coord3=1.; coord4=0.; break;
+
+				case  9: coord1=0.; coord2=.5; coord3=.5; coord4=-1.;break;
+				case 10: coord1=.5; coord2=0.; coord3=.5; coord4=-1.;break;
+				case 11: coord1=.5; coord2=.5; coord3=0.; coord4=-1.;break;
+				case 12: coord1=0.; coord2=.5; coord3=.5; coord4=+1.;break;
+				case 13: coord1=.5; coord2=0.; coord3=.5; coord4=+1.;break;
+				case 14: coord1=.5; coord2=.5; coord3=0.; coord4=+1.;break;
+				default: _error_("node index should be in [0 5]");
+			}
+			break;
+		default: _error_("supported number of nodes are 6 and 15");
 	}
 
Index: /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.h	(revision 15624)
+++ /issm/trunk-jpl/src/c/classes/gauss/GaussPenta.h	(revision 15625)
@@ -44,4 +44,5 @@
 		void GaussPoint(int ig);
 		void GaussVertex(int iv);
+		void GaussNode(int numnodes,int iv);
 		void GaussFaceTria(int index1, int index2, int index3, int order);
 		void GaussCenter(void);