function md=setflowequation(md,varargin) %SETELEMENTSTYPE - associate a solution type to each element % % This routine works like plotmodel: it works with an even number of inputs % 'hutter','macayeal','l1l2','pattyn','stokes' and 'fill' are the possible options % that must be followed by the corresponding exp file or flags list % It can either be a domain file (argus type, .exp extension), or an array of element flags. % If user wants every element outside the domain to be % setflowequationd, add '~' to the name of the domain file (ex: '~Pattyn.exp'); % an empty string '' will be considered as an empty domain % a string 'all' will be considered as the entire domain % You can specify the type of coupling, 'penalties' or 'tiling', to use with the input 'coupling' % NB: l1l2 cannot currently be coupled to any other ice flow model % % Usage: % md=setflowequation(md,varargin) % % Example: % md=setflowequation(md,'pattyn','Pattyn.exp','macayeal',md.mask.elementonfloatingice,'fill','hutter'); % md=setflowequation(md,'pattyn','Pattyn.exp',fill','hutter','coupling','tiling'); %some checks on list of arguments if ((nargin<2) | (nargout~=1)), error('setflowequation error message'); end %Process options options=pairoptions(varargin{:}); options=deleteduplicates(options,1); %Find_out what kind of coupling to use coupling_method=getfieldvalue(options,'coupling','tiling'); if (~strcmpi(coupling_method,'tiling') & ~strcmpi(coupling_method,'penalties')), error('coupling type can only be: tiling or penalties'); end %recover elements distribution hutterflag = FlagElements(md,getfieldvalue(options,'hutter','')); macayealflag = FlagElements(md,getfieldvalue(options,'macayeal','')); pattynflag = FlagElements(md,getfieldvalue(options,'pattyn','')); l1l2flag = FlagElements(md,getfieldvalue(options,'l1l2','')); stokesflag = FlagElements(md,getfieldvalue(options,'stokes','')); filltype = getfieldvalue(options,'fill','none'); %Flag the elements that have not been flagged as filltype if strcmpi(filltype,'hutter'), hutterflag(find(~(macayealflag | pattynflag)))=1; elseif strcmpi(filltype,'macayeal'), macayealflag(find(~(hutterflag | pattynflag | stokesflag)))=1; elseif strcmpi(filltype,'pattyn'), pattynflag(find(~(hutterflag | macayealflag | stokesflag)))=1; end %check that each element has at least one flag if any(hutterflag+macayealflag+pattynflag+l1l2flag+stokesflag==0), error('elements type not assigned, must be specified') end %check that each element has only one flag if any(hutterflag+macayealflag+pattynflag+l1l2flag+stokesflag>1), disp('setflowequation warning message: some elements have several types, higher order type is used for them') hutterflag(find(hutterflag & macayealflag))=0; hutterflag(find(hutterflag & pattynflag))=0; macayealflag(find(macayealflag & pattynflag))=0; end %check that l1l2 is not coupled to any other model for now if any(l1l2flag) & any(hutterflag | macayealflag | pattynflag | stokesflag) error('l1l2 cannot be coupled to any other model'); end %Check that no l1l2 or pattyn or stokes for 2d mesh if (md.mesh.dimension==2), if any(l1l2flag | stokesflag | pattynflag) error('stokes and pattyn elements not allowed in 2d mesh, extrude it first') end end %Stokes can only be used alone for now: if any(stokesflag) &any(hutterflag), error('stokes cannot be used with any other model for now, put stokes everywhere') end %Initialize node fields nodeonhutter=zeros(md.mesh.numberofvertices,1); nodeonhutter(md.mesh.elements(find(hutterflag),:))=1; nodeonmacayeal=zeros(md.mesh.numberofvertices,1); nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1; nodeonpattyn=zeros(md.mesh.numberofvertices,1); nodeonpattyn(md.mesh.elements(find(pattynflag),:))=1; nodeonl1l2=zeros(md.mesh.numberofvertices,1); nodeonl1l2(md.mesh.elements(find(l1l2flag),:))=1; nodeonstokes=zeros(md.mesh.numberofvertices,1); noneflag=zeros(md.mesh.numberofelements,1); %First modify stokesflag to get rid of elements contrained everywhere (spc + border with pattyn or macayeal) if any(stokesflag), 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 fullspcelems=double(sum(fullspcnodes(md.mesh.elements),2)==6); %find all the nodes on the boundary of the domain without icefront stokesflag(find(fullspcelems))=0; nodeonstokes(md.mesh.elements(find(stokesflag),:))=1; end %Then complete with NoneApproximation or the other model used if there is no stokes if any(stokesflag), if any(pattynflag), %fill with pattyn pattynflag(~stokesflag)=1; nodeonpattyn(md.mesh.elements(find(pattynflag),:))=1; elseif any(macayealflag), %fill with macayeal macayealflag(~stokesflag)=1; nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1; else %fill with none noneflag(find(~stokesflag))=1; end end %Now take care of the coupling between MacAyeal and Pattyn md.diagnostic.vertex_pairing=[]; nodeonmacayealpattyn=zeros(md.mesh.numberofvertices,1); nodeonpattynstokes=zeros(md.mesh.numberofvertices,1); nodeonmacayealstokes=zeros(md.mesh.numberofvertices,1); macayealpattynflag=zeros(md.mesh.numberofelements,1); macayealstokesflag=zeros(md.mesh.numberofelements,1); pattynstokesflag=zeros(md.mesh.numberofelements,1); if strcmpi(coupling_method,'penalties'), %Create the border nodes between Pattyn and MacAyeal and extrude them numnodes2d=md.mesh.numberofvertices2d; numlayers=md.mesh.numberoflayers; bordernodes2d=find(nodeonpattyn(1:numnodes2d) & nodeonmacayeal(1:numnodes2d)); %Nodes connected to two different types of elements %initialize and fill in penalties structure if ~isnan(bordernodes2d), penalties=[]; for i=1:numlayers-1, penalties=[penalties; [bordernodes2d bordernodes2d+md.mesh.numberofvertices2d*(i)]]; end md.diagnostic.vertex_pairing=penalties; end elseif strcmpi(coupling_method,'tiling'), if any(macayealflag) & any(pattynflag), %coupling macayeal pattyn %Find node at the border nodeonmacayealpattyn(find(nodeonmacayeal & nodeonpattyn))=1; %Macayeal elements in contact with this layer become MacAyealPattyn elements matrixelements=ismember(md.mesh.elements,find(nodeonmacayealpattyn)); commonelements=sum(matrixelements,2)~=0; commonelements(find(pattynflag))=0; %only one layer: the elements previously in macayeal macayealflag(find(commonelements))=0; %these elements are now macayealpattynelements macayealpattynflag(find(commonelements))=1; nodeonmacayeal(:)=0; nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1; %rule out elements that don't touch the 2 boundaries pos=find(macayealpattynflag); elist=zeros(length(pos),1); elist = elist + any(sum(nodeonmacayeal(md.mesh.elements(pos,:)),2),2); elist = elist - any(sum(nodeonpattyn(md.mesh.elements(pos,:)) ,2),2); pos1=find(elist==1); macayealflag(pos(pos1))=1; macayealpattynflag(pos(pos1))=0; pos2=find(elist==-1); pattynflag(pos(pos2))=1; macayealpattynflag(pos(pos2))=0; %Recompute nodes associated to these elements nodeonmacayeal(:)=0; nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1; nodeonpattyn(:)=0; nodeonpattyn(md.mesh.elements(find(pattynflag),:))=1; nodeonmacayealpattyn(:)=0; nodeonmacayealpattyn(md.mesh.elements(find(macayealpattynflag),:))=1; elseif any(pattynflag) & any(stokesflag), %coupling pattyn stokes %Find node at the border nodeonpattynstokes(find(nodeonpattyn & nodeonstokes))=1; %Stokes elements in contact with this layer become PattynStokes elements matrixelements=ismember(md.mesh.elements,find(nodeonpattynstokes)); commonelements=sum(matrixelements,2)~=0; commonelements(find(pattynflag))=0; %only one layer: the elements previously in macayeal stokesflag(find(commonelements))=0; %these elements are now macayealpattynelements pattynstokesflag(find(commonelements))=1; nodeonstokes=zeros(md.mesh.numberofvertices,1); nodeonstokes(md.mesh.elements(find(stokesflag),:))=1; %rule out elements that don't touch the 2 boundaries pos=find(pattynstokesflag); elist=zeros(length(pos),1); elist = elist + any(sum(nodeonstokes(md.mesh.elements(pos,:)),2),2); elist = elist - any(sum(nodeonpattyn(md.mesh.elements(pos,:)),2),2); pos1=find(elist==1); stokesflag(pos(pos1))=1; pattynstokesflag(pos(pos1))=0; pos2=find(elist==-1); pattynflag(pos(pos2))=1; pattynstokesflag(pos(pos2))=0; %Recompute nodes associated to these elements nodeonstokes(:)=0; nodeonstokes(md.mesh.elements(find(stokesflag),:))=1; nodeonpattyn(:)=0; nodeonpattyn(md.mesh.elements(find(pattynflag),:))=1; nodeonpattynstokes(:)=0; nodeonpattynstokes(md.mesh.elements(find(pattynstokesflag),:))=1; elseif any(stokesflag) & any(macayealflag), %Find node at the border nodeonmacayealstokes(find(nodeonmacayeal & nodeonstokes))=1; %Stokes elements in contact with this layer become MacAyealStokes elements matrixelements=ismember(md.mesh.elements,find(nodeonmacayealstokes)); commonelements=sum(matrixelements,2)~=0; commonelements(find(macayealflag))=0; %only one layer: the elements previously in macayeal stokesflag(find(commonelements))=0; %these elements are now macayealmacayealelements macayealstokesflag(find(commonelements))=1; nodeonstokes=zeros(md.mesh.numberofvertices,1); nodeonstokes(md.mesh.elements(find(stokesflag),:))=1; %rule out elements that don't touch the 2 boundaries pos=find(macayealstokesflag); elist=zeros(length(pos),1); elist = elist + any(sum(nodeonmacayeal(md.mesh.elements(pos,:)),2),2); elist = elist - any(sum(nodeonstokes(md.mesh.elements(pos,:)) ,2),2); pos1=find(elist==1); macayealflag(pos(pos1))=1; macayealstokesflag(pos(pos1))=0; pos2=find(elist==-1); stokesflag(pos(pos2))=1; macayealstokesflag(pos(pos2))=0; %Recompute nodes associated to these elements nodeonmacayeal(:)=0; nodeonmacayeal(md.mesh.elements(find(macayealflag),:))=1; nodeonstokes(:)=0; nodeonstokes(md.mesh.elements(find(stokesflag),:))=1; nodeonmacayealstokes(:)=0; nodeonmacayealstokes(md.mesh.elements(find(macayealstokesflag),:))=1; elseif any(stokesflag) & any(hutterflag), error('type of coupling not supported yet'); end end %Create MacaAyealPattynApproximation where needed md.flowequation.element_equation=zeros(md.mesh.numberofelements,1); md.flowequation.element_equation(find(noneflag))=0; md.flowequation.element_equation(find(hutterflag))=1; md.flowequation.element_equation(find(macayealflag))=2; md.flowequation.element_equation(find(l1l2flag))=8; md.flowequation.element_equation(find(pattynflag))=3; md.flowequation.element_equation(find(stokesflag))=4; md.flowequation.element_equation(find(macayealpattynflag))=5; md.flowequation.element_equation(find(macayealstokesflag))=6; md.flowequation.element_equation(find(pattynstokesflag))=7; %border md.flowequation.borderpattyn=nodeonpattyn; md.flowequation.bordermacayeal=nodeonmacayeal; md.flowequation.borderstokes=nodeonstokes; %Create vertices_type md.flowequation.vertex_equation=zeros(md.mesh.numberofvertices,1); pos=find(nodeonmacayeal); md.flowequation.vertex_equation(pos)=2; pos=find(nodeonl1l2); md.flowequation.vertex_equation(pos)=8; pos=find(nodeonpattyn); md.flowequation.vertex_equation(pos)=3; pos=find(nodeonhutter); md.flowequation.vertex_equation(pos)=1; pos=find(nodeonmacayealpattyn); md.flowequation.vertex_equation(pos)=5; pos=find(nodeonstokes); md.flowequation.vertex_equation(pos)=4; if any(stokesflag), pos=find(~nodeonstokes); if(~any(pattynflag) & ~any(macayealflag)), md.flowequation.vertex_equation(pos)=0; end end pos=find(nodeonpattynstokes); md.flowequation.vertex_equation(pos)=7; pos=find(nodeonmacayealstokes); md.flowequation.vertex_equation(pos)=6; %figure out solution types md.flowequation.ishutter=double(any(md.flowequation.element_equation==1)); md.flowequation.ismacayealpattyn=double(any(md.flowequation.element_equation==2 | md.flowequation.element_equation==3)); md.flowequation.isstokes=double(any(md.flowequation.element_equation==4)); md.flowequation.isl1l2=double(any(md.flowequation.element_equation==8)); return %Check that tiling can work: if any(md.flowequation.bordermacayeal) & any(md.flowequation.borderpattyn) & any(md.flowequation.borderpattyn + md.flowequation.bordermacayeal ~=1), error('error coupling domain too irregular'); end if any(md.flowequation.bordermacayeal) & any(md.flowequation.borderstokes) & any(md.flowequation.borderstokes + md.flowequation.bordermacayeal ~=1), error('error coupling domain too irregular'); end if any(md.flowequation.borderstokes) & any(md.flowequation.borderpattyn) & any(md.flowequation.borderpattyn + md.flowequation.borderstokes~=1), error('error coupling domain too irregular'); end