source: issm/trunk/src/m/classes/public/partition/adjacency.jschierm.m@ 4660

%
%  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