%--------------------------------------------------------------------------
% Mesh a L-shape using 4 noded Elements        
%--------------------------------------------------------------------------
function [nodes, eles, eleNum, nodeNum,plotIdx, iK, jK] = GenerateL(nelx, nely)

% dbstop if error
%     nelx = 20;
%     nely = 20;

    % left - rec
    [nodes1, eles1, eleNum1, nodeNum1] = GenerateMesh(nelx/2, nely, nelx/2, nely);
    % right - square
    [nodes2, eles2, eleNum2, nodeNum2] = GenerateMesh(nelx/2, nely/2, nelx/2, nely/2);

    nodes1(:,3) = 1:nodeNum1;
    borderNodes = nodes1(end-nely:end-nely/2,3); % get border nodes index for eles2

    % offset + nodes2
    nodes2(:,1) = nodes2(:,1)+nelx/2;
    nodes2(:,2) = nodes2(:,2)+nely/2;

    % eles2(1:nely/2+1) -> borderNodes
    [~,cols] = intersect(eles2,1:nely/2+1);

    % offset + eles2
    eles2 = eles2+nodeNum1-(nely/2+1);
    
    % eles2(1:nely/2+1) -> borderNodes
    tempeles2 = eles2(:);
    tempeles2(cols,1) = borderNodes;
    eles2 = reshape(tempeles2,eleNum2,4);    

    nodes = [nodes1;nodes2(nely/2+2:end,:)];
    eles = [eles1;eles2];
    eleNum = eleNum1+eleNum2;
    nodeNum = nodeNum1+nodeNum2-(nely/2+1);
    nodes(:,3) = 1:nodeNum;
    
    % compute the ele IDX of L shape, include the right-down empty square  
    plotIdx = 1:(nely*nelx/2); % left - rec
    for i = 0:(nelx/2-1)
        temp = (i*nely+1) : (nely/2+(i*nely));
        plotIdx = [plotIdx,temp+(nely*nelx/2)];
    end

    edofMat = zeros(eleNum, 8);
    for i = 1:eleNum
        for j = 1:4
            edofMat(i, [2*j-1,2*j]) = [2 * eles(i, j) - 1, 2 * eles(i, j)];
        end    
    end
    iK = reshape(kron(edofMat, ones(8, 1)).', 64 * eleNum, 1);
    jK = reshape(kron(edofMat, ones(1, 8)).', 64 * eleNum, 1);

    % test
%     scatter(nodes(:,1),nodes(:,2));
%     figure;
%     for i = 1:eleNum
%         node = eles(i,:);
%         scatter(nodes(node,1),nodes(node,2));hold on;
%     end
end