Index: /issm/trunk/src/m/kml/kml_mesh_write.m
===================================================================
--- /issm/trunk/src/m/kml/kml_mesh_write.m (revision 6414)
+++ /issm/trunk/src/m/kml/kml_mesh_write.m (revision 6415)
@@ -93,4 +93,17 @@
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,4 +206,158 @@
% 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,' \n');
+fprintf(fid,'\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,' \n');
if exist('cdata','var') && ~isempty(cdata)
@@ -204,17 +371,4 @@
% 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,5 +410,69 @@
fprintf(fid,' \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,5 +510,69 @@
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,10 +616,118 @@
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,' \n');
- fprintf(fid,' Partitions\n');
+ fprintf(fid,' Partition Elements\n');
fprintf(fid,' 1\n');
fprintf(fid,' Partitions=%d, Nodes=%d\n',...
@@ -415,5 +805,113 @@
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,5 +1008,5 @@
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,9 +1019,10 @@
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,5 +1045,5 @@
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,5 +1065,5 @@
% 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,6 +1096,6 @@
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,5 +1115,5 @@
% 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,12 +1123,10 @@
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,5 +1137,5 @@
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,5 +1165,5 @@
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,5 +1177,5 @@
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,5 +1183,5 @@
% 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,5 +1199,5 @@
(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,5 +1215,5 @@
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,7 +1225,7 @@
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,5 +1243,5 @@
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,10 +1282,476 @@
end
-% write trailer data
-
-fprintf(fid,' \n');
-fprintf(fid,'\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,' \n');
+ fprintf(fid,' Partitions\n');
+ fprintf(fid,' 1\n');
+ fprintf(fid,' Partitions=%d, Nodes=%d\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,' \n');
+ if (length(iloop)-1 > 1)
+ fprintf(fid,' Partition %d, Loop %d\n',k,i);
+ else
+ fprintf(fid,' Partition %d\n',k);
+ end
+ fprintf(fid,' 1\n');
+ if exist('pdata','var')
+ fprintf(fid,' Partition data: %g\n',pdata(k));
+ imap = fix((pdata(k)-cmin)/(cmax-cmin)*size(cmap,1))+1;
+ if (imap >= 1) && (imap <= size(cmap,1))
+ fprintf(fid,' #MatlabColor%d\n',imap);
+ elseif (pdata(k) == cmax)
+ fprintf(fid,' #MatlabColor%d\n',size(cmap,1));
+ else
+ fprintf(fid,' #BlackLineEmptyPoly\n');
+ end
+ else
+ fprintf(fid,' #BlackLineRandomPoly\n');
+ end
+ fprintf(fid,' \n');
+ fprintf(fid,' 1\n');
+ fprintf(fid,' relativeToGround\n');
+ fprintf(fid,' \n');
+ fprintf(fid,' \n');
+ fprintf(fid,' \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,' \n');
+ fprintf(fid,' \n');
+ fprintf(fid,' \n');
+ fprintf(fid,' \n');
+ fprintf(fid,' \n');
+ end
+ end
+ fprintf(fid,' \n');
+end
+
+end
+