1

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.

Ander Biguri
  • 35,140
  • 11
  • 74
  • 120
  • The `delaunayTriangulation` object you pass to your function has its own `neighbors` method that will be hard to beat on performance to solve the problem you use to introduce your question. Does the `struct` you describe take the form it does for any other necessary purpose? – Will Sep 07 '18 at 10:36
  • @Will there is no `delaunayTriangulation` object passed, but I'll check that version of the triangulation in MATLAB. Yes, indeed I need this information for other own code stuff. – Ander Biguri Sep 07 '18 at 10:37
  • @Will to be clear, I had a typo from copy pasting from somewhere else, I am using `delaunay()`. But maybe I should use that one – Ander Biguri Sep 07 '18 at 10:46
  • 1
    That clears things up! Anything the `delaunayTriangulation` class can do should be used directly when performance is the objective. There aren't many triangulation problems left that are likely to need to be solved from first-principles. – Will Sep 07 '18 at 10:54
  • @Will indeed there are not! I am exploring, that will probably have the answer ;) thanks – Ander Biguri Sep 07 '18 at 10:55
  • Are all your tetrahedra the same size (i.e. the only reaon why it can have less than 4 neighbours is because it is close to an edge)? Otherwise you might be able to define neighbours by calculating their centres and the distance between these centre points. – Floris Sep 07 '18 at 11:05
  • @FlorisSA ah, yes, but indeed this is not the case. arbitrary-sized and shaped ones, as this is not for a FEM, so I dont care if they are ugly. Found the solution though ;) – Ander Biguri Sep 07 '18 at 11:08

1 Answers1

0

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.

Ander Biguri
  • 35,140
  • 11
  • 74
  • 120