Why does 'pdegplot' give the error "Error using checkSelfIntersections ... Edges overlap or intersect at non-end points" when trying to plot a geometry with circles in contact with each other

I am trying to create a geometry with circles in contact with each other. This is meant to simulate heat transfer between two pipes in physical contact. When I plot it with 'pdegplot', I get the following error.

model = createpde(1);

circle3 = struct('x', 0, 'y', -0.05, 'radius', 0.045);
circle4 = struct('x', 0, 'y', -0.05, 'radius', 0.05);
circle5 = struct('x', 0, 'y', 0, 'radius', 0.10);
circle6 = struct('x', 0, 'y', 0, 'radius', 0.11);

% pipes C4-C3 and C6-C5
C3 = create_csg_circle(circle3);
C4 = create_csg_circle(circle4);
C5 = create_csg_circle(circle5);
C6 = create_csg_circle(circle6);

% combine pipes
geometry_description = [C3, C4, C5, C6];
names_of_shapes      = char('C3','C4','C5','C6')';
set_formula          = 'C6-C5+C4-C3';

% create geometry
regions = decsg(geometry_description, set_formula, names_of_shapes);

figure
pdegplot(regions, 'EdgeLabels', 'on', 'FaceLabels', 'on');
axis equal

geometryFromEdges(model, regions);
generateMesh(model, 'Hmin', 0.0005);

figure
pdeplot(model);
axis equal

function csg = create_csg_circle(c)
CSG_CIRCLE = 1;
csg = [CSG_CIRCLE; c.x; c.y; c.radius];
end

>> error_script

Error using checkSelfIntersections (line 59)

Invalid geometry detected. Edges overlap or intersect at non-end points.

Error in initmesh (line 119)

checkSelfIntersections(g)

Error in pdegplot>plotTwoDGeometry (line 139)

[p1,~,t1]=initmesh(g,'init','on', 'MesherVersion','R2013a');

Error in pdegplot (line 102)

hh = plotTwoDGeometry(g, plotVertexLabels, plotEdgeLabels, plotSubLabels);

Error in tst2 (line 24)

pdegplot(regions, 'EdgeLabels', 'on', 'FaceLabels', 'on');

Answer 1

[MathWorks Support Team]

The reason for this error is that there are several issues with the geometry.

Note that 'pdegplot' is a red herring here. Even if you remove the call to 'pdegplot', you will see the same error when calling 'generateMesh'.

(1) In general, circles in a geometry will not be true circles (i.e. they are approximated with a discrete number of edges). The edges of the circles may then overlap when discretized. This is a likely culprit of the error.

(2) This geometry is challenging to mesh.

Because the angle between two circles that touch at a point is technically 0, the triangles in that region will have very small angles. A triangle with such a morbidly acute angle will not yield good results in the simulation.

This issue is the result of a limitation of most meshing tools and is industry-wide.

(3) This geometry would not simulate heat transfer between two pipes as expected.

For one, if the region of contact between the pipes in the simulation is a point, it will have zero area. This means no heat will be conducted between the pipes in the simulation.

In reality, the pipes will not be perfectly round. Imperfections will create a region of contact between the pipes with nonzero area.

The workaround is to move the circles apart slightly and add a tiny region that overlaps the two circles to model the contact between the pipes. This simplifies the geometry to be meshed and allows the simulation to more realistically model conduction between two pipes (see the revised script below).

model = createpde(1);

circle3 = struct('x', 0, 'y', -0.05, 'radius', 0.045);
circle4 = struct('x', 0, 'y', -0.05, 'radius', 0.05);
circle5 = struct('x', 0, 'y', 0, 'radius', 0.10);
circle6 = struct('x', 0, 'y', 0, 'radius', 0.11);

% pipes C4-C3 and C6-C5
C3 = create_csg_circle(circle3);
C4 = create_csg_circle(circle4);
C5 = create_csg_circle(circle5);
C6 = create_csg_circle(circle6);

% patch
R1 = [3; 4; -0.005; 0.005; 0.005; -0.005; -0.1001; -0.1001; -0.0997; -0.0997];

C3 = [C3;zeros(length(R1)-length(C3),1)];
C4 = [C4;zeros(length(R1)-length(C4),1)];
C5 = [C5;zeros(length(R1)-length(C5),1)];
C6 = [C6;zeros(length(R1)-length(C6),1)];

% combine pipes and patch
geometry_description = [C3, C4, C5, C6, R1];
names_of_shapes      = char('C3','C4','C5','C6','R1')';
set_formula          = 'C6-C5+C4-C3+R1';

% create geometry and remove junk inner boundaries
[regions,bt] = decsg(geometry_description, set_formula, names_of_shapes);
[domain,bt] = csgdel(regions,bt, [15, 16, 22, 23]);

figure
pdegplot(domain, 'EdgeLabels', 'on', 'FaceLabels', 'on');
axis equal

geometryFromEdges(model, domain);
generateMesh(model, 'Hmin', 0.0005);

figure
pdeplot(model);
axis equal

function csg = create_csg_circle(c)
CSG_CIRCLE = 1;
csg = [CSG_CIRCLE; c.x; c.y; c.radius];
end