%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2020-10-29 15:04:15
 % @LastEditTime: 2021-06-17 11:29:39
%--------------------------
function [localx, globalx] = generateStruct(microx, microy, nelx, nely, lx, ly, opt, scalar)

    [nodes, eles, ~, ~] = GenerateMesh(nelx*lx, nely*ly, nelx, nely);
    [globalOffsetx,globalOffsety] = offsetSIMP(eles,nodes, microx, microy, nelx, nely);

    x = zeros(microy, microx);
    globalx = zeros(nely*microy, nelx*microx);
    localx = cell(nely, nelx);

    mina = microx/10; maxa = round((microx-microx/4)/2);
    minb = microx/10; maxb = round((microx-microx/4)/2);  
%     mina = microx/10; maxa = round((microx-microx/3)/2);
%     minb = microx/10; maxb = round((microx-microx/3)/2); 

    switch opt
        case 'interval2'
%             globalx(1:microx, round(microx/3):round(microx*2/3)) = 1;
%             globalx(microx+round(microx/3):microx+round(microx*2/3), :) = 1;
%             globalx(1:microx, microx+round(microx/3):microx+round(microx*2/3)) = 1;

            globalx(1:microx, round(microx/3):round(microx*2/3)) = 1;   
            globalx(microx+round(microx/3):microx+round(microx*2/3), :) = 1;
            globalx(1:microx, microx+round(microx/3):microx+round(microx*2/3)) = 1;
            globalx = rot90(globalx);
            globalx = fliplr(globalx);
            globalx = 1-globalx;
            
        case 'checkboard2'  
            tmp_interval = round(microx/scalar);
            x(:) = 0;
            for i = 1:2:scalar
                x((i-1)*tmp_interval+1 : i*tmp_interval, :) = 1;
            end
            for i = 1:2:scalar
                x(:,(i-1)*tmp_interval+1 : i*tmp_interval) = 1;
            end
            for i = 1:2:scalar
                for j = 1:2:scalar
                    x((j-1)*tmp_interval+1 : j*tmp_interval, (i-1)*tmp_interval+1 : i*tmp_interval) = 0;
                end
            end
             figure; imagesc(x); axis equal; axis off; colormap(gray);
            
            for i = 1:nelx
                for j = 1:nely                                      
                    globalx(globalOffsety(j,i,:), globalOffsetx(j,i,:)) = 1-x;
                    localx{j,i} = 1-x;
                end
            end
            
        case 'checkboard'
            for i = 1:nelx
                for j = 1:nely
                    x(:) = 0;
                    x(1:round(microx/2), 1:round(microx/2)) = 1;
                    x(round(microx/2)+1:end, round(microx/2)+1:end) = 1;
                    
                    globalx(globalOffsety(j,i,:), globalOffsetx(j,i,:)) = 1-x;
                    localx{j,i} = 1-x;
                end
            end
            
        case 'multi_shape'
            a = maxa;
            theta = 0;            
            b = maxb;            
            i = 1; j = 1;
            x(:) = 0;
            for ii = 1:microx
                for jj = 1:microy
                    ey = jj-microy/2;
                    ex = ii-microx/2;
                    if (ey*sin(theta)+ex*cos(theta))^2/a^2 + (ey*cos(theta)-ex*sin(theta))^2/b^2<1
                        x(jj,ii) = 1;
                    end
                end
            end
            globalx(globalOffsety(j,i,:), globalOffsetx(j,i,:)) = x;
            localx{j,i} = x;
            
            i = 2; j = 1;
            w = 1/5;
            x(:) = 0;
            x(1:round(w*microy), :) = 1;
            x(round(microy-w*microy+1):microy, :) = 1;
            x(:, 1:round(w*microx)) = 1;
            x(:, round(microy-w*microx+1):microy) = 1;
            x = 1-x;
            globalx(globalOffsety(j,i,:), globalOffsetx(j,i,:)) = x;
            localx{j,i} = x;
            
            i = 2; j = 2;
            w = 1/3;
            x(:) = 0;
            x(1:round(w*microy), :) = 1;
            x(round(microy-w*microy+1):microy, :) = 1;
            x(:, 1:round(w*microx)) = 1;
            x(:, round(microy-w*microx+1):microy) = 1;
            x = 1-x;
            globalx(globalOffsety(j,i,:), globalOffsetx(j,i,:)) = x;
            localx{j,i} = x;

            a = mina;
            theta = pi/4;            
            b = maxb;            
            i = 1; j = 2;
            x(:) = 0;
            for ii = 1:microx
                for jj = 1:microy
                    ey = jj-microy/2;
                    ex = ii-microx/2;
                    if (ey*sin(theta)+ex*cos(theta))^2/a^2 + (ey*cos(theta)-ex*sin(theta))^2/b^2<1
                        x(jj,ii) = 1;
                    end
                end
            end
            globalx(globalOffsety(j,i,:), globalOffsetx(j,i,:)) = x;
            localx{j,i} = x;

            globalx = 1-globalx; 
        case 'diamond'
            maxa = round((microx-microx/5)/2);
            mina = microx/15;
            a = maxa;
            [gridx, gridy] = meshgrid(-microx/2+0.5:microx/2-0.5, -microx/2+0.5:microx/2-0.5);
            gridy = flipud(gridy);
            
            for i = 1:nelx
                for j = 1:nely
                    P = [0, a;
                        a, 0;
                        0, -a;
                        -a, 0];
                    pgon = polyshape(P);

                    x = isinterior(pgon, gridx(:), gridy(:));
                    x = reshape(1-x, microx, microx);

                    globalx(globalOffsety(j,i,:), globalOffsetx(j,i,:)) = x;
                    localx{j,i} = x;
%                     a = a - (maxa-mina)/nelx;
                end
            end
            
        case 'interval_dense'
            for i = 1:nelx
                for j = 1:nely
                    x(:) = 0;
                    x(:, 1:round(microx/4)) = 1;
                    x(:, round(microx/2)+1:end-round(microx/4)) = 1;
                    
                    globalx(globalOffsety(j,i,:), globalOffsetx(j,i,:)) = x;
                    localx{j,i} = x;
                end
            end

        case 'interval'
            for i = 1:nelx
                for j = 1:nely
                    x(:) = 0;
                    x(:, 1:round(microx/3)) = 1;
                    x(:, end-round(microx/3):end) = 1;
                    
                    globalx(globalOffsety(j,i,:), globalOffsetx(j,i,:)) = x;
                    localx{j,i} = x;
                end
            end
            
        case 'ellipse2'
            a = maxa;
            theta = 0;            
            for i = 1:nelx    
                b = maxb;
                for j = 1:nely    
                    x(:) = 0;

                    for ii = 1:microx
                        for jj = 1:microy
                            ey = jj-microy/2;
                            ex = ii-microx/2;
                            if (ey*sin(theta)+ex*cos(theta))^2/a^2 + (ey*cos(theta)-ex*sin(theta))^2/b^2<1
                                x(jj,ii) = 1;
                            end
                        end
                    end
                    globalx(globalOffsety(j,i,:), globalOffsetx(j,i,:)) = x;
                    localx{j,i} = 1-x;
                end
            end
            globalx = 1-globalx;  

        case 'ellipse_test'
            % for section of testing shape functions
            a = maxa;
            theta = 0;            
            for i = 1:nelx    
                b = maxb;
                for j = 1:nely    
                    x(:) = 1;                                        

                    for ii = 1:microx
                        for jj = 1:microy
                            ey = jj-microy/2;
                            ex = ii-microx/2;
                            if (ey*sin(theta)+ex*cos(theta))^2/a^2 + (ey*cos(theta)-ex*sin(theta))^2/b^2<1
                                x(jj,ii) = 0;
                            end
                        end
                    end
                    
                    if i==1 && j==1
%                         m0 = 0;
                    else
                        x = x*2;
                    end 
%                     x = x*m0;
%                     if i==1 && j==1
%                     else
%                         x(:) = 0;
%                     end
                    globalx(globalOffsety(j,i,:), globalOffsetx(j,i,:)) = x;
                    localx{j,i} = x;                                        
                end
            end
%             globalx = 1-globalx;  

        case 'ellipse_MBB'
            a = maxa;
            theta = 0;            
            for i = 1:nelx    
                b = maxb;
                for j = 1:nely    
                    x(:) = 0;

                    for ii = 1:microx
                        for jj = 1:microy
                            ey = jj-microy/2;
                            ex = ii-microx/2;
                            if (ey*sin(theta)+ex*cos(theta))^2/a^2 + (ey*cos(theta)-ex*sin(theta))^2/b^2<1
                                x(jj,ii) = 1;
                            end
                        end
                    end
                    globalx(globalOffsety(j,i,:), globalOffsetx(j,i,:)) = x;
                    localx{j,i} = 1-x;

                    b = b - (maxb-minb)/nely;
                end
                a = a - (maxa-mina)/nelx;
                theta = theta + (pi/2)/(nelx-1);
            end
            globalx = 1-globalx;       
            
        case 'ellipse_bending'
            a = maxa;
            theta = 0;            
            for i = 1:nelx    
                b = maxb;
                for j = 1:nely    
                    x(:) = 0;

                    for ii = 1:microx
                        for jj = 1:microy
                            ey = jj-microy/2;
                            ex = ii-microx/2;
                            if (ey*sin(theta)+ex*cos(theta))^2/a^2 + (ey*cos(theta)-ex*sin(theta))^2/b^2<1
                                x(jj,ii) = 1;
                            end
                        end
                    end
                    globalx(globalOffsety(j,i,:), globalOffsetx(j,i,:)) = x;
                    localx{j,i} = 1-x;
                end
                a = a - (maxa-mina)/nelx;
                theta = theta + (pi/2)/(nelx-1);
            end
            globalx = 1-globalx;   

        case 'ellipse'                      
            a = maxa;
            theta = 0;            
            for i = 1:nelx    
                b = maxb;
                for j = 1:nely    
                    x(:) = 0;

                    for ii = 1:microx
                        for jj = 1:microy
                            ey = jj-microy/2;
                            ex = ii-microx/2;
                            if (ey*sin(theta)+ex*cos(theta))^2/a^2 + (ey*cos(theta)-ex*sin(theta))^2/b^2<1
                                x(jj,ii) = 1;
                            end
                        end
                    end
                    globalx(globalOffsety(j,i,:), globalOffsetx(j,i,:)) = x;
                    localx{j,i} = 1-x;

                    b = b - (maxb-minb)/nely;
                end
                a = a - (maxa-mina)/nelx;
                theta = theta + (pi/2)/(nelx-1);
            end
            globalx = 1-globalx;            
            
        case 'bar'        
            w = mina/microx;
            for i = 1:nelx
                for j = 1:nely
                    x(:) = 0;

                    x(1:round(w*microy), :) = 1;
                    x(round(microy-w*microy+1):microy, :) = 1;
                    x(:, 1:round(w*microx)) = 1;
                    x(:, round(microy-w*microx+1):microy) = 1;

                    globalx(globalOffsety(j,i,:), globalOffsetx(j,i,:)) = x;
                    localx{j,i} = x;

                    w = w + (maxa-mina)/microx/nelx/nely;
                end
            end
        case 'bar-size'
            w = mina/microx;
            w = w + (maxa-mina)/microx/2/2 * scalar;
            for i = 1:nelx
                for j = 1:nely
                    x(:) = 0;

                    x(1:round(w*microy), :) = 1;
                    x(round(microy-w*microy+1):microy, :) = 1;
                    x(:, 1:round(w*microx)) = 1;
                    x(:, round(microy-w*microx+1):microy) = 1;

                    globalx(globalOffsety(j,i,:), globalOffsetx(j,i,:)) = x;
                    localx{j,i} = x;
                end
            end
    end

%     imshow(1-globalx);

end