source: issm/trunk-jpl/src/m/kml/kml_part_edges.m@ 13646

Last change on this file since 13646 was 13646, checked in by Mathieu Morlighem, 12 years ago

CHG: cosmetics, removed multiple blank lines
File size: 20.2 KB
Line 
1%%
2% create kml linestrings for the partition edges.
3%
4% [kfold]=kml_part_edges(md,params)
5%
6% where the required input is:
7% md (model, model class object)
8%
9% the optional input is:
10% params (string/numeric, parameter names and values)
11%
12% and the optional input is:
13% latsgn (numeric, +1/-1 for north/south latitude)
14% data (numeric, element or nodal results data)
15% alt (numeric, altitude for polygons, default 10000)
16% cmin (numeric, minimum of color map)
17% cmax (numeric, maximum of color map)
18% cmap (char or numeric, colormap definition)
19% prtplt (char, 'off'/'no' for partition segment plot)
20%
21% and the required output is:
22% kfold (kml_folder, folder of linestring placemarks)
23%
24function [kfold]=kml_part_edges(varargin)
25
26if ~nargin
27 help kml_part_edges
28 return
29end
30
31%% process input data
32
33iarg=1;
34if (nargin >= 1)
35 md=varargin{1};
36end
37if ~exist('md','var') || isempty(md) || ~isa(md,'model')
38 error(['Model ''' inputname(iarg) ''' is unrecognized class ''' class(md) '''.']);
39end
40
41% parameters
42
43iarg=iarg+1;
44while (iarg <= nargin-1)
45 if ischar(varargin{iarg})
46 if ~isempty(strmatch(varargin{iarg},...
47 {'latsgn','data','alt',...
48 'cmin','cmax','cmap','prtplt'},...
49 'exact'))
50 eval([varargin{iarg} '=varargin{iarg+1};']);
51 disp([varargin{iarg} '=' any2str(varargin{iarg+1},20) ';']);
52 else
53 warning([varargin{iarg} '=' any2str(varargin{iarg+1},20) ' is not recognized.']);
54 end
55 else
56 error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']);
57 end
58 if strcmpi(varargin{iarg},'data')
59 cdata=inputname(iarg+1);
60 end
61 iarg=iarg+2;
62end
63
64if isempty(md.mesh.lat) || ((numel(md.mesh.lat) == 1) && isnan(md.mesh.lat)) || ...
65 isempty(md.mesh.long) || ((numel(md.mesh.long) == 1) && isnan(md.mesh.long))
66 if ~exist('latsgn','var')
67 error(['Missing ''latsgn'' parameter to calculate missing lat/long data.']);
68 elseif (abs(latsgn) ~= 1)
69 error(['Incorrect latsgn=' num2str(latsgn) ' parameter to calculate missing lat/long data.']);
70 else
71 display('Converting x/y data to lat/long data.');
72 [md.mesh.lat,md.mesh.long]=xy2ll(md.x,md.y,latsgn);
73 end
74end
75
76if exist('data','var') && ~isempty(data)
77 if (numel(data)==md.mesh.numberofelements)
78 edata=data;
79 elseif (numel(data)==md.mesh.numberofvertices)
80 ndata=data;
81 display('Averaging nodal data to element data.');
82 edata=zeros(1,md.mesh.numberofelements);
83 for i=1:size(md.mesh.elements,1)
84 for j=1:size(md.mesh.elements,2)
85 edata(i)=edata(i)+ndata(md.mesh.elements(i,j));
86 end
87 edata(i)=edata(i)/size(md.mesh.elements,2);
88 end
89 else
90 error(['Data has incorrect number of ' num2str(numel(data)) ' values.']);
91 end
92end
93
94% colormap command operates on a figure, so create an invisible one
95% (could also directly call colormaps, e.g. jet(64), but risky)
96
97hfig=figure('Visible','off');
98if exist('cmap','var')
99 colormap(cmap)
100end
101cmap=colormap;
102close(hfig)
103
104if exist('edata','var')
105 if ~exist('cmin','var')
106 cmin=min(min(edata));
107 end
108 if ~exist('cmax','var')
109 cmax=max(max(edata));
110 end
111end
112
113if ~exist('alt','var')
114 alt=10000;
115end
116
117%% write folder for partition edges
118
119if (~exist('prtplt','var') || strncmpi(prtplt,'on' ,2) || strncmpi(prtplt,'y',1)) && ...
120 md.qmu.numberofpartitions
121 kfold=kml_folder();
122 kfold.name ='Partition Edges';
123 kfold.visibility=1;
124 kfold.descript =sprintf('Partitions=%d, Nodes=%d',...
125 md.qmu.numberofpartitions,md.mesh.numberofvertices);
126 kfold.feature ={repmat(kml_placemark(),1,md.qmu.numberofpartitions)};
127
128% write each partition as a linestring multigeometry placemark
129
130 disp(['Writing ' num2str(md.qmu.numberofpartitions) ' partitions as KML linestrings.']);
131 epart=md.qmu.partition(md.mesh.elements)+1;
132 if exist('ndata','var') || exist('edata','var')
133 pdata=zeros(1,md.qmu.numberofpartitions);
134 pdata(:)=NaN;
135 end
136
137% loop over each partition
138
139 for k=1:md.qmu.numberofpartitions
140% disp(['partition k=' int2str(k)])
141
142% for each partition, find all the included elements and determine the
143% perimeter (including those shared by another partition)
144
145 [icol,irow]=find(epart'==k);
146 if isempty(irow)
147 continue;
148 end
149 irow=unique(irow);
150 elemp=md.mesh.elements(irow,:);
151 epartp=epart(irow,:);
152 nodeconp=kmlnodeconnectivity(elemp,md.mesh.numberofvertices);
153 [edgeadjp]=edgeadjacency(elemp,nodeconp);
154 [edgeper,elemper,iloop]=edgeperimeter(elemp,nodeconp,edgeadjp);
155 iloop(end+1)=size(edgeper,1)+1;
156
157% determine the data to be used for the colors (if any)
158
159 if exist('ndata','var')
160 pdata(k)=ndata(find(md.qmu.partition+1==k,1));
161 elseif exist('edata','var')
162 for i=1:size(epartp,1)
163 if isempty(find(epart(i,:)~=k,1))
164 pdata(k)=edata(irow(i));
165 break
166 end
167 end
168 if isnan(pdata(k))
169 warning('Data for Partition %d is not defined.\n',k)
170 end
171 end
172
173% set up the placemark with multigeometry
174
175 kplace=kml_placemark();
176 if (length(iloop)-1 > 1)
177 kplace.name =sprintf('Partition %d (%d loops)',k,length(iloop)-1);
178 else
179 kplace.name =sprintf('Partition %d',k);
180 end
181 kplace.visibility=1;
182 if exist('pdata','var')
183 kplace.descript =sprintf('Partition data: %g',pdata(k));
184 imap = fix((pdata(k)-cmin)/(cmax-cmin)*size(cmap,1))+1;
185 if (imap >= 1) && (imap <= size(cmap,1))
186 kplace.styleurl =sprintf('#MatlabColor%d',imap);
187 elseif (pdata(k) == cmax)
188 kplace.styleurl =sprintf('#MatlabColor%d',size(cmap,1));
189 else
190 kplace.styleurl =sprintf('#BlackLineEmptyPoly');
191 end
192 else
193 kplace.styleurl =sprintf('#BlackLineRandomPoly');
194 end
195
196 kmgeom=kml_multigeometry();
197 kmgeom.geometry ={repmat(kml_linestring(),1,length(iloop)-1)};
198
199% loop over each loop of the perimeter for the given partition
200
201 for i=1:length(iloop)-1
202 kline=kml_linestring();
203 kline.extrude =1;
204 kline.tessellate=1;
205 kline.altmode ='relativeToGround';
206 kline.coords =zeros(0,3);
207
208 elast=0;
209 nlast=0;
210 slast=0;
211 lat=[];
212 long=[];
213
214% loop over the element edges on the loop of the partition
215
216 j=iloop(i);
217 while (j < iloop(i+1))
218% find which side of element is referenced in perimeter list
219 for l=1:size(elemp,2)
220 if ((elemp(elemper(j),l) == edgeper(j,1)) && ...
221 (elemp(elemper(j),mod(l,3)+1) == edgeper(j,2))) || ...
222 ((elemp(elemper(j),l) == edgeper(j,2)) && ...
223 (elemp(elemper(j),mod(l,3)+1) == edgeper(j,1)))
224 jedge=l;
225 break
226 end
227 end
228
229% check if element side connects nodes in partition
230 if (epartp(elemper(j),jedge) == k) && ...
231 (epartp(elemper(j),mod(jedge,3)+1) == k)
232% write out specified element side
233% disp(['segment j=' int2str(j) ' unshared edge ' int2str(edgeper(j,1)) ' to ' int2str(edgeper(j,2)) ' on side ' int2str(jedge) ' from element ' int2str(elemper(j)) ' written.'])
234% if first edge, write out first node
235 if ~elast
236 kline.coords(end+1,:)=[md.mesh.long(edgeper(j,1)) md.mesh.lat(edgeper(j,1)) alt];
237 end
238 kline.coords(end+1,:)=[md.mesh.long(edgeper(j,2)) md.mesh.lat(edgeper(j,2)) alt];
239 elast=elemper(j);
240 nlast=edgeper(j,2);
241 slast=0;
242 j=j+1;
243
244% element not entirely within partition, so figure out boundary
245 else
246% disp(['segment j=' int2str(j) ' from element ' int2str(elemper(j)) ' shared by other partitions.'])
247 ielem=elemper(j);
248
249% follow partition boundary through elements not wholly in partition
250% (may include elements not in perimeter list)
251
252 while 1
253% if first edge, figure out direction from perimeter edge direction
254 if ~nlast && ~slast
255 nlast=find(elemp(ielem,:)==edgeper(j,1));
256 nnext=find(elemp(ielem,:)==edgeper(j,2));
257 if (nlast+nnext == 3)
258 slast=1;
259 elseif (nlast+nnext == 5)
260 slast=2;
261 elseif (nlast+nnext == 4)
262 slast=3;
263 end
264 if (nnext+(6-nlast-nnext) == 3)
265 snext=1;
266 elseif (nnext+(6-nlast-nnext) == 5)
267 snext=2;
268 elseif (nnext+(6-nlast-nnext) == 4)
269 snext=3;
270 end
271
272% find how many nodes of current element are in current partition
273% (1 or 2, not 3) and handle each permutation separately
274 ipart=find(epartp(ielem,:)==k);
275% two nodes are in current partition, so cut off other node
276 if (length(ipart) == 2)
277 switch 6-sum(ipart)
278 case nlast
279 slast=6-slast-snext;
280 nlast=0;
281 case nnext
282 if (epartp(ielem,nnext) == k)
283 nlast=nnext;
284 end
285 otherwise
286 slast=6-slast-snext;
287 nlast=0;
288 end
289% one node is in current partition
290 else
291% all different, so cut through centroid
292 if (epartp(ielem,1) ~= epartp(ielem,2)) && ...
293 (epartp(ielem,2) ~= epartp(ielem,3)) && ...
294 (epartp(ielem,3) ~= epartp(ielem,1))
295 switch ipart
296 case {nlast,nnext}
297 if (epartp(ielem,nnext) == k)
298 nlast=nnext;
299 end
300 otherwise
301 slast=6-slast-snext;
302 nlast=0;
303 end
304% other two are in the same partition, so cut them off
305 else
306 switch ipart
307 case nlast
308 if (epartp(ielem,nnext) == k)
309 nlast=nnext;
310 end
311 case nnext
312 slast=snext;
313 nlast=0;
314 otherwise
315 slast=6-slast-snext;
316 nlast=0;
317 end
318 end
319 end
320
321% last edge exited last element at node
322 if nlast
323% write out first node of first side for half-edge to midpoint
324% disp(['segment j=' int2str(j) ' unshared half edge from node ' int2str(elemp(ielem,nlast)) ' (node ' int2str(nlast) ') on side ' int2str(slast) ' from element ' int2str(ielem) ' written.'])
325 kline.coords(end+1,:)=[md.mesh.long(elemp(ielem,nlast)) ...
326 md.mesh.lat(elemp(ielem,nlast)) alt];
327 end
328 nlast=0;
329
330% write out midpoint of first side
331 kline.coords(end+1,:)=[(md.mesh.long(elemp(ielem,slast))...
332 +md.mesh.long(elemp(ielem,mod(slast,3)+1)))/2. ...
333 (md.mesh.lat(elemp(ielem,slast))...
334 +md.mesh.lat(elemp(ielem,mod(slast,3)+1)))/2. alt];
335 end
336
337% last edge exited last element at node
338 if nlast
339 if elast
340% find where last node on previous element occurs on current element
341 nlast=find(elemp(ielem,:)==nlast,1);
342 end
343% half-edge occurs on unshared side from current node (unique unless mesh
344% is only attached at node)
345 switch nlast
346 case 1
347 if ~edgeadjp(ielem,1)
348 nnext=2;
349 slast=1;
350 else
351 nnext=3;
352 slast=3;
353 end
354 case 2
355 if ~edgeadjp(ielem,2)
356 nnext=3;
357 slast=2;
358 else
359 nnext=1;
360 slast=1;
361 end
362 case 3
363 if ~edgeadjp(ielem,3)
364 nnext=1;
365 slast=3;
366 else
367 nnext=2;
368 slast=2;
369 end
370 end
371% write out half-edge from current node to midpoint of unshared side
372% disp(['segment j=' int2str(j) ' unshared half edge from node ' int2str(elemp(ielem,nlast)) ' (node ' int2str(nlast) ') on side ' int2str(slast) ' from element ' int2str(ielem) ' written.'])
373 kline.coords(end+1,:)=[(md.mesh.long(elemp(ielem,nlast))...
374 +md.mesh.long(elemp(ielem,nnext)))/2. ...
375 (md.mesh.lat(elemp(ielem,nlast))...
376 +md.mesh.lat(elemp(ielem,nnext)))/2. alt];
377 nlast=0;
378
379% last edge exited last element at midpoint of side
380 elseif slast
381 if elast
382% find where last side on previous element occurs on current element
383 slast=find(edgeadjp(ielem,:)==elast,1);
384 end
385 end
386
387% find how many nodes of current element are in current partition
388% (1 or 2, not 3) and handle each permutation separately
389 ipart=find(epartp(ielem,:)==k);
390 if (length(ipart) == 2)
391% two nodes are in current partition, so cut off other node
392 switch 6-sum(ipart)
393 case 1
394 snext=3+1-slast;
395 case 2
396 snext=1+2-slast;
397 case 3
398 snext=2+3-slast;
399 end
400 else
401 if (epartp(ielem,1) ~= epartp(ielem,2)) && ...
402 (epartp(ielem,2) ~= epartp(ielem,3)) && ...
403 (epartp(ielem,3) ~= epartp(ielem,1))
404% all different, so cut through centroid
405% disp(['element ielem=' int2str(ielem) ' centroid written.'])
406 kline.coords(end+1,:)=[sum(md.mesh.long(elemp(ielem,:)))/3. ...
407 sum(md.mesh.lat(elemp(ielem,:)))/3. alt];
408 end
409% one node is in current partition, so cut off other two nodes
410 switch ipart
411 case 1
412 snext=3+1-slast;
413 case 2
414 snext=1+2-slast;
415 case 3
416 snext=2+3-slast;
417 end
418 end
419% write out midpoint of opposite side
420% disp(['segment j=' int2str(j) ' internal edge from side ' int2str(slast) ' to side ' int2str(snext) ' from element ' int2str(ielem) ' written.'])
421 kline.coords(end+1,:)=[(md.mesh.long(elemp(ielem,snext))...
422 +md.mesh.long(elemp(ielem,mod(snext,3)+1)))/2. ...
423 (md.mesh.lat(elemp(ielem,snext))...
424 +md.mesh.lat(elemp(ielem,mod(snext,3)+1)))/2. alt];
425 elast=ielem;
426 nlast=0;
427 slast=snext;
428% find adjacent element to opposite side
429 ielem=edgeadjp(elast,slast);
430% if opposite side is unshared, find it in edge perimeter list
431 if ~ielem
432 jlast=find(elemper(j:end)==elast)+j-1;
433 j=0;
434 for l=1:length(jlast)
435 if ((elemp(elast,slast) == edgeper(jlast(l),1)) && ...
436 (elemp(elast,mod(slast,3)+1) == edgeper(jlast(l),2))) || ...
437 ((elemp(elast,slast) == edgeper(jlast(l),2)) && ...
438 (elemp(elast,mod(slast,3)+1) == edgeper(jlast(l),1)))
439 j=jlast(l);
440 break
441 end
442 end
443 if ~j
444 j=iloop(i+1)-1;
445 end
446% write out half-edge from midpoint of unshared side to node
447 if (epartp(elast,slast) == k)
448 nnext=slast;
449 else
450 nnext=mod(slast,3)+1;
451 end
452% disp(['segment j=' int2str(j) ' unshared half edge on side ' int2str(slast) ' to node ' int2str(elemp(elast,nnext)) ' (node ' int2str(nnext) ') from element ' int2str(elast) ' written.'])
453 kline.coords(end+1,:)=[md.mesh.long(elemp(elast,nnext)) ...
454 md.mesh.lat(elemp(elast,nnext)) alt];
455 break
456% if not unshared, advance perimeter list and watch for end
457 else
458 if (elast == elemper(j))
459 if (j+1 < iloop(i+1)) && ...
460 ~isempty(find(elemper(j+1:end)~=elast,1))
461 j=j+find(elemper(j+1:end)~=elast,1);
462 else
463 break
464 end
465 end
466 end
467 end
468 j=j+1;
469 end
470 end
471
472 kmgeom.geometry{1}(i)=kline;
473 clear kline
474 end
475
476 kplace.geometry=kmgeom;
477 kfold.feature{1}(k)=kplace;
478 clear kmgeom kplace
479 end
480end
481
482end
Note: See TracBrowser for help on using the repository browser.