1 | function segments=contourenvelope(md,varargin)
|
---|
2 | %CONTOURENVELOPE - build a set of segments enveloping a contour .exp
|
---|
3 | %
|
---|
4 | % Usage:
|
---|
5 | % segments=contourenvelope(md,varargin)
|
---|
6 | %
|
---|
7 | % Example:
|
---|
8 | % segments=contourenvelope(md,'Stream.exp');
|
---|
9 | % segments=contourenvelope(md,md.mask.elementonfloatingice)
|
---|
10 | % segments=contourenvelope(md);
|
---|
11 |
|
---|
12 | %some checks
|
---|
13 | if nargin>2,
|
---|
14 | help contourenvelope
|
---|
15 | error('contourenvelope error message: bad usage');
|
---|
16 | end
|
---|
17 | if nargin==2,
|
---|
18 | flags=varargin{1};
|
---|
19 |
|
---|
20 | if ischar(flags),
|
---|
21 | file=flags;
|
---|
22 | file=varargin{1};
|
---|
23 | if ~exist(file),
|
---|
24 | error(['contourenvelope error message: file ' file ' not found']);
|
---|
25 | end
|
---|
26 | isfile=1;
|
---|
27 | elseif isnumeric(flags),
|
---|
28 | %do nothing for now
|
---|
29 | isfile=0;
|
---|
30 | else
|
---|
31 | error('contourenvelope error message: second argument should a file or an elements flag');
|
---|
32 | end
|
---|
33 | end
|
---|
34 |
|
---|
35 | %Now, build the connectivity tables for this mesh.
|
---|
36 | %Computing connectivity
|
---|
37 | if (size(md.mesh.vertexconnectivity,1)~=md.mesh.numberofvertices & size(md.mesh.vertexconnectivity,1)~=md.mesh.numberofvertices2d),
|
---|
38 | md.mesh.vertexconnectivity=NodeConnectivity(md.mesh.elements,md.mesh.numberofvertices);
|
---|
39 | end
|
---|
40 | if (size(md.mesh.elementconnectivity,1)~=md.mesh.numberofelements & size(md.mesh.elementconnectivity,1)~=md.mesh.numberofelements2d),
|
---|
41 | md.mesh.elementconnectivity=ElementConnectivity(md.mesh.elements,md.mesh.vertexconnectivity);
|
---|
42 | end
|
---|
43 |
|
---|
44 | %get nodes inside profile
|
---|
45 | mesh.elementconnectivity=md.mesh.elementconnectivity;
|
---|
46 | if md.mesh.dimension==2;
|
---|
47 | mesh.elements=md.mesh.elements;
|
---|
48 | mesh.x=md.mesh.x;
|
---|
49 | mesh.y=md.mesh.y;
|
---|
50 | mesh.numberofvertices=md.mesh.numberofvertices;
|
---|
51 | mesh.numberofelements=md.mesh.numberofelements;
|
---|
52 | else
|
---|
53 | mesh.elements=md.mesh.elements2d;
|
---|
54 | mesh.x=md.mesh.x2d;
|
---|
55 | mesh.y=md.mesh.y2d;
|
---|
56 | mesh.numberofvertices=md.mesh.numberofvertices2d;
|
---|
57 | mesh.numberofelements=md.mesh.numberofelements2d;
|
---|
58 | end
|
---|
59 |
|
---|
60 | if nargin==2,
|
---|
61 |
|
---|
62 | if isfile,
|
---|
63 | %get flag list of elements and nodes inside the contour
|
---|
64 | nodein=ContourToMesh(mesh.elements,mesh.x,mesh.y,file,'node',1);
|
---|
65 | elemin=(sum(nodein(mesh.elements),2)==size(mesh.elements,2));
|
---|
66 | %modify element connectivity
|
---|
67 | elemout=find(~elemin);
|
---|
68 | mesh.elementconnectivity(elemout,:)=0;
|
---|
69 | mesh.elementconnectivity(find(ismember(mesh.elementconnectivity,elemout)))=0;
|
---|
70 | else
|
---|
71 | %get flag list of elements and nodes inside the contour
|
---|
72 | nodein=zeros(mesh.numberofvertices,1);
|
---|
73 | elemin=zeros(mesh.numberofelements,1);
|
---|
74 |
|
---|
75 | pos=find(flags);
|
---|
76 | elemin(pos)=1;
|
---|
77 | nodein(mesh.elements(pos,:))=1;
|
---|
78 |
|
---|
79 | %modify element connectivity
|
---|
80 | elemout=find(~elemin);
|
---|
81 | mesh.elementconnectivity(elemout,:)=0;
|
---|
82 | mesh.elementconnectivity(find(ismember(mesh.elementconnectivity,elemout)))=0;
|
---|
83 | end
|
---|
84 | end
|
---|
85 |
|
---|
86 | %Find element on boundary
|
---|
87 | %First: find elements on the boundary of the domain
|
---|
88 | flag=mesh.elementconnectivity;
|
---|
89 | if nargin==2,
|
---|
90 | flag(find(flag))=elemin(flag(find(flag)));
|
---|
91 | end
|
---|
92 | elementonboundary=double(prod(flag,2)==0 & sum(flag,2)>0);
|
---|
93 |
|
---|
94 | %Find segments on boundary
|
---|
95 | pos=find(elementonboundary);
|
---|
96 | num_segments=length(pos);
|
---|
97 | segments=zeros(num_segments,3);
|
---|
98 | count=1;
|
---|
99 |
|
---|
100 | for i=1:num_segments,
|
---|
101 | el1=pos(i);
|
---|
102 | els2=mesh.elementconnectivity(el1,find(mesh.elementconnectivity(el1,:)));
|
---|
103 | if length(els2)>1,
|
---|
104 | flag=intersect(mesh.elements(els2(1),:),mesh.elements(els2(2),:));
|
---|
105 | nods1=mesh.elements(el1,:);
|
---|
106 | nods1(find(nods1==flag))=[];
|
---|
107 | segments(count,:)=[nods1 el1];
|
---|
108 |
|
---|
109 | ord1=find(nods1(1)==mesh.elements(el1,:));
|
---|
110 | ord2=find(nods1(2)==mesh.elements(el1,:));
|
---|
111 |
|
---|
112 | %swap segment nodes if necessary
|
---|
113 | if ( (ord1==1 & ord2==2) | (ord1==2 & ord2==3) | (ord1==3 & ord2==1) ),
|
---|
114 | temp=segments(count,1);
|
---|
115 | segments(count,1)=segments(count,2);
|
---|
116 | segments(count,2)=temp;
|
---|
117 | end
|
---|
118 | segments(count,1:2)=fliplr(segments(count,1:2));
|
---|
119 | count=count+1;
|
---|
120 | else
|
---|
121 | nods1=mesh.elements(el1,:);
|
---|
122 | flag=setdiff(nods1,mesh.elements(els2,:));
|
---|
123 | for j=1:3,
|
---|
124 | nods=nods1; nods(j)=[];
|
---|
125 | if any(ismember(flag,nods)),
|
---|
126 | segments(count,:)=[nods el1];
|
---|
127 | ord1=find(nods(1)==mesh.elements(el1,:));
|
---|
128 | ord2=find(nods(2)==mesh.elements(el1,:));
|
---|
129 | if ( (ord1==1 & ord2==2) | (ord1==2 & ord2==3) | (ord1==3 & ord2==1) ),
|
---|
130 | temp=segments(count,1);
|
---|
131 | segments(count,1)=segments(count,2);
|
---|
132 | segments(count,2)=temp;
|
---|
133 | end
|
---|
134 | segments(count,1:2)=fliplr(segments(count,1:2));
|
---|
135 | count=count+1;
|
---|
136 | end
|
---|
137 | end
|
---|
138 | end
|
---|
139 | end
|
---|