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 
%   modelextract2d, 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
%   add an argument 0 if you do not want the elements to be checked (faster)
%
%   Usage:
%      md2=modelextract(md1,area);
%
%   Examples:
%      md2=modelextract(md,'Domain.exp');
%      md2=modelextract(md,md.elementoniceshelf);
%
%   See also: EXTRUDE, COLLAPSE

%some checks
if ((nargin~=2) | (nargout~=1)),
	help modelextract
	error('modelextract error message: bad usage');
end

%get check option
if (nargin==3 & varargin{1}==0),
	checkoutline=0;
else
	checkoutline=1;
end

%first
flag_elem=FlagElements(md1,area);

%kick out all elements with 3 dirichlets
spc_elem=find(~flag_elem);
spc_grid=sort(unique(md1.elements(spc_elem,:)));
flag=ones(md1.numberofgrids,1);
flag(spc_grid)=0;
pos=find(sum(flag(md1.elements),2)==0);
flag_elem(pos)=0;

%extracted elements and nodes lists
pos_elem=find(flag_elem);
pos_grid=sort(unique(md1.elements(pos_elem,:)));

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

%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 md1.dim==3,
	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;

	%automatically modify fields

	%loop over model fields
	model_fields=fields(md1);
	for i=1:length(model_fields),

		%get field
		field=md1.(model_fields(i));
		fieldsize=size(field);

		%size = number of grids * n
		if fieldsize(1)==numberofgrids1
			md2.(model_fields(i))=field(pos_grid,:);

		%size = number of elements * n
		elseif fieldsize(1)==numberofelements1
			md2.(model_fields(i))=field(pos_elem,:);
		end
	end

	%modify some specific fields

	%Mesh
	md2.numberofelements=numberofelements2;
	md2.numberofgrids=numberofgrids2;
	md2.elements=elements_2;

	%uppernodes lowernodes
	if md1.dim==3
		md2.uppergrids=md1.uppergrids(pos_grid);
		pos=find(~isnan(md2.uppergrids));
		md2.uppergrids(pos)=Pgrid(md2.uppergrids(pos));

		md2.lowergrids=md1.lowergrids(pos_grid);
		pos=find(~isnan(md2.lowergrids));
		md2.lowergrids(pos)=Pgrid(md2.lowergrids(pos));

		md2.upperelements=md1.upperelements(pos_elem);
		pos=find(~isnan(md2.upperelements));
		md2.upperelements(pos)=Pelem(md2.upperelements(pos));

		md2.lowerelements=md1.lowerelements(pos_elem);
		pos=find(~isnan(md2.lowerelements));
		md2.lowerelements(pos)=Pelem(md2.lowerelements(pos));
	end

	%Initial 2d mesh 
	if md1.dim==3
		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

	%Edges
	if size(md2.edges,2)>1, %do not use ~isnan because there are some NaNs...
		%renumber first two columns
		pos=find(~isnan(md2.edges(:,4)));
		md2.edges(:  ,1)=Pgrid(md2.edges(:,1)); 
		md2.edges(:  ,2)=Pgrid(md2.edges(:,2)); 
		md2.edges(:  ,3)=Pelem(md2.edges(:,3));
		md2.edges(pos,4)=Pelem(md2.edges(pos,4));
		%remove edges when the 2 vertices are not in the domain.
		md2.edges=md2.edges(find(md2.edges(:,1) & md2.edges(:,2)),:);
		%Replace all zeros by NaN in the last two columns;
		pos=find(md2.edges(:,3)==0);
		md2.edges(pos,3)=NaN;
		pos=find(md2.edges(:,4)==0);
		md2.edges(pos,4)=NaN;
		%Invert NaN of the third column with last column (Also invert first two columns!!)
		pos=find(isnan(md2.edges(:,3)));
		md2.edges(pos,3)=md2.edges(pos,4);
		md2.edges(pos,4)=NaN;
		values=md2.edges(pos,2);
		md2.edges(pos,2)=md2.edges(pos,1);
		md2.edges(pos,1)=values;
		%Finally remove edges that do not belong to any element
		pos=find(isnan(md2.edges(:,3)) & isnan(md2.edges(:,4)));
		md2.edges(pos,:)=[];
	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

	%recreate segments
	if md1.dim==2
		md2.nodeconnectivity=NodeConnectivity(md2.elements,md2.numberofgrids);
		md2.elementconnectivity=ElementConnectivity(md2.elements,md2.nodeconnectivity);
		md2.segments=contourenvelope(md2);
		md2.gridonboundary=zeros(numberofgrids2,1); md2.gridonboundary(md2.segments(:,1:2))=1;
	end

	%Boundary conditions: Dirichlets on new boundary
	%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);
	if ~isnan(md1.spcvelocity),
		md2.spcvelocity(gridstoflag2,1:3)=1;
		if ~isnan(md1.vx_obs) & ~isnan(md1.vy_obs)
			md2.spcvelocity(gridstoflag2,4)=md2.vx_obs(gridstoflag2); 
			md2.spcvelocity(gridstoflag2,5)=md2.vy_obs(gridstoflag2);
		else
			md2.spcvelocity(gridstoflag2,4:5)=zeros(length(gridstoflag2),2);
			disp(' ')
			disp('!! modelextract warning: spc values should be checked !!')
			disp(' ')
		end
	end
	if ~isnan(md1.spctemperature),
		md2.spctemperature(gridstoflag2,1)=1;
		if ~isnan(md1.observed_temperature)
			md2.spctemperature(gridstoflag2,2)=md2.observed_temperature(gridstoflag2); 
		else
			md2.spctemperature(gridstoflag2,2)=zeros(length(gridstoflag2),2);
			disp(' ')
			disp('!! modelextract warning: spc values should be checked !!')
			disp(' ')
		end
	end

	%Diagnostic
	if ~isnan(md2.pressureload)
		md2.pressureload(:,1)=Pgrid(md1.pressureload(:,1)); 
		md2.pressureload(:,2)=Pgrid(md1.pressureload(:,2)); 
		md2.pressureload(:,end-1)=Pelem(md1.pressureload(:,end-1));
		if md1.dim==3
			md2.pressureload(:,3)=Pgrid(md1.pressureload(:,3)); 
			md2.pressureload(:,4)=Pgrid(md1.pressureload(:,4)); 
		end
		md2.pressureload=md2.pressureload(find(md2.pressureload(:,1) & md2.pressureload(:,2) & md2.pressureload(:,end)),:);
	end

	%Results fields
	if isstruct(md1.results),
		md2.results=struct();
		solutionfields=fields(md1.results);
		for i=1:length(solutionfields),
			%get subfields
			solutionsubfields=fields(md1.results.(solutionfields(i)));
			for j=1:length(solutionsubfields),
				field=md1.results.(solutionfields(i)).(solutionsubfields(j));
				if length(field)==numberofgrids1,
					md2.results.(solutionfields(i)).(solutionsubfields(j))=field(pos_grid);
				elseif length(field)==numberofelements1,
					md2.results.(solutionfields(i)).(solutionsubfields(j))=field(pos_elem);
				else
					md2.results.(solutionfields(i)).(solutionsubfields(j))=field;
				end
			end
		end
	end

	%reinitialize 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;

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