Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5372)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5373)
@@ -773,5 +773,4 @@
 		}
 		else if(approximation==MacAyealPattynApproximationEnum){
-			ISSMERROR("coupling MacAyeal/Pattyn not implemented yet");
 			CreateKMatrixDiagnosticMacAyeal( Kgg);
 			CreateKMatrixDiagnosticPattyn( Kgg);
@@ -2377,4 +2376,5 @@
 	/* local element matrices: */
 	double Ke_gg[numdofp][numdofm]={0.0}; //local element stiffness matrix 
+	double Ke_gg_transp[numdofm][numdofp]={0.0}; //local element stiffness matrix 
 	double Ke_gg_gaussian[numdofp][numdofm]; //stiffness matrix evaluated at the gaussian point.
 	double Jdet;
@@ -2490,8 +2490,8 @@
 	} //for (ig1=0; ig1<num_area_gauss; ig1++)
 
-	/*Add Ke_gg to global matrix Kgg: */
-	// one need to be transposed
+	/*Add Ke_gg and its transpose to global matrix Kgg: */
+	MatrixTranspose(&Ke_gg_transp[0][0],&Ke_gg[0][0],12,6);
 	MatSetValues(Kgg,numdofp,doflistp,numdofm,doflistm,(const double*)Ke_gg,ADD_VALUES);
-	MatSetValues(Kgg,numdofm,doflistm,numdofp,doflistp,(const double*)Ke_gg,ADD_VALUES);
+	MatSetValues(Kgg,numdofm,doflistm,numdofp,doflistp,(const double*)Ke_gg_transp,ADD_VALUES);
 
 	//Deal with 2d friction at the bedrock interface
@@ -2505,6 +2505,5 @@
 	}
 
-	delete tria;
-	delete pentabase;
+	delete tria->matice; delete tria;
 
 	xfree((void**)&first_gauss_area_coord);
@@ -4973,7 +4972,17 @@
 	double       vx;
 	double       vy;
-
-	/*Get dof list: */
-	GetDofList(&doflist);
+	int          approximation;
+
+	/*Get approximation enum and dof list: */
+	inputs->GetParameterValue(&approximation,ApproximationEnum);
+
+	/*If the element is a coupling, do nothing: every grid is also on an other elements 
+	 * (as coupling is between MacAyeal and Pattyn) so the other element will take care of it*/
+	if(approximation==MacAyealPattynApproximationEnum){
+	 return;
+	}
+	else{
+		GetDofList(&doflist,approximation);
+	}
 
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
@@ -5458,5 +5467,5 @@
 
 	/*Get dof list: */
-	GetDofList(&doflist);
+	GetDofList(&doflist,MacAyealApproximationEnum);
 
 	/*Use the dof list to index into the solution vector: */
@@ -5529,9 +5538,10 @@
 
 	const int    numvertices=6;
-	const int    numdofpervertex=4;
-	const int    solutionnumdof=2;
-	const int    numdof=solutionnumdof*numvertices;
-	int*         doflist=NULL;
-	int*         doflistbase=NULL;
+	const int    numvertices2d=3;
+	const int    numdofpervertex=2;
+	const int    numdof=numdofpervertex*numvertices;
+	const int    numdof2d=numdofpervertex*numvertices2d;
+	int*         doflistp=NULL;
+	int*         doflistm=NULL;
 	double       macayeal_values[numdof];
 	double       pattyn_values[numdof];
@@ -5556,6 +5566,6 @@
 
 	/*Get dof listof this element (pattyn dofs) and of the penta at base (macayeal dofs): */
-	GetDofList(&doflist);
-	penta->GetDofList(&doflistbase);
+	GetDofList(&doflistp,PattynApproximationEnum);
+	penta->GetDofList(&doflistm,MacAyealApproximationEnum);
 
 	/*Get node data: */
@@ -5563,7 +5573,11 @@
 
 	/*Use the dof list to index into the solution vector: */
-	for(i=0;i<numdof;i++){
-		pattyn_values[i]=solution[doflist[i]];
-		macayeal_values[i]=solution[doflistbase[i]];
+	for(i=0;i<numdof2d;i++){
+		pattyn_values[i]=solution[doflistp[i]];
+		macayeal_values[i]=solution[doflistm[i]];
+	}
+	for(i=numdof2d;i<numdof;i++){
+		pattyn_values[i]=solution[doflistp[i]];
+		macayeal_values[i]=macayeal_values[i-numdof2d];
 	}
 
@@ -5588,5 +5602,7 @@
 
 	/*Now Compute vel*/
-	for(i=0;i<numvertices;i++) vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
+	for(i=0;i<numvertices;i++) {
+		vel[i]=pow( pow(vx[i],2.0) + pow(vy[i],2.0) + pow(vz[i],2.0) , 0.5);
+	}
 
 	/*For pressure: we have not computed pressure in this analysis, for this element. We are in 3D, 
@@ -5612,6 +5628,6 @@
 
 	/*Free ressources:*/
-	xfree((void**)&doflist);
-	xfree((void**)&doflistbase);
+	xfree((void**)&doflistp);
+	xfree((void**)&doflistm);
 }
 /*}}}*/
@@ -5641,5 +5657,5 @@
 
 	/*Get dof list: */
-	GetDofList(&doflist);
+	GetDofList(&doflist,PattynApproximationEnum);
 
 	/*Get node data: */
Index: /issm/trunk/src/c/objects/Loads/Icefront.cpp
===================================================================
--- /issm/trunk/src/c/objects/Loads/Icefront.cpp	(revision 5372)
+++ /issm/trunk/src/c/objects/Loads/Icefront.cpp	(revision 5373)
@@ -448,5 +448,5 @@
 
 	/* Get dof list and node coordinates: */
-	GetDofList(&doflist);
+	GetDofList(&doflist,MacAyealApproximationEnum);
 	GetVerticesCoordinates(&xyz_list[0][0],nodes,numgrids);
 	
@@ -571,5 +571,5 @@
 
 	/* Get dof list and node coordinates: */
-	GetDofList(&doflist);
+	GetDofList(&doflist,PattynApproximationEnum);
 	GetVerticesCoordinates(&xyz_list[0][0], nodes, numgrids);
 	
@@ -641,5 +641,4 @@
 	normal_norm=norm(normal3);normal3[0]=normal3[0]/normal_norm;normal3[1]=normal3[1]/normal_norm;normal3[2]=normal3[2]/normal_norm;
 	normal_norm=norm(normal4);normal4[0]=normal4[0]/normal_norm;normal4[1]=normal4[1]/normal_norm;normal4[2]=normal4[2]/normal_norm;
-
 
 	//Compute load contribution for this quad:
@@ -796,5 +795,5 @@
 /*}}}*/
 /*FUNCTION Icefront::GetDofList {{{1*/
-void  Icefront::GetDofList(int** pdoflist){
+void  Icefront::GetDofList(int** pdoflist,int approximation_enum){
 
 	int i,j;
@@ -825,5 +824,5 @@
 	/*Figure out size of doflist: */
 	for(i=0;i<numberofnodes;i++){
-		numberofdofs+=nodes[i]->GetNumberOfDofs();
+		numberofdofs+=nodes[i]->GetNumberOfDofs(approximation_enum);
 	}
 
@@ -834,6 +833,6 @@
 	count=0;
 	for(i=0;i<numberofnodes;i++){
-		nodes[i]->GetDofList(doflist+count);
-		count+=nodes[i]->GetNumberOfDofs();
+		nodes[i]->GetDofList(doflist+count,approximation_enum);
+		count+=nodes[i]->GetNumberOfDofs(approximation_enum);
 	}
 
Index: /issm/trunk/src/c/objects/Loads/Icefront.h
===================================================================
--- /issm/trunk/src/c/objects/Loads/Icefront.h	(revision 5372)
+++ /issm/trunk/src/c/objects/Loads/Icefront.h	(revision 5373)
@@ -75,5 +75,5 @@
 		void  CreatePVectorDiagnosticHorizQuad( Vec pg);
 		void  CreatePVectorDiagnosticStokes( Vec pg);
-		void  GetDofList(int** pdoflist);
+		void  GetDofList(int** pdoflist,int approximation_enum=0);
 		void  QuadPressureLoad(double* pe_g,double rho_water,double rho_ice,double gravity, double* thickness_list, double* bed_list, 
 		                              double* normal1,double* normal2,double* normal3,double* normal4,double* xyz_list);
Index: /issm/trunk/src/c/objects/Materials/Matice.cpp
===================================================================
--- /issm/trunk/src/c/objects/Materials/Matice.cpp	(revision 5372)
+++ /issm/trunk/src/c/objects/Materials/Matice.cpp	(revision 5373)
@@ -378,5 +378,5 @@
 		if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) && 
 				(epsilon[3]==0) && (epsilon[4]==0)){
-			viscosity3d=pow((double)10,(double)14);
+			viscosity3d=0.5*pow((double)10,(double)14);
 		}
 		else{
@@ -450,5 +450,5 @@
 		if((epsilon[0]==0) && (epsilon[1]==0) && (epsilon[2]==0) && 
 				(epsilon[3]==0) && (epsilon[4]==0) && (epsilon[5]==0)){
-			viscosity3d=pow((double)10,(double)14);
+			viscosity3d=0.5*pow((double)10,(double)14);
 		}
 		else{