%--------------------------
 % @Author: Jingqiao Hu
 % @Date: 2022-01-12 20:53:52
 % @LastEditTime: 2022-01-12 23:27:18

 % triangulation mesh
%--------------------------
function xPhys = top88_nonregular(volfrac, e1, nu, F, fixeddofs, nodes, faces)
    penal=2;        
    Emin = 1e-9;
    
    [sK0, area] = intKE_triangle(nodes, faces, e1, nu); % sK: [36,m]
    area0 = sum(area);
    
    [iK, jK, edofMat, ~, ~] = forAssemble_trimesh(faces);
    
    alldofs = 2*size(nodes, 1);
    U = zeros(alldofs, 1);
    freedofs = setdiff(1:alldofs,fixeddofs);
    
    %% INITIALIZE ITERATION
    nele = size(faces, 1);
    x = repmat(volfrac, nele, 1);
    xPhys = x;
    loop = 0;
    change = 1;
    
    x1 = nodes(:, 1); % [n, 1]
    x1 = x1(faces);    % [m, 3]
    y1 = nodes(:, 2); % [n, 1]
    y1 = y1(faces);    % [m, 3]
    
    m     = 1;                % The number of general constraints. 绾︽潫鏉′欢涓暟
    n     = nele;      % The number of design variables x_j. 璁捐鍙橀噺x涓暟nelx * nely
    xmin  = zeros(n,1);       % Column vector with the lower bounds for the variables x_j. x鐨勪笅鐣岀殑鍒楀悜閲?
    xmax  = ones(n,1);        % Column vector with the upper bounds for the variables x_j. x鐨勪笂鐣岀殑鍒楀悜閲?
    xold1 = x(:);             % xval, one iteration ago (provided that iter>1). 涓?釜杩唬涔嬪墠鐨剎val
    xold2 = x(:);             % xval, two iterations ago (provided that iter>2). 涓や釜杩唬涔嬪墠鐨剎val
    low   = ones(n,1);        % Column vector with the lower asymptotes from the previous iteration (provided that iter>1). 涓婁竴涓凯浠ｇ殑涓嬫笎杩戠嚎鐨勫垪鍚戦噺
    upp   = ones(n,1);        % Column vector with the upper asymptotes from the previous iteration (provided that iter>1). 涓婁竴涓凯浠ｇ殑涓婃笎杩戠嚎鐨勫垪鍚戦噺
    a0    = 1;                % The constants a_0 in the term a_0*z.
    a     = zeros(m,1);       % Column vector with the constants a_i in the terms a_i*z.
    c_MMA = 1000*ones(m,1);   % Column vector with the constants c_i in the terms c_i*y_i.
    d     = zeros(m,1);       % Column vector with the constants d_i in the terms 0.5*d_i*(y_i)^2.
    
    %% START ITERATION
    while change > 0.01 && loop < 100
        loop = loop + 1;
        
        %% FE-ANALYSIS
        sK = sK0 .* (Emin+xPhys(:)'.^penal*(e1-Emin)); % [36,m]
        K = sparse(iK,jK,sK(:)); K = (K+K')/2;
        U(freedofs) = K(freedofs,freedofs)\F(freedofs);
        
        %% OBJECTIVE FUNCTION AND SENSITIVITY ANALYSIS
        U1(1, :, :) = U(edofMat)'; % [1,6,n]
        U2(:, 1, :) = U(edofMat)'; % [6,1,n]
        sK1 = sK0 .* (penal*(e1-Emin)*xPhys'.^(penal-1)); % [36,m]
        sK2 = reshape(sK1, 6, 6, []);
        dc = mtimesx(mtimesx(U1, sK2), U2);
        dc = squeeze(dc);
        
        c = F' * U;
        
        v = sum(area .* xPhys(:));
        dv = area;
        
        %% MMA
        xval  = x(:);
        f0val = c;
        df0dx = - dc(:);   
    
        fval  = v/(volfrac*area0) - 1;
        dfdx  = dv(:)';
    
        [xnew, ~, ~, ~, ~, ~, ~, ~, ~, low,upp] = mmasub(m, n, loop, xval, ...
            xmin, xmax, xold1, xold2, f0val,df0dx,fval,dfdx,low,upp,a0,a,c_MMA,d);

        % Update MMA Variables
        xPhys    = xnew(:);
        xold2    = xold1(:); xold1    = x(:);

        change = max(abs(xnew(:)-x(:)));
        x = xnew;
        
        fprintf(' It.:%5i Obj.:%11.4f Vol.:%7.3f ch.:%7.3f\n',loop,c, v/area0, change);
        figure(1); colormap(gray); patch(x1', y1', 1-x); caxis([0,1]); axis equal;
    end
end