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