%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2021-02-25 15:13:41
 % @LastEditTime: 2021-02-26 19:43:34

 % Generate 3d hex mesh for FEM;
 %   y              4-----3
 %   |             /|    /|
 %   0---x        8-----7 2   
 %  /             |     |/
 % z              5-----6
 % 1. COORDINATES:  first y(1-0) -> x(0-1) -> z(0-1)
 % 2. CONNECTIVITY: 1->8 TODO:
 % 3. PARAMETERS:   nelx * nely * nelz macro-elements with microx^3 micro-elements, 
 %                  size of each element is lx^3
 % 4. RETRUN:       nodes: nodes coordinates as (cx, cy, cz), size of [nodes_num, 3]
%--------------------------
function nodes = mesh3D(nelx, nely, nelz, microx, num_b, lx, opt_scale)

    if opt_scale==1
        %% 1-scale
        
        dx = lx / microx;
        globaly = nely * microx;
        globalx = nelx * microx;
        globalz = nelz * microx;
        % Number of points
        npx = globalx + 1;
        npy = globaly + 1;
        npz = globalz + 1;
        % Discretizing
        nx = linspace(0, dx*globalx, npx);
        ny = linspace(0, dx*globaly, npy);
        nz = linspace(0, dx*globalz, npz);
        [gridx, gridy, gridz] = meshgrid(nx, fliplr(ny), nz);       
%         [gridx, gridy, gridz] = meshgrid(nx, (ny), nz);       

        nodes = [gridx(:), gridy(:), gridz(:)];

    else
        %% 2-scale
        nodes = nodes_2scale(nelx, nely, nelz, num_b, lx);
    end
end

function nodes = nodes_2scale(nelx, nely, nelz, num_b, lx)
    
    new_nelx = nelx * (num_b + 1);
    new_nely = nely * (num_b + 1);
    new_nelz = nelz * (num_b + 1);
    Num_node = (1+new_nelx) * (new_nely+1) * (new_nelz+1);
    nodenrs = reshape(1:Num_node, 1+new_nely, 1+new_nelx, 1+new_nelz);

    dx = lx / (num_b+1);
    nx = linspace(0, lx*nelx, new_nelx+1);
    ny = linspace(0, lx*nely, new_nely+1);
    nz = linspace(0, lx*nelz, new_nelz+1);
    [gridx, gridy, gridz] = meshgrid(nx, fliplr(ny), nz);
    gridx = gridx(:);
    gridy = gridy(:);
    gridz = gridz(:);
    
    [up, front, down, back, left, right, dofs_num] = surface_dofs_id(num_b+1);
    surface_nodes = zeros(nelx, nely, nelz, dofs_num/3);
    
    for i = 1 : nelx
        for j = 1 : nely
            for k = 1 : nelz
                nodes_ele = nodenrs((j-1)*(num_b+1)+1 : j*(num_b+1)+1, ...
                                    (i-1)*(num_b+1)+1 : i*(num_b+1)+1, ...
                                    (k-1)*(num_b+1)+1 : k*(num_b+1)+1);
                inner = nodes_ele(2:end-1, 2:end-1, 2:end-1);
                surface_nodes(i,j,k,:) = setdiff(nodes_ele(:), inner(:));
            end
        end
    end
    global_nodes = unique(surface_nodes(:)); % remove interior nodes of every macro-ele

    nodes = [gridx(global_nodes), gridy(global_nodes), gridz(global_nodes)];
end

function nodes_2scale_way2()
    dx = lx / (num_b+1);
    new_nelx = nelx * (num_b + 1);
    new_nely = nely * (num_b + 1);
    new_nelz = nelz * (num_b + 1);
    
    % there are two nodes distribution of y0x, n is orignal node & p is bezier pnt
    % 1. n-p-p-n    2. p-p-p-p
    %    p-p-p-p       p     p
    %    p-p-p-p       p     p
    %    n-p-p-n       p-p-p-p
    xy_nodes1_num = (new_nelx+1) * (new_nely+1); 
    xy_nodes2_num = (nelx+1)*(nely+1) + num_b*nely*(nelx+1) + num_b*nelx*(nely+1); 

    nodes = zeros(xy_nodes1_num*(nelz+1) + xy_nodes2_num*nelz*num_b, 3);
    nodes_y0x = mesh_y0x(lx, nelx, nely, num_b);

    nz  = linspace(0, dx*new_nelz, new_nelz+1);
    nid = 0; 
    % figure;
    for k = 1 : nelz*(num_b+1) + 1
        % now nid is the last node-id of this y0x
        if mod(k-1, num_b+1)==0
            nx = linspace(0, lx*nelx, new_nelx+1);
            ny = linspace(0, lx*nely, new_nely+1);                
            [gridx, gridy, gridz] = meshgrid(nx, fliplr(ny), nz(k));

            nodes(nid+1:nid+xy_nodes1_num, :) = [gridx(:), gridy(:), gridz(:)];
            % scatter3(nodes(nid+1:nid+xy_nodes1_num,1),nodes(nid+1:nid+xy_nodes1_num,2),nodes(nid+1:nid+xy_nodes1_num,3),'filled'); axis equal; hold on
            nid = nid + xy_nodes1_num;                                
        else
            gridz = repmat(nz(k), size(nodes_y0x(:,1)));
            nodes(nid+1:nid+xy_nodes2_num, :) = [nodes_y0x, gridz(:)];
            % scatter3(nodes(nid+1:nid+xy_nodes2_num,1),nodes(nid+1:nid+xy_nodes2_num,2),nodes(nid+1:nid+xy_nodes2_num,3),'filled'); axis equal;
            nid = nid + xy_nodes2_num;                                
        end
    end
%         figure; scatter3(nodes(:,1),nodes(:,2),nodes(:,3),'filled'); axis equal;
end

% mesh nodes distribution-2 of y0x
% addNum: add control points num, could be zero
% nodes order: col by col
function [nodes] = mesh_y0x(lx, nelx, nely, addNum)
    nodeNum = (nelx+1)*(nely+1) + addNum*nely*(nelx+1) + addNum*nelx*(nely+1);
    nodes = zeros(nodeNum, 2);

    left_pnum   = nely+1 + addNum*nely;
    center_pnum = (nely+1) * addNum;
    dx = lx / (addNum+1);

    for i = 0:nelx
        % left col
        cy = linspace(lx*nely, 0, left_pnum)';
        cx = repmat(lx * i, left_pnum, 1);

        idx = (1:left_pnum) + i * (left_pnum + center_pnum);
        nodes(idx, :) = [cx, cy];

        if i>nelx-1
            break;
        end

        % center col
        cy = linspace(lx*nely, 0, nely+1)';
        cy = repmat(cy, addNum, 1);
        cx = linspace(lx * i+dx, lx * (i+1)-dx, addNum);
        cx = repmat(cx, nely+1, 1);

        idx = (1:center_pnum)+left_pnum + i * (left_pnum + center_pnum);
        nodes(idx, :) = [cx(:), cy];
        
        % NOTE: no right col, as it is next left col
    end    
end