demo18 of Im2mesh package

demo18 - 2d image to tetrahedral mesh
Cite as: Ma, J., & Li, Y. (2025). Im2mesh: A MATLAB/Octave package for generating finite element mesh based on 2D multi-phase image (2.1.5). Zenodo. https://doi.org/10.5281/zenodo.14847059
Jiexian Ma, mjx0799@gmail.com. Project website.
Table of Contents

Note

I suggest familiarizing yourself with Im2mesh_GUI before learning Im2mesh package. With graphical user interface, Im2mesh_GUI will help you better understand the workflow and parameters of Im2mesh package.
Im2mesh_GUI: https://www.mathworks.com/matlabcentral/fileexchange/179684-im2mesh_gui-2d-image-to-finite-element-mesh
If you are using Im2mesh package in MATLAB, you need to install MATLAB Image Processing Toolbox and Mapping Toolbox because Im2mesh package use a few functions in these toolboxes.

Overview

There are 4 methods to generate tetrahedral or hexahedral mesh based on a 2D image.
Method 1. Use MATLAB built-in function generateMesh to create tetrahedral mesh. We will demostrate this method in demo18.
Method 2. CAD file. Import the dxf file created by Im2mesh into finite element software.
Method 3. Gmsh. Modify the geo file created by Im2mesh and then generate 3D mesh via Gmsh.
Method 4. fTetWild. See demo20.
To run demo18, you need to have MATLAB Partial Differential Equation Toolbox.

Setup

Before we start, please set folder "Im2mesh_Matlab" as your current folder of MATLAB.
clearvars
Set default image size.
x = 250; y = 250; width = 400; height = 400;
set(groot, 'DefaultFigurePosition', [x,y,width,height])
% To reset:
% set(groot, 'DefaultFigurePosition', 'factory')
We add the folder 'mesh2d-master' to the path of MATLAB.
addpath(genpath('mesh2d-master'))

Import image

im = imread('logo.tif');
if size(im,3) == 3; im = rgb2gray( im ); end
imshow( im,'InitialMagnification','fit' );

Get boundries from image

We have demostrate how to do this in example3 of demo01.mlx
opt = []; % reset opt
opt.tolerance = 1;
opt.tf_mesh = false;
bounds = im2mesh( im, opt );
% bounds is the simplified 2d polygonal boundaries
plotBounds( bounds );

Create 3d pde model object

We can use function bounds2pde3d to create Matlab 3d pde model object.
'height' is the extrusion height in the Z direction.
'scale_factor' is the scale factor for the node coordinites. We can set it as the physical dimensions of a pixel in the image.
height = 50;
scale_factor = 1;
model3d = bounds2pde3d( bounds, height, scale_factor )
model3d =
PDEModel with properties: PDESystemSize: 1 IsTimeDependent: 0 Geometry: [1×1 DiscreteGeometry] EquationCoefficients: [] BoundaryConditions: [] InitialConditions: [] Mesh: [] SolverOptions: [1×1 pde.PDESolverOptions]
We can visualize the geometry.
pdegplot(model3d)
Awesome!

Generate mesh

We can use Matlab built-in function generateMesh to generate tetrahedral mesh.
generateMesh( model3d, 'Hmax',500, 'Hmin', 5, 'Hgrad', 1.2, 'GeometricOrder', 'quadratic' )
ans =
FEMesh with properties: Nodes: [3×22776 double] Elements: [10×14642 double] MaxElementSize: 500 MinElementSize: 5 MeshGradation: 1.2000 GeometricOrder: 'quadratic'

Find elements in different phases

Before plotting the mesh, we need to find elements in different phases, also known as cells or parts.
numCell = model3d.Geometry.NumCells;
elementIDsCell = cell( [1, numCell] );
for i = 1: numCell
elementIDsCell{i} = findElements( model3d.Mesh, "region", Cell=i );
end

Plot mesh

Now we can plot the mesh. Different phases are labelled with different colors.
figure
colors = lines(numCell); % colormap
for i = 1: numCell
pdemesh( model3d.Mesh.Nodes, ...
model3d.Mesh.Elements(:,elementIDsCell{i}), ...
'FaceColor', colors(i,:))
hold on
end
hold off
Awesome! We saw 3 phases in the mesh.

Change view angle

figure
colors = lines(numCell); % colormap
for i = 1: numCell
pdemesh(model3d.Mesh.Nodes, ...
model3d.Mesh.Elements(:,elementIDsCell{i}), ...
'FaceColor', colors(i,:))
hold on
end
view([-0 90])
hold off
Please check the following link about mesh refinement.
https://www.mathworks.com/help/pde/ug/pde.pdemodel.generatemesh.html

Evaluate mesh quality

We use Matlab built-in function meshQuality to check mesh qualtiy.
Q = meshQuality( model3d.Mesh );
mean(Q)
ans = 0.8509
min(Q)
ans = 0.1072
We can plot the histogram of mesh qualtiy.
figure
histogram(Q)
xlabel( "Element Shape Quality", "fontweight", "b" )
ylabel( "Number of Elements", "fontweight", "b" )

Extract mesh data

We can use the following code to extract mesh data from the pde model object for other purposes.
vert = model3d.Mesh.Nodes';
ele = model3d.Mesh.Elements';
numCell = model3d.Geometry.NumCells;
tnum = zeros( size(ele,1), 1 );
for i = 1: numCell
index = findElements( model3d.Mesh, "region", Cell=i );
tnum(index) = i;
end
We got 3 variables: vert, ele, tnum
vert: Mesh nodes. It’s a Nn-by-3 matrix, where Nn is the number of nodes in the mesh. Each row of vert contains the x, y, z coordinates for that mesh node.
ele: Mesh elements. Each row in tria contains the indices of the nodes for that mesh element.
tnum: Label of phase. Ne-by-1 array, where Ne is the number of elements.
tnum(j,1) = k; means the j-th element belongs to the k-th phase.

Plot mesh

We can also use function plotMesh3d in Im2mesh package to plot 3d mesh.
color_code = 1;
plotMesh3d( vert, ele, tnum, color_code );
Define cutting plane to check interior. You can ask Gemini to define more setting.
color_code = 1;
opt = []; % reset
opt.cutPlaneYZ = 180; % plot elements where x <= 180
plotMesh3d(vert, ele, tnum, color_code, opt);

Export as inp file (Abaqus)

We use function printInp3d to write 3d mesh to inp file. Please refer to the following link to learn more about exporting mesh: https://github.com/mjx888/writeMesh/blob/main/README.md
writeMesh demo05, demo06, and demo09 could be useful to you.
When using function printInp3d, we can specify element type, precision, file name, boundary nodes, and format.
Since the generated mesh consists of quadratic elements, we can set 'ele_type' as 'C3D10' Abaqus element. If the mesh consists of linear elements, we can set 'ele_type' as 'C3D4' instead.
ele_type = 'C3D10';
precision = 8;
file_name = 'test.inp';
opt = []; % reset
opt.tf_printMaxMinNode = 0;
opt.tf_printInterfNode = 0;
printInp3d( vert, ele, tnum, ele_type, precision, file_name, opt );
printInp3d Done! Check the inp file!

Export as bdf file (Nastran bulk data)

We use function printBdf3d to write 3d mesh to bdf file (Nastran bulk data, compatible with COMSOL). Note that function printBdf3d only works for linear element.
% convert to linear element before using function printBdf3d
if size(ele, 2) == 10
eleL = ele(:, 1:4);
[ vertL, eleL ] = delRedundantVertex( vert, eleL );
elseif size(ele, 2) == 4
% do nothing
else
error('Wrong dimension.')
end
precision = 8;
file_name = 'test_3d.bdf';
printBdf3d( vertL, eleL, tnum, [], precision, file_name );
printBdf3d Done! Check the bdf file!
% reset image size
set(groot, 'DefaultFigurePosition', 'factory')
% end of demo