generatemesh not respecting boundaries in 2023

I am having some issues with the generatemesh function. I used the following code to mesh a 2D plate with a crack in it. It works perfectly in v2022a, but does not produce the expected result in v2023b.

% boundary polyshape from vertices
x = [0 150 150 75 150 150 0];
y = [0 0 72.5 75 77.5 150 150];
pgon = polyshape({x}, {y});

% triangulation of geometry
T = triangulation(pgon);
tnodes = T.Points';
telements = T.ConnectivityList';
model = createpde;
geometryFromMesh(model, tnodes, telements);

% generate mesh
h = 20;
mesh = generateMesh(model, 'GeometricOrder', 'linear', 'Hmax', h);
pdemesh(model)

In v2022a, it produces the following mesh. The crack boundary is respected.

However, in v2023b, it does not respect the geometry of the crack, and meshes over it!

I'm aware that I can set a finer element size at vertices and edges, for example, at the crack tip, but this is not working properly either. For example, using

mesh = generateMesh(model, 'GeometricOrder', 'linear', 'Hmax', h, 'Hvertex', {4,h/20});

to set a finer mesh at the crack tip, improves the situation, but the crack tip is still not in the right place, as shown in the zoom in on the crack tip in the image below (it should be at coordinate 75,75)!

Is there any way I can use generatemesh in v2023b and have it still respect the boundaries?

Answer 1

[Christopher]

I have found a workaround for this specific problem. It is not pretty. The principle is to mesh the geometry in sections and then merge the meshes. To get generatemesh to respect the boundaries, this required four quadrants around the crack tip. You also have to be careful to generate conforming elements along shared edges. Then I look for duplicate nodes and merge them, making appropriate updates to the node numbering in the connectivity matrix.

% boundary polyshape from vertices
pgon1 = polyshape({[0 75 75 0]}, {[0 0 75 75]}); % lower left quadrant
pgon2 = polyshape({[75 150 150 75]}, {[0 0 72.5 75]}); % lower right quadrant
pgon3 = polyshape({[0 75 75 0]}, {[75 75 150 150]}); % upper left quadrant
pgon4 = polyshape({[75 150 150 75]}, {[75 77.5 150 150]}); % upper right quadrant

figure(1)
plot(pgon1); hold on
plot(pgon2);
plot(pgon3);
plot(pgon4); axis equal

% triangulation of geometry, quadrant by quadrant
T1 = triangulation(pgon1);
T2 = triangulation(pgon2);
T3 = triangulation(pgon3);
T4 = triangulation(pgon4);
model1 = createpde;
model2 = createpde;
model3 = createpde;
model4 = createpde;
geometryFromMesh(model1, T1.Points', T1.ConnectivityList'); % lower left quadrant
geometryFromMesh(model2, T2.Points', T2.ConnectivityList'); % lower right quadrant
geometryFromMesh(model3, T3.Points', T3.ConnectivityList'); % upper left quadrant
geometryFromMesh(model4, T4.Points', T4.ConnectivityList'); % upper right quadrant

% generate mesh, quadrant by quadrant
% **BE CAREFUL TO GENERATE CONFORMING ELEMENTS ALONG SHARED EDGES**
h = 20;
mesh1 = generateMesh(model1, 'GeometricOrder', 'linear', 'Hmax', h, 'Hvertex', {4, h/10}, 'Hgrad', 1.2);
mesh2 = generateMesh(model2, 'GeometricOrder', 'linear', 'Hmax', h, 'Hvertex', {2, h/10}, 'Hgrad', 1.2);
mesh3 = generateMesh(model3, 'GeometricOrder', 'linear', 'Hmax', h, 'Hvertex', {2, h/10}, 'Hgrad', 1.2);
mesh4 = generateMesh(model4, 'GeometricOrder', 'linear', 'Hmax', h, 'Hvertex', {1, h/10}, 'Hgrad', 1.2);

figure(2)
pdemesh(model1); hold on
pdemesh(model2)
pdemesh(model3)
pdemesh(model4); axis equal

% merge meshes by equivalencing duplicate nodes
[nodes, elements] = mergeMeshes(mesh1.Nodes', mesh2.Nodes', mesh1.Elements', mesh2.Elements', 1e-4);
[nodes, elements] = mergeMeshes(nodes, mesh3.Nodes', elements, mesh3.Elements', 1e-4);
[nodes, elements] = mergeMeshes(nodes, mesh4.Nodes', elements, mesh4.Elements', 1e-4);

figure(3)
triplot(elements, nodes(:,1), nodes(:,2), 'k');
axis equal

% merge mesh function
function [nodes, elements] = mergeMeshes(nodes1, nodes2, elements1, elements2, tolerance)
    nodes = [nodes1; nodes2];
    elements = [elements1; elements2 + size(nodes1,1)];

    % find duplicate nodes (within tolerance)
    duplicateNodes = [];
    mergeWith = [];
    for i = 1:size(nodes,1)
        lia = ismembertol(nodes, nodes(i,:), tolerance, 'ByRows', true);
        if(sum(lia) > 1)
            duplicates = find(lia);
            duplicateNodes = [duplicateNodes duplicates(2:end)];
            mergeWith = [mergeWith i];
        end
    end

    duplicateNodes = duplicateNodes(1:length(duplicateNodes)/2);
    mergeWith = mergeWith(1:length(mergeWith)/2);

    % merge duplicate nodes
    nodes(duplicateNodes',:) = [];

    % renumber nodes in elements matrix
    [duplicateNodes, sortOrder] = sort(duplicateNodes, 'descend');
    mergeWith = mergeWith(sortOrder);
    for i = 1:length(duplicateNodes)
        elements(elements == duplicateNodes(i)) = mergeWith(i);
        elements(elements >= duplicateNodes(i)) = elements(elements >= duplicateNodes(i)) - 1;
    end
end

Answer 2

[Jiexian Ma]

Hi Christopher,

I encoutered the same issue when using function generateMesh.

You could use Im2mesh package to generate mesh instead.

Matlab code

You could check demo14 to demo17 in Im2mesh package. I provide a few examples.