[9664] | 1 | function 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
|
---|
| 23 | if ((nargin<2) | (nargout~=1)),
|
---|
[9664] | 24 | error('setflowequation error message');
|
---|
[1] | 25 | end
|
---|
| 26 |
|
---|
[13028] | 27 | %Process options
|
---|
| 28 | options=pairoptions(varargin{:});
|
---|
| 29 | options=deleteduplicates(options,1);
|
---|
| 30 |
|
---|
[5379] | 31 | %Find_out what kind of coupling to use
|
---|
| 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');
|
---|
| 35 | end
|
---|
| 36 |
|
---|
[13028] | 37 | %recover elements distribution
|
---|
| 38 | hutterflag = FlagElements(md,getfieldvalue(options,'hutter',''));
|
---|
| 39 | macayealflag = FlagElements(md,getfieldvalue(options,'macayeal',''));
|
---|
| 40 | pattynflag = FlagElements(md,getfieldvalue(options,'pattyn',''));
|
---|
| 41 | l1l2flag = FlagElements(md,getfieldvalue(options,'l1l2',''));
|
---|
| 42 | stokesflag = FlagElements(md,getfieldvalue(options,'stokes',''));
|
---|
| 43 | filltype = getfieldvalue(options,'fill','none');
|
---|
[1] | 44 |
|
---|
[12888] | 45 | %Flag the elements that have not been flagged as filltype
|
---|
[1] | 46 | if strcmpi(filltype,'hutter'),
|
---|
[12949] | 47 | hutterflag(find(~(macayealflag | pattynflag)))=1;
|
---|
[1] | 48 | elseif strcmpi(filltype,'macayeal'),
|
---|
[12950] | 49 | macayealflag(find(~(hutterflag | pattynflag | stokesflag)))=1;
|
---|
[1] | 50 | elseif strcmpi(filltype,'pattyn'),
|
---|
[12950] | 51 | pattynflag(find(~(hutterflag | macayealflag | stokesflag)))=1;
|
---|
[1] | 52 | end
|
---|
| 53 |
|
---|
| 54 | %check that each element has at least one flag
|
---|
[13048] | 55 | if any(hutterflag+macayealflag+pattynflag+l1l2flag+stokesflag==0),
|
---|
| 56 | error('elements type not assigned, must be specified')
|
---|
[1] | 57 | end
|
---|
| 58 |
|
---|
| 59 | %check that each element has only one flag
|
---|
[13048] | 60 | if 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;
|
---|
| 65 | end
|
---|
| 66 |
|
---|
[13048] | 67 | %check that l1l2 is not coupled to any other model for now
|
---|
| 68 | if any(l1l2flag) & any(hutterflag | macayealflag | pattynflag | stokesflag)
|
---|
| 69 | error('l1l2 cannot be coupled to any other model');
|
---|
| 70 | end
|
---|
| 71 |
|
---|
| 72 | %Check that no l1l2 or pattyn or stokes for 2d mesh
|
---|
[9719] | 73 | if (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
|
---|
| 77 | end
|
---|
| 78 |
|
---|
[5433] | 79 | %Stokes can only be used alone for now:
|
---|
[6705] | 80 | if any(stokesflag) &any(hutterflag),
|
---|
[13048] | 81 | error('stokes cannot be used with any other model for now, put stokes everywhere')
|
---|
[5433] | 82 | end
|
---|
[1] | 83 |
|
---|
[11004] | 84 | %Initialize node fields
|
---|
[9725] | 85 | nodeonhutter=zeros(md.mesh.numberofvertices,1);
|
---|
[9733] | 86 | nodeonhutter(md.mesh.elements(find(hutterflag),:))=1;
|
---|
[9725] | 87 | nodeonmacayeal=zeros(md.mesh.numberofvertices,1);
|
---|
[9733] | 88 | nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1;
|
---|
[9725] | 89 | nodeonpattyn=zeros(md.mesh.numberofvertices,1);
|
---|
[9733] | 90 | nodeonpattyn(md.mesh.elements(find(pattynflag),:))=1;
|
---|
[13048] | 91 | nodeonl1l2=zeros(md.mesh.numberofvertices,1);
|
---|
| 92 | nodeonl1l2(md.mesh.elements(find(l1l2flag),:))=1;
|
---|
[11004] | 93 | nodeonstokes=zeros(md.mesh.numberofvertices,1);
|
---|
| 94 | noneflag=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)
|
---|
| 97 | if 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] | 102 | end
|
---|
[1] | 103 |
|
---|
[5613] | 104 | %Then complete with NoneApproximation or the other model used if there is no stokes
|
---|
[5433] | 105 | if 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] | 115 | end
|
---|
[1] | 116 |
|
---|
[5379] | 117 | %Now take care of the coupling between MacAyeal and Pattyn
|
---|
[9679] | 118 | md.diagnostic.vertex_pairing=[];
|
---|
[9725] | 119 | nodeonmacayealpattyn=zeros(md.mesh.numberofvertices,1);
|
---|
| 120 | nodeonpattynstokes=zeros(md.mesh.numberofvertices,1);
|
---|
| 121 | nodeonmacayealstokes=zeros(md.mesh.numberofvertices,1);
|
---|
[11003] | 122 | macayealpattynflag=zeros(md.mesh.numberofelements,1);
|
---|
| 123 | macayealstokesflag=zeros(md.mesh.numberofelements,1);
|
---|
| 124 | pattynstokesflag=zeros(md.mesh.numberofelements,1);
|
---|
[5379] | 125 | if 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
|
---|
| 139 | elseif 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] | 239 | end
|
---|
| 240 |
|
---|
[11003] | 241 | %Create MacaAyealPattynApproximation where needed
|
---|
| 242 | md.flowequation.element_equation=zeros(md.mesh.numberofelements,1);
|
---|
| 243 | md.flowequation.element_equation(find(noneflag))=0;
|
---|
| 244 | md.flowequation.element_equation(find(hutterflag))=1;
|
---|
| 245 | md.flowequation.element_equation(find(macayealflag))=2;
|
---|
[13048] | 246 | md.flowequation.element_equation(find(l1l2flag))=8;
|
---|
[11003] | 247 | md.flowequation.element_equation(find(pattynflag))=3;
|
---|
| 248 | md.flowequation.element_equation(find(stokesflag))=4;
|
---|
| 249 | md.flowequation.element_equation(find(macayealpattynflag))=5;
|
---|
| 250 | md.flowequation.element_equation(find(macayealstokesflag))=6;
|
---|
| 251 | md.flowequation.element_equation(find(pattynstokesflag))=7;
|
---|
| 252 |
|
---|
| 253 | %border
|
---|
| 254 | md.flowequation.borderpattyn=nodeonpattyn;
|
---|
| 255 | md.flowequation.bordermacayeal=nodeonmacayeal;
|
---|
| 256 | md.flowequation.borderstokes=nodeonstokes;
|
---|
| 257 |
|
---|
[5071] | 258 | %Create vertices_type
|
---|
[9725] | 259 | md.flowequation.vertex_equation=zeros(md.mesh.numberofvertices,1);
|
---|
[8298] | 260 | pos=find(nodeonmacayeal);
|
---|
[9661] | 261 | md.flowequation.vertex_equation(pos)=2;
|
---|
[13048] | 262 | pos=find(nodeonl1l2);
|
---|
| 263 | md.flowequation.vertex_equation(pos)=8;
|
---|
[8298] | 264 | pos=find(nodeonpattyn);
|
---|
[9661] | 265 | md.flowequation.vertex_equation(pos)=3;
|
---|
[13058] | 266 | pos=find(nodeonhutter);
|
---|
| 267 | md.flowequation.vertex_equation(pos)=1;
|
---|
[8298] | 268 | pos=find(nodeonmacayealpattyn);
|
---|
[9661] | 269 | md.flowequation.vertex_equation(pos)=5;
|
---|
[8298] | 270 | pos=find(nodeonstokes);
|
---|
[9661] | 271 | md.flowequation.vertex_equation(pos)=4;
|
---|
[5433] | 272 | if 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] | 277 | end
|
---|
[8298] | 278 | pos=find(nodeonpattynstokes);
|
---|
[9661] | 279 | md.flowequation.vertex_equation(pos)=7;
|
---|
[8298] | 280 | pos=find(nodeonmacayealstokes);
|
---|
[9661] | 281 | md.flowequation.vertex_equation(pos)=6;
|
---|
[5071] | 282 |
|
---|
[202] | 283 | %figure out solution types
|
---|
[9661] | 284 | md.flowequation.ishutter=double(any(md.flowequation.element_equation==1));
|
---|
| 285 | md.flowequation.ismacayealpattyn=double(any(md.flowequation.element_equation==2 | md.flowequation.element_equation==3));
|
---|
| 286 | md.flowequation.isstokes=double(any(md.flowequation.element_equation==4));
|
---|
[13048] | 287 | md.flowequation.isl1l2=double(any(md.flowequation.element_equation==8));
|
---|
[202] | 288 |
|
---|
[11003] | 289 | return
|
---|
| 290 |
|
---|
| 291 | %Check that tiling can work:
|
---|
| 292 | if any(md.flowequation.bordermacayeal) & any(md.flowequation.borderpattyn) & any(md.flowequation.borderpattyn + md.flowequation.bordermacayeal ~=1),
|
---|
| 293 | error('error coupling domain too irregular');
|
---|
[1] | 294 | end
|
---|
[11003] | 295 | if any(md.flowequation.bordermacayeal) & any(md.flowequation.borderstokes) & any(md.flowequation.borderstokes + md.flowequation.bordermacayeal ~=1),
|
---|
| 296 | error('error coupling domain too irregular');
|
---|
| 297 | end
|
---|
| 298 | if any(md.flowequation.borderstokes) & any(md.flowequation.borderpattyn) & any(md.flowequation.borderpattyn + md.flowequation.borderstokes~=1),
|
---|
| 299 | error('error coupling domain too irregular');
|
---|
| 300 | end
|
---|