function md2=mesh2(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:
%      md2=mesh2(md,domainname,resolution)
%   or md2=mesh2(md,domainname,riftname, resolution)
%
%   Examples:
%      md2=mesh2(md,'DOmainOutline.exp',1000);
%      md2=mesh2(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;
%Initialize return model;
md2=model;

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

md2=subsasgn(md2,struct('type','.','subs','elements'),elements); 
md2=subsasgn(md2,struct('type','.','subs','segmentmarkers'),segmentmarkers); 
md2=subsasgn(md2,struct('type','.','subs','segments'),segments); 
md2=subsasgn(md2,struct('type','.','subs','y'),y);
md2=subsasgn(md2,struct('type','.','subs','x'),x); 

%Fill in rest of fields:
md2=subsasgn(md2,struct('type','.','subs','numberofelements'),length(elements));
md2=subsasgn(md2,struct('type','.','subs','numberofgrids'),length(x));
md2=subsasgn(md2,struct('type','.','subs','z'),zeros(length(x),1));
gridonboundary=zeros(length(x),1); gridonboundary(segments(:,1:2))=1;
md2=subsasgn(md2,struct('type','.','subs','gridonboundary'),gridonboundary);
md2=subsasgn(md2,struct('type','.','subs','elements_type'),3*ones(md2.numberofelements,1)); % type determined by number of grids per element

%outline names
md2=subsasgn(md2,struct('type','.','subs','domainoutline'),readfile(domainname));
if strcmp(riftname,''),
	md2=subsasgn(md2,struct('type','.','subs','riftoutline'),'');
else
	md2=subsasgn(md2,struct('type','.','subs','riftoutline'),readfile(riftname));
end

%type of model
md2=subsasgn(md2,struct('type','.','subs','type'),'2d');
	
%augment counter  keeping track of what has been done to this model
md2=subsasgn(md2,struct('type','.','subs','counter'),1);
