4 node quadrilateral element for non symmetrical rectangles

Below you may find the Matlab function for 4 node quadrilateral element which is capable of meshing non symmetric rectangles. I will try to do the same for triangles.


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


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


ld=abs(xcorner(2)-xcorner(1));
lu=abs(xcorner(3)-xcorner(4));
hr=abs(ycorner(3)-ycorner(2));
hl=abs(ycorner(4)-ycorner(1));
for i=1:ny
        for j=1:nx
            dl=(i-1)*(nx+1)+j;
            ul=dl+(nx+1); 
            dr=dl+1;
            ur=ul+1;
            l1=(i-1)*(lu-ld)/ny+ld; l2=i*(lu-ld)/ny+ld; 
            h1=(j-1)*(hr-hl)/nx+hl; h2=j*(hr-hl)/nx+hl;
            x_cd=(i-1)*(xcorner(4)-xcorner(1))/ny+xcorner(1);
            x_cu=i*(xcorner(4)-xcorner(1))/ny+xcorner(1);
            y_cl=(j-1)*(ycorner(2)-ycorner(1))/nx+ycorner(1);
            y_cr=j*(ycorner(2)-ycorner(1))/nx+ycorner(1);
            xy_dl=[(j-1)*(l1/nx)+x_cd (i-1)*(h1/ny)+y_cl];
            xy_dr=[j*(l1/nx)+x_cd     (i-1)*(h2/ny)+y_cr];
            xy_ul=[(j-1)*(l2/nx)+x_cu  i*(h1/ny)+y_cl];
            xy_ur=[j*(l2/nx)+x_cu      i*(h2/ny)+y_cr];
            el=(i-1)*nx+j;
            Edof(el,:)=[el,2*dl-1,2*dl,2*dr-1,2*dr,2*ur-1,2*ur,2*ul-1,2*ul];
            Ex(el,:)=[xy_dl(1) xy_dr(1) xy_ur(1) xy_ul(1)];
            Ey(el,:)=[xy_dl(2) xy_dr(2) xy_ur(2) xy_ul(2)];
        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