Prevent nodes from seeing each other across a circle.

Hello,

I am trying to create a mesh using a given set of coordinates, but I am having a slight issue. Basically, I want to triangulate the entire domain but exclude a cut-out circle inside of it. The cut-out circle is excluded from the given set of coordinates, but it seems that the nodes still see other nodes across the circle. I was wonering if anyone has any ideas on how to prevent the nodes from seeing each other across the cricle.

n = length(ax); 

% Perform Delaunay triangulation 
e2g = delaunay(ax,ay);        %element to global index & Delaunay function
Ne = length(ei2gi);             %total number of elements


%change number of elements to reflect domain elements
Ne = length(e2g);
%plot mesh of triangulation
figure(1)
trimesh(e2g,ax,ay)
axis equal; % Plot the mesh

[genevibes]

great question! would love to know this as well.

Answer 1

[John D'Errico]

Your problem is you are trying to use a delaunay triangulation tool to create the result. A general delaunay triangulation describes a CONVEX domain. A region with a hole inside is not convex by definition, so it will fail.

(It WOULD have helped if you provided the data, so I could easily give you some examples. Instead, I will need to cook some up. And, of course, I won't spend the time to make a very good example.) I'll even try constraining the edges of the circle to be edges in the tringulation.

nc = 31;
t = linspace(0,2*pi,nc)';
ps = polyshape(cos(t)+1,sin(t)+1); % used to initially exclude generated points in the circle

t(end) = [];
xyc = [cos(t),sin(t)] + 1; % A unit radius circle centered at [1,1]
tric = (1:nc-1)';
tric = [tric,tric+1];
tric(end,2) = 1;

xyr = rand(500,2)*6 - 3;
k = isinterior(ps,xyr(:,1),xyr(:,2));
xyr(k,:) = [];

xy = [xyc;xyr];
tri = delaunayTriangulation(xy,tric)

triplot(tri)
axis equal

So, did it work? Of course not. A Delaunay trinagulation describes a convex set. So even though I constrained the edges of the circle to be edges of the triangulation, the interior of the circle is still triangulated too.

Could I have done it differently? Well, yes, that is not difficult. For example, start with the delaunay triangulation, as above, wherein I constrained it to have edges on the circle. Then remove any triangles from the triangulation that have the property that ALL 3 vertices of the triangle lie on the circular boundary. This will be sufficient to create a triangulation with the desired property - an interior "circular" hole.

Thus, starting with the constrained Delaunay triangulation, do this:

triCL = tri.ConnectivityList;
triCL(all(triCL <= 30,2),:) = [];
triplot(triCL,tri.Points(:,1),tri.Points(:,2))
axis equal

While that was not so difficult, we can use a simpler scheme yet - an alpha shape, even though it lacks an explicit constraint option.

stri = alphaShape(xy(:,1),xy(:,2),.75,'HoleThreshhold',1);
plot(stri)
axis equal

This did work of couse, since an alpha shape is able to describe a non-convex domain. Note the difference of the two triangulations around the perimeter of the domain. A delaunay triangulation will sometimes tend to have many very long thin triangles around the perimeter. You can see them in the first two figures I show. This does depends of the way the points were generated, of course. However, an alpha shape will not have that failing.

[matlabpasta]

Thanks it worked! I didn't realize that the triangulation is only used for convex sets.

[John D'Errico]

I just edited my post, showing how to remove the triangles that cross the circular hole from the constrained delaunay triangulation. Either way will work. So you can choose your preference. In fact, an alphs shape also starts with a delaunay triangulation, then modifying it by deleting specific triangles.