faceneighbors.m

function facenb = faceneighbors(t, opt)
%
% facenb=faceneighbors(t,opt)
%
% to find 4 face-neighboring elements of a tetrahedron
%
% author: Qianqian Fang, <q.fang at neu.edu>
%
% input:
%     t: tetrahedron element list, 4 columns of integers
%     opt: if opt='surface', return boundary triangle list
%          (should be the same as the face output from v2m)
%          if opt='rowmajor', same as 'surface', except the
%          order of the triangles are in the row-major order
%%
%          otherwise, return the element list for each element:
%          each row contains 4 numbers, representing the element
%          indices sharing triangular faces [1 2 3],[1 2 4],[1 3 4]
%          and [2 3 4] in order, where 1~4 is the node local index.
%          if the index is 0, indicating the face has no neighbor
%          (i.e. a boundary face)
%
% output:
%     facenb: see opt
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

faces = [t(:, [1, 2, 3])
         t(:, [1, 2, 4])
         t(:, [1, 3, 4])
         t(:, [2, 3, 4])];

faces = sort(faces, 2);
[foo, ix, jx] = unique(faces, 'rows');
if (isoctavemesh)
    u = unique(jx);
    qx = u(hist(jx, u) == 2);
else
    vec = histc(jx, 1:max(jx));
    qx = find(vec == 2);
end

nn = max(t(:));
ne = size(t, 1);
facenb = zeros(size(t));

% now I need to find all repeatitive elements
% that share a face, to do this, unique('first')
% will give me the 1st element, and 'last' will
% give me the second. There will be no more than 2

% doing this is 60 times faster than doing find(jx==qx(i))
% inside a loop

if (isoctavemesh || datenum(version('-date')) > datenum('January 27 2006')) % compare to matlab 7.2
    [ujx, ii] = unique(jx, 'first');
    [ujx, ii2] = unique(jx, 'last');
else
    ujx = unique(jx);
    [t1, ii2] = ismember(ujx, jx);
    [t1, ii] = ismember(ujx, flipwd(jx(:)));
    ii = length(jx) - ii + 1;
end

% iddup is the list of all pairs that share a common face

iddup = [ii(qx) ii2(qx)];
faceid = ceil(iddup / ne);
eid = mod(iddup, ne);
eid(eid == 0) = ne;

% now rearrange this list into an element format

for i = 1:length(qx)
    facenb(eid(i, 1), faceid(i, 1)) = eid(i, 2);
    facenb(eid(i, 2), faceid(i, 2)) = eid(i, 1);
end

% facenb may contain 0s, that just means the corresponding
% face is a boundary face and has no neighbor.

% if the second option is 'surface', I am going to find
% and return surface patches only

if (nargin == 2)
    if (strcmp(opt, 'surface'))
        facenb = faces(find(facenb == 0), :);
    elseif (strcmp(opt, 'rowmajor'))
        index = [1:length(faces)];
        index = (reshape(index, [], 4))';
        faces = faces(index(:), :);
        facenb = faces(find(facenb' == 0), :);
    else
        error(['supplied option "' opt '" is not supported.']);
    end
end