Index: /issm/trunk/src/c/objects/Elements/Beam.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Beam.cpp	(revision 4911)
+++ /issm/trunk/src/c/objects/Elements/Beam.cpp	(revision 4912)
@@ -249,22 +249,4 @@
 double Beam::CostFunction(void){
 	ISSMERROR(" not supported yet!");
-}
-/*}}}*/
-/*FUNCTION Beam::CreateKMatrix{{{1*/
-void  Beam::CreateKMatrix(Mat Kgg){
-
-	int analysis_type;
-
-	/*retrive parameters: */
-	parameters->FindParam(&analysis_type,AnalysisTypeEnum);
-
-	/*Just branch to the correct element stiffness matrix generator, according to the type of analysis we are carrying out: */
-	if (analysis_type==DiagnosticHutterAnalysisEnum) {
-		CreateKMatrixDiagnosticHutter( Kgg);
-	}
-	else{
-		ISSMERROR("analysis %i (%s) not supported yet",analysis_type,EnumAsString(analysis_type));
-	}
-
 }
 /*}}}*/
@@ -840,64 +822,4 @@
 
 /*Beam specific routines: */
-/*FUNCTION Beam::CreateKMatrixDiagnosticHutter{{{1*/
-
-void  Beam::CreateKMatrixDiagnosticHutter(Mat Kgg){
-	
-	
-	const int numgrids=2;
-	const int NDOF2=2;
-	const int numdofs=NDOF2*numgrids;
-	int       doflist[numdofs];
-	double    Ke_gg[numdofs][numdofs]={{0,0,0,0},{0,0,0,0},{0,0,0,0},{0,0,0,0}};
-	int       numberofdofspernode;
-	bool onbed;
-	bool onsurface;
-	int connectivity[2];
-	double one0,one1;
-	
-	connectivity[0]=nodes[0]->GetConnectivity();
-	connectivity[1]=nodes[1]->GetConnectivity();
-	
-	one0=1/(double)connectivity[0];
-	one1=1/(double)connectivity[1];
-	inputs->GetParameterValue(&onbed,ElementOnBedEnum);
-	inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
-
-	GetDofList(&doflist[0],&numberofdofspernode);
-
-	if (onbed && onsurface){
-		Ke_gg[0][0]=one0;
-		Ke_gg[1][1]=one0;
-		Ke_gg[2][0]=-one1;
-		Ke_gg[2][2]=one1;
-		Ke_gg[3][1]=-one1;
-		Ke_gg[3][3]=one1;
-	}
-	else if (onbed){
-		Ke_gg[0][0]=one0;
-		Ke_gg[1][1]=one0;
-		Ke_gg[2][0]=-2*one1;
-		Ke_gg[2][2]=2*one1;
-		Ke_gg[3][1]=-2*one1;
-		Ke_gg[3][3]=2*one1;
-	}
-	else if (onsurface){
-		Ke_gg[2][0]=-one1;
-		Ke_gg[2][2]=one1;
-		Ke_gg[3][1]=-one1;
-		Ke_gg[3][3]=one1;
-	}
-	else{ //node is on two horizontal layers and beams include the values only once, so the have to use half of the connectivity
-		Ke_gg[2][0]=-2*one1;
-		Ke_gg[2][2]=2*one1;
-		Ke_gg[3][1]=-2*one1;
-		Ke_gg[3][3]=2*one1;
-	}
-
-	/*Add Ke_gg to global matrix Kgg: */
-	MatSetValues(Kgg,numdofs,doflist,numdofs,doflist,(const double*)Ke_gg,ADD_VALUES);
-
-}
-/*}}}*/
 /*FUNCTION Beam::CreatePVectorDiagnosticHutter{{{1*/
 
Index: /issm/trunk/src/c/objects/Elements/Beam.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Beam.h	(revision 4911)
+++ /issm/trunk/src/c/objects/Elements/Beam.h	(revision 4912)
@@ -69,5 +69,5 @@
 		void	   SetCurrentConfiguration(Elements* elements,Loads* loads,DataSet* nodes,Materials* materials,Parameters* parameters);
 		double	CostFunction(void);
-		void	   CreateKMatrix(Mat Kgg);
+		void	   CreateKMatrix(Mat Kgg){ISSMERROR("Not implemented yet");};
 		void	   CreatePVector(Vec pg);
 		void     DeleteResults(void);
@@ -112,5 +112,4 @@
 		/*}}}*/
 		/*Beam specific routines: {{{1*/
-		void	  CreateKMatrixDiagnosticHutter(Mat Kgg);
 		void	  CreatePVectorDiagnosticHutter(Vec pg);
 		void	  GetDofList(int* doflist,int* pnumberofdofs);
Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 4911)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 4912)
@@ -2240,24 +2240,79 @@
 
 	/*Collapsed formulation: */
-	Beam*  beam=NULL;
-	int    i;
-
+	int       i;
+	const int numgrids=6;
+	const int NDOF2=2;
+	const int numdofs=NDOF2*numgrids;
+	int       doflist[numdofs];
+	double    Ke_gg[numdofs][numdofs]={0.0};
+	int       numberofdofspernode;
+	bool      onbed;
+	bool      onsurface;
+	int       connectivity[2];
+	double    one0,one1;
+	int       i0,i1,j0,j1;
+	
 	/*flags: */
 	bool onwater;
-
+	
 	/*recover some inputs: */
 	inputs->GetParameterValue(&onwater,ElementOnWaterEnum);
-
+	
 	/*If on water, skip: */
 	if(onwater)return;
 
+	inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+	inputs->GetParameterValue(&onsurface,ElementOnSurfaceEnum);
+
+	GetDofList(&doflist[0],&numberofdofspernode);
+
 	/*Spawn 3 beam elements: */
 	for(i=0;i<3;i++){
-		beam=(Beam*)SpawnBeam(i,i+3); //[0 3], [1 4] and [2 5] are the four vertical edges of the Penta
-		beam->CreateKMatrix(Kgg);
-	}
-
-	/*clean up*/
-	delete beam;
+		/*2 dofs of first node*/
+		i0=2*i;
+		i1=2*i+1;
+		/*2 dofs of second node*/
+		j0=2*(i+3);
+		j1=2*(i+3)+1;
+
+		/*Find connectivity for the two nodes*/
+		connectivity[0]=nodes[i]->GetConnectivity();
+		connectivity[1]=nodes[i+3]->GetConnectivity();
+		one0=1/(double)connectivity[0];
+		one1=1/(double)connectivity[1];
+
+		/*Create matrix for these two nodes*/
+		if (onbed && onsurface){
+			Ke_gg[i0][i0]=one0;
+			Ke_gg[i1][i1]=one0;
+			Ke_gg[j0][i0]=-one1;
+			Ke_gg[j0][j0]=one1;
+			Ke_gg[j1][i1]=-one1;
+			Ke_gg[j1][j1]=one1;
+		}
+		else if (onbed){
+			Ke_gg[i0][i0]=one0;
+			Ke_gg[i1][i1]=one0;
+			Ke_gg[j0][i0]=-2*one1;
+			Ke_gg[j0][j0]=2*one1;
+			Ke_gg[j1][i1]=-2*one1;
+			Ke_gg[j1][j1]=2*one1;
+		}
+		else if (onsurface){
+			Ke_gg[j0][i0]=-one1;
+			Ke_gg[j0][j0]=one1;
+			Ke_gg[j1][i1]=-one1;
+			Ke_gg[j1][j1]=one1;
+		}
+		else{ //node is on two horizontal layers and beams include the values only once, so the have to use half of the connectivity
+			Ke_gg[j0][i0]=-2*one1;
+			Ke_gg[j0][j0]=2*one1;
+			Ke_gg[j1][i1]=-2*one1;
+			Ke_gg[j1][j1]=2*one1;
+		}
+	}
+
+	/*Add Ke_gg to global matrix Kgg: */
+	MatSetValues(Kgg,numdofs,doflist,numdofs,doflist,(const double*)Ke_gg,ADD_VALUES);
 
 }