Height = 0.5;
Width = 1;

Nx = 120; %Elements in x-directon.
Ny = 60; %Elements in y-direction.

DX  = 0.35;
DH = Height/10;
eps=1e-6;

coord = zeros(Nx+Ny,2);
nnodes = 0;
nodes = zeros(Ny+1,Nx+1);
for i=1:Nx+1;
    for j=1:Ny+1;
        
        xcoord = (i-1)/Nx*Width;
        ycoord = (j-1)/Ny*Height;
        
        
        
        if ( (xcoord > DX-eps) && (ycoord < DH + eps) )
            nodes(j,i) = 0;
        else
            nnodes = nnodes+1;
            coord(nnodes,1) = xcoord;
            coord(nnodes,2) = ycoord;
            nodes(j,i) = nnodes;
        end
        
    end
end

nelmH = Nx*Ny;
Enod = [];
el = 0;
for y = 1:size(nodes,1)-1;
    for x =1:size(nodes,2)-1;
        
        if ((nodes(y,x) == 0) || (nodes(y,x+1) == 0) || (nodes(y+1,x+1) == 0) || (nodes(y+1,x) == 0 ) )
            
        else
            el = el + 1;
            Enod(el,:) = [el, nodes(y,x), nodes(y,x+1), nodes(y+1,x+1), nodes(y+1,x)];
        end
    end
end

Edof = [Enod(:,1), Enod(:,2)*2-1, Enod(:,2)*2, Enod(:,3)*2-1,Enod(:,3)*2, Enod(:,4)*2-1, Enod(:,4)*2, Enod(:,5)*2-1, Enod(:,5)*2];
ndof = size(coord,1)*2;
nelm = size(Edof,1);
dof = [1:2:ndof-1; 2:2:ndof]';
[Ex,Ey] = coordxtr(Edof,coord,dof,4);

plotpar = [1 2 3];
figure(1), eldraw2(Ex,Ey,plotpar)

%If the number of elements on vertical side is even it will be exactly in
%middle. Otherwise on the node below the midpoint.
eps = 1e-5;

bc = [];
yc = Height;
xc = 0;
for i=1:ndof/2
    
    %LoadDof
    if (coord(i,1) > Width - eps) && (coord(i,1) < Width + eps)
        
        if ((coord(i,2)  < yc))
            yc = coord(i,2);
            dofLoad = i*2;
        end
        
    end
    
    
    %   BotMiddle
    if ((coord(i,2)  > 0- eps) && (coord(i,2)  < 0 + eps) )
        
        if coord(i,1) > xc
            xc = coord(i,1);
            dofBotMid = i*2-1;
        end
    end
    
    
    %Symmetry boundary
    if ((coord(i,2)  > 0- eps) && (coord(i,2)  < 0 + eps) )
        
        if  (coord(i,1)  < DX + eps)
            
            bc = [bc; i*2 0;];
        end
    end
    
    
    
end
   

bc = [bc; dofBotMid 0;];


F = zeros(ndof,1);
F(dofLoad) = 0.2e-3;
