Changeset 5379
- Timestamp:
- 08/18/10 15:19:12 (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/m/classes/public/setelementstype.m
r5204 r5379 10 10 % an empty string '' will be considered as an empty domain 11 11 % 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' 12 13 % 13 14 % Usage: … … 16 17 % Example: 17 18 % md=setelementstype(md,'pattyn','Pattyn.exp','macayeal',md.elementoniceshelf,'fill','hutter'); 19 % md=setelementstype(md,'pattyn','Pattyn.exp',fill','hutter','coupling','tiling'); 18 20 19 21 %some checks on list of arguments … … 24 26 if md.counter<3, 25 27 error('only fully parameterized 2d models can be setelementstyped'); 28 end 29 30 %Find_out what kind of coupling to use 31 options=pairoptions(varargin{:}); 32 coupling_method=getfieldvalue(options,'coupling','tiling'); 33 if (~strcmpi(coupling_method,'tiling') & ~strcmpi(coupling_method,'penalties')), 34 error('coupling type can only be: tiling or penalties'); 26 35 end 27 36 … … 98 107 md.elements_type(find(~stokesflag),2)=NoneApproximationEnum(); 99 108 109 %Now take care of the coupling between MacAyeal and Pattyn 110 md.penalties=[]; 111 if 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 126 elseif 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; 147 end 148 100 149 %Create vertices_type 101 150 md.vertices_type=zeros(md.numberofgrids,2); … … 110 159 pos=find(gridonpattyn & gridonmacayeal); 111 160 md.vertices_type(pos,1)=PattynApproximationEnum(); 161 pos=find(gridonmacayealpattyn); 162 md.vertices_type(pos,1)=MacAyealPattynApproximationEnum(); 112 163 pos=find(gridonstokes); 113 164 md.vertices_type(pos,2)=StokesApproximationEnum(); 114 165 pos=find(~gridonstokes); 115 166 md.vertices_type(pos,2)=NoneApproximationEnum(); 116 117 %Create the border grids between Pattyn and MacAyeal and extrude them118 numgrids2d=md.numberofgrids2d;119 numlayers=md.numlayers;120 bordergrids2d=find(gridonpattyn(1:numgrids2d) & gridonmacayeal(1:numgrids2d)); %Grids connected to two different types of elements121 122 %initialize and fill in penalties structure123 md.penalties=[];124 if ~isnan(bordergrids2d),125 penalties=[];126 for i=1:numlayers-1,127 penalties=[penalties; [bordergrids2d bordergrids2d+md.numberofgrids2d*(i)]];128 end129 md.penalties=penalties;130 end131 167 132 168 %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.