function [ fea ] = optexec()
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  Topology Optimization example using FEATool Multiphysics v1.7
%
%  Copyright 2017-2018 Prof. Shintaro Yamasaki, Osaka University
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if( isOct )
  if( isempty(pkg('list','optim')) )
    error( ['Please install the Octave optimization package from:',char([10,10]), ...
            '    https://octave.sourceforge.io/optim/',char(10)] )
  end
  pkg load optim
end
if( ~exist('linprog') )
  error( 'MATLAB or Octave with Optimization toolbox required.' )
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       Construct fea structure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fea.sdim = {'x', 'y'};
fea = addphys(fea, @customeqn, {'des', 'phi'});
fea = addphys(fea, @customeqn, {'rho'});
fea = addphys(fea, @customeqn, {'u', 'v'});

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       Define constants and expressions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

curitr = 1;     % current iteration number
maxitr = 100;   % maximum iteration number

fea.expr = { ...
        'g1max',    0.75;                  % upper limit of the constraint function
        'radius',   0.05;                  % filter parameter
        'E0',       1e-6;                  % young's modulus for the air domain
        'Ec',       1;                     % young's modulus for the solid domain
        'nu',       0.3;                   % poisson's ratio
        'ttx',      0;                     % x-directional load
        'tty',      -1;                    % y-directional load
        'E',        '(Ec-E0)*rho^3+E0' ;   % modulus of elasticity
        'dEdrho',   '3*(Ec-E0)*rho^2' };   % derivative of E with respect to density

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       Grid generation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

grid1    = rectgrid(  52, 40, [0 1.3; 0 1] );
grid2    = rectgrid(  8,  40, [1.3 1.5; 0 1] );
fea.grid = gridmerge( grid1, 2, grid2, 4 );

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       Set fields for the design variables (des),
%       and filtered node-wise material distribution (phi)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fea.phys.ce.dvar = {'des', 'phi'};
fea.phys.ce.sfun = {'sflag1', 'sflag1'};
fea.phys.ce.eqn.seqn = { ...
        '0 = 0', ...
        '-radius*radius*(phix_x+phiy_y)+phi_t = des' ...
};
fea.phys.ce.eqn.coef = { ...
        'des0_ce', 'des_0', 'Initial condition for des', {0.4, 0.4};
        'phi0_ce', 'phi_0', 'Initial condition for phi', {0.3, 0.3} ...
};
fea.phys.ce.bdr.sel  = [1, 1, 1, 1, 1, 1];
fea.phys.ce.bdr.coef = { ...
                         'bcnd_ce', ...
                         'Dirichlet/Neumann boundary conditions for des', ...
                         'Dirichlet/Neumann boundary conditions for phi', ...
                         {'Dirichlet, r_des', 'Neumann, g_des';
                          'Dirichlet, r_phi', 'Neumann, g_phi'}, ...
                         {0, 0, 0, 0, 0, 0;
                          0, 0, 0, 0, 0, 0}, ...
                         [ ], ...
                         {0, 0, 0, 0, 0, 0;
                          0, 0, 0, 0, 0, 0} ...
                       };

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    Set a field for filtered element-wise material distribution (rho)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fea.phys.ce2.dvar = {'rho'};
fea.phys.ce2.sfun = {'sf_disc0'};
fea.phys.ce2.eqn.seqn = {'rho_t = phi'};
fea.phys.ce2.eqn.coef = {'rho0_ce2', 'rho_0', ...
                         'Initial condition for rho', {0.2, 0.2}};

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       Set the displacement field
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fea.phys.ce3.dvar = {'u', 'v'};
fea.phys.ce3.sfun = {'sflag1', 'sflag1'};
fea.phys.ce3.eqn.seqn = { ...
        '-E/(1-nu^2)*( (ux + nu*vy)_x + ((1-nu)/2*(uy+vx))_y ) = 0', ...
        '-E/(1-nu^2)*( (nu*ux + vy)_y + ((1-nu)/2*(uy+vx))_x ) = 0' ...
};
fea.phys.ce3.eqn.coef = { ...
        'u0_ce3', 'u_0', 'Initial condition for u', {0, 0};
        'v0_ce3', 'v_0', 'Initial condition for v', {0, 0} ...
};
fea.phys.ce3.bdr.sel = [1, 1, 1, 1, 1, 1];
fea.phys.ce3.bdr.coef = { ...
        'bcnd_ce', ...
        'Dirichlet/Neumann boundary conditions', ...
        'Dirichlet/Neumann boundary conditions', ...
        {'Dirichlet, r_u', 'Neumann, g_u'; 'Dirichlet, r_v', 'Neumann, g_v'}, ...
        {0, 0, 1, 0, 0, 0; 0, 0, 1, 0, 0, 0}, ...
        [], ...
        {0, 0, 0, 'ttx', 0, 0; 0, 0, 0, 'tty', 0, 0} ...
};

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       Check and parse fea struct, and evaluate initial conditions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fea = parsephys(fea);
fea = parseprob(fea);
fea.sol.u = evalinit(fea);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       Construct param structure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

param.maxitr    = maxitr;
param.gmax(1,1) = fea.expr{find(strcmp(fea.expr,'g1max')),end};

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       Make index for des, phi, rho, and uv
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tmp.ind.des = find(fea.sol.u == 0.4);
tmp.ind.phi = find(fea.sol.u == 0.3);
tmp.ind.rho = find(fea.sol.u == 0.2);
tmp.ind.uv  = find(fea.sol.u == 0);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       Construct fea_sens structure for the sensitivity analysis
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

fea_sens = fea;
fea_sens.phys.ce.eqn.seqn{2}  = '-radius*radius*(phix_x+phiy_y)+phi_t = des_t';
fea_sens.phys.ce2.eqn.seqn{1} = ' rho_t = phi_t';
fea_sens.phys.ce3.eqn.seqn{1} = ...
    '-dEdrho/(1-nu^2)*( rho_x*(ux+nu*vy) + rho_y*((1-nu)/2*(uy+vx)) ) = 0';
fea_sens.phys.ce3.eqn.seqn{2} = ...
    '-dEdrho/(1-nu^2)*( rho_y*(nu*ux+vy) + rho_x*((1-nu)/2*(uy+vx)) ) = 0';
fea_sens = rmfield(fea_sens, 'eqn');
fea_sens = parsephys(fea_sens);
fea_sens = parseprob(fea_sens);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Construct slp structure for the sequential linear programming (SLP)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

slp.Gmax   = param.gmax./abs(param.gmax);
slp.xl     = zeros(size(tmp.ind.des));
slp.xu     = ones(size(tmp.ind.des));
slp.mvlmt  = 0.05;
slp.tolfun = 1e-8;
slp.alpha  = 0.1;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       Optimization loops
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

it = 0;
u  = [];
while ( true )
  it = it + 1;

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %       Compute phi, rho, u, and v
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        fea.sol.u = solvestat(fea, 'fid', [], 'init', fea.sol.u, ...
                                   'solcomp', {'phi', 'rho', 'u', 'v'});
        u = [ u, fea.sol.u ];

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %       Compute the objective function f
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        tmp.f = intbdr('ttx*u+tty*v', fea, [4]);

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %       Compute the constraint function g
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        tmp.g(1, 1) = intsubd('rho', fea, [1 2]);

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %       Compute the sensitivties of f and g w.r.t. 'rho'
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        fea_sens.sol.u = fea.sol.u;
        [tmp.M, tmp.K, tmp.F] = assembleprob(fea_sens, 'f_a', 1, 'f_sparse', 1 );
        tmp.dF_drho = -(fea.sol.u(tmp.ind.uv).'*tmp.K(tmp.ind.uv, tmp.ind.rho)).';
        tmp.dG_drho = (ones(size(tmp.ind.rho)).'*tmp.K(tmp.ind.rho, tmp.ind.rho)).';

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %       Compute the sensitivties of f and g w.r.t. 'phi'
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        tmp.dF_dphi = (tmp.dF_drho.'*inv(-tmp.K(tmp.ind.rho, tmp.ind.rho))* ...
                       tmp.K(tmp.ind.rho, tmp.ind.phi)).';
        tmp.dG_dphi = (tmp.dG_drho.'*inv(-tmp.K(tmp.ind.rho, tmp.ind.rho))* ...
                       tmp.K(tmp.ind.rho, tmp.ind.phi)).';

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %       Compute the sensitivties of f and g w.r.t. 'des'
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        tmp.dF_ddes = -tmp.K(tmp.ind.phi, tmp.ind.des).'* ...
                       (tmp.K(tmp.ind.phi, tmp.ind.phi).'\tmp.dF_dphi);
        tmp.dG_ddes = -tmp.K(tmp.ind.phi, tmp.ind.des).'*...
                       (tmp.K(tmp.ind.phi, tmp.ind.phi).'\tmp.dG_dphi);

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %       Plot the current material distribution
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        clf
        plotbdr( fea, 'labels', 'off', 'colors', {'k'} )
        postplot( fea, 'surfexpr', 'rho', 'evalstyle', 'exact' )
        caxis([0 1])
        drawnow

        tmp.disp = sprintf('%8d,%14.6e,%14.6e\n', curitr, tmp.f, ...
                           tmp.g(1,1)/abs(param.gmax(1,1)));

        tmp.Fid = fopen('result_featool.log', 'a');
        fprintf(tmp.disp);
        fprintf(tmp.Fid, tmp.disp);
        fclose(tmp.Fid);

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %       Termination decision
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        if ( curitr >= maxitr )
          break
        end

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %       Update the design variables using SLP
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        slp.dFdx = tmp.dF_ddes;
        slp.dGdx(1, :) = tmp.dG_ddes.'/abs(param.gmax(1, 1));
        slp.Gcur = tmp.g([1], 1)./abs(param.gmax);
        slp.x = fea.sol.u(tmp.ind.des);
        slp.x(find(slp.x > 1)) = 1;
        slp.x(find(slp.x < 0)) = 0;

        [slp.solx,slp.exitflag] = opt_slp(slp.dFdx, slp.dGdx, slp.Gmax, slp.Gcur, ...
                                          slp.xl, slp.xu, slp.x, slp.mvlmt, slp.tolfun);

        if ( 1 ~= slp.exitflag )
                fprintf('Error in SLP!\n');
                slp.w = 2*( slp.Gcur-slp.Gmax >= 0 ).*(slp.Gcur-slp.Gmax);
                slp.dFdx = (slp.w.'*slp.dGdx).';
                slp.solx = slp.x-slp.alpha/max(abs(slp.dFdx))*slp.dFdx;
                slp.solx(find(slp.solx > 1)) = 1;
                slp.solx(find(slp.solx < 0)) = 0;
        end

        fea.sol.u(tmp.ind.des) = slp.solx;

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %       Update the iteration number
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

        curitr = curitr + 1;
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%       Perform Sequential Linear Programming (SLP)
%
%       arguments: dFdx      sensitivities w.r.t. the objective function(n x 1)
%                  dGdx      sensitivities w.r.t. the constraint functions(m x n)
%                  Gmax      allowable upper limit of the constraints(m x 1)
%                  Gcur      current values of the constraints(m x 1)
%                  xl        lower limits of the design variables(n x 1)
%                  xu        upper limits of the design variables(n x 1)
%                  x         design variables(n x 1)
%                  mvlmt     move limit
%                  tolfun    small value for convergence decision
%
%       returns:   solx      updated design variables
%                  exitflag  exit flag
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [solx, exitflag] = opt_slp(dFdx, dGdx, Gmax, Gcur, xl, xu, x, mvlmt, tolfun);

xl1 = xl;
xu1 = xu;

for ii=1:1:size(x, 1)
  if ( x(ii)-mvlmt >= xl(ii) )
    xl1(ii) = x(ii)-mvlmt;
  end

  if ( x(ii)+mvlmt <= xu(ii) )
    xu1(ii) = x(ii)+mvlmt;
  end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%       solve the following linear programming
%       minimize:               dFdx.'*solx
%       subject to:             dGdx*solx <= Gmax-Gcur+dGdx*x
%                                       xl(ii) <= solx(ii) <= xu(ii)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if( isOct )

  options = optimset('Maxiter', 1e+5, 'Display', 'off', 'TolFun', tolfun);
  if ( size(dGdx, 1) ~= 0 )
    [solx, fval] = linprog(dFdx, dGdx, Gmax-Gcur+dGdx*x, [], [], xl1, xu1);
  else
    [solx, fval] = linprog(dFdx, [], [], [], [], xl1, xu1);
  end
  exitflag = 1;

else

  options = optimset('Maxiter', 1e+5, 'Display', 'off', ...
                     'LargeScale', 'on', 'Diagnostics', 'off', 'TolFun', tolfun);
  if ( size(dGdx, 1) ~= 0 )
    [solx, fval, exitflag] = linprog(dFdx, dGdx, Gmax-Gcur+dGdx*x, [], [], xl1, xu1, x, options);
  else
    [solx, fval, exitflag] = linprog(dFdx, [], [], [], [], xl1, xu1, x, options);
  end

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
