Index: ../trunk-jpl/src/c/classes/Elements/PentaRef.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 15791)
+++ ../trunk-jpl/src/c/classes/Elements/PentaRef.cpp	(revision 15792)
@@ -2059,3 +2059,61 @@
 	return -1;
 }
 /*}}}*/
+/*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: ../trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15791)
+++ ../trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15792)
@@ -2973,6 +2973,8 @@
 void  Penta::ResetCoordinateSystem(void){
 
 	int        approximation;
+	int        numindices;
+	int       *indices = NULL;
 	IssmDouble slopex,slopey;
 	IssmDouble xz_plane[6];
 
@@ -2981,7 +2983,7 @@
 	if(IsFloating() || !IsOnBed() || (approximation!=FSApproximationEnum && approximation!=SSAFSApproximationEnum &&  approximation!=HOFSApproximationEnum)) return;
 
 	/*Get number of nodes for velocity only and base*/
-	int  vnumnodes = this->NumberofNodesVelocity();
+	BasalNodeIndices(&numindices,&indices,this->VelocityInterpolation());
 
 	/*Get inputs*/
 	Input* slopex_input=inputs->GetInput(BedSlopeXEnum); _assert_(slopex_input);
@@ -2989,8 +2991,8 @@
 
 	/*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);
 		slopey_input->GetInputValue(&slopey,gauss);
@@ -3000,10 +3002,11 @@
 		xz_plane[1]=0.;       xz_plane[4]=-slopey;  
 		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: ../trunk-jpl/src/c/classes/Elements/PentaRef.h
===================================================================
--- ../trunk-jpl/src/c/classes/Elements/PentaRef.h	(revision 15791)
+++ ../trunk-jpl/src/c/classes/Elements/PentaRef.h	(revision 15792)
@@ -68,6 +68,7 @@
 		void GetInputDerivativeValue(IssmDouble* pvalues, IssmDouble* plist,IssmDouble* xyz_list, GaussPenta* gauss);
 		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);
 		int  NumberofNodesPressure(void);