Large, mostly empty data structure for tetrahedron

Question

I'm trying to link the node id of every face in a tetrahedron with it's corresponding tetra id.

tetras = [1 2 3 4  % Tetra 1
          5 6 7 8] % Tetra 2

For tetra 1, there are four faces:

faces = [1 2 3; 1 2 4; 1 3 4; 2 3 4] % Notice these are sorted

Then I'd like to store these in a data structure:

tet_for_face = cell(8,8,8) % 8 allows for the maximum node id

tet_for_face{1,2,3} = 1;
tet_for_face{1,2,4} = 1;
tet_for_face{1,3,4} = 1;
tet_for_face{2,3,4} = 1;

This means that I can find the tetra ID of any particular face in O(1):

tet_for_face{2,3,3}
ans = []
tet_for_face{2,3,4}
ans = 1

The problem with this approach is that it requires contiguous memory. As my mesh gets larger, I run out of memory:

cell(1000, 1000, 1000)
??? Error using ==> cell
Out of memory. Type HELP MEMORY for your options.

I've also played around with using nested cells:

tet = cell(num_nodes, 1);
tet2 = cellfun(@(x) cell(num_nodes, 1), tet, 'UniformOutput', 0);
tet3 = cellfun(@(x) cellfun(@(y) cell(num_nodes, 1), x, 'UniformOutput', 0), tet2, 'UniformOutput', 0);

tet3{2}{3}{4} = 1;
...

Although this works for small meshes, and doesn't require contiguous memory (AFAIK), it has a nasty habit of crashing MATLAB with N=1000.

Any ideas?

Answer 1

After a bit of playing with sparse arrays (which can only be 1D or 2D, not 3D), and not getting anywhere, I decided to go with containers.Map (HashMap).

I used string keys, and the fastest way I found of producing them I found was using sprintf (rather than int2str or mat2str)

Sample code:

tet = containers.Map;
for tetra_id in tetras 
    for face in faces
        face_key = sprintf('%d ', face);
        tet(face_key) = tetra_id;

This gives me a map like so:

tet('1 2 3') = 1

Answer 2

You can use sparse matrices to deal with many problems arising with meshes. It depends on what you want to do with this data structure in practice, but here is one example:

% tetras and faces are transposed - column-wise storage
tetras = [1 2 3 4; 5 6 7 8]';
faces = [1 2 3; 1 2 4; 1 3 4; 2 3 4]';

ntetras    = size(tetras, 2);
nfaces     = size(faces, 2);
nfacenodes = size(faces, 1);

% construct face definitions for all tetras
tetras_faces = reshape(tetras(faces, :), nfacenodes, ntetras*nfaces);

% assign the faces to tetras keeping the face id within the tetra, if you need it
enum_faces  = repmat(1:ntetras*nfaces, nfacenodes, 1);

% create a sparse matrix connecting tetra faces to tetras.
% Every column contains 3 non-zeros - 1 for every node in a face
% The number of matrix columns is ntetras*nfaces - 4 columns for every element.
A = sparse(tetras_faces, enum_faces, 1);

Now to extract the information you need you can multiply A by a vector holding the information about the face you are looking for:

v = sparse(ntetras*nfaces, 1);
v([1 2 3]) = 1;
tetra_id = ceil(find(A*v==nfacenodes)/nfaces)

Note that this is just an example. You can extract much more useful information this way, and you could perform more sophisticated searches using matrix-matrix multiplication instead of matrix-vector multiplication.