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