function md2=modelextract(md1,area)
%MODELEXTRACT - extract a model according to an Argus contour or flag list
%
%   This routine extracts a submodel from a bigger model with respect to a given contour
%   md 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 
%   modelextractd, 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:
%      md2=modelextract(md1,area);
%
%   Examples:
%      md2=modelextract(md,'Domain.exp');
%      md2=modelextract(md,md.elementoniceshelf);
%
%   See also: EXTRUDE, COLLAPSE

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

%first
if ischar(area),
	if isempty(area),
		flag_elem=zeros(md1.numberofelements,1);
		invert=0;
	elseif strcmpi(area,'all')
		flag_elem=ones(md1.numberofelements,1);
		invert=0;
	else
		%make sure that we actually don't want the elements outside the domain outline!
		if strcmpi(area(1),'~'),
			area=area(2:length(area));
			invert=1;
		else
			invert=0;
		end
		%ok, flag_elem elements
		A=expread(area,1);
		flag_elem=ArgusContourToMesh(md1.elements(:,1:3),md1.x,md1.y,A,'element',1);
		%check that the elements do not cross the contour -> flag as 0
		pos=find(flag_elem);
		for contour=1:length(A)
			for i=1: length(pos)
				test=cross_extract(md1.x(md1.elements(pos(i),:)),md1.y(md1.elements(pos(i),:)),A(contour).x,A(contour).y);
				if test
					flag_elem(pos(i))=0;
				end
			end
		end
	end
	if invert,
		flag_elem=~flag_elem;
	end
elseif isfloat(area),
	if size(area,1)~=md1.numberofelements,
		setelementstypeusage();
		error('Flags must be of same size as number of elements in model');
	end
	flag_elem=area;
else
	error('Invalide option');
end
pos_elem=find(flag_elem);

%get the flags of grids
flag_grid=ismember([1:md1.numberofgrids]',sort(unique(md1.elements(pos_elem,:)))');
pos_grid=find(flag_grid);

%keep track of some fields
numberofgrids1=md1.numberofgrids;
numberofelements1=md1.numberofelements;
numberofgrids2=length(pos_grid);
numberofelements2=length(pos_elem);

%Create Pelem and Pgrid (transform old grids in new grids and same thing for the elements)
Pelem=zeros(numberofelements1,1);
Pelem(pos_elem)=[1:numberofelements2]';
Pgrid=zeros(numberofgrids1,1);
Pgrid(pos_grid)=[1:numberofgrids2]';

%renumber the elements (some grid won't exist anymore)
elements_1=md1.elements;
elements_2=elements_1(pos_elem,:);

elements_2(:,1)=Pgrid(elements_2(:,1));
elements_2(:,2)=Pgrid(elements_2(:,2));
elements_2(:,3)=Pgrid(elements_2(:,3));
if strcmpi(md1.type,'3d'),
	elements_2(:,4)=Pgrid(elements_2(:,4));
	elements_2(:,5)=Pgrid(elements_2(:,5));
	elements_2(:,6)=Pgrid(elements_2(:,6));
end

%OK, now create the new model !

	%take every fields from model
	md2=md1;

	%Mesh
	md2.numberofelements=numberofelements2;
	md2.numberofgrids=numberofgrids2;
	md2.elements=elements_2;
	if ~isnan(md2.elements_type), md2.elements_type=md1.elements_type(pos_elem,:);end
	md2.x=md1.x(pos_grid);
	md2.y=md1.y(pos_grid);
	md2.z=md1.z(pos_grid);
	if ~isnan(md2.bed_slopex), md2.bed_slopex=md1.bed_slopex(pos_grid); end;
	if ~isnan(md2.bed_slopey), md2.bed_slopex=md1.bed_slopey(pos_grid); end;
	if ~isnan(md2.surface_slopex), md2.surface_slopex=md1.surface_slopex(pos_grid); end;
	if ~isnan(md2.surface_slopey), md2.surface_slopex=md1.surface_slopey(pos_grid); end;

	%Initial 2d mesh 
	if strcmpi(md1.type,'3d')
		flag_elem_2d=flag_elem(1:md1.numberofelements2d);
		pos_elem_2d=find(flag_elem_2d);
		flag_grid_2d=flag_grid(1:md1.numberofgrids2d);
		pos_grid_2d=find(flag_grid_2d);

		md2.numberofelements2d=length(pos_elem_2d);
		md2.numberofgrids2d=length(pos_grid_2d);
		md2.elements2d=md1.elements2d(pos_elem_2d,:);
		md2.elements2d(:,1)=Pgrid(md2.elements2d(:,1));
		md2.elements2d(:,2)=Pgrid(md2.elements2d(:,2));
		md2.elements2d(:,3)=Pgrid(md2.elements2d(:,3));

		if ~isnan(md2.elements_type2d), md2.elements_type2d=md1.elements_type2d(pos_elem_2d,:); end;
		md2.x2d=md1.x(pos_grid_2d);
		md2.y2d=md1.y(pos_grid_2d);
		md2.z2d=md1.z(pos_grid_2d);
	end
		
	%Elements
	if ~isnan(md2.elementonhutter), md2.elementonhutter=md1.elementonhutter(pos_elem); end;
	if ~isnan(md2.elementonmacayeal), md2.elementonmacayeal=md1.elementonmacayeal(pos_elem); end;
	if ~isnan(md2.elementonpattyn), md2.elementonpattyn=md1.elementonpattyn(pos_elem); end;
	if ~isnan(md2.elementonstokes), md2.elementonstokes=md1.elementonstokes(pos_elem); end;

	%Nodes
	if ~isnan(md2.gridonhutter), md2.gridonhutter=md1.gridonhutter(pos_grid); end;
	if ~isnan(md2.gridonmacayeal), md2.gridonmacayeal=md1.gridonmacayeal(pos_grid); end;
	if ~isnan(md2.gridonpattyn), md2.gridonpattyn=md1.gridonpattyn(pos_grid); end;
	if ~isnan(md2.gridonstokes), md2.gridonstokes=md1.gridonstokes(pos_grid); end;
	
	%Penalties
	if ~isnan(md2.penalties),
		for i=1:size(md1.penalties,1);
			md2.penalties(i,:)=Pgrid(md1.penalties(i,:));
		end
		md2.penalties=md2.penalties(find(md2.penalties(:,1)),:);
	end
	md2.rifts=NaN; %????????????????????????????
	md2.numrifts=0;%????????????????????????????

	%Projections
	if ~isnan(md2.uppergrids), md2.uppergrids=Pgrid(md1.uppergrids(pos_grid)); end;
	if ~isnan(md2.lowergrids), md2.lowergrids=Pgrid(md1.lowergrids(pos_grid)); end;
	if ~isnan(md2.deadgrids), md2.deadgrids=md1.deadgrids(pos_grid);end;
	
	%Extrusion
	md2.elementonbed=md1.elementonbed(pos_elem);
	md2.elementonsurface=md1.elementonsurface(pos_elem);
	md2.gridonbed=md1.gridonbed(pos_grid);
	md2.gridonsurface=md1.gridonsurface(pos_grid);
	if ~isnan(md2.firn_layer), md2.firn_layer=md1.firn_layer(pos_grid); end;
	
	%Physical parameters
	if ~isnan(md2.drag), md2.drag=md1.drag(pos_grid);end;
	if ~isnan(md2.p), md2.p=md1.p(pos_elem);end;
	if ~isnan(md2.q), md2.q=md1.q(pos_elem);end;
	if ~isnan(md2.n), md2.n=md1.n(pos_elem);end;
	if ~isnan(md2.B), md2.B=md1.B(pos_grid);end;

	%Geometrical parameters
	if ~isnan(md2.elementoniceshelf), md2.elementoniceshelf=md1.elementoniceshelf(pos_elem);end;
	if ~isnan(md2.elementonicesheet), md2.elementonicesheet=md1.elementonicesheet(pos_elem);end;
	if ~isnan(md2.gridoniceshelf), md2.gridoniceshelf=md1.gridoniceshelf(pos_grid);end;
	if ~isnan(md2.gridonicesheet), md2.gridonicesheet=md1.gridonicesheet(pos_grid);end;
	if ~isnan(md2.surface), md2.surface=md1.surface(pos_grid); end;
	if ~isnan(md2.thickness), md2.thickness=md1.thickness(pos_grid); end;
	if ~isnan(md2.new_thickness), md2.new_thickness=md1.new_thickness(pos_grid); end;
	if ~isnan(md2.bed), md2.bed=md1.bed(pos_grid); end;

	%Boundary conditions
	md2.gridonboundary=md1.gridonboundary(pos_grid);

	%Diagnostic
	if ~isnan(md2.segmentonneumann_diag)
		md2.segmentonneumann_diag(:,1)=Pgrid(md1.segmentonneumann_diag(:,1)); 
		md2.segmentonneumann_diag(:,2)=Pgrid(md1.segmentonneumann_diag(:,2)); 
		md2.segmentonneumann_diag(:,end)=Pelem(md1.segmentonneumann_diag(:,end)); 
		if strcmpi(md1.type,'3d')
			md2.segmentonneumann_diag(:,3)=Pgrid(md1.segmentonneumann_diag(:,3)); 
			md2.segmentonneumann_diag(:,4)=Pgrid(md1.segmentonneumann_diag(:,4)); 
		end
		md2.segmentonneumann_diag=md2.segmentonneumann_diag(find(md2.segmentonneumann_diag(:,1) & md2.segmentonneumann_diag(:,2)),:);
	end
	md2.neumannvalues_diag=NaN; %???????????????????????????
	if ~isnan(md2.gridondirichlet_diag), md2.gridondirichlet_diag=md1.gridondirichlet_diag(pos_grid); end;
	if ~isnan(md2.dirichletvalues_diag),
		md2.dirichletvalues_diag=[];
		for i=1:size(md1.dirichletvalues_diag,2)
			md2.dirichletvalues_diag(:,i)=md1.dirichletvalues_diag(pos_grid,i);
		end
	end

	%Thermal
	if ~isnan(md2.gridondirichlet_thermal), md2.gridondirichlet_thermal=md1.gridondirichlet_thermal(pos_grid); end;
	if ~isnan(md2.dirichletvalues_thermal), md2.dirichletvalues_thermal=md1.dirichletvalues_thermal(pos_grid); end;

	%Transient
	if ~isnan(md2.segmentonneumann_prog)
		md2.segmentonneumann_prog(:,1)=Pgrid(md1.segmentonneumann_prog(:,1)); 
		md2.segmentonneumann_prog(:,2)=Pgrid(md1.segmentonneumann_prog(:,2)); 
		md2.segmentonneumann_prog(:,end)=Pelem(md1.segmentonneumann_prog(:,end)); 
		md2.segmentonneumann_prog=md2.segmentonneumann_prog(find(md2.segmentonneumann_prog(:,1) & md2.segmentonneumann_prog(:,2)),:);
	end
	md2.neumannvalues_prog=NaN; %???????????????????????????
	if ~isnan(md2.segmentonneumann_prog2)
		md2.segmentonneumann_prog2(:,1)=Pgrid(md1.segmentonneumann_prog2(:,1)); 
		md2.segmentonneumann_prog2(:,2)=Pgrid(md1.segmentonneumann_prog2(:,2)); 
		md2.segmentonneumann_prog2(:,end)=Pelem(md1.segmentonneumann_prog2(:,end)); 
		md2.segmentonneumann_prog2=md2.segmentonneumann_prog2(find(md2.segmentonneumann_prog2(:,1) & md2.segmentonneumann_prog2(:,2)),:);
	end
	md2.neumannvalues_prog=NaN; %???????????????????????????
	if ~isnan(md2.gridondirichlet_prog), md2.gridondirichlet_prog=md1.gridondirichlet_prog(pos_grid); end;
	if ~isnan(md2.dirichletvalues_prog), md2.dirichletvalues_prog=md1.dirichletvalues_prog(pos_grid); end;

	%Observations
	if ~isnan(md2.vx_obs), md2.vx_obs=md1.vx_obs(pos_grid); end;
	if ~isnan(md2.vy_obs), md2.vy_obs=md1.vy_obs(pos_grid); end;
	if ~isnan(md2.vel_obs), md2.vel_obs=md1.vel_obs(pos_grid); end;
	if ~isnan(md2.accumulation), md2.accumulation=md1.accumulation(pos_grid); end;
	if ~isnan(md2.geothermalflux), md2.geothermalflux=md1.geothermalflux(pos_grid); end;
	if ~isnan(md2.observed_temperature), md2.observed_temperature=md1.observed_temperature(pos_grid); end;
	
	%Output parameters
	md2.viscousheating=NaN;
	md2.pressure_elem=NaN;
	md2.stress=NaN;
	md2.stress_surface=NaN;
	md2.stress_bed=NaN;
	md2.deviatoricstress=NaN;
	md2.strainrate=NaN;

	%Results fields
	if ~isnan(md2.vx), md2.vx=md1.vx(pos_grid); end;
	if ~isnan(md2.vy), md2.vy=md1.vy(pos_grid); end;
	if ~isnan(md2.vz), md2.vz=md1.vz(pos_grid); end;
	if ~isnan(md2.vel), md2.vel=md1.vel(pos_grid); end;
	if ~isnan(md2.temperature), md2.temperature=md1.temperature(pos_grid); end;
	if ~isnan(md2.melting), md2.melting=md1.melting(pos_grid); end;
	if ~isnan(md2.pressure), md2.pressure=md1.pressure(pos_grid); end;

%Now take care of the boudaries (Dirichlet).

%Catch the elements that have not been extracted
orphans_elem=find(~flag_elem);
orphans_grid=unique(md1.elements(orphans_elem,:))';

%Figure out which grid are on the boundary between md2 and md1
gridstoflag1=intersect(orphans_grid,pos_grid);
gridstoflag2=Pgrid(gridstoflag1);

%Same thing for the elements
flaglist=ismember(md1.elements,gridstoflag1);
elementstoflag1=find((flaglist(:,1) | flaglist(:,2) | flaglist(:,3)) & ~flag_elem);
flaglist=ismember(md2.elements,gridstoflag2);
elementstoflag2=find(flaglist(:,1) | flaglist(:,2) | flaglist(:,3));

%These grids are "on boundary" and hence "on dirichlet"
test=0;
if ~isnan(md2.gridonboundary), md2.gridonboundary(gridstoflag2)=1; end;
if ~isnan(md2.gridondirichlet_diag), md2.gridondirichlet_diag(gridstoflag2)=1; end;
if ~isnan(md2.gridondirichlet_prog), md2.gridondirichlet_prog(gridstoflag2)=1; end;
if ~isnan(md2.gridondirichlet_thermal), md2.gridondirichlet_thermal(gridstoflag2)=1; end;
if ~isnan(md2.dirichletvalues_diag), 
	test=1; 
	if strcmpi(md1.type,'3d')
		if ~isnan(md2.vx) & ~isnan(md2.vy)
			md2.dirichletvalues_diag(gridstoflag2,1)=md2.vx(gridstoflag2); 
			md2.dirichletvalues_diag(gridstoflag2,2)=md2.vy(gridstoflag2);
		end
	else
		if ~isnan(md2.vx_obs) & ~isnan(md2.vy_obs)
			md2.dirichletvalues_diag(gridstoflag2,1)=md2.vx_obs(gridstoflag2); 
			md2.dirichletvalues_diag(gridstoflag2,2)=md2.vy_obs(gridstoflag2);
		end
	end
end
if ~isnan(md2.dirichletvalues_prog), test=1; end;
if ~isnan(md2.dirichletvalues_thermal), test=1; end;

if test
	disp(' ')
	disp('!! modelextract warning: dirichlet values should be checked !!')
	disp(' ')
end

%Keep track of pos_grid and pos_elem
md2.extractedgrids=pos_grid;
md2.extractedelements=pos_elem;

%get segments
elementsonboundary=find(sum(md2.gridonboundary(md2.elements),2)>=2);
md2.segments=[];
for i=1:length(elementsonboundary),
	el=elementsonboundary(i);
	for j=1:3,
		grid1=md2.elements(el,j);
		if (j==3),
			grid2=md2.elements(el,1);
		else
			grid2=md2.elements(el,j+1);
		end
		%ok, [grid1,grid2] are a segment only if they belong to one element
		pos=find(    ((md2.elements(:,1)==grid1) | (md2.elements(:,2)==grid1) | (md2.elements(:,3)==grid1)) & ...
		            ((md2.elements(:,1)==grid2) | (md2.elements(:,2)==grid2) | (md2.elements(:,3)==grid2))...
				);
	    if length(pos)==1,
			md2.segments=[md2.segments; [grid1 grid2 pos]];
		end
	end
end
md2.segmentmarkers=ones(size(md2.segments,1));


end %function

function bool=cross_extract(x1,y1,x2,y2);
%return 1 if the polygon defined by (x1,y1) crosses the polygon defined by (x2,y2)
	x1 = x1(:); y1 = y1(:);
	x2 = x2(:); y2 = y2(:);
	l1 = length(x1);
	l2 = length(x2);

	% Indexing
	[i11,i12] = meshgrid(1:l1,1:l2);
	[i21,i22] = meshgrid([2:l1 1],[2:l2 1]);
	i11 = i11(:); i12 = i12(:);
	i21 = i21(:); i22 = i22(:);

	% Calculate matrix of intersections
	S = iscross([x1(i11) x1(i21)]',[y1(i11) y1(i21)]',...
		[x2(i12) x2(i22)]',[y2(i12) y2(i22)]')';
	S(find(S==0.5))=0;

	bool=any(any(S));
end
