Index: sm/trunk/src/m/classes/public/partition/adjacency.jschierm.m =================================================================== --- /issm/trunk/src/m/classes/public/partition/adjacency.jschierm.m (revision 4660) +++ (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