Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5210)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 5211)
@@ -3321,6 +3321,8 @@
 					D_scalar=D_scalar*dt;
 				}
-				K[0][0]=D_scalar*pow(u,2);       K[0][1]=D_scalar*fabs(u)*fabs(v);
-				K[1][0]=D_scalar*fabs(u)*fabs(v);K[1][1]=D_scalar*pow(v,2);
+				K[0][0]=D_scalar*pow(u,2);       K[0][1]=D_scalar*(u)*(v);
+				K[1][0]=D_scalar*(u)*(v);        K[1][1]=D_scalar*pow(v,2);
+				//K[0][0]=D_scalar*pow(u,2);       K[0][1]=D_scalar*fabs(u)*fabs(v);
+				//K[1][0]=D_scalar*fabs(u)*fabs(v);K[1][1]=D_scalar*pow(v,2);
 
 				/*Get B_artdiff: */
@@ -5191,6 +5193,9 @@
 		}
 	}
-	else{
+	else if (approximation==PattynApproximationEnum){
 		InputUpdateFromSolutionDiagnosticPattyn(solution);
+	}
+	else if (approximation==MacAyealPattynApproximationEnum){
+		InputUpdateFromSolutionDiagnosticMacAyealPattyn(solution);
 	}
 }
@@ -5288,4 +5293,97 @@
 	/*Free ressources:*/
 	xfree((void**)&doflist);
+}
+/*}}}*/
+/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticMacAyealPattyn {{{1*/
+void  Penta::InputUpdateFromSolutionDiagnosticMacAyealPattyn(double* solution){
+
+	int i;
+
+	const int    numvertices=6;
+	const int    numdofpervertex=4;
+	const int    solutionnumdof=2;
+	const int    numdof=solutionnumdof*numvertices;
+	int*         doflist=NULL;
+	int*         doflistbase=NULL;
+	double       macayeal_values[numdof];
+	double       pattyn_values[numdof];
+	double       vx[numvertices];
+	double       vy[numvertices];
+	double       vz[numvertices];
+	double       vel[numvertices];
+	int          dummy;
+	double       pressure[numvertices];
+	double       surface[numvertices];
+	double       rho_ice,g;
+	double       xyz_list[numvertices][3];
+	double       gauss[numvertices][4]={{1,0,0,-1},{0,1,0,-1},{0,0,1,-1},{1,0,0,1},{0,1,0,1},{0,0,1,1}};
+
+	Input  *vz_input       = NULL;
+	double *vz_ptr         = NULL;
+	Penta  *penta          = NULL;
+
+	/*OK, we have to add results of this element for pattyn 
+	 * and results from the penta at base for macayeal. Now recover results*/
+	penta=GetBasalElement();
+
+	/*Get dof listof this element (pattyn dofs) and of the penta at base (macayeal dofs): */
+	GetDofList(&doflist);
+	penta->GetDofList(&doflistbase);
+
+	/*Get node data: */
+	GetVerticesCoordinates(&xyz_list[0][0], nodes, numvertices);
+
+	/*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]];
+	}
+
+	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
+	for(i=0;i<numvertices;i++){
+		vx[i]=macayeal_values[i*numdofpervertex+0]+pattyn_values[i*numdofpervertex+0];
+		vy[i]=macayeal_values[i*numdofpervertex+1]+pattyn_values[i*numdofpervertex+1];
+	}
+
+	/*Get Vz*/
+	vz_input=inputs->GetInput(VzEnum);
+	if (vz_input){
+		if (vz_input->Enum()!=PentaVertexInputEnum){
+			ISSMERROR("Cannot compute Vel as Vz is of type %s",EnumToString(vz_input->Enum()));
+		}
+		vz_input->GetValuesPtr(&vz_ptr,&dummy);
+		for(i=0;i<numvertices;i++) vz[i]=vz_ptr[i];
+	}
+	else{
+		for(i=0;i<numvertices;i++) vz[i]=0.0;
+	}
+
+	/*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 pressure: we have not computed pressure in this analysis, for this element. We are in 3D, 
+	 *so the pressure is just the pressure at the z elevation: */
+	rho_ice=matpar->GetRhoIce();
+	g=matpar->GetG();
+	inputs->GetParameterValues(&surface[0],&gauss[0][0],6,SurfaceEnum);
+
+	for(i=0;i<numvertices;i++){
+		pressure[i]=rho_ice*g*(surface[i]-xyz_list[i][2]);
+	}
+	/*Now, we have to move the previous Vx and Vy inputs  to old 
+	 * status, otherwise, we'll wipe them off: */
+	this->inputs->ChangeEnum(VxEnum,VxOldEnum);
+	this->inputs->ChangeEnum(VyEnum,VyOldEnum);
+	this->inputs->ChangeEnum(PressureEnum,PressureOldEnum);
+
+	/*Add vx and vy as inputs to the tria element: */
+	this->inputs->AddInput(new PentaVertexInput(VxEnum,vx));
+	this->inputs->AddInput(new PentaVertexInput(VyEnum,vy));
+	this->inputs->AddInput(new PentaVertexInput(VelEnum,vel));
+	this->inputs->AddInput(new PentaVertexInput(PressureEnum,pressure));
+
+	/*Free ressources:*/
+	xfree((void**)&doflist);
+	xfree((void**)&doflistbase);
 }
 /*}}}*/
Index: /issm/trunk/src/c/objects/Elements/Penta.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.h	(revision 5210)
+++ /issm/trunk/src/c/objects/Elements/Penta.h	(revision 5211)
@@ -167,4 +167,5 @@
 		void    InputUpdateFromSolutionDiagnosticHoriz( double* solutiong);
 		void    InputUpdateFromSolutionDiagnosticMacAyeal( double* solutiong);
+		void    InputUpdateFromSolutionDiagnosticMacAyealPattyn( double* solutiong);
 		void    InputUpdateFromSolutionDiagnosticPattyn( double* solutiong);
 		void    InputUpdateFromSolutionDiagnosticHutter( double* solutiong);