Building a connected graph from a triangulation

Question

TL;DR: I have a bunch of tetrahedra and I want to know which are its 4 (or less) neighboring (face-sharing) tetrahedra and I do that by checking the points in 3D that they share. This is slow.

---

I am building a connected graph from a triangulation. The graph has a structure as[1]:

Graph.element.nodeId
             .neighbours
     .nodes.positions

And the output of the triangulation has 2 matrices, TRI and points, the first one is a Ntri x 4 array, in each row having the node Ids for each tetrahedron, while the secodn one is a list of point Npoint x 3 size.

I am currently building the graph as in the below code, but it is absurdly slow for any decent sized mesh. The single line that takes almost all time is the find line (marked in the code), the part where the neighbors of the current element are being found.

The current algorithm takes for each tetrahedron all its nodes, and then finds all other tetrahedra that also contain these same nodes. Then it filters out all of the tetrahedra that does not contain 3 same nodes as the current one, leaving only the current tetrahedron's neighbors.

function graph=trimesh2graph(TRI,points)
nD=size(points,2);
% preallocate.
graph.elements(size(TRI,1)).neighbours=[];


% For each element, find its node Ids and neighbouring elements
for ii=1:size(TRI,1)
    nodeids=TRI(ii,:);
    elem=[];
    for jj=1:(nD+1)
        [iind,~]=find(nodeids(jj)==TRI); % this is 80% of the time
        elem=[elem; iind];
    end
    u = unique(elem);

    % of all elements that have shared nodes with the current element, 
    % get only the ones that have 3 shared nodes.
    graph.elements(ii).neighbours = int32((u(histc(elem,u)==nD)));


   % some other code
end
% some other code

Run this script with demo data:

x = gallery('uniformdata',[30 1],0);
y = gallery('uniformdata',[30 1],1);
z = gallery('uniformdata',[30 1],2);
vertices=[x,y,z];
TRI= delaunay(x,y,z)
trimesh2graph(TRI,vertices);

How can I improve the performance of this code? I am expecting that it will require a different approach, rather than just faster commands. I looked a bit to voronoi diagrams, but it seems that the finding (find) needs to be done anyway.

Note: I don't necessarily need MATLAB code, if you know an algorithm to solve the problem please answer, I'll code it in MATLAB.

---

[1] yes, it is better to store this info in 1D arrays. I will later, but easier to understand now with the current structure.

Answer 1

Ah, of course there is an inbuilt for this

if using delaunay() instead of the newer version with its own class, just do

neighbors(triangulation(TRI,points))

To get the neighbors. Boundary elements will have NaNs to fill in the matrix.