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=setmesh(md,domainname)
|
---|
12 | % or md=setmesh(md,domainname,riftname);
|
---|
13 | %
|
---|
14 | % Examples:
|
---|
15 | % md=setmesh(md,'DomainOutline.exp');
|
---|
16 | % md=setmesh(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.mesh.segments,md.mesh.segmentmarkers]=TriMeshNoDensity(domainname);
|
---|
28 | else
|
---|
29 | [elements,x,y,segments,segmentmarkers]=TriMeshNoDensity(domainname,riftname);
|
---|
30 |
|
---|
31 | %check that all the created nodes belong to at least one element
|
---|
32 | orphan=find(~ismember([1:length(x)],sort(unique(elements(:)))));
|
---|
33 | for i=1:length(orphan),
|
---|
34 | %get rid of the orphan node i
|
---|
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.mesh.segments=segments;
|
---|
53 | md.mesh.segmentmarkers=segmentmarkers;
|
---|
54 | end
|
---|
55 |
|
---|
56 | %Fill in rest of fields:
|
---|
57 | md.numberofelements=length(md.elements);
|
---|
58 | md.numberofnodes=length(md.x);
|
---|
59 | md.z=zeros(md.numberofnodes,1);
|
---|
60 | md.mesh.vertexonboundary=zeros(md.numberofnodes,1); md.mesh.vertexonboundary(md.mesh.segments(:,1:2))=1;
|
---|
61 | md.nodeonbed=ones(md.numberofnodes,1);
|
---|
62 | md.nodeonsurface=ones(md.numberofnodes,1);
|
---|
63 | md.elementonbed=ones(md.numberofelements,1);
|
---|
64 | md.elementonsurface=ones(md.numberofelements,1);
|
---|
65 |
|
---|
66 | %Now, build the connectivity tables for this mesh.
|
---|
67 | md.nodeconnectivity=NodeConnectivity(md.elements,md.numberofnodes);
|
---|
68 | md.elementconnectivity=ElementConnectivity(md.elements,md.nodeconnectivity);
|
---|
69 |
|
---|
70 | %type of model
|
---|
71 | md.dim=2;
|
---|