% Q: the elastic tensor % dQ: the derivatives of elastic tensor for microstructures function [Q,dQ,xPhys,xold] = inverseHomog3d_sym8_opx8() % 优化八分之一结构 E = 0.1; nu = 0.15; % E = 0.5; nu = 0.35; % E = 0.25; nu = 0.35; 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) nelxyz = 40; % 整个结构的尺寸 nelx = nelxyz; nely = nelxyz; nelz = nelxyz; nele = nelx*nely*nelz; volfrac = 0.5; rmin = 2; E0 = 1; Emin = 1e-9; nu = 0.3; penal = 5; nelx_sym8 = nelx/2; nely_sym8 = nely/2; nelz_sym8 = nelz/2; nele_sym8 = nelx_sym8*nely_sym8*nelz_sym8; % the initial definitions of the PUC 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;]; disp(D0) % dx = lx/nelx; dy = ly/nely; dz = lz/nelz; % KE = elementMatVec3D(1.0/2, 1.0/2, 1.0/2, D0); KE = elementMatVec3D(1.0, 1.0, 1.0, D0); % [KE] = lk_H8(E0,nu); Num_node = (1+nely)*(1+nelx)*(1+nelz); nodenrs = reshape(1:Num_node,1+nely,1+nelx,1+nelz); edofVec = reshape(3*nodenrs(1:end-1,1:end-1,1:end-1)+1,nelx*nely*nelz,1); edofMat = repmat(edofVec,1,24)+repmat([0 1 2 3*nely+[3 4 5 0 1 2] -3 -2 -1 3*(nelx+1)*(nely+1)+[0 1 2 3*nely+[3 4 5 0 1 2] -3 -2 -1]], nele, 1); iK = reshape(kron(edofMat,ones(24,1))',24*24*nele,1); jK = reshape(kron(edofMat,ones(1,24))',24*24*nele,1); [H,Hs] = filter(nelx_sym8,nely_sym8,nelz_sym8,rmin); % periodic boundary conditions n1 = [nodenrs(end, [1 end], 1) nodenrs(1, [end 1], 1) nodenrs(end, [1 end], end) nodenrs(1, [end 1], end)]; d1 = reshape([3*n1-2; 3*n1-1; 3*n1],3*numel(n1),1); n3 = [reshape(squeeze(nodenrs(end,1,2:end-1)),1,numel(squeeze(nodenrs(end,1,2:end-1))))... reshape(squeeze(nodenrs(1, 1, 2:end-1)),1,numel(squeeze(nodenrs(1, 1, 2:end-1))))... reshape(squeeze(nodenrs(end,2:end-1,1)),1,numel(squeeze(nodenrs(end,2:end-1,1))))... reshape(squeeze(nodenrs(1, 2:end-1, 1)),1,numel(squeeze(nodenrs(1, 2:end-1, 1))))... reshape(squeeze(nodenrs(2:end-1, 1, 1)),1,numel(squeeze(nodenrs(2:end-1, 1, 1))))... reshape(squeeze(nodenrs(2:end-1,1,end)),1,numel(squeeze(nodenrs(2:end-1,1,end))))... reshape(squeeze(nodenrs(2:end-1, 2:end-1, 1)),1,numel(squeeze(nodenrs(2:end-1, 2:end-1, 1))))... reshape(squeeze(nodenrs(2:end-1, 1, 2:end-1)),1,numel(squeeze(nodenrs(2:end-1, 1, 2:end-1))))... reshape(squeeze(nodenrs(end,2:end-1,2:end-1)),1,numel(squeeze(nodenrs(end,2:end-1,2:end-1))))]; d3 = reshape([3*n3-2; 3*n3-1; 3*n3],3*numel(n3),1); n4 = [reshape(squeeze(nodenrs(1, end, 2:end-1)),1,numel(squeeze(nodenrs(1, end, 2:end-1))))... reshape(squeeze(nodenrs(end,end,2:end-1)),1,numel(squeeze(nodenrs(end,end,2:end-1))))... reshape(squeeze(nodenrs(1, 2:end-1, end)),1,numel(squeeze(nodenrs(1, 2:end-1, end))))... reshape(squeeze(nodenrs(end,2:end-1,end)),1,numel(squeeze(nodenrs(end,2:end-1,end))))... reshape(squeeze(nodenrs(2:end-1,end,end)),1,numel(squeeze(nodenrs(2:end-1,end,end))))... reshape(squeeze(nodenrs(2:end-1, end, 1)),1,numel(squeeze(nodenrs(2:end-1, end, 1))))... reshape(squeeze(nodenrs(2:end-1,2:end-1,end)),1,numel(squeeze(nodenrs(2:end-1,2:end-1,end))))... reshape(squeeze(nodenrs(2:end-1,end,2:end-1)),1,numel(squeeze(nodenrs(2:end-1,end,2:end-1))))... reshape(squeeze(nodenrs(1, 2:end-1, 2:end-1)),1,numel(squeeze(nodenrs(1, 2:end-1, 2:end-1))))]; d4 = reshape([3*n4-2; 3*n4-1; 3*n4],3*numel(n4),1); n2 = setdiff(nodenrs(:),[n1(:);n3(:);n4(:)]); d2 = reshape([3*n2-2; 3*n2- 1; 3*n2],3*numel(n2),1); % the imposing of six linearly independent unit test strains e = eye(6); ufixed = zeros(24,6); vert_cor = [0, nelx, nelx, 0, 0, nelx, nelx, 0; 0, 0, nely, nely, 0, 0, nely, nely; 0, 0, 0, 0, nelz, nelz, nelz, nelz]; for i = 1:6 epsilon = [ e(i,1), e(i,4)/2, e(i,6)/2; e(i,4)/2, e(i,2), e(i,5)/2; e(i,6)/2, e(i,5)/2, e(i,3)]; ufixed(:,i) = reshape(epsilon*vert_cor,24,1); end % 3D boundary constraint equations wfixed = [repmat(ufixed(7:9,:),numel(squeeze(nodenrs(end,1,2:end-1))),1); repmat(ufixed( 4:6,:)-ufixed(10:12,:),numel(squeeze(nodenrs(1, 1, 2:end-1))),1); repmat(ufixed(22:24,:),numel(squeeze(nodenrs(end,2:end-1,1))),1); repmat(ufixed(13:15,:)-ufixed(10:12,:),numel(squeeze(nodenrs(1, 2:end-1, 1))),1); repmat(ufixed(16:18,:),numel(squeeze(nodenrs(2:end-1, 1, 1))),1); repmat(ufixed( 4:6,:)-ufixed(13:15,:),numel(squeeze(nodenrs(2:end-1,1,end))),1); repmat(ufixed(13:15,:),numel(squeeze(nodenrs(2:end-1, 2:end-1,1))),1); repmat(ufixed(4:6,:),numel(squeeze(nodenrs(2:end-1, 1, 2:end- 1))),1); repmat(ufixed(10:12,:),numel(squeeze(nodenrs(end,2:end-1,2:end-1))),1)]; % intial qe = cell(6,6); Q = zeros(6,6); dQ = cell(6,6); dQ_sym8 = cell(6,6); % 1/8 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; % beta = 1; % xTilde = x; % xPhys = 1-exp(-beta*xTilde)+xTilde*exp(-beta); loopbeta = 0; loop = 0; maxloop = 14; change = 1; OBJ = zeros(maxloop,1); VOL = zeros(maxloop,1); BWR = zeros(maxloop,1); % figure; while loop < maxloop && change > 0.01 loop = loop+1; loopbeta = loopbeta+1; % FE-ANALYSIS sK = reshape(KE(:)*(Emin+xPhys(:)'.^penal*(1-Emin)),24*24*nele,1); K = sparse(iK,jK,sK); K = (K+K')/2; Kr = [K(d2,d2), K(d2,d3)+K(d2,d4); K(d3,d2)+K(d4,d2),K(d3,d3)+K(d4,d3)+K(d3,d4)+K(d4,d4)]; U(d1,:) = ufixed; U([d2;d3],:) = Kr\(-[K(d2,d1);K(d3,d1)+K(d4,d1)]*ufixed-[K(d2,d4);K(d3,d4)+K(d4,d4)]*wfixed); % for i_load=1:6 % U([d2;d3],i_load) = pcg(Kr, (-[K(d2,d1);K(d3,d1)+K(d4,d1)]*ufixed(:,i_load)-[K(d2,d4);K(d3,d4)+K(d4,d4)]*wfixed(:,i_load)),1e-3,200); % end U(d4,:) = U(d3,:)+wfixed; % homog cellVolume = nele; for i = 1:6 for j = 1:6 U1 = U(:,i); U2 = U(:,j); qe{i,j} = reshape(sum((U1(edofMat)*KE).*U2(edofMat),2),nely,nelx,nelz); Q(i,j) = 1/cellVolume*sum(sum(sum((Emin+xPhys.^penal*(1-Emin)).*qe{i,j}))); dQ{i,j} = 1/cellVolume*(penal*(1-Emin)*xPhys.^(penal-1).*qe{i,j}); % compute dQ_sym8 dQx = dQ{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 disp(Q) % size(dQ_sym8{1,1}) 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 c = sum((Q-D_expect).^2,'all'); dc = zeros(nelxyz/2,nelxyz/2,nelxyz/2); for ix = 1:nelxyz/2 for iy = 1:nelxyz/2 for iz = 1:nelxyz/2 dc(iy,ix,iz) = sum((Q-D_expect).'.*squeeze(dQ_sym8_ele(iy,ix,iz,:,:)).*2,'all'); end end end % % max shear modulus % c = -(Q(4,4)+Q(5,5)+Q(6,6)); % dc = -(dQ_sym8{4,4}+dQ_sym8{5,5}+dQ_sym8{6,6}); disp('sensitivity:') disp([min(dc,[],'all'), max(dc,[],'all')]) % % c = -(Q(1,1)+Q(1,2)+Q(1,3)+Q(2,1)+Q(2,2)+Q(2,3)+Q(3,1)+Q(3,2)+Q(3,3)); % dc = -(dQ_sym8{1,1}+dQ_sym8{1,2}+dQ_sym8{1,3}+dQ_sym8{2,1}+dQ_sym8{2,2}+dQ_sym8{2,3}... % +dQ_sym8{3,1}+dQ_sym8{3,2}+dQ_sym8{3,3}); % % c = -(Q(1,1)+Q(2,2)+Q(3,3)); % dc = -(dQ_sym8{1,1}+dQ_sym8{2,2}+dQ_sym8{3,3}); % 优化出实心的 % c = -(Q(1,2)+Q(1,3)+Q(2,1)+Q(2,3)+Q(3,1)+Q(3,2)); % dc = -(dQ_sym8{1,2}+dQ_sym8{1,3}+dQ_sym8{2,1}+dQ_sym8{2,3}+dQ_sym8{3,1}+dQ_sym8{3,2}); % c = -(Q(6,6)); % 平板，不连续 % dc = -(dQ_sym8{6,6}); % c = -(Q(4,4)+Q(5,5)); % dc = -(dQ_sym8{4,4}+dQ_sym8{5,5}); dv = ones(nely_sym8,nelx_sym8,nelz_sym8).*8; % FILTERING AND MODIFICATION OF SENSITIVITIES dc(:) = H*(dc(:)./Hs); dv(:) = H*(dv(:)./Hs); % OPTIMALITY CRITERIA UPDATE l1 = 0; l2 = 1e9; move = 0.2;move=0.1; 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); xold = xPhys; % xPhys = x; xPhys = expandX_3Dsym8(xPhys_sym8); % PRINT RESULTS fprintf(' It.:%5i Obj.:%11.4f Vol.:%7.3f ch.:%7.3f\n',loop,c,mean(xPhys_sym8(:)),change); % PLOT DENSITIES if mod(loop,30)==1 && loop>1 % figure(1) % clf; display_3D(xPhys_sym8);pause(1e-6); figure(2) clf; display_3D(xPhys);pause(1e-6); end OBJ(loop) = c; VOL(loop) = mean(xPhys(:)); BWR(loop) = max(xPhys,[],'all')-min(xPhys,[],'all'); plot_convergence_all(loop,maxloop,OBJ,VOL,BWR); end figure() clf; display_3D(xPhys); %% % plot: 让微观结构看的更清楚一点，从里到外密度递增 center_x = nelx/2+0.5; center_y = nely/2+0.5; center_z = nelz/2+0.5; len_diag = sqrt(nelx^2+nely^2+nelz^2)/2; x_plot = zeros(nely,nelx,nelz); for i=1:nelx for j=1:nely for k=1:nelz dist_x = sqrt((i-center_x)^2+(j-center_y)^2+(k-center_z)^2); x_plot(j,i,k) = xPhys(j,i,k)*dist_x/len_diag; end end end figure(6) clf; display_3D(0.9*x_plot);pause(1e-6); 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 % === DISPLAY 3D TOPOLOGY (ISO-VIEW) === function display_3D(rho) [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 % 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 function Ke = elementMatVec3D(a, b, c, DH) GN_x=[-1/sqrt(3),1/sqrt(3)]; GN_y=GN_x; GN_z=GN_x; GaussWeigh=[1,1]; Ke = zeros(24,24); L = zeros(6,9); L(1,1) = 1; L(2,5) = 1; L(3,9) = 1; L(4,2) = 1; L(4,4) = 1; L(5,6) = 1; L(5,8) = 1; L(6,3) = 1; L(6,7) = 1; for i=1:length(GN_x) for j=1:length(GN_y) for k=1:length(GN_z) x = GN_x(i);y = GN_y(j);z = GN_z(k); dNx = 1/8*[-(1-y)*(1-z), (1-y)*(1-z), (1+y)*(1-z), -(1+y)*(1-z), -(1-y)*(1+z), (1-y)*(1+z), (1+y)*(1+z) -(1+y)*(1+z)]; dNy = 1/8*[-(1-x)*(1-z) -(1+x)*(1-z), (1+x)*(1-z), (1-x)*(1-z), -(1-x)*(1+z), -(1+x)*(1+z), (1+x)*(1+z), (1-x)*(1+z)]; dNz = 1/8*[-(1-x)*(1-y) -(1+x)*(1-y) -(1+x)*(1+y) -(1-x)*(1+y), (1-x)*(1-y), (1+x)*(1-y), (1+x)*(1+y), (1-x)*(1+y)]; J = [dNx;dNy;dNz]*[ -a a a -a -a a a -a ; -b -b b b -b -b b b; -c -c -c -c c c c c]'; G = [inv(J) zeros(3) zeros(3);zeros(3) inv(J) zeros(3);zeros(3) zeros(3) inv(J)]; dN(1,1:3:24) = dNx; dN(2,1:3:24) = dNy; dN(3,1:3:24) = dNz; dN(4,2:3:24) = dNx; dN(5,2:3:24) = dNy; dN(6,2:3:24) = dNz; dN(7,3:3:24) = dNx; dN(8,3:3:24) = dNy; dN(9,3:3:24) = dNz; Be = L*G*dN; Ke = Ke + GaussWeigh(i)*GaussWeigh(j)*GaussWeigh(k)*det(J)*(Be'*DH*Be); end end end end % === GENERATE ELEMENT STIFFNESS MATRIX === function [KE] = lk_H8(E0,nu) A = [32 6 -8 6 -6 4 3 -6 -10 3 -3 -3 -4 -8; -48 0 0 -24 24 0 0 0 12 -12 0 12 12 12]; k = 1/144*A'*[1; nu]; K1 = [k(1) k(2) k(2) k(3) k(5) k(5); k(2) k(1) k(2) k(4) k(6) k(7); k(2) k(2) k(1) k(4) k(7) k(6); k(3) k(4) k(4) k(1) k(8) k(8); k(5) k(6) k(7) k(8) k(1) k(2); k(5) k(7) k(6) k(8) k(2) k(1)]; K2 = [k(9) k(8) k(12) k(6) k(4) k(7); k(8) k(9) k(12) k(5) k(3) k(5); k(10) k(10) k(13) k(7) k(4) k(6); k(6) k(5) k(11) k(9) k(2) k(10); k(4) k(3) k(5) k(2) k(9) k(12) k(11) k(4) k(6) k(12) k(10) k(13)]; K3 = [k(6) k(7) k(4) k(9) k(12) k(8); k(7) k(6) k(4) k(10) k(13) k(10); k(5) k(5) k(3) k(8) k(12) k(9); k(9) k(10) k(2) k(6) k(11) k(5); k(12) k(13) k(10) k(11) k(6) k(4); k(2) k(12) k(9) k(4) k(5) k(3)]; K4 = [k(14) k(11) k(11) k(13) k(10) k(10); k(11) k(14) k(11) k(12) k(9) k(8); k(11) k(11) k(14) k(12) k(8) k(9); k(13) k(12) k(12) k(14) k(7) k(7); k(10) k(9) k(8) k(7) k(14) k(11); k(10) k(8) k(9) k(7) k(11) k(14)]; K5 = [k(1) k(2) k(8) k(3) k(5) k(4); k(2) k(1) k(8) k(4) k(6) k(11); k(8) k(8) k(1) k(5) k(11) k(6); k(3) k(4) k(5) k(1) k(8) k(2); k(5) k(6) k(11) k(8) k(1) k(8); k(4) k(11) k(6) k(2) k(8) k(1)]; K6 = [k(14) k(11) k(7) k(13) k(10) k(12); k(11) k(14) k(7) k(12) k(9) k(2); k(7) k(7) k(14) k(10) k(2) k(9); k(13) k(12) k(10) k(14) k(7) k(11); k(10) k(9) k(2) k(7) k(14) k(7); k(12) k(2) k(9) k(11) k(7) k(14)]; KE = E0/((nu+1)*(1-2*nu))*... [ K1 K2 K3 K4; K2' K5 K6 K3'; K3' K6 K5' K2'; K4 K3 K2 K1']; end