source: issm/trunk/src/m/model/mesh/rifts/meshprocessrifts.m@ 9733

function md=meshprocessrifts(md,domainoutline)
%MESHPROCESSRIFTS - process mesh when rifts are present
%
%   split rifts inside mesh (rifts are defined by presence of
%   segments inside the domain outline)
%   if domain outline is provided, check for rifts that could touch it, and open them up.
%
%   Usage:
%      md=meshprocessrifts(md,domainoutline)
%
%   Ex: 
%      md=meshprocessrifts(md,'DomainOutline.exp');
%

%some checks on arguments: 
if nargout~=1,
        help meshprocessrifts
        error('meshprocessrifts usage error:');
end

if nargin~=2,
        help meshprocessrifts
        error('meshprocessrifts usage error:');
end

%Call MEX file
[md.mesh.elements,md.x,md.y,md.mesh.segments,md.mesh.segmentmarkers,md.rifts.riftstruct]=TriMeshProcessRifts(md.mesh.elements,md.x,md.y,md.mesh.segments,md.mesh.segmentmarkers);
md.rifts
md.rifts.riftstruct
if ~isstruct(md.rifts.riftstruct),
        error('TriMeshProcessRifts did not find any rift');
end

%Fill in rest of fields:
md.mesh.numberofelements=length(md.mesh.elements);
md.mesh.numberofvertices=length(md.x);
md.z=zeros(md.mesh.numberofvertices,1);
md.mesh.vertexonboundary=zeros(length(md.x),1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
md.rifts.numrifts=length(md.rifts.riftstruct);
md.flowequation.element_equation=3*ones(md.mesh.numberofelements,1);
md.mesh.vertexonbed=ones(md.mesh.numberofvertices,1);
md.mesh.vertexonsurface=ones(md.mesh.numberofvertices,1);
md.mesh.elementonbed=ones(md.mesh.numberofelements,1);
md.mesh.elementonsurface=ones(md.mesh.numberofelements,1);

%get coordinates of rift tips
for i=1:md.rifts.numrifts,
        md.rifts.riftstruct(i).tip1coordinates=[md.x(md.rifts.riftstruct(i).tips(1)) md.y(md.rifts.riftstruct(i).tips(1))];
        md.rifts.riftstruct(i).tip2coordinates=[md.x(md.rifts.riftstruct(i).tips(2)) md.y(md.rifts.riftstruct(i).tips(2))];
end

%In case we have rifts that open up the domain outline, we need to open them: 
flags=ContourToMesh(md.mesh.elements,md.x,md.y,domainoutline,'node',0);
found=0;
for i=1:md.rifts.numrifts,
        if flags(md.rifts.riftstruct(i).tips(1))==0,
                found=1;
                break;
        end
        if flags(md.rifts.riftstruct(i).tips(2))==0,
                found=1;
                break;
        end
end
if found,
        md=meshprocessoutsiderifts(md,domainoutline);
end

%get elements that are not correctly oriented in the correct direction:
aires=GetAreas(md.mesh.elements,md.x,md.y);
pos=find(aires<0);
md.mesh.elements(pos,:)=[md.mesh.elements(pos,2) md.mesh.elements(pos,1) md.mesh.elements(pos,3)];