demo21 of Im2mesh package

demo21 - 2D mesh for periodic boundary conditions

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

Overview Setup Create geometry Polyshape Create matrix Create fiber Cell array of polyshape Polygonal Boundaries Convert to polygonal boundaries Add intersect points (optional) Simplify boundaries Refine outermost boundary Generate mesh A bad case

Overview

In this demo, we will demostrate creating a 2D mesh for periodic boundary conditions (PBC). PBC requires node matching on the opposite edges of a Representative Volume Element (RVE).

Node matching for PBC:

Left and right boundaries: For every node on the left boundary, there must be a corresponding node on the right boundary with the exact same y-coordinate
Top and bottom boundaries: For every node on the bottom boundary, there must be a corresponding node on the top boundary with the exact same x-coordinate.

We will use function refineOuterBoundary, function checkPBCNode, and the opt.tf_preserve option of function bounds2mesh.

Setup

Before we start, please set folder "Im2mesh_Matlab" as your current folder of MATLAB.

clearvars

Set default image size (optional).

x = 250; y = 250; width = 300; height = 300;
set(groot, 'DefaultFigurePosition', [x,y,width,height])
% To reset:
% set(groot, 'DefaultFigurePosition', 'factory')

Function bounds2mesh use a mesh generator called MESH2D (developed by Darren Engwirda). We can use the following command to add the folder 'mesh2d-master' to the path of MATLAB.

addpath(genpath('mesh2d-master'))

Create geometry

Polyshape

We use Matlab polyshape objects to create a simple RVE for demostration. Suppose we are interested in fiber composites.

Start with a square.

vertex = 10*[ 0 0; 1 0; 1 1; 0 1 ];
psSq = polyshape(vertex);

Create a circle.

t = 0.05:0.1:2*pi;
x1 = 5 + 1.5*cos(t);
y1 = 5 + 1.5*sin(t);
psC0 = polyshape(x1,y1);

Translate

psC1 = translate( psC0, [5.1, 1.5] );
psC2 = translate( psC0, [-4.9, 1.5] );
psC3 = translate( psC0, [1, 5.6] );
psC4 = translate( psC0, [1, -4.4] );

Plot them together.

figure

plot( [psSq; psC0; psC1; psC2; psC3; psC4] );

axis equal

We saw that they are overlapped.

Create matrix

psMatrix = subtract( psSq, psC0 );

psMatrix = subtract( psMatrix, psC1 );

psMatrix = subtract( psMatrix, psC2 );

psMatrix = subtract( psMatrix, psC3 );

psMatrix = subtract( psMatrix, psC4 );

figure

plot( psMatrix );

axis equal

Create fiber

psC1 = intersect( psC1, psSq );
psC2 = intersect( psC2, psSq );
psC3 = intersect( psC3, psSq );
psC4 = intersect( psC4, psSq );

Plot them together.

figure

plot( [psC0; psC1; psC2; psC3; psC4] );

axis equal

Union of fibers

psFiber = psC0;

psFiber = union( psFiber, psC1);

psFiber = union( psFiber, psC2);

psFiber = union( psFiber, psC3);

psFiber = union( psFiber, psC4);

figure

plot( psFiber );

axis equal

Cell array of polyshape

We put them into a cell array - psCell. Each element in psCell represent different parts. It's like labelling.

psCell = { psMatrix; psFiber };

figure

for i = 1: length(psCell)

subplot(2,1,i);

plot( psCell{i} ); axis equal;

title(['psCell\{' num2str(i), '\}']);

end

Polygonal Boundaries

Convert to polygonal boundaries

Convert the cell array of polyshape 'psCell' to a cell array of polygonal boundaries 'bounds'

bounds = polyshape2bound(psCell);

% plot boundaries and show all vertices

plotBounds(bounds,false,'ko-')

Add intersect points (optional)

Polyshape objects are not able to find intersect points automatically. Therefore, we use function addIntersectPnts to add intersect points to 'bounds'.

Note that if we don't add intersect points, mesh generation may fail.

tol_intersect = 1e-6;   % distance tolerance for intersect
bounds = addIntersectPnts( bounds, tol_intersect );

Simplify boundaries

We want to reduce the number of vertices on the circle.

boundsCtrlP = getCtrlPnts( bounds );

tolerance = 0.05; % for Douglas-Peucker polyline simplification

boundsReduced= simplifyBounds( boundsCtrlP, tolerance );

% show all vertices

plotBounds(boundsReduced,false,'ko-')

Refine outermost boundary

We use function refineOuterBoundary to add mesh seed to the outermost boundary.

In the function refineOuterBoundary, the seeds are inserted into each edge according to the following equation. 'len' is the length of an edge. 'targetSpace' is the target spacing between seeds. 'actualSpace' is the actual spacing between seeds.

numSegment = round( len / targetSpace );

actualSpace = len / numSegment;

targetSpace = 0.6;

newB = refineOuterBoundary( boundsReduced, targetSpace );

plotBounds(newB,false,'ko-')

Generate mesh

To generate mesh, we use function bounds2mesh. Function bounds2mesh uses MESH2D (developed by Darren Engwirda) as mesh generator.

Here, we set opt.tf_preserve as true, which enables MESH2D to retain the initial edges without further subdivision. Therefore, we can get 2d mesh with node matching.

hmax = 500; 
grad_limit = 0.25;
opt = [];
opt.disp = inf;     % silence verbosity
opt.tf_preserve = 1;
[vert,tria,tnum,vert2,tria2] = bounds2mesh( newB, hmax, grad_limit, opt );

Plot mesh.

color = 1;

plotMeshes( vert,tria,tnum, color);

We use function checkPBCNode to check whether the mesh satisfies node matching for PBC.

checkPBCNode( vert );
PBC Check Passed: All opposite boundary nodes are aligned.

We use function tricost to check the mesh quality.

figure

tricost( vert, tria );

A bad case

When creating geometry using polyshape objects, we need to be careful. Here, a bad case is demostrated.

clearvars

Start with a square.

vertex = 10*[ 0 0; 1 0; 1 1; 0 1 ];
psSq = polyshape(vertex);

Create a circle.

t = 0.05:0.1:2*pi;
x1 = 5 + 1.5*cos(t);
y1 = 5 + 1.5*sin(t);
psC0 = polyshape(x1,y1);

Translate

psC1 = translate( psC0, [5.1, 1.5] );
psC2 = translate( psC0, [-4.9, 1.5] );
psC3 = translate( psC0, [1, 5.6] );
psC4 = translate( psC0, [1, -4.4] );

Plot them together.

figure

plot( [psSq; psC0; psC1; psC2; psC3; psC4] );

axis equal

Create fiber

psC1 = intersect( psC1, psSq );

psC2 = intersect( psC2, psSq );

psC3 = intersect( psC3, psSq );

psC4 = intersect( psC4, psSq );

psFiber = psC0;

psFiber = union( psFiber, psC1);

psFiber = union( psFiber, psC2);

psFiber = union( psFiber, psC3);

psFiber = union( psFiber, psC4);

figure

plot( psFiber );

axis equal

Create matrix

psMatrix = subtract( psSq, psFiber );

figure

plot( psMatrix );

axis equal

We saw that the outer edges in the matrix are wierd. Some edges are not supposed to there. If we proceed with this geometry, the 2d mesh generation would fail.

% % reset image size
% set(groot, 'DefaultFigurePosition', 'factory')
% 
% % end of demo