source: issm/trunk/src/m/model/plot/plot_contour.m@ 5901

Last change on this file since 5901 was 4867, checked in by seroussi, 15 years ago

fixed problem in plot_contour

File size: 6.2 KB
Line 
1function plot_contour(md,datain,options);
2%PLOT_CONTOUR - plot contours of a given field
3%
4% Usage:
5% plot_contour(md,data,options);
6%
7% See also: PLOTMODEL
8
9%process data and model
10[x y z index is2d]=processmesh(md,[],options);
11[data datatype]=processdata(md,datain,options);
12
13%check is2d
14if ~is2d,
15 error('plot_contour error message: contour not supported for 3d meshes, project on a layer');
16end
17
18%first, process data: must be on grids
19if datatype==1,
20 %elements -> take average
21 data=averaging(md,data,0);
22elseif datatype==2,
23 %nodes -> do nothing
24elseif datatype==3,
25 %quiver -> take norm
26 data=sqrt(sum(datain.*datain,2));
27else
28 error('datatype not supported yet');
29end
30
31%prepare colors
32if exist(options,'contouronly')
33 %remove the previous plots
34 cla
35end
36color=getfieldvalue(options,'contourcolor','y');
37
38%get contours levels
39contourlevels=getfieldvalue(options,'contourlevels');
40if isnumeric(contourlevels),
41 levels=round_ice(linspace(max(data),min(data),contourlevels),2);
42else
43 levels=[];
44 for i=1:length(contourlevels)
45 levels(end+1)=contourlevels{i};
46 end
47 levels=sort(unique(levels),'descend');
48end
49numlevels=length(levels);
50
51%initialization of some variables
52numberofelements=size(index,1);
53elementslist=1:numberofelements;
54c=[];
55h=[];
56
57%get unique edges in mesh
58%1: list of edges
59edges=[index(:,[1,2]); index(:,[2,3]); index(:,[3,1])];
60%2: find unique edges
61[edges,I,J]=unique(sort(edges,2),'rows');
62%3: unique edge numbers
63vec=J;
64%4: unique edges numbers in each triangle (2 triangles sharing the same edge will have
65% the same edge number)
66edges_tria=[vec(elementslist), vec(elementslist+numberofelements), vec(elementslist+2*numberofelements)];
67
68%segments [grids1 gruids2]
69Seg1=index(:,[1 2]);
70Seg2=index(:,[2 3]);
71Seg3=index(:,[3 1]);
72
73%segment numbers [1;4;6;...]
74Seg1_num=edges_tria(:,1);
75Seg2_num=edges_tria(:,2);
76Seg3_num=edges_tria(:,3);
77
78%value of data on each tips of the segments
79Data1=data(Seg1);
80Data2=data(Seg2);
81Data3=data(Seg3);
82
83%get the ranges for each segment
84Range1=sort(Data1,2);
85Range2=sort(Data2,2);
86Range3=sort(Data3,2);
87
88for i=1:numlevels
89
90 level=levels(i);
91
92 %find the segments that contain this value
93 pos1=(Range1(:,1)<level & Range1(:,2)>level);
94 pos2=(Range2(:,1)<level & Range2(:,2)>level);
95 pos3=(Range3(:,1)<level & Range3(:,2)>level);
96
97 %get elements
98 poselem12=(pos1 & pos2);
99 poselem13=(pos1 & pos3);
100 poselem23=(pos2 & pos3);
101 poselem=find(poselem12 | poselem13 | poselem23);
102 numelems=length(poselem);
103
104 %if no element has been flagged, skip to the next level
105 if numelems==0,
106 continue,
107 end
108
109 %go through the elements and build the coordinates for each segment (1 by element)
110 x1=zeros(numelems,1);
111 x2=zeros(numelems,1);
112 y1=zeros(numelems,1);
113 y2=zeros(numelems,1);
114 edge_l=zeros(numelems,2);
115
116 for j=1:numelems,
117
118 weight1=(level-Data1(poselem(j),1))/(Data1(poselem(j),2)-Data1(poselem(j),1));
119 weight2=(level-Data2(poselem(j),1))/(Data2(poselem(j),2)-Data2(poselem(j),1));
120 weight3=(level-Data3(poselem(j),1))/(Data3(poselem(j),2)-Data3(poselem(j),1));
121
122 if poselem12(poselem(j));
123
124 x1(j)=x(Seg1(poselem(j),1))+weight1*(x(Seg1(poselem(j),2))-x(Seg1(poselem(j),1)));
125 x2(j)=x(Seg2(poselem(j),1))+weight2*(x(Seg2(poselem(j),2))-x(Seg2(poselem(j),1)));
126 y1(j)=y(Seg1(poselem(j),1))+weight1*(y(Seg1(poselem(j),2))-y(Seg1(poselem(j),1)));
127 y2(j)=y(Seg2(poselem(j),1))+weight2*(y(Seg2(poselem(j),2))-y(Seg2(poselem(j),1)));
128 edge_l(j,1)=Seg1_num(poselem(j));
129 edge_l(j,2)=Seg2_num(poselem(j));
130
131 elseif poselem13(poselem(j)),
132
133 x1(j)=x(Seg1(poselem(j),1))+weight1*(x(Seg1(poselem(j),2))-x(Seg1(poselem(j),1)));
134 x2(j)=x(Seg3(poselem(j),1))+weight3*(x(Seg3(poselem(j),2))-x(Seg3(poselem(j),1)));
135 y1(j)=y(Seg1(poselem(j),1))+weight1*(y(Seg1(poselem(j),2))-y(Seg1(poselem(j),1)));
136 y2(j)=y(Seg3(poselem(j),1))+weight3*(y(Seg3(poselem(j),2))-y(Seg3(poselem(j),1)));
137 edge_l(j,1)=Seg1_num(poselem(j));
138 edge_l(j,2)=Seg3_num(poselem(j));
139
140 elseif poselem23(poselem(j)),
141
142 x1(j)=x(Seg2(poselem(j),1))+weight2*(x(Seg2(poselem(j),2))-x(Seg2(poselem(j),1)));
143 x2(j)=x(Seg3(poselem(j),1))+weight3*(x(Seg3(poselem(j),2))-x(Seg3(poselem(j),1)));
144 y1(j)=y(Seg2(poselem(j),1))+weight2*(y(Seg2(poselem(j),2))-y(Seg2(poselem(j),1)));
145 y2(j)=y(Seg3(poselem(j),1))+weight3*(y(Seg3(poselem(j),2))-y(Seg3(poselem(j),1)));
146 edge_l(j,1)=Seg2_num(poselem(j));
147 edge_l(j,2)=Seg3_num(poselem(j));
148 else
149 %it shoud not go here
150 end
151 end
152
153 %now that we have the segments, we must try to connect them...
154
155 %loop over the subcontours
156 while ~isempty(edge_l),
157
158 %take the right edge of the second segment and connect it to the next segments if any
159 e1=edge_l(1,1); e2=edge_l(1,2);
160 xc=[x1(1);x2(1)]; yc=[y1(1);y2(1)];
161
162 %erase the lines corresponding to this edge
163 edge_l(1,:)=[];
164 x1(1)=[]; x2(1)=[];
165 y1(1)=[]; y2(1)=[];
166
167 [ro1,co1]=find(edge_l==e1);
168
169 while ~isempty(ro1)
170
171 if co1==1,
172 xc=[x2(ro1);xc]; yc=[y2(ro1);yc];
173
174 %next edge:
175 e1=edge_l(ro1,2);
176
177 else
178 xc=[x1(ro1);xc]; yc=[y1(ro1);yc];
179
180 %next edge:
181 e1=edge_l(ro1,1);
182 end
183
184 %erase the lines of this
185 edge_l(ro1,:)=[];
186 x1(ro1)=[]; x2(ro1)=[];
187 y1(ro1)=[]; y2(ro1)=[];
188
189 %next connection
190 [ro1,co1]=find(edge_l==e1);
191 end
192
193 %same thing the other way (to the right)
194 [ro2,co2]=find(edge_l==e2);
195
196 while ~isempty(ro2)
197
198 if co2==1,
199 xc=[xc;x2(ro2)]; yc=[yc;y2(ro2)];
200
201 %next edge:
202 e2=edge_l(ro2,2);
203 else
204 xc=[xc;x1(ro2)]; yc=[yc;y1(ro2)];
205
206 %next edge:
207 e2=edge_l(ro2,1);
208 end
209
210 %erase the lines of this
211 edge_l(ro2,:)=[];
212 x1(ro2)=[]; x2(ro2)=[];
213 y1(ro2)=[]; y2(ro2)=[];
214
215 %next connection
216 [ro2,co2]=find(edge_l==e2);
217 end
218
219 %we now have one subcontour ready to be plotted
220 zc=level*ones(length(xc)+1,1);
221 if getfieldvalue(options,'contouronly',0),
222 h=[h;patch('Xdata',[xc;NaN],'Ydata',[yc;NaN],'Zdata',zc,'Cdata',zc,'facecolor','none','edgecolor','flat')];
223 hold on
224 else
225 h=[h;patch('Xdata',[xc;NaN],'Ydata',[yc;NaN],'facecolor','none','edgecolor',color)];
226 hold on
227 end
228
229 % Update the CS data structure as per "contours.m"
230 % so that clabel works
231 c = horzcat(c,[level, xc'; length(xc), yc']);
232
233 end
234end
235
236%labels?
237if (~strcmpi(getfieldvalue(options,'contourticks','on'),'off') & ~isempty(c) & ~isempty(h))
238 if exist(options,'contouronly')
239 clabel(c,h);
240 else
241 clabel(c,h,'color',color,'FontSize',14);
242 end
243end
Note: See TracBrowser for help on using the repository browser.