Transfer your mesh from GMSH to Matlab

I believe examples are on of the best ways to learn, so I start from an example.

1- build your geometry in GMSH (you can easily do that by studying the first tutorial of the software)

2- Mesh the geometry by going to the Mesh tab and clicking on the 2D. This will build a triangular non structured mesh

Homework 8 solution

Homework 8 solution (click to download from Scribd)

HW#8

HW8 problem description (click to download from Scribd)

GMSH

Here is the link to a very useful free mesh generator and post processor software (GMSH) which is distributed under the terms of the GNU General Public License (GPL). http://geuz.org/gmsh/ There are versions available for Win, Mac and Linux.
You can find tutorials in the tutorial folder of the downloaded file.

Below is an example of unstructured mesh for a circular surface with a squared hole inside.

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

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

3 node triangular element

Below you may find a Matlab function to generate 3 node triangular elements. The output will not be symmetric, so I will try to write one for symmetric triangular elements soon.
Note that nx and ny are the number of rectangular elements in each direction.

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

% function [Edof,Ex,Ey]=trimesh(xcorner,ycorner,nx,ny)
%----------- Description ---------
% Generates triangular 3 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
%    3         4
%    o---------o
%    | +       |
%    |   +     |
%    |     +   | 
%    |       + |
%    o---------o
%    1         2
% ----------- 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)*(nx+1)+j;

            ul=dl+(nx+1); 

            dr=dl+1;

            ur=ul+1;

            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_ul=[(j-1)*(l/nx)+xcorner(1) i*(h/ny)+ycorner(1)];

            xy_ur=[j*(l/nx)+xcorner(1) i*(h/ny)+ycorner(1)];

            el=2*(i-1)*nx+2*j-1;

            Edof(el,:)=[el,2*dl-1,2*dl,2*dr-1,2*dr,2*ul-1,2*ul];

            Edof(el+1,:)=[el+1,2*dr-1,2*dr,2*ur-1,2*ur,2*ul-1,2*ul];

            Ex(el,:)=[xy_dl(1) xy_dr(1) xy_ul(1)];

            Ex(el+1,:)=[xy_dr(1) xy_ur(1) xy_ul(1)];

            Ey(el,:)=[xy_dl(2) xy_dr(2) xy_ul(2)];

            Ey(el+1,:)=[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

4 node quadrilateral element

Below you may find a Matlab function written by me to generate 4 node quadrilateral elements. It is very easy and straight forward.


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


% function [Edof,Ex,Ey]=quadmesh(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, 10 November 2011


l=xcorner(2)-xcorner(1);
h=ycorner(3)-ycorner(2);
for i=1:ny
        for j=1:nx
            dl=(i-1)*(nx+1)+j;
            ul=dl+(nx+1); 
            dr=dl+1;
            ur=ul+1;
            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_ul=[(j-1)*(l/nx)+xcorner(1) i*(h/ny)+ycorner(1)];
            xy_ur=[j*(l/nx)+xcorner(1) i*(h/ny)+ycorner(1)];
            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

Homework 7 solution

Homework 7 solution (click to download from Scribd)

HW#7

HW7 problem description (click to download from Scribd)