Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15729)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15730)
@@ -954,5 +954,5 @@
 	/*First, figure out size of doflist and create it: */
 	int numberofdofs=0;
-	for(int i=0;i<numnodes;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(FSApproximationEnum,setenum);
+	for(int i=0;i<numnodes;i++) numberofdofs+=nodes[i]->GetNumberOfDofs(FSvelocityEnum,setenum);
 
 	/*Allocate output*/
@@ -962,6 +962,6 @@
 	int count=0;
 	for(int i=0;i<numnodes;i++){
-		nodes[i]->GetDofList(doflist+count,FSApproximationEnum,setenum);
-		count+=nodes[i]->GetNumberOfDofs(FSApproximationEnum,setenum);
+		nodes[i]->GetDofList(doflist+count,FSvelocityEnum,setenum);
+		count+=nodes[i]->GetNumberOfDofs(FSvelocityEnum,setenum);
 	}
 
@@ -7090,4 +7090,5 @@
 	int       cs_list[numnodes];
 	int       i,j;
+	int       init;
 
 	/*Prepare node list*/
@@ -7099,4 +7100,7 @@
 	}
 
+	/*Some parameters needed*/
+	init=this->element_type;
+
 	/*compute all stiffness matrices for this element*/
 	ElementMatrix* Ke1=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,HOApproximationEnum);
@@ -7105,6 +7109,9 @@
 	delete Ke1;
 	delete Ke2;
-	Ke1=CreateKMatrixDiagnosticHO(); TransformInvStiffnessMatrixCoord(Ke1,this->nodes,NUMVERTICES,XYEnum);
-	Ke2=CreateKMatrixDiagnosticFS(); TransformInvStiffnessMatrixCoord(Ke2,this->nodes,NUMVERTICES,XYZEnum);
+	
+	this->element_type=P1Enum;
+	Ke1=CreateKMatrixDiagnosticHO(); 
+	this->element_type=init;
+	Ke2=CreateKMatrixDiagnosticFS(); 
 
 	for(i=0;i<numdofs;i++) for(j=0;j<NUMVERTICES;j++){
@@ -7131,5 +7138,4 @@
 	int approximation;
 	inputs->GetInputValue(&approximation,ApproximationEnum);
-
 	switch(approximation){
 		case SSAApproximationEnum:
@@ -7664,5 +7670,8 @@
 
 	/*compute all stiffness matrices for this element*/
+	int init = this->element_type;
+	this->element_type=P1Enum;
 	ElementMatrix* Ke1=CreateKMatrixDiagnosticHO();
+	this->element_type=init;
 	ElementMatrix* Ke2=CreateKMatrixDiagnosticFS();
 	ElementMatrix* Ke3=CreateKMatrixCouplingHOFS();
@@ -9530,5 +9539,4 @@
 	for(int i=0;i<vnumnodes;i++){
 		gauss->GaussNode(this->VelocityInterpolation(),i);
-
 		vx_input->GetInputValue(&vx,gauss);
 		vy_input->GetInputValue(&vy,gauss);
@@ -9670,7 +9678,8 @@
 void  Penta::InputUpdateFromSolutionDiagnosticSSA(IssmDouble* solution){
 
-	const int    numdof=NDOF2*NUMVERTICES;
-
-	int     i;
+	int         numnodes = this->NumberofNodes();
+	int         numdof=NDOF2*numnodes;
+
+	int         i;
 	IssmDouble  rho_ice,g;
 	IssmDouble  values[numdof];
@@ -9861,5 +9870,5 @@
 	/*Get dof listof this element (SSA dofs) and of the penta at base (SSA dofs): */
 	penta->GetDofList(&doflistm,SSAApproximationEnum,GsetEnum);
-	GetDofList(&doflists,FSApproximationEnum,GsetEnum);
+	GetDofList(&doflists,FSvelocityEnum,GsetEnum);
 	this->parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
 
@@ -10090,9 +10099,11 @@
 
 	const int    numdofp=NDOF2*NUMVERTICES;
-	const int    numdofs=NDOF4*NUMVERTICES;
-
-	int    i;
+	const int    numdofs=NDOF3*NUMVERTICES;
+	const int    numdofpressure=NDOF1*NUMVERTICES;
+
+	int        i;
 	IssmDouble HO_values[numdofp];
 	IssmDouble FS_values[numdofs];
+	IssmDouble Pressure_values[numdofpressure];
 	IssmDouble vx[NUMVERTICES];
 	IssmDouble vy[NUMVERTICES];
@@ -10104,7 +10115,8 @@
 	IssmDouble xyz_list[NUMVERTICES][3];
 	IssmDouble FSreconditioning;
-	int*   doflistp      = NULL;
-	int*   doflists      = NULL;
-	Penta  *penta        = NULL;
+	int*       doflistp        = NULL;
+	int*       doflists        = NULL;
+	int*       doflistpressure = NULL;
+	Penta      *penta          = NULL;
 
 	/*OK, we have to add results of this element for HO 
@@ -10114,5 +10126,6 @@
 	/*Get dof listof this element (HO dofs) and of the penta at base (SSA dofs): */
 	GetDofList(&doflistp,HOApproximationEnum,GsetEnum);
-	GetDofList(&doflists,FSApproximationEnum,GsetEnum);
+	GetDofList(&doflists,FSvelocityEnum,GsetEnum);
+	GetDofListPressure(&doflistpressure,GsetEnum);
 	this->parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
 
@@ -10123,4 +10136,5 @@
 	for(i=0;i<numdofp;i++) HO_values[i]=solution[doflistp[i]];
 	for(i=0;i<numdofs;i++) FS_values[i]=solution[doflists[i]];
+	for(i=0;i<numdofpressure;i++) Pressure_values[i]=solution[doflistpressure[i]];
 
 	/*Transform solution in Cartesian Space*/
@@ -10130,8 +10144,8 @@
 	/*Ok, we have vx and vy in values, fill in vx and vy arrays: */
 	for(i=0;i<NUMVERTICES;i++){
-		vx[i]=FS_values[i*NDOF4+0]+HO_values[i*NDOF2+0];
-		vy[i]=FS_values[i*NDOF4+1]+HO_values[i*NDOF2+1];
-		vzFS[i]=FS_values[i*NDOF4+2];
-		pressure[i]=FS_values[i*NDOF4+3]*FSreconditioning;
+		vx[i]=FS_values[i*NDOF3+0]+HO_values[i*NDOF2+0];
+		vy[i]=FS_values[i*NDOF3+1]+HO_values[i*NDOF2+1];
+		vzFS[i]=FS_values[i*NDOF3+2];
+		pressure[i]=FS_values[i*NDOF1]*FSreconditioning;
 
 		/*Check solution*/
@@ -10178,4 +10192,5 @@
 	xDelete<int>(doflistp);
 	xDelete<int>(doflists);
+	xDelete<int>(doflistpressure);
 }
 /*}}}*/
@@ -10183,5 +10198,6 @@
 void  Penta::InputUpdateFromSolutionDiagnosticSIA(IssmDouble* solution){
 
-	const int    numdof=NDOF2*NUMVERTICES;
+	int         numnodes = this->NumberofNodes();
+	int         numdof=NDOF2*numnodes;
 
 	int     i;
@@ -10246,8 +10262,9 @@
 void  Penta::InputUpdateFromSolutionDiagnosticVert(IssmDouble* solution){
 
-	const int numdof=NDOF1*NUMVERTICES;
-
-	int      i;
-	int      approximation;
+	int          numnodes = this->NumberofNodes();
+	int          numdof=NDOF1*numnodes;
+
+	int          i;
+	int          approximation;
 	IssmDouble   rho_ice,g;
 	IssmDouble   values[numdof];
@@ -10262,5 +10279,5 @@
 	IssmDouble   surface[NUMVERTICES];
 	IssmDouble   xyz_list[NUMVERTICES][3];
-	int*     doflist      = NULL;
+	int*         doflist      = NULL;
 
 	/*Get the approximation and do nothing if the element in FS or None*/