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

Last change on this file since 20500 was 20500, checked in by Mathieu Morlighem, 10 years ago

merged trunk-jpl and trunk for revision 20497
File size: 11.4 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
118md.stressbalance.vertex_pairing=[];
119nodeonSSAHO=zeros(md.mesh.numberofvertices,1);
120nodeonHOFS=zeros(md.mesh.numberofvertices,1);
121nodeonSSAFS=zeros(md.mesh.numberofvertices,1);
122SSAHOflag=zeros(md.mesh.numberofelements,1);
123SSAFSflag=zeros(md.mesh.numberofelements,1);
124HOFSflag=zeros(md.mesh.numberofelements,1);
125if strcmpi(coupling_method,'penalties'),
126 %Create the border nodes between HO and SSA and extrude them
127 numnodes2d=md.mesh.numberofvertices2d;
128 numlayers=md.mesh.numberoflayers;
129 bordernodes2d=find(nodeonHO(1:numnodes2d) & nodeonSSA(1:numnodes2d)); %Nodes connected to two different types of elements
130
131 %initialize and fill in penalties structure
132 if ~isnan(bordernodes2d),
133 penalties=[];
134 for i=1:numlayers-1,
135 penalties=[penalties; [bordernodes2d bordernodes2d+md.mesh.numberofvertices2d*(i)]];
136 end
137 md.stressbalance.vertex_pairing=penalties;
138 end
139elseif strcmpi(coupling_method,'tiling'),
140 if any(SSAflag) & any(HOflag), %coupling SSA HO
141 %Find node at the border
142 nodeonSSAHO(find(nodeonSSA & nodeonHO))=1;
143 %SSA elements in contact with this layer become SSAHO elements
144 matrixelements=ismember(md.mesh.elements,find(nodeonSSAHO));
145 commonelements=sum(matrixelements,2)~=0;
146 commonelements(find(HOflag))=0; %only one layer: the elements previously in SSA
147 SSAflag(find(commonelements))=0; %these elements are now SSAHOelements
148 SSAHOflag(find(commonelements))=1;
149 nodeonSSA(:)=0;
150 nodeonSSA(md.mesh.elements(find(SSAflag),:))=1;
151
152 %rule out elements that don't touch the 2 boundaries
153 pos=find(SSAHOflag);
154 elist=zeros(length(pos),1);
155 elist = elist + any(sum(nodeonSSA(md.mesh.elements(pos,:)),2),2);
156 elist = elist - any(sum(nodeonHO(md.mesh.elements(pos,:)) ,2),2);
157 pos1=find(elist==1);
158 SSAflag(pos(pos1))=1;
159 SSAHOflag(pos(pos1))=0;
160 pos2=find(elist==-1);
161 HOflag(pos(pos2))=1;
162 SSAHOflag(pos(pos2))=0;
163
164 %Recompute nodes associated to these elements
165 nodeonSSA(:)=0;
166 nodeonSSA(md.mesh.elements(find(SSAflag),:))=1;
167 nodeonHO(:)=0;
168 nodeonHO(md.mesh.elements(find(HOflag),:))=1;
169 nodeonSSAHO(:)=0;
170 nodeonSSAHO(md.mesh.elements(find(SSAHOflag),:))=1;
171
172 elseif any(HOflag) & any(FSflag), %coupling HO FS
173 %Find node at the border
174 nodeonHOFS(find(nodeonHO & nodeonFS))=1;
175 %FS elements in contact with this layer become HOFS elements
176 matrixelements=ismember(md.mesh.elements,find(nodeonHOFS));
177 commonelements=sum(matrixelements,2)~=0;
178 commonelements(find(HOflag))=0; %only one layer: the elements previously in SSA
179 FSflag(find(commonelements))=0; %these elements are now SSAHOelements
180 HOFSflag(find(commonelements))=1;
181 nodeonFS=zeros(md.mesh.numberofvertices,1);
182 nodeonFS(md.mesh.elements(find(FSflag),:))=1;
183
184 %rule out elements that don't touch the 2 boundaries
185 pos=find(HOFSflag);
186 elist=zeros(length(pos),1);
187 elist = elist + any(sum(nodeonFS(md.mesh.elements(pos,:)),2),2);
188 elist = elist - any(sum(nodeonHO(md.mesh.elements(pos,:)),2),2);
189 pos1=find(elist==1);
190 FSflag(pos(pos1))=1;
191 HOFSflag(pos(pos1))=0;
192 pos2=find(elist==-1);
193 HOflag(pos(pos2))=1;
194 HOFSflag(pos(pos2))=0;
195
196 %Recompute nodes associated to these elements
197 nodeonFS(:)=0;
198 nodeonFS(md.mesh.elements(find(FSflag),:))=1;
199 nodeonHO(:)=0;
200 nodeonHO(md.mesh.elements(find(HOflag),:))=1;
201 nodeonHOFS(:)=0;
202 nodeonHOFS(md.mesh.elements(find(HOFSflag),:))=1;
203
204 elseif any(FSflag) & any(SSAflag),
205 %Find node at the border
206 nodeonSSAFS(find(nodeonSSA & nodeonFS))=1;
207 %FS elements in contact with this layer become SSAFS elements
208 matrixelements=ismember(md.mesh.elements,find(nodeonSSAFS));
209 commonelements=sum(matrixelements,2)~=0;
210 commonelements(find(SSAflag))=0; %only one layer: the elements previously in SSA
211 FSflag(find(commonelements))=0; %these elements are now SSASSAelements
212 SSAFSflag(find(commonelements))=1;
213 nodeonFS=zeros(md.mesh.numberofvertices,1);
214 nodeonFS(md.mesh.elements(find(FSflag),:))=1;
215
216 %rule out elements that don't touch the 2 boundaries
217 pos=find(SSAFSflag);
218 elist=zeros(length(pos),1);
219 elist = elist + any(sum(nodeonSSA(md.mesh.elements(pos,:)),2),2);
220 elist = elist - any(sum(nodeonFS(md.mesh.elements(pos,:)) ,2),2);
221 pos1=find(elist==1);
222 SSAflag(pos(pos1))=1;
223 SSAFSflag(pos(pos1))=0;
224 pos2=find(elist==-1);
225 FSflag(pos(pos2))=1;
226 SSAFSflag(pos(pos2))=0;
227
228 %Recompute nodes associated to these elements
229 nodeonSSA(:)=0;
230 nodeonSSA(md.mesh.elements(find(SSAflag),:))=1;
231 nodeonFS(:)=0;
232 nodeonFS(md.mesh.elements(find(FSflag),:))=1;
233 nodeonSSAFS(:)=0;
234 nodeonSSAFS(md.mesh.elements(find(SSAFSflag),:))=1;
235
236 elseif any(FSflag) & any(SIAflag),
237 error('type of coupling not supported yet');
238 end
239end
240
241%Create element equations
242md.flowequation.element_equation=zeros(md.mesh.numberofelements,1);
243md.flowequation.element_equation(find(noneflag))=0;
244md.flowequation.element_equation(find(SIAflag))=1;
245md.flowequation.element_equation(find(SSAflag))=2;
246md.flowequation.element_equation(find(L1L2flag))=3;
247md.flowequation.element_equation(find(HOflag))=4;
248md.flowequation.element_equation(find(FSflag))=5;
249md.flowequation.element_equation(find(SSAHOflag))=6;
250md.flowequation.element_equation(find(SSAFSflag))=7;
251md.flowequation.element_equation(find(HOFSflag))=8;
252
253%border
254md.flowequation.borderHO=nodeonHO;
255md.flowequation.borderSSA=nodeonSSA;
256md.flowequation.borderFS=nodeonFS;
257
258%Create vertices_type
259md.flowequation.vertex_equation=zeros(md.mesh.numberofvertices,1);
260pos=find(nodeonSSA);
261md.flowequation.vertex_equation(pos)=2;
262pos=find(nodeonL1L2);
263md.flowequation.vertex_equation(pos)=3;
264pos=find(nodeonHO);
265md.flowequation.vertex_equation(pos)=4;
266pos=find(nodeonFS);
267md.flowequation.vertex_equation(pos)=5;
268%DO SIA LAST! Otherwise spcs might not be set up correctly (SIA should have priority)
269pos=find(nodeonSIA);
270md.flowequation.vertex_equation(pos)=1;
271if any(FSflag),
272 pos=find(~nodeonFS);
273 if(~any(HOflag) & ~any(SSAflag)),
274 md.flowequation.vertex_equation(pos)=0;
275 end
276end
277pos=find(nodeonSSAHO);
278md.flowequation.vertex_equation(pos)=6;
279pos=find(nodeonHOFS);
280md.flowequation.vertex_equation(pos)=7;
281pos=find(nodeonSSAFS);
282md.flowequation.vertex_equation(pos)=8;
283
284%figure out solution types
285md.flowequation.isSIA = double(any(md.flowequation.element_equation == 1));
286md.flowequation.isSSA = double(any(md.flowequation.element_equation == 2));
287md.flowequation.isL1L2 = double(any(md.flowequation.element_equation == 3));
288md.flowequation.isHO = double(any(md.flowequation.element_equation == 4));
289md.flowequation.isFS = double(any(md.flowequation.element_equation == 5));
290
291return
292
293%Check that tiling can work:
294if any(md.flowequation.borderSSA) & any(md.flowequation.borderHO) & any(md.flowequation.borderHO + md.flowequation.borderSSA ~=1),
295 error('error coupling domain too irregular');
296end
297if any(md.flowequation.borderSSA) & any(md.flowequation.borderFS) & any(md.flowequation.borderFS + md.flowequation.borderSSA ~=1),
298 error('error coupling domain too irregular');
299end
300if any(md.flowequation.borderFS) & any(md.flowequation.borderHO) & any(md.flowequation.borderHO + md.flowequation.borderFS~=1),
301 error('error coupling domain too irregular');
302end
Note: See TracBrowser for help on using the repository browser.