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.
454 lines
16 KiB
454 lines
16 KiB
function ani_opt_3d_case2_fitting()
|
|
nx=60;
|
|
ny=20;
|
|
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;
|
|
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 * nele + 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.
|
|
% c_MMA = 1e3.*ones(m,1); % Column vector with the constants c_i in the terms c_i*y_i.
|
|
c_MMA = [0.5e3.*ones(m-1,1);1e3]; % Column vector with the constants c_i in the terms c_i*y_i.
|
|
d = zeros(m,1); % Column vector with the constants d_i in the terms 0.5*d_i*(y_i)^2.
|
|
|
|
|
|
% aold1 = a(:); % xval, one iteration ago (provided that iter>1).
|
|
% aold2 = a(:); % xval, two iterations ago (provided that iter>2).
|
|
% bold1 = b(:);
|
|
% bold2 = b(:);
|
|
% cold1 = c(:);
|
|
% cold2 = c(:);
|
|
%
|
|
% 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).
|
|
% alow=low;aupp=upp;
|
|
% blow=low;bupp=upp;
|
|
% clow=low;cupp=upp;
|
|
|
|
maxloop = 300;
|
|
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,ch_da,ch_db,ch_dc]=fitting_val(aPhys,bPhys,cPhys);
|
|
ke_all=lk_all(ch);
|
|
U=FEM(nx,ny,nz,ke_all);
|
|
|
|
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-aPhys(:)-bPhys(:)-cPhys(:);
|
|
cnstrt2 = -(-1+aPhys(:)+bPhys(:)-cPhys(:));
|
|
cnstrt3 = -(1-aPhys(:)+bPhys(:)-cPhys(:));
|
|
cnstrt4 = -(1+aPhys(:)-bPhys(:)-cPhys(:));
|
|
cnstrt1 =cnstrt1./6;
|
|
cnstrt2 =cnstrt2./8;
|
|
cnstrt3 =cnstrt3./6;
|
|
cnstrt4 =cnstrt4./6;
|
|
fval = [cnstrt1;cnstrt2;cnstrt3;cnstrt4;...
|
|
sum(vol)/(volfrac*nele)-1];
|
|
% EYES=eye(nele);
|
|
EYES=sparse(1:nele,1:nele,1);
|
|
dfdx = [EYES.*(-1)./6,EYES.*(-1)./6,EYES.*(-1)./6;
|
|
EYES.*(-1)./8,EYES.*(-1)./8,EYES.*( 1)./8;
|
|
EYES.*( 1)./6,EYES.*(-1)./6,EYES.*( 1)./6;
|
|
EYES.*(-1)./6,EYES.*( 1)./6,EYES.*( 1)./6;
|
|
dv(:)' / (volfrac*nele)]; % m X n
|
|
[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)));
|
|
|
|
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;
|
|
|
|
plot_convergence_all(loop,maxloop,sumCompliance,volfraction);
|
|
|
|
% [c,cPhys,cold1,cold2,change,compliance,vol,clow,cupp,comp]=...
|
|
% optimizeC(a,b,c,aPhys,bPhys,cPhys,cold1,cold2,Hs,H,volfrac,clow,cupp,nx,ny,nz,loop);
|
|
% sumCompliance(loop,1)=compliance;
|
|
% volfraction(loop,1) = mean(vol);
|
|
% disp([' c-It.: ' sprintf('%4i',loop)...
|
|
% ' compliance.: ' sprintf('%10.10f',compliance)...
|
|
% ' volfrac:' sprintf('%10.10f',mean(vol))...
|
|
% ' change.: ' sprintf('%6.3f',change )]);
|
|
% c_log{loop,1}=c;c_log{loop,2}=cPhys;
|
|
% % plot_comp(comp);
|
|
%
|
|
% [a,aPhys,aold1,aold2,change,compliance,vol,alow,aupp,comp]=...
|
|
% optimizeA(a,b,c,aPhys,bPhys,cPhys,aold1,aold2,Hs,H,volfrac,alow,aupp,nx,ny,nz,loop);
|
|
% sumCompliance(loop,2)=compliance;
|
|
% volfraction(loop,2) = mean(vol);
|
|
% disp([' a-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;
|
|
% % plot_comp(comp);
|
|
%
|
|
% [b,bPhys,bold1,bold2,change,compliance,vol,blow,bupp,comp]=...
|
|
% optimizeB(a,b,c,aPhys,bPhys,cPhys,bold1,bold2,Hs,H,volfrac,blow,bupp,nx,ny,nz,loop);
|
|
% sumCompliance(loop,3)=compliance;
|
|
% volfraction(loop,3) = mean(vol);
|
|
% disp([' b-It.: ' sprintf('%4i',loop)...
|
|
% ' compliance.: ' sprintf('%10.10f',compliance)...
|
|
% ' volfrac:' sprintf('%10.10f',mean(vol))...
|
|
% ' change.: ' sprintf('%6.3f',change )]);
|
|
% b_log{loop,1}=b;b_log{loop,2}=bPhys;
|
|
% plot_comp(comp);
|
|
|
|
% plot_distributions(aPhys,bPhys,cPhys);
|
|
% plot_convergence(loop,maxloop,sumCompliance,volfraction);
|
|
|
|
tim(loop)=toc;
|
|
save('ani_rmin1.5_Beam3d_Case2_fitting.mat','aPhys','bPhys','cPhys',...
|
|
'sumCompliance','volfraction','a_log','b_log','c_log','tim');
|
|
end
|
|
end
|
|
|
|
% function [c,cPhys,cold1,cold2,change,compliance,vol,clow,cupp,comp]=...
|
|
% optimizeC(a,b,c,aPhys,bPhys,cPhys,cold1,cold2,Hs,H,volfrac,clow,cupp,nx,ny,nz,loop)
|
|
% m = 1; % The number of general constraints.
|
|
% n = nx * ny * nz; % The number of design variables x_j.
|
|
% nele = nx * ny * nz;
|
|
% 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.
|
|
% x_MMA = 1e1.*ones(m,1); % Column vector with the constants c_i in the terms c_i*y_i.
|
|
% d = zeros(m,1); % Column vector with the constants d_i in the terms 0.5*d_i*(y_i)^2.
|
|
%
|
|
% ch=interp_ch(aPhys,bPhys,cPhys);
|
|
% ke_all=lk_all(ch);
|
|
% U=FEM(nx,ny,nz,ke_all);
|
|
%
|
|
%
|
|
% %OPTIMIZE c
|
|
% % c<=min(-1+a+b,1-a+b,1+a-b)
|
|
% % c_min = max(1-a+c,-1+a+c);
|
|
% c_min = -1-aPhys-bPhys;
|
|
% c_max = min(min(-1+aPhys+bPhys,1-aPhys+bPhys),1+aPhys-bPhys);
|
|
% xmin = max(c_min(:),-5.00);% Column vector with the lower bounds for the variables x_j.
|
|
% xmax = min(c_max(:),1.00);
|
|
%
|
|
% [ch_dc]=interp_ch_dc(aPhys,bPhys,cPhys);
|
|
% ke_dc=lk_all(ch_dc);
|
|
% [comp,compliance,dcdc] = COMP(U,ke_all,ke_dc,nx,ny,nz);
|
|
% [dvdc]=interp_dvdc(aPhys(:),bPhys(:),cPhys(:));
|
|
% dv=dvdc;
|
|
% xval=cPhys(:);
|
|
% f0val = compliance;
|
|
% df0dx = dcdc(:);
|
|
% vol=interp_v(aPhys(:),bPhys(:),cPhys(:));
|
|
% fval = sum(vol)/(volfrac*nele)-1;
|
|
% dfdx = dv(:)' / (volfrac*nele);
|
|
% [xmma, ~, ~, ~, ~, ~, ~, ~, ~, clow,cupp] = ...
|
|
% mmasub(m, n, loop, xval, xmin, xmax, cold1, cold2, ...
|
|
% f0val,df0dx,fval,dfdx,clow,cupp,a0,a1,x_MMA,d);
|
|
% c = reshape(xmma,nz,ny,nx);
|
|
% cPhys(:) = (double(H)*double(c(:)))./double(Hs);
|
|
% cold2 = cold1(:);
|
|
% cold1 = xval(:);
|
|
% change = max(max(abs(c(:)-xval)));
|
|
% vol=interp_v(aPhys(:),bPhys(:),cPhys(:));
|
|
% end
|
|
|
|
% function [b,bPhys,bold1,bold2,change,compliance,vol,blow,bupp,comp]=...
|
|
% optimizeB(a,b,c,aPhys,bPhys,cPhys,bold1,bold2,Hs,H,volfrac,blow,bupp,nx,ny,nz,loop)
|
|
% m = 1; % The number of general constraints.
|
|
% n = nx * ny * nz; % The number of design variables x_j.
|
|
% nele = nx * ny * nz;
|
|
% 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.
|
|
% x_MMA = 0.8e1.*ones(m,1); % Column vector with the constants c_i in the terms c_i*y_i.
|
|
% d = zeros(m,1); % Column vector with the constants d_i in the terms 0.5*d_i*(y_i)^2.
|
|
%
|
|
% ch=interp_ch(aPhys,bPhys,cPhys);
|
|
% ke_all=lk_all(ch);
|
|
% U=FEM(nx,ny,nz,ke_all);
|
|
%
|
|
% %OPTIMIZE b
|
|
% % b>=max(1-a+c,-1+a+c)
|
|
% % a<=1+a-c
|
|
% b_min = max(max(1-aPhys+cPhys,-1+aPhys+cPhys),-1-aPhys-cPhys);
|
|
% b_max = 1+aPhys-cPhys;
|
|
% xmin = max(b_min(:),0.52);% Column vector with the lower bounds for the variables x_j.
|
|
% xmax = min(b_max(:),1.98);
|
|
%
|
|
% [ch_db]=interp_ch_db(aPhys,bPhys,cPhys);
|
|
% ke_db=lk_all(ch_db);
|
|
% [comp,compliance,dcdb] = COMP(U,ke_all,ke_db,nx,ny,nz);
|
|
% [dvdb]=interp_dvdb(aPhys(:),bPhys(:),cPhys(:));
|
|
% dv=dvdb;
|
|
% xval=bPhys(:);
|
|
% f0val = compliance;
|
|
% df0dx = dcdb(:);
|
|
% vol=interp_v(aPhys(:),bPhys(:),cPhys(:));
|
|
% fval = sum(vol)/(volfrac*nele)-1;
|
|
% dfdx = dv(:)' / (volfrac*nele);
|
|
% [xmma, ~, ~, ~, ~, ~, ~, ~, ~, blow,bupp] = ...
|
|
% mmasub(m, n, loop, xval, xmin, xmax, bold1, bold2, ...
|
|
% f0val,df0dx,fval,dfdx,blow,bupp,a0,a1,x_MMA,d);
|
|
% b = reshape(xmma,nz,ny,nx);
|
|
% bPhys(:) = (double(H)*double(b(:)))./double(Hs);
|
|
% bold2 = bold1(:);
|
|
% bold1 = xval(:);
|
|
% change = max(max(abs(b(:)-xval)));
|
|
% vol=interp_v(aPhys(:),bPhys(:),cPhys(:));
|
|
% end
|
|
|
|
% function [a,aPhys,aold1,aold2,change,compliance,vol,alow,aupp,comp]=...
|
|
% optimizeA(a,b,c,aPhys,bPhys,cPhys,aold1,aold2,Hs,H,volfrac,alow,aupp,nx,ny,nz,loop)
|
|
% m = 1; % The number of general constraints.
|
|
% n = nx * ny * nz; % The number of design variables x_j.
|
|
% nele = nx * ny * nz;
|
|
%
|
|
% 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.
|
|
% x_MMA = 0.8e1.*ones(m,1); % Column vector with the constants c_i in the terms c_i*y_i.
|
|
% d = zeros(m,1); % Column vector with the constants d_i in the terms 0.5*d_i*(y_i)^2.
|
|
%
|
|
% ch=interp_ch(aPhys,bPhys,cPhys);
|
|
% ke_all=lk_all(ch);
|
|
% U=FEM(nx,ny,nz,ke_all);
|
|
%
|
|
% a_min = max(max(1-bPhys+cPhys,-1+bPhys+cPhys),-1-bPhys-cPhys);
|
|
% a_max = 1+bPhys-cPhys;
|
|
% xmin = max(a_min(:),0.52);% Column vector with the lower bounds for the variables x_j.
|
|
% xmax = min(a_max(:),1.98);
|
|
%
|
|
% [ch_da]=interp_ch_da(aPhys,bPhys,cPhys);
|
|
% ke_da=lk_all(ch_da);
|
|
% [comp,compliance,dcda] = COMP(U,ke_all,ke_da,nx,ny,nz);
|
|
%
|
|
% [dvda]=interp_dvda(aPhys(:),bPhys(:),cPhys(:));
|
|
% dv=dvda;
|
|
% xval=aPhys(:);
|
|
% f0val = compliance;
|
|
% df0dx = dcda(:);
|
|
% vol=interp_v(aPhys(:),bPhys(:),cPhys(:));
|
|
% fval = sum(vol)/(volfrac*nele)-1;
|
|
% dfdx = dv(:)' / (volfrac*nele);
|
|
% [xmma, ~, ~, ~, ~, ~, ~, ~, ~, alow,aupp] = ...
|
|
% mmasub(m, n, loop, xval, xmin, xmax, aold1, aold2, ...
|
|
% f0val,df0dx,fval,dfdx,alow,aupp,a0,a1,x_MMA,d);
|
|
% a = reshape(xmma,nz,ny,nx);
|
|
% aPhys(:) = (double(H)*double(a(:)))./double(Hs);
|
|
% aold2 = aold1(:);
|
|
% aold1 = xval(:);
|
|
% change = max(max(abs(a(:)-xval)));
|
|
% vol=interp_v(aPhys(:),bPhys(:),cPhys(:));
|
|
% end
|
|
|
|
|
|
%%%%%%%%%% 3D FEM with brick element %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
function [U]=FEM(nx,ny,nz,ke_all)
|
|
N = 3*(nx+1)*(ny+1)*(nz+1);
|
|
% K = sparse(N, N);
|
|
nele=nx*ny*nz;
|
|
si=zeros(1,nele*24^2);
|
|
sj=zeros(1,nele*24^2);
|
|
sv=zeros(1,nele*24^2);
|
|
|
|
for elx = 1:nx
|
|
for ely = 1:ny
|
|
for elz = 1:nz
|
|
% n1 = (nz+1)*(ny+1)*elx + (nz+1)*(ely-1) + elz;
|
|
% n2 = n1 - (nz+1)*(ny+1);
|
|
n2 = (nz+1)*(ny+1)*elx + (nz+1)*(ely-1) + elz;
|
|
n1 = n2-(nz+1)*(ny+1);
|
|
n3 = n1 + (nz+1);
|
|
n4 = n2 + (nz+1);
|
|
n5 = n1 + 1;
|
|
n6 = n2 + 1;
|
|
n7 = n3 + 1;
|
|
n8 = n4 + 1;
|
|
edof = [3*n1-2; 3*n1-1; 3*n1;
|
|
3*n2-2; 3*n2-1; 3*n2;
|
|
3*n3-2; 3*n3-1; 3*n3;
|
|
3*n4-2; 3*n4-1; 3*n4;
|
|
3*n5-2; 3*n5-1; 3*n5;
|
|
3*n6-2; 3*n6-1; 3*n6;
|
|
3*n7-2; 3*n7-1; 3*n7;
|
|
3*n8-2; 3*n8-1; 3*n8];
|
|
ne=(elx-1)*ny*nz+(ely-1)*nz+elz;
|
|
si((ne-1)*24^2+1:ne*24^2)=kron(ones(1,24),edof.');
|
|
sj((ne-1)*24^2+1:ne*24^2)=kron(edof.',ones(1,24));
|
|
ke=ke_all(:,:,elz,ely,elx);
|
|
sv((ne-1)*24^2+1:ne*24^2)=ke(:);
|
|
% K(edof,edof) = K(edof,edof) + ke_all(:,:,elz,ely,elx);
|
|
end
|
|
end
|
|
end
|
|
K = sparse(si,sj,sv,N,N);
|
|
|
|
%loadnid = nz+1:nz+1:(nx+1)*(ny+1)*(nz+1);
|
|
% loadnid = (nx/2+1)*(ny+1)*(nz+1) - nz/2;
|
|
loadnid = nx*(ny+1)*(nz+1) + (nz+1)*ny/2 + nz/2 + 1;
|
|
loaddof = 3.*loadnid(:);
|
|
F = sparse(N,1); F(loaddof)=-1;
|
|
|
|
fixednid=[1:(nz+1)*6,(ny+1-6)*(nz+1)+1:(ny+1)*(nz+1)];
|
|
%fixednid=[1,(ny+1)*(nz+1)-nz,nx*(ny+1)*(nz+1)+1,(nx+1)*(ny+1)*(nz+1)-nz];
|
|
fixeddofs = [3*fixednid(:); 3*fixednid(:)-1; 3*fixednid(:)-2];
|
|
alldofs = [1:N];
|
|
freedofs = setdiff(alldofs,fixeddofs);
|
|
U = zeros(N,1);
|
|
U(freedofs,:) = K(freedofs,freedofs) \ F(freedofs,:);
|
|
end
|
|
|
|
|
|
function plot_comp(comp)
|
|
figure(1)
|
|
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')
|
|
pause(1e-6);
|
|
end
|
|
|
|
function plot_distributions(aPhys,bPhys,cPhys)
|
|
figure(2)
|
|
set(gcf,'position',[410 650 1200 300])
|
|
subplot(1,3,1)
|
|
colormap(gray); imagesc(squeeze(aPhys(:,1,:)));
|
|
axis equal; axis tight; axis off;
|
|
title('a Distribution')
|
|
set(gca,'YDir','normal')
|
|
subplot(1,3,2)
|
|
colormap(gray); imagesc(squeeze(bPhys(:,1,:)));
|
|
axis equal; axis tight; axis off;
|
|
title('b Distribution')
|
|
set(gca,'YDir','normal')
|
|
subplot(1,3,3)
|
|
colormap(gray); imagesc(squeeze(cPhys(:,1,:)));
|
|
axis equal; axis tight; axis off;
|
|
title('c Distribution')
|
|
set(gca,'YDir','normal')
|
|
pause(1e-6);
|
|
end
|
|
|
|
% function [comp,compliance,dcdx] = COMP(U,ke_all,ke_dx,nx,ny,nz)
|
|
% comp =zeros(nz,ny,nx);
|
|
% compliance = 0.;
|
|
% dcdx=zeros(nz,ny,nx);
|
|
% for elz = 1:nz
|
|
% % z = nz*h - (elz-1)*h;
|
|
% for ely = 1:ny
|
|
% % y = (ely-1)*h + 0.5*h;
|
|
% for elx = 1:nx
|
|
% % x = (elx-1)*h + 0.5*h;
|
|
% n2 = (nz+1)*(ny+1)*elx + (nz+1)*(ely-1) + elz;
|
|
% n1 = n2-(nz+1)*(ny+1);
|
|
% n3 = n1 + (nz+1);
|
|
% n4 = n2 + (nz+1);
|
|
% n5 = n1 + 1;
|
|
% n6 = n2 + 1;
|
|
% n7 = n3 + 1;
|
|
% n8 = n4 + 1;
|
|
% Ue = U([3*n1-2; 3*n1-1; 3*n1;
|
|
% 3*n2-2; 3*n2-1; 3*n2;
|
|
% 3*n3-2; 3*n3-1; 3*n3;
|
|
% 3*n4-2; 3*n4-1; 3*n4;
|
|
% 3*n5-2; 3*n5-1; 3*n5;
|
|
% 3*n6-2; 3*n6-1; 3*n6;
|
|
% 3*n7-2; 3*n7-1; 3*n7;
|
|
% 3*n8-2; 3*n8-1; 3*n8]);
|
|
%
|
|
% comp(elz,ely,elx) = Ue'* ke_all(:,:,elz,ely,elx)*Ue;
|
|
% compliance = compliance + comp(elz,ely,elx);
|
|
% dcdx(elz,ely,elx) =-Ue'*ke_dx(:,:,elz,ely,elx)*Ue;
|
|
% end
|
|
% end
|
|
% 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))
|
|
xlim([0 maxloop])
|
|
subplot(2,1,2)
|
|
plotV=volfraction.';
|
|
plot([1:loop].',plotV(1:loop))
|
|
xlim([0 maxloop])
|
|
pause(1e-6);
|
|
end
|