Index: /issm/trunk/src/m/model/setelementstype.m
===================================================================
--- /issm/trunk/src/m/model/setelementstype.m	(revision 6704)
+++ /issm/trunk/src/m/model/setelementstype.m	(revision 6705)
@@ -67,5 +67,5 @@
 
 %Stokes can only be used alone for now:
-if any(stokesflag) &any(hutterflag+macayealflag),
+if any(stokesflag) &any(hutterflag),
 	error('setelementstype error message: stokes cannot be used with any other model for now, put stokes everywhere')
 end
@@ -123,4 +123,5 @@
 gridonmacayealpattyn=zeros(md.numberofgrids,1);
 gridonpattynstokes=zeros(md.numberofgrids,1);
+gridonmacayealstokes=zeros(md.numberofgrids,1);
 if strcmpi(coupling_method,'penalties'),
 	%Create the border grids between Pattyn and MacAyeal and extrude them
@@ -177,5 +178,25 @@
 		gridonpattynstokes=zeros(md.numberofgrids,1);
 		gridonpattynstokes(md.elements(find(pattynstokesflag),:))=1;
-	elseif any(stokesflag) & any(macayealflag+hutterflag),
+	elseif any(stokesflag) & any(macayealflag),
+		%Find grid at the border
+		gridonmacayealstokes(find(gridonmacayeal & gridonstokes))=1;
+		%Stokes elements in contact with this layer become MacAyealStokes elements
+		matrixelements=ismember(md.elements,find(gridonmacayealstokes));
+		commonelements=sum(matrixelements,2)~=0;
+		commonelements(find(macayealflag))=0; %only one layer: the elements previously in macayeal
+		stokesflag(find(commonelements))=0; %these elements are now macayealmacayealelements
+		macayealstokesflag=zeros(md.numberofelements,1);
+		macayealstokesflag(find(commonelements))=1;
+		gridonstokes=zeros(md.numberofgrids,1);
+		gridonstokes(md.elements(find(stokesflag),:))=1;
+		md.gridonstokes=gridonstokes;
+
+		%Create MacaAyealMacAyealApproximation where needed
+		md.elements_type(find(macayealstokesflag))=MacAyealStokesApproximationEnum();
+
+		%Now recreate gridonmacayealstokes
+		gridonmacayealstokes=zeros(md.numberofgrids,1);
+		gridonmacayealstokes(md.elements(find(macayealstokesflag),:))=1;
+	elseif any(stokesflag) & any(hutterflag),
 		error('type of coupling not supported yet');
 	end
@@ -206,4 +227,6 @@
 pos=find(gridonpattynstokes);
 md.vertices_type(pos)=PattynStokesApproximationEnum();
+pos=find(gridonmacayealstokes);
+md.vertices_type(pos)=MacAyealStokesApproximationEnum();
 
 %figure out solution types