function BdrSides = computeBoundarySides(c4n, n4e)

allSides = [n4e(:,[2 3]); n4e(:,[3 1]); n4e(:,[1 2])];
[~, FirstInd, Back] = unique(sort(allSides, 2), 'rows', 'first');
[~, sortInd] = sort(FirstInd);
S = length(sortInd);
sideNo(sortInd) = 1:S;
s4e = reshape(sideNo(Back), size(n4e));

% ...



end

