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;