Index: /issm/trunk/src/m/classes/public/partition/adjacency.jschierm.m
===================================================================
--- /issm/trunk/src/m/classes/public/partition/adjacency.jschierm.m	(revision 4660)
+++ /issm/trunk/src/m/classes/public/partition/adjacency.jschierm.m	(revision 4660)
@@ -0,0 +1,216 @@
+%
+%  function to create the adjacency matrix from the connectivity table.
+%
+%  [adj_mat,vlist,vwgt]=adjacency_matrix(elem_con,xyz)
+%
+%  where the required input is:
+%    elem_con    (double [ne x mvpe], element connectivity table)
+%    xyz         (double [nv x 3   ], vertex coordinates)
+%
+%  the required output is:
+%    adj_mat     (double [sparse nv x nv], vertex adjacency matrix)
+%    vlist       (double [nv], vertex labels)
+%    vwgt        (double [nv], vertex weights)
+%
+function [adj_mat,vlist,vwgt]=adjacency_matrix(elem_con,xyz)
+
+if ~nargin
+    help adjacency_matrix
+    return
+end
+
+%  can modify local elem_con, since it's not being returned
+
+elem_con(~isfinite(elem_con))=0;
+
+%%  create unique sorted vertex list, eliminating NaN, Inf, and zero
+
+disp('Creating unique sorted vertex list.');
+vlist=unique(nonzeros(elem_con));
+fprintf('  Found %d vertices numbered from %d to %d.\n',...
+        length(vlist),vlist(1),vlist(end));
+
+%%  create edge list and vertex weights
+
+disp('Creating edge list and vertex weights.');
+nedge=0;
+elist=zeros(numel(elem_con),2);
+asum= 0.;
+amin= Inf;
+amax=-Inf;
+vsum= 0.;
+vmin= Inf;
+vmax=-Inf;
+vwgt =zeros(size(vlist));
+
+%  loop over elements
+
+fprintf('  Processing %d elements with a maximum of %d vertices.\n',...
+        size(elem_con,1),size(elem_con,2));
+%hwbar=waitbar(0);
+for i=1:size(elem_con,1)
+    elem=nonzeros(elem_con(i,:));
+    
+%  replace vertex numbers by indices (if necessary)
+
+    if (vlist(1) ~= 1) || (vlist(end) ~= length(vlist))
+        elem=bin_search(elem,vlist);
+    end
+
+%  accumulate edge list for each element type
+
+    switch length(elem)
+        
+%  tria
+
+        case 3
+            elist(nedge+1,1)=elem(1);
+            elist(nedge+1,2)=elem(2);
+            elist(nedge+2,1)=elem(2);
+            elist(nedge+2,2)=elem(3);
+            elist(nedge+3,1)=elem(3);
+            elist(nedge+3,2)=elem(1);
+            nedge=nedge+3;
+            
+            if exist('xyz','var')
+                area=mag(cross(xyz(elem(2),:)-xyz(elem(1),:),...
+                               xyz(elem(3),:)-xyz(elem(1),:)))/2;
+                asum=asum+area;
+                if (area < amin)
+                    amin=area;
+                end
+                if (area > amax)
+                    amax=area;
+                end
+                for j=1:3
+                    vwgt(elem(j))=vwgt(elem(j))+area/3;
+                end
+            end
+            
+%  penta
+
+        case 6
+            elist(nedge+1,1)=elem(1);
+            elist(nedge+1,2)=elem(2);
+            elist(nedge+2,1)=elem(2);
+            elist(nedge+2,2)=elem(3);
+            elist(nedge+3,1)=elem(3);
+            elist(nedge+3,2)=elem(1);
+            elist(nedge+4,1)=elem(1);
+            elist(nedge+4,2)=elem(4);
+            elist(nedge+5,1)=elem(2);
+            elist(nedge+5,2)=elem(5);
+            elist(nedge+6,1)=elem(3);
+            elist(nedge+6,2)=elem(6);
+            elist(nedge+7,1)=elem(7);
+            elist(nedge+7,2)=elem(8);
+            elist(nedge+8,1)=elem(8);
+            elist(nedge+8,2)=elem(9);
+            elist(nedge+9,1)=elem(9);
+            elist(nedge+9,2)=elem(7);
+            nedge=nedge+9;
+            
+            if exist('xyz','var')
+                vol =mag(cross(xyz(elem(2),:)-xyz(elem(1),:),...
+                               xyz(elem(3),:)-xyz(elem(1),:)))/2*...
+                     mag(xyz(elem(4),:)-xyz(elem(1),:));
+                vsum=vsum+vol;
+                if (vol  < vmin)
+                    vmin=vol;
+                end
+                if (vol  > vmax)
+                    vmax=vol;
+                end
+                for j=1:6
+                    vwgt(elem(j))=vwgt(elem(j))+vol/6;
+                end
+            end
+            
+        otherwise
+            error(['Unrecognized element of length' length(elem) '.']);
+    end
+    if (i/100 == floor(i/100))
+         fprintf('  %d elements processed.\n',i);
+        %waitbar(i/size(elem_con,1),hwbar,sprintf('%d elements processed.',i));
+    end
+end
+fprintf('  %d total elements processed.\n\n',i);
+%waitbar(1,hwbar,sprintf('%d total elements processed.',i));
+%close(hwbar)
+
+if (asum > 0)
+    fprintf('Total area=%f; min area=%f, max area=%f, ratio=%f.\n',...
+            asum,amin,amax,amax/amin);
+end
+if (vsum > 0)
+    fprintf('Total volume=%f; min volume=%f, max volume=%f, ratio=%f.\n',...
+            vsum,vmin,vmax,vmax/vmin);
+end
+fprintf('Total weight=%f; min weight=%f, max weight=%f, ratio=%f.\n',...
+        sum(vwgt),min(vwgt),max(vwgt),max(vwgt)/min(vwgt));
+
+elist=elist(1:nedge,:);
+
+%%  replace vertex numbers by indices (if necessary)
+
+% if (vlist(1) ~= 1) || (vlist(end) ~= length(vlist))
+%     elist=bin_search(elist,vlist);
+% end
+
+%%  create adjacency matrix and make symmetric
+
+adj_mat=sparse(elist(:,1),elist(:,2),1,length(vlist),length(vlist));
+adj_mat=double(adj_mat | adj_mat');
+
+end
+
+%
+%  function to perform a recursive binary search for a matrix of values
+%  in an ordered vector (loop separately outside of recursion for
+%  efficiency purposes)
+%
+%  function [ind]=bin_search(val,vect)
+%
+function [ind]=bin_search(val,vect)
+
+ind=zeros(size(val));
+
+for i=1:numel(val)
+    ind(i)=bin_search_val(val(i),vect);
+end
+
+end
+
+%
+%  function to perform a recursive binary search in an ordered vector,
+%  returning NaN if value does not exist (more efficient than find or
+%  ismember, which must use linear searches and/or sort)
+%
+%  function [ind]=bin_search_val(val,vect)
+%
+function [ind]=bin_search_val(val,vect)
+
+imid=floor((1+length(vect))/2);
+
+if (val == vect(imid))
+    ind=imid;
+elseif (val < vect(imid)) && (imid > 1)
+    ind=     bin_search(val,vect(1:imid-1));
+elseif (val > vect(imid)) && (imid < length(vect))
+    ind=imid+bin_search(val,vect(imid+1:length(vect)));
+else
+    ind=NaN;
+end
+
+end
+
+%
+%  function to calculate the magnitude of a vector
+%
+%  function [vmag]=mag(vect)
+%
+function [vmag]=mag(vect)
+
+vmag=sqrt(dot(vect,vect));
+
+end
Index: /issm/trunk/src/m/classes/public/partition/adjacency.m
===================================================================
--- /issm/trunk/src/m/classes/public/partition/adjacency.m	(revision 4660)
+++ /issm/trunk/src/m/classes/public/partition/adjacency.m	(revision 4660)
@@ -0,0 +1,24 @@
+function md=adjacency(md);
+%adjacency:  compute adjacency matrix, list of vertices and list of weights.
+%
+%  function to create the adjacency matrix from the connectivity table.
+%
+%  the required output is:
+%    md.adj_mat     (double [sparse nv x nv], vertex adjacency matrix)
+%    md.vwgt        (double [nv], vertex weights)
+
+
+
+indi=[md.elements(:,1);md.elements(:,2);md.elements(:,3)];
+indj=[md.elements(:,2);md.elements(:,3);md.elements(:,1)];
+values=1;
+
+md.adjacency=sparse(indi,indj,values,md.numberofgrids,md.numberofgrids);
+md.adjacency=double([md.adjacency | md.adjacency']);
+
+%now, build vwgt:
+areas=GetAreas(md.elements,md.x,md.y);
+connectivity=md.nodeconnectivity(:,1:end-1);
+pos=find(connectivity);
+connectivity(pos)=areas(connectivity(pos))/3;
+md.vwgt=sum(connectivity,2);
Index: /issm/trunk/src/m/classes/public/partition/flagedges.m
===================================================================
--- /issm/trunk/src/m/classes/public/partition/flagedges.m	(revision 4660)
+++ /issm/trunk/src/m/classes/public/partition/flagedges.m	(revision 4660)
@@ -0,0 +1,46 @@
+function [xsegments ysegments]=flagedges(elements,x,y,partition);
+%FLAGEDGES return pairs of x,y segments, delimiting partitions.
+%
+%         Usage: [xsegments ysegments]=flagedges(elements,x,y,partition)
+%
+
+xsegments=[];
+ysegments=[];
+
+for i=1:size(elements,1),
+	m1=partition(elements(i,1));
+	m2=partition(elements(i,2));
+	m3=partition(elements(i,3));
+	x1=x(elements(i,1));
+	x2=x(elements(i,2));
+	x3=x(elements(i,3));
+	y1=y(elements(i,1));
+	y2=y(elements(i,2));
+	y3=y(elements(i,3));
+
+	if (m1~=m2) & (m1~=m3) & (m2~=m3),
+		xmiddle=(x1+x2+x3)/3;
+		ymiddle=(y1+y2+y3)/3;
+		xsegments=[xsegments; (x1+x2)/2 xmiddle];
+		xsegments=[xsegments; (x1+x3)/2 xmiddle];
+		xsegments=[xsegments; (x2+x3)/2 xmiddle];
+		ysegments=[ysegments; (y1+y2)/2 ymiddle];
+		ysegments=[ysegments; (y1+y3)/2 ymiddle];
+		ysegments=[ysegments; (y2+y3)/2 ymiddle];
+	end
+
+	if (m1==m2) & (m1~=m3),
+		xsegments=[xsegments; (x1+x3)/2 (x2+x3)/2];
+		ysegments=[ysegments; (y1+y3)/2 (y2+y3)/2];
+	end
+	if (m1==m3) & (m2~=m3),
+		xsegments=[xsegments; (x1+x2)/2 (x2+x3)/2];
+		ysegments=[ysegments; (y1+y2)/2 (y2+y3)/2];
+	end
+
+	if (m2==m3) & (m1~=m3),
+		xsegments=[xsegments; (x1+x2)/2 (x1+x3)/2];
+		ysegments=[ysegments; (y1+y2)/2 (y1+y3)/2];
+	end
+end
+
Index: /issm/trunk/src/m/classes/public/partition/part_hist.m
===================================================================
--- /issm/trunk/src/m/classes/public/partition/part_hist.m	(revision 4660)
+++ /issm/trunk/src/m/classes/public/partition/part_hist.m	(revision 4660)
@@ -0,0 +1,35 @@
+function []=part_hist(maptab,wgt);
+
+if (size(maptab,1) > size(maptab,2))
+    map=maptab(:,end);
+else
+    map=maptab(end,:)';
+end
+
+imin=min(map);
+imax=max(map);
+
+part=zeros(imax-imin+1,2);
+
+for i=imin:imax
+    ind=find(map == i);
+    part(i-imin+1,1)=length(ind);
+    if exist('wgt','var')
+        part(i-imin+1,2)=sum(wgt(ind));
+    else
+        part(i-imin+1,2)=part(i-imin+1,1);
+    end
+end
+
+figure
+subplot(2,1,1)
+bar(imin:imax,part(:,1));
+xlim([imin imax])
+title('Number of Nodes in Each Partition')
+
+subplot(2,1,2)
+bar(imin:imax,part(:,2));
+xlim([imin imax])
+title('Weight of Nodes in Each Partition')
+
+end
Index: /issm/trunk/src/m/classes/public/partition/partitioner.m
===================================================================
--- /issm/trunk/src/m/classes/public/partition/partitioner.m	(revision 4660)
+++ /issm/trunk/src/m/classes/public/partition/partitioner.m	(revision 4660)
@@ -0,0 +1,80 @@
+function part=partitioner(md,varargin)
+%PARTITIONER - partition mesh 
+%
+%   List of options to partitioner: 
+%
+%   package: 'chaco', 'metis' or 'scotch'
+%   npart: number of partitions.
+%   weighting: 'on' or 'off': default off
+%   section:  1 by defaults(1=bisection, 2=quadrisection, 3=octasection)
+%   
+%   Usage:
+%      md.part=partitioner(md,'package','chaco','npart',100,'weighting','on');
+%
+
+
+%get options: 
+options=pairoptions(varargin{:});
+
+%set defaults
+options=addfielddefault(options,'package','chaco');
+options=addfielddefault(options,'npart',10);
+options=addfielddefault(options,'weighting','on');
+options=addfielddefault(options,'section',1);
+
+%get package: 
+package=getfieldvalue(options,'package');
+npart=getfieldvalue(options,'npart');
+
+%adjacency matrix if needed:
+if (strcmpi(package,'chaco') || strcmpi(package,'scotch')),
+	md=adjacency(md);
+end
+
+
+if strcmpi(package,'chaco'),
+
+	%  default method (from chaco.m)
+	method=[1 1 0 0 1 1 50 0 .001 7654321]';
+	method(1)=3;    %  global method (3=inertial (geometric))
+	method(3)=0;    %  vertex weights (0=off, 1=on)
+	
+	%specify bisection
+	method(6)=getfieldvalue(options,'section');%  ndims (1=bisection, 2=quadrisection, 3=octasection)
+
+	%are we using weights? 
+	if strcmpi(getfieldvalue(options,'weighting'),'on'),
+		weights=floor(md.vwgt/min(md.vwgt));
+		method(3)=1;
+	else 
+		weights=[];
+	end
+	
+	%  partition into nparts
+	part=Chaco(md.adjacency,weights,[],md.x, md.y ,md.z,method,npart,[])'+1; %index partitions from 1 up. like metis.
+
+elseif strcmpi(package,'scotch'),
+
+	%are we using weights? 
+	if strcmpi(getfieldvalue(options,'weighting'),'on'),
+		weights=floor(md.vwgt/min(md.vwgt));
+	end
+	maptab=Scotch(md.adjacency,[],weights,[],'cmplt',[npart]);
+	
+	part=maptab(:,2)+1;%index partitions from 1 up. like metis.
+
+
+elseif strcmpi(package,'linear'),
+
+	part=1:1:md.numberofgrids;
+
+elseif strcmpi(package,'metis'),
+
+	[element_partitioning,part]=MeshPartition(md,npart);
+
+else
+
+	error(['partitioner error message: could not find ' package ' partitioner']);
+	help partitioner
+
+end