Index: issm/trunk/src/m/model/setflowequation.m
===================================================================
--- issm/trunk/src/m/model/setflowequation.m	(revision 9733)
+++ issm/trunk/src/m/model/setflowequation.m	(revision 11027)
@@ -67,37 +67,21 @@
 end
 
-%add in model who is who
-md.flowequation.element_equation=zeros(md.mesh.numberofelements,1);
-
-%1: Hutter elements
+%Initialize node fields
 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
+nodeonstokes=zeros(md.mesh.numberofvertices,1);
+noneflag=zeros(md.mesh.numberofelements,1);
+
 %First modify stokesflag to get rid of elements contrained everywhere (spc + border with pattyn or macayeal)
 if any(stokesflag),
-	nodeonstokes=zeros(md.mesh.numberofvertices,1);
-	nodeonstokes(md.mesh.elements(find(stokesflag),:))=1;
 	fullspcnodes=double((~isnan(md.diagnostic.spcvx)+~isnan(md.diagnostic.spcvy)+~isnan(md.diagnostic.spcvz))==3 | (nodeonpattyn & nodeonstokes));         %find all the nodes on the boundary of the domain without icefront
 	fullspcelems=double(sum(fullspcnodes(md.mesh.elements),2)==6);         %find all the nodes on the boundary of the domain without icefront
 	stokesflag(find(fullspcelems))=0;
-end
-nodeonstokes=zeros(md.mesh.numberofvertices,1);
-nodeonstokes(md.mesh.elements(find(stokesflag),:))=1;
-md.flowequation.borderstokes=nodeonstokes;
-md.flowequation.element_equation(find(stokesflag))=4;
+	nodeonstokes(md.mesh.elements(find(stokesflag),:))=1;
+end
 
 %Then complete with NoneApproximation or the other model used if there is no stokes
@@ -106,14 +90,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 +103,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 +129,28 @@
 		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
+		nodeonmacayeal(:)=0;
+		nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1;
+
+		%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);
+		macayealflag(pos(pos1))=1;
+		macayealpattynflag(pos(pos1))=0;
+		pos2=find(list==-1);
+		pattynflag(pos(pos2))=1;
+		macayealpattynflag(pos(pos2))=0;
+
+		%Recompute nodes associated to these elements
+		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,16 +161,28 @@
 		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
-		nodeonpattynstokes=zeros(md.mesh.numberofvertices,1);
+
+		%rule out elements that don't touch the 2 boundaries
+		pos=find(pattynstokesflag);
+		list=zeros(length(pos),1);
+		list = list + any(sum(nodeonstokes(md.mesh.elements(pos,:)),2),2);
+		list = list - any(sum(nodeonpattyn(md.mesh.elements(pos,:)),2),2);
+		pos1=find(list==1);
+		stokesflag(pos(pos1))=1;
+		pattynstokesflag(pos(pos1))=0;
+		pos2=find(list==-1);
+		pattynflag(pos(pos2))=1;
+		pattynstokesflag(pos(pos2))=0;
+
+		%Recompute nodes associated to these elements
+		nodeonstokes(:)=0;
+		nodeonstokes(md.mesh.elements(find(stokesflag),:))=1;
+		nodeonpattyn(:)=0;
+		nodeonpattyn(md.mesh.elements(find(pattynflag),:))=1;
+		nodeonpattynstokes(:)=0;
 		nodeonpattynstokes(md.mesh.elements(find(pattynstokesflag),:))=1;
+
 	elseif any(stokesflag) & any(macayealflag),
 		%Find node at the border
@@ -186,20 +193,48 @@
 		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
-		nodeonmacayealstokes=zeros(md.mesh.numberofvertices,1);
+
+		%rule out elements that don't touch the 2 boundaries
+		pos=find(macayealstokesflag);
+		list=zeros(length(pos),1);
+		list = list + any(sum(nodeonmacayeal(md.mesh.elements(pos,:)),2),2);
+		list = list - any(sum(nodeonstokes(md.mesh.elements(pos,:))  ,2),2);
+		pos1=find(list==1);
+		macayealflag(pos(pos1))=1;
+		macayealstokesflag(pos(pos1))=0;
+		pos2=find(list==-1);
+		stokesflag(pos(pos2))=1;
+		macayealstokesflag(pos(pos2))=0;
+
+		%Recompute nodes associated to these elements
+		nodeonmacayeal(:)=0;
+		nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1;
+		nodeonstokes(:)=0;
+		nodeonstokes(md.mesh.elements(find(stokesflag),:))=1;
+		nodeonmacayealstokes(:)=0;
 		nodeonmacayealstokes(md.mesh.elements(find(macayealstokesflag),:))=1;
+
 	elseif any(stokesflag) & any(hutterflag),
 		error('type of coupling not supported yet');
 	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 +248,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 +268,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
