source: issm/trunk/src/m/model/mesh/meshnodensity.m@ 8298

function md=meshnodensity(md,domainname,varargin)
%MESH - create model mesh
%
%   This routine creates a model mesh using TriMeshNoDensity and a domain outline
%   where md is a @model object, domainname is the name of an Argus domain outline file, 
%   Riftname is an optional argument (Argus domain outline) describing rifts.
%   The  difference with mesh.m is that the resolution of the mesh follows that of the domain 
%   outline and the riftoutline
%
%   Usage:
%      md=mesh(md,domainname)
%   or md=mesh(md,domainname,riftname);
%
%   Examples:
%      md=mesh(md,'DomainOutline.exp');
%      md=mesh(md,'DomainOutline.exp','Rifts.exp');

if (nargin==2),
        riftname='';
end
if (nargin==3),
        riftname=varargin{1};
end

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

        %check that all the created nodes 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 node 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.numberofnodes=length(md.x);
md.z=zeros(md.numberofnodes,1);
md.nodeonboundary=zeros(md.numberofnodes,1); md.nodeonboundary(md.segments(:,1:2))=1;
md.nodeonbed=ones(md.numberofnodes,1);
md.nodeonsurface=ones(md.numberofnodes,1);
md.elementonbed=ones(md.numberofelements,1);
md.elementonsurface=ones(md.numberofelements,1);

%Now, build the connectivity tables for this mesh.
md.nodeconnectivity=NodeConnectivity(md.elements,md.numberofnodes);
md.elementconnectivity=ElementConnectivity(md.elements,md.nodeconnectivity);

%outline names
md.domainoutline=char(textread(domainname,'%s','delimiter','\n'));
if strcmp(riftname,''),
        md.riftoutline='';
else
        md.riftoutline=readfile(riftname);
end

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