Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5849)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5850)
@@ -2212,5 +2212,7 @@
 			break;
 		case MacAyealPattynApproximationEnum:
-			CreateKMatrixDiagnosticMacAyeal3d( Kgg);
+			Ke=CreateKMatrixDiagnosticMacAyeal3d();
+			if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL);
+			delete Ke;
 			Ke=CreateKMatrixDiagnosticPattyn();
 			if(Ke) Ke->AddToGlobal(Kgg,NULL,NULL);
@@ -2324,10 +2326,24 @@
 /*}}}*/
 /*FUNCTION Penta::CreateKMatrixDiagnosticMacAyeal3d{{{1*/
-void Penta::CreateKMatrixDiagnosticMacAyeal3d( Mat Kgg){
+ElementMatrix* Penta::CreateKMatrixDiagnosticMacAyeal3d(void){
+
+
+	/*compute all stiffness matrices for this element*/
+	ElementMatrix* Ke1=CreateKMatrixDiagnosticMacAyeal3dViscous();
+	ElementMatrix* Ke2=CreateKMatrixDiagnosticMacAyeal3dFriction();
+	ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2);
+
+	/*clean-up and return*/
+	delete Ke1;
+	delete Ke2;
+	return Ke;
+}
+/*}}}*/
+/*FUNCTION Penta::CreateKMatrixDiagnosticMacAyeal3dViscous{{{1*/
+ElementMatrix* Penta::CreateKMatrixDiagnosticMacAyeal3dViscous(void){
 
 	const int    numdof2d=2*NUMVERTICES2D;
 	int             i,j,ig;
 	double       xyz_list[NUMVERTICES][3];
-	int*         doflist=NULL;
 	GaussPenta *gauss=NULL;
 	GaussTria  *gauss_tria=NULL;
@@ -2342,5 +2358,4 @@
 	double DL[2][2]={0.0}; //for basal drag
 	double DL_scalar;
-	double Ke_gg[numdof2d][numdof2d]={0.0}; //local element stiffness matrix 
 	double Ke_gg_gaussian[numdof2d][numdof2d]; //stiffness matrix evaluated at the gaussian point.
 	double Jdet;
@@ -2354,13 +2369,11 @@
 	Penta* pentabase=NULL;
 
-	/*If on water, skip stiffness: */
-	if(IsOnWater())return;
-
 	/*Find penta on bed as this is a macayeal elements: */
 	pentabase=GetBasalElement();
 	tria=pentabase->SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
 
-	/* Get dof list: */
-	tria->GetDofList(&doflist,MacAyealApproximationEnum,GsetEnum);  //Pattyn dof list
+	/*Initialize Element matrix and return if necessary*/
+	if(IsOnWater()) return NULL;
+	ElementMatrix* Ke=tria->NewElementMatrix(MacAyealApproximationEnum);
 
 	/*Retrieve all inputs and parameters*/
@@ -2398,16 +2411,5 @@
 					&Ke_gg_gaussian[0][0],0);
 
-		for(i=0;i<numdof2d;i++) for(j=0;j<numdof2d;j++) Ke_gg[i][j]+=Ke_gg_gaussian[i][j];
-	}
-
-	MatSetValues(Kgg,numdof2d,doflist,numdof2d,doflist,(const double*)Ke_gg,ADD_VALUES);
-
-	//Deal with 2d friction at the bedrock interface
-	if(IsOnBed() && !IsOnShelf()){
-		/*Build a tria element using the 3 grids of the base of the penta. Then use 
-		 * the tria functionality to build a friction stiffness matrix on these 3
-		 * grids: */
-
-		tria->CreateKMatrixDiagnosticMacAyealFriction(Kgg,NULL,NULL);
+		for(i=0;i<numdof2d;i++) for(j=0;j<numdof2d;j++) Ke->values[i*numdof2d+j]+=Ke_gg_gaussian[i][j];
 	}
 
@@ -2417,5 +2419,22 @@
 	delete gauss_tria;
 	delete gauss;
-	xfree((void**)&doflist);
+	return Ke;
+}
+/*}}}*/
+/*FUNCTION Penta::CreateKMatrixDiagnosticMacAyeal3dFriction{{{1*/
+ElementMatrix* Penta::CreateKMatrixDiagnosticMacAyeal3dFriction(void){
+
+	/*Initialize Element matrix and return if necessary*/
+	if(IsOnWater() || IsOnShelf() || !IsOnBed()) return NULL;
+
+	/*Build a tria element using the 3 grids of the base of the penta. Then use 
+	 * the tria functionality to build a friction stiffness matrix on these 3
+	 * grids: */
+	Tria* tria=(Tria*)SpawnTria(0,1,2); //grids 0, 1 and 2 make the new tria.
+	ElementMatrix* Ke=tria->CreateKMatrixDiagnosticMacAyealFriction();
+	delete tria->matice; delete tria;
+
+	/*clean-up and return*/
+	return Ke;
 }
 /*}}}*/
Index: /issm/trunk/src/c/objects/Elements/Penta.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.h	(revision 5849)
+++ /issm/trunk/src/c/objects/Elements/Penta.h	(revision 5850)
@@ -128,5 +128,7 @@
 		void	  CreateKMatrixDiagnosticHutter( Mat Kgg);
 		ElementMatrix* CreateKMatrixDiagnosticMacAyeal2d(void);
-		void	  CreateKMatrixDiagnosticMacAyeal3d(Mat Kgg);
+		ElementMatrix* CreateKMatrixDiagnosticMacAyeal3d(void);
+		ElementMatrix* CreateKMatrixDiagnosticMacAyeal3dViscous(void);
+		ElementMatrix* CreateKMatrixDiagnosticMacAyeal3dFriction(void);
 		ElementMatrix* CreateKMatrixDiagnosticPattyn(void);
 		ElementMatrix* CreateKMatrixDiagnosticPattynViscous(void);