Index: /issm/trunk/src/m/classes/public/setelementstype.m
===================================================================
--- /issm/trunk/src/m/classes/public/setelementstype.m	(revision 5378)
+++ /issm/trunk/src/m/classes/public/setelementstype.m	(revision 5379)
@@ -10,4 +10,5 @@
 %   an empty string '' will be considered as an empty domain
 %   a string 'all' will be considered as the entire domain
+%   You can specify the type of coupling, 'penalties' or 'tiling', to use with the input 'coupling'
 %
 %   Usage:
@@ -16,4 +17,5 @@
 %   Example:
 %      md=setelementstype(md,'pattyn','Pattyn.exp','macayeal',md.elementoniceshelf,'fill','hutter');
+%      md=setelementstype(md,'pattyn','Pattyn.exp',fill','hutter','coupling','tiling');
 
 %some checks on list of arguments
@@ -24,4 +26,11 @@
 if md.counter<3,
 	error('only fully parameterized 2d models can be setelementstyped');
+end
+
+%Find_out what kind of coupling to use
+options=pairoptions(varargin{:});
+coupling_method=getfieldvalue(options,'coupling','tiling');
+if (~strcmpi(coupling_method,'tiling') & ~strcmpi(coupling_method,'penalties')),
+	error('coupling type can only be: tiling or penalties');
 end
 
@@ -98,4 +107,44 @@
 md.elements_type(find(~stokesflag),2)=NoneApproximationEnum();
 
+%Now take care of the coupling between MacAyeal and Pattyn
+md.penalties=[];
+if strcmpi(coupling_method,'penalties'),
+	gridonmacayealpattyn=zeros(md.numberofgrids,1);
+	%Create the border grids between Pattyn and MacAyeal and extrude them
+	numgrids2d=md.numberofgrids2d;
+	numlayers=md.numlayers;
+	bordergrids2d=find(gridonpattyn(1:numgrids2d) & gridonmacayeal(1:numgrids2d)); %Grids connected to two different types of elements
+
+	%initialize and fill in penalties structure
+	if ~isnan(bordergrids2d),
+		penalties=[];
+		for	i=1:numlayers-1,
+			penalties=[penalties; [bordergrids2d bordergrids2d+md.numberofgrids2d*(i)]];
+		end
+		md.penalties=penalties;
+	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))=1; %these elements are now macayealpattynelements
+	macayealpattynflag=zeros(md.numberofelements);
+	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),1)=MacAyealPattynApproximationEnum();
+
+	%Now recreate gridonmacayeal and gridonmacayealpattyn
+	gridonmacayeal(md.elements(find(macayealflag),:))=1;
+	gridonmacayealpattyn(md.elements(find(macayealpattynflag),:))=1;
+end
+
 %Create vertices_type
 md.vertices_type=zeros(md.numberofgrids,2);
@@ -110,23 +159,10 @@
 pos=find(gridonpattyn & gridonmacayeal);
 md.vertices_type(pos,1)=PattynApproximationEnum();
+pos=find(gridonmacayealpattyn);
+md.vertices_type(pos,1)=MacAyealPattynApproximationEnum();
 pos=find(gridonstokes);
 md.vertices_type(pos,2)=StokesApproximationEnum();
 pos=find(~gridonstokes);
 md.vertices_type(pos,2)=NoneApproximationEnum();
-
-%Create the border grids between Pattyn and MacAyeal and extrude them
-numgrids2d=md.numberofgrids2d;
-numlayers=md.numlayers;
-bordergrids2d=find(gridonpattyn(1:numgrids2d) & gridonmacayeal(1:numgrids2d)); %Grids connected to two different types of elements
-
-%initialize and fill in penalties structure
-md.penalties=[];
-if ~isnan(bordergrids2d),
-	penalties=[];
-	for	i=1:numlayers-1,
-		penalties=[penalties; [bordergrids2d bordergrids2d+md.numberofgrids2d*(i)]];
-	end
-	md.penalties=penalties;
-end
 
 %flag dead grids (strictly in MacAyeal and not on bed -> not used in diagnostic horiz)