Changeset 5379


Ignore:
Timestamp:
08/18/10 15:19:12 (15 years ago)
Author:
seroussi
Message:

added tiling capability in setelements_type

File:
1 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk/src/m/classes/public/setelementstype.m

    r5204 r5379  
    1010%   an empty string '' will be considered as an empty domain
    1111%   a string 'all' will be considered as the entire domain
     12%   You can specify the type of coupling, 'penalties' or 'tiling', to use with the input 'coupling'
    1213%
    1314%   Usage:
     
    1617%   Example:
    1718%      md=setelementstype(md,'pattyn','Pattyn.exp','macayeal',md.elementoniceshelf,'fill','hutter');
     19%      md=setelementstype(md,'pattyn','Pattyn.exp',fill','hutter','coupling','tiling');
    1820
    1921%some checks on list of arguments
     
    2426if md.counter<3,
    2527        error('only fully parameterized 2d models can be setelementstyped');
     28end
     29
     30%Find_out what kind of coupling to use
     31options=pairoptions(varargin{:});
     32coupling_method=getfieldvalue(options,'coupling','tiling');
     33if (~strcmpi(coupling_method,'tiling') & ~strcmpi(coupling_method,'penalties')),
     34        error('coupling type can only be: tiling or penalties');
    2635end
    2736
     
    98107md.elements_type(find(~stokesflag),2)=NoneApproximationEnum();
    99108
     109%Now take care of the coupling between MacAyeal and Pattyn
     110md.penalties=[];
     111if strcmpi(coupling_method,'penalties'),
     112        gridonmacayealpattyn=zeros(md.numberofgrids,1);
     113        %Create the border grids between Pattyn and MacAyeal and extrude them
     114        numgrids2d=md.numberofgrids2d;
     115        numlayers=md.numlayers;
     116        bordergrids2d=find(gridonpattyn(1:numgrids2d) & gridonmacayeal(1:numgrids2d)); %Grids connected to two different types of elements
     117
     118        %initialize and fill in penalties structure
     119        if ~isnan(bordergrids2d),
     120                penalties=[];
     121                for     i=1:numlayers-1,
     122                        penalties=[penalties; [bordergrids2d bordergrids2d+md.numberofgrids2d*(i)]];
     123                end
     124                md.penalties=penalties;
     125        end
     126elseif strcmpi(coupling_method,'tiling'),
     127        %Find grid at the border
     128        gridonmacayealpattyn=zeros(md.numberofgrids,1);
     129        gridonmacayealpattyn(find(gridonmacayeal & gridonpattyn))=1;
     130        %Macayeal elements in contact with this layer become MacAyealPattyn elements
     131        matrixelements=ismember(md.elements,find(gridonmacayealpattyn));
     132        commonelements=sum(matrixelements,2)~=0;
     133        commonelements(find(pattynflag))=0; %only one layer: the elements previously in macayeal
     134        macayealflag(find(commonelements))=1; %these elements are now macayealpattynelements
     135        macayealpattynflag=zeros(md.numberofelements);
     136        macayealpattynflag(find(commonelements))=1;
     137        gridonmacayeal=zeros(md.numberofgrids,1);
     138        gridonmacayeal(md.elements(find(macayealflag),:))=1;
     139        md.gridonmacayeal=gridonmacayeal;
     140
     141        %Create MacaAyealPattynApproximation where needed
     142        md.elements_type(find(macayealpattynflag),1)=MacAyealPattynApproximationEnum();
     143
     144        %Now recreate gridonmacayeal and gridonmacayealpattyn
     145        gridonmacayeal(md.elements(find(macayealflag),:))=1;
     146        gridonmacayealpattyn(md.elements(find(macayealpattynflag),:))=1;
     147end
     148
    100149%Create vertices_type
    101150md.vertices_type=zeros(md.numberofgrids,2);
     
    110159pos=find(gridonpattyn & gridonmacayeal);
    111160md.vertices_type(pos,1)=PattynApproximationEnum();
     161pos=find(gridonmacayealpattyn);
     162md.vertices_type(pos,1)=MacAyealPattynApproximationEnum();
    112163pos=find(gridonstokes);
    113164md.vertices_type(pos,2)=StokesApproximationEnum();
    114165pos=find(~gridonstokes);
    115166md.vertices_type(pos,2)=NoneApproximationEnum();
    116 
    117 %Create the border grids between Pattyn and MacAyeal and extrude them
    118 numgrids2d=md.numberofgrids2d;
    119 numlayers=md.numlayers;
    120 bordergrids2d=find(gridonpattyn(1:numgrids2d) & gridonmacayeal(1:numgrids2d)); %Grids connected to two different types of elements
    121 
    122 %initialize and fill in penalties structure
    123 md.penalties=[];
    124 if ~isnan(bordergrids2d),
    125         penalties=[];
    126         for     i=1:numlayers-1,
    127                 penalties=[penalties; [bordergrids2d bordergrids2d+md.numberofgrids2d*(i)]];
    128         end
    129         md.penalties=penalties;
    130 end
    131167
    132168%flag dead grids (strictly in MacAyeal and not on bed -> not used in diagnostic horiz)
Note: See TracChangeset for help on using the changeset viewer.