1 | function 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
|
---|
14 | if ~is2d,
|
---|
15 | error('plot_contour error message: contour not supported for 3d meshes, project on a layer');
|
---|
16 | end
|
---|
17 |
|
---|
18 | %first, process data: must be on grids
|
---|
19 | if datatype==1,
|
---|
20 | %elements -> take average
|
---|
21 | data=averaging(md,data,0);
|
---|
22 | elseif datatype==2,
|
---|
23 | %nodes -> do nothing
|
---|
24 | elseif datatype==3,
|
---|
25 | %quiver -> take norm
|
---|
26 | data=sqrt(sum(datain.*datain,2));
|
---|
27 | else
|
---|
28 | error('datatype not supported yet');
|
---|
29 | end
|
---|
30 |
|
---|
31 | %prepare colors
|
---|
32 | if exist(options,'contouronly')
|
---|
33 | %remove the previous plots
|
---|
34 | cla
|
---|
35 | end
|
---|
36 | color=getfieldvalue(options,'contourcolor','y');
|
---|
37 |
|
---|
38 | %get contours levels
|
---|
39 | contourlevels=getfieldvalue(options,'contourlevels');
|
---|
40 | if isnumeric(contourlevels),
|
---|
41 | levels=round_ice(linspace(max(data),min(data),contourlevels),2);
|
---|
42 | else
|
---|
43 | levels=[];
|
---|
44 | for i=1:length(contourlevels)
|
---|
45 | levels(end+1)=contourlevels{i};
|
---|
46 | end
|
---|
47 | levels=sort(unique(levels),'descend');
|
---|
48 | end
|
---|
49 | numlevels=length(levels);
|
---|
50 |
|
---|
51 | %initialization of some variables
|
---|
52 | numberofelements=size(index,1);
|
---|
53 | elementslist=1:numberofelements;
|
---|
54 | c=[];
|
---|
55 | h=[];
|
---|
56 |
|
---|
57 | %get unique edges in mesh
|
---|
58 | %1: list of edges
|
---|
59 | edges=[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
|
---|
63 | vec=J;
|
---|
64 | %4: unique edges numbers in each triangle (2 triangles sharing the same edge will have
|
---|
65 | % the same edge number)
|
---|
66 | edges_tria=[vec(elementslist), vec(elementslist+numberofelements), vec(elementslist+2*numberofelements)];
|
---|
67 |
|
---|
68 | %segments [grids1 gruids2]
|
---|
69 | Seg1=index(:,[1 2]);
|
---|
70 | Seg2=index(:,[2 3]);
|
---|
71 | Seg3=index(:,[3 1]);
|
---|
72 |
|
---|
73 | %segment numbers [1;4;6;...]
|
---|
74 | Seg1_num=edges_tria(:,1);
|
---|
75 | Seg2_num=edges_tria(:,2);
|
---|
76 | Seg3_num=edges_tria(:,3);
|
---|
77 |
|
---|
78 | %value of data on each tips of the segments
|
---|
79 | Data1=data(Seg1);
|
---|
80 | Data2=data(Seg2);
|
---|
81 | Data3=data(Seg3);
|
---|
82 |
|
---|
83 | %get the ranges for each segment
|
---|
84 | Range1=sort(Data1,2);
|
---|
85 | Range2=sort(Data2,2);
|
---|
86 | Range3=sort(Data3,2);
|
---|
87 |
|
---|
88 | for 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
|
---|
234 | end
|
---|
235 |
|
---|
236 | %labels?
|
---|
237 | if (~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
|
---|
243 | end
|
---|