Index: /issm/trunk/src/m/partition/adjacency_matrix.m
===================================================================
--- /issm/trunk/src/m/partition/adjacency_matrix.m	(revision 2915)
+++ /issm/trunk/src/m/partition/adjacency_matrix.m	(revision 2915)
@@ -0,0 +1,201 @@
+%
+%  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));
+
+%%  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
+
+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
+end
+
+disp(sprintf('Total area=%f; min area=%f, max area=%f, ratio=%f.',...
+             asum,amin,amax,amax/amin));
+disp(sprintf('Total volume=%f; min volume=%f, max volume=%f, ratio=%f.',...
+             vsum,vmin,vmax,vmax/vmin));
+disp(sprintf('Total weight=%f; min weight=%f, max weight=%f, ratio=%f.',...
+             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
+
+disp('Creating adjacency matrix.');
+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/partition/eric_100126.m.bak
===================================================================
--- /issm/trunk/src/m/partition/eric_100126.m.bak	(revision 2915)
+++ /issm/trunk/src/m/partition/eric_100126.m.bak	(revision 2915)
@@ -0,0 +1,15 @@
+addpath '/u/wilkes-r1b/jschierm/Libs/meshpart'
+addpath '/u/wilkes-r1b/jschierm/Libs/meshpart/chaco'
+addpath /u/wilkes-r1b/jschierm/Libs/scotch_5.1
+addpath /home/jschierm/Libs/scotch_5.1/bin
+load Pig_mesh1M  % just for me -- create your own mesh
+%  create adjacency matrix, vertex list, and vertex weights
+[adj_mat,vlist,vwgt]=adjacency_matrix(md.elements,[md.x md.y md.z]);
+%  partition into 100 parts, but ignore vertex and edge weights
+[status,maptab]=gmap(adj_mat,vlist,[],[],'cmplt',[100],'-vm','-vs','-vt');
+%  plot partitions
+figure
+plotmodel(md,'data',maptab(:,2))
+%  plot histogram of number of nodes in each partition
+figure
+hist(maptab(:,2),min(maptab(:,2)):1:max(maptab(:,2)))
Index: /issm/trunk/src/m/partition/gmap.m
===================================================================
--- /issm/trunk/src/m/partition/gmap.m	(revision 2915)
+++ /issm/trunk/src/m/partition/gmap.m	(revision 2915)
@@ -0,0 +1,64 @@
+%
+%  function to call the gmap module of the scotch partitioner.
+%
+%  [status,maptab]=gmap(adj_mat,vlist,vwgt,ewgt,atype,apar,...
+%                       options)
+%
+%  where the required input is:
+%    adj_mat    (double [sparse nv x nv], vertex adjacency matrix)
+%    vlist      (double [nv], vertex labels or [])
+%    vwgt       (double [nv], vertex weights (integers) or [])
+%    ewgt       (double [sparse nv x nv], edge weights (integers) or [])
+%    atype      (character, architecture type)
+%                 'cmplt'      complete graph
+%                 'cmpltw'     weighted complete graph
+%                 'hcub'       binary hypercube
+%                 'leaf'       tree-leaf architecture
+%                 'mesh2d'     bidimensional array
+%                 'mesh3d'     tridimensional array
+%                 'torus2d'    bidimensional array with wraparound edges
+%                 'torus3d'    tridimensional array with wraparound edges
+%    apars      (double, architecture params (corresponding to atype))
+%                 [size]                     cmplt
+%                 [size load0 load1 ...]     cmpltw
+%                 [dim]                      hcub
+%                 [height cluster weight]    leaf
+%                 [dimX dimY]                mesh2d
+%                 [dimX dimY dimZ]           mesh3d
+%                 [dimX dimY]                torus2d
+%                 [dimX dimY dimZ]           torus3d
+%
+%  the required output is:
+%    status     (double, return code from gmap)
+%    maptab     (double [nv x 2], vertex labels and partitions)
+%
+%  the optional input is:
+%    options    (character, options to gmap)
+%               "  -h         : Display this help"
+%               "  -m<strat>  : Set mapping strategy (see user's manual)"
+%               "  -s<obj>    : Force unity weights on <obj>:"
+%               "                 e  : edges"
+%               "                 v  : vertices"
+%               "  -V         : Print program version and copyright"
+%               "  -v<verb>   : Set verbose mode to <verb>:"
+%               "                 m  : mapping information"
+%               "                 s  : strategy information"
+%               "                 t  : timing information"
+%               ""
+%               "See default strategy with option '-vs'"
+%
+function [status,maptab]=gmap(adj_mat,vlist,vwgt,ewgt,atype,apars,...
+                              varargin)
+
+if ~nargin
+    help gmap
+    return
+end
+
+%  gmap_mex uses static variables, so clear those out before every run
+clear gmap_mex
+
+[status,maptab]=gmap_mex(adj_mat,vlist,vwgt,ewgt,atype,apars,...
+                         varargin{:});
+
+end
Index: /issm/trunk/src/m/partition/part_hist.m
===================================================================
--- /issm/trunk/src/m/partition/part_hist.m	(revision 2915)
+++ /issm/trunk/src/m/partition/part_hist.m	(revision 2915)
@@ -0,0 +1,29 @@
+function []=part_hist(maptab,wgt);
+
+imin=min(maptab(:,2));
+imax=max(maptab(:,2));
+
+part=zeros(imax-imin+1,2);
+
+for i=imin:imax
+    ind=find(maptab(:,2) == 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 Elements in Each Partition')
+
+subplot(2,1,2)
+bar(imin:imax,part(:,2));
+xlim([imin imax])
+title('Weight of Elements in Each Partition')
+
+end
Index: /issm/trunk/src/m/partition/partition.m
===================================================================
--- /issm/trunk/src/m/partition/partition.m	(revision 2915)
+++ /issm/trunk/src/m/partition/partition.m	(revision 2915)
@@ -0,0 +1,13 @@
+%  create adjacency matrix, vertex list, and vertex weights
+[adj_mat,vlist,vwgt]=adjacency_matrix(md.elements,[md.x md.y md.z]);
+%  partition into 100 parts, but ignore vertex and edge weights
+%[status,maptab]=gmap(adj_mat,vlist,[],[],'cmplt',[100],'-vm','-vs','-vt');
+[status,maptab]=gmap(adj_mat,vlist,[],[],'cmplt',[1000]);
+
+%  plot partitions
+%figure
+%plotmodel(md,'data',maptab(:,2))
+%  plot histogram of number of nodes in each partition
+figure
+%hist(maptab(:,2),min(maptab(:,2)):1:max(maptab(:,2)))
+part_hist(maptab,vwgt);
Index: /issm/trunk/src/m/partition/plotedges.m
===================================================================
--- /issm/trunk/src/m/partition/plotedges.m	(revision 2915)
+++ /issm/trunk/src/m/partition/plotedges.m	(revision 2915)
@@ -0,0 +1,18 @@
+flags=zeros(md.numberofgrids,1);
+for i=1:md.numberofgrids,
+
+	if mod(i,100),
+		disp([num2str(i) 'vs' num2str(md.numberofgrids)]);
+	end
+
+
+	grids=adj_mat(i,:);
+	grids=grids(find(grids));
+	part_i=maptab(i,2);
+	part_connected=maptab(grids,2);
+
+	index=find(part_i~=part_connected);
+	grids_on_edge=grids(index);
+	flags(grids_on_edge)=1;
+
+end