function md=mesh(md,domainname,varargin)
%MESH - create model mesh
%
%   This routine creates a model mesh using TriMesh and a domain outline, to within a certain resolution
%   where md is a @model object, domainname is the name of an Argus domain outline file, 
%   and resolution is a characteristic length for the mesh (same unit as the domain outline
%   unit). Riftname is an optional argument (Argus domain outline) describing rifts.
%
%   Usage:
%      md=mesh(md,domainname,resolution)
%   or md=mesh(md,domainname,riftname, resolution)
%
%   Examples:
%      md=mesh(md,'DOmainOutline.exp',1000);
%      md=mesh(md,'DOmainOutline.exp','Rifts.exp',1500);

%Figure out a characteristic area. Resolution is a grid oriented concept (ex a 1000m  resolution grid would 
%be made of 1000*1000 area squares). 
if (nargin==3),
	resolution=varargin{1};
	riftname='';
end
if (nargin==4),
	riftname=varargin{1};
	resolution=varargin{2};
end

%Check that mesh was not already run, and warn user: 
if subsref(md,struct('type','.','subs','counter'))>=1,
	choice=input('This model already has a mesh. Are you sure you want to go ahead? (y/n)','s');
	if ~strcmp(choice,'y')
		error('no meshing done ... exiting');
	end
end

area=resolution^2;

%Mesh using TriMesh
if strcmp(riftname,''),
	[md.elements,md.x,md.y,md.segments,md.segmentmarkers]=TriMesh(domainname,area,'yes');
else
	[elements,x,y,segments,segmentmarkers]=TriMesh(domainname,riftname,area,'yes');

	%check that all the created grids belong to at least one element
	orphan=find(~ismember([1:length(x)],sort(unique(elements(:)))));
	for i=1:length(orphan),
		%get rid of the orphan grid i
		%update x and y
		x=[x(1:orphan(i)-(i-1)-1); x(orphan(i)-(i-1)+1:end)];
		y=[y(1:orphan(i)-(i-1)-1); y(orphan(i)-(i-1)+1:end)];
		%update elements
		pos=find(elements>orphan(i)-(i-1));
		elements(pos)=elements(pos)-1;
		%update segments
		pos1=find(segments(:,1)>orphan(i)-(i-1));
		pos2=find(segments(:,2)>orphan(i)-(i-1));
		segments(pos1,1)=segments(pos1,1)-1;
		segments(pos2,2)=segments(pos2,2)-1;
	end

	%plug into md
	md.x=x;
	md.y=y;
	md.elements=elements;
	md.segments=segments;
	md.segmentmarkers=segmentmarkers;
end

%Fill in rest of fields:
md.numberofelements=length(md.elements);
md.numberofgrids=length(md.x);
md.z=zeros(md.numberofgrids,1);
md.gridonboundary=zeros(md.numberofgrids,1); md.gridonboundary(md.segments(:,1:2))=1;
md.gridonbed=ones(md.numberofgrids,1);
md.gridonsurface=ones(md.numberofgrids,1);
md.elementonbed=ones(md.numberofelements,1);
md.elementonsurface=ones(md.numberofelements,1);

%outline names
md.domainoutline=readfile(domainname);
if strcmp(riftname,''),
	md.riftoutline='';
else
	md.riftoutline=readfile(riftname);
end

%type of model
md.type='2d';
	
%augment counter  keeping track of what has been done to this model
md.counter=1;
