function ani_opt_cab()
nx=60;
ny=1;
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)

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(:);
n     = nx * ny * nz;
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 = 173;
a_log = cell(maxloop,2);
b_log = cell(maxloop,2);
c_log = cell(maxloop,2);
tim=zeros(maxloop,1);
sumCompliance=zeros(maxloop,3);
volfraction = zeros(maxloop,3);
loop = 0;
change = 1.;
while change > 0.001 && loop < maxloop
    loop = loop + 1;
    tic
    
    [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_cab173.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 [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 = 2e3.*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_Beam(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-a-b;
%     c_max = min(min(-1+a+b,1-a+b),1+a-b);
%     c_min = -1-aPhys-bPhys;
%     c_max = min(min(-1+aPhys+bPhys,1-aPhys+bPhys),1+aPhys-bPhys);
%     xmin  = max(c_min(:),-4.98);% Column vector with the lower bounds for the variables x_j.
%     xmax  = min(c_max(:),0.98);
    
    xmin = -2.6.*ones(size(aPhys(:)));
    xmax =  0.6.*ones(size(aPhys(:)));
    
    [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 = 1e2.*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_Beam(nx,ny,nz,ke_all);
    
    %OPTIMIZE b
    % b>=max(1-a+c,-1+a+c)
    % a<=1+a-c
%     b_min = max(max(1-a+c,-1+a+c),-1-a-c);
%     b_max = 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);
      
    xmin = 0.8.*ones(size(aPhys(:)));
    xmax = 1.2.*ones(size(aPhys(:)));
    
    [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 = 1e2.*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_Beam(nx,ny,nz,ke_all);
    
%     a_min = max(max(1-b+c,-1+b+c),-1-b-c);
%     a_max = 1+b-c;
%     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);

    xmin = 0.8.*ones(size(aPhys(:)));
    xmax = 1.2.*ones(size(aPhys(:)));
    
    [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

function plot_convergence(loop,maxloop,sumCompliance,volfraction)
figure(3)
set(gcf,'position',[100 100 1200 400])
subplot(2,1,1)
plotC=sumCompliance.';
plot([1:3*loop].',plotC(1:3*loop))
xlim([0 3*maxloop])
subplot(2,1,2)
plotV=volfraction.';
plot([1:3*loop].',plotV(1:3*loop))
xlim([0 3*maxloop])
pause(1e-6);
end