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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%