Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15791)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15792)
@@ -2974,4 +2974,6 @@
 
 	int        approximation;
+	int        numindices;
+	int       *indices = NULL;
 	IssmDouble slopex,slopey;
 	IssmDouble xz_plane[6];
@@ -2982,5 +2984,5 @@
 
 	/*Get number of nodes for velocity only and base*/
-	int  vnumnodes = this->NumberofNodesVelocity();
+	BasalNodeIndices(&numindices,&indices,this->VelocityInterpolation());
 
 	/*Get inputs*/
@@ -2990,6 +2992,6 @@
 	/*Loop over basal nodes and update their CS*/
 	GaussPenta* gauss = new GaussPenta();
-	for(int i=0;i<3;i++){//FIXME
-		gauss->GaussNode(this->VelocityInterpolation(),i);
+	for(int i=0;i<numindices;i++){//FIXME
+		gauss->GaussNode(this->VelocityInterpolation(),indices[i]);
 
 		slopex_input->GetInputValue(&slopex,gauss);
@@ -3001,8 +3003,9 @@
 		xz_plane[2]=slopex;   xz_plane[5]=1.;          
 
-		XZvectorsToCoordinateSystem(&this->nodes[i]->coord_system[0][0],&xz_plane[0]);
+		XZvectorsToCoordinateSystem(&this->nodes[indices[i]]->coord_system[0][0],&xz_plane[0]);
 	}
 
 	/*cleanup*/
+	xDelete<int>(indices);
 	delete gauss;
 }
Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 15791)
+++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 15792)
@@ -2060,2 +2060,60 @@
 }
 /*}}}*/
+/*FUNCTION PentaRef::BasalNodeIndices{{{*/
+void PentaRef::BasalNodeIndices(int* pnumindices,int** pindices,int finiteelement){
+
+	/*Output*/
+	int  numindices;
+	int* indices = NULL;
+
+	switch(finiteelement){
+		case P1Enum: case P1DGEnum:
+			numindices = 3;
+			indices    = xNew<int>(numindices);
+			indices[0] = 0;
+			indices[1] = 1;
+			indices[2] = 2;
+			break;
+		case P1bubbleEnum: case P1bubblecondensedEnum:
+			numindices = 3;
+			indices    = xNew<int>(numindices);
+			indices[0] = 0;
+			indices[1] = 1;
+			indices[2] = 2;
+			break;
+		case P2xP1Enum:
+			numindices = 6;
+			indices    = xNew<int>(numindices);
+			indices[0] = 0;
+			indices[1] = 1;
+			indices[2] = 2;
+			indices[3] = 6;
+			indices[4] = 7;
+			indices[5] = 8;
+			break;
+		case P1xP2Enum:
+			numindices = 3;
+			indices    = xNew<int>(numindices);
+			indices[0] = 0;
+			indices[1] = 1;
+			indices[2] = 2;
+			return;
+		case P2Enum:
+			numindices = 6;
+			indices    = xNew<int>(numindices);
+			indices[0] = 0;
+			indices[1] = 1;
+			indices[2] = 2;
+			indices[3] = 9;
+			indices[4] = 10;
+			indices[5] = 11;
+			break;
+		default:
+			_error_("Element type "<<EnumToStringx(this->element_type)<<" not supported yet");
+	}
+
+	/*Assign output pointer*/
+	*pnumindices = numindices;
+	*pindices    = indices;
+}
+/*}}}*/
Index: /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h	(revision 15791)
+++ /issm/trunk-jpl/src/c/classes/Elements/PentaRef.h	(revision 15792)
@@ -69,4 +69,5 @@
 		void GetInputDerivativeValue(IssmDouble* pvalues, IssmDouble* plist,IssmDouble* xyz_list, GaussTria* gauss){_error_("only PentaGauss are supported");};
 
+		int  BasalNodeIndices(int* pnumindices,int** pindices,int finiteelement);
 		int  NumberofNodes(void);
 		int  NumberofNodesVelocity(void);