3 node symmetric triangular elements

Below you may find a Matlab function to generate 3 node symmetric triangular elements as I promised. 
Note that nx and ny are the number of rectangular elements in each direction.



function [Edof,Ex,Ey]=symtrimesh(xcorner,ycorner,nx,ny)

% function [Edof,Ex,Ey]=symtrimesh(xcorner,ycorner,nx,ny)
%----------- Description ---------
% Generates triangular 3 node symmetric elements
% ------------ INPUT -------------
% xcorner=[lower left, lower right, upper right, upper left] 
% ycorner=[lower left, lower right, upper right, upper left]
% nx= No. of rectangulat elements in x direction
% ny= No. of rectangulat elements in y direction
%    7         8         9        
%    o---------o---------o
%    | +       |       + |
%    |   +     |     +   |
%    |     +   |   +     |
%    |       + | +       |
%   4o---------5---------o6
%    |       + | +       |
%    |     +   |   +     |  
%    |   +     |     +   |
%    | +       |       + |
%    o---------o---------o
%    1         2         3
% ----------- OUTPUT -------------
% Ex= element x coordinates matrix [nrelem X 3]
% Ey= element y coordinates matrix [nrelem X 3]
% Edof= element DOF [nrelem X 7]
% Written by 
% Shah, 11 November 2011

l=xcorner(2)-xcorner(1);

h=ycorner(3)-ycorner(2);

for i=1:ny
        for j=1:nx

            dl=(i-1)*(2*nx+1)+2*j-1;
            dm=dl+1;
            dr=dl+2;
            ul=dl+(2*nx+1); 
            um=ul+1;
            ur=ul+2;
            xy_dl=[(j-1)*(l/nx)+xcorner(1) (i-1)*(h/ny)+ycorner(1)];
            xy_dr=[j*(l/nx)+xcorner(1) (i-1)*(h/ny)+ycorner(1)];
            xy_dm=(xy_dl+xy_dr)./2;
            xy_ul=[(j-1)*(l/nx)+xcorner(1) i*(h/ny)+ycorner(1)];
            xy_ur=[j*(l/nx)+xcorner(1) i*(h/ny)+ycorner(1)];
            xy_um=(xy_ul+xy_ur)./2;
            el=4*nx*(i-1)+4*(j-1)+1;
            if floor(i/2)-i/2~=0
                Edof(el,:)=[el,2*dl-1,2*dl,2*um-1,2*um,2*ul-1,2*ul];
                Edof(el+1,:)=[el+1,2*dl-1,2*dl,2*dm-1,2*dm,2*um-1,2*um];
                Edof(el+2,:)=[el+2,2*dm-1,2*dm,2*dr-1,2*dr,2*um-1,2*um];
                Edof(el+3,:)=[el+3,2*dr-1,2*dr,2*ur-1,2*ur,2*um-1,2*um];
                Ex(el,:)=[xy_dl(1) xy_um(1) xy_ul(1)];
                Ex(el+1,:)=[xy_dl(1) xy_dm(1) xy_um(1)];
                Ex(el+2,:)=[xy_dm(1) xy_dr(1) xy_um(1)];
                Ex(el+3,:)=[xy_dr(1) xy_ur(1) xy_um(1)];
                Ey(el,:)=[xy_dl(2) xy_um(2) xy_ul(2)];
                Ey(el+1,:)=[xy_dl(2) xy_dm(2) xy_um(2)];
                Ey(el+2,:)=[xy_dm(2) xy_dr(2) xy_um(2)];
                Ey(el+3,:)=[xy_dr(2) xy_ur(2) xy_um(2)];
            else
                Edof(el,:)=[el,2*dl-1,2*dl,2*dm-1,2*dm,2*ul-1,2*ul];
                Edof(el+1,:)=[el+1,2*dm-1,2*dm,2*um-1,2*um,2*ul-1,2*ul];
                Edof(el+2,:)=[el+2,2*dm-1,2*dm,2*ur-1,2*ur,2*um-1,2*um];
                Edof(el+3,:)=[el+3,2*dm-1,2*dm,2*dr-1,2*dr,2*ur-1,2*ur];
                Ex(el,:)=[xy_dl(1) xy_dm(1) xy_ul(1)];
                Ex(el+1,:)=[xy_dm(1) xy_um(1) xy_ul(1)];
                Ex(el+2,:)=[xy_dm(1) xy_ur(1) xy_um(1)];
                Ex(el+3,:)=[xy_dm(1) xy_dr(1) xy_ur(1)];
                Ey(el,:)=[xy_dl(2) xy_dm(2) xy_ul(2)];
                Ey(el+1,:)=[xy_dm(2) xy_um(2) xy_ul(2)];
                Ey(el+2,:)=[xy_dm(2) xy_ur(2) xy_um(2)];
                Ey(el+3,:)=[xy_dm(2) xy_dr(2) xy_ur(2)];       

             end            
        end
end

figure(1);

hold on;
for i=1:length(Edof(:,1));

    plot(Ex(i,[1:end 1]),Ey(i,[1:end 1]));

end
axis equal