function md2=removeholes(md,field)
%REMOVEHOLES - interpolate a field on the a mesh without any hole
%
%   as its name indicates, this routine takes a model, a field (value of some 
%   physical quantity evaluated at the mesh grids) of the model, and interpolates this field 
%   on the model mesh, without any hole.
%
%   Usage:
%      md=removeholes(md,field)

%Check that model is complete
if md.counter<3,
	error('removeholes error message: model is incomplete ... exiting');
end

if nargin~=2,
	removeholesusage;
	error('removeholes error message');
end

if ~ischar(field), 
	removeholesusage;
	error('removeholes error message');
end


%Ok, retrieve and write domain outline without holes, to disk.
domainoutline_string=md.domainoutline;
name_slots=findstr(domainoutline_string,'## Name');
domainoutline_string=domainoutline_string(1:(name_slots(2)-1)); %only keep first outline
writefile('DomainOutlineTemp.exp',domainoutline_string);

%Now create new model with mesh based on DomainOutlineTemp: 
%get average resolution
resolution=mean(sqrt(2*area(md)));
md2=model;
md2=mesh(md2,'DomainOutlineTemp.exp',resolution);

%Ok, now interpolate field onto this new mesh
fieldvalue=getfield(md,field);
md2=setfield(md2,field,griddata(md.x,md.y,fieldvalue,md2.x,md2.y));

end %end of function

function removeholesusage(),
disp('usage: md2=removeholes(md,field)');
end
