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.
276 lines
9.1 KiB
276 lines
9.1 KiB
function ani_opt_concurrent_pNorm()
|
|
nx=60;
|
|
ny=2;
|
|
nz=20;
|
|
volfrac=0.5;
|
|
rmin=1.5;
|
|
[Hs,H]=filter_coef_3D(nx,ny,nz,rmin);
|
|
|
|
a(1:nz,1:ny,1:nx) = 1.0;
|
|
b(1:nz,1:ny,1:nx) = 1.0;
|
|
c(1:nz,1:ny,1:nx) = 0.0;
|
|
aPhys = a;
|
|
bPhys = b;
|
|
cPhys = c;
|
|
% -1+a+b-c>=0
|
|
% 1-a+b-c>=0
|
|
% 1+a-b-c>=0
|
|
|
|
%OPTIMIZE a
|
|
% a>=max(1-b+c,-1+b+c)
|
|
% a<=1+b-c
|
|
%OPTIMIZE b
|
|
% b>=max(1-a+c,-1+a+c)
|
|
% b<=1+a-c
|
|
%OPTIMIZE c
|
|
% c<=min(-1+a+b,1-a+b,1+a-b)
|
|
nele = nx * ny * nz;
|
|
|
|
m = 4 + 1; % The number of general constraints.
|
|
n = 3 * nx * ny * nz; % The number of design variables x_j.
|
|
|
|
xmin = [ones(nele,1).*0.6;...
|
|
ones(nele,1).*0.6;...
|
|
ones(nele,1).*(-4.9)];
|
|
% Column vector with the lower bounds for the variables x_j.
|
|
xmax = [ones(nele,1).*1.9;...
|
|
ones(nele,1).*1.9;...
|
|
ones(nele,1).*0.9];
|
|
xold1 = [ a(:); b(:);c(:)]; % xval, one iteration ago (provided that iter>1).
|
|
xold2 = [ a(:); b(:);c(:)]; % xval, two iterations ago (provided that iter>2).
|
|
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.
|
|
% p=4:
|
|
% c_MMA = [4e-5;1.5e-5;4e-5;2e-5;1e2].*5; % Column vector with the constants c_i in the terms c_i*y_i.
|
|
% p=6
|
|
% c_MMA = [4e-9;1.5e-9;4e-9;2e-9;1e2].*1e1;
|
|
% 1/p
|
|
% c_MMA = [4e-1;1.5e-1;2e4;2e-1;7e1].*1e1;
|
|
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.
|
|
|
|
maxloop = 100;
|
|
a_log = cell(maxloop,2);
|
|
b_log = cell(maxloop,2);
|
|
c_log = cell(maxloop,2);
|
|
tim=zeros(maxloop,1);
|
|
sumCompliance=zeros(maxloop,1);
|
|
volfraction = zeros(maxloop,1);
|
|
loop = 0;
|
|
change = 1.;
|
|
while change > 0.001 && loop < maxloop
|
|
loop = loop + 1;
|
|
tic
|
|
|
|
ch=interp_ch(aPhys,bPhys,cPhys);
|
|
ke_all=lk_all(ch);
|
|
U=FEM_Beam(nx,ny,nz,ke_all);
|
|
|
|
[ch_da]=interp_ch_da(aPhys,bPhys,cPhys);
|
|
[ch_db]=interp_ch_db(aPhys,bPhys,cPhys);
|
|
[ch_dc]=interp_ch_dc(aPhys,bPhys,cPhys);
|
|
ke_da=lk_all(ch_da);
|
|
ke_db=lk_all(ch_db);
|
|
ke_dc=lk_all(ch_dc);
|
|
[~,~,dcdb] = COMP(U,ke_all,ke_db,nx,ny,nz);
|
|
[~,~,dcda] = COMP(U,ke_all,ke_da,nx,ny,nz);
|
|
[comp,compliance,dcdc] = COMP(U,ke_all,ke_dc,nx,ny,nz);
|
|
|
|
[dvda]=interp_dvda(aPhys(:),bPhys(:),cPhys(:));
|
|
[dvdb]=interp_dvdb(aPhys(:),bPhys(:),cPhys(:));
|
|
[dvdc]=interp_dvdc(aPhys(:),bPhys(:),cPhys(:));
|
|
|
|
dv=[dvda;dvdb;dvdc];
|
|
xval=[aPhys(:);bPhys(:);cPhys(:)];
|
|
f0val = compliance;
|
|
df0dx = [dcda(:);dcdb(:);dcdc(:)];
|
|
|
|
vol=interp_v(aPhys(:),bPhys(:),cPhys(:));
|
|
% cnstrt1 = -1-a(:)-b(:)-c(:);
|
|
% cnstrt2 = -(-1+a(:)+b(:)-c(:));
|
|
% cnstrt3 = -(1-a(:)+b(:)-c(:));
|
|
% cnstrt4 = -(1+a(:)-b(:)-c(:));
|
|
% cnstrt1 = -1-aPhys(:)-bPhys(:)-cPhys(:);
|
|
% cnstrt2 = -(-1+aPhys(:)+bPhys(:)-cPhys(:));
|
|
% cnstrt3 = -(1-aPhys(:)+bPhys(:)-cPhys(:));
|
|
% cnstrt4 = -(1+aPhys(:)-bPhys(:)-cPhys(:));
|
|
p=100;
|
|
% cnstrt1 = sum((5-aPhys(:)-bPhys(:)-cPhys(:)).^p)-nele*6^p;
|
|
% cnstrt2 = sum((9-aPhys(:)-bPhys(:)+cPhys(:)).^p)-nele*8^p;
|
|
% cnstrt3 = sum((5+aPhys(:)-bPhys(:)+cPhys(:)).^p)-nele*6^p;
|
|
% cnstrt4 = sum((5-aPhys(:)+bPhys(:)+cPhys(:)).^p)-nele*6^p;
|
|
cnstrt1 = mean((5-aPhys(:)-bPhys(:)-cPhys(:)).^p).^(1/p)-6;
|
|
cnstrt2 = mean((9-aPhys(:)-bPhys(:)+cPhys(:)).^p).^(1/p)-8;
|
|
cnstrt3 = mean((5+aPhys(:)-bPhys(:)+cPhys(:)).^p).^(1/p)-6;
|
|
cnstrt4 = mean((5-aPhys(:)+bPhys(:)+cPhys(:)).^p).^(1/p)-6;
|
|
|
|
fval = [cnstrt1;cnstrt2;cnstrt3;cnstrt4;...
|
|
sum(vol)/(volfrac*nele)-1];
|
|
disp('fval: ')
|
|
% disp([cnstrt1,cnstrt2,cnstrt3,cnstrt4])
|
|
disp(fval.')
|
|
disp((fval.*c_MMA).')
|
|
% fval = sum(vol)/(volfrac*nele)-1;
|
|
% dc1da=-p.*(5-aPhys(:)-bPhys(:)-cPhys(:)).^(p-1);
|
|
%
|
|
% dc2da=-p.*(9-aPhys(:)-bPhys(:)+cPhys(:)).^(p-1);
|
|
% dc2dc= p.*(9-aPhys(:)-bPhys(:)+cPhys(:)).^(p-1);
|
|
%
|
|
% dc3da= p.*(5+aPhys(:)-bPhys(:)+cPhys(:)).^(p-1);
|
|
% dc3db=-p.*(5+aPhys(:)-bPhys(:)+cPhys(:)).^(p-1);
|
|
%
|
|
% dc4da=-p.*(5-aPhys(:)+bPhys(:)+cPhys(:)).^(p-1);
|
|
% dc4db= p.*(5-aPhys(:)+bPhys(:)+cPhys(:)).^(p-1);
|
|
|
|
dc1da=-sum((5-aPhys(:)-bPhys(:)-cPhys(:)).^p).^(1/p-1).*((5-aPhys(:)-bPhys(:)-cPhys(:)).^(p-1));
|
|
|
|
dc2da=-sum((9-aPhys(:)-bPhys(:)+cPhys(:)).^p).^(1/p-1).*(9-aPhys(:)-bPhys(:)+cPhys(:)).^(p-1);
|
|
dc2dc= -dc2da;
|
|
|
|
dc3da= sum((5+aPhys(:)-bPhys(:)+cPhys(:)).^p).^(1/p-1).*(5+aPhys(:)-bPhys(:)+cPhys(:)).^(p-1);
|
|
dc3db=-dc3da;
|
|
|
|
dc4da=-sum((5-aPhys(:)+bPhys(:)+cPhys(:)).^p).^(1/p-1).*(5-aPhys(:)+bPhys(:)+cPhys(:)).^(p-1);
|
|
dc4db= -dc4da;
|
|
|
|
% EYES=eye(nele);
|
|
dfdx = [dc1da.',dc1da.',dc1da.';
|
|
dc2da.',dc2da.',dc2dc.';
|
|
dc3da.',dc3db.',dc3da.';
|
|
dc4da.',dc4db.',dc4db.';
|
|
dv(:)' / (volfrac*nele)]; % m X n
|
|
% dfdx = dv(:)' / (volfrac*nele);
|
|
[xmma, ~, ~, ~, ~, ~, ~, ~, ~, low,upp] = ...
|
|
mmasub(m, n, loop, xval, xmin, xmax, xold1, xold2, ...
|
|
f0val,df0dx,fval,dfdx,low,upp,a0,a1,c_MMA,d);
|
|
a = reshape(xmma(1:nele),nz,ny,nx);
|
|
b = reshape(xmma((nele+1):2*nele),nz,ny,nx);
|
|
c = reshape(xmma((2*nele+1):n),nz,ny,nx);
|
|
aPhys(:) = (double(H)*double(a(:)))./double(Hs);
|
|
bPhys(:) = (double(H)*double(b(:)))./double(Hs);
|
|
cPhys(:) = (double(H)*double(c(:)))./double(Hs);
|
|
xold2 = xold1(:);
|
|
xold1 = xval(:);
|
|
change = max(max(abs([a(:);b(:);c(:)]-xval)));
|
|
% [a,b,c,aPhys,bPhys,cPhys,xold1,xold2,change,compliance,vol,low,upp,comp]=...
|
|
% optimizeABC(a,b,c,aPhys,bPhys,cPhys,xold1,xold2,Hs,H,volfrac,low,upp,nx,ny,nz,loop,a0,a1,x_MMA,d,...
|
|
% n,m,xmin,xmax);
|
|
sumCompliance(loop,1)=compliance;
|
|
volfraction(loop,1) = mean(vol);
|
|
disp([' It.: ' sprintf('%4i',loop)...
|
|
' compliance.: ' sprintf('%10.10f',compliance)...
|
|
' volfrac:' sprintf('%10.10f',mean(vol))...
|
|
' change.: ' sprintf('%6.3f',change )]);
|
|
a_log{loop,1}=a;a_log{loop,2}=aPhys;
|
|
b_log{loop,1}=b;b_log{loop,2}=bPhys;
|
|
c_log{loop,1}=c;c_log{loop,2}=cPhys;
|
|
% num_unconnect=connect_check(a,b,c).'
|
|
disp('connect_check: ')
|
|
disp([sum(cnstrt1),sum(cnstrt2),sum(cnstrt3),sum(cnstrt4)])
|
|
[num_unconnect,cnstrt1,cnstrt2,cnstrt3,cnstrt4]=connect_check(a,b,c);
|
|
disp(num_unconnect.')
|
|
% disp([sum(cnstrt1),sum(cnstrt2),sum(cnstrt3),sum(cnstrt4)])
|
|
% plot_comp(comp);
|
|
|
|
% plot_distributions(aPhys,bPhys,cPhys);
|
|
plot_convergence_all(loop,maxloop,sumCompliance,volfraction);
|
|
|
|
tim(loop)=toc;
|
|
save('ani_rmin1.5_concurrent.mat','aPhys','bPhys','cPhys',...
|
|
'sumCompliance','volfraction','a_log','b_log','c_log','tim');
|
|
% save ani_DATA_cab.mat
|
|
end
|
|
disp(['Time = ',num2str(sum(tim))])
|
|
% ch=interp_ch(aPhys,bPhys,cPhys);
|
|
% [ch_dc]=interp_ch_dc(aPhys,bPhys,cPhys);
|
|
% ke_all=lk_all(ch);ke_dc=lk_all(ch_dc);
|
|
% U=FEM_Beam(nx,ny,nz,ke_all);
|
|
% [comp,compliance,dcdc] = COMP(U,ke_all,ke_dc,nx,ny,nz);
|
|
|
|
disp(['C = ',num2str(compliance)])
|
|
figure()
|
|
set(gcf,'position',[100 650 300 300])
|
|
cc = squeeze(comp(:,1,:));
|
|
colormap(1-gray); imagesc(cc);
|
|
axis equal; axis tight; axis off;
|
|
% title('Compliance Distribution')
|
|
set(gca,'YDir','normal')
|
|
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))
|
|
xlim([0 maxloop])
|
|
subplot(2,1,2)
|
|
plotV=volfraction.';
|
|
plot([1:loop].',plotV(1:loop))
|
|
xlim([0 maxloop])
|
|
pause(1e-6);
|
|
end
|
|
|
|
function [a,b,c,aPhys,bPhys,cPhys,xold1,xold2,change,compliance,vol,low,upp,comp]=...
|
|
optimizeABC(a,b,c,aPhys,bPhys,cPhys,xold1,xold2,Hs,H,volfrac,low,upp,nx,ny,nz,loop,a0,a1,x_MMA,d,...
|
|
n,m,xmin,xmax)
|
|
nele = nx * ny * nz;
|
|
ch=interp_ch(aPhys,bPhys,cPhys);
|
|
ke_all=lk_all(ch);
|
|
U=FEM_Beam(nx,ny,nz,ke_all);
|
|
|
|
[ch_da]=interp_ch_da(aPhys,bPhys,cPhys);
|
|
[ch_db]=interp_ch_db(aPhys,bPhys,cPhys);
|
|
[ch_dc]=interp_ch_dc(aPhys,bPhys,cPhys);
|
|
ke_da=lk_all(ch_da);
|
|
ke_db=lk_all(ch_db);
|
|
ke_dc=lk_all(ch_dc);
|
|
[~,~,dcdb] = COMP(U,ke_all,ke_db,nx,ny,nz);
|
|
[~,~,dcda] = COMP(U,ke_all,ke_da,nx,ny,nz);
|
|
[comp,compliance,dcdc] = COMP(U,ke_all,ke_dc,nx,ny,nz);
|
|
|
|
[dvda]=interp_dvda(aPhys(:),bPhys(:),cPhys(:));
|
|
[dvdb]=interp_dvdb(aPhys(:),bPhys(:),cPhys(:));
|
|
[dvdc]=interp_dvdc(aPhys(:),bPhys(:),cPhys(:));
|
|
|
|
dv=[dvda;dvdb;dvdc];
|
|
xval=[aPhys(:);bPhys(:);cPhys(:)];
|
|
f0val = compliance;
|
|
df0dx = [dcda(:);dcdb(:);dcdc(:)];
|
|
vol=interp_v(aPhys(:),bPhys(:),cPhys(:));
|
|
fval = sum(vol)/(volfrac*nele)-1;
|
|
dfdx = dv(:)' / (volfrac*nele);
|
|
|
|
|
|
[xmma, ~, ~, ~, ~, ~, ~, ~, ~, low,upp] = ...
|
|
mmasub(m, n, loop, xval, xmin, xmax, xold1, xold2, ...
|
|
f0val,df0dx,fval,dfdx,low,upp,a0,a1,x_MMA,d);
|
|
|
|
a = reshape(xmma(1:nele),nz,ny,nx);
|
|
b = reshape(xmma((nele+1):2*nele),nz,ny,nx);
|
|
c = reshape(xmma((2*nele+1):n),nz,ny,nx);
|
|
aPhys(:) = (double(H)*double(a(:)))./double(Hs);
|
|
bPhys(:) = (double(H)*double(b(:)))./double(Hs);
|
|
cPhys(:) = (double(H)*double(c(:)))./double(Hs);
|
|
xold2 = xold1(:);
|
|
xold1 = xval(:);
|
|
|
|
change = max(max(abs([a(:);b(:);c(:)]-xval)));
|
|
|
|
end
|
|
|
|
function [num_unconnect,cnstrt1,cnstrt2,cnstrt3,cnstrt4]=connect_check(a,b,c)
|
|
num_unconnect=zeros(4,1);
|
|
cnstrt1 = -1-a(:)-b(:)-c(:);
|
|
cnstrt2 = -(-1+a(:)+b(:)-c(:));
|
|
cnstrt3 = -(1-a(:)+b(:)-c(:));
|
|
cnstrt4 = -(1+a(:)-b(:)-c(:));
|
|
num_unconnect(1)=sum(cnstrt1>0);
|
|
num_unconnect(2)=sum(cnstrt2>0);
|
|
num_unconnect(3)=sum(cnstrt3>0);
|
|
num_unconnect(4)=sum(cnstrt4>0);
|
|
end
|
|
|