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

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
%   'SIA','SSA','L1L2','HO','FS' 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: '~HO.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,'HO','HO.exp',fill','SIA','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
SIAflag  = FlagElements(md,getfieldvalue(options,'SIA',''));
SSAflag  = FlagElements(md,getfieldvalue(options,'SSA',''));
HOflag   = FlagElements(md,getfieldvalue(options,'HO',''));
L1L2flag = FlagElements(md,getfieldvalue(options,'L1L2',''));
FSflag   = FlagElements(md,getfieldvalue(options,'FS',''));
filltype = getfieldvalue(options,'fill','none');
displayunused(options);

%Flag the elements that have not been flagged as filltype
if strcmpi(filltype,'SIA'),
        SIAflag(find(~(SSAflag | HOflag)))=1;
elseif strcmpi(filltype,'SSA'),
        SSAflag(find(~(SIAflag | HOflag | FSflag)))=1;
elseif strcmpi(filltype,'HO'),
        HOflag(find(~(SIAflag | SSAflag | FSflag)))=1;
end

%check that each element has at least one flag
if any(SIAflag+SSAflag+HOflag+L1L2flag+FSflag==0),
        error('elements type not assigned, supported models are ''SIA'',''SSA'',''HO'' and ''FS''')
end

%check that each element has only one flag
if any(SIAflag+SSAflag+HOflag+L1L2flag+FSflag>1),
        disp('setflowequation warning message: some elements have several types, higher order type is used for them')
        SIAflag(find(SIAflag & SSAflag))=0;
        SIAflag(find(SIAflag & HOflag))=0;
        SSAflag(find(SSAflag & HOflag))=0;
end

%check that L1L2 is not coupled to any other model for now
if any(L1L2flag) & any(SIAflag | SSAflag | HOflag | FSflag)
        error('L1L2 cannot be coupled to any other model');
end

%Check that no L1L2 or HO or FS for 2d mesh
if strcmp(meshtype(md.mesh),'2Dhorizontal')
        if any(L1L2flag | FSflag | HOflag)
                error('FS and HO elements not allowed in 2d mesh, extrude it first')
        end
end

%FS can only be used alone for now:
if any(FSflag) &any(SIAflag),
        error('FS cannot be used with any other model for now, put FS everywhere')
end

%Initialize node fields
nodeonSIA=zeros(md.mesh.numberofvertices,1);
nodeonSIA(md.mesh.elements(find(SIAflag),:))=1;
nodeonSSA=zeros(md.mesh.numberofvertices,1);
nodeonSSA(md.mesh.elements(find(SSAflag),:))=1;
nodeonHO=zeros(md.mesh.numberofvertices,1);
nodeonHO(md.mesh.elements(find(HOflag),:))=1;
nodeonL1L2=zeros(md.mesh.numberofvertices,1);
nodeonL1L2(md.mesh.elements(find(L1L2flag),:))=1;
nodeonFS=zeros(md.mesh.numberofvertices,1);
noneflag=zeros(md.mesh.numberofelements,1);

%First modify FSflag to get rid of elements contrained everywhere (spc + border with HO or SSA)
if any(FSflag),
        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
        fullspcelems=double(sum(fullspcnodes(md.mesh.elements),2)==6);         %find all the nodes on the boundary of the domain without icefront
        FSflag(find(fullspcelems))=0;
        nodeonFS(md.mesh.elements(find(FSflag),:))=1;
end

%Then complete with NoneApproximation or the other model used if there is no FS
if any(FSflag), 
        if any(HOflag), %fill with HO
                HOflag(~FSflag)=1;
                nodeonHO(md.mesh.elements(find(HOflag),:))=1;
        elseif any(SSAflag), %fill with SSA
                SSAflag(~FSflag)=1;
                nodeonSSA(md.mesh.elements(find(SSAflag),:))=1;
        else %fill with none 
                noneflag(find(~FSflag))=1;
        end
end

%Now take care of the coupling between SSA and HO
md.stressbalance.vertex_pairing=[];
nodeonSSAHO=zeros(md.mesh.numberofvertices,1);
nodeonHOFS=zeros(md.mesh.numberofvertices,1);
nodeonSSAFS=zeros(md.mesh.numberofvertices,1);
SSAHOflag=zeros(md.mesh.numberofelements,1);
SSAFSflag=zeros(md.mesh.numberofelements,1);
HOFSflag=zeros(md.mesh.numberofelements,1);
if strcmpi(coupling_method,'penalties'),
        %Create the border nodes between HO and SSA and extrude them
        numnodes2d=md.mesh.numberofvertices2d;
        numlayers=md.mesh.numberoflayers;
        bordernodes2d=find(nodeonHO(1:numnodes2d) & nodeonSSA(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.stressbalance.vertex_pairing=penalties;
        end
elseif strcmpi(coupling_method,'tiling'),
        if any(SSAflag) & any(HOflag), %coupling SSA HO
                %Find node at the border
                nodeonSSAHO(find(nodeonSSA & nodeonHO))=1;
                %SSA elements in contact with this layer become SSAHO elements
                matrixelements=ismember(md.mesh.elements,find(nodeonSSAHO));
                commonelements=sum(matrixelements,2)~=0;
                commonelements(find(HOflag))=0; %only one layer: the elements previously in SSA
                SSAflag(find(commonelements))=0; %these elements are now SSAHOelements
                SSAHOflag(find(commonelements))=1;
                nodeonSSA(:)=0;
                nodeonSSA(md.mesh.elements(find(SSAflag),:))=1;

                %rule out elements that don't touch the 2 boundaries
                pos=find(SSAHOflag);
                elist=zeros(length(pos),1);
                elist = elist + any(sum(nodeonSSA(md.mesh.elements(pos,:)),2),2);
                elist = elist - any(sum(nodeonHO(md.mesh.elements(pos,:))  ,2),2);
                pos1=find(elist==1);
                SSAflag(pos(pos1))=1;
                SSAHOflag(pos(pos1))=0;
                pos2=find(elist==-1);
                HOflag(pos(pos2))=1;
                SSAHOflag(pos(pos2))=0;

                %Recompute nodes associated to these elements
                nodeonSSA(:)=0;
                nodeonSSA(md.mesh.elements(find(SSAflag),:))=1;
                nodeonHO(:)=0;
                nodeonHO(md.mesh.elements(find(HOflag),:))=1;
                nodeonSSAHO(:)=0;
                nodeonSSAHO(md.mesh.elements(find(SSAHOflag),:))=1;

        elseif any(HOflag) & any(FSflag), %coupling HO FS
                %Find node at the border
                nodeonHOFS(find(nodeonHO & nodeonFS))=1;
                %FS elements in contact with this layer become HOFS elements
                matrixelements=ismember(md.mesh.elements,find(nodeonHOFS));
                commonelements=sum(matrixelements,2)~=0;
                commonelements(find(HOflag))=0; %only one layer: the elements previously in SSA
                FSflag(find(commonelements))=0; %these elements are now SSAHOelements
                HOFSflag(find(commonelements))=1;
                nodeonFS=zeros(md.mesh.numberofvertices,1);
                nodeonFS(md.mesh.elements(find(FSflag),:))=1;

                %rule out elements that don't touch the 2 boundaries
                pos=find(HOFSflag);
                elist=zeros(length(pos),1);
                elist = elist + any(sum(nodeonFS(md.mesh.elements(pos,:)),2),2);
                elist = elist - any(sum(nodeonHO(md.mesh.elements(pos,:)),2),2);
                pos1=find(elist==1);
                FSflag(pos(pos1))=1;
                HOFSflag(pos(pos1))=0;
                pos2=find(elist==-1);
                HOflag(pos(pos2))=1;
                HOFSflag(pos(pos2))=0;

                %Recompute nodes associated to these elements
                nodeonFS(:)=0;
                nodeonFS(md.mesh.elements(find(FSflag),:))=1;
                nodeonHO(:)=0;
                nodeonHO(md.mesh.elements(find(HOflag),:))=1;
                nodeonHOFS(:)=0;
                nodeonHOFS(md.mesh.elements(find(HOFSflag),:))=1;

        elseif any(FSflag) & any(SSAflag),
                %Find node at the border
                nodeonSSAFS(find(nodeonSSA & nodeonFS))=1;
                %FS elements in contact with this layer become SSAFS elements
                matrixelements=ismember(md.mesh.elements,find(nodeonSSAFS));
                commonelements=sum(matrixelements,2)~=0;
                commonelements(find(SSAflag))=0; %only one layer: the elements previously in SSA
                FSflag(find(commonelements))=0; %these elements are now SSASSAelements
                SSAFSflag(find(commonelements))=1;
                nodeonFS=zeros(md.mesh.numberofvertices,1);
                nodeonFS(md.mesh.elements(find(FSflag),:))=1;

                %rule out elements that don't touch the 2 boundaries
                pos=find(SSAFSflag);
                elist=zeros(length(pos),1);
                elist = elist + any(sum(nodeonSSA(md.mesh.elements(pos,:)),2),2);
                elist = elist - any(sum(nodeonFS(md.mesh.elements(pos,:))  ,2),2);
                pos1=find(elist==1);
                SSAflag(pos(pos1))=1;
                SSAFSflag(pos(pos1))=0;
                pos2=find(elist==-1);
                FSflag(pos(pos2))=1;
                SSAFSflag(pos(pos2))=0;

                %Recompute nodes associated to these elements
                nodeonSSA(:)=0;
                nodeonSSA(md.mesh.elements(find(SSAflag),:))=1;
                nodeonFS(:)=0;
                nodeonFS(md.mesh.elements(find(FSflag),:))=1;
                nodeonSSAFS(:)=0;
                nodeonSSAFS(md.mesh.elements(find(SSAFSflag),:))=1;

        elseif any(FSflag) & any(SIAflag),
                error('type of coupling not supported yet');
        end
end

%Create MacaAyealHOApproximation where needed
md.flowequation.element_equation=zeros(md.mesh.numberofelements,1);
md.flowequation.element_equation(find(noneflag))=0;
md.flowequation.element_equation(find(SIAflag))=1;
md.flowequation.element_equation(find(SSAflag))=2;
md.flowequation.element_equation(find(L1L2flag))=8;
md.flowequation.element_equation(find(HOflag))=3;
md.flowequation.element_equation(find(FSflag))=4;
md.flowequation.element_equation(find(SSAHOflag))=5;
md.flowequation.element_equation(find(SSAFSflag))=6;
md.flowequation.element_equation(find(HOFSflag))=7;

%border
md.flowequation.borderHO=nodeonHO;
md.flowequation.borderSSA=nodeonSSA;
md.flowequation.borderFS=nodeonFS;

%Create vertices_type
md.flowequation.vertex_equation=zeros(md.mesh.numberofvertices,1);
pos=find(nodeonSSA);
md.flowequation.vertex_equation(pos)=2;
pos=find(nodeonL1L2);
md.flowequation.vertex_equation(pos)=8;
pos=find(nodeonHO);
md.flowequation.vertex_equation(pos)=3;
pos=find(nodeonSIA);
md.flowequation.vertex_equation(pos)=1;
pos=find(nodeonSSAHO);
md.flowequation.vertex_equation(pos)=5;
pos=find(nodeonFS);
md.flowequation.vertex_equation(pos)=4;
if any(FSflag),
        pos=find(~nodeonFS);
        if(~any(HOflag) & ~any(SSAflag)),
                md.flowequation.vertex_equation(pos)=0;
        end
end
pos=find(nodeonHOFS);
md.flowequation.vertex_equation(pos)=7;
pos=find(nodeonSSAFS);
md.flowequation.vertex_equation(pos)=6;

%figure out solution types
md.flowequation.isSIA  = double(any(md.flowequation.element_equation == 1));
md.flowequation.isSSA  = double(any(md.flowequation.element_equation == 2));
md.flowequation.isHO   = double(any(md.flowequation.element_equation == 3));
md.flowequation.isL1L2 = double(any(md.flowequation.element_equation == 8));
md.flowequation.isFS   = double(any(md.flowequation.element_equation == 4));

return

%Check that tiling can work:
if any(md.flowequation.borderSSA) & any(md.flowequation.borderHO) & any(md.flowequation.borderHO + md.flowequation.borderSSA ~=1),
        error('error coupling domain too irregular');
end
if any(md.flowequation.borderSSA) & any(md.flowequation.borderFS) & any(md.flowequation.borderFS + md.flowequation.borderSSA ~=1),
        error('error coupling domain too irregular');
end
if any(md.flowequation.borderFS) & any(md.flowequation.borderHO) & any(md.flowequation.borderHO + md.flowequation.borderFS~=1),
        error('error coupling domain too irregular');
end