Index: /issm/trunk/src/c/objects/Elements/Penta.cpp
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 4831)
+++ /issm/trunk/src/c/objects/Elements/Penta.cpp	(revision 4832)
@@ -2125,4 +2125,121 @@
 /*FUNCTION Penta::InputUpdateFromSolutionDiagnosticHoriz {{{1*/
 void  Penta::InputUpdateFromSolutionDiagnosticHoriz(double* solution){
+
+	/*Inputs*/
+	bool    collapsed,onbed;
+
+	/*Recover inputs*/
+	inputs->GetParameterValue(&collapsed,CollapseEnum);
+	inputs->GetParameterValue(&onbed,ElementOnBedEnum);
+
+	/*MacAyeal, everything is done by the element on bed*/
+	if (collapsed){
+		if (!onbed){
+			return;
+		}
+		else{
+			InputUpdateFromSolutionDiagnosticMacAyeal(solution);
+			return;
+		}
+	}
+	else{
+		InputUpdateFromSolutionDiagnosticPattyn(solution);
+	}
+}
+
+/*}}}*/
+/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticMacAyeal {{{1*/
+void  Penta::InputUpdateFromSolutionDiagnosticMacAyeal(double* solution){
+
+	int i;
+
+	const int    numvertices=6;
+	const int    numdofpervertex=2;
+	const int    numdof=numdofpervertex*numvertices;
+
+	int          doflist[numdof];
+	double       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  *VzInput        = NULL;
+	double *VzPtr          = NULL;
+	Penta  *penta          = NULL;
+
+	/*OK, we are on bed. Now recover results*/
+
+	/*Get dof list: */
+	GetDofList(&doflist[0],&dummy);
+
+	/*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++){
+		values[i]=solution[doflist[i]];
+	}
+
+	/*Ok, we have vx and vy in values, fill in vx and vy arrays and extrude */
+	for(i=0;i<3;i++){
+		vx[i]  =values[i*numdofpervertex+0];
+		vy[i]  =values[i*numdofpervertex+1];
+		vx[i+3]=vx[i];
+		vy[i+3]=vy[i];
+	}
+
+	/*Get parameters fro pressure computation*/
+	rho_ice=matpar->GetRhoIce();
+	g=matpar->GetG();
+
+	/*Start looping over all elements above current element and update all inputs*/
+	penta=this;
+	for(;;){
+
+		/*Now Compute vel*/
+		VzInput=inputs->GetInput(VzEnum);
+		if (VzInput){
+			if (VzInput->Enum()!=PentaVertexInputEnum) ISSMERROR("Cannot compute Vel as Vz is of type %s",EnumAsString(VzInput->Enum()));
+			VzInput->GetValuesPtr(&VzPtr,&dummy);
+			for(i=0;i<numvertices;i++) vz[i]=VzPtr[i];
+		}
+		else{for(i=0;i<numvertices;i++) vz[i]=0.0;}
+		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);
+
+		/*Now compute pressure*/
+		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: */
+		penta->inputs->ChangeEnum(VxEnum,VxOldEnum);
+		penta->inputs->ChangeEnum(VyEnum,VyOldEnum);
+		penta->inputs->ChangeEnum(PressureEnum,PressureOldEnum);
+
+		/*Add vx and vy as inputs to the tria element: */
+		penta->inputs->AddInput(new PentaVertexInput(VxEnum,vx));
+		penta->inputs->AddInput(new PentaVertexInput(VyEnum,vy));
+		penta->inputs->AddInput(new PentaVertexInput(VelEnum,vel));
+		penta->inputs->AddInput(new PentaVertexInput(PressureEnum,pressure));
+
+		/*Stop if we have reached the surface*/
+		if (penta->IsOnSurface()) break;
+
+		/* get upper Penta*/
+		penta=penta->GetUpperElement(); ISSMASSERT(penta->Id()!=this->id);
+	}
+}
+/*}}}*/
+/*FUNCTION Penta::InputUpdateFromSolutionDiagnosticPattyn {{{1*/
+void  Penta::InputUpdateFromSolutionDiagnosticPattyn(double* solution){
 	
 	int i;
@@ -2145,6 +2262,6 @@
 	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*       VzInput=NULL;
-	double*      VzPtr=NULL;
+	Input  *VzInput        = NULL;
+	double *VzPtr          = NULL;
 
 	/*Get dof list: */
@@ -2186,5 +2303,5 @@
 	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]);
@@ -2199,5 +2316,5 @@
 	this->inputs->AddInput(new PentaVertexInput(VxEnum,vx));
 	this->inputs->AddInput(new PentaVertexInput(VyEnum,vy));
-	this->inputs->AddInput(new TriaVertexInput(VelEnum,vel));
+	this->inputs->AddInput(new PentaVertexInput(VelEnum,vel));
 	this->inputs->AddInput(new PentaVertexInput(PressureEnum,pressure));
 }
@@ -4916,5 +5033,5 @@
 
 	/*Are we on the base, not on the surface?:*/
-	if(onbed==1){
+	if(onbed){
 
 		/*OK, we are on bed. we will follow the steps:
Index: /issm/trunk/src/c/objects/Elements/Penta.h
===================================================================
--- /issm/trunk/src/c/objects/Elements/Penta.h	(revision 4831)
+++ /issm/trunk/src/c/objects/Elements/Penta.h	(revision 4832)
@@ -177,4 +177,6 @@
 		void    InputUpdateFromSolutionBalancedvelocities( double* solutiong);
 		void    InputUpdateFromSolutionDiagnosticHoriz( double* solutiong);
+		void    InputUpdateFromSolutionDiagnosticMacAyeal( double* solutiong);
+		void    InputUpdateFromSolutionDiagnosticPattyn( double* solutiong);
 		void    InputUpdateFromSolutionDiagnosticHutter( double* solutiong);
 		void    InputUpdateFromSolutionDiagnosticVert( double* solutiong);
Index: /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp
===================================================================
--- /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp	(revision 4831)
+++ /issm/trunk/src/c/solvers/solver_diagnostic_nonlinear.cpp	(revision 4832)
@@ -128,12 +128,4 @@
 	}
 	
-
-	/*extrude if we are in 3D: */
-	if (dim==3){
-		if(verbose)_printf_("%s\n","extruding velocity of collapsed elements...");
-		InputExtrudex( femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum,true); //true means that the velocity is only extruded if collapsed
-		InputExtrudex( femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum,true);
-	}
-
 	//more output might be needed, when running in control
 	if(pKff0){