Changeset 6415

Timestamp:

10/25/10 11:53:59 (14 years ago)

Author:

jschierm

Message:

Addition of partition polygons to kml_mesh_write.m and reorganization into functions.

File:

• issm/trunk/src/m/kml/kml_mesh_write.m (modified) (23 diffs)

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

@@ -93 +93 @@
         error(['Data has incorrect number of ' num2str(numel(data)) ' values.']);
     end
+end
+
+if exist('edata','var')
+    if ~exist('cmin','var')
+        cmin=min(min(edata));
+    end
+    if ~exist('cmax','var')
+        cmax=max(max(edata));
+    end
+end
+
+if ~exist('alt','var')
+    alt=10000;
 end
 
@@ -193 +206 @@
 %  write folder for mesh
 
+kml_mesh_elem(fid,md,varargin{3:end});
+
+%  write folder for partition segments
+
+if (~exist('prtplt','var') || strncmpi(prtplt,'on' ,2) || strncmpi(prtplt,'y',1)) && ...
+    md.npart
+    kml_part_flagedges(fid,md,varargin{3:end});
+end
+
+%  write folder for unshared edges
+
+if (~exist('prtplt','var') || strncmpi(prtplt,'on' ,2) || strncmpi(prtplt,'y',1)) && ...
+    md.npart
+    kml_unsh_edges(fid,md,varargin{3:end});
+end
+
+%  write folder for partition elements
+
+if (~exist('prtplt','var') || strncmpi(prtplt,'on' ,2) || strncmpi(prtplt,'y',1)) && ...
+    md.npart
+    kml_part_elems(fid,md,varargin{3:end});
+end
+
+%  write folder for partition edges
+
+if (~exist('prtplt','var') || strncmpi(prtplt,'on' ,2) || strncmpi(prtplt,'y',1)) && ...
+    md.npart
+    kml_part_edges(fid,md,varargin{3:end});
+end
+
+%  write folder for partitions
+
+if (~exist('prtplt','var') || strncmpi(prtplt,'on' ,2) || strncmpi(prtplt,'y',1)) && ...
+    md.npart
+    kml_partitions(fid,md,varargin{3:end});
+end
+
+%  write trailer data
+
+fprintf(fid,'  </Document>\n');
+fprintf(fid,'</kml>\n');
+
+fclose(fid);
+display('End of file successfully written.');
+
+end
+
+%%
+%  create kml polygons for the element mesh.
+%
+%  []=kml_mesh_elem(fid,md,params)
+%
+%  where the required input is:
+%    fid           (numeric, file ID of .kml file)
+%    md            (model, model class object)
+%
+%  the optional input is:
+%    params        (string/numeric, parameter names and values)
+%
+%  and the optional input is:
+%    data          (numeric, element or nodal results data)
+%    alt           (numeric, altitude for polygons, default 10000)
+%    cmin          (numeric, minimum of color map)
+%    cmax          (numeric, maximum of color map)
+%    cmap          (char or numeric, colormap definition)
+%
+function []=kml_mesh_elem(varargin)
+
+if ~nargin
+    help kml_mesh_elem
+    return
+end
+
+%%  process input data
+
+iarg=1;
+if (nargin >= 1)
+    fid=varargin{1};
+end
+if ~exist('fid','var') || isempty(fid) || (fid < 0)
+    error('File ID ''%d'' must be open.',fid);
+end
+
+iarg=iarg+1;
+if (nargin >= 2)
+    md=varargin{2};
+end
+if ~exist('md','var') || isempty(md) || ~isa(md,'model')
+    error(['Model ''' inputname(iarg) ''' is unrecognized class ''' class(md) '''.']);
+end
+
+%  parameters
+
+iarg=iarg+1;
+while (iarg <= nargin-1)
+    if ischar(varargin{iarg})
+        eval([varargin{iarg} '=varargin{iarg+1};']);
+        if (numel(varargin{iarg+1}) <= 20)
+            disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']);
+        else
+            disp([varargin{iarg} '=' string_size(varargin{iarg+1}) ' ' class(varargin{iarg+1}) ';']);
+        end
+        if strcmpi(varargin{iarg},'data')
+            cdata=inputname(iarg+1);
+        end
+    else
+        error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']);
+    end
+    iarg=iarg+2;
+end
+
+if exist('data','var') && ~isempty(data)
+    if     (numel(data)==md.numberofelements)
+        edata=data;
+    elseif (numel(data)==md.numberofgrids)
+        ndata=data;
+        display('Averaging nodal data to element data.');
+        edata=zeros(1,md.numberofelements);
+        for i=1:size(md.elements,1)
+            for j=1:size(md.elements,2)
+                edata(i)=edata(i)+ndata(md.elements(i,j));
+            end
+            edata(i)=edata(i)/size(md.elements,2);
+        end
+    else
+        error(['Data has incorrect number of ' num2str(numel(data)) ' values.']);
+    end
+end
+
+%  colormap command operates on a figure, so create an invisible one
+%  (could also directly call colormaps, e.g. jet(64), but risky)
+
+hfig=figure('Visible','off');
+if exist('cmap','var')
+    colormap(cmap)
+end
+cmap=colormap;
+close(hfig)
+    
+if exist('edata','var')
+    if ~exist('cmin','var')
+        cmin=min(min(edata));
+    end
+    if ~exist('cmax','var')
+        cmax=max(max(edata));
+    end
+end
+
+if ~exist('alt','var')
+    alt=10000;
+end
+
+%%  write folder for mesh
+
 fprintf(fid,'    <Folder>\n');
 if exist('cdata','var') && ~isempty(cdata)
@@ -204 +371 @@
 
 %  write each element as a polygon
-
-if exist('edata','var')
-    if ~exist('cmin','var')
-        cmin=min(min(edata));
-    end
-    if ~exist('cmax','var')
-        cmax=max(max(edata));
-    end
-end
-
-if ~exist('alt','var')
-    alt=10000;
-end
 
 disp(['Writing ' num2str(size(md.elements,1)) ' tria elements as KML polygons.']);
@@ -256 +410 @@
 fprintf(fid,'    </Folder>\n');
 
-%  write folder for partition segments
+end
+
+%%
+%  create kml linestrings for the flagged partition edges.
+%
+%  []=kml_part_flagedges(fid,md,params)
+%
+%  where the required input is:
+%    fid           (char, file ID of .kml file)
+%    md            (model, model class object)
+%
+%  the optional input is:
+%    params        (string/numeric, parameter names and values)
+%
+%  and the optional input is:
+%    alt           (numeric, altitude for polygons, default 10000)
+%    prtplt        (char, 'off'/'no' for partition segment plot)
+%
+function []=kml_part_flagedges(varargin)
+
+if ~nargin
+    help kml_part_flagedges
+    return
+end
+
+%%  process input data
+
+iarg=1;
+if (nargin >= 1)
+    fid=varargin{1};
+end
+if ~exist('fid','var') || isempty(fid) || (fid < 0)
+    error('File ID ''%d'' must be open.',fid);
+end
+
+iarg=iarg+1;
+if (nargin >= 2)
+    md=varargin{2};
+end
+if ~exist('md','var') || isempty(md) || ~isa(md,'model')
+    error(['Model ''' inputname(iarg) ''' is unrecognized class ''' class(md) '''.']);
+end
+
+%  parameters
+
+iarg=iarg+1;
+while (iarg <= nargin-1)
+    if ischar(varargin{iarg})
+        eval([varargin{iarg} '=varargin{iarg+1};']);
+        if (numel(varargin{iarg+1}) <= 20)
+            disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']);
+        else
+            disp([varargin{iarg} '=' string_size(varargin{iarg+1}) ' ' class(varargin{iarg+1}) ';']);
+        end
+    else
+        error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']);
+    end
+    iarg=iarg+2;
+end
+
+if ~exist('alt','var')
+    alt=10000;
+end
+
+%%  write folder for partition segments
 
 if (~exist('prtplt','var') || strncmpi(prtplt,'on' ,2) || strncmpi(prtplt,'y',1)) && ...
@@ -292 +510 @@
 end
 
-%  write folder for unshared edges
+end
+
+%%
+%  create kml linestrings for the unshared element edges.
+%
+%  []=kml_unsh_edges(fid,md,params)
+%
+%  where the required input is:
+%    fid           (char, file ID of .kml file)
+%    md            (model, model class object)
+%
+%  the optional input is:
+%    params        (string/numeric, parameter names and values)
+%
+%  and the optional input is:
+%    alt           (numeric, altitude for polygons, default 10000)
+%    prtplt        (char, 'off'/'no' for partition segment plot)
+%
+function []=kml_unsh_edges(varargin)
+
+if ~nargin
+    help kml_unsh_edges
+    return
+end
+
+%%  process input data
+
+iarg=1;
+if (nargin >= 1)
+    fid=varargin{1};
+end
+if ~exist('fid','var') || isempty(fid) || (fid < 0)
+    error('File ID ''%d'' must be open.',fid);
+end
+
+iarg=iarg+1;
+if (nargin >= 2)
+    md=varargin{2};
+end
+if ~exist('md','var') || isempty(md) || ~isa(md,'model')
+    error(['Model ''' inputname(iarg) ''' is unrecognized class ''' class(md) '''.']);
+end
+
+%  parameters
+
+iarg=iarg+1;
+while (iarg <= nargin-1)
+    if ischar(varargin{iarg})
+        eval([varargin{iarg} '=varargin{iarg+1};']);
+        if (numel(varargin{iarg+1}) <= 20)
+            disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']);
+        else
+            disp([varargin{iarg} '=' string_size(varargin{iarg+1}) ' ' class(varargin{iarg+1}) ';']);
+        end
+    else
+        error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']);
+    end
+    iarg=iarg+2;
+end
+
+if ~exist('alt','var')
+    alt=10000;
+end
+
+%%  write folder for unshared edges
 
 if (~exist('prtplt','var') || strncmpi(prtplt,'on' ,2) || strncmpi(prtplt,'y',1)) && ...
@@ -334 +616 @@
 end
 
-%  write folder for partitions
+end
+
+%%
+%  create kml polygons for the partition elements.
+%
+%  []=kml_part_elems(fid,md,params)
+%
+%  where the required input is:
+%    fid           (char, file ID of .kml file)
+%    md            (model, model class object)
+%
+%  the optional input is:
+%    params        (string/numeric, parameter names and values)
+%
+%  and the optional input is:
+%    data          (numeric, element or nodal results data)
+%    alt           (numeric, altitude for polygons, default 10000)
+%    cmin          (numeric, minimum of color map)
+%    cmax          (numeric, maximum of color map)
+%    cmap          (char or numeric, colormap definition)
+%    prtplt        (char, 'off'/'no' for partition segment plot)
+%
+function []=kml_part_elems(varargin)
+
+if ~nargin
+    help kml_part_elems
+    return
+end
+
+%%  process input data
+
+iarg=1;
+if (nargin >= 1)
+    fid=varargin{1};
+end
+if ~exist('fid','var') || isempty(fid) || (fid < 0)
+    error('File ID ''%d'' must be open.',fid);
+end
+
+iarg=iarg+1;
+if (nargin >= 2)
+    md=varargin{2};
+end
+if ~exist('md','var') || isempty(md) || ~isa(md,'model')
+    error(['Model ''' inputname(iarg) ''' is unrecognized class ''' class(md) '''.']);
+end
+
+%  parameters
+
+iarg=iarg+1;
+while (iarg <= nargin-1)
+    if ischar(varargin{iarg})
+        eval([varargin{iarg} '=varargin{iarg+1};']);
+        if (numel(varargin{iarg+1}) <= 20)
+            disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']);
+        else
+            disp([varargin{iarg} '=' string_size(varargin{iarg+1}) ' ' class(varargin{iarg+1}) ';']);
+        end
+        if strcmpi(varargin{iarg},'data')
+            cdata=inputname(iarg+1);
+        end
+    else
+        error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']);
+    end
+    iarg=iarg+2;
+end
+
+if exist('data','var') && ~isempty(data)
+    if     (numel(data)==md.numberofelements)
+        edata=data;
+    elseif (numel(data)==md.numberofgrids)
+        ndata=data;
+        display('Averaging nodal data to element data.');
+        edata=zeros(1,md.numberofelements);
+        for i=1:size(md.elements,1)
+            for j=1:size(md.elements,2)
+                edata(i)=edata(i)+ndata(md.elements(i,j));
+            end
+            edata(i)=edata(i)/size(md.elements,2);
+        end
+    else
+        error(['Data has incorrect number of ' num2str(numel(data)) ' values.']);
+    end
+end
+
+%  colormap command operates on a figure, so create an invisible one
+%  (could also directly call colormaps, e.g. jet(64), but risky)
+
+hfig=figure('Visible','off');
+if exist('cmap','var')
+    colormap(cmap)
+end
+cmap=colormap;
+close(hfig)
+    
+if exist('edata','var')
+    if ~exist('cmin','var')
+        cmin=min(min(edata));
+    end
+    if ~exist('cmax','var')
+        cmax=max(max(edata));
+    end
+end
+
+if ~exist('alt','var')
+    alt=10000;
+end
+
+%  write folder for partition elements
 
 if (~exist('prtplt','var') || strncmpi(prtplt,'on' ,2) || strncmpi(prtplt,'y',1)) && ...
     md.npart
     fprintf(fid,'    <Folder>\n');
-    fprintf(fid,'      <name>Partitions</name>\n');
+    fprintf(fid,'      <name>Partition Elements</name>\n');
     fprintf(fid,'      <visibility>1</visibility>\n');
     fprintf(fid,'      <description>Partitions=%d, Nodes=%d</description>\n',...
@@ -415 +805 @@
 end
 
-%  write folder for partition edges
+end
+
+%%
+%  create kml linestrings for the partition edges.
+%
+%  []=kml_part_edges(fid,md,params)
+%
+%  where the required input is:
+%    fid           (char, file ID of .kml file)
+%    md            (model, model class object)
+%
+%  the optional input is:
+%    params        (string/numeric, parameter names and values)
+%
+%  and the optional input is:
+%    data          (numeric, element or nodal results data)
+%    alt           (numeric, altitude for polygons, default 10000)
+%    cmin          (numeric, minimum of color map)
+%    cmax          (numeric, maximum of color map)
+%    cmap          (char or numeric, colormap definition)
+%    prtplt        (char, 'off'/'no' for partition segment plot)
+%
+function []=kml_part_edges(varargin)
+
+if ~nargin
+    help kml_part_edges
+    return
+end
+
+%%  process input data
+
+iarg=1;
+if (nargin >= 1)
+    fid=varargin{1};
+end
+if ~exist('fid','var') || isempty(fid) || (fid < 0)
+    error('File ID ''%d'' must be open.',fid);
+end
+
+iarg=iarg+1;
+if (nargin >= 2)
+    md=varargin{2};
+end
+if ~exist('md','var') || isempty(md) || ~isa(md,'model')
+    error(['Model ''' inputname(iarg) ''' is unrecognized class ''' class(md) '''.']);
+end
+
+%  parameters
+
+iarg=iarg+1;
+while (iarg <= nargin-1)
+    if ischar(varargin{iarg})
+        eval([varargin{iarg} '=varargin{iarg+1};']);
+        if (numel(varargin{iarg+1}) <= 20)
+            disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']);
+        else
+            disp([varargin{iarg} '=' string_size(varargin{iarg+1}) ' ' class(varargin{iarg+1}) ';']);
+        end
+        if strcmpi(varargin{iarg},'data')
+            cdata=inputname(iarg+1);
+        end
+    else
+        error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']);
+    end
+    iarg=iarg+2;
+end
+
+if exist('data','var') && ~isempty(data)
+    if     (numel(data)==md.numberofelements)
+        edata=data;
+    elseif (numel(data)==md.numberofgrids)
+        ndata=data;
+        display('Averaging nodal data to element data.');
+        edata=zeros(1,md.numberofelements);
+        for i=1:size(md.elements,1)
+            for j=1:size(md.elements,2)
+                edata(i)=edata(i)+ndata(md.elements(i,j));
+            end
+            edata(i)=edata(i)/size(md.elements,2);
+        end
+    else
+        error(['Data has incorrect number of ' num2str(numel(data)) ' values.']);
+    end
+end
+
+%  colormap command operates on a figure, so create an invisible one
+%  (could also directly call colormaps, e.g. jet(64), but risky)
+
+hfig=figure('Visible','off');
+if exist('cmap','var')
+    colormap(cmap)
+end
+cmap=colormap;
+close(hfig)
+    
+if exist('edata','var')
+    if ~exist('cmin','var')
+        cmin=min(min(edata));
+    end
+    if ~exist('cmax','var')
+        cmax=max(max(edata));
+    end
+end
+
+if ~exist('alt','var')
+    alt=10000;
+end
+
+%%  write folder for partition edges
 
 if (~exist('prtplt','var') || strncmpi(prtplt,'on' ,2) || strncmpi(prtplt,'y',1)) && ...
@@ -510 +1008 @@
             j=iloop(i);
             while (j < iloop(i+1))
-%  find which edge of element is referenced in perimeter list
+%  find which side of element is referenced in perimeter list
                 for l=1:size(elemp,2)
                     if ((elemp(elemper(j),l)          == edgeper(j,1)) && ...
@@ -521 +1019 @@
                 end
 
-%  check if element edge connects nodes in partition
+%  check if element side connects nodes in partition
                 if (epartp(elemper(j),jedge)          == k) && ...
                    (epartp(elemper(j),mod(jedge,3)+1) == k)
-%  write out specified element edge
+%  write out specified element side
                     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.'])
+%  if first edge, write out first node
                     if ~elast
                         [lat(end+1),long(end+1)]=mapxy(md.x(edgeper(j,1)),md.y(edgeper(j,1)),'s');
@@ -546 +1045 @@
 
                     while 1
-%  if first edge, figure out direction from perimeter direction
+%  if first edge, figure out direction from perimeter edge direction
                         if ~nlast && ~slast
                             nlast=find(elemp(ielem,:)==edgeper(j,1));
@@ -566 +1065 @@
 
 %  find how many nodes of current element are in current partition
-%  (1 or 2, not 3)
+%  (1 or 2, not 3) and handle each permutation separately
                             ipart=find(epartp(ielem,:)==k);
 %  two nodes are in current partition, so cut off other node
@@ -597 +1096 @@
                                             nlast=0;
                                     end
+%  other two are in the same partition, so cut them off
                                 else
-%  other two are in the same partition, so cut them off
                                     switch ipart
                                         case nlast
@@ -616 +1115 @@
 %  last edge exited last element at node
                             if nlast
-%  write out first node of first edge for half-edge to midpoint
+%  write out first node of first side for half-edge to midpoint
                                 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.'])
                                 [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,nlast))),...
@@ -624 +1123 @@
                             nlast=0;
                             
-                            if ~elast
-%  write out midpoint of first edge
-                                [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,slast))...
-                                                               +md.x(elemp(ielem,mod(slast,3)+1)))/2.,...
-                                                               (md.y(elemp(ielem,slast))...
-                                                               +md.y(elemp(ielem,mod(slast,3)+1)))/2.,'s');
-                                fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
-                            end
+%  write out midpoint of first side
+                            [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,slast))...
+                                                           +md.x(elemp(ielem,mod(slast,3)+1)))/2.,...
+                                                           (md.y(elemp(ielem,slast))...
+                                                           +md.y(elemp(ielem,mod(slast,3)+1)))/2.,'s');
+                            fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
                         end
 
@@ -640 +1137 @@
                                 nlast=find(elemp(ielem,:)==nlast,1);
                             end
-%  half-edge occurs on unshared edge from current node (unique unless mesh
+%  half-edge occurs on unshared side from current node (unique unless mesh
 %  is only attached at node)
                             switch nlast
@@ -668 +1165 @@
                                     end
                             end
-%  write out half-edge from current node to midpoint of unshared edge
+%  write out half-edge from current node to midpoint of unshared side
                             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.'])
                             [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,nlast))...
@@ -680 +1177 @@
                         elseif slast
                             if elast
-%  find where last edge on previous element occurs on current element
+%  find where last side on previous element occurs on current element
                                 slast=find(edgeadjp(ielem,:)==elast,1);
                             end
@@ -686 +1183 @@
 
 %  find how many nodes of current element are in current partition
-%  (1 or 2, not 3)
+%  (1 or 2, not 3) and handle each permutation separately
                         ipart=find(epartp(ielem,:)==k);
                         if (length(ipart) == 2)
@@ -702 +1199 @@
                                (epartp(ielem,2) ~= epartp(ielem,3)) && ...
                                (epartp(ielem,3) ~= epartp(ielem,1))
-%  all different, so write out centroid
+%  all different, so cut through centroid
                                 disp(['element ielem=' int2str(ielem) ' centroid written.'])
                                 [lat(end+1),long(end+1)]=mapxy(sum(md.x(elemp(ielem,:)))/3.,...
@@ -718 +1215 @@
                             end
                         end
-%  write out midpoint of opposite edge
+%  write out midpoint of opposite side
                         disp(['segment j=' int2str(j) ' internal edge from side ' int2str(slast) ' to side ' int2str(snext) ' from element ' int2str(ielem) ' written.'])
                         [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,snext))...
@@ -728 +1225 @@
                         nlast=0;
                         slast=snext;
-%  find adjacent element to opposite edge
+%  find adjacent element to opposite side
                         ielem=edgeadjp(elast,slast);
-%  if opposite edge is unshared, find it in edge perimeter list
+%  if opposite side is unshared, find it in edge perimeter list
                         if ~ielem
                             jlast=find(elemper(j:end)==elast)+j-1;
@@ -746 +1243 @@
                                 j=iloop(i+1)-1;
                             end
-%  write out half-edge from midpoint of unshared edge to node
+%  write out half-edge from midpoint of unshared side to node
                             if (epartp(elast,slast) == k)
                                 nnext=slast;
@@ -785 +1282 @@
 end
 
-%  write trailer data
-
-fprintf(fid,'  </Document>\n');
-fprintf(fid,'</kml>\n');
-
-fclose(fid);
-display('End of file successfully written.');
-
+end
+
+%%
+%  create kml polygons for the partitions.
+%
+%  []=kml_partitions(fid,md,params)
+%
+%  where the required input is:
+%    fid           (char, file ID of .kml file)
+%    md            (model, model class object)
+%
+%  the optional input is:
+%    params        (string/numeric, parameter names and values)
+%
+%  and the optional input is:
+%    data          (numeric, element or nodal results data)
+%    alt           (numeric, altitude for polygons, default 10000)
+%    cmin          (numeric, minimum of color map)
+%    cmax          (numeric, maximum of color map)
+%    cmap          (char or numeric, colormap definition)
+%    prtplt        (char, 'off'/'no' for partition segment plot)
+%
+function []=kml_partitions(varargin)
+
+if ~nargin
+    help kml_part_edges
+    return
+end
+
+%%  process input data
+
+iarg=1;
+if (nargin >= 1)
+    fid=varargin{1};
+end
+if ~exist('fid','var') || isempty(fid) || (fid < 0)
+    error('File ID ''%d'' must be open.',fid);
+end
+
+iarg=iarg+1;
+if (nargin >= 2)
+    md=varargin{2};
+end
+if ~exist('md','var') || isempty(md) || ~isa(md,'model')
+    error(['Model ''' inputname(iarg) ''' is unrecognized class ''' class(md) '''.']);
+end
+
+%  parameters
+
+iarg=iarg+1;
+while (iarg <= nargin-1)
+    if ischar(varargin{iarg})
+        eval([varargin{iarg} '=varargin{iarg+1};']);
+        if (numel(varargin{iarg+1}) <= 20)
+            disp([varargin{iarg} '=' any2str(varargin{iarg+1}) ';']);
+        else
+            disp([varargin{iarg} '=' string_size(varargin{iarg+1}) ' ' class(varargin{iarg+1}) ';']);
+        end
+        if strcmpi(varargin{iarg},'data')
+            cdata=inputname(iarg+1);
+        end
+    else
+        error(['''' any2str(varargin{iarg}) ''' is not a parameter name.']);
+    end
+    iarg=iarg+2;
+end
+
+if exist('data','var') && ~isempty(data)
+    if     (numel(data)==md.numberofelements)
+        edata=data;
+    elseif (numel(data)==md.numberofgrids)
+        ndata=data;
+        display('Averaging nodal data to element data.');
+        edata=zeros(1,md.numberofelements);
+        for i=1:size(md.elements,1)
+            for j=1:size(md.elements,2)
+                edata(i)=edata(i)+ndata(md.elements(i,j));
+            end
+            edata(i)=edata(i)/size(md.elements,2);
+        end
+    else
+        error(['Data has incorrect number of ' num2str(numel(data)) ' values.']);
+    end
+end
+
+%  colormap command operates on a figure, so create an invisible one
+%  (could also directly call colormaps, e.g. jet(64), but risky)
+
+hfig=figure('Visible','off');
+if exist('cmap','var')
+    colormap(cmap)
+end
+cmap=colormap;
+close(hfig)
+    
+if exist('edata','var')
+    if ~exist('cmin','var')
+        cmin=min(min(edata));
+    end
+    if ~exist('cmax','var')
+        cmax=max(max(edata));
+    end
+end
+
+if ~exist('alt','var')
+    alt=10000;
+end
+
+%%  write folder for partitions
+
+if (~exist('prtplt','var') || strncmpi(prtplt,'on' ,2) || strncmpi(prtplt,'y',1)) && ...
+    md.npart
+    fprintf(fid,'    <Folder>\n');
+    fprintf(fid,'      <name>Partitions</name>\n');
+    fprintf(fid,'      <visibility>1</visibility>\n');
+    fprintf(fid,'      <description>Partitions=%d, Nodes=%d</description>\n',...
+        md.npart,md.numberofgrids);
+
+%  write each partition loop as linestrings
+
+    disp(['Writing ' num2str(md.npart) ' partitions as KML polygons.']);
+    epart=md.part(md.elements)+1;
+    if exist('ndata','var') || exist('edata','var')
+        pdata=zeros(1,md.npart);
+        pdata(:)=NaN;
+    end
+
+%  loop over each partition
+
+    for k=1:md.npart
+        disp(['partition k=' int2str(k)])
+        
+%  for each partition, find all the included elements and determine the
+%  perimeter (including those shared by another partition)
+
+        [icol,irow]=find(epart'==k);
+        irow=unique(irow);
+        elemp=md.elements(irow,:);
+        epartp=epart(irow,:);
+        nodeconp=nodeconnectivity(elemp,md.numberofgrids);
+        [edgeadjp]=edgeadjacency(elemp,nodeconp);
+        [edgeper,elemper,iloop]=edgeperimeter(elemp,nodeconp,edgeadjp);
+        iloop(end+1)=size(edgeper,1)+1;
+        
+%  determine the data to be used for the colors (if any)
+
+        if exist('ndata','var')
+            pdata(k)=ndata(find(md.part+1==k,1));
+        elseif exist('edata','var')
+            for i=1:size(epartp,1)
+                if isempty(find(epart(i,:)~=k,1))
+                    pdata(k)=edata(irow(i));
+                    break
+                end
+            end
+            if isnan(pdata(k))
+                warning('Data for Partition %d is not defined.\n',k)
+            end
+        end
+        
+%  loop over each loop of the perimeter for the given partition
+
+        for i=1:length(iloop)-1
+            fprintf(fid,'      <Placemark>\n');
+            if (length(iloop)-1 > 1)
+                fprintf(fid,'        <name>Partition %d, Loop %d</name>\n',k,i);
+            else
+                fprintf(fid,'        <name>Partition %d</name>\n',k);
+            end
+            fprintf(fid,'        <visibility>1</visibility>\n');
+            if exist('pdata','var')
+                fprintf(fid,'        <description>Partition data: %g</description>\n',pdata(k));
+                imap = fix((pdata(k)-cmin)/(cmax-cmin)*size(cmap,1))+1;
+                if     (imap >= 1) && (imap <= size(cmap,1))
+                    fprintf(fid,'        <styleUrl>#MatlabColor%d</styleUrl>\n',imap);
+                elseif (pdata(k) == cmax)
+                    fprintf(fid,'        <styleUrl>#MatlabColor%d</styleUrl>\n',size(cmap,1));
+                else
+                    fprintf(fid,'        <styleUrl>#BlackLineEmptyPoly</styleUrl>\n');
+                end
+            else
+                fprintf(fid,'        <styleUrl>#BlackLineRandomPoly</styleUrl>\n');
+            end
+            fprintf(fid,'        <Polygon>\n');
+            fprintf(fid,'          <extrude>1</extrude>\n');
+            fprintf(fid,'          <altitudeMode>relativeToGround</altitudeMode>\n');
+            fprintf(fid,'          <outerBoundaryIs>\n');
+            fprintf(fid,'            <LinearRing>\n');
+            fprintf(fid,'              <coordinates>\n');
+            elast=0;
+            nlast=0;
+            slast=0;
+            lat=[];
+            long=[];
+            
+%  loop over the element edges on the loop of the partition
+
+            j=iloop(i);
+            while (j < iloop(i+1))
+%  find which side of element is referenced in perimeter list
+                for l=1:size(elemp,2)
+                    if ((elemp(elemper(j),l)          == edgeper(j,1)) && ...
+                        (elemp(elemper(j),mod(l,3)+1) == edgeper(j,2))) || ...
+                       ((elemp(elemper(j),l)          == edgeper(j,2)) && ...
+                        (elemp(elemper(j),mod(l,3)+1) == edgeper(j,1)))
+                        jedge=l;
+                        break
+                    end
+                end
+
+%  check if element side connects nodes in partition
+                if (epartp(elemper(j),jedge)          == k) && ...
+                   (epartp(elemper(j),mod(jedge,3)+1) == k)
+%  write out specified element side
+                    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.'])
+%  if first edge, write out first node
+                    if ~elast
+                        [lat(end+1),long(end+1)]=mapxy(md.x(edgeper(j,1)),md.y(edgeper(j,1)),'s');
+                        fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
+                    end
+                    [lat(end+1),long(end+1)]=mapxy(md.x(edgeper(j,2)),md.y(edgeper(j,2)),'s');
+                    fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
+                    elast=elemper(j);
+                    nlast=edgeper(j,2);
+                    slast=0;
+                    j=j+1;
+                    
+%  element not entirely within partition, so figure out boundary
+                else
+                    disp(['segment j=' int2str(j) ' from element ' int2str(elemper(j)) ' shared by other partitions.'])
+                    ielem=elemper(j);
+
+%  follow partition boundary through elements not wholly in partition
+%  (may include elements not in perimeter list)
+
+                    while 1
+%  if first edge, figure out direction from perimeter edge direction
+                        if ~nlast && ~slast
+                            nlast=find(elemp(ielem,:)==edgeper(j,1));
+                            nnext=find(elemp(ielem,:)==edgeper(j,2));
+                            if     (nlast+nnext == 3)
+                                slast=1;
+                            elseif (nlast+nnext == 5)
+                                slast=2;
+                            elseif (nlast+nnext == 4)
+                                slast=3;
+                            end
+                            if     (nnext+(6-nlast-nnext) == 3)
+                                snext=1;
+                            elseif (nnext+(6-nlast-nnext) == 5)
+                                snext=2;
+                            elseif (nnext+(6-nlast-nnext) == 4)
+                                snext=3;
+                            end
+
+%  find how many nodes of current element are in current partition
+%  (1 or 2, not 3) and handle each permutation separately
+                            ipart=find(epartp(ielem,:)==k);
+%  two nodes are in current partition, so cut off other node
+                            if (length(ipart) == 2)
+                                switch 6-sum(ipart)
+                                    case nlast
+                                        slast=6-slast-snext;
+                                        nlast=0;
+                                    case nnext
+                                        if (epartp(ielem,nnext) == k)
+                                            nlast=nnext;
+                                        end
+                                    otherwise
+                                        slast=6-slast-snext;
+                                        nlast=0;
+                                end
+%  one node is in current partition
+                            else
+%  all different, so cut through centroid
+                                if (epartp(ielem,1) ~= epartp(ielem,2)) && ...
+                                   (epartp(ielem,2) ~= epartp(ielem,3)) && ...
+                                   (epartp(ielem,3) ~= epartp(ielem,1))
+                                    switch ipart
+                                        case {nlast,nnext}
+                                            if (epartp(ielem,nnext) == k)
+                                                nlast=nnext;
+                                            end
+                                        otherwise
+                                            slast=6-slast-snext;
+                                            nlast=0;
+                                    end
+%  other two are in the same partition, so cut them off
+                                else
+                                    switch ipart
+                                        case nlast
+                                            if (epartp(ielem,nnext) == k)
+                                                nlast=nnext;
+                                            end
+                                        case nnext
+                                            slast=snext;
+                                            nlast=0;
+                                        otherwise
+                                            slast=6-slast-snext;
+                                            nlast=0;
+                                    end
+                                end
+                            end
+
+%  last edge exited last element at node
+                            if nlast
+%  write out first node of first side for half-edge to midpoint
+                                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.'])
+                                [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,nlast))),...
+                                                               (md.y(elemp(ielem,nlast))),'s');
+                                fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
+                            end
+                            nlast=0;
+                            
+%  write out midpoint of first side
+                            [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,slast))...
+                                                           +md.x(elemp(ielem,mod(slast,3)+1)))/2.,...
+                                                           (md.y(elemp(ielem,slast))...
+                                                           +md.y(elemp(ielem,mod(slast,3)+1)))/2.,'s');
+                            fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
+                        end
+
+%  last edge exited last element at node
+                        if nlast
+                            if elast
+%  find where last node on previous element occurs on current element
+                                nlast=find(elemp(ielem,:)==nlast,1);
+                            end
+%  half-edge occurs on unshared side from current node (unique unless mesh
+%  is only attached at node)
+                            switch nlast
+                                case 1
+                                    if ~edgeadjp(ielem,1)
+                                        nnext=2;
+                                        slast=1;
+                                    else
+                                        nnext=3;
+                                        slast=3;
+                                    end
+                                case 2
+                                    if ~edgeadjp(ielem,2)
+                                        nnext=3;
+                                        slast=2;
+                                    else
+                                        nnext=1;
+                                        slast=1;
+                                    end
+                                case 3
+                                    if ~edgeadjp(ielem,3)
+                                        nnext=1;
+                                        slast=3;
+                                    else
+                                        nnext=2;
+                                        slast=2;
+                                    end
+                            end
+%  write out half-edge from current node to midpoint of unshared side
+                            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.'])
+                            [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,nlast))...
+                                                           +md.x(elemp(ielem,nnext)))/2.,...
+                                                           (md.y(elemp(ielem,nlast))...
+                                                           +md.y(elemp(ielem,nnext)))/2.,'s');
+                            fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
+                            nlast=0;
+
+%  last edge exited last element at midpoint of side
+                        elseif slast
+                            if elast
+%  find where last side on previous element occurs on current element
+                                slast=find(edgeadjp(ielem,:)==elast,1);
+                            end
+                        end
+
+%  find how many nodes of current element are in current partition
+%  (1 or 2, not 3) and handle each permutation separately
+                        ipart=find(epartp(ielem,:)==k);
+                        if (length(ipart) == 2)
+%  two nodes are in current partition, so cut off other node
+                            switch 6-sum(ipart)
+                                case 1
+                                    snext=3+1-slast;
+                                case 2
+                                    snext=1+2-slast;
+                                case 3
+                                    snext=2+3-slast;
+                            end
+                        else
+                            if (epartp(ielem,1) ~= epartp(ielem,2)) && ...
+                               (epartp(ielem,2) ~= epartp(ielem,3)) && ...
+                               (epartp(ielem,3) ~= epartp(ielem,1))
+%  all different, so cut through centroid
+                                disp(['element ielem=' int2str(ielem) ' centroid written.'])
+                                [lat(end+1),long(end+1)]=mapxy(sum(md.x(elemp(ielem,:)))/3.,...
+                                                               sum(md.y(elemp(ielem,:)))/3.,'s');
+                                fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
+                            end
+%  one node is in current partition, so cut off other two nodes
+                            switch ipart
+                                case 1
+                                    snext=3+1-slast;
+                                case 2
+                                    snext=1+2-slast;
+                                case 3
+                                    snext=2+3-slast;
+                            end
+                        end
+%  write out midpoint of opposite side
+                        disp(['segment j=' int2str(j) ' internal edge from side ' int2str(slast) ' to side ' int2str(snext) ' from element ' int2str(ielem) ' written.'])
+                        [lat(end+1),long(end+1)]=mapxy((md.x(elemp(ielem,snext))...
+                                                       +md.x(elemp(ielem,mod(snext,3)+1)))/2.,...
+                                                       (md.y(elemp(ielem,snext))...
+                                                       +md.y(elemp(ielem,mod(snext,3)+1)))/2.,'s');
+                        fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
+                        elast=ielem;
+                        nlast=0;
+                        slast=snext;
+%  find adjacent element to opposite side
+                        ielem=edgeadjp(elast,slast);
+%  if opposite side is unshared, find it in edge perimeter list
+                        if ~ielem
+                            jlast=find(elemper(j:end)==elast)+j-1;
+                            j=0;
+                            for l=1:length(jlast)
+                                if ((elemp(elast,slast)          == edgeper(jlast(l),1)) && ...
+                                    (elemp(elast,mod(slast,3)+1) == edgeper(jlast(l),2))) || ...
+                                   ((elemp(elast,slast)          == edgeper(jlast(l),2)) && ...
+                                    (elemp(elast,mod(slast,3)+1) == edgeper(jlast(l),1)))
+                                    j=jlast(l);
+                                    break
+                                end
+                            end
+                            if ~j
+                                j=iloop(i+1)-1;
+                            end
+%  write out half-edge from midpoint of unshared side to node
+                            if (epartp(elast,slast) == k)
+                                nnext=slast;
+                            else
+                                nnext=mod(slast,3)+1;
+                            end
+                            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.'])
+                            [lat(end+1),long(end+1)]=mapxy(md.x(elemp(elast,nnext)),...
+                                                           md.y(elemp(elast,nnext)),'s');
+                            fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(end),lat(end),alt);
+                            break
+%  if not unshared, advance perimeter list and watch for end
+                        else
+                            if (elast == elemper(j))
+                                if (j+1 < iloop(i+1)) && ...
+                                   ~isempty(find(elemper(j+1:end)~=elast,1))
+                                    j=j+find(elemper(j+1:end)~=elast,1);
+                                else
+                                    break
+                                end
+                            end
+                        end
+                    end
+                    j=j+1;
+                end
+            end
+%            fprintf(fid,'                %0.16g,%0.16g,%0.16g\n',long(iloop(1)),lat(iloop(1)),alt);
+
+            fprintf(fid,'              </coordinates>\n');
+            fprintf(fid,'            </LinearRing>\n');
+            fprintf(fid,'          </outerBoundaryIs>\n');
+            fprintf(fid,'        </Polygon>\n');
+            fprintf(fid,'      </Placemark>\n');
+        end
+    end
+    fprintf(fid,'    </Folder>\n');
+end
+
+end
+