%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2020-10-29 15:04:15
 % @LastEditTime: 2022-01-11 17:21:56

 % NOTE: nely & nelx is globaly & globalx
%--------------------------
function [freedofs,fixeddofs,F,loadnid] = designDomain(opt,nelx,nely,vF, microx, x_offset)
  
    alldofs = 2*(nelx+1)*(nely+1);
    
    switch opt
    case 'bending'
        fixeddofs = [1, 2*(nely+1), 2*(nely+1)*(nelx+1)];        

        % mimicking contact force
        [F,loadnid] = mimicking_contact_F(vF, nelx, nely, microx, x_offset);

        % 3-point concentrate force
        loadnid = [(nely+1) * (microx*x_offset+1) * 2, ...
                    (nely+1) * ((nelx+1+1)/2-1) * 2 + 2, ...
                    (nely+1) * (nelx+1-microx*x_offset) * 2];
        % vF = [vF/2, -vF, vF/2];

    case 'MBB'
        fixeddofs = [1:2:2*(nely+1),2*(nely+1)*(nelx+1)];
        loadnid = 2;
        vF = -vF;
        
    case 'cantilever'
        fixeddofs = 1:2*(nely+1);        
        loadnid = 2*(nely+1)*(nelx+1);
        vF = -vF;
        
    case 'cant2'
        fixeddofs = 1:2*(nely+1);        
        loadnid = [2*(nely+1)*(nelx+1)-1, 2*(nely+1)*(nelx+1)];
        vF = -vF;

    case 'Lshape'
%         fixeddofs = [2*(nely+1):2*(nely+1):(2*(nely+1)*nelx/2),...
%                     (2*(nely+1)-1):2*(nely+1):(2*(nely+1)*nelx/2-1)];
%         loadnid = 2*nelx*(nely+1)+nely+2;
        fixeddofs = [2:2*(nely+1):(2*(nely+1)*nelx/2),1:2*(nely+1):(2*(nely+1)*nelx/2-1)];
        loadnid = 2*nelx*(nely+1)+nely+2;       
    
    case 'michell'
        fixeddofs = 1:2*(nely+1);
        loadnid = [((nely+1)*(nelx+1)-nely/2)*2, ((nely+1)*(nelx+1)-nely/2)*2-1];
        vF = [-vF,vF]*sqrt(2)/2;
        
    case 'miche2'
        fixeddofs = 1:2*(nely+1);
        loadnid = (nelx*(nely+1)+1) : (nelx+1)*(nely+1);
        loadnid = [loadnid*2-1];   
        vF = vF/length(loadnid(:));
        
    case 'bridge1'
        loadnid = [2*(nely+1)];
        vF = -vF;
        fixeddofs = union(1:2:2*(nely+1),[2*(nelx+1)*(nely+1),2*(nelx+1)*(nely+1)-1]);
        
    case 'bridge2'
        loadnid = [2*(nely+1)*(nelx/2+1),2*(nely+1)];
        vF = [-vF, -2*vF];
        fixeddofs = union(1:2:2*(nely+1),[2*(nelx+1)*(nely+1),2*(nelx+1)*(nely+1)-1]);

    case 'internal'
        f1 = 1:2:2*(nely+1)-1; % left
        f2 = 2:2*(nely+1):2+2*(nely+1)*nelx; % up
        f3 = 1+2*(nely+1)*nelx:2:2*(nely+1)*(nelx+1)-1; % right
        f4 = 2*(nely+1):2*nely+2:2*(nelx+1)*(nely+1); % down
        loadnid = [f2,f3];
        vF = repmat(vF*-1, size(loadnid));
        
        fixeddofs = [1:2:2*(nely+1), 2*(nely+1):2*(nely+1):2*(nely+1)*(nelx+1)];
    
    case 'two-loads'
%         vF = 2*vF;
        fixeddofs = union([1:2:2*(nely+1)],2*(nely+1)*(nelx+1)-[1:1:2*(nely+1)]+1);    
        loadnid = [2,2*(nely+1)];
        vF = [-vF,vF];
        
%     case 'bending'
%         fixeddofs = [1:2:2*(nely+1),2*(nely+1)*(nelx+1)];
%         loadnid = [2*(nely+1), 2*(nely+1)*(nelx+1), 2*(nely+1)*(nelx)+1];
%         vF = [vF, -vF, -vF];
        
    case 'stretch'
        fixeddofs = 1:2*(nely+1);
        loadnid = 2*(nely+1)*nelx+1:2:alldofs;
    end
    
    freedofs = setdiff(1:alldofs,fixeddofs);
    
    if ~strcmp(opt, 'bending')
        F = sparse(loadnid,1,vF,alldofs,1);
    end

%     % nonregular mesh
%     if nargin > 4
%         [fixeddofs, freedofs, F] = FEA_init_nonregular(nodes, fixeddofs, loadnid,vF);
%     end

%     figure; scatter(nodes(:,1),nodes(:,2)); axis equal; hold on
%     scatter(nodes(round(fixeddofs/2), 1),nodes(round(fixeddofs/2), 2), 'filled'); axis equal; hold on
%     scatter(nodes(loadnid/2, 1),nodes(loadnid/2, 2), 'filled'); axis equal; hold on
%     scatter(nodes(loadnid, 1),nodes(loadnid, 2), 'filled'); axis equal; hold on
end

% only located on that macro-nodes
function [F,loadnid] = mimicking_contact_F(vF, globalx, globaly, microx, x_offset)
    nelx = globalx / microx;
    nely = globaly / microx;

    fj1 = nely+1;
    fi1 = x_offset+1;
    fv1 = vF/2;
    F1 = neighbors_F(fi1, fj1, fv1, 1, nelx, nely);
    
    fj2 = 1;
    fi2 = (nelx+1+1) / 2;
    fv2 = -vF;
    F2  = neighbors_F(fi2, fj2, fv2, 1, nelx, nely);

    fj3 = nely+1;
    fi3 = nelx+1 - x_offset;
    fv3 = vF/2;
    F3  = neighbors_F(fi3, fj3, fv3, 1, nelx, nely);

    F = F1 + F2 + F3;
    F_macro = reshape(F', [], 1);
    
%     [nodes, eles, eleNum, nodeNum] = GenerateMesh(nelx, nely, nelx, nely);
%     figure; scatter(nodes(:,1),nodes(:,2)); axis equal; hold on;
%     scatter(nodes(loadnid/2, 1),nodes(loadnid/2, 2), 'filled'); axis equal; hold on

    loadnid = [];

    % now F is only on macro-nodes, expand to full-scale
    [~,~,edofMat_ma,~]   = forAssemble(nelx, nely);
    F_full = zeros((globalx+1)*(globaly+1), 2);
    
    for i = 1:nelx
        for j = 1:nely
            edof = edofMat_ma(j+(i-1)*nely, :);
            Fe   = reshape(F_macro(edof), 2, [])';

            % the macro-nodes corresponding to micro-nodes: [x,1],[x,x],[1,x],[1,1]            
            micro_nodes = [microx+1, 1;
                            microx+1, microx+1;
                            1, microx+1;
                            1, 1];

            for im = 1:4
                jj = micro_nodes(im, 1);
                ii = micro_nodes(im, 2);

                fulli = (i-1)*microx + ii;
                fullj = (j-1)*microx + jj;
                full_nid = (fulli - 1) * (globaly+1) + fullj;
                
                F_full(full_nid, :) = Fe(im, :);

%                 loadnid = [loadnid, full_nid];
            end
        end
    end
    F = reshape(F_full', [], 1);
end

% using pressure to mimicking contact force
function F = mimicking_contact_F2(vF, nelx, nely, microx)
    
    fj1 = nely;
    fi1 = microx;
    fv1 = vF/2;
    F1  = neighbors_F(fi1, fj1, fv1, microx, nelx, nely);
    
    fj2 = 1;
    fi2 = (nelx+1+1) / 2;
    fv2 = -vF;
    F2  = neighbors_F(fi2, fj2, fv2, microx, nelx, nely);

    fj3 = nely;
    fi3 = nelx - microx;
    fv3 = vF/2;
    F3  = neighbors_F(fi3, fj3, fv3, microx, nelx, nely);

    F = F1 + F2 + F3;
    F = reshape(F', [], 1);
end

function F = neighbors_F(f_i, f_j, fv, step, nelx, nely)

    nodes_num = (nelx+1)*(nely+1);
    F = zeros(nodes_num, 2);
    
    for i = max(f_i-step, 1) : min(f_i+step, nelx+1)
        for j = max(f_j-step, 1) : min(f_j+step, nely+1)
            fx1 = fv * (1 - ((i-f_i)/step)^2);
            fy1 = fv * (1 - ((j-f_j)/step)^2);

            nid = j + (i-1) * (nely+1);
            F(nid, :) = [0, fy1];
        end
    end
end