Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5862)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5863)
@@ -702,5 +702,7 @@
 			break;
 		case DiagnosticHutterAnalysisEnum:
-			CreateKMatrixDiagnosticHutter( Kgg);
+			Ke=CreateKMatrixDiagnosticHutter();
+			if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL);
+			delete Ke;
 			break;
 		case DiagnosticVertAnalysisEnum:
@@ -2248,19 +2250,15 @@
 /*}}}*/
 /*FUNCTION Penta::CreateKMatrixDiagnosticHutter{{{1*/
-void  Penta::CreateKMatrixDiagnosticHutter(Mat Kgg){
-
-	/*Collapsed formulation: */
-	int       i;
-	const int numdofs=NDOF2*NUMVERTICES;
-	int*      doflist=NULL;
-	double    Ke_gg[numdofs][numdofs]={0.0};
+ElementMatrix* Penta::CreateKMatrixDiagnosticHutter(void){
+
+	/*Intermediaries*/
+	const int numdof=NDOF2*NUMVERTICES;
 	int       connectivity[2];
 	double    one0,one1;
-	int       i0,i1,j0,j1;
-	
-	/*If on water, skip: */
-	if(IsOnWater())return;
-
-	GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
+	int       i,i0,i1,j0,j1;
+
+	/*Initialize Element matrix and return if necessary*/
+	if(IsOnWater()) return NULL;
+	ElementMatrix* Ke=this->NewElementMatrix(NoneApproximationEnum);
 
 	/*Spawn 3 beam elements: */
@@ -2281,39 +2279,35 @@
 		/*Create matrix for these two nodes*/
 		if (IsOnBed() && IsOnSurface()){
-			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;
+			Ke->values[i0*numdof+i0]=one0;
+			Ke->values[i1*numdof+i1]=one0;
+			Ke->values[j0*numdof+i0]=-one1;
+			Ke->values[j0*numdof+j0]=one1;
+			Ke->values[j1*numdof+i1]=-one1;
+			Ke->values[j1*numdof+j1]=one1;
 		}
 		else if (IsOnBed()){
-			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;
+			Ke->values[i0*numdof+i0]=one0;
+			Ke->values[i1*numdof+i1]=one0;
+			Ke->values[j0*numdof+i0]=-2*one1;
+			Ke->values[j0*numdof+j0]=2*one1;
+			Ke->values[j1*numdof+i1]=-2*one1;
+			Ke->values[j1*numdof+j1]=2*one1;
 		}
 		else if (IsOnSurface()){
-			Ke_gg[j0][i0]=-one1;
-			Ke_gg[j0][j0]=one1;
-			Ke_gg[j1][i1]=-one1;
-			Ke_gg[j1][j1]=one1;
+			Ke->values[j0*numdof+i0]=-one1;
+			Ke->values[j0*numdof+j0]=one1;
+			Ke->values[j1*numdof+i1]=-one1;
+			Ke->values[j1*numdof+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);
-
-	/*Free ressources:*/
-	xfree((void**)&doflist);
-
+			Ke->values[j0*numdof+i0]=-2*one1;
+			Ke->values[j0*numdof+j0]=2*one1;
+			Ke->values[j1*numdof+i1]=-2*one1;
+			Ke->values[j1*numdof+j1]=2*one1;
+		}
+	}
+
+	/*Clean up and return*/
+	return Ke;
 }
 /*FUNCTION Penta::CreateKMatrixDiagnosticMacAyeal2d{{{1*/
Index: /issm/trunk/src/c/objects/Elements/Penta.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.h	(revision 5862)
+++ /issm/trunk/src/c/objects/Elements/Penta.h	(revision 5863)
@@ -128,5 +128,5 @@
 		void	  CreateKMatrixCouplingPattynStokes( Mat Kgg);
 		void	  CreateKMatrixDiagnosticHoriz( Mat Kgg);
-		void	  CreateKMatrixDiagnosticHutter( Mat Kgg);
+		ElementMatrix* CreateKMatrixDiagnosticHutter(void);
 		ElementMatrix* CreateKMatrixDiagnosticMacAyeal2d(void);
 		ElementMatrix* CreateKMatrixDiagnosticMacAyeal3d(void);