You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
413 lines
12 KiB
413 lines
12 KiB
function xPhys = myInvHom3D_sym8_MMA()
|
|
|
|
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 = 5;
|
|
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 = 150;
|
|
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-6
|
|
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
|
|
|
|
xval = xPhys_sym8(:);
|
|
f0val = c;
|
|
df0dx = dc(:);
|
|
fval = sum(xPhys(:))/(volfrac*nele)-1;
|
|
dfdx = dv(:).'./(volfrac*nele);
|
|
|
|
[xnew_sym8, ~, ~, ~, ~, ~, ~, ~, ~, low,upp] = ...
|
|
mmasub(m, n, loop, xval, xmin, xmax, xold1, xold2, ...
|
|
f0val,df0dx,fval,dfdx,low,upp,a0,a1,c_MMA,d);
|
|
xPhys_sym8(:) = (H*xnew_sym8(:))./Hs;
|
|
xold2 = xold1(:);
|
|
xold1 = x_sym8(:);
|
|
|
|
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
|