Changeset 6417


Ignore:
Timestamp:
10/25/10 14:25:10 (14 years ago)
Author:
jschierm
Message:

Separated kml_mesh_write.m into multiple functions in different files.

Location:
issm/trunk/src/m/kml
Files:
6 added
1 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/m/kml/kml_mesh_write.m

    r6415 r6417  
    1 %
     1%%
    22%  write a kml file of the mesh from the model.
    33%
    … …  
    248248fprintf(fid,'</kml>\n');
    249249
    250 fclose(fid);
    251 display('End of file successfully written.');
    252 
    253 end
    254 
    255 %%
    256 %  create kml polygons for the element mesh.
    257 %
    258 %  []=kml_mesh_elem(fid,md,params)
    259 %
    260 %  where the required input is:
    261 %    fid           (numeric, file ID of .kml file)
    262 %    md            (model, model class object)
    263 %
    264 %  the optional input is:
    265 %    params        (string/numeric, parameter names and values)
    266 %
    267 %  and the optional input is:
    268 %    data          (numeric, element or nodal results data)
    269 %    alt           (numeric, altitude for polygons, default 10000)
    270 %    cmin          (numeric, minimum of color map)
    271 %    cmax          (numeric, maximum of color map)
    272 %    cmap          (char or numeric, colormap definition)
    273 %
    274 function []=kml_mesh_elem(varargin)
    275 
    276 if ~nargin
    277     help kml_mesh_elem
    278     return
    279 end
    280 
    281 %%  process input data
    282 
    283 iarg=1;
    284 if (nargin >= 1)
    285     fid=varargin{1};
    286 end
    287 if ~exist('fid','var') || isempty(fid) || (fid < 0)
    288     error('File ID ''%d'' must be open.',fid);
    289 end
    290 
    291 iarg=iarg+1;
    292 if (nargin >= 2)
    293     md=varargin{2};
    294 end
    295 if ~exist('md','var') || isempty(md) || ~isa(md,'model')
    296     error(['Model ''' inputname(iarg) ''' is unrecognized class ''' class(md) '''.']);
    297 end
    298 
    299 %  parameters
    300 
    301 iarg=iarg+1;
    302 while (iarg <= nargin-1)
    303     if ischar(varargin{iarg})
    304         eval([varargin{iarg} '=varargin{iarg+1};']);
    305         if (numel(varargin{iarg+1}) <= 20)
    306             disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']);
    307         else
    308             disp([varargin{iarg} '=' string_size(varargin{iarg+1}) ' ' class(varargin{iarg+1}) ';']);
    309         end
    310         if strcmpi(varargin{iarg},'data')
    311             cdata=inputname(iarg+1);
    312         end
    313     else
    314         error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']);
    315     end
    316     iarg=iarg+2;
    317 end
    318 
    319 if exist('data','var') && ~isempty(data)
    320     if     (numel(data)==md.numberofelements)
    321         edata=data;
    322     elseif (numel(data)==md.numberofgrids)
    323         ndata=data;
    324         display('Averaging nodal data to element data.');
    325         edata=zeros(1,md.numberofelements);
    326         for i=1:size(md.elements,1)
    327             for j=1:size(md.elements,2)
    328                 edata(i)=edata(i)+ndata(md.elements(i,j));
    329             end
    330             edata(i)=edata(i)/size(md.elements,2);
    331         end
    332     else
    333         error(['Data has incorrect number of ' num2str(numel(data)) ' values.']);
    334     end
    335 end
    336 
    337 %  colormap command operates on a figure, so create an invisible one
    338 %  (could also directly call colormaps, e.g. jet(64), but risky)
    339 
    340 hfig=figure('Visible','off');
    341 if exist('cmap','var')
    342     colormap(cmap)
    343 end
    344 cmap=colormap;
    345 close(hfig)
    346    
    347 if exist('edata','var')
    348     if ~exist('cmin','var')
    349         cmin=min(min(edata));
    350     end
    351     if ~exist('cmax','var')
    352         cmax=max(max(edata));
    353     end
    354 end
    355 
    356 if ~exist('alt','var')
    357     alt=10000;
    358 end
    359 
    360 %%  write folder for mesh
    361 
    362 fprintf(fid,'    <Folder>\n');
    363 if exist('cdata','var') && ~isempty(cdata)
    364     fprintf(fid,'      <name>Data: %s</name>\n',cdata);
     250if (fclose(fid) < 0)
     251    error('File ''%s'' could not be closed.',filek2);
    365252else
    366     fprintf(fid,'      <name>Mesh</name>\n');
    367 end
    368 fprintf(fid,'      <visibility>1</visibility>\n');
    369 fprintf(fid,'      <description>Elements=%d, Grids=%d</description>\n',...
    370     md.numberofelements,md.numberofgrids);
    371 
    372 %  write each element as a polygon
    373 
    374 disp(['Writing ' num2str(size(md.elements,1)) ' tria elements as KML polygons.']);
    375 for i=1:size(md.elements,1)
    376     fprintf(fid,'      <Placemark>\n');
    377     fprintf(fid,'        <name>Element %d</name>\n',i);
    378     fprintf(fid,'        <visibility>1</visibility>\n');
    379     if exist('edata','var')
    380         fprintf(fid,'        <description>Element data: %g</description>\n',edata(i));
    381         imap = fix((edata(i)-cmin)/(cmax-cmin)*size(cmap,1))+1;
    382         if     (imap >= 1) && (imap <= size(cmap,1))
    383             fprintf(fid,'        <styleUrl>#MatlabColor%d</styleUrl>\n',imap);
    384         elseif (edata(i) == cmax)
    385             fprintf(fid,'        <styleUrl>#MatlabColor%d</styleUrl>\n',size(cmap,1));
    386         else
    387             fprintf(fid,'        <styleUrl>#BlackLineEmptyPoly</styleUrl>\n');
    388         end
    389     else
    390         fprintf(fid,'        <styleUrl>#BlackLineRandomPoly</styleUrl>\n');
    391     end
    392     fprintf(fid,'        <Polygon>\n');
    393     fprintf(fid,'          <extrude>1</extrude>\n');
    394     fprintf(fid,'          <altitudeMode>relativeToGround</altitudeMode>\n');
    395     fprintf(fid,'          <outerBoundaryIs>\n');
    396     fprintf(fid,'            <LinearRing>\n');
    397     fprintf(fid,'              <coordinates>\n');
    398     for j=1:size(md.elements,2)
    399         [lat(j),long(j)]=mapxy(md.x(md.elements(i,j)),md.y(md.elements(i,j)),'s');
    400         fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(j),lat(j),alt);
    401     end
    402     fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(1),lat(1),alt);
    403 
    404     fprintf(fid,'              </coordinates>\n');
    405     fprintf(fid,'            </LinearRing>\n');
    406     fprintf(fid,'          </outerBoundaryIs>\n');
    407     fprintf(fid,'        </Polygon>\n');
    408     fprintf(fid,'      </Placemark>\n');
    409 end
    410 fprintf(fid,'    </Folder>\n');
    411 
    412 end
    413 
    414 %%
    415 %  create kml linestrings for the flagged partition edges.
    416 %
    417 %  []=kml_part_flagedges(fid,md,params)
    418 %
    419 %  where the required input is:
    420 %    fid           (char, file ID of .kml file)
    421 %    md            (model, model class object)
    422 %
    423 %  the optional input is:
    424 %    params        (string/numeric, parameter names and values)
    425 %
    426 %  and the optional input is:
    427 %    alt           (numeric, altitude for polygons, default 10000)
    428 %    prtplt        (char, 'off'/'no' for partition segment plot)
    429 %
    430 function []=kml_part_flagedges(varargin)
    431 
    432 if ~nargin
    433     help kml_part_flagedges
    434     return
    435 end
    436 
    437 %%  process input data
    438 
    439 iarg=1;
    440 if (nargin >= 1)
    441     fid=varargin{1};
    442 end
    443 if ~exist('fid','var') || isempty(fid) || (fid < 0)
    444     error('File ID ''%d'' must be open.',fid);
    445 end
    446 
    447 iarg=iarg+1;
    448 if (nargin >= 2)
    449     md=varargin{2};
    450 end
    451 if ~exist('md','var') || isempty(md) || ~isa(md,'model')
    452     error(['Model ''' inputname(iarg) ''' is unrecognized class ''' class(md) '''.']);
    453 end
    454 
    455 %  parameters
    456 
    457 iarg=iarg+1;
    458 while (iarg <= nargin-1)
    459     if ischar(varargin{iarg})
    460         eval([varargin{iarg} '=varargin{iarg+1};']);
    461         if (numel(varargin{iarg+1}) <= 20)
    462             disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']);
    463         else
    464             disp([varargin{iarg} '=' string_size(varargin{iarg+1}) ' ' class(varargin{iarg+1}) ';']);
    465         end
    466     else
    467         error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']);
    468     end
    469     iarg=iarg+2;
    470 end
    471 
    472 if ~exist('alt','var')
    473     alt=10000;
    474 end
    475 
    476 %%  write folder for partition segments
    477 
    478 if (~exist('prtplt','var') || strncmpi(prtplt,'on' ,2) || strncmpi(prtplt,'y',1)) && ...
    479     md.npart
    480     [xseg,yseg]=flagedges(md.elements,md.x,md.y,md.part);
    481     fprintf(fid,'    <Folder>\n');
    482     fprintf(fid,'      <name>Partition Segments</name>\n');
    483     fprintf(fid,'      <visibility>1</visibility>\n');
    484     fprintf(fid,'      <description>Partitions=%d, Segments=%d</description>\n',...
    485         md.npart,size(xseg,1));
    486 
    487 %  write each segment as a linestring
    488 
    489     disp(['Writing ' num2str(size(xseg,1)) ' partition segments as KML linestrings.']);
    490     for i=1:size(xseg,1)
    491         fprintf(fid,'      <Placemark>\n');
    492         fprintf(fid,'        <name>Segment %d</name>\n',i);
    493         fprintf(fid,'        <visibility>1</visibility>\n');
    494         fprintf(fid,'        <styleUrl>#RedLineRedPoly</styleUrl>\n');
    495         fprintf(fid,'        <LineString>\n');
    496         fprintf(fid,'          <extrude>1</extrude>\n');
    497         fprintf(fid,'          <tessellate>1</tessellate>\n');
    498         fprintf(fid,'          <altitudeMode>relativeToGround</altitudeMode>\n');
    499         fprintf(fid,'          <coordinates>\n');
    500         for j=1:2
    501             [lat(j),long(j)]=mapxy(xseg(i,j),yseg(i,j),'s');
    502             fprintf(fid,'            %0.16g,%0.16g,%0.16g\n',long(j),lat(j),alt);
    503         end
    504 
    505         fprintf(fid,'          </coordinates>\n');
    506         fprintf(fid,'        </LineString>\n');
    507         fprintf(fid,'      </Placemark>\n');
    508     end
    509     fprintf(fid,'    </Folder>\n');
    510 end
    511 
    512 end
    513 
    514 %%
    515 %  create kml linestrings for the unshared element edges.
    516 %
    517 %  []=kml_unsh_edges(fid,md,params)
    518 %
    519 %  where the required input is:
    520 %    fid           (char, file ID of .kml file)
    521 %    md            (model, model class object)
    522 %
    523 %  the optional input is:
    524 %    params        (string/numeric, parameter names and values)
    525 %
    526 %  and the optional input is:
    527 %    alt           (numeric, altitude for polygons, default 10000)
    528 %    prtplt        (char, 'off'/'no' for partition segment plot)
    529 %
    530 function []=kml_unsh_edges(varargin)
    531 
    532 if ~nargin
    533     help kml_unsh_edges
    534     return
    535 end
    536 
    537 %%  process input data
    538 
    539 iarg=1;
    540 if (nargin >= 1)
    541     fid=varargin{1};
    542 end
    543 if ~exist('fid','var') || isempty(fid) || (fid < 0)
    544     error('File ID ''%d'' must be open.',fid);
    545 end
    546 
    547 iarg=iarg+1;
    548 if (nargin >= 2)
    549     md=varargin{2};
    550 end
    551 if ~exist('md','var') || isempty(md) || ~isa(md,'model')
    552     error(['Model ''' inputname(iarg) ''' is unrecognized class ''' class(md) '''.']);
    553 end
    554 
    555 %  parameters
    556 
    557 iarg=iarg+1;
    558 while (iarg <= nargin-1)
    559     if ischar(varargin{iarg})
    560         eval([varargin{iarg} '=varargin{iarg+1};']);
    561         if (numel(varargin{iarg+1}) <= 20)
    562             disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']);
    563         else
    564             disp([varargin{iarg} '=' string_size(varargin{iarg+1}) ' ' class(varargin{iarg+1}) ';']);
    565         end
    566     else
    567         error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']);
    568     end
    569     iarg=iarg+2;
    570 end
    571 
    572 if ~exist('alt','var')
    573     alt=10000;
    574 end
    575 
    576 %%  write folder for unshared edges
    577 
    578 if (~exist('prtplt','var') || strncmpi(prtplt,'on' ,2) || strncmpi(prtplt,'y',1)) && ...
    579     md.npart
    580     [edgeadj]=edgeadjacency(md.elements,md.nodeconnectivity);
    581     [icol,irow]=find(edgeadj'==0);
    582     edgeuns=zeros(length(irow),2);
    583     for i=1:length(irow)
    584         edgeuns(i,1)=md.elements(irow(i),icol(i));
    585         edgeuns(i,2)=md.elements(irow(i),mod(icol(i),size(md.elements,2))+1);
    586     end
    587     fprintf(fid,'    <Folder>\n');
    588     fprintf(fid,'      <name>Unshared Edges</name>\n');
    589     fprintf(fid,'      <visibility>1</visibility>\n');
    590     fprintf(fid,'      <description>Partitions=%d, Edges=%d</description>\n',...
    591         md.npart,size(edgeuns,1));
    592 
    593 %  write each edge as a linestring
    594 
    595     disp(['Writing ' num2str(size(edgeuns,1)) ' unshared edges as KML linestrings.']);
    596     for i=1:size(edgeuns,1)
    597         fprintf(fid,'      <Placemark>\n');
    598         fprintf(fid,'        <name>Edge %d</name>\n',i);
    599         fprintf(fid,'        <visibility>1</visibility>\n');
    600         fprintf(fid,'        <styleUrl>#RedLineRedPoly</styleUrl>\n');
    601         fprintf(fid,'        <LineString>\n');
    602         fprintf(fid,'          <extrude>1</extrude>\n');
    603         fprintf(fid,'          <tessellate>1</tessellate>\n');
    604         fprintf(fid,'          <altitudeMode>relativeToGround</altitudeMode>\n');
    605         fprintf(fid,'          <coordinates>\n');
    606         for j=1:2
    607             [lat(j),long(j)]=mapxy(md.x(edgeuns(i,j)),md.y(edgeuns(i,j)),'s');
    608             fprintf(fid,'            %0.16g,%0.16g,%0.16g\n',long(j),lat(j),alt);
    609         end
    610 
    611         fprintf(fid,'          </coordinates>\n');
    612         fprintf(fid,'        </LineString>\n');
    613         fprintf(fid,'      </Placemark>\n');
    614     end
    615     fprintf(fid,'    </Folder>\n');
    616 end
    617 
    618 end
    619 
    620 %%
    621 %  create kml polygons for the partition elements.
    622 %
    623 %  []=kml_part_elems(fid,md,params)
    624 %
    625 %  where the required input is:
    626 %    fid           (char, file ID of .kml file)
    627 %    md            (model, model class object)
    628 %
    629 %  the optional input is:
    630 %    params        (string/numeric, parameter names and values)
    631 %
    632 %  and the optional input is:
    633 %    data          (numeric, element or nodal results data)
    634 %    alt           (numeric, altitude for polygons, default 10000)
    635 %    cmin          (numeric, minimum of color map)
    636 %    cmax          (numeric, maximum of color map)
    637 %    cmap          (char or numeric, colormap definition)
    638 %    prtplt        (char, 'off'/'no' for partition segment plot)
    639 %
    640 function []=kml_part_elems(varargin)
    641 
    642 if ~nargin
    643     help kml_part_elems
    644     return
    645 end
    646 
    647 %%  process input data
    648 
    649 iarg=1;
    650 if (nargin >= 1)
    651     fid=varargin{1};
    652 end
    653 if ~exist('fid','var') || isempty(fid) || (fid < 0)
    654     error('File ID ''%d'' must be open.',fid);
    655 end
    656 
    657 iarg=iarg+1;
    658 if (nargin >= 2)
    659     md=varargin{2};
    660 end
    661 if ~exist('md','var') || isempty(md) || ~isa(md,'model')
    662     error(['Model ''' inputname(iarg) ''' is unrecognized class ''' class(md) '''.']);
    663 end
    664 
    665 %  parameters
    666 
    667 iarg=iarg+1;
    668 while (iarg <= nargin-1)
    669     if ischar(varargin{iarg})
    670         eval([varargin{iarg} '=varargin{iarg+1};']);
    671         if (numel(varargin{iarg+1}) <= 20)
    672             disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']);
    673         else
    674             disp([varargin{iarg} '=' string_size(varargin{iarg+1}) ' ' class(varargin{iarg+1}) ';']);
    675         end
    676         if strcmpi(varargin{iarg},'data')
    677             cdata=inputname(iarg+1);
    678         end
    679     else
    680         error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']);
    681     end
    682     iarg=iarg+2;
    683 end
    684 
    685 if exist('data','var') && ~isempty(data)
    686     if     (numel(data)==md.numberofelements)
    687         edata=data;
    688     elseif (numel(data)==md.numberofgrids)
    689         ndata=data;
    690         display('Averaging nodal data to element data.');
    691         edata=zeros(1,md.numberofelements);
    692         for i=1:size(md.elements,1)
    693             for j=1:size(md.elements,2)
    694                 edata(i)=edata(i)+ndata(md.elements(i,j));
    695             end
    696             edata(i)=edata(i)/size(md.elements,2);
    697         end
    698     else
    699         error(['Data has incorrect number of ' num2str(numel(data)) ' values.']);
    700     end
    701 end
    702 
    703 %  colormap command operates on a figure, so create an invisible one
    704 %  (could also directly call colormaps, e.g. jet(64), but risky)
    705 
    706 hfig=figure('Visible','off');
    707 if exist('cmap','var')
    708     colormap(cmap)
    709 end
    710 cmap=colormap;
    711 close(hfig)
    712    
    713 if exist('edata','var')
    714     if ~exist('cmin','var')
    715         cmin=min(min(edata));
    716     end
    717     if ~exist('cmax','var')
    718         cmax=max(max(edata));
    719     end
    720 end
    721 
    722 if ~exist('alt','var')
    723     alt=10000;
    724 end
    725 
    726 %  write folder for partition elements
    727 
    728 if (~exist('prtplt','var') || strncmpi(prtplt,'on' ,2) || strncmpi(prtplt,'y',1)) && ...
    729     md.npart
    730     fprintf(fid,'    <Folder>\n');
    731     fprintf(fid,'      <name>Partition Elements</name>\n');
    732     fprintf(fid,'      <visibility>1</visibility>\n');
    733     fprintf(fid,'      <description>Partitions=%d, Nodes=%d</description>\n',...
    734         md.npart,md.numberofgrids);
    735 
    736 %  write each partition loop as a polygon
    737 
    738     disp(['Writing ' num2str(md.npart) ' partitions as KML polygons.']);
    739     epart=md.part(md.elements)+1;
    740     if exist('ndata','var') || exist('edata','var')
    741         pdata=zeros(1,md.npart);
    742         pdata(:)=NaN;
    743     end
    744     for k=1:md.npart
    745         [icol,irow]=find(epart'==k);
    746         irow=unique(irow);
    747         elem=md.elements(irow,:);
    748         nodecon=nodeconnectivity(elem,md.numberofgrids);
    749         [edgeper,elemper,iloop]=edgeperimeter(elem,nodecon);
    750         iloop(end+1)=size(edgeper,1)+1;
    751         if exist('ndata','var')
    752             pdata(k)=ndata(find(md.part+1==k,1));
    753         elseif exist('edata','var')
    754             for i=1:size(epart,1)
    755                 if isempty(find(epart(i,:)~=k,1))
    756                     pdata(k)=edata(i);
    757                     break
    758                 end
    759             end
    760             if isnan(pdata(k))
    761                 warning('Data for Partition %d is not defined.\n',k)
    762             end
    763         end
    764         for i=1:length(iloop)-1
    765             fprintf(fid,'      <Placemark>\n');
    766             if (length(iloop)-1 > 1)
    767                 fprintf(fid,'        <name>Partition %d, Loop %d</name>\n',k,i);
    768             else
    769                 fprintf(fid,'        <name>Partition %d</name>\n',k);
    770             end
    771             fprintf(fid,'        <visibility>1</visibility>\n');
    772             if exist('pdata','var')
    773                 fprintf(fid,'        <description>Partition data: %g</description>\n',pdata(k));
    774                 imap = fix((pdata(k)-cmin)/(cmax-cmin)*size(cmap,1))+1;
    775                 if     (imap >= 1) && (imap <= size(cmap,1))
    776                     fprintf(fid,'        <styleUrl>#MatlabColor%d</styleUrl>\n',imap);
    777                 elseif (pdata(k) == cmax)
    778                     fprintf(fid,'        <styleUrl>#MatlabColor%d</styleUrl>\n',size(cmap,1));
    779                 else
    780                     fprintf(fid,'        <styleUrl>#BlackLineEmptyPoly</styleUrl>\n');
    781                 end
    782             else
    783                 fprintf(fid,'        <styleUrl>#BlackLineRandomPoly</styleUrl>\n');
    784             end
    785             fprintf(fid,'        <Polygon>\n');
    786             fprintf(fid,'          <extrude>1</extrude>\n');
    787             fprintf(fid,'          <altitudeMode>relativeToGround</altitudeMode>\n');
    788             fprintf(fid,'          <outerBoundaryIs>\n');
    789             fprintf(fid,'            <LinearRing>\n');
    790             fprintf(fid,'              <coordinates>\n');
    791             for j=iloop(i):iloop(i+1)-1
    792                 [lat(j),long(j)]=mapxy(md.x(edgeper(j,1)),md.y(edgeper(j,1)),'s');
    793                 fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(j),lat(j),alt);
    794             end
    795             fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(iloop(i)),lat(iloop(i)),alt);
    796 
    797             fprintf(fid,'              </coordinates>\n');
    798             fprintf(fid,'            </LinearRing>\n');
    799             fprintf(fid,'          </outerBoundaryIs>\n');
    800             fprintf(fid,'        </Polygon>\n');
    801             fprintf(fid,'      </Placemark>\n');
    802         end
    803     end
    804     fprintf(fid,'    </Folder>\n');
    805 end
    806 
    807 end
    808 
    809 %%
    810 %  create kml linestrings for the partition edges.
    811 %
    812 %  []=kml_part_edges(fid,md,params)
    813 %
    814 %  where the required input is:
    815 %    fid           (char, file ID of .kml file)
    816 %    md            (model, model class object)
    817 %
    818 %  the optional input is:
    819 %    params        (string/numeric, parameter names and values)
    820 %
    821 %  and the optional input is:
    822 %    data          (numeric, element or nodal results data)
    823 %    alt           (numeric, altitude for polygons, default 10000)
    824 %    cmin          (numeric, minimum of color map)
    825 %    cmax          (numeric, maximum of color map)
    826 %    cmap          (char or numeric, colormap definition)
    827 %    prtplt        (char, 'off'/'no' for partition segment plot)
    828 %
    829 function []=kml_part_edges(varargin)
    830 
    831 if ~nargin
    832     help kml_part_edges
    833     return
    834 end
    835 
    836 %%  process input data
    837 
    838 iarg=1;
    839 if (nargin >= 1)
    840     fid=varargin{1};
    841 end
    842 if ~exist('fid','var') || isempty(fid) || (fid < 0)
    843     error('File ID ''%d'' must be open.',fid);
    844 end
    845 
    846 iarg=iarg+1;
    847 if (nargin >= 2)
    848     md=varargin{2};
    849 end
    850 if ~exist('md','var') || isempty(md) || ~isa(md,'model')
    851     error(['Model ''' inputname(iarg) ''' is unrecognized class ''' class(md) '''.']);
    852 end
    853 
    854 %  parameters
    855 
    856 iarg=iarg+1;
    857 while (iarg <= nargin-1)
    858     if ischar(varargin{iarg})
    859         eval([varargin{iarg} '=varargin{iarg+1};']);
    860         if (numel(varargin{iarg+1}) <= 20)
    861             disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']);
    862         else
    863             disp([varargin{iarg} '=' string_size(varargin{iarg+1}) ' ' class(varargin{iarg+1}) ';']);
    864         end
    865         if strcmpi(varargin{iarg},'data')
    866             cdata=inputname(iarg+1);
    867         end
    868     else
    869         error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']);
    870     end
    871     iarg=iarg+2;
    872 end
    873 
    874 if exist('data','var') && ~isempty(data)
    875     if     (numel(data)==md.numberofelements)
    876         edata=data;
    877     elseif (numel(data)==md.numberofgrids)
    878         ndata=data;
    879         display('Averaging nodal data to element data.');
    880         edata=zeros(1,md.numberofelements);
    881         for i=1:size(md.elements,1)
    882             for j=1:size(md.elements,2)
    883                 edata(i)=edata(i)+ndata(md.elements(i,j));
    884             end
    885             edata(i)=edata(i)/size(md.elements,2);
    886         end
    887     else
    888         error(['Data has incorrect number of ' num2str(numel(data)) ' values.']);
    889     end
    890 end
    891 
    892 %  colormap command operates on a figure, so create an invisible one
    893 %  (could also directly call colormaps, e.g. jet(64), but risky)
    894 
    895 hfig=figure('Visible','off');
    896 if exist('cmap','var')
    897     colormap(cmap)
    898 end
    899 cmap=colormap;
    900 close(hfig)
    901    
    902 if exist('edata','var')
    903     if ~exist('cmin','var')
    904         cmin=min(min(edata));
    905     end
    906     if ~exist('cmax','var')
    907         cmax=max(max(edata));
    908     end
    909 end
    910 
    911 if ~exist('alt','var')
    912     alt=10000;
    913 end
    914 
    915 %%  write folder for partition edges
    916 
    917 if (~exist('prtplt','var') || strncmpi(prtplt,'on' ,2) || strncmpi(prtplt,'y',1)) && ...
    918     md.npart
    919     fprintf(fid,'    <Folder>\n');
    920     fprintf(fid,'      <name>Partition Edges</name>\n');
    921     fprintf(fid,'      <visibility>1</visibility>\n');
    922     fprintf(fid,'      <description>Partitions=%d, Nodes=%d</description>\n',...
    923         md.npart,md.numberofgrids);
    924 
    925 %  write each partition loop as linestrings
    926 
    927     disp(['Writing ' num2str(md.npart) ' partitions as KML linestrings.']);
    928     epart=md.part(md.elements)+1;
    929     if exist('ndata','var') || exist('edata','var')
    930         pdata=zeros(1,md.npart);
    931         pdata(:)=NaN;
    932     end
    933 
    934 %  loop over each partition
    935 
    936     for k=1:md.npart
    937         disp(['partition k=' int2str(k)])
    938        
    939 %  for each partition, find all the included elements and determine the
    940 %  perimeter (including those shared by another partition)
    941 
    942         [icol,irow]=find(epart'==k);
    943         irow=unique(irow);
    944         elemp=md.elements(irow,:);
    945         epartp=epart(irow,:);
    946         nodeconp=nodeconnectivity(elemp,md.numberofgrids);
    947         [edgeadjp]=edgeadjacency(elemp,nodeconp);
    948         [edgeper,elemper,iloop]=edgeperimeter(elemp,nodeconp,edgeadjp);
    949         iloop(end+1)=size(edgeper,1)+1;
    950        
    951 %  determine the data to be used for the colors (if any)
    952 
    953         if exist('ndata','var')
    954             pdata(k)=ndata(find(md.part+1==k,1));
    955         elseif exist('edata','var')
    956             for i=1:size(epartp,1)
    957                 if isempty(find(epart(i,:)~=k,1))
    958                     pdata(k)=edata(irow(i));
    959                     break
    960                 end
    961             end
    962             if isnan(pdata(k))
    963                 warning('Data for Partition %d is not defined.\n',k)
    964             end
    965         end
    966        
    967 %  loop over each loop of the perimeter for the given partition
    968 
    969         for i=1:length(iloop)-1
    970             fprintf(fid,'      <Placemark>\n');
    971             if (length(iloop)-1 > 1)
    972                 fprintf(fid,'        <name>Partition %d, Loop %d</name>\n',k,i);
    973             else
    974                 fprintf(fid,'        <name>Partition %d</name>\n',k);
    975             end
    976             fprintf(fid,'        <visibility>1</visibility>\n');
    977             if exist('pdata','var')
    978                 fprintf(fid,'        <description>Partition data: %g</description>\n',pdata(k));
    979                 imap = fix((pdata(k)-cmin)/(cmax-cmin)*size(cmap,1))+1;
    980                 if     (imap >= 1) && (imap <= size(cmap,1))
    981                     fprintf(fid,'        <styleUrl>#MatlabColor%d</styleUrl>\n',imap);
    982                 elseif (pdata(k) == cmax)
    983                     fprintf(fid,'        <styleUrl>#MatlabColor%d</styleUrl>\n',size(cmap,1));
    984                 else
    985                     fprintf(fid,'        <styleUrl>#BlackLineEmptyPoly</styleUrl>\n');
    986                 end
    987             else
    988                 fprintf(fid,'        <styleUrl>#BlackLineRandomPoly</styleUrl>\n');
    989             end
    990 %            fprintf(fid,'        <Polygon>\n');
    991 %            fprintf(fid,'          <extrude>1</extrude>\n');
    992 %            fprintf(fid,'          <altitudeMode>relativeToGround</altitudeMode>\n');
    993 %            fprintf(fid,'          <outerBoundaryIs>\n');
    994 %            fprintf(fid,'            <LinearRing>\n');
    995         fprintf(fid,'        <LineString>\n');
    996         fprintf(fid,'          <extrude>1</extrude>\n');
    997         fprintf(fid,'          <tessellate>1</tessellate>\n');
    998         fprintf(fid,'          <altitudeMode>relativeToGround</altitudeMode>\n');
    999             fprintf(fid,'              <coordinates>\n');
    1000             elast=0;
    1001             nlast=0;
    1002             slast=0;
    1003             lat=[];
    1004             long=[];
    1005            
    1006 %  loop over the element edges on the loop of the partition
    1007 
    1008             j=iloop(i);
    1009             while (j < iloop(i+1))
    1010 %  find which side of element is referenced in perimeter list
    1011                 for l=1:size(elemp,2)
    1012                     if ((elemp(elemper(j),l)          == edgeper(j,1)) && ...
    1013                         (elemp(elemper(j),mod(l,3)+1) == edgeper(j,2))) || ...
    1014                        ((elemp(elemper(j),l)          == edgeper(j,2)) && ...
    1015                         (elemp(elemper(j),mod(l,3)+1) == edgeper(j,1)))
    1016                         jedge=l;
    1017                         break
    1018                     end
    1019                 end
    1020 
    1021 %  check if element side connects nodes in partition
    1022                 if (epartp(elemper(j),jedge)          == k) && ...
    1023                    (epartp(elemper(j),mod(jedge,3)+1) == k)
    1024 %  write out specified element side
    1025                     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.'])
    1026 %  if first edge, write out first node
    1027                     if ~elast
    1028                         [lat(end+1),long(end+1)]=mapxy(md.x(edgeper(j,1)),md.y(edgeper(j,1)),'s');
    1029                         fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
    1030                     end
    1031                     [lat(end+1),long(end+1)]=mapxy(md.x(edgeper(j,2)),md.y(edgeper(j,2)),'s');
    1032                     fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
    1033                     elast=elemper(j);
    1034                     nlast=edgeper(j,2);
    1035                     slast=0;
    1036                     j=j+1;
    1037                    
    1038 %  element not entirely within partition, so figure out boundary
    1039                 else
    1040                     disp(['segment j=' int2str(j) ' from element ' int2str(elemper(j)) ' shared by other partitions.'])
    1041                     ielem=elemper(j);
    1042 
    1043 %  follow partition boundary through elements not wholly in partition
    1044 %  (may include elements not in perimeter list)
    1045 
    1046                     while 1
    1047 %  if first edge, figure out direction from perimeter edge direction
    1048                         if ~nlast && ~slast
    1049                             nlast=find(elemp(ielem,:)==edgeper(j,1));
    1050                             nnext=find(elemp(ielem,:)==edgeper(j,2));
    1051                             if     (nlast+nnext == 3)
    1052                                 slast=1;
    1053                             elseif (nlast+nnext == 5)
    1054                                 slast=2;
    1055                             elseif (nlast+nnext == 4)
    1056                                 slast=3;
    1057                             end
    1058                             if     (nnext+(6-nlast-nnext) == 3)
    1059                                 snext=1;
    1060                             elseif (nnext+(6-nlast-nnext) == 5)
    1061                                 snext=2;
    1062                             elseif (nnext+(6-nlast-nnext) == 4)
    1063                                 snext=3;
    1064                             end
    1065 
    1066 %  find how many nodes of current element are in current partition
    1067 %  (1 or 2, not 3) and handle each permutation separately
    1068                             ipart=find(epartp(ielem,:)==k);
    1069 %  two nodes are in current partition, so cut off other node
    1070                             if (length(ipart) == 2)
    1071                                 switch 6-sum(ipart)
    1072                                     case nlast
    1073                                         slast=6-slast-snext;
    1074                                         nlast=0;
    1075                                     case nnext
    1076                                         if (epartp(ielem,nnext) == k)
    1077                                             nlast=nnext;
    1078                                         end
    1079                                     otherwise
    1080                                         slast=6-slast-snext;
    1081                                         nlast=0;
    1082                                 end
    1083 %  one node is in current partition
    1084                             else
    1085 %  all different, so cut through centroid
    1086                                 if (epartp(ielem,1) ~= epartp(ielem,2)) && ...
    1087                                    (epartp(ielem,2) ~= epartp(ielem,3)) && ...
    1088                                    (epartp(ielem,3) ~= epartp(ielem,1))
    1089                                     switch ipart
    1090                                         case {nlast,nnext}
    1091                                             if (epartp(ielem,nnext) == k)
    1092                                                 nlast=nnext;
    1093                                             end
    1094                                         otherwise
    1095                                             slast=6-slast-snext;
    1096                                             nlast=0;
    1097                                     end
    1098 %  other two are in the same partition, so cut them off
    1099                                 else
    1100                                     switch ipart
    1101                                         case nlast
    1102                                             if (epartp(ielem,nnext) == k)
    1103                                                 nlast=nnext;
    1104                                             end
    1105                                         case nnext
    1106                                             slast=snext;
    1107                                             nlast=0;
    1108                                         otherwise
    1109                                             slast=6-slast-snext;
    1110                                             nlast=0;
    1111                                     end
    1112                                 end
    1113                             end
    1114 
    1115 %  last edge exited last element at node
    1116                             if nlast
    1117 %  write out first node of first side for half-edge to midpoint
    1118                                 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.'])
    1119                                 [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,nlast))),...
    1120                                                                (md.y(elemp(ielem,nlast))),'s');
    1121                                 fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
    1122                             end
    1123                             nlast=0;
    1124                            
    1125 %  write out midpoint of first side
    1126                             [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,slast))...
    1127                                                            +md.x(elemp(ielem,mod(slast,3)+1)))/2.,...
    1128                                                            (md.y(elemp(ielem,slast))...
    1129                                                            +md.y(elemp(ielem,mod(slast,3)+1)))/2.,'s');
    1130                             fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
    1131                         end
    1132 
    1133 %  last edge exited last element at node
    1134                         if nlast
    1135                             if elast
    1136 %  find where last node on previous element occurs on current element
    1137                                 nlast=find(elemp(ielem,:)==nlast,1);
    1138                             end
    1139 %  half-edge occurs on unshared side from current node (unique unless mesh
    1140 %  is only attached at node)
    1141                             switch nlast
    1142                                 case 1
    1143                                     if ~edgeadjp(ielem,1)
    1144                                         nnext=2;
    1145                                         slast=1;
    1146                                     else
    1147                                         nnext=3;
    1148                                         slast=3;
    1149                                     end
    1150                                 case 2
    1151                                     if ~edgeadjp(ielem,2)
    1152                                         nnext=3;
    1153                                         slast=2;
    1154                                     else
    1155                                         nnext=1;
    1156                                         slast=1;
    1157                                     end
    1158                                 case 3
    1159                                     if ~edgeadjp(ielem,3)
    1160                                         nnext=1;
    1161                                         slast=3;
    1162                                     else
    1163                                         nnext=2;
    1164                                         slast=2;
    1165                                     end
    1166                             end
    1167 %  write out half-edge from current node to midpoint of unshared side
    1168                             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.'])
    1169                             [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,nlast))...
    1170                                                            +md.x(elemp(ielem,nnext)))/2.,...
    1171                                                            (md.y(elemp(ielem,nlast))...
    1172                                                            +md.y(elemp(ielem,nnext)))/2.,'s');
    1173                             fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
    1174                             nlast=0;
    1175 
    1176 %  last edge exited last element at midpoint of side
    1177                         elseif slast
    1178                             if elast
    1179 %  find where last side on previous element occurs on current element
    1180                                 slast=find(edgeadjp(ielem,:)==elast,1);
    1181                             end
    1182                         end
    1183 
    1184 %  find how many nodes of current element are in current partition
    1185 %  (1 or 2, not 3) and handle each permutation separately
    1186                         ipart=find(epartp(ielem,:)==k);
    1187                         if (length(ipart) == 2)
    1188 %  two nodes are in current partition, so cut off other node
    1189                             switch 6-sum(ipart)
    1190                                 case 1
    1191                                     snext=3+1-slast;
    1192                                 case 2
    1193                                     snext=1+2-slast;
    1194                                 case 3
    1195                                     snext=2+3-slast;
    1196                             end
    1197                         else
    1198                             if (epartp(ielem,1) ~= epartp(ielem,2)) && ...
    1199                                (epartp(ielem,2) ~= epartp(ielem,3)) && ...
    1200                                (epartp(ielem,3) ~= epartp(ielem,1))
    1201 %  all different, so cut through centroid
    1202                                 disp(['element ielem=' int2str(ielem) ' centroid written.'])
    1203                                 [lat(end+1),long(end+1)]=mapxy(sum(md.x(elemp(ielem,:)))/3.,...
    1204                                                                sum(md.y(elemp(ielem,:)))/3.,'s');
    1205                                 fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
    1206                             end
    1207 %  one node is in current partition, so cut off other two nodes
    1208                             switch ipart
    1209                                 case 1
    1210                                     snext=3+1-slast;
    1211                                 case 2
    1212                                     snext=1+2-slast;
    1213                                 case 3
    1214                                     snext=2+3-slast;
    1215                             end
    1216                         end
    1217 %  write out midpoint of opposite side
    1218                         disp(['segment j=' int2str(j) ' internal edge from side ' int2str(slast) ' to side ' int2str(snext) ' from element ' int2str(ielem) ' written.'])
    1219                         [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,snext))...
    1220                                                        +md.x(elemp(ielem,mod(snext,3)+1)))/2.,...
    1221                                                        (md.y(elemp(ielem,snext))...
    1222                                                        +md.y(elemp(ielem,mod(snext,3)+1)))/2.,'s');
    1223                         fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
    1224                         elast=ielem;
    1225                         nlast=0;
    1226                         slast=snext;
    1227 %  find adjacent element to opposite side
    1228                         ielem=edgeadjp(elast,slast);
    1229 %  if opposite side is unshared, find it in edge perimeter list
    1230                         if ~ielem
    1231                             jlast=find(elemper(j:end)==elast)+j-1;
    1232                             j=0;
    1233                             for l=1:length(jlast)
    1234                                 if ((elemp(elast,slast)          == edgeper(jlast(l),1)) && ...
    1235                                     (elemp(elast,mod(slast,3)+1) == edgeper(jlast(l),2))) || ...
    1236                                    ((elemp(elast,slast)          == edgeper(jlast(l),2)) && ...
    1237                                     (elemp(elast,mod(slast,3)+1) == edgeper(jlast(l),1)))
    1238                                     j=jlast(l);
    1239                                     break
    1240                                 end
    1241                             end
    1242                             if ~j
    1243                                 j=iloop(i+1)-1;
    1244                             end
    1245 %  write out half-edge from midpoint of unshared side to node
    1246                             if (epartp(elast,slast) == k)
    1247                                 nnext=slast;
    1248                             else
    1249                                 nnext=mod(slast,3)+1;
    1250                             end
    1251                             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.'])
    1252                             [lat(end+1),long(end+1)]=mapxy(md.x(elemp(elast,nnext)),...
    1253                                                            md.y(elemp(elast,nnext)),'s');
    1254                             fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
    1255                             break
    1256 %  if not unshared, advance perimeter list and watch for end
    1257                         else
    1258                             if (elast == elemper(j))
    1259                                 if (j+1 < iloop(i+1)) && ...
    1260                                    ~isempty(find(elemper(j+1:end)~=elast,1))
    1261                                     j=j+find(elemper(j+1:end)~=elast,1);
    1262                                 else
    1263                                     break
    1264                                 end
    1265                             end
    1266                         end
    1267                     end
    1268                     j=j+1;
    1269                 end
    1270             end
    1271 %            fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(iloop(i)),lat(iloop(i)),alt);
    1272 
    1273             fprintf(fid,'              </coordinates>\n');
    1274         fprintf(fid,'        </LineString>\n');
    1275 %            fprintf(fid,'            </LinearRing>\n');
    1276 %            fprintf(fid,'          </outerBoundaryIs>\n');
    1277 %            fprintf(fid,'        </Polygon>\n');
    1278             fprintf(fid,'      </Placemark>\n');
    1279         end
    1280     end
    1281     fprintf(fid,'    </Folder>\n');
    1282 end
    1283 
    1284 end
    1285 
    1286 %%
    1287 %  create kml polygons for the partitions.
    1288 %
    1289 %  []=kml_partitions(fid,md,params)
    1290 %
    1291 %  where the required input is:
    1292 %    fid           (char, file ID of .kml file)
    1293 %    md            (model, model class object)
    1294 %
    1295 %  the optional input is:
    1296 %    params        (string/numeric, parameter names and values)
    1297 %
    1298 %  and the optional input is:
    1299 %    data          (numeric, element or nodal results data)
    1300 %    alt           (numeric, altitude for polygons, default 10000)
    1301 %    cmin          (numeric, minimum of color map)
    1302 %    cmax          (numeric, maximum of color map)
    1303 %    cmap          (char or numeric, colormap definition)
    1304 %    prtplt        (char, 'off'/'no' for partition segment plot)
    1305 %
    1306 function []=kml_partitions(varargin)
    1307 
    1308 if ~nargin
    1309     help kml_part_edges
    1310     return
    1311 end
    1312 
    1313 %%  process input data
    1314 
    1315 iarg=1;
    1316 if (nargin >= 1)
    1317     fid=varargin{1};
    1318 end
    1319 if ~exist('fid','var') || isempty(fid) || (fid < 0)
    1320     error('File ID ''%d'' must be open.',fid);
    1321 end
    1322 
    1323 iarg=iarg+1;
    1324 if (nargin >= 2)
    1325     md=varargin{2};
    1326 end
    1327 if ~exist('md','var') || isempty(md) || ~isa(md,'model')
    1328     error(['Model ''' inputname(iarg) ''' is unrecognized class ''' class(md) '''.']);
    1329 end
    1330 
    1331 %  parameters
    1332 
    1333 iarg=iarg+1;
    1334 while (iarg <= nargin-1)
    1335     if ischar(varargin{iarg})
    1336         eval([varargin{iarg} '=varargin{iarg+1};']);
    1337         if (numel(varargin{iarg+1}) <= 20)
    1338             disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']);
    1339         else
    1340             disp([varargin{iarg} '=' string_size(varargin{iarg+1}) ' ' class(varargin{iarg+1}) ';']);
    1341         end
    1342         if strcmpi(varargin{iarg},'data')
    1343             cdata=inputname(iarg+1);
    1344         end
    1345     else
    1346         error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']);
    1347     end
    1348     iarg=iarg+2;
    1349 end
    1350 
    1351 if exist('data','var') && ~isempty(data)
    1352     if     (numel(data)==md.numberofelements)
    1353         edata=data;
    1354     elseif (numel(data)==md.numberofgrids)
    1355         ndata=data;
    1356         display('Averaging nodal data to element data.');
    1357         edata=zeros(1,md.numberofelements);
    1358         for i=1:size(md.elements,1)
    1359             for j=1:size(md.elements,2)
    1360                 edata(i)=edata(i)+ndata(md.elements(i,j));
    1361             end
    1362             edata(i)=edata(i)/size(md.elements,2);
    1363         end
    1364     else
    1365         error(['Data has incorrect number of ' num2str(numel(data)) ' values.']);
    1366     end
    1367 end
    1368 
    1369 %  colormap command operates on a figure, so create an invisible one
    1370 %  (could also directly call colormaps, e.g. jet(64), but risky)
    1371 
    1372 hfig=figure('Visible','off');
    1373 if exist('cmap','var')
    1374     colormap(cmap)
    1375 end
    1376 cmap=colormap;
    1377 close(hfig)
    1378    
    1379 if exist('edata','var')
    1380     if ~exist('cmin','var')
    1381         cmin=min(min(edata));
    1382     end
    1383     if ~exist('cmax','var')
    1384         cmax=max(max(edata));
    1385     end
    1386 end
    1387 
    1388 if ~exist('alt','var')
    1389     alt=10000;
    1390 end
    1391 
    1392 %%  write folder for partitions
    1393 
    1394 if (~exist('prtplt','var') || strncmpi(prtplt,'on' ,2) || strncmpi(prtplt,'y',1)) && ...
    1395     md.npart
    1396     fprintf(fid,'    <Folder>\n');
    1397     fprintf(fid,'      <name>Partitions</name>\n');
    1398     fprintf(fid,'      <visibility>1</visibility>\n');
    1399     fprintf(fid,'      <description>Partitions=%d, Nodes=%d</description>\n',...
    1400         md.npart,md.numberofgrids);
    1401 
    1402 %  write each partition loop as linestrings
    1403 
    1404     disp(['Writing ' num2str(md.npart) ' partitions as KML polygons.']);
    1405     epart=md.part(md.elements)+1;
    1406     if exist('ndata','var') || exist('edata','var')
    1407         pdata=zeros(1,md.npart);
    1408         pdata(:)=NaN;
    1409     end
    1410 
    1411 %  loop over each partition
    1412 
    1413     for k=1:md.npart
    1414         disp(['partition k=' int2str(k)])
    1415        
    1416 %  for each partition, find all the included elements and determine the
    1417 %  perimeter (including those shared by another partition)
    1418 
    1419         [icol,irow]=find(epart'==k);
    1420         irow=unique(irow);
    1421         elemp=md.elements(irow,:);
    1422         epartp=epart(irow,:);
    1423         nodeconp=nodeconnectivity(elemp,md.numberofgrids);
    1424         [edgeadjp]=edgeadjacency(elemp,nodeconp);
    1425         [edgeper,elemper,iloop]=edgeperimeter(elemp,nodeconp,edgeadjp);
    1426         iloop(end+1)=size(edgeper,1)+1;
    1427        
    1428 %  determine the data to be used for the colors (if any)
    1429 
    1430         if exist('ndata','var')
    1431             pdata(k)=ndata(find(md.part+1==k,1));
    1432         elseif exist('edata','var')
    1433             for i=1:size(epartp,1)
    1434                 if isempty(find(epart(i,:)~=k,1))
    1435                     pdata(k)=edata(irow(i));
    1436                     break
    1437                 end
    1438             end
    1439             if isnan(pdata(k))
    1440                 warning('Data for Partition %d is not defined.\n',k)
    1441             end
    1442         end
    1443        
    1444 %  loop over each loop of the perimeter for the given partition
    1445 
    1446         for i=1:length(iloop)-1
    1447             fprintf(fid,'      <Placemark>\n');
    1448             if (length(iloop)-1 > 1)
    1449                 fprintf(fid,'        <name>Partition %d, Loop %d</name>\n',k,i);
    1450             else
    1451                 fprintf(fid,'        <name>Partition %d</name>\n',k);
    1452             end
    1453             fprintf(fid,'        <visibility>1</visibility>\n');
    1454             if exist('pdata','var')
    1455                 fprintf(fid,'        <description>Partition data: %g</description>\n',pdata(k));
    1456                 imap = fix((pdata(k)-cmin)/(cmax-cmin)*size(cmap,1))+1;
    1457                 if     (imap >= 1) && (imap <= size(cmap,1))
    1458                     fprintf(fid,'        <styleUrl>#MatlabColor%d</styleUrl>\n',imap);
    1459                 elseif (pdata(k) == cmax)
    1460                     fprintf(fid,'        <styleUrl>#MatlabColor%d</styleUrl>\n',size(cmap,1));
    1461                 else
    1462                     fprintf(fid,'        <styleUrl>#BlackLineEmptyPoly</styleUrl>\n');
    1463                 end
    1464             else
    1465                 fprintf(fid,'        <styleUrl>#BlackLineRandomPoly</styleUrl>\n');
    1466             end
    1467             fprintf(fid,'        <Polygon>\n');
    1468             fprintf(fid,'          <extrude>1</extrude>\n');
    1469             fprintf(fid,'          <altitudeMode>relativeToGround</altitudeMode>\n');
    1470             fprintf(fid,'          <outerBoundaryIs>\n');
    1471             fprintf(fid,'            <LinearRing>\n');
    1472             fprintf(fid,'              <coordinates>\n');
    1473             elast=0;
    1474             nlast=0;
    1475             slast=0;
    1476             lat=[];
    1477             long=[];
    1478            
    1479 %  loop over the element edges on the loop of the partition
    1480 
    1481             j=iloop(i);
    1482             while (j < iloop(i+1))
    1483 %  find which side of element is referenced in perimeter list
    1484                 for l=1:size(elemp,2)
    1485                     if ((elemp(elemper(j),l)          == edgeper(j,1)) && ...
    1486                         (elemp(elemper(j),mod(l,3)+1) == edgeper(j,2))) || ...
    1487                        ((elemp(elemper(j),l)          == edgeper(j,2)) && ...
    1488                         (elemp(elemper(j),mod(l,3)+1) == edgeper(j,1)))
    1489                         jedge=l;
    1490                         break
    1491                     end
    1492                 end
    1493 
    1494 %  check if element side connects nodes in partition
    1495                 if (epartp(elemper(j),jedge)          == k) && ...
    1496                    (epartp(elemper(j),mod(jedge,3)+1) == k)
    1497 %  write out specified element side
    1498                     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.'])
    1499 %  if first edge, write out first node
    1500                     if ~elast
    1501                         [lat(end+1),long(end+1)]=mapxy(md.x(edgeper(j,1)),md.y(edgeper(j,1)),'s');
    1502                         fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
    1503                     end
    1504                     [lat(end+1),long(end+1)]=mapxy(md.x(edgeper(j,2)),md.y(edgeper(j,2)),'s');
    1505                     fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
    1506                     elast=elemper(j);
    1507                     nlast=edgeper(j,2);
    1508                     slast=0;
    1509                     j=j+1;
    1510                    
    1511 %  element not entirely within partition, so figure out boundary
    1512                 else
    1513                     disp(['segment j=' int2str(j) ' from element ' int2str(elemper(j)) ' shared by other partitions.'])
    1514                     ielem=elemper(j);
    1515 
    1516 %  follow partition boundary through elements not wholly in partition
    1517 %  (may include elements not in perimeter list)
    1518 
    1519                     while 1
    1520 %  if first edge, figure out direction from perimeter edge direction
    1521                         if ~nlast && ~slast
    1522                             nlast=find(elemp(ielem,:)==edgeper(j,1));
    1523                             nnext=find(elemp(ielem,:)==edgeper(j,2));
    1524                             if     (nlast+nnext == 3)
    1525                                 slast=1;
    1526                             elseif (nlast+nnext == 5)
    1527                                 slast=2;
    1528                             elseif (nlast+nnext == 4)
    1529                                 slast=3;
    1530                             end
    1531                             if     (nnext+(6-nlast-nnext) == 3)
    1532                                 snext=1;
    1533                             elseif (nnext+(6-nlast-nnext) == 5)
    1534                                 snext=2;
    1535                             elseif (nnext+(6-nlast-nnext) == 4)
    1536                                 snext=3;
    1537                             end
    1538 
    1539 %  find how many nodes of current element are in current partition
    1540 %  (1 or 2, not 3) and handle each permutation separately
    1541                             ipart=find(epartp(ielem,:)==k);
    1542 %  two nodes are in current partition, so cut off other node
    1543                             if (length(ipart) == 2)
    1544                                 switch 6-sum(ipart)
    1545                                     case nlast
    1546                                         slast=6-slast-snext;
    1547                                         nlast=0;
    1548                                     case nnext
    1549                                         if (epartp(ielem,nnext) == k)
    1550                                             nlast=nnext;
    1551                                         end
    1552                                     otherwise
    1553                                         slast=6-slast-snext;
    1554                                         nlast=0;
    1555                                 end
    1556 %  one node is in current partition
    1557                             else
    1558 %  all different, so cut through centroid
    1559                                 if (epartp(ielem,1) ~= epartp(ielem,2)) && ...
    1560                                    (epartp(ielem,2) ~= epartp(ielem,3)) && ...
    1561                                    (epartp(ielem,3) ~= epartp(ielem,1))
    1562                                     switch ipart
    1563                                         case {nlast,nnext}
    1564                                             if (epartp(ielem,nnext) == k)
    1565                                                 nlast=nnext;
    1566                                             end
    1567                                         otherwise
    1568                                             slast=6-slast-snext;
    1569                                             nlast=0;
    1570                                     end
    1571 %  other two are in the same partition, so cut them off
    1572                                 else
    1573                                     switch ipart
    1574                                         case nlast
    1575                                             if (epartp(ielem,nnext) == k)
    1576                                                 nlast=nnext;
    1577                                             end
    1578                                         case nnext
    1579                                             slast=snext;
    1580                                             nlast=0;
    1581                                         otherwise
    1582                                             slast=6-slast-snext;
    1583                                             nlast=0;
    1584                                     end
    1585                                 end
    1586                             end
    1587 
    1588 %  last edge exited last element at node
    1589                             if nlast
    1590 %  write out first node of first side for half-edge to midpoint
    1591                                 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.'])
    1592                                 [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,nlast))),...
    1593                                                                (md.y(elemp(ielem,nlast))),'s');
    1594                                 fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
    1595                             end
    1596                             nlast=0;
    1597                            
    1598 %  write out midpoint of first side
    1599                             [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,slast))...
    1600                                                            +md.x(elemp(ielem,mod(slast,3)+1)))/2.,...
    1601                                                            (md.y(elemp(ielem,slast))...
    1602                                                            +md.y(elemp(ielem,mod(slast,3)+1)))/2.,'s');
    1603                             fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
    1604                         end
    1605 
    1606 %  last edge exited last element at node
    1607                         if nlast
    1608                             if elast
    1609 %  find where last node on previous element occurs on current element
    1610                                 nlast=find(elemp(ielem,:)==nlast,1);
    1611                             end
    1612 %  half-edge occurs on unshared side from current node (unique unless mesh
    1613 %  is only attached at node)
    1614                             switch nlast
    1615                                 case 1
    1616                                     if ~edgeadjp(ielem,1)
    1617                                         nnext=2;
    1618                                         slast=1;
    1619                                     else
    1620                                         nnext=3;
    1621                                         slast=3;
    1622                                     end
    1623                                 case 2
    1624                                     if ~edgeadjp(ielem,2)
    1625                                         nnext=3;
    1626                                         slast=2;
    1627                                     else
    1628                                         nnext=1;
    1629                                         slast=1;
    1630                                     end
    1631                                 case 3
    1632                                     if ~edgeadjp(ielem,3)
    1633                                         nnext=1;
    1634                                         slast=3;
    1635                                     else
    1636                                         nnext=2;
    1637                                         slast=2;
    1638                                     end
    1639                             end
    1640 %  write out half-edge from current node to midpoint of unshared side
    1641                             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.'])
    1642                             [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,nlast))...
    1643                                                            +md.x(elemp(ielem,nnext)))/2.,...
    1644                                                            (md.y(elemp(ielem,nlast))...
    1645                                                            +md.y(elemp(ielem,nnext)))/2.,'s');
    1646                             fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
    1647                             nlast=0;
    1648 
    1649 %  last edge exited last element at midpoint of side
    1650                         elseif slast
    1651                             if elast
    1652 %  find where last side on previous element occurs on current element
    1653                                 slast=find(edgeadjp(ielem,:)==elast,1);
    1654                             end
    1655                         end
    1656 
    1657 %  find how many nodes of current element are in current partition
    1658 %  (1 or 2, not 3) and handle each permutation separately
    1659                         ipart=find(epartp(ielem,:)==k);
    1660                         if (length(ipart) == 2)
    1661 %  two nodes are in current partition, so cut off other node
    1662                             switch 6-sum(ipart)
    1663                                 case 1
    1664                                     snext=3+1-slast;
    1665                                 case 2
    1666                                     snext=1+2-slast;
    1667                                 case 3
    1668                                     snext=2+3-slast;
    1669                             end
    1670                         else
    1671                             if (epartp(ielem,1) ~= epartp(ielem,2)) && ...
    1672                                (epartp(ielem,2) ~= epartp(ielem,3)) && ...
    1673                                (epartp(ielem,3) ~= epartp(ielem,1))
    1674 %  all different, so cut through centroid
    1675                                 disp(['element ielem=' int2str(ielem) ' centroid written.'])
    1676                                 [lat(end+1),long(end+1)]=mapxy(sum(md.x(elemp(ielem,:)))/3.,...
    1677                                                                sum(md.y(elemp(ielem,:)))/3.,'s');
    1678                                 fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
    1679                             end
    1680 %  one node is in current partition, so cut off other two nodes
    1681                             switch ipart
    1682                                 case 1
    1683                                     snext=3+1-slast;
    1684                                 case 2
    1685                                     snext=1+2-slast;
    1686                                 case 3
    1687                                     snext=2+3-slast;
    1688                             end
    1689                         end
    1690 %  write out midpoint of opposite side
    1691                         disp(['segment j=' int2str(j) ' internal edge from side ' int2str(slast) ' to side ' int2str(snext) ' from element ' int2str(ielem) ' written.'])
    1692                         [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,snext))...
    1693                                                        +md.x(elemp(ielem,mod(snext,3)+1)))/2.,...
    1694                                                        (md.y(elemp(ielem,snext))...
    1695                                                        +md.y(elemp(ielem,mod(snext,3)+1)))/2.,'s');
    1696                         fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
    1697                         elast=ielem;
    1698                         nlast=0;
    1699                         slast=snext;
    1700 %  find adjacent element to opposite side
    1701                         ielem=edgeadjp(elast,slast);
    1702 %  if opposite side is unshared, find it in edge perimeter list
    1703                         if ~ielem
    1704                             jlast=find(elemper(j:end)==elast)+j-1;
    1705                             j=0;
    1706                             for l=1:length(jlast)
    1707                                 if ((elemp(elast,slast)          == edgeper(jlast(l),1)) && ...
    1708                                     (elemp(elast,mod(slast,3)+1) == edgeper(jlast(l),2))) || ...
    1709                                    ((elemp(elast,slast)          == edgeper(jlast(l),2)) && ...
    1710                                     (elemp(elast,mod(slast,3)+1) == edgeper(jlast(l),1)))
    1711                                     j=jlast(l);
    1712                                     break
    1713                                 end
    1714                             end
    1715                             if ~j
    1716                                 j=iloop(i+1)-1;
    1717                             end
    1718 %  write out half-edge from midpoint of unshared side to node
    1719                             if (epartp(elast,slast) == k)
    1720                                 nnext=slast;
    1721                             else
    1722                                 nnext=mod(slast,3)+1;
    1723                             end
    1724                             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.'])
    1725                             [lat(end+1),long(end+1)]=mapxy(md.x(elemp(elast,nnext)),...
    1726                                                            md.y(elemp(elast,nnext)),'s');
    1727                             fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
    1728                             break
    1729 %  if not unshared, advance perimeter list and watch for end
    1730                         else
    1731                             if (elast == elemper(j))
    1732                                 if (j+1 < iloop(i+1)) && ...
    1733                                    ~isempty(find(elemper(j+1:end)~=elast,1))
    1734                                     j=j+find(elemper(j+1:end)~=elast,1);
    1735                                 else
    1736                                     break
    1737                                 end
    1738                             end
    1739                         end
    1740                     end
    1741                     j=j+1;
    1742                 end
    1743             end
    1744 %            fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(iloop(1)),lat(iloop(1)),alt);
    1745 
    1746             fprintf(fid,'              </coordinates>\n');
    1747             fprintf(fid,'            </LinearRing>\n');
    1748             fprintf(fid,'          </outerBoundaryIs>\n');
    1749             fprintf(fid,'        </Polygon>\n');
    1750             fprintf(fid,'      </Placemark>\n');
    1751         end
    1752     end
    1753     fprintf(fid,'    </Folder>\n');
    1754 end
    1755 
    1756 end
    1757 
     253    disp(['End of file ''' filek2 ''' successfully written.']);
     254end
     255
     256end
     257
Note: See TracChangeset for help on using the changeset viewer.