%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-02-19 22:39:04
 % @LastEditTime: 2021-02-20 11:46:56

 % PARAMATERS: macro:[nely,nelx,nelz], micro:[microx^3], lx: macro-ele length 
 % opt: type of microstructure
%--------------------------
function [localx, globalx] = generateStruct(microx, nelx, nely, nelz, opt)

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

    switch opt         
        case 'interval'
%             x(microx/2 : end, :) = 1;
%             x(:, round(microx/2) : end, :) = 1;
            x(:) = 1;
            quartx = round(microx / 4);
            x(:, 1+quartx:microx-quartx, :) = 0;
%             x = 1-x;
            
            for i = 1:nelx    
                for j = 1:nely
                    for k = 1:nelz
                        [offsetx, offsety, offsetz] = micro_offset(i, j, k, microx);
                        globalx(offsety, offsetx, offsetz) = x;
                        localx{j,i,k} = x;
                    end
                end
            end
            
        case 'ellipse'    
            mina = 0.1; maxa = 0.4;
            minb = 0.1; maxb = 0.4;
            minc = 0.1; maxc = 0.4;
    
            a = maxa;
%             theta = 0;            
            for i = 1:nelx    
                b = maxb;
                for j = 1:nely
                    c = maxc;
                    for k = 1:nelz
                        x(:) = 0;

                        for ii = 1:microx
                            for jj = 1:microx
                                for kk = 1:microx
                                    ey = jj-microx/2;
                                    ex = ii-microx/2;
                                    ez = kk-microx/2; % TODO: ADD rotation
                                    if (ex/microx/a)^2 + (ey/microx/b)^2 + (ez/microx/c)^2 < 1
                                        x(jj,ii,kk) = 1;
                                    end
                                end
                            end
                        end
                        [offsetx, offsety, offsetz] = micro_offset(i, j, k, microx);
                        globalx(offsety, offsetx, offsetz) = x;
                        localx{j,i,k} = x;

                        c = c - (maxc-minc)/nelz;
                    end
                    b = b - (maxb-minb)/nely;
                end
                a = a - (maxa-mina)/nelx;
                % theta = theta + (pi/2)/nelx;
            end
            globalx = 1-globalx;            
            
        case 'bar'   
            mina = 0.1; maxa = 0.4;
            w = mina;
            for i = 1:nelx
                for j = 1:nely
                    for k = 1:nelz
                        x(:) = 0;

                        % FIXME: Z-axis don't need to be changed
                        x(1:round(w*microx), :, :) = 1;
                        x(round(microx-w*microx+1):microx, :, :) = 1;
                        x(:, 1:round(w*microx), :) = 1;
                        x(:, round(microx-w*microx+1):microx, :) = 1;
                        x(:, :, 1:round(w*microx)) = 1;
                        x(:, :, round(microx-w*microx+1):microx) = 1;

                        [offsetx, offsety, offsetz] = micro_offset(i, j, k, microx);
                        globalx(offsety, offsetx, offsetz) = x;
                        localx{j,i,k} = x;

                        w = w + (maxa-mina)/nelx/nely/nelz;
                    end
                end
            end
    end
    
%     for i = 1:nelz*microx
%         figure; imagesc(1-globalx(:,:,i)); axis equal;
%     end
end

% compute the index of each micro-ele in full-scale
function [offsetx, offsety, offsetz] = micro_offset(i, j, k, microx)
    offsetx = (i-1)*microx+1 : i*microx;
    offsety = (j-1)*microx+1 : j*microx;
    offsetz = (k-1)*microx+1 : k*microx;
end