function xPhys = myInvHom3D_sym8()

E = 0.1; nu = 0.15;
D_expect = E/(1+nu)/(1-2*nu)* [ 1-nu nu nu 0 0 0;
                            nu 1-nu nu 0 0 0;
                            nu nu 1-nu 0 0 0;
                            0 0 0 (1-2*nu)/2 0 0;
                            0 0 0 0 (1-2*nu)/2 0;
                            0 0 0 0 0 (1-2*nu)/2;];
disp(D_expect)
num=20; h = 1.0/num;
nelx = num; nely = num; nelz = num;
nelxyz = num;

% nele = nelx*nely*nelz;
volfrac = 0.5;
rmin = 2;    

E0 = 1; nu = 0.3;
Emin = 1e-9;
penal = 9;
I=eye(6);

nelx_sym8 = nelx/2;
nely_sym8 = nely/2;
nelz_sym8 = nelz/2;
nele_sym8 = nelx_sym8*nely_sym8*nelz_sym8;


D0 = E0/(1+nu)/(1-2*nu)* [ 1-nu nu nu 0 0 0;
                            nu 1-nu nu 0 0 0;
                            nu nu 1-nu 0 0 0;
                            0 0 0 (1-2*nu)/2 0 0;
                            0 0 0 0 (1-2*nu)/2 0;
                            0 0 0 0 0 (1-2*nu)/2;];

% load KEFB_20.mat KE Fe intB
% load KEFBE_40.mat KE Fe intB
datafile = ['homKeFe_',num2str(num),'.mat'];
tic
if exist(datafile,'file')
    load(datafile,'Ke','Fe','intB','intDB')
else
    [Ke,Fe,intB,intDB] = symbolicKeFe(num,E0,nu);%========================================可预计算                 
    save(datafile,'Ke','Fe','intB','intDB')
end
toc
[eleidx,mesh,VE]=periodicMesh(num);
nele=size(mesh,1);nnode=nele;
edofMat = zeros(size(mesh,1),24);
edofMat(:,1:3:24) = mesh.*3-2;
edofMat(:,2:3:24) = mesh.*3-1;
edofMat(:,3:3:24) = mesh.*3;
iK = reshape(kron(edofMat,ones(24,1))',24*24*nele,1);
jK = reshape(kron(edofMat,ones(1,24))',24*24*nele,1);

iF = reshape(kron(edofMat,ones(6,1))',24*6*nele,1);
jF = reshape(kron(repmat([1 2 3 4 5 6],[nele,1]),ones(1,24))',24*6*nele,1);

% [H,Hs] = filter(nelx,nely,nelz,rmin);
[H,Hs] = filter(nelx_sym8,nely_sym8,nelz_sym8,rmin);
% x = ones(nely,nelx,nelz).*volfrac;
% x(nely/2:nely/2+1,nelx/2:nelx/2+1,nelz /2:nelz/2+1) = Emin;
% x = rand(nely,nelx,nelz);
% xPhys = x;

x_sym8 = ones(nely_sym8,nelx_sym8,nelz_sym8).*volfrac;
x_sym8(nely_sym8/2:nely_sym8/2+1,nelx_sym8/2:nelx_sym8/2+1,nelz_sym8 /2:nelz_sym8/2+1) = Emin;
x = expandX_3Dsym8(x_sym8);
xPhys = x;
xPhys_sym8=x_sym8;

loopbeta = 0; loop = 0; maxloop = 50;
change = 1;


n = nele_sym8;
m = 1;
xmin  = ones(n,1).*Emin;
xmax  = ones(n,1).*1.0;
xold1 = x(:);
xold2 = x(:);
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.
a1    = zeros(m,1);       % Column vector with the constants a_i in the terms a_i*z.
c_MMA = 1e3.*ones(m,1);
d     = zeros(m,1);       % Column vector with the constants d_i in the terms 0.5*d_i*(y_i)^2.

OBJ = zeros(maxloop,1);
VOL = zeros(maxloop,1);
BWR = zeros(maxloop,1);
% figure;
dQ_sym8 = cell(6,6);
while loop < maxloop && change > 1e-4
    loop = loop+1; loopbeta = loopbeta+1;
    sK = reshape(Ke(:)*((x(:).^penal).'),24*24*nele,1);
    K = sparse(iK,jK,sK); K = (K+K')/2;
    sF = reshape(Fe(:)*((x(:).^penal).'),24*6*nele,1);
    F = sparse(iF,jF,sF);
    u = zeros(nnode*3,6);
%     for i=1:6
%         u(4:end,i)=pcg(K(4:end,4:end),F(4:end,i),1e-6,300);
%     end
    u(4:end,:)=K(4:end,4:end)\F(4:end,:);
    x_vec=x(:);%x = x(:);

    DH=zeros(6,6);
    dDHdx = zeros(6,6,nele);
    for iele=1:nele
        element=mesh(iele,:);
        edof=zeros(1,24);
        edof(1:3:24)=3*element(:)-2;
        edof(2:3:24)=3*element(:)-1;
        edof(3:3:24)=3*element(:);
        uie=u(edof,:);
        De = D0*(x_vec(iele).^penal);
        DH=DH+De*(I*h^3-intB*uie);
        dDHdx(:,:,iele) = penal*x_vec(iele)^(penal-1).*(D0.*h^3 - 2.*intDB*uie + uie.'*Ke*uie);
    end
    for i=1:6
        for j=1:6
            dQx = reshape(squeeze(dDHdx(i,j,:)),nely,nelx,nelz);
            dQx_z1 = dQx(:,:,1:nelz_sym8) + flip(dQx(:,:,nelz_sym8+1:end),3);
            dQx_z1_y1 = dQx_z1(1:nely_sym8,:,:) + flipud(dQx_z1(nely_sym8+1:end,:,:));
            dQx_z1_y1_x1 = dQx_z1_y1(:,1:nelx_sym8,:) + fliplr(dQx_z1_y1(:,nelx_sym8+1:end,:));
            
            dQ_sym8{i,j} = dQx_z1_y1_x1;
        end
    end
    dQ_sym8_ele = zeros(nelxyz/2,nelxyz/2,nelxyz/2,6,6);
    for i=1:6
        for j=1:6
            dQ_sym8_ele(:,:,:,i,j) = dQ_sym8{i,j};
        end
    end
%     disp(DH)
    c = sum((DH-D_expect).^2,'all');
%     dc = zeros(nele,1);
%     for iele = 1:nele
%         dc(iele) = sum((DH-D_expect).*dDHdx(:,:,iele).*2,'all');
%     end

    dc = zeros(nely_sym8,nelx_sym8,nelz_sym8);
    for ix = 1:nelxyz/2
        for iy = 1:nelxyz/2
            for iz = 1:nelxyz/2
                dc(iy,ix,iz) = sum((DH-D_expect).'.*squeeze(dQ_sym8_ele(iy,ix,iz,:,:)).*2,'all');
            end
        end
    end
%     dc = dc(:);
    disp('sensitivity:')    
disp([min(dc,[],'all'),        max(dc,[],'all')])
%     dc = dc.*1e1;

%     c = -(DH(4,4)+DH(5,5)+DH(6,6));
%     dc = zeros(nele,1);
%     for iele = 1:nele
%         dc(iele) = -(dDHdx(4,4,iele)+dDHdx(5,5,iele)+dDHdx(6,6,iele));
%     end

%     dv = ones(nele,1);
    dv = ones(nely_sym8,nelx_sym8,nelz_sym8).*8;
    dc(:) = H*(dc(:)./Hs);  
    dv(:) = H*(dv(:)./Hs);
    
    l1 = 0; l2 = 1e9; move = 0.2;
    
%     while (l2-l1)/(l1+l2) > 1e-3
%         lmid = 0.5*(l2+l1);
%         xnew = max(Emin,max(x-move,min(1,min(x+move,x.*sqrt(-dc./dv/lmid)))));
%         xPhys(:) = (H*xnew(:))./Hs;
% %         xPhys(:) = xnew(:);
%         if sum(xPhys(:)) > volfrac*nele
%             l1 = lmid; 
%         else
%             l2 = lmid; 
%         end
%     end
%     x(:) = xPhys(:);
    
    while (l2-l1)/(l1+l2) > 1e-3
        lmid = 0.5*(l2+l1);
%         xnew_sym8 = max(Emin,max(x_sym8-move,min(1,min(x_sym8+move,x_sym8.*sqrt(-dc./dv/lmid)))));
xnew_sym8 = max(Emin,max(x_sym8-move,min(1,min(x_sym8+move,x_sym8.*(-dc./dv/lmid)))));
        xPhys_sym8(:) = (H*xnew_sym8(:))./Hs;
        if sum(xPhys_sym8(:)) > volfrac*nele_sym8, l1 = lmid; else l2 = lmid; end
    end
    change = max(abs(xnew_sym8(:)-x_sym8(:)));
    x_sym8 = xnew_sym8;
    x = expandX_3Dsym8(x_sym8);
%     xPhys = x;
    xPhys = expandX_3Dsym8(xPhys_sym8);
% xval = x;
% f0val = c;
% df0dx = dc;
% fval = sum(xPhys(:))/(volfrac*nele)-1;
% dfdx = dv.'./(volfrac*nele);
% 
% 
% [xnew, ~, ~, ~, ~, ~, ~, ~, ~, low,upp] = ...
%     mmasub(m, n, loop, xval, xmin, xmax, xold1, xold2, ...
%     f0val,df0dx,fval,dfdx,low,upp,a0,a1,c_MMA,d);
%     change = max(abs(xnew(:)-x(:)));
%     x = xnew;
%     xPhys(:) = x(:);
%     xPhys(:) = (H*xnew(:))./Hs;
    % PRINT RESULTS
    fprintf(' It.:%5i Obj.:%11.4f Vol.:%7.3f ch.:%7.3f\n',loop,c,mean(xPhys(:)),change);
%     clf; 
% display_3D(xPhys); 
disp(['black-white ratio : ',num2str(max(xPhys,[],'all')-min(xPhys,[],'all'))])
OBJ(loop) = c;
VOL(loop) = mean(xPhys(:));
BWR(loop) = max(xPhys,[],'all')-min(xPhys,[],'all');
plot_convergence_all(loop,maxloop,OBJ,VOL,BWR);
end
end

% function plot_convergence_all(loop,maxloop,sumCompliance,volfraction)
% figure(3)
% % figure()
% % loop=75;maxloop=loop;
% set(gcf,'position',[100 100 1200 400])
% subplot(2,1,1)
% plotC=sumCompliance.';
% % plot([1:loop].',plotC(1:loop))
% semilogy([1:loop].',plotC(1:loop))
% xlim([0 maxloop])
% subplot(2,1,2)
% plotV=volfraction.';
% plot([1:loop].',plotV(1:loop))
% xlim([0 maxloop])
% pause(1e-6);
% end

function display_3D(rho)
figure(1)
    [nely,nelx,nelz] = size(rho);
    hx = 1; hy = 1; hz = 1;            % User-defined unit element size
    face = [1 2 3 4; 2 6 7 3; 4 3 7 8; 1 5 8 4; 1 2 6 5; 5 6 7 8];
%     set(gcf,'Name','ISO display','NumberTitle','off');
    for k = 1:nelz
        z = (k-1)*hz;
        for i = 1:nelx
            x = (i-1)*hx;
            for j = 1:nely
                y = nely*hy - (j-1)*hy;
                if (rho(j,i,k) > 0.5)  % User-defined display density threshold
                    vert = [x y z; x y-hx z; x+hx y-hx z; x+hx y z; x y z+hx;x y-hx z+hx; x+hx y-hx z+hx;x+hx y z+hx];
                    vert(:,[2 3]) = vert(:,[3 2]); vert(:,2,:) = -vert(:,2,:);
                    patch('Faces',face,'Vertices',vert,'FaceColor',[0.2+0.8*(1-rho(j,i,k)),0.2+0.8*(1-rho(j,i,k)),0.2+0.8*(1-rho(j,i,k))]);
                    hold on;
                end
            end
        end
    end
    axis equal; axis tight; axis off; box on; view([30,30]); pause(1e-6);
end

function [eleidx,mesh,VE]=periodicMesh(num)
% [0,1]^3 
nele=num^3;
nnod=(num+1)^3;
nodenrs=reshape(1:nnod,1+num,1+num,1+num);
edofvec=reshape(nodenrs(1:end-1,1:end-1,1:end-1),nele,1);
edofMat=repmat(edofvec,1,8)+repmat([0 1 num+1 num+2 (num+1)^2+[0 1 num+1 num+2]],nele,1);

nnp=nele;%total number of unique nodes
nnpArray=reshape(1:nnp,num,num,num);
eleidx=nnpArray;%element indices

nnpArray(end+1,:,:)=nnpArray(1,:,:);
nnpArray(:,end+1,:)=nnpArray(:,1,:);
nnpArray(:,:,end+1)=nnpArray(:,:,1);

dofvector=zeros(nnod,1);
dofvector(:)=nnpArray(:);
% dofvector(2:3:end)=3*nnpArray(:)-1;
% dofvector(3:3:end)=3*nnpArray(:);

mesh=dofvector(edofMat);

% coordinates of the 1st node of each element
h=1.0/num;
%V=zeros(num,num,num,3);
VE=zeros(nele,3);
for i=1:num
    for j=1:num
        for k=1:num
            %V(i,j,k,:)=[(i-1)*h,(j-1)*h,(k-1)*h];
            VE(eleidx(i,j,k),:)=[(i-1)*h,(j-1)*h,(k-1)*h];
        end
    end
end
end

function [Ke,Fe,intB,intDB] = symbolicKeFe(num,E,m)
h=1.0/num;
syms x y z
x0=0; y0=0; z0=0;
N(1) = (h+x0-x) * (h+y0-y) * (h+z0-z) / (h^3);
N(2) = (x-x0) * (h+y0-y) * (h+z0-z) / (h^3);
N(3) = (h+x0-x) * (y-y0) * (h+z0-z) / (h^3);
N(4) = (x-x0) * (y-y0) * (h+z0-z) / (h^3);
N(5) = (h+x0-x) * (h+y0-y) * (z-z0) / (h^3);
N(6) = (x-x0) * (h+y0-y) * (z-z0) / (h^3);
N(7) = (h+x0-x) * (y-y0) * (z-z0) / (h^3);
N(8) = (x-x0) * (y-y0) * (z-z0) / (h^3);

for j=1:8
    for jj=1:3
        Dif(3*jj-2,j*3+jj-3) = diff(N(j),x);
        Dif(3*jj-1,j*3+jj-3) = diff(N(j),y);
        Dif(3*jj,  j*3+jj-3) = diff(N(j),z);
    end
end
B9x24 = Dif;
for j=1:24
    B(1,j) = B9x24(1,j);
    B(2,j) = B9x24(5,j);
    B(3,j) = B9x24(9,j);
    B(4,j) = B9x24(2,j) + B9x24(4,j);
    B(5,j) = B9x24(6,j) + B9x24(8,j);
    B(6,j) = B9x24(3,j) + B9x24(7,j);
end
% Young's modulus & Poisson's ratio
%E = 2e11; m = 0.33;
a = E / ( (1+m) * (1-2*m) );
% 6 X 6 空间问题的弹性系数矩阵
D = [a*(1-m) a*m     a*m     0         0         0
    a*m      a*(1-m) a*m     0         0         0
    a*m      a*m     a*(1-m) 0         0         0
    0        0       0       a*(0.5-m) 0         0
    0        0       0       0         a*(0.5-m) 0
    0        0       0       0         0         a*(0.5-m)];

K=B.'*D*B;
Ke=zeros(24,24);
for i=1:24
    for j=i:24
        Ke(i,j)=int(int(int(K(i,j),'x',0,h),'y',0,h),'z',0,h);
    end
end
Ke=Ke+Ke.'-diag(diag(Ke));

F=B.'*D;
Fe=zeros(24,6);
for i=1:24
    for j=1:6
        Fe(i,j)=int(int(int(F(i,j),'x',0,h),'y',0,h),'z',0,h);
    end
end
intB=zeros(6,24);
for i=1:6
    for j=1:24
        intB(i,j)=double(int(int(int(B(i,j),'x',0,h),'y',0,h),'z',0,h));
    end
end

intDB=zeros(6,24);
DB = D*B;
for i=1:6
    for j=1:24
        intB(i,j)=double(int(int(int(DB(i,j),'x',0,h),'y',0,h),'z',0,h));
    end
end

end

% PREPARE FILTER
function [H,Hs] = filter(nelx,nely,nelz,rmin)  
    nele = nelx*nely*nelz;
    iH = ones(nele*(2*(ceil(rmin)-1)+1)^2,1);
    jH = ones(size(iH));
    sH = zeros(size(iH));
    k = 0;
    for k1 = 1:nelz
        for i1 = 1:nelx
            for j1 = 1:nely
                e1 = (k1-1)*nelx*nely + (i1-1)*nely+j1;
                for k2 = max(k1-(ceil(rmin)-1),1):min(k1+(ceil(rmin)-1),nelz)
                    for i2 = max(i1-(ceil(rmin)-1),1):min(i1+(ceil(rmin)-1),nelx)
                        for j2 = max(j1-(ceil(rmin)-1),1):min(j1+(ceil(rmin)-1),nely)
                            e2 = (k2-1)*nelx*nely + (i2-1)*nely+j2;
                            k = k+1;
                            iH(k) = e1;
                            jH(k) = e2;
                            sH(k) = max(0,rmin-sqrt((i1-i2)^2+(j1-j2)^2+(k1-k2)^2));
                        end
                    end
                end
            end
        end
    end
    H = sparse(iH,jH,sH);
    Hs = sum(H,2);
end