Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp
===================================================================
--- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15743)
+++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp	(revision 15744)
@@ -401,17 +401,24 @@
 
 	if(Ke){
-		/*Condense if requested*/
-		if(this->element_type==MINIcondensedEnum){
-			int indices[3]={18,19,20};
-			Ke->StaticCondensation(3,&indices[0]);
-		}
-		else if(this->element_type==P1bubblecondensedEnum){
-			int size   = nodes[6]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
-			int offset = 0;
-			for(int i=0;i<6;i++) offset+=nodes[i]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
-			int* indices=xNew<int>(size);
-			for(int i=0;i<size;i++) indices[i] = offset+i;
-			Ke->StaticCondensation(size,indices);
-			xDelete<int>(indices);
+		int approximation;
+		inputs->GetInputValue(&approximation,ApproximationEnum);
+		if(approximation==HOFSApproximationEnum){
+			//Do nothing condensatino already done for Stokes part
+		}
+		else{
+			/*Condense if requested*/
+			if(this->element_type==MINIcondensedEnum){
+				int indices[3]={18,19,20};
+				Ke->StaticCondensation(3,&indices[0]);
+			}
+			else if(this->element_type==P1bubblecondensedEnum){
+				int size   = nodes[6]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
+				int offset = 0;
+				for(int i=0;i<6;i++) offset+=nodes[i]->GetNumberOfDofs(NoneApproximationEnum,GsetEnum);
+				int* indices=xNew<int>(size);
+				for(int i=0;i<size;i++) indices[i] = offset+i;
+				Ke->StaticCondensation(size,indices);
+				xDelete<int>(indices);
+			}
 		}
 
@@ -565,12 +572,20 @@
 		/*StaticCondensation if requested*/
 		if(this->element_type==MINIcondensedEnum){
-			int indices[3]={18,19,20};
-
-			this->element_type=MINIEnum;
-			ElementMatrix* Ke = CreateKMatrixDiagnosticFS();
-			this->element_type=MINIcondensedEnum;
-
-			pe->StaticCondensation(Ke,3,&indices[0]);
-			delete Ke;
+			int approximation;
+			inputs->GetInputValue(&approximation,ApproximationEnum);
+			if(approximation==HOFSApproximationEnum){
+				//Do nothing, condensation already done in PVectorCoupling
+			}
+			else{
+				int indices[3]={18,19,20};
+
+				this->element_type=MINIEnum;
+				ElementMatrix* Ke = CreateKMatrixDiagnosticFS();
+				this->element_type=MINIcondensedEnum;
+
+				pe->StaticCondensation(Ke,3,&indices[0]);
+				delete Ke;
+			}
+
 		}
 		else if(this->element_type==P1bubblecondensedEnum){
@@ -6624,5 +6639,5 @@
 	int pnumnodes = this->NumberofNodesPressure();
 
-	De=new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSApproximationEnum);
+	De=new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
 
 	for(int i=0;i<vnumnodes;i++){
@@ -7068,22 +7083,11 @@
 
 	/*Constants*/
-	const int numnodes  = 2 *NUMVERTICES;
-	const int numdofp     = NDOF2 *NUMVERTICES;
-	const int numdofs     = NDOF4 *NUMVERTICES;
-	const int numdoftotal = (NDOF2+NDOF4) *NUMVERTICES;
+	const int numnodes       = 2 *NUMVERTICES;
+	const int numdofp        = NDOF2 *NUMVERTICES;
+	const int numdofs        = NDOF4 * 6 + NDOF3;
+	const int numdoftotal    = (NDOF2+NDOF4) *NUMVERTICES + NDOF3;
 
 	/*Intermediaries*/
-	Node     *node_list[numnodes];
-	int       cs_list[numnodes];
-	int       i,j;
-	int       init;
-
-	/*Prepare node list*/
-	for(i=0;i<NUMVERTICES;i++){
-		node_list[i+0*NUMVERTICES] = this->nodes[i];
-		node_list[i+1*NUMVERTICES] = this->nodes[i];
-		cs_list[i+0*NUMVERTICES] = XYEnum;
-		cs_list[i+1*NUMVERTICES] = XYZEnum;
-	}
+	int       i,j,init;
 
 	/*Some parameters needed*/
@@ -7092,25 +7096,29 @@
 	/*compute all stiffness matrices for this element*/
 	ElementMatrix* Ke1=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,HOApproximationEnum);
-	ElementMatrix* Ke2=new ElementMatrix(this->nodes,NUMVERTICES,this->parameters,FSvelocityEnum);
+	ElementMatrix* Ke2=new ElementMatrix(this->nodes,2*NUMVERTICES+1,this->parameters,FSvelocityEnum);
 	ElementMatrix* Ke=new ElementMatrix(Ke1,Ke2);
 	delete Ke1;
 	delete Ke2;
 	
+	/*Compute HO Matrix with P1 element type\n");*/
 	this->element_type=P1Enum;
-	Ke1=CreateKMatrixDiagnosticHO(); 
+	Ke1=CreateKMatrixDiagnosticHO();
 	this->element_type=init;
+	/*Compute FS Matrix and condense it \n");*/
 	Ke2=CreateKMatrixDiagnosticFS(); 
+	int indices[3]={18,19,20};
+	Ke2->StaticCondensation(3,&indices[0]);
 
 	for(i=0;i<numdofs;i++) for(j=0;j<NUMVERTICES;j++){
-		Ke->values[(i+numdofp)*numdoftotal+NDOF2*j+0]+=Ke2->values[i*numdofs+NDOF4*j+0];
-		Ke->values[(i+numdofp)*numdoftotal+NDOF2*j+1]+=Ke2->values[i*numdofs+NDOF4*j+1];
+		Ke->values[(i+numdofp)*numdoftotal+NDOF2*j+0]+=Ke2->values[i*numdofs+NDOF3*j+0];
+		Ke->values[(i+numdofp)*numdoftotal+NDOF2*j+1]+=Ke2->values[i*numdofs+NDOF3*j+1];
 	}
 	for(i=0;i<numdofp;i++) for(j=0;j<NUMVERTICES;j++){
-		Ke->values[i*numdoftotal+numdofp+NDOF4*j+0]+=Ke1->values[i*numdofp+NDOF2*j+0];
-		Ke->values[i*numdoftotal+numdofp+NDOF4*j+1]+=Ke1->values[i*numdofp+NDOF2*j+1];
-	}
-
-	/*Transform Coordinate System*/
-	TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list);
+		Ke->values[i*numdoftotal+numdofp+NDOF3*j+0]+=Ke1->values[i*numdofp+NDOF2*j+0];
+		Ke->values[i*numdoftotal+numdofp+NDOF3*j+1]+=Ke1->values[i*numdofp+NDOF2*j+1];
+	}
+
+	/*Transform Coordinate System*/ //Do not transform, already sone in the matrices
+	//TransformStiffnessMatrixCoord(Ke,node_list,numnodes,cs_list);
 
 	/*clean-up and return*/
@@ -7119,5 +7127,5 @@
 	return Ke;
 }
-/*}}}*/
+//*}}}*/
 /*FUNCTION Penta::CreateKMatrixDiagnosticHoriz {{{*/
 ElementMatrix* Penta::CreateKMatrixDiagnosticHoriz(void){
@@ -7657,9 +7665,9 @@
 
 	/*compute all stiffness matrices for this element*/
+	ElementMatrix* Ke1=CreateKMatrixDiagnosticFS();
 	int init = this->element_type;
-	this->element_type=P1Enum;
-	ElementMatrix* Ke1=CreateKMatrixDiagnosticHO();
+	this->element_type=P1Enum; //P1 needed for HO 
+	ElementMatrix* Ke2=CreateKMatrixDiagnosticHO();
 	this->element_type=init;
-	ElementMatrix* Ke2=CreateKMatrixDiagnosticFS();
 	ElementMatrix* Ke3=CreateKMatrixCouplingHOFS();
 	ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2,Ke3);
@@ -8230,12 +8238,21 @@
 	inputs->GetInputValue(&approximation,ApproximationEnum);
 	if(approximation!=HOFSApproximationEnum) return NULL;
-	ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters,FSvelocityEnum);
+	int vnumnodes = this->NumberofNodesVelocity();
+	int pnumnodes = this->NumberofNodesPressure();
+
+	/*Prepare coordinate system list*/
+	int* cs_list = xNew<int>(vnumnodes+pnumnodes);
+	for(i=0;i<vnumnodes;i++) cs_list[i]           = XYZEnum;
+	for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
+
+	/*Initialize Element matrix and vectors*/
+	ElementVector* pe=new ElementVector(nodes,vnumnodes+pnumnodes,this->parameters,FSvelocityEnum);
 
 	/*Retrieve all inputs and parameters*/
 	GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
 	this->parameters->FindParam(&FSreconditioning,DiagnosticFSreconditioningEnum);
-	Input* vx_input=inputs->GetInput(VxEnum);               _assert_(vx_input);
-	Input* vy_input=inputs->GetInput(VyEnum);               _assert_(vy_input);
-	Input* vz_input=inputs->GetInput(VzEnum);               _assert_(vz_input);
+	Input* vx_input=inputs->GetInput(VxEnum);       _assert_(vx_input);
+	Input* vy_input=inputs->GetInput(VyEnum);       _assert_(vy_input);
+	Input* vz_input=inputs->GetInput(VzEnum);       _assert_(vz_input);
 	Input* vzHO_input=inputs->GetInput(VzHOEnum);   _assert_(vzHO_input);
 
@@ -8256,15 +8273,16 @@
 
 		for(i=0;i<NUMVERTICES;i++){
-			pe->values[i*NDOF4+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dbasis[2][i];
-			pe->values[i*NDOF4+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dbasis[2][i];
-			pe->values[i*NDOF4+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dbasis[0][i]+dw[1]*dbasis[1][i]+2*dw[2]*dbasis[2][i]);
-			pe->values[i*NDOF4+3]+=Jdet*gauss->weight*FSreconditioning*dw[2]*basis[i];
+			pe->values[i*NDOF3+0]+=-Jdet*gauss->weight*viscosity*dw[0]*dbasis[2][i];
+			pe->values[i*NDOF3+1]+=-Jdet*gauss->weight*viscosity*dw[1]*dbasis[2][i];
+			pe->values[i*NDOF3+2]+=-Jdet*gauss->weight*viscosity*(dw[0]*dbasis[0][i]+dw[1]*dbasis[1][i]+2*dw[2]*dbasis[2][i]);
+			pe->values[NDOF3*vnumnodes+i]+=Jdet*gauss->weight*FSreconditioning*dw[2]*basis[i];
 		}
 	}
 
 	/*Transform coordinate system*/
-	TransformLoadVectorCoord(pe,nodes,NUMVERTICES,XYZEnum);
+	TransformLoadVectorCoord(pe,nodes,vnumnodes+pnumnodes,cs_list);
 
 	/*Clean up and return*/
+	xDelete<int>(cs_list);
 	delete gauss;
 	return pe;
@@ -8407,8 +8425,17 @@
 
 	/*compute all load vectors for this element*/
+	int init = this->element_type;
+	this->element_type=P1Enum;
 	ElementVector* pe1=CreatePVectorDiagnosticHO();
+	this->element_type=init;
 	ElementVector* pe2=CreatePVectorDiagnosticFS();
+	int indices[3]={18,19,20};
+	this->element_type=MINIcondensedEnum;
+	ElementMatrix* Ke = CreateKMatrixDiagnosticFS();
+	this->element_type=init;
+	pe2->StaticCondensation(Ke,3,&indices[0]);
+	delete Ke;
 	ElementVector* pe3=CreatePVectorCouplingHOFS();
-	ElementVector* pe =new ElementVector(pe1,pe2,pe3);
+	ElementVector* pe =new ElementVector(pe1,pe2);
 
 	/*clean-up and return*/
Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp	(revision 15743)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateNodesDiagnosticHoriz.cpp	(revision 15744)
@@ -100,5 +100,7 @@
 			for(int i=0;i<iomodel->numberofelements;i++){
 				if(iomodel->my_elements[i]){
-					nodes->AddObject(new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,0,iomodel,DiagnosticHorizAnalysisEnum,FSvelocityEnum));
+					Node* node = new Node(iomodel->nodecounter+iomodel->numberofvertices+i+1,iomodel->numberofvertices+i,0,iomodel,DiagnosticHorizAnalysisEnum,FSvelocityEnum);
+					node->Deactivate();
+					nodes->AddObject(node);
 				}
 			}