Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp
===================================================================
--- /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp	(revision 15623)
+++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/DiagnosticHoriz/CreateConstraintsDiagnosticHoriz.cpp	(revision 15624)
@@ -6,4 +6,5 @@
 #include "../../../classes/classes.h"
 #include "../../../shared/shared.h"
+#include "../../IoModelToConstraintsx/IoModelToConstraintsx.h"
 #include "../ModelProcessorx.h"
 
@@ -12,10 +13,9 @@
 	/*Intermediary*/
 	int        i,j;
-	int        count;
+	int        count,temp,finiteelement;
 	IssmDouble g;
 	IssmDouble rho_ice;
 	IssmDouble FSreconditioning;
-	bool       isSSA,isL1L2,isHO,isFS;
-	int        fe_ssa,fe_ho;
+	bool       isSSA,isL1L2,isHO,isFS,iscoupling;
 	bool       spcpresent = false;
 	int        Mx,Nx;
@@ -50,6 +50,4 @@
 	iomodel->Constant(&isHO,FlowequationIsHOEnum);
 	iomodel->Constant(&isFS,FlowequationIsFSEnum);
-	iomodel->Constant(&fe_ssa,FlowequationFeSSAEnum);
-	iomodel->Constant(&fe_ho,FlowequationFeHOEnum);
 
 	/*Recover pointer: */
@@ -58,4 +56,43 @@
 	/*Now, is the flag macayaealHO on? otherwise, do nothing: */
 	if(!isSSA & !isHO & !isFS & !isL1L2){
+		*pconstraints=constraints;
+		return;
+	}
+
+	/*Do we have coupling*/
+	if( (isSSA?1.:0.) + (isL1L2?1.:0.) + (isHO?1.:0.) + (isFS?1.:0.) >1.)
+	 iscoupling = true;
+	else
+	 iscoupling = false;
+
+	/*If no coupling, call Regular IoModelToConstraintsx, else, use P1 elements only*/
+	if(!iscoupling && !isFS){
+
+		/*Get finite element type*/
+		if(isSSA){
+			iomodel->Constant(&temp,FlowequationFeSSAEnum);
+			switch(temp){
+				case 0 : finiteelement = P1Enum; break;
+				case 1 : finiteelement = P2Enum; break;
+				default: _error_("finite element "<<temp<<" not supported");
+			}
+		}
+		else if(isL1L2){
+			finiteelement = P1Enum;
+		}
+		else if(isHO){
+			iomodel->Constant(&temp,FlowequationFeHOEnum);
+			switch(temp){
+				case 0 : finiteelement = P1Enum; break;
+				case 1 : finiteelement = P2Enum; break;
+				default: _error_("finite element "<<temp<<" not supported");
+			}
+		}
+		else if(isFS){
+			finiteelement = P1Enum;
+		}
+		IoModelToConstraintsx(constraints,iomodel,DiagnosticSpcvxEnum,DiagnosticHorizAnalysisEnum,finiteelement,1);
+		IoModelToConstraintsx(constraints,iomodel,DiagnosticSpcvyEnum,DiagnosticHorizAnalysisEnum,finiteelement,2);
+
 		*pconstraints=constraints;
 		return;
@@ -99,5 +136,5 @@
 
 			/*Start with adding spcs of coupling: zero at the border SSA/HO for the appropriate dofs*/
-			if (reCast<int,IssmDouble>(vertices_type[i]==SSAHOApproximationEnum)){
+			if(reCast<int,IssmDouble>(vertices_type[i]==SSAHOApproximationEnum)){
 				/*If grionSSA, spc HO dofs: 3 & 4*/
 					if (reCast<int,IssmDouble>(nodeonHO[i])){
@@ -309,44 +346,4 @@
 			}
 		}
-	}
-
-	/*SPC Quadratic elements*/
-	if((isSSA && fe_ssa==1) || (isHO && fe_ho==1)){
-
-		int   v1,v2;
-		bool *my_edges = NULL;
-
-		if(Mx!=iomodel->numberofvertices) _error_("transient spc not supported yet");
-		EdgesPartitioning(&my_edges,iomodel);
-
-		for(i=0;i<iomodel->numberofedges;i++){
-
-			if(my_edges[i]){
-
-				v1 = iomodel->edges[2*i+0]-1;
-				v2 = iomodel->edges[2*i+1]-1;
-
-				if(!xIsNan<IssmDouble>(spcvx[v1]) && !xIsNan<IssmDouble>(spcvx[v2])){
-					constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1,
-									1,(spcvx[v1]+spcvx[v2])/(2.),DiagnosticHorizAnalysisEnum));
-					count++;
-				}
-				if(!xIsNan<IssmDouble>(spcvy[v1]) && !xIsNan<IssmDouble>(spcvy[v2])){
-					constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1,
-									2,(spcvy[v1]+spcvy[v2])/(2.),DiagnosticHorizAnalysisEnum));
-					count++;
-				}
-				if (reCast<int,IssmDouble>(vertices_type[v1])==FSApproximationEnum ||  (reCast<int,IssmDouble>(vertices_type[v1])==NoneApproximationEnum)){
-					if(!xIsNan<IssmDouble>(spcvz[v1]) && !xIsNan<IssmDouble>(spcvz[v2])){
-						constraints->AddObject(new SpcStatic(iomodel->constraintcounter+count+1,iomodel->nodecounter+iomodel->numberofvertices+i+1,
-										3,(spcvz[v1]+spcvz[v2])/(2.),DiagnosticHorizAnalysisEnum));
-						count++;
-					}
-				}
-			}
-		}
-
-		/*Clean up*/
-		xDelete<bool>(my_edges);
 	}