Changeset 5379

Timestamp:

08/18/10 15:19:12 (15 years ago)

Author:

seroussi

Message:

added tiling capability in setelements_type

File:

• issm/trunk/src/m/classes/public/setelementstype.m (modified) (5 diffs)

issm/trunk/src/m/classes/public/setelementstype.m

@@ -10 +10 @@
 %   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 +17 @@
 %   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 +26 @@
 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 +107 @@
 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 +159 @@
 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)