function md=meshprocessrifts(md)
%MESHPROCESSRIFTS - process mesh when rifts are present
%
%   split rifts inside mesh (rifts are defined by presence of
%   segments inside the domain outline)
%
%   Usage:
%      md=meshprocessrifts(md)

%some checks on list of arguments
if ((nargin~=1) | (nargout~=1)),
	meshprocessriftsusage();
	error('meshprocessrifts error message');
end

[md.elements,md.x,md.y,md.segments,md.segmentmarkers,md.rifts]=TriMeshProcessRifts(md.elements,md.x,md.y,md.segments,md.segmentmarkers);

%Fill in rest of fields:
md.numberofelements=length(md.elements);
md.numberofgrids=length(md.x);
md.z=zeros(md.numberofgrids,1);
md.gridonboundary=zeros(length(md.x),1); md.gridonboundary(md.segments(:,1:2))=1;
md.numrifts=length(md.rifts);
md.elements_type=3*ones(md.numberofelements,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);

function meshprocessriftsusage(),
disp('usage: md=meshprocessrifts(md)');
