Index: /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp
===================================================================
--- /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16907)
+++ /issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp	(revision 16908)
@@ -3096,5 +3096,5 @@
 
 	/*Constants*/
-	int numnodes = 2*element->GetNumberOfNodes();
+	int numnodes    = element->GetNumberOfNodes();
 	int numdofm     = 1 *numnodes; //*2/2
 	int numdofp     = 2 *numnodes;
@@ -3111,6 +3111,6 @@
 	IssmDouble* Ke_gg          = xNewZeroInit<IssmDouble>(numdofp*numdofm);
 	IssmDouble* Ke_gg_gaussian = xNew<IssmDouble>(numdofp*numdofm);
-	Node        *node_list[numnodes];
-	int         cs_list[numnodes];
+	Node        *node_list[2*numnodes];
+	int*         cs_list= xNew<int>(2*numnodes);
 
 	/*Find penta on bed as HO must be coupled to the dofs on the bed: */
@@ -3150,5 +3150,4 @@
 		element->JacobianDeterminant(&Jdet, xyz_list,gauss);
 		this->GetBSSAHO(B, element,xyz_list, gauss);
-		//basaltria->GetBprimeSSA(Bprime, xyz_list, gauss_tria); /FIXME
 		this->GetBSSAprime(Bprime, basaltria,xyz_list, gauss_tria); 
 		element->ViscosityHO(&viscosity,xyz_list,gauss,vx_input,vy_input);
@@ -3164,16 +3163,20 @@
 					Ke_gg_gaussian,0);
 
-		for( i=0; i<numdofp; i++) for(j=0;j<numdofm; j++) Ke_gg[i*numdofp+j]+=Ke_gg_gaussian[i*numdofp+j];
+		for( i=0; i<numdofp; i++) for(j=0;j<numdofm; j++) Ke_gg[i*numdofm+j]+=Ke_gg_gaussian[i*numdofm+j];
 	} 
-	for(i=0;i<numdofp;i++) for(j=0;j<numdofm;j++) Ke->values[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i*numdofp+j];
-	for(i=0;i<numdofm;i++) for(j=0;j<numdofp;j++) Ke->values[i*numdoftotal+(j+2*numdofm)]+=Ke_gg[j*numdofp+i];
+	for(i=0;i<numdofp;i++) for(j=0;j<numdofm;j++) Ke->values[(i+2*numdofm)*numdoftotal+j]+=Ke_gg[i*numdofm+j];
+	for(i=0;i<numdofm;i++) for(j=0;j<numdofp;j++) Ke->values[i*numdoftotal+(j+2*numdofm)]+=Ke_gg[j*numdofm+i];
 
 	/*Transform Coordinate System*/
-	element->TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list);
+	element->TransformStiffnessMatrixCoord(Ke,node_list,2*numnodes,cs_list);
 
 	/*Clean-up and return*/
 	basaltria->DeleteMaterials(); delete basaltria;
+	
 	delete gauss;
 	delete gauss_tria;
+	xDelete<IssmDouble>(Ke_gg);
+	xDelete<IssmDouble>(Ke_gg_gaussian);
+	xDelete<int>(cs_list);
 	return Ke;
 
@@ -3211,10 +3214,10 @@
 	/*Build B: */
 	for(int i=0;i<numnodes;i++){
-		B[2*numnodes*0+2*i+0] = dbasis[0*3+i];
+		B[2*numnodes*0+2*i+0] = dbasis[0*numnodes+i];
 		B[2*numnodes*0+2*i+1] = 0.;
 		B[2*numnodes*1+2*i+0] = 0.;
-		B[2*numnodes*1+2*i+1] = dbasis[1*3+i];
-		B[2*numnodes*2+2*i+0] = .5*dbasis[1*3+i];
-		B[2*numnodes*2+2*i+1] = .5*dbasis[0*3+i];
+		B[2*numnodes*1+2*i+1] = dbasis[1*numnodes+i];
+		B[2*numnodes*2+2*i+0] = .5*dbasis[1*numnodes+i];
+		B[2*numnodes*2+2*i+1] = .5*dbasis[0*numnodes+i];
 	}
 
Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16907)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 16908)
@@ -1328,4 +1328,6 @@
 /*FUNCTION Penta::GetNode(int node_number) {{{*/
 Node* Penta::GetNode(int node_number){
+	_assert_(node_number>=0); 
+	_assert_(node_number<this->NumberofNodes()); 
 	return this->nodes[node_number];
 }