function [DH]=homogenization3D_me(x)
Emin = 1e-9;
x = max(x,Emin);
% x: num X num X num, 
% xmin <= x <= 1 density of the material
% E0 = 1; 
% nu = 0.33;
% penal=1;

% 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;];
num=size(x,1);
h=1/num;

% load KEFB_20.mat KE Fe intB
% load KEFBE_40.mat KE Fe intB
% KE = symbolicKe(num,E0,nu);%========================================可预计算
[D0,Ke,Fe] = LocalSMRT(40);
KE=Ke;
intB = GaussIntB(h);

[eleidx,mesh,VE]=periodicMesh(num);
nele=size(mesh,1);
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);
% sK = reshape(KE(:)*(Emin+x(:)'.^penal*(1-Emin)),24*24*nele,1);
sK = reshape(KE(:)*(x(:).'),24*24*nele,1);
K = sparse(iK,jK,sK); K = (K+K')/2;

% Fe=symbolicFe_new(num,E0,nu);%======================================可预计算
iK = reshape(kron(edofMat,ones(6,1))',24*6*nele,1);
jK = reshape(kron(repmat([1 2 3 4 5 6],[nele,1]),ones(1,24))',24*6*nele,1);
sK = reshape(Fe(:)*(x(:).'),24*6*nele,1);
F = sparse(iK,jK,sK);

nnode=nele;% because of the periodic B.C.

% tic
u = zeros(nnode*3,6);
% u(4:end,:)=K(4:end,4:end)\F(4:end,:);
% toc
for i=1:6
    u(4:end,i)=pcg(K(4:end,4:end),F(4:end,i),1e-4,300);
end

% tic
% u = zeros(nnode*3,6);
% for i = 1:6
%     u(4:end,i)=pcg(K(4:end,4:end),F(4:end,i));
% end
% toc
% uRE = (u-u0)./u0;
% mean(mean(abs(uRE(4:end,:))))

% u = K\F;
% intB=compIntB(h);%==================================================可预计算

x_vec=x(:);

I=eye(6);
DH=zeros(6,6);

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);
    DH=DH+De*(I*h^3-intB*uie);
%     未定义函数或变量 'I'�?
end
% RE = (DH-DH0)./DH0;
% mean(mean(abs(RE)))
% Imask = [[1,1,1,0,0,0];[1,1,1,0,0,0];[1,1,1,0,0,0];[0,0,0,1,0,0];[0,0,0,0,1,0];[0,0,0,0,0,1]];
% mean(mean(abs(RE.*Imask)))
end