function md=setelementstype(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','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 
%   setelementstyped, 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
%
%   Usage:
%      md=setelementstype(md,varargin)
%
%   Example:
%      md=setelementstype(md,'pattyn','Pattyn.exp','macayeal',md.elementoniceshelf,'fill','hutter');

%some checks on list of arguments
if ((nargin<2) | (nargout~=1)),
	error('setelementstype error message');
end

if md.counter<2,
	error('only fully parameterized 2d models can be setelementstyped');
end

[hutterflag macayealflag pattynflag stokesflag filltype]=recover_areas(md,varargin);

%Flag the elements that has not been flagged as filltype
if strcmpi(filltype,'hutter'),
	hutterflag(find(~macayealflag & ~pattynflag))=1;
elseif strcmpi(filltype,'macayeal'),
	macayealflag(find(~hutterflag & ~pattynflag))=1;
elseif strcmpi(filltype,'pattyn'),
	pattynflag(find(~hutterflag & ~macayealflag))=1;
end

%check that each element has at least one flag
if any(hutterflag+ macayealflag+pattynflag==0),
	error('setelementstype error message: elements type not assigned, must be specified')
end

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

%Check that no pattyn or stokes for 2d mesh
if strcmpi(md.type,'2d'),
	if any(stokesflag | pattynflag)
		error('setelementstype error message: stokes and pattyn elements no allowed in 2d mesh, extrude it first')
	end
end

%Modify stokesflag to get rid of elements contrained everywhere
nonstokesflag=find(~stokesflag);                   %non stokes elements
gridonstokes=ones(md.numberofgrids,1);
gridonstokes(md.elements(nonstokesflag,:))=0;      %non stokes grids
pos=find(md.gridondirichlet_diag);                 %find all the grids on the boundary of the domain without icefront
gridonstokes(pos)=0;                               %we have to constrain all the boundary of the domain without icefront

%Rule out elements where all grids are on the border of the stokes grids (because those elements are completely 
%constrained 
pos=find(stokesflag);
constrained_list=pos(find(sum(gridonstokes(md.elements(pos,:)),2)==0));
stokesflag(constrained_list)=0;

%add in model who is who
md.elements_type=zeros(md.numberofelements,2);

%1: Hutter elements
gridonhutter=zeros(md.numberofgrids,1);
gridonhutter(md.elements(find(hutterflag),:))=1;
md.gridonhutter=gridonhutter;
md.elements_type(find(hutterflag),1)=hutterenum();

%2: MacAyeal elements
gridonmacayeal=zeros(md.numberofgrids,1);
gridonmacayeal(md.elements(find(macayealflag),:))=1;
md.gridonmacayeal=gridonmacayeal;
md.elements_type(find(macayealflag),1)=macayealenum();

%3: Pattyn elements
gridonpattyn=zeros(md.numberofgrids,1);
gridonpattyn(md.elements(find(pattynflag),:))=1;
md.gridonpattyn=gridonpattyn;
md.elements_type(find(pattynflag),1)=pattynenum();

%4: Stokes elements
md.gridonstokes=gridonstokes;
md.elements_type(find(stokesflag),2)=stokesenum();

%5: None elements (non Stokes)
md.elements_type(find(~stokesflag),2)=noneenum();

%Create the border grids between Pattyn and MacAyeal and extrude them
numgrids2d=md.numberofgrids2d;
numlayers=md.numlayers;
bordergrids2d=find(gridonpattyn(1:numgrids2d) & gridonmacayeal(1:numgrids2d)); %Grids connected to two different types of elements

%initialize and fill in penalties structure
md.penalties=[];
if ~isnan(bordergrids2d),
	penalties=[];
	for	i=1:numlayers,
		penalties=[penalties [bordergrids2d+md.numberofgrids2d*(i-1)]];
	end
	md.penalties=penalties;
end

%flag dead grids (strictly in MacAyeal and not on bed -> not used in diagnostic horiz)
nonmacayeal_el=find(~macayealflag);                   %non macayeal elements
deadgrids=ones(md.numberofgrids,1);
deadgrids(md.elements(nonmacayeal_el,:))=0;          %non macayeal grids are not dead
deadgrids(find(md.gridonbed))=0;                      %grid on bed are not dead
md.deadgrids=deadgrids;
end
