Index: /issm/trunk-jpl/src/m/model/setflowequation.m
===================================================================
--- /issm/trunk-jpl/src/m/model/setflowequation.m	(revision 11002)
+++ /issm/trunk-jpl/src/m/model/setflowequation.m	(revision 11003)
@@ -67,23 +67,15 @@
 end
 
-%add in model who is who
-md.flowequation.element_equation=zeros(md.mesh.numberofelements,1);
-
 %1: Hutter elements
 nodeonhutter=zeros(md.mesh.numberofvertices,1);
 nodeonhutter(md.mesh.elements(find(hutterflag),:))=1;
-md.flowequation.element_equation(find(hutterflag))=1;
 
 %2: MacAyeal elements
 nodeonmacayeal=zeros(md.mesh.numberofvertices,1);
 nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1;
-md.flowequation.bordermacayeal=nodeonmacayeal;
-md.flowequation.element_equation(find(macayealflag))=2;
 
 %3: Pattyn elements
 nodeonpattyn=zeros(md.mesh.numberofvertices,1);
 nodeonpattyn(md.mesh.elements(find(pattynflag),:))=1;
-md.flowequation.borderpattyn=nodeonpattyn;
-md.flowequation.element_equation(find(pattynflag))=3;
 
 %4: Stokes elements
@@ -98,6 +90,7 @@
 nodeonstokes=zeros(md.mesh.numberofvertices,1);
 nodeonstokes(md.mesh.elements(find(stokesflag),:))=1;
-md.flowequation.borderstokes=nodeonstokes;
-md.flowequation.element_equation(find(stokesflag))=4;
+
+%5: None elements
+noneflag=zeros(md.mesh.numberofelements,1);
 
 %Then complete with NoneApproximation or the other model used if there is no stokes
@@ -106,14 +99,9 @@
 		pattynflag(~stokesflag)=1;
 		nodeonpattyn(md.mesh.elements(find(pattynflag),:))=1;
-		md.flowequation.borderpattyn=nodeonpattyn;
-		md.flowequation.element_equation(find(~stokesflag))=3;
 	elseif any(macayealflag), %fill with macayeal
 		macayealflag(~stokesflag)=1;
 		nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1;
-		md.flowequation.bordermacayeal=nodeonmacayeal;
-		md.flowequation.element_equation(find(~stokesflag))=2;
 	else %fill with none 
-	%5: None elements (non Stokes)
-		md.flowequation.element_equation(find(~stokesflag))=0;
+		noneflag(find(~stokesflag))=1;
 	end
 end
@@ -124,4 +112,7 @@
 nodeonpattynstokes=zeros(md.mesh.numberofvertices,1);
 nodeonmacayealstokes=zeros(md.mesh.numberofvertices,1);
+macayealpattynflag=zeros(md.mesh.numberofelements,1);
+macayealstokesflag=zeros(md.mesh.numberofelements,1);
+pattynstokesflag=zeros(md.mesh.numberofelements,1);
 if strcmpi(coupling_method,'penalties'),
 	%Create the border nodes between Pattyn and MacAyeal and extrude them
@@ -147,15 +138,29 @@
 		commonelements(find(pattynflag))=0; %only one layer: the elements previously in macayeal
 		macayealflag(find(commonelements))=0; %these elements are now macayealpattynelements
-		macayealpattynflag=zeros(md.mesh.numberofelements,1);
 		macayealpattynflag(find(commonelements))=1;
 		nodeonmacayeal=zeros(md.mesh.numberofvertices,1);
 		nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1;
-		md.flowequation.bordermacayeal=nodeonmacayeal;
-
-		%Create MacaAyealPattynApproximation where needed
-		md.flowequation.element_equation(find(macayealpattynflag))=5;
-
-		%Now recreate nodeonmacayeal and nodeonmacayealpattyn
+
+		%rule out elements that don't touch the 2 boundaries
+		pos=find(macayealpattynflag);
+		list=zeros(length(pos),1);
+		list = list + any(sum(nodeonmacayeal(md.mesh.elements(pos,:)),2),2);
+		list = list - any(sum(nodeonpattyn(md.mesh.elements(pos,:))  ,2),2);
+		pos1=find(list==1);
+		pattynflag(pos(pos1))=0;
+		macayealflag(pos(pos1))=1;
+		macayealpattynflag(pos(pos1))=0;
+		pos2=find(list==-1);
+		pattynflag(pos(pos2))=1;
+		macayealflag(pos(pos2))=0;
+		macayealpattynflag(pos(pos2))=0;
+
+		nodeonmacayeal(:)=0;
+		nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1;
+		nodeonpattyn(:)=0;
+		nodeonpattyn(md.mesh.elements(find(pattynflag),:))=1;
+		nodeonmacayealpattyn(:)=0;
 		nodeonmacayealpattyn(md.mesh.elements(find(macayealpattynflag),:))=1;
+
 	elseif any(pattynflag) & any(stokesflag), %coupling pattyn stokes
 		%Find node at the border
@@ -166,12 +171,7 @@
 		commonelements(find(pattynflag))=0; %only one layer: the elements previously in macayeal
 		stokesflag(find(commonelements))=0; %these elements are now macayealpattynelements
-		pattynstokesflag=zeros(md.mesh.numberofelements,1);
 		pattynstokesflag(find(commonelements))=1;
 		nodeonstokes=zeros(md.mesh.numberofvertices,1);
 		nodeonstokes(md.mesh.elements(find(stokesflag),:))=1;
-		md.flowequation.borderstokes=nodeonstokes;
-
-		%Create MacaAyealPattynApproximation where needed
-		md.flowequation.element_equation(find(pattynstokesflag))=7;
 
 		%Now recreate nodeonpattynstokes
@@ -186,12 +186,7 @@
 		commonelements(find(macayealflag))=0; %only one layer: the elements previously in macayeal
 		stokesflag(find(commonelements))=0; %these elements are now macayealmacayealelements
-		macayealstokesflag=zeros(md.mesh.numberofelements,1);
 		macayealstokesflag(find(commonelements))=1;
 		nodeonstokes=zeros(md.mesh.numberofvertices,1);
 		nodeonstokes(md.mesh.elements(find(stokesflag),:))=1;
-		md.flowequation.borderstokes=nodeonstokes;
-
-		%Create MacaAyeal Approximation where needed
-		md.flowequation.element_equation(find(macayealstokesflag))=6;
 
 		%Now recreate nodeonmacayealstokes
@@ -202,4 +197,20 @@
 	end
 end
+
+%Create MacaAyealPattynApproximation where needed
+md.flowequation.element_equation=zeros(md.mesh.numberofelements,1);
+md.flowequation.element_equation(find(noneflag))=0;
+md.flowequation.element_equation(find(hutterflag))=1;
+md.flowequation.element_equation(find(macayealflag))=2;
+md.flowequation.element_equation(find(pattynflag))=3;
+md.flowequation.element_equation(find(stokesflag))=4;
+md.flowequation.element_equation(find(macayealpattynflag))=5;
+md.flowequation.element_equation(find(macayealstokesflag))=6;
+md.flowequation.element_equation(find(pattynstokesflag))=7;
+
+%border
+md.flowequation.borderpattyn=nodeonpattyn;
+md.flowequation.bordermacayeal=nodeonmacayeal;
+md.flowequation.borderstokes=nodeonstokes;
 
 %Create vertices_type
@@ -213,6 +224,4 @@
 pos=find(nodeonhutter);
 md.flowequation.vertex_equation(pos)=1;
-pos=find(nodeonpattyn & nodeonmacayeal);
-md.flowequation.vertex_equation(pos)=3;
 pos=find(nodeonmacayealpattyn);
 md.flowequation.vertex_equation(pos)=5;
@@ -235,3 +244,14 @@
 md.flowequation.isstokes=double(any(md.flowequation.element_equation==4));
 
-end
+return
+
+%Check that tiling can work:
+if any(md.flowequation.bordermacayeal) & any(md.flowequation.borderpattyn) & any(md.flowequation.borderpattyn + md.flowequation.bordermacayeal ~=1),
+	error('error coupling domain too irregular');
+end
+if any(md.flowequation.bordermacayeal) & any(md.flowequation.borderstokes) & any(md.flowequation.borderstokes + md.flowequation.bordermacayeal ~=1),
+	error('error coupling domain too irregular');
+end
+if any(md.flowequation.borderstokes) & any(md.flowequation.borderpattyn) & any(md.flowequation.borderpattyn + md.flowequation.borderstokes~=1),
+	error('error coupling domain too irregular');
+end