Index: /issm/trunk-jpl/src/m/model/setflowequation.m
===================================================================
--- /issm/trunk-jpl/src/m/model/setflowequation.m	(revision 11003)
+++ /issm/trunk-jpl/src/m/model/setflowequation.m	(revision 11004)
@@ -67,30 +67,22 @@
 end
 
-%1: Hutter elements
+%Initialize node fields
 nodeonhutter=zeros(md.mesh.numberofvertices,1);
 nodeonhutter(md.mesh.elements(find(hutterflag),:))=1;
-
-%2: MacAyeal elements
 nodeonmacayeal=zeros(md.mesh.numberofvertices,1);
 nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1;
-
-%3: Pattyn elements
 nodeonpattyn=zeros(md.mesh.numberofvertices,1);
 nodeonpattyn(md.mesh.elements(find(pattynflag),:))=1;
-
-%4: Stokes elements
+nodeonstokes=zeros(md.mesh.numberofvertices,1);
+nodeonstokes(md.mesh.elements(find(stokesflag),:))=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;
-
-%5: None elements
-noneflag=zeros(md.mesh.numberofelements,1);
+	nodeonstokes(md.mesh.elements(find(stokesflag),:))=1;
+end
 
 %Then complete with NoneApproximation or the other model used if there is no stokes
@@ -139,5 +131,5 @@
 		macayealflag(find(commonelements))=0; %these elements are now macayealpattynelements
 		macayealpattynflag(find(commonelements))=1;
-		nodeonmacayeal=zeros(md.mesh.numberofvertices,1);
+		nodeonmacayeal(:)=0;
 		nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1;
 
@@ -148,12 +140,11 @@
 		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;
 
+		%Recompute nodes associated to these elements
 		nodeonmacayeal(:)=0;
 		nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1;
@@ -175,7 +166,24 @@
 		nodeonstokes(md.mesh.elements(find(stokesflag),:))=1;
 
-		%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
@@ -190,7 +198,24 @@
 		nodeonstokes(md.mesh.elements(find(stokesflag),:))=1;
 
-		%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');