Using generateMesh to refine a triangulation, without losing boundary points

I'm trying to refine a constrained Delaunay triangulation mesh automatically along the lines of Chew's algorithm. I was hoping to use MATLAB's generateMesh function, which does a nice job refining the mesh -- except that the mesh boundaries are not respected. Is there a way of forcing the generateMesh function to preserve the edge bounds?

In the example below, I have a smooth curve whose shape is distorted by generateMesh. It's not horrible in this case, but this is just an toy model for coastline polygons, where it's important to not lose any of the shape.

N= 100;

t= (1:N)/N*2*pi;

X= 1*cos(t)+ .15*cos(3*t)+ .1*sin(5*t);
X= X/range(X)*.9;
X= X- min(X) -.1;

Y= 1*sin(t)+ .5*cos(4*t)+ .4*sin(6*t);
Y= Y/range(Y)*.8;
Y= Y-min(Y)+.1;


% Create polyshape
P= polyshape(X,Y);

% Constrained triangulation of the polyshape
TRI= triangulation(P);

% Create mesh from original triangulation
model= createpde;
geometryFromMesh(model,TRI.Points',TRI.ConnectivityList');


figure(1)
tiledlayout('flow')

nexttile
plot(P)
title('Polyshape')

nexttile
pdemesh(model)
axis normal
title('Original Triangulation')



% Refine mesh
generateMesh(model,'Hmax',.3);

nexttile
pdemesh(model)
title('Distorted Shape')

Answer 1

[Samyuktha]

Hi Zel,

I understand that you wish to use the 'generateMesh' function to preserve the edge bounds.

Please refer to 'Refine Mesh on Specified Edges and Vertices' in the 'Examples' section of the documentation generateMesh, in order to generate a 2-D mesh with finer spots around the specified edges and vertices.

Hope it helps!

Answer 2

[Jiexian Ma]

I encountered the same issue when using function generateMesh. Setting 'Hmin' to a smaller value could improve the result, but some boundary points were still missing. I think the root cause may be numerical roundoff error of function delaunayTriangulation.

% --------------------------------------
N= 100;
t= (1:N)/N*2*pi;

X= 1*cos(t)+ .15*cos(3*t)+ .1*sin(5*t);
X= X/(max(X)-min(X)) *.9;
X= X- min(X) -.1;

Y= 1*sin(t)+ .5*cos(4*t)+ .4*sin(6*t);
Y= Y/(max(Y)-min(Y)) *.8;
Y= Y-min(Y)+.1;

% Create polyshape
P= polyshape(X,Y);

% Constrained triangulation of the polyshape
TRI= triangulation(P);

% Create mesh from original triangulation
model= createpde;
geometryFromMesh(model,TRI.Points',TRI.ConnectivityList');

% --------------------------------------
figure(1)
tiledlayout('flow')

% plot
nexttile
generateMesh(model,'GeometricOrder','linear','Hmax',.3);
pdemesh(model)
xlim([-0.2 0.9])
ylim([0 1])
numV = size(model.Mesh.Nodes,2);
title(['Set Hmax. Num of node:',num2str(numV)])

% plot
nexttile
generateMesh(model,'GeometricOrder','linear','Hmax',.3, 'Hmin', 0.0001);
pdemesh(model)
xlim([-0.2 0.9])
ylim([0 1])
numV = size(model.Mesh.Nodes,2);
title(['Set Hmax, Hmin. Num of node:',num2str(numV)])

% --------------------------------------
% Set scale to avoid numeric error
scale_factor = 100;
X = X*scale_factor;
Y = Y*scale_factor;

% Create polyshape
P= polyshape(X,Y);

% Constrained triangulation of the polyshape
TRI= triangulation(P);

% Create mesh from original triangulation
model2= createpde;
geometryFromMesh(model2,TRI.Points',TRI.ConnectivityList');

% --------------------------------------
% plot
nexttile
generateMesh(model2,'GeometricOrder','linear','Hmax',30, 'Hmin', 0.01);
pdemesh(model2)
xlim([-20 90])
ylim([0 100])
numV = size(model2.Mesh.Nodes,2);
title(['Try to fix. Num of node:',num2str(numV)])

●