Index: /issm/trunk/src/m/classes/public/ismodelselfconsistent.m
===================================================================
--- /issm/trunk/src/m/classes/public/ismodelselfconsistent.m	(revision 5612)
+++ /issm/trunk/src/m/classes/public/ismodelselfconsistent.m	(revision 5613)
@@ -58,5 +58,5 @@
 checksize(md,fields,[md.numberofelements 1]);
 %Check the values of elements_type
-checkvalues(md,{'elements_type'},[MacAyealApproximationEnum() HutterApproximationEnum() PattynApproximationEnum() MacAyealPattynApproximationEnum() StokesApproximationEnum() NoneApproximationEnum()]);
+checkvalues(md,{'elements_type'},[MacAyealApproximationEnum() HutterApproximationEnum() PattynApproximationEnum() MacAyealPattynApproximationEnum() PattynStokesApproximationEnum() StokesApproximationEnum() NoneApproximationEnum()]);
 if (md.dim==2),
 	checkvalues(md,{'elements_type'},[MacAyealApproximationEnum() HutterApproximationEnum()]);
Index: /issm/trunk/src/m/classes/public/plot/plot_elementstype.m
===================================================================
--- /issm/trunk/src/m/classes/public/plot/plot_elementstype.m	(revision 5612)
+++ /issm/trunk/src/m/classes/public/plot/plot_elementstype.m	(revision 5613)
@@ -69,4 +69,12 @@
 	patch( 'Faces', [B E F C ], 'Vertices', [x y z],'CData', PattynApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
 	patch( 'Faces', [C A D F ], 'Vertices', [x y z],'CData', PattynApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
+	%Stokes elements
+	posS=find(data==StokesApproximationEnum);
+	A=elements(posS,1); B=elements(posS,2); C=elements(posS,3); D=elements(posS,4); E=elements(posS,5); F=elements(posS,6);
+	p4=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData', StokesApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
+	patch( 'Faces', [D E F], 'Vertices', [x y z],'CData', StokesApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
+	patch( 'Faces', [A B E D], 'Vertices', [x y z],'CData', StokesApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
+	patch( 'Faces', [B E F C ], 'Vertices', [x y z],'CData', StokesApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
+	patch( 'Faces', [C A D F ], 'Vertices', [x y z],'CData', StokesApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
 	%MacAyealPattyn elements
 	posP=find(data==MacAyealPattynApproximationEnum);
@@ -77,18 +85,13 @@
 	patch( 'Faces', [B E F C ], 'Vertices', [x y z],'CData', MacAyealPattynApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
 	patch( 'Faces', [C A D F ], 'Vertices', [x y z],'CData', MacAyealPattynApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
-	%Stokes elements
-	alpha=0.35;
-	posS=find(data==StokesApproximationEnum);
-	if ~isempty(posS)
-		A=elements(posS,1); B=elements(posS,2); C=elements(posS,3); D=elements(posS,4); E=elements(posS,5); F=elements(posS,6);
-		p4=patch( 'Faces', [A B C], 'Vertices', [x y z],'FaceColor','none','EdgeColor','black');
-		patch( 'Faces', [D E F], 'Vertices', [x y z],'FaceColor','none','EdgeColor','black');
-		patch( 'Faces', [A B E D], 'Vertices', [x y z],'FaceColor','none','EdgeColor','black');
-		patch( 'Faces', [B E F C ], 'Vertices', [x y z],'FaceColor','none','EdgeColor','black');
-		patch( 'Faces', [C A D F ], 'Vertices', [x y z],'FaceColor','none','EdgeColor','black');
-		legend([p1 p2 p3 p5 p4],'Hutter''s elements','MacAyeal''s elements','Pattyn''s elements','MacAyealPattyn''s elements','Stokes''s elements');
-	else
-		legend([p1 p2 p3 p5],'Hutter''s elements','MacAyeal''s elements','Pattyn''s elements','MacAyealPattyn''s elements');
-	end
+	%PattynStokes elements
+	PosPS=find(data==PattynStokesApproximationEnum);
+	A=elements(PosPS,1); B=elements(PosPS,2); C=elements(PosPS,3); D=elements(PosPS,4); E=elements(PosPS,5); F=elements(PosPS,6);
+	p6=patch( 'Faces', [A B C], 'Vertices', [x y z],'CData', PattynStokesApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
+	patch( 'Faces', [D E F], 'Vertices', [x y z],'CData', PattynStokesApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
+	patch( 'Faces', [A B E D], 'Vertices', [x y z],'CData', PattynStokesApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
+	patch( 'Faces', [B E F C ], 'Vertices', [x y z],'CData', PattynStokesApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
+	patch( 'Faces', [C A D F ], 'Vertices', [x y z],'CData', PattynStokesApproximationEnum,'FaceColor','flat','EdgeColor',edgecolor);
+	legend([p1 p2 p3 p4 p5 p6],'Hutter''s elements','MacAyeal''s elements','Pattyn''s elements','Stokes''s elements','MacAyealPattyn''s elements','PattynStokes''s elements');
 end
 
Index: /issm/trunk/src/m/classes/public/setelementstype.m
===================================================================
--- /issm/trunk/src/m/classes/public/setelementstype.m	(revision 5612)
+++ /issm/trunk/src/m/classes/public/setelementstype.m	(revision 5613)
@@ -43,5 +43,5 @@
 	macayealflag(find(~hutterflag & ~pattynflag))=1;
 elseif strcmpi(filltype,'pattyn'),
-	pattynflag(find(~hutterflag & ~macayealflag))=1;
+	pattynflag(find(~hutterflag & ~macayealflag & ~stokesflag))=1;
 end
 
@@ -52,9 +52,9 @@
 
 %check that each element has only one flag
-if any(hutterflag+ macayealflag+pattynflag>1),
+if any(hutterflag+ macayealflag+pattynflag+stokesflag>1),
 	disp('setelementstype warning message: some elements have several types, higher order type is used for them')
 	hutterflag(find(hutterflag & macayealflag))=0;
+	hutterflag(find(hutterflag & pattynflag))=0;
 	macayealflag(find(macayealflag & pattynflag))=0;
-	pattynflag(find(pattynflag & macayealflag))=0;
 end
 
@@ -67,15 +67,8 @@
 
 %Stokes can only be used alone for now:
-if any(stokesflag) &any(hutterflag+pattynflag+macayealflag),
+if any(stokesflag) &any(hutterflag+macayealflag),
 	error('setelementstype error message: stokes cannot be used with any other model for now, put stokes everywhere')
 end
 
-if any(stokesflag),
-	%Modify stokesflag to get rid of elements contrained everywhere
-	fullspcnodes=double(sum(md.spcvelocity(:,1:3),2)==3);         %find all the grids on the boundary of the domain without icefront
-	fullspcelems=double(sum(fullspcnodes(md.elements),2)==6);         %find all the grids on the boundary of the domain without icefront
-	stokesflag(find(fullspcelems))=0;
-end
-	
 %add in model who is who
 md.elements_type=zeros(md.numberofelements,1);
@@ -100,4 +93,12 @@
 
 %4: Stokes elements
+%First modify stokesflag to get rid of elements contrained everywhere (spc + border with pattyn or macayeal)
+if any(stokesflag),
+	gridonstokes=zeros(md.numberofgrids,1);
+	gridonstokes(md.elements(find(stokesflag),:))=1;
+	fullspcnodes=double(sum(md.spcvelocity(:,1:3),2)==3 | (gridonpattyn & gridonstokes));         %find all the grids on the boundary of the domain without icefront
+	fullspcelems=double(sum(fullspcnodes(md.elements),2)==6);         %find all the grids on the boundary of the domain without icefront
+	stokesflag(find(fullspcelems))=0;
+end
 gridonstokes=zeros(md.numberofgrids,1);
 gridonstokes(md.elements(find(stokesflag),:))=1;
@@ -105,14 +106,22 @@
 md.elements_type(find(stokesflag))=StokesApproximationEnum();
 
-%Complete with NoneApproximation if there is no stokes
+%Then complete with NoneApproximation or the other model used if there is no stokes
 if any(stokesflag), 
+	if any(pattynflag), %fill with pattyn
+		pattynflag(~stokesflag)=1;
+		gridonpattyn(md.elements(find(pattynflag),:))=1;
+		md.gridonpattyn=gridonpattyn;
+		md.elements_type(find(~stokesflag))=PattynApproximationEnum();
+	else %fill with none 
 	%5: None elements (non Stokes)
-	md.elements_type(find(~stokesflag))=NoneApproximationEnum();
+		md.elements_type(find(~stokesflag))=NoneApproximationEnum();
+	end
 end
 
 %Now take care of the coupling between MacAyeal and Pattyn
 md.penalties=[];
+gridonmacayealpattyn=zeros(md.numberofgrids,1);
+gridonpattynstokes=zeros(md.numberofgrids,1);
 if strcmpi(coupling_method,'penalties'),
-	gridonmacayealpattyn=zeros(md.numberofgrids,1);
 	%Create the border grids between Pattyn and MacAyeal and extrude them
 	numgrids2d=md.numberofgrids2d;
@@ -129,23 +138,45 @@
 	end
 elseif strcmpi(coupling_method,'tiling'),
-	%Find grid at the border
-	gridonmacayealpattyn=zeros(md.numberofgrids,1);
-	gridonmacayealpattyn(find(gridonmacayeal & gridonpattyn))=1;
-	%Macayeal elements in contact with this layer become MacAyealPattyn elements
-	matrixelements=ismember(md.elements,find(gridonmacayealpattyn));
-	commonelements=sum(matrixelements,2)~=0;
-	commonelements(find(pattynflag))=0; %only one layer: the elements previously in macayeal
-	macayealflag(find(commonelements))=0; %these elements are now macayealpattynelements
-	macayealpattynflag=zeros(md.numberofelements,1);
-	macayealpattynflag(find(commonelements))=1;
-	gridonmacayeal=zeros(md.numberofgrids,1);
-	gridonmacayeal(md.elements(find(macayealflag),:))=1;
-	md.gridonmacayeal=gridonmacayeal;
-
-	%Create MacaAyealPattynApproximation where needed
-	md.elements_type(find(macayealpattynflag))=MacAyealPattynApproximationEnum();
-
-	%Now recreate gridonmacayeal and gridonmacayealpattyn
-	gridonmacayealpattyn(md.elements(find(macayealpattynflag),:))=1;
+	if any(macayealflag) & any(pattynflag), %coupling macayeal pattyn
+		%Find grid at the border
+		gridonmacayealpattyn(find(gridonmacayeal & gridonpattyn))=1;
+		%Macayeal elements in contact with this layer become MacAyealPattyn elements
+		matrixelements=ismember(md.elements,find(gridonmacayealpattyn));
+		commonelements=sum(matrixelements,2)~=0;
+		commonelements(find(pattynflag))=0; %only one layer: the elements previously in macayeal
+		macayealflag(find(commonelements))=0; %these elements are now macayealpattynelements
+		macayealpattynflag=zeros(md.numberofelements,1);
+		macayealpattynflag(find(commonelements))=1;
+		gridonmacayeal=zeros(md.numberofgrids,1);
+		gridonmacayeal(md.elements(find(macayealflag),:))=1;
+		md.gridonmacayeal=gridonmacayeal;
+
+		%Create MacaAyealPattynApproximation where needed
+		md.elements_type(find(macayealpattynflag))=MacAyealPattynApproximationEnum();
+
+		%Now recreate gridonmacayeal and gridonmacayealpattyn
+		gridonmacayealpattyn(md.elements(find(macayealpattynflag),:))=1;
+	elseif any(pattynflag) & any(stokesflag), %coupling pattyn stokes
+		%Find grid at the border
+		gridonpattynstokes(find(gridonpattyn & gridonstokes))=1;
+		%Stokes elements in contact with this layer become PattynStokes elements
+		matrixelements=ismember(md.elements,find(gridonpattynstokes));
+		commonelements=sum(matrixelements,2)~=0;
+		commonelements(find(pattynflag))=0; %only one layer: the elements previously in macayeal
+		stokesflag(find(commonelements))=0; %these elements are now macayealpattynelements
+		pattynstokesflag=zeros(md.numberofelements,1);
+		pattynstokesflag(find(commonelements))=1;
+		gridonstokes=zeros(md.numberofgrids,1);
+		gridonstokes(md.elements(find(stokesflag),:))=1;
+		md.gridonstokes=gridonstokes;
+
+		%Create MacaAyealPattynApproximation where needed
+		md.elements_type(find(pattynstokesflag))=PattynStokesApproximationEnum();
+
+		%Now recreate gridonmacayeal and gridonmacayealpattyn
+		gridonpattynstokes(md.elements(find(pattynstokesflag),:))=1;
+	elseif any(stokesflag) & any(macayealflag+hutterflag),
+		error('type of coupling not supported yet');
+	end
 end
 
@@ -164,4 +195,6 @@
 pos=find(gridonmacayealpattyn);
 md.vertices_type(pos)=MacAyealPattynApproximationEnum();
+pos=find(gridonpattynstokes);
+md.vertices_type(pos)=PattynStokesApproximationEnum();
 pos=find(gridonstokes);
 md.vertices_type(pos)=StokesApproximationEnum();