MATLAB: Prevent nodes from seeing each other across a circle.

2dfem

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

Best Answer

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.
Related Question