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

Last change on this file since 13395 was 13395, checked in by Mathieu Morlighem, 12 years ago

merged trunk-jpl and trunk for revision 13393

File size: 12.6 KB
RevLine 
[9664]1function md=setflowequation(md,varargin)
[1]2%SETELEMENTSTYPE - associate a solution type to each element
3%
4% This routine works like plotmodel: it works with an even number of inputs
[13027]5% 'hutter','macayeal','l1l2','pattyn','stokes' and 'fill' are the possible options
[1]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
[9664]9% setflowequationd, add '~' to the name of the domain file (ex: '~Pattyn.exp');
[1]10% an empty string '' will be considered as an empty domain
11% a string 'all' will be considered as the entire domain
[5379]12% You can specify the type of coupling, 'penalties' or 'tiling', to use with the input 'coupling'
[13048]13% NB: l1l2 cannot currently be coupled to any other ice flow model
[1]14%
15% Usage:
[9664]16% md=setflowequation(md,varargin)
[1]17%
18% Example:
[9664]19% md=setflowequation(md,'pattyn','Pattyn.exp','macayeal',md.mask.elementonfloatingice,'fill','hutter');
20% md=setflowequation(md,'pattyn','Pattyn.exp',fill','hutter','coupling','tiling');
[1]21
22%some checks on list of arguments
23if ((nargin<2) | (nargout~=1)),
[9664]24 error('setflowequation error message');
[1]25end
26
[13028]27%Process options
28options=pairoptions(varargin{:});
29options=deleteduplicates(options,1);
30
[5379]31%Find_out what kind of coupling to use
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');
35end
36
[13028]37%recover elements distribution
38hutterflag = FlagElements(md,getfieldvalue(options,'hutter',''));
39macayealflag = FlagElements(md,getfieldvalue(options,'macayeal',''));
40pattynflag = FlagElements(md,getfieldvalue(options,'pattyn',''));
41l1l2flag = FlagElements(md,getfieldvalue(options,'l1l2',''));
42stokesflag = FlagElements(md,getfieldvalue(options,'stokes',''));
43filltype = getfieldvalue(options,'fill','none');
[1]44
[12888]45%Flag the elements that have not been flagged as filltype
[1]46if strcmpi(filltype,'hutter'),
[12949]47 hutterflag(find(~(macayealflag | pattynflag)))=1;
[1]48elseif strcmpi(filltype,'macayeal'),
[12950]49 macayealflag(find(~(hutterflag | pattynflag | stokesflag)))=1;
[1]50elseif strcmpi(filltype,'pattyn'),
[12950]51 pattynflag(find(~(hutterflag | macayealflag | stokesflag)))=1;
[1]52end
53
54%check that each element has at least one flag
[13048]55if any(hutterflag+macayealflag+pattynflag+l1l2flag+stokesflag==0),
56 error('elements type not assigned, must be specified')
[1]57end
58
59%check that each element has only one flag
[13048]60if any(hutterflag+macayealflag+pattynflag+l1l2flag+stokesflag>1),
[9664]61 disp('setflowequation warning message: some elements have several types, higher order type is used for them')
[1]62 hutterflag(find(hutterflag & macayealflag))=0;
[5613]63 hutterflag(find(hutterflag & pattynflag))=0;
[1]64 macayealflag(find(macayealflag & pattynflag))=0;
65end
66
[13048]67%check that l1l2 is not coupled to any other model for now
68if any(l1l2flag) & any(hutterflag | macayealflag | pattynflag | stokesflag)
69 error('l1l2 cannot be coupled to any other model');
70end
71
72%Check that no l1l2 or pattyn or stokes for 2d mesh
[9719]73if (md.mesh.dimension==2),
[13048]74 if any(l1l2flag | stokesflag | pattynflag)
75 error('stokes and pattyn elements not allowed in 2d mesh, extrude it first')
[1]76 end
77end
78
[5433]79%Stokes can only be used alone for now:
[6705]80if any(stokesflag) &any(hutterflag),
[13048]81 error('stokes cannot be used with any other model for now, put stokes everywhere')
[5433]82end
[1]83
[11004]84%Initialize node fields
[9725]85nodeonhutter=zeros(md.mesh.numberofvertices,1);
[9733]86nodeonhutter(md.mesh.elements(find(hutterflag),:))=1;
[9725]87nodeonmacayeal=zeros(md.mesh.numberofvertices,1);
[9733]88nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1;
[9725]89nodeonpattyn=zeros(md.mesh.numberofvertices,1);
[9733]90nodeonpattyn(md.mesh.elements(find(pattynflag),:))=1;
[13048]91nodeonl1l2=zeros(md.mesh.numberofvertices,1);
92nodeonl1l2(md.mesh.elements(find(l1l2flag),:))=1;
[11004]93nodeonstokes=zeros(md.mesh.numberofvertices,1);
94noneflag=zeros(md.mesh.numberofelements,1);
[1]95
[5613]96%First modify stokesflag to get rid of elements contrained everywhere (spc + border with pattyn or macayeal)
97if any(stokesflag),
[9679]98 fullspcnodes=double((~isnan(md.diagnostic.spcvx)+~isnan(md.diagnostic.spcvy)+~isnan(md.diagnostic.spcvz))==3 | (nodeonpattyn & nodeonstokes)); %find all the nodes on the boundary of the domain without icefront
[9733]99 fullspcelems=double(sum(fullspcnodes(md.mesh.elements),2)==6); %find all the nodes on the boundary of the domain without icefront
[5613]100 stokesflag(find(fullspcelems))=0;
[11004]101 nodeonstokes(md.mesh.elements(find(stokesflag),:))=1;
[5613]102end
[1]103
[5613]104%Then complete with NoneApproximation or the other model used if there is no stokes
[5433]105if any(stokesflag),
[5613]106 if any(pattynflag), %fill with pattyn
107 pattynflag(~stokesflag)=1;
[9733]108 nodeonpattyn(md.mesh.elements(find(pattynflag),:))=1;
[6712]109 elseif any(macayealflag), %fill with macayeal
110 macayealflag(~stokesflag)=1;
[9733]111 nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1;
[5613]112 else %fill with none
[11003]113 noneflag(find(~stokesflag))=1;
[5613]114 end
[5433]115end
[1]116
[5379]117%Now take care of the coupling between MacAyeal and Pattyn
[9679]118md.diagnostic.vertex_pairing=[];
[9725]119nodeonmacayealpattyn=zeros(md.mesh.numberofvertices,1);
120nodeonpattynstokes=zeros(md.mesh.numberofvertices,1);
121nodeonmacayealstokes=zeros(md.mesh.numberofvertices,1);
[11003]122macayealpattynflag=zeros(md.mesh.numberofelements,1);
123macayealstokesflag=zeros(md.mesh.numberofelements,1);
124pattynstokesflag=zeros(md.mesh.numberofelements,1);
[5379]125if strcmpi(coupling_method,'penalties'),
[8298]126 %Create the border nodes between Pattyn and MacAyeal and extrude them
[9725]127 numnodes2d=md.mesh.numberofvertices2d;
128 numlayers=md.mesh.numberoflayers;
[8298]129 bordernodes2d=find(nodeonpattyn(1:numnodes2d) & nodeonmacayeal(1:numnodes2d)); %Nodes connected to two different types of elements
[5379]130
131 %initialize and fill in penalties structure
[8298]132 if ~isnan(bordernodes2d),
[5379]133 penalties=[];
134 for i=1:numlayers-1,
[9725]135 penalties=[penalties; [bordernodes2d bordernodes2d+md.mesh.numberofvertices2d*(i)]];
[5379]136 end
[9679]137 md.diagnostic.vertex_pairing=penalties;
[5379]138 end
139elseif strcmpi(coupling_method,'tiling'),
[5613]140 if any(macayealflag) & any(pattynflag), %coupling macayeal pattyn
[8298]141 %Find node at the border
142 nodeonmacayealpattyn(find(nodeonmacayeal & nodeonpattyn))=1;
[5613]143 %Macayeal elements in contact with this layer become MacAyealPattyn elements
[9733]144 matrixelements=ismember(md.mesh.elements,find(nodeonmacayealpattyn));
[5613]145 commonelements=sum(matrixelements,2)~=0;
146 commonelements(find(pattynflag))=0; %only one layer: the elements previously in macayeal
147 macayealflag(find(commonelements))=0; %these elements are now macayealpattynelements
148 macayealpattynflag(find(commonelements))=1;
[11004]149 nodeonmacayeal(:)=0;
[9733]150 nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1;
[5379]151
[11003]152 %rule out elements that don't touch the 2 boundaries
153 pos=find(macayealpattynflag);
[12959]154 elist=zeros(length(pos),1);
155 elist = elist + any(sum(nodeonmacayeal(md.mesh.elements(pos,:)),2),2);
156 elist = elist - any(sum(nodeonpattyn(md.mesh.elements(pos,:)) ,2),2);
157 pos1=find(elist==1);
[11003]158 macayealflag(pos(pos1))=1;
159 macayealpattynflag(pos(pos1))=0;
[12959]160 pos2=find(elist==-1);
[11003]161 pattynflag(pos(pos2))=1;
162 macayealpattynflag(pos(pos2))=0;
[5379]163
[11004]164 %Recompute nodes associated to these elements
[11003]165 nodeonmacayeal(:)=0;
166 nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1;
167 nodeonpattyn(:)=0;
168 nodeonpattyn(md.mesh.elements(find(pattynflag),:))=1;
169 nodeonmacayealpattyn(:)=0;
[9733]170 nodeonmacayealpattyn(md.mesh.elements(find(macayealpattynflag),:))=1;
[11003]171
[5613]172 elseif any(pattynflag) & any(stokesflag), %coupling pattyn stokes
[8298]173 %Find node at the border
174 nodeonpattynstokes(find(nodeonpattyn & nodeonstokes))=1;
[5613]175 %Stokes elements in contact with this layer become PattynStokes elements
[9733]176 matrixelements=ismember(md.mesh.elements,find(nodeonpattynstokes));
[5613]177 commonelements=sum(matrixelements,2)~=0;
178 commonelements(find(pattynflag))=0; %only one layer: the elements previously in macayeal
179 stokesflag(find(commonelements))=0; %these elements are now macayealpattynelements
180 pattynstokesflag(find(commonelements))=1;
[9725]181 nodeonstokes=zeros(md.mesh.numberofvertices,1);
[9733]182 nodeonstokes(md.mesh.elements(find(stokesflag),:))=1;
[5613]183
[11004]184 %rule out elements that don't touch the 2 boundaries
185 pos=find(pattynstokesflag);
[12959]186 elist=zeros(length(pos),1);
187 elist = elist + any(sum(nodeonstokes(md.mesh.elements(pos,:)),2),2);
188 elist = elist - any(sum(nodeonpattyn(md.mesh.elements(pos,:)),2),2);
189 pos1=find(elist==1);
[11004]190 stokesflag(pos(pos1))=1;
191 pattynstokesflag(pos(pos1))=0;
[12959]192 pos2=find(elist==-1);
[11004]193 pattynflag(pos(pos2))=1;
194 pattynstokesflag(pos(pos2))=0;
195
196 %Recompute nodes associated to these elements
197 nodeonstokes(:)=0;
198 nodeonstokes(md.mesh.elements(find(stokesflag),:))=1;
199 nodeonpattyn(:)=0;
200 nodeonpattyn(md.mesh.elements(find(pattynflag),:))=1;
201 nodeonpattynstokes(:)=0;
[9733]202 nodeonpattynstokes(md.mesh.elements(find(pattynstokesflag),:))=1;
[11004]203
[6705]204 elseif any(stokesflag) & any(macayealflag),
[8298]205 %Find node at the border
206 nodeonmacayealstokes(find(nodeonmacayeal & nodeonstokes))=1;
[6705]207 %Stokes elements in contact with this layer become MacAyealStokes elements
[9733]208 matrixelements=ismember(md.mesh.elements,find(nodeonmacayealstokes));
[6705]209 commonelements=sum(matrixelements,2)~=0;
210 commonelements(find(macayealflag))=0; %only one layer: the elements previously in macayeal
211 stokesflag(find(commonelements))=0; %these elements are now macayealmacayealelements
212 macayealstokesflag(find(commonelements))=1;
[9725]213 nodeonstokes=zeros(md.mesh.numberofvertices,1);
[9733]214 nodeonstokes(md.mesh.elements(find(stokesflag),:))=1;
[6705]215
[11004]216 %rule out elements that don't touch the 2 boundaries
217 pos=find(macayealstokesflag);
[12959]218 elist=zeros(length(pos),1);
219 elist = elist + any(sum(nodeonmacayeal(md.mesh.elements(pos,:)),2),2);
220 elist = elist - any(sum(nodeonstokes(md.mesh.elements(pos,:)) ,2),2);
221 pos1=find(elist==1);
[11004]222 macayealflag(pos(pos1))=1;
223 macayealstokesflag(pos(pos1))=0;
[12959]224 pos2=find(elist==-1);
[11004]225 stokesflag(pos(pos2))=1;
226 macayealstokesflag(pos(pos2))=0;
227
228 %Recompute nodes associated to these elements
229 nodeonmacayeal(:)=0;
230 nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1;
231 nodeonstokes(:)=0;
232 nodeonstokes(md.mesh.elements(find(stokesflag),:))=1;
233 nodeonmacayealstokes(:)=0;
[9733]234 nodeonmacayealstokes(md.mesh.elements(find(macayealstokesflag),:))=1;
[11004]235
[6705]236 elseif any(stokesflag) & any(hutterflag),
[5613]237 error('type of coupling not supported yet');
238 end
[5379]239end
240
[11003]241%Create MacaAyealPattynApproximation where needed
242md.flowequation.element_equation=zeros(md.mesh.numberofelements,1);
243md.flowequation.element_equation(find(noneflag))=0;
244md.flowequation.element_equation(find(hutterflag))=1;
245md.flowequation.element_equation(find(macayealflag))=2;
[13048]246md.flowequation.element_equation(find(l1l2flag))=8;
[11003]247md.flowequation.element_equation(find(pattynflag))=3;
248md.flowequation.element_equation(find(stokesflag))=4;
249md.flowequation.element_equation(find(macayealpattynflag))=5;
250md.flowequation.element_equation(find(macayealstokesflag))=6;
251md.flowequation.element_equation(find(pattynstokesflag))=7;
252
253%border
254md.flowequation.borderpattyn=nodeonpattyn;
255md.flowequation.bordermacayeal=nodeonmacayeal;
256md.flowequation.borderstokes=nodeonstokes;
257
[5071]258%Create vertices_type
[9725]259md.flowequation.vertex_equation=zeros(md.mesh.numberofvertices,1);
[8298]260pos=find(nodeonmacayeal);
[9661]261md.flowequation.vertex_equation(pos)=2;
[13048]262pos=find(nodeonl1l2);
263md.flowequation.vertex_equation(pos)=8;
[8298]264pos=find(nodeonpattyn);
[9661]265md.flowequation.vertex_equation(pos)=3;
[13058]266pos=find(nodeonhutter);
267md.flowequation.vertex_equation(pos)=1;
[8298]268pos=find(nodeonmacayealpattyn);
[9661]269md.flowequation.vertex_equation(pos)=5;
[8298]270pos=find(nodeonstokes);
[9661]271md.flowequation.vertex_equation(pos)=4;
[5433]272if any(stokesflag),
[8298]273 pos=find(~nodeonstokes);
[6712]274 if(~any(pattynflag) & ~any(macayealflag)),
[9661]275 md.flowequation.vertex_equation(pos)=0;
[5969]276 end
[5433]277end
[8298]278pos=find(nodeonpattynstokes);
[9661]279md.flowequation.vertex_equation(pos)=7;
[8298]280pos=find(nodeonmacayealstokes);
[9661]281md.flowequation.vertex_equation(pos)=6;
[5071]282
[202]283%figure out solution types
[9661]284md.flowequation.ishutter=double(any(md.flowequation.element_equation==1));
285md.flowequation.ismacayealpattyn=double(any(md.flowequation.element_equation==2 | md.flowequation.element_equation==3));
286md.flowequation.isstokes=double(any(md.flowequation.element_equation==4));
[13048]287md.flowequation.isl1l2=double(any(md.flowequation.element_equation==8));
[202]288
[11003]289return
290
291%Check that tiling can work:
292if any(md.flowequation.bordermacayeal) & any(md.flowequation.borderpattyn) & any(md.flowequation.borderpattyn + md.flowequation.bordermacayeal ~=1),
293 error('error coupling domain too irregular');
[1]294end
[11003]295if any(md.flowequation.bordermacayeal) & any(md.flowequation.borderstokes) & any(md.flowequation.borderstokes + md.flowequation.bordermacayeal ~=1),
296 error('error coupling domain too irregular');
297end
298if any(md.flowequation.borderstokes) & any(md.flowequation.borderpattyn) & any(md.flowequation.borderpattyn + md.flowequation.borderstokes~=1),
299 error('error coupling domain too irregular');
300end
Note: See TracBrowser for help on using the repository browser.