How to plot solid tetrahedrons in MATLAB?

Hi all,

I have information of a mesh including coordinates of nodes and connectivity, how can I plot solid tetrahedrons in MATLAB? The nodes and connectivity are like this (giving only one tetrahedron as an example):

coordinate =
  121.9864   -4.9838  223.7500
  121.9864   -4.9838   13.7500
  121.9864    4.9838   13.7500
  121.9864    4.9838  223.7500
  160.0000         0  223.7500
  160.0000         0   13.7500

connectivity =
     6     5     4     1

I tried this function:

function PlotMesh(coordinates,nodes,show)

if nargin == 2
    show = 0 ;
end
dimension = size(coordinates,2) ;  % Dimension of the mesh
nel = size(nodes,1) ;                  % number of elements
nnode = length(coordinates) ;          % total number of nodes in system
nnel = size(nodes,2);                % number of nodes per element
%
% Initialization of the required matrices
X = zeros(nnel,nel) ;
Y = zeros(nnel,nel) ;
Z = zeros(nnel,nel) ;

if dimension == 3   % For 3D plots
    if nnel==4 % surface in 3D
        for iel=1:nel
            nd = nodes(iel,:) ;

              X(:,iel)=coordinates(nd,1);    % extract x value of the node
              Y(:,iel)=coordinates(nd,2);    % extract y value of the node
              Z(:,iel)=coordinates(nd,3) ;   % extract z value of the node
          end
          % Plotting the FEM mesh
          figure
          fill3(X,Y,Z,'w')
          rotate3d ;
          title('Finite Element Mesh') ;
          axis off ;
          % display Node numbers and Element numbers
          if show ~= 0
              k = 1:nnode ;
              nd = k' ;
              for i = 1:nel
                  text(X(:,i),Y(:,i),Z(:,i),int2str(nd(i)),....
                      'fontsize',8,'color','k');
                  text(sum(X(:,i))/4,sum(Y(:,i))/4,sum(Z(:,i))/4,int2str(i),.....
                      'fontsize',10,'color','r') ;
              end
          end

      elseif nnel==8  % solid in 3D
          fm = [1 2 6 5; 2 3 7 6; 3 4 8 7; 4 1 5 8; 1 2 3 4; 5 6 7 8];
          XYZ = cell(1,nel) ;
          for e=1:nel
              nd=nodes(e,:);
              X(:,e) = coordinates(nd,1) ;
              Y(:,e) = coordinates(nd,2) ;
              Z(:,e) = coordinates(nd,3) ;
              XYZ{e} = [X(:,e)  Y(:,e) Z(:,e)] ;
          end
          % Plot FEM mesh
          figure
          set(gcf,'color','w')
          axis off
          cellfun(@patch,repmat({'Vertices'},1,nel),XYZ,.......
              repmat({'Faces'},1,nel),repmat({fm},1,nel),......
              repmat({'FaceColor'},1,nel),repmat({'w'},1,nel));
          view(3)
          set(gca,'XTick',[]) ; set(gca,'YTick',[]); set(gca,'ZTick',[]) ;
          % display Node numbers and Element numbers
          if show ~= 0
              k = 1:nnode ;
              nd = k' ;
              for i = 1:nel
                  text(X(:,i),Y(:,i),Z(:,i),int2str(nd(i)),....
                      'fontsize',8,'color','k');
                  text(sum(X(:,i))/8,sum(Y(:,i))/8,sum(Z(:,i))/8,int2str(i),.....
                      'fontsize',10,'color','r') ;
              end
          end
      end

elseif dimension == 2           % For 2D plots
    disp('2d case skipped')
end

However it plots 2 triangular surfaces instead of 1 tetrahedron solid. Can anyone help?

Thanks!

[Walter Roberson]

Did you look at tetramesh() ?

[Xh Du]

Hi Walter,

Yeah I checked tetramesh and it works, but I think the above function would work for all element types? Also I'd like to understand why it plots 2 surfaces instead of 1 tetrahedron. Thanks!

Answer 1

[Precise Simulation]

With the coordinates and connectivities in your example, and the Matlab FEM Toolbox, you can create a grid struct as follows

    grid.p = coordinate';   
    grid.c = connectivity';  

    % Reconstruct adjacency, boundary, and subdomain fields.
    grid.a = gridadj( grid.c, size(grid.p,1) );   
    grid.b = gridbdr( grid.p, grid.c, grid.a );   
    grid.s = ones(1,size(grid.c,2));

And simply plot and visualize the grid with

    plotgrid( grid )

If you expand the data struct to a full FEA struct you can also utilize all postprocessing functionality to plot iso-surface, slice, arrow, and contour plots for unstructured Matlab meshes and grids.

Answer 2

[John D'Errico]

You don't plot a solid tetrahedron. You plot 4 facets of a tetrahedron. Done. WTP? Use patch.

There are some things you can do to improve on that, since if you have a large tessellation, there will be MANY facets, so 4 times as many facets as you have tetrahedra.

First, you can recognize that there is no need to plot a shared facet in a tessellation twice. If done carefully, you will remove all replicated facets. Since in a tessellation, most of the facets will be internal, and EVERY internal facet will be duplicated exactly twice, this will cut the overall plotting task by almost 50%.

Next if the facets are solid color, not transparent, then you cannot see the inside anyway. So all that matters is to completely remove ALL of the shared facets! If you cannot see them, then don't plot them at all. So the trick is to identify only the list of what I might call exterior facets. An exterior fact appears in the list of all facets only one time, not twice.

There are a couple of tricks in all of this, but the task is really far easier than it may sound.

XYZ = randn(1000,3);
tess = delaunayn(XYZ);

facets = [tess(:,[1 2 3]);tess(:,[1 2 4]);tess(:,[1 3 4]);tess(:,[2 3 4])];

size(tess)
ans =
        6464           4

size(facets)
ans =
       25856           3

So 1000 vertices, which resulted in 6464 tetrahedra, but 23856 facets.

If we were using transparency, then we could see into the inside of the tessellation. So we would just use unique, with the rows option to drop out the duplicates of those facets. So this would be correct:

% do this if the facets are translucent
facets = unique(sort(facets,2),'rows');

We would now just plot the facets using a call to patch.

But with opaque facets, we cannot see into the inside at all, so then we would just drop all facets that appeared twice in the list. How do we find the replicated facets? This is essentially a careful pair of sorts. This does nothing to the facets themselves, but makes them identifiable as dups.

% do this instead if the facets are opaque
facets = sortrows(sort(facets,2));

So now a simple test using diff and then a find can identify the dups, to then be removed completely

df = diff(facets,[],1);
dups = find(all(df == 0,2));
facets([dups,dups+1],:) = [];

size(facets)
ans =
    68     3

So even though there were tens of thousands of facets in the overall list, only 68 of them matter if the facets are opaque.

In either case, once you have the list of facets to plot, then a single call to patch will plot them all. Which of the sorts you use above, thus either the sort rows or the unique, will depend on whether you can see inside.

[Hugh Collett]

Hi John, I'm not sure where to ask this but I'm hoping you may be able to help me please. I have a set of 3d points and I did delaunay(x,y) to give a surface triangulation. I didn't do delaunay(x,y,z) because I wanted triangular not tetrahedral mesh. So I got how the points connected by the connectivity list this triangulation gave me. I then worked out the centre point of the triangles because I want to add a node at an offset of these centre points.

I'm now hoping to do a tetrahedral mesh of the original triangles connecting to my central points. I know all the points, but I'm not sure how to connect them. Could you advise me how to do this please, I would be very grateful I'm really stuck. If I haven't explained myself well let me know. I'm basically trying to do something similar to the image I've attached.

Thanks

Hugh

[John D'Errico]

I think you are confused about several things concerning a mesh, about delaunay, etc. Or possibly it may just be you have not explained what you do know. But when an explanation is confusing, that almost always means you are yourself confused.

First, if you have a point cloud in 3-d, then using delaunay(x,y) does not create a mesh in 3-d, since you never told it about z. MATLAB is just humming along, doing what you told it to do. It triangulates the set of points you gave it, in TWO DIMENSIONS, thus the (x,y) plane. z does not exist as far as MATLAB is concerned.

A surface triangulation in 3-d can be arrived at using a convex hull.

x = rand(1,10);
y = rand(1,10);
z = rand(1,10);
tri = convhull(x,y,z)
tri =
     1     3     7
     1     7    10
     1     9     3
     1    10     9
     2     3     4
     2     4     8
     2     7     3
     2     8     7
     3     9     4
     4     9     8
     7     8    10
     8     9    10

The result is a set of triangles as a surface, referencing the point cloud (x,y,z).

Note that these triangles will be a convex representation of that surface. Had we used delaunay

tet = delaunay(x,y,z)
tet =
     9     5     4     8
    10     6     5     2
     8     5     4     2
     5     3     4     2
     7     6    10     2
     1     7     6    10
     4     3     5     9
     9     3     5     2
     6     3     9     2
     6     9     5     2
     7     3     6     2
     6     3     7     1
     9     3     6     1
    10     5     8     2
     7    10     8     2
    10     5     9     8
     6     5     9    10
     9     1     6    10

this would be a volume tessellation of that point cloud, thus now simplexes in 3 dimensions, tetrahedrae.

The surface of that tessellation is the same as the convex hull. since a delaunay tessellation of a point cloud will always produce a convex result.

So, we COULD have used delaunay to compute a delaunay tessellation in 3-d, then we would be forced to extract the surface triangles from that mesh. (This is a fine mesh I've gotten us into.) And while it is not difficult to extract that surface from the 3-d tessellation, it is easier to just use convhull, if that is your end point goal.

In either case, you need to recognizer that what will come out is a convex representation of some point cloud. That may or may not be an issue, since it appears you MAY be looking to deal with non-convex domains. I say that because the pictures drawn were not of convex domains.

Finally, you ask about appending new points to each triangle on the surface. I'm not exactly sure what your goal really is here.

First, I'll suggest it is better to combine your point cloud into ONE array.

xyz = [x;y;z]';

Rather than store the cloud as three independent vectors, where one of them can at some point be corrupted by some later accident, store it in one array. Learn to use arrays in MATLAB.

So I've created one 10x3 array here, each row of which represents the (x,y,z) coordiantes of one point in your cloud. I can again form the convex hull, now using convhulln.

tri = convhulln(xyz);

Next, you ask how to convert those triangles into tetrahedrae, by now connecting each facet to a different point?

I might connect each surface triangular facet with the interior of the mesh like this:

nxyz = size(xyz,1);
ntri = size(tri,1);
xyznew = [xyz;mean(xyz,1)];
tet = [tri,repmat(nxyz + 1,[ntri,1])];

tet
tet =
     8     4     2    11
     4     3     2    11
     7     8     2    11
     8     7    10    11
     3     7     2    11
     9     3     4    11
     9     8    10    11
     9     4     8    11
     1     7     3    11
     9     1     3    11
     7     1    10    11
     1     9    10    11

So you can see that here I have created a volume tessellation from the convex surface triangulation. Each triangular facet is now connected to that same new central point in the interior of that cdomain.

I am not sure this is truly what you are asking, based on those pictures. Your words are confusing though. However, the general idea of converting a set of triangular facets into tetrahedrae is the same as what I have shown.