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

Last change on this file since 24313 was 24313, checked in by Mathieu Morlighem, 5 years ago

merged trunk-jpl and trunk for revision 24310
File size: 11.5 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','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 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',''));
41FSflag = FlagElements(md,getfieldvalue(options,'FS',''));
42filltype = getfieldvalue(options,'fill','none');
43displayunused(options);
44
45%Flag the elements that have not been flagged as filltype
46if strcmpi(filltype,'SIA'),
47 SIAflag(find(~(SSAflag | HOflag)))=1;
48elseif strcmpi(filltype,'SSA'),
49 SSAflag(find(~(SIAflag | HOflag | FSflag)))=1;
50elseif strcmpi(filltype,'HO'),
51 HOflag(find(~(SIAflag | SSAflag | FSflag)))=1;
52end
53
54%check that each element has at least one flag
55if any(SIAflag+SSAflag+HOflag+L1L2flag+FSflag==0),
56 error('elements type not assigned, supported models are ''SIA'',''SSA'',''HO'' and ''FS''')
57end
58
59%check that each element has only one flag
60if any(SIAflag+SSAflag+HOflag+L1L2flag+FSflag>1),
61 disp('setflowequation warning message: some elements have several types, higher order type is used for them')
62 SIAflag(find(SIAflag & SSAflag))=0;
63 SIAflag(find(SIAflag & HOflag))=0;
64 SSAflag(find(SSAflag & HOflag))=0;
65end
66
67%check that L1L2 is not coupled to any other model for now
68if any(L1L2flag) & any(SIAflag | SSAflag | HOflag | FSflag)
69 error('L1L2 cannot be coupled to any other model');
70end
71
72%Check that no HO or FS for 2d mesh
73if strcmp(domaintype(md.mesh),'2Dhorizontal')
74 if any(FSflag | HOflag)
75 error('FS and HO elements not allowed in 2d mesh, extrude it first')
76 end
77end
78
79%FS can only be used alone for now:
80if any(FSflag) &any(SIAflag),
81 error('FS cannot be used with any other model for now, put FS everywhere')
82end
83
84%Initialize node fields
85nodeonSIA=zeros(md.mesh.numberofvertices,1);
86nodeonSIA(md.mesh.elements(find(SIAflag),:))=1;
87nodeonSSA=zeros(md.mesh.numberofvertices,1);
88nodeonSSA(md.mesh.elements(find(SSAflag),:))=1;
89nodeonHO=zeros(md.mesh.numberofvertices,1);
90nodeonHO(md.mesh.elements(find(HOflag),:))=1;
91nodeonL1L2=zeros(md.mesh.numberofvertices,1);
92nodeonL1L2(md.mesh.elements(find(L1L2flag),:))=1;
93nodeonFS=zeros(md.mesh.numberofvertices,1);
94noneflag=zeros(md.mesh.numberofelements,1);
95
96%First modify FSflag to get rid of elements contrained everywhere (spc + border with HO or SSA)
97if any(FSflag),
98 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
99 fullspcelems=double(sum(fullspcnodes(md.mesh.elements),2)==6); %find all the nodes on the boundary of the domain without icefront
100 FSflag(find(fullspcelems))=0;
101 nodeonFS(md.mesh.elements(find(FSflag),:))=1;
102end
103
104%Then complete with NoneApproximation or the other model used if there is no FS
105if any(FSflag),
106 if any(HOflag), %fill with HO
107 HOflag(~FSflag)=1;
108 nodeonHO(md.mesh.elements(find(HOflag),:))=1;
109 elseif any(SSAflag), %fill with SSA
110 SSAflag(~FSflag)=1;
111 nodeonSSA(md.mesh.elements(find(SSAflag),:))=1;
112 else %fill with none
113 noneflag(find(~FSflag))=1;
114 end
115end
116
117%Now take care of the coupling between SSA and HO
118if strcmpi(coupling_method,'penalties'),
119 md.stressbalance.vertex_pairing=[];
120end
121nodeonSSAHO=zeros(md.mesh.numberofvertices,1);
122nodeonHOFS=zeros(md.mesh.numberofvertices,1);
123nodeonSSAFS=zeros(md.mesh.numberofvertices,1);
124SSAHOflag=zeros(md.mesh.numberofelements,1);
125SSAFSflag=zeros(md.mesh.numberofelements,1);
126HOFSflag=zeros(md.mesh.numberofelements,1);
127if strcmpi(coupling_method,'penalties'),
128 %Create the border nodes between HO and SSA and extrude them
129 numnodes2d=md.mesh.numberofvertices2d;
130 numlayers=md.mesh.numberoflayers;
131 bordernodes2d=find(nodeonHO(1:numnodes2d) & nodeonSSA(1:numnodes2d)); %Nodes connected to two different types of elements
132
133 %initialize and fill in penalties structure
134 if ~isnan(bordernodes2d),
135 penalties=[];
136 for i=1:numlayers-1,
137 penalties=[penalties; [bordernodes2d bordernodes2d+md.mesh.numberofvertices2d*(i)]];
138 end
139 md.stressbalance.vertex_pairing=penalties;
140 end
141elseif strcmpi(coupling_method,'tiling'),
142 if any(SSAflag) & any(HOflag), %coupling SSA HO
143 %Find node at the border
144 nodeonSSAHO(find(nodeonSSA & nodeonHO))=1;
145 %SSA elements in contact with this layer become SSAHO elements
146 matrixelements=ismember(md.mesh.elements,find(nodeonSSAHO));
147 commonelements=sum(matrixelements,2)~=0;
148 commonelements(find(HOflag))=0; %only one layer: the elements previously in SSA
149 SSAflag(find(commonelements))=0; %these elements are now SSAHOelements
150 SSAHOflag(find(commonelements))=1;
151 nodeonSSA(:)=0;
152 nodeonSSA(md.mesh.elements(find(SSAflag),:))=1;
153
154 %rule out elements that don't touch the 2 boundaries
155 pos=find(SSAHOflag);
156 elist=zeros(length(pos),1);
157 elist = elist + any(sum(nodeonSSA(md.mesh.elements(pos,:)),2),2);
158 elist = elist - any(sum(nodeonHO(md.mesh.elements(pos,:)) ,2),2);
159 pos1=find(elist==1);
160 SSAflag(pos(pos1))=1;
161 SSAHOflag(pos(pos1))=0;
162 pos2=find(elist==-1);
163 HOflag(pos(pos2))=1;
164 SSAHOflag(pos(pos2))=0;
165
166 %Recompute nodes associated to these elements
167 nodeonSSA(:)=0;
168 nodeonSSA(md.mesh.elements(find(SSAflag),:))=1;
169 nodeonHO(:)=0;
170 nodeonHO(md.mesh.elements(find(HOflag),:))=1;
171 nodeonSSAHO(:)=0;
172 nodeonSSAHO(md.mesh.elements(find(SSAHOflag),:))=1;
173
174 elseif any(HOflag) & any(FSflag), %coupling HO FS
175 %Find node at the border
176 nodeonHOFS(find(nodeonHO & nodeonFS))=1;
177 %FS elements in contact with this layer become HOFS elements
178 matrixelements=ismember(md.mesh.elements,find(nodeonHOFS));
179 commonelements=sum(matrixelements,2)~=0;
180 commonelements(find(HOflag))=0; %only one layer: the elements previously in SSA
181 FSflag(find(commonelements))=0; %these elements are now SSAHOelements
182 HOFSflag(find(commonelements))=1;
183 nodeonFS=zeros(md.mesh.numberofvertices,1);
184 nodeonFS(md.mesh.elements(find(FSflag),:))=1;
185
186 %rule out elements that don't touch the 2 boundaries
187 pos=find(HOFSflag);
188 elist=zeros(length(pos),1);
189 elist = elist + any(sum(nodeonFS(md.mesh.elements(pos,:)),2),2);
190 elist = elist - any(sum(nodeonHO(md.mesh.elements(pos,:)),2),2);
191 pos1=find(elist==1);
192 FSflag(pos(pos1))=1;
193 HOFSflag(pos(pos1))=0;
194 pos2=find(elist==-1);
195 HOflag(pos(pos2))=1;
196 HOFSflag(pos(pos2))=0;
197
198 %Recompute nodes associated to these elements
199 nodeonFS(:)=0;
200 nodeonFS(md.mesh.elements(find(FSflag),:))=1;
201 nodeonHO(:)=0;
202 nodeonHO(md.mesh.elements(find(HOflag),:))=1;
203 nodeonHOFS(:)=0;
204 nodeonHOFS(md.mesh.elements(find(HOFSflag),:))=1;
205
206 elseif any(FSflag) & any(SSAflag),
207 %Find node at the border
208 nodeonSSAFS(find(nodeonSSA & nodeonFS))=1;
209 %FS elements in contact with this layer become SSAFS elements
210 matrixelements=ismember(md.mesh.elements,find(nodeonSSAFS));
211 commonelements=sum(matrixelements,2)~=0;
212 commonelements(find(SSAflag))=0; %only one layer: the elements previously in SSA
213 FSflag(find(commonelements))=0; %these elements are now SSASSAelements
214 SSAFSflag(find(commonelements))=1;
215 nodeonFS=zeros(md.mesh.numberofvertices,1);
216 nodeonFS(md.mesh.elements(find(FSflag),:))=1;
217
218 %rule out elements that don't touch the 2 boundaries
219 pos=find(SSAFSflag);
220 elist=zeros(length(pos),1);
221 elist = elist + any(sum(nodeonSSA(md.mesh.elements(pos,:)),2),2);
222 elist = elist - any(sum(nodeonFS(md.mesh.elements(pos,:)) ,2),2);
223 pos1=find(elist==1);
224 SSAflag(pos(pos1))=1;
225 SSAFSflag(pos(pos1))=0;
226 pos2=find(elist==-1);
227 FSflag(pos(pos2))=1;
228 SSAFSflag(pos(pos2))=0;
229
230 %Recompute nodes associated to these elements
231 nodeonSSA(:)=0;
232 nodeonSSA(md.mesh.elements(find(SSAflag),:))=1;
233 nodeonFS(:)=0;
234 nodeonFS(md.mesh.elements(find(FSflag),:))=1;
235 nodeonSSAFS(:)=0;
236 nodeonSSAFS(md.mesh.elements(find(SSAFSflag),:))=1;
237
238 elseif any(FSflag) & any(SIAflag),
239 error('type of coupling not supported yet');
240 end
241end
242
243%Create element equations
244md.flowequation.element_equation=zeros(md.mesh.numberofelements,1);
245md.flowequation.element_equation(find(noneflag))=0;
246md.flowequation.element_equation(find(SIAflag))=1;
247md.flowequation.element_equation(find(SSAflag))=2;
248md.flowequation.element_equation(find(L1L2flag))=3;
249md.flowequation.element_equation(find(HOflag))=4;
250md.flowequation.element_equation(find(FSflag))=5;
251md.flowequation.element_equation(find(SSAHOflag))=6;
252md.flowequation.element_equation(find(SSAFSflag))=7;
253md.flowequation.element_equation(find(HOFSflag))=8;
254
255%border
256md.flowequation.borderHO=nodeonHO;
257md.flowequation.borderSSA=nodeonSSA;
258md.flowequation.borderFS=nodeonFS;
259
260%Create vertices_type
261md.flowequation.vertex_equation=zeros(md.mesh.numberofvertices,1);
262pos=find(nodeonSSA);
263md.flowequation.vertex_equation(pos)=2;
264pos=find(nodeonL1L2);
265md.flowequation.vertex_equation(pos)=3;
266pos=find(nodeonHO);
267md.flowequation.vertex_equation(pos)=4;
268pos=find(nodeonFS);
269md.flowequation.vertex_equation(pos)=5;
270%DO SIA LAST! Otherwise spcs might not be set up correctly (SIA should have priority)
271pos=find(nodeonSIA);
272md.flowequation.vertex_equation(pos)=1;
273if any(FSflag),
274 pos=find(~nodeonFS);
275 if(~any(HOflag) & ~any(SSAflag)),
276 md.flowequation.vertex_equation(pos)=0;
277 end
278end
279pos=find(nodeonSSAHO);
280md.flowequation.vertex_equation(pos)=6;
281pos=find(nodeonHOFS);
282md.flowequation.vertex_equation(pos)=7;
283pos=find(nodeonSSAFS);
284md.flowequation.vertex_equation(pos)=8;
285
286%figure out solution types
287md.flowequation.isSIA = double(any(md.flowequation.element_equation == 1));
288md.flowequation.isSSA = double(any(md.flowequation.element_equation == 2));
289md.flowequation.isL1L2 = double(any(md.flowequation.element_equation == 3));
290md.flowequation.isHO = double(any(md.flowequation.element_equation == 4));
291md.flowequation.isFS = double(any(md.flowequation.element_equation == 5));
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.