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

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