[2737] | 1 | function md=meshnodensity(md,domainname,varargin)
|
---|
| 2 | %MESH - create model mesh
|
---|
| 3 | %
|
---|
| 4 | % This routine creates a model mesh using TriMeshNoDensity and a domain outline
|
---|
| 5 | % where md is a @model object, domainname is the name of an Argus domain outline file,
|
---|
| 6 | % Riftname is an optional argument (Argus domain outline) describing rifts.
|
---|
| 7 | % The difference with mesh.m is that the resolution of the mesh follows that of the domain
|
---|
| 8 | % outline and the riftoutline
|
---|
| 9 | %
|
---|
| 10 | % Usage:
|
---|
| 11 | % md=mesh(md,domainname)
|
---|
| 12 | % or md=mesh(md,domainname,riftname);
|
---|
| 13 | %
|
---|
| 14 | % Examples:
|
---|
| 15 | % md=mesh(md,'DomainOutline.exp');
|
---|
| 16 | % md=mesh(md,'DomainOutline.exp','Rifts.exp');
|
---|
| 17 |
|
---|
| 18 | if (nargin==2),
|
---|
| 19 | riftname='';
|
---|
| 20 | end
|
---|
| 21 | if (nargin==3),
|
---|
| 22 | riftname=varargin{1};
|
---|
| 23 | end
|
---|
| 24 |
|
---|
| 25 | %Mesh using TriMeshNoDensity
|
---|
| 26 | if strcmp(riftname,''),
|
---|
| 27 | [md.elements,md.x,md.y,md.segments,md.segmentmarkers]=TriMeshNoDensity(domainname);
|
---|
| 28 | else
|
---|
| 29 | [elements,x,y,segments,segmentmarkers]=TriMeshNoDensity(domainname,riftname);
|
---|
| 30 |
|
---|
[8298] | 31 | %check that all the created nodes belong to at least one element
|
---|
[2737] | 32 | orphan=find(~ismember([1:length(x)],sort(unique(elements(:)))));
|
---|
| 33 | for i=1:length(orphan),
|
---|
[8298] | 34 | %get rid of the orphan node i
|
---|
[2737] | 35 | %update x and y
|
---|
| 36 | x=[x(1:orphan(i)-(i-1)-1); x(orphan(i)-(i-1)+1:end)];
|
---|
| 37 | y=[y(1:orphan(i)-(i-1)-1); y(orphan(i)-(i-1)+1:end)];
|
---|
| 38 | %update elements
|
---|
| 39 | pos=find(elements>orphan(i)-(i-1));
|
---|
| 40 | elements(pos)=elements(pos)-1;
|
---|
| 41 | %update segments
|
---|
| 42 | pos1=find(segments(:,1)>orphan(i)-(i-1));
|
---|
| 43 | pos2=find(segments(:,2)>orphan(i)-(i-1));
|
---|
| 44 | segments(pos1,1)=segments(pos1,1)-1;
|
---|
| 45 | segments(pos2,2)=segments(pos2,2)-1;
|
---|
| 46 | end
|
---|
| 47 |
|
---|
| 48 | %plug into md
|
---|
| 49 | md.x=x;
|
---|
| 50 | md.y=y;
|
---|
| 51 | md.elements=elements;
|
---|
| 52 | md.segments=segments;
|
---|
| 53 | md.segmentmarkers=segmentmarkers;
|
---|
| 54 | end
|
---|
| 55 |
|
---|
| 56 | %Fill in rest of fields:
|
---|
| 57 | md.numberofelements=length(md.elements);
|
---|
[8298] | 58 | md.numberofnodes=length(md.x);
|
---|
| 59 | md.z=zeros(md.numberofnodes,1);
|
---|
| 60 | md.nodeonboundary=zeros(md.numberofnodes,1); md.nodeonboundary(md.segments(:,1:2))=1;
|
---|
| 61 | md.nodeonbed=ones(md.numberofnodes,1);
|
---|
| 62 | md.nodeonsurface=ones(md.numberofnodes,1);
|
---|
[2737] | 63 | md.elementonbed=ones(md.numberofelements,1);
|
---|
| 64 | md.elementonsurface=ones(md.numberofelements,1);
|
---|
| 65 |
|
---|
| 66 | %Now, build the connectivity tables for this mesh.
|
---|
[8298] | 67 | md.nodeconnectivity=NodeConnectivity(md.elements,md.numberofnodes);
|
---|
[2737] | 68 | md.elementconnectivity=ElementConnectivity(md.elements,md.nodeconnectivity);
|
---|
| 69 |
|
---|
| 70 | %outline names
|
---|
| 71 | md.domainoutline=char(textread(domainname,'%s','delimiter','\n'));
|
---|
| 72 | if strcmp(riftname,''),
|
---|
| 73 | md.riftoutline='';
|
---|
| 74 | else
|
---|
| 75 | md.riftoutline=readfile(riftname);
|
---|
| 76 | end
|
---|
| 77 |
|
---|
| 78 | %type of model
|
---|
[3994] | 79 | md.dim=2;
|
---|
[2737] | 80 |
|
---|
| 81 | %augment counter keeping track of what has been done to this model
|
---|
| 82 | md.counter=1;
|
---|