Index: /issm/trunk/src/m/partition/adjacency.jschierm.m =================================================================== --- /issm/trunk/src/m/partition/adjacency.jschierm.m (revision 4659) +++ /issm/trunk/src/m/partition/adjacency.jschierm.m (revision 4659) @@ -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: sm/trunk/src/m/partition/adjacency_matrix.m =================================================================== --- /issm/trunk/src/m/partition/adjacency_matrix.m (revision 4658) +++ (revision ) @@ -1,216 +1,0 @@ -% -% 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: sm/trunk/src/m/partition/part_chaco.m =================================================================== --- /issm/trunk/src/m/partition/part_chaco.m (revision 4658) +++ (revision ) @@ -1,26 +1,0 @@ -% create adjacency matrix, vertex list, and vertex weights -[adj_mat,vlist,vwgt]=adjacency_matrix(md.elements,[md.x md.y md.z]); - -% add integer scaled vertex weights into adjacency matrix for chaco (on diagonal) -adj_mat2=adj_mat+sparse(1:size(adj_mat,1),1:size(adj_mat,2),floor(vwgt/min(vwgt))); - -% 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)=1; % vertex weights (0=off, 1=on) -method(6)=1; % ndims (1=bisection, 2=quadrisection, 3=octasection) - -% partition into 100 parts -[map]=chaco(adj_mat2,[md.x md.y],method,100); -%adj_mat: don't use the weights -%adj_mat2: use the weights - -error; - -% plot partitions -figure -plotmodel(md,'data',map) - -% plot histograms of number and weight of vertices in each partition -figure -part_hist([(1:length(map))' map'],vwgt) Index: sm/trunk/src/m/partition/part_scotch.m =================================================================== --- /issm/trunk/src/m/partition/part_scotch.m (revision 4658) +++ (revision ) @@ -1,15 +1,0 @@ -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 histograms of number and weight of vertices in each partition -figure -part_hist(maptab(:,2),vwgt) Index: /issm/trunk/src/m/partition/partitioner.m =================================================================== --- /issm/trunk/src/m/partition/partitioner.m (revision 4658) +++ /issm/trunk/src/m/partition/partitioner.m (revision 4659) @@ -27,15 +27,16 @@ npart=getfieldvalue(options,'npart'); +%adjacency matrix if needed: +if (strcmpi(package,'chaco') || strcmpi(package,'scotch')), + md=adjacency(md); +end + + if strcmpi(package,'chaco'), - - % create adjacency matrix, vertex list, and vertex weights - if isempty(md.adjacency) - [md.adjacency,vlist,md.vwgt]=adjacency_matrix(md.elements,[md.x md.y md.z]); - end % 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)=1; % vertex weights (0=off, 1=on) + method(3)=0; % vertex weights (0=off, 1=on) %specify bisection @@ -44,34 +45,23 @@ %are we using weights? if strcmpi(getfieldvalue(options,'weighting'),'on'), + weights=floor(md.vwgt/min(md.vwgt)); + method(3)=1; + else + weights=[]; + end - % add integer scaled vertex weights into adjacency matrix for chaco (on diagonal) -% adj_mat2=md.adjacency+sparse(1:size(md.adjacency,1),1:size(md.adjacency,2),floor(md.vwgt/min(md.vwgt))); - - % partition into nparts - part=Chaco(md.adjacency,floor(md.vwgt/min(md.vwgt)),[],md.x, md.y ,md.z,method,npart,[])'; - else - % partition into nparts - method(3)=0; % vertex weights (0=off, 1=on) - part=Chaco(md.adjacency,[],[],md.x, md.y, md.z,method,npart,[])'; - end - part=part+1; %index partitions from 1 up. like metis. + % 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'), - % create adjacency matrix, vertex list, and vertex weights - if isempty(md.adjacency) - [md.adjacency,vlist,md.vwgt]=adjacency_matrix(md.elements,[md.x md.y md.z]); + %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. - % partition into 100 parts, but ignore vertex and edge weights - %maptab=gmap(md.adjacency,[],[],[],'cmplt',[npart],'-vm','-vs','-vt'); - - if strcmpi(getfieldvalue(options,'weighting'),'on'), - maptab=gmap(md.adjacency,[],floor(md.vwgt/min(md.vwgt)),[],'cmplt',[npart]); - else - maptab=gmap(md.adjacency,[],[],[],'cmplt',[npart]); - end - - part=maptab(:,2)+1; elseif strcmpi(package,'linear'),