source: issm/trunk/src/m/parameterization/setflowequation.m@ 26744

Last change on this file since 26744 was 26744, checked in by Mathieu Morlighem, 4 years ago

merged trunk-jpl and trunk for revision 26742
File size: 12.0 KB
Line 
1function md=setflowequation(md,varargin)
2%SETFLOWEQUATION - associate a solution type to each element
3%
4% This routine works like plotmodel: it works with an even number of inputs
5% 'SIA','SSA','L1L2','MLHO','HO','FS' and 'fill' are the possible options
6% that must be followed by the corresponding exp file or flags list
7% It can either be a domain file (argus type, .exp extension), or an array of element flags.
8% If user wants every element outside the domain to be
9% setflowequationd, add '~' to the name of the domain file (ex: '~HO.exp');
10% an empty string '' will be considered as an empty domain
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'
13% NB: L1L2 and MLHO cannot currently be coupled to any other ice flow model
14%
15% Usage:
16% md=setflowequation(md,varargin)
17%
18% Example:
19% md=setflowequation(md,'HO','HO.exp',fill','SIA','coupling','tiling');
20
21%some checks on list of arguments
22if ((nargin<2) | (nargout~=1)),
23 error('setflowequation error message');
24end
25
26%Process options
27options=pairoptions(varargin{:});
28options=deleteduplicates(options,1);
29
30%Find_out what kind of coupling to use
31coupling_method=getfieldvalue(options,'coupling','tiling');
32if (~strcmpi(coupling_method,'tiling') & ~strcmpi(coupling_method,'penalties')),
33 error('coupling type can only be: tiling or penalties');
34end
35
36%recover elements distribution
37SIAflag = FlagElements(md,getfieldvalue(options,'SIA',''));
38SSAflag = FlagElements(md,getfieldvalue(options,'SSA',''));
39HOflag = FlagElements(md,getfieldvalue(options,'HO',''));
40L1L2flag = FlagElements(md,getfieldvalue(options,'L1L2',''));
41MLHOflag = FlagElements(md,getfieldvalue(options,'MLHO',''));
42FSflag = FlagElements(md,getfieldvalue(options,'FS',''));
43filltype = getfieldvalue(options,'fill','none');
44displayunused(options);
45
46%Flag the elements that have not been flagged as filltype
47if strcmpi(filltype,'SIA'),
48 SIAflag(find(~(SSAflag | HOflag)))=1;
49elseif strcmpi(filltype,'SSA'),
50 SSAflag(find(~(SIAflag | HOflag | FSflag)))=1;
51elseif strcmpi(filltype,'HO'),
52 HOflag(find(~(SIAflag | SSAflag | FSflag)))=1;
53end
54
55%check that each element has at least one flag
56if any(SIAflag+SSAflag+HOflag+L1L2flag+MLHOflag+FSflag==0),
57 error('elements type not assigned, supported models are ''SIA'',''SSA'',''HO'' and ''FS''')
58end
59
60%check that each element has only one flag
61if any(SIAflag+SSAflag+HOflag+L1L2flag+MLHOflag+FSflag>1),
62 disp('setflowequation.m: Warning: some elements have several types, higher order type is used for them')
63 SIAflag(find(SIAflag & SSAflag))=0;
64 SIAflag(find(SIAflag & HOflag))=0;
65 SSAflag(find(SSAflag & HOflag))=0;
66end
67
68%check that L1L2 is not coupled to any other model for now
69if any(L1L2flag) & any(SIAflag | SSAflag | HOflag | FSflag)
70 error('L1L2 cannot be coupled to any other model');
71end
72if any(MLHOflag) & any(SIAflag | SSAflag | HOflag | FSflag)
73 error('MLHO cannot be coupled to any other model');
74end
75
76%Check that no HO or FS for 2d mesh
77if strcmp(domaintype(md.mesh),'2Dhorizontal')
78 if any(FSflag | HOflag)
79 error('FS and HO elements not allowed in 2d mesh, extrude it first')
80 end
81end
82
83%FS can only be used alone for now:
84if any(FSflag) &any(SIAflag),
85 error('FS cannot be used with any other model for now, put FS everywhere')
86end
87
88%Initialize node fields
89nodeonSIA=zeros(md.mesh.numberofvertices,1); nodeonSIA(md.mesh.elements(find(SIAflag),:))=1;
90nodeonSSA=zeros(md.mesh.numberofvertices,1); nodeonSSA(md.mesh.elements(find(SSAflag),:))=1;
91nodeonHO=zeros(md.mesh.numberofvertices,1); nodeonHO(md.mesh.elements(find(HOflag),:))=1;
92nodeonL1L2=zeros(md.mesh.numberofvertices,1); nodeonL1L2(md.mesh.elements(find(L1L2flag),:))=1;
93nodeonMLHO=zeros(md.mesh.numberofvertices,1); nodeonMLHO(md.mesh.elements(find(MLHOflag),:))=1;
94nodeonFS=zeros(md.mesh.numberofvertices,1);
95noneflag=zeros(md.mesh.numberofelements,1);
96
97%First modify FSflag to get rid of elements contrained everywhere (spc + border with HO or SSA)
98if any(FSflag),
99 fullspcnodes=double((~isnan(md.stressbalance.spcvx)+~isnan(md.stressbalance.spcvy)+~isnan(md.stressbalance.spcvz))==3 | (nodeonHO & nodeonFS)); %find all the nodes on the boundary of the domain without icefront
100 fullspcelems=double(sum(fullspcnodes(md.mesh.elements),2)==6); %find all the nodes on the boundary of the domain without icefront
101 FSflag(find(fullspcelems))=0;
102 nodeonFS(md.mesh.elements(find(FSflag),:))=1;
103end
104
105%Then complete with NoneApproximation or the other model used if there is no FS
106if any(FSflag),
107 if any(HOflag), %fill with HO
108 HOflag(~FSflag)=1;
109 nodeonHO(md.mesh.elements(find(HOflag),:))=1;
110 elseif any(SSAflag), %fill with SSA
111 SSAflag(~FSflag)=1;
112 nodeonSSA(md.mesh.elements(find(SSAflag),:))=1;
113 else %fill with none
114 noneflag(find(~FSflag))=1;
115 end
116end
117
118%Now take care of the coupling between SSA and HO
119if strcmpi(coupling_method,'penalties'),
120 md.stressbalance.vertex_pairing=[];
121end
122nodeonSSAHO=zeros(md.mesh.numberofvertices,1);
123nodeonHOFS=zeros(md.mesh.numberofvertices,1);
124nodeonSSAFS=zeros(md.mesh.numberofvertices,1);
125SSAHOflag=zeros(md.mesh.numberofelements,1);
126SSAFSflag=zeros(md.mesh.numberofelements,1);
127HOFSflag=zeros(md.mesh.numberofelements,1);
128if strcmpi(coupling_method,'penalties'),
129 %Create the border nodes between HO and SSA and extrude them
130 numnodes2d=md.mesh.numberofvertices2d;
131 numlayers=md.mesh.numberoflayers;
132 bordernodes2d=find(nodeonHO(1:numnodes2d) & nodeonSSA(1:numnodes2d)); %Nodes connected to two different types of elements
133
134 %initialize and fill in penalties structure
135 if ~isnan(bordernodes2d),
136 penalties=[];
137 for i=1:numlayers-1,
138 penalties=[penalties; [bordernodes2d bordernodes2d+md.mesh.numberofvertices2d*(i)]];
139 end
140 md.stressbalance.vertex_pairing=penalties;
141 end
142elseif strcmpi(coupling_method,'tiling'),
143 if any(SSAflag) & any(HOflag), %coupling SSA HO
144 %Find node at the border
145 nodeonSSAHO(find(nodeonSSA & nodeonHO))=1;
146 %SSA elements in contact with this layer become SSAHO elements
147 matrixelements=ismember(md.mesh.elements,find(nodeonSSAHO));
148 commonelements=sum(matrixelements,2)~=0;
149 commonelements(find(HOflag))=0; %only one layer: the elements previously in SSA
150 SSAflag(find(commonelements))=0; %these elements are now SSAHOelements
151 SSAHOflag(find(commonelements))=1;
152 nodeonSSA(:)=0;
153 nodeonSSA(md.mesh.elements(find(SSAflag),:))=1;
154
155 %rule out elements that don't touch the 2 boundaries
156 pos=find(SSAHOflag);
157 elist=zeros(length(pos),1);
158 elist = elist + any(sum(nodeonSSA(md.mesh.elements(pos,:)),2),2);
159 elist = elist - any(sum(nodeonHO(md.mesh.elements(pos,:)) ,2),2);
160 pos1=find(elist==1);
161 SSAflag(pos(pos1))=1;
162 SSAHOflag(pos(pos1))=0;
163 pos2=find(elist==-1);
164 HOflag(pos(pos2))=1;
165 SSAHOflag(pos(pos2))=0;
166
167 %Recompute nodes associated to these elements
168 nodeonSSA(:)=0;
169 nodeonSSA(md.mesh.elements(find(SSAflag),:))=1;
170 nodeonHO(:)=0;
171 nodeonHO(md.mesh.elements(find(HOflag),:))=1;
172 nodeonSSAHO(:)=0;
173 nodeonSSAHO(md.mesh.elements(find(SSAHOflag),:))=1;
174
175 elseif any(HOflag) & any(FSflag), %coupling HO FS
176 %Find node at the border
177 nodeonHOFS(find(nodeonHO & nodeonFS))=1;
178 %FS elements in contact with this layer become HOFS elements
179 matrixelements=ismember(md.mesh.elements,find(nodeonHOFS));
180 commonelements=sum(matrixelements,2)~=0;
181 commonelements(find(HOflag))=0; %only one layer: the elements previously in SSA
182 FSflag(find(commonelements))=0; %these elements are now SSAHOelements
183 HOFSflag(find(commonelements))=1;
184 nodeonFS=zeros(md.mesh.numberofvertices,1);
185 nodeonFS(md.mesh.elements(find(FSflag),:))=1;
186
187 %rule out elements that don't touch the 2 boundaries
188 pos=find(HOFSflag);
189 elist=zeros(length(pos),1);
190 elist = elist + any(sum(nodeonFS(md.mesh.elements(pos,:)),2),2);
191 elist = elist - any(sum(nodeonHO(md.mesh.elements(pos,:)),2),2);
192 pos1=find(elist==1);
193 FSflag(pos(pos1))=1;
194 HOFSflag(pos(pos1))=0;
195 pos2=find(elist==-1);
196 HOflag(pos(pos2))=1;
197 HOFSflag(pos(pos2))=0;
198
199 %Recompute nodes associated to these elements
200 nodeonFS(:)=0;
201 nodeonFS(md.mesh.elements(find(FSflag),:))=1;
202 nodeonHO(:)=0;
203 nodeonHO(md.mesh.elements(find(HOflag),:))=1;
204 nodeonHOFS(:)=0;
205 nodeonHOFS(md.mesh.elements(find(HOFSflag),:))=1;
206
207 elseif any(FSflag) & any(SSAflag),
208 %Find node at the border
209 nodeonSSAFS(find(nodeonSSA & nodeonFS))=1;
210 %FS elements in contact with this layer become SSAFS elements
211 matrixelements=ismember(md.mesh.elements,find(nodeonSSAFS));
212 commonelements=sum(matrixelements,2)~=0;
213 commonelements(find(SSAflag))=0; %only one layer: the elements previously in SSA
214 FSflag(find(commonelements))=0; %these elements are now SSASSAelements
215 SSAFSflag(find(commonelements))=1;
216 nodeonFS=zeros(md.mesh.numberofvertices,1);
217 nodeonFS(md.mesh.elements(find(FSflag),:))=1;
218
219 %rule out elements that don't touch the 2 boundaries
220 pos=find(SSAFSflag);
221 elist=zeros(length(pos),1);
222 elist = elist + any(sum(nodeonSSA(md.mesh.elements(pos,:)),2),2);
223 elist = elist - any(sum(nodeonFS(md.mesh.elements(pos,:)) ,2),2);
224 pos1=find(elist==1);
225 SSAflag(pos(pos1))=1;
226 SSAFSflag(pos(pos1))=0;
227 pos2=find(elist==-1);
228 FSflag(pos(pos2))=1;
229 SSAFSflag(pos(pos2))=0;
230
231 %Recompute nodes associated to these elements
232 nodeonSSA(:)=0;
233 nodeonSSA(md.mesh.elements(find(SSAflag),:))=1;
234 nodeonFS(:)=0;
235 nodeonFS(md.mesh.elements(find(FSflag),:))=1;
236 nodeonSSAFS(:)=0;
237 nodeonSSAFS(md.mesh.elements(find(SSAFSflag),:))=1;
238
239 elseif any(FSflag) & any(SIAflag),
240 error('type of coupling not supported yet');
241 end
242end
243
244%Create element equations
245md.flowequation.element_equation=zeros(md.mesh.numberofelements,1);
246md.flowequation.element_equation(find(noneflag))=0;
247md.flowequation.element_equation(find(SIAflag))=1;
248md.flowequation.element_equation(find(SSAflag))=2;
249md.flowequation.element_equation(find(L1L2flag))=3;
250md.flowequation.element_equation(find(MLHOflag))=4;
251md.flowequation.element_equation(find(HOflag))=5;
252md.flowequation.element_equation(find(FSflag))=6;
253md.flowequation.element_equation(find(SSAHOflag))=7;
254md.flowequation.element_equation(find(SSAFSflag))=8;
255md.flowequation.element_equation(find(HOFSflag))=9;
256
257%border
258md.flowequation.borderHO=nodeonHO;
259md.flowequation.borderSSA=nodeonSSA;
260md.flowequation.borderFS=nodeonFS;
261
262%Create vertices_type
263md.flowequation.vertex_equation=zeros(md.mesh.numberofvertices,1);
264pos=find(nodeonSSA); md.flowequation.vertex_equation(pos)=2;
265pos=find(nodeonL1L2); md.flowequation.vertex_equation(pos)=3;
266pos=find(nodeonMLHO); md.flowequation.vertex_equation(pos)=4;
267pos=find(nodeonHO); md.flowequation.vertex_equation(pos)=5;
268pos=find(nodeonFS); md.flowequation.vertex_equation(pos)=6;
269%DO SIA LAST! Otherwise spcs might not be set up correctly (SIA should have priority)
270pos=find(nodeonSIA);
271md.flowequation.vertex_equation(pos)=1;
272if any(FSflag),
273 pos=find(~nodeonFS);
274 if(~any(HOflag) & ~any(SSAflag)),
275 md.flowequation.vertex_equation(pos)=0;
276 end
277end
278pos=find(nodeonSSAHO);
279md.flowequation.vertex_equation(pos)=7;
280pos=find(nodeonHOFS);
281md.flowequation.vertex_equation(pos)=8;
282pos=find(nodeonSSAFS);
283md.flowequation.vertex_equation(pos)=9;
284
285%figure out solution types
286md.flowequation.isSIA = double(any(md.flowequation.element_equation == 1));
287md.flowequation.isSSA = double(any(md.flowequation.element_equation == 2));
288md.flowequation.isL1L2 = double(any(md.flowequation.element_equation == 3));
289md.flowequation.isMLHO = double(any(md.flowequation.element_equation == 4));
290md.flowequation.isHO = double(any(md.flowequation.element_equation == 5));
291md.flowequation.isFS = double(any(md.flowequation.element_equation == 6));
292
293return
294
295%Check that tiling can work:
296if any(md.flowequation.borderSSA) & any(md.flowequation.borderHO) & any(md.flowequation.borderHO + md.flowequation.borderSSA ~=1),
297 error('error coupling domain too irregular');
298end
299if any(md.flowequation.borderSSA) & any(md.flowequation.borderFS) & any(md.flowequation.borderFS + md.flowequation.borderSSA ~=1),
300 error('error coupling domain too irregular');
301end
302if any(md.flowequation.borderFS) & any(md.flowequation.borderHO) & any(md.flowequation.borderHO + md.flowequation.borderFS~=1),
303 error('error coupling domain too irregular');
304end
Note: See TracBrowser for help on using the repository browser.