source: issm/trunk/src/m/model/contourenvelope.m@ 12706

Last change on this file since 12706 was 12706, checked in by Mathieu Morlighem, 13 years ago

merged trunk-jpl and trunk for revision 12703

File size: 4.1 KB
Line 
1function 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
13if nargin>2,
14 help contourenvelope
15 error('contourenvelope error message: bad usage');
16end
17if 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
33end
34
35%Now, build the connectivity tables for this mesh.
36%Computing connectivity
37if (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);
39end
40if (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);
42end
43
44%get nodes inside profile
45mesh.elementconnectivity=md.mesh.elementconnectivity;
46if 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;
52else
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;
58end
59
60if 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
84end
85
86%Find element on boundary
87%First: find elements on the boundary of the domain
88flag=mesh.elementconnectivity;
89if nargin==2,
90 flag(find(flag))=elemin(flag(find(flag)));
91end
92elementonboundary=double(prod(flag,2)==0 & sum(flag,2)>0);
93
94%Find segments on boundary
95pos=find(elementonboundary);
96num_segments=length(pos);
97segments=zeros(num_segments,3);
98count=1;
99
100for 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
139end
Note: See TracBrowser for help on using the repository browser.