MATLAB: Delaunay Triangulation removal of a box toward the center

delaunay triangulationMATLAB

I am having trouble removing triangulation from 2-5 on the x axis and 1.8 to 3.8 on the y axis. I managed to remove most it but I cannot figure out why excess triangulation was removed at the top left and bottom right. I have attachd my code.
A1 = 2;
B1 =1.8;
ElementSize1 = .33;
nx1 = ceil(A1/ElementSize1);
ny1 = ceil(B1/ElementSize1);
A2 = 2.6;
B2 =1.8;
ElementSize2 = .288;
nx2 = ceil(A2/ElementSize2);
ny2 = ceil(B2/ElementSize1);
% temperature at the bound1ary
TemperatureAtTop = 100;
TemperatureAtBottom = 0;
TemperatureAtLeft = 0;
TemperatureAtRight = 0;
%% Define nodes
% use the meshgrid and linspace functions to greate a rectangular grid of nodes
%8 meshgrids which make up the whole grid
[x1,y1] = meshgrid(linspace(0,A1,nx1),linspace(0,1.8,ny1+1));
[x2,y2] = meshgrid(linspace(2.2,4.8,nx2),linspace(0,1.8,ny2+1));
[x3,y3] = meshgrid(linspace(5,7,nx1),linspace(0,1.8,ny1+1));
[x4,y4] = meshgrid(linspace(5,7,nx1),linspace(2,3.6,ny2+1));
[x5,y5] = meshgrid(linspace(5,7,nx1),linspace(3.8,5.6,ny2+1));
[x6,y6] = meshgrid(linspace(2.2,4.8,nx2),linspace(3.8,5.6,ny2+1));
[x7,y7] = meshgrid(linspace(0,2,nx1),linspace(3.8,5.6,ny2+1));
[x8,y8] = meshgrid(linspace(0,2,nx1),linspace(2,3.6,ny2+1));
x=[x1 x2 x3 x4 x5 x6 x7 x8];
y=[y1 y2 y3 y4 y5 y6 y7 y8];
n = length(x(:)); % Total number of nodes
% Perform Dulanay triangulation to construct the elements from this set of
% nodes and build the matrix maping the element node index to the global node index.
ei2gi = delaunay(x,y); % Element index (EI) to global index (GI)
Ne = length(ei2gi); % Total number of elements
%% Identify boundary and free nodes
%
% Define the list of indicies for the fixed nodes and the list of their
% corresponding boundary temperatures. Use this to find the list of free
% nodes at which we are solving the temperature.
% Define the list of indices for all the nodes (these are the GI)
eta = 1:n;
etain=1:n;
% Build lists of indices for each of the 4 edges and the 4 inner edges
etaL = []; etaR = []; etaT = []; etaB = [];
etainL = []; etainR = []; etainT = []; etainB = [];
for i = eta,
if x(i) == 0, etaL(end+1) = i;
elseif x(i) == 7, etaR(end+1) = i;
else
if y(i) == 0, etaB(end+1) = i;
elseif y(i) == 5.6, etaT(end+1) = i;
end
end
end
for i= etain,
if x(i) ==2, etainL(end+1) = i;
elseif x(i)== 5 etainR(end+1)= i;
else
if y(i) == 1.8, etainB(end+1) = i;
elseif y(i) == 3.8, etainT(end+1) = i;
end
end
end
etain=[etainL etainR etainT etainB]
ElementsInHole = [];
for el = 1:Ne,
gi = ei2gi(el,:);
membership = ismember(gi,etain);
if sum(membership)== 3.0,
ElementsInHole(end+1) = el;
end
end
ElementsInDomain = setdiff(1:Ne,ElementsInHole);
% Update the list of elements and plot it to check it looks OK
ei2gi = ei2gi(ElementsInDomain,:);
Ne = length(ei2gi);
figure(1), trimesh(ei2gi,x,y), axis equal;

Best Answer

Were it my problem to solve, I'd try one of several things.
Simplest seems to be to start with a coarse quadtrilateration, thus 8 large rectangles. Split each rect into 2 triangles. Then refine each triangle as deeply as I wanted. 
[x,y] = meshgrid([0 2 5 7],[0 1.8 3.8 5.5]);
xy = [x(:),y(:)];
q = [1;2;3;5;7;9;10;11];
quads = [q,q+1,q+5,q+4];
% split each quad into a pair of triangles
tri = [quads(:,[1 2 3]);quads(:,[4 1 3])];
triplot(tri,xy(:,1),xy(:,2))
You should see this has the desired coarse geometry. Next, can you refine the triangulation? For example, split each triangle into 4 smaller, similar triangles? Use the midpoint of each edge. This is not too difficult, and I'd rather not do all of your work for you, so why not try it?
Here, for example, is the result of 3 such iterative refinements, splitting each triangle into 4 children 3 times.
This should not take much more than a couple of lines of code. So why not take a shot? (Ok, while I admit, that after looking at some code I wrote long ago to refine a triangulation, it took more than a couple of lines, but my code was far more sophisticated than a simple similar triangle splitting. Not that hard though, and you are halfway there now.)
Related Question