Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15769)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15770)
@@ -6990,5 +6990,4 @@
 
 	/*Constants*/
-	const int numnodes  = 2 *NUMVERTICES;
 	const int numdof    = NUMVERTICES *NDOF4;
 	const int numdofm   = NUMVERTICES *NDOF2;
@@ -7016,17 +7015,23 @@
 	Friction*  friction=NULL;
 	GaussPenta *gauss=NULL;
-	Node       *node_list[numnodes];
-	int         cs_list[numnodes];
 
 	/*If on water or not FS, skip stiffness: */
 	inputs->GetInputValue(&approximation,ApproximationEnum);
 	if(IsFloating() || !IsOnBed()) return NULL;
+
+	int vnumnodes = this->NumberofNodesVelocity();
+	int pnumnodes = this->NumberofNodesPressure();
+	int numnodes  = 2*vnumnodes-1+pnumnodes;
+
+	Node       *node_list[numnodes];
+	int         cs_list[numnodes];
+
 	ElementMatrix* Ke1=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,SSAApproximationEnum);
-	ElementMatrix* Ke2=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,FSvelocityEnum);
+	ElementMatrix* Ke2=new ElementMatrix(this->nodes,numnodes,this->parameters,FSvelocityEnum);
 	ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
 	delete Ke1; delete Ke2;
 
 	/*Prepare node list*/
-	for(i=0;i<NUMVERTICES;i++){
+	for(i=0;i<numnodes+NUMVERTICES;i++){
 		node_list[i+0*NUMVERTICES] = this->nodes[i];
 		node_list[i+1*NUMVERTICES] = this->nodes[i];
@@ -7097,4 +7102,5 @@
 
 	/*Clean up and return*/
+	xDelete<int>(cs_list);
 	delete gauss;
 	delete friction;
@@ -8352,5 +8358,15 @@
 	inputs->GetInputValue(&approximation,ApproximationEnum);
 	if(approximation!=HOFSApproximationEnum) return NULL;
-	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,FSvelocityEnum);
+
+	int vnumnodes = this->NumberofNodesVelocity();
+	int pnumnodes = this->NumberofNodesPressure();
+	int numnodes  = vnumnodes+pnumnodes;
+
+	/*Prepare coordinate system list*/
+	int* cs_list = xNew<int>(vnumnodes+pnumnodes);
+	for(i=0;i<vnumnodes;i++) cs_list[i]           = XYZEnum;
+	for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
+
+	ElementVector* pe=new ElementVector(nodes,numnodes,this->parameters,FSvelocityEnum);
 
 	/*Retrieve all inputs and parameters*/
@@ -8386,14 +8402,15 @@
 
 		for(i=0;i<NUMVERTICES2D;i++){
-			pe->values[i*NDOF4+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[0])*basis[i];
-			pe->values[i*NDOF4+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[1])*basis[i];
-			pe->values[i*NDOF4+2]+=Jdet2d*gauss->weight*2*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*basis[i];
+			pe->values[i*NDOF3+0]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[0]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[0])*basis[i];
+			pe->values[i*NDOF3+1]+=Jdet2d*gauss->weight*(alpha2_gauss*w*bed_normal[1]*bed_normal[2]+2*viscosity*dw[2]*bed_normal[1])*basis[i];
+			pe->values[i*NDOF3+2]+=Jdet2d*gauss->weight*2*viscosity*(dw[0]*bed_normal[0]+dw[1]*bed_normal[1]+dw[2]*bed_normal[2])*basis[i];
 		}
 	}
 
 	/*Transform coordinate system*/
-	TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZEnum);
+	TransformLoadVectorCoord(pe,nodes,vnumnodes+pnumnodes,cs_list);
 
 	/*Clean up and return*/
+	xDelete<int>(cs_list);
 	delete gauss;
 	delete friction;