% FEMLAB CODE FOR THE 4-TERMINAL DEVICE EXAMPLE OF SECTION 4.2
clear fem femadj
% DEFINE REYNOLDS NUMBER, DARCY NUMBER, LENGTH OF DESIGN DOMAIN, AND VOLUME FRACTION
Re = 50; 
Da = 1e-4; 
L0 = 3.0; 
beta = 0.4;
% DEFINE GEOMETRY, MESH, AND SUBDOMAIN/BOUNDARY GROUPS [SEE FIGURE 6]
fem.geom = rect2(0,L0,0,5) + rect2(-2,0,1,2) + rect2(-2,0,3,4) +rect2(L0,L0+2,1,2)+ rect2(L0,L0+2,3,4);
fem.mesh =meshinit(fem,'Hmaxsub',[3 0.125]);
figure;
meshplot(fem.mesh);
% subdomain groups 1:design domain 2:inlet/outlet leads
fem.equ.ind = {[3] [1 2 4 5]};
% boundary groups 1:walls 2:inlets 3:outlets 4:interior
fem.bnd.ind = {[2:3 5:8 10 12:14 16:18 20:22] [4 23] [1 24] [9 11 15 19]};
% DEFINE SPACE CO-ORDINATES, DEPENDENT VARIABLES, AND SHAPE FUNCTIONS
fem.sdim = {'x' 'y'}; 
fem.dim = {'u', 'v', 'p', 'gamma'}; 
fem.shape =[2 2 1 1];
% DEFINE CONSTANTS
fem.const.rho = 1; 
fem.const.eta = 1; 
fem.const.umax = Re;
fem.const.alphamin = 0; 
fem.const.alphamax = 1/Da; 
fem.const.q =0.1; 
Phi0 = 96*fem.const.eta*(L0+4)*fem.const.umax^2/9;
% DEFINE EXPRESSIONS ON SUBDOMAIN AND BOUNDARY GROUPS
fem.equ.expr = {'A' 'eta*(2*ux*ux+2*vy*vy+(uy+vx)*(uy+vx))+alpha*(u*u+v*v)' 'alpha' ...
{'alphamin+(alphamax-alphamin)*q*(1-gamma)/(q+gamma)'  '0'}};
fem.bnd.expr = {'B' '0'};
% DEFINE GOVERNING EQUATIONS AND INITIAL CONDITIONS [SEE EQUATIONS (8) AND (9)]
fem.form = 'general';
fem.equ.shape = {[1:4] [1:3]}; % only define gamma on subdomain group 1
fem.equ.ga = {{{'-p+2*eta*ux' 'eta*(uy+vx)'} {'eta*(uy+vx)' '-p+2*eta*vy'} {0 0} {0 0}}}; 
fem.equ.f = {{'rho*(u*ux+v*uy)+alpha*u' 'rho*(u*vx+v*vy)+alpha*v' 'ux+vy' 1}};
fem.equ.init = {{0 0 0 beta}};
% DEFINE BOUNDARY CONDITIONS
fem.bnd.shape = {[1:3]}; % do not define gamma on any boundaries
fem.bnd.r = {{'u' 'v' 0 0} ... % walls: no-slip
{'u*nx+4*umax*s*(1-s)' 'v' 0 0} ... % inlets: parabolic profile
{0 'v' 0 0} ... % outlets: normal flow
{0 0 0 0}}; % interior: nothing
fem.bnd.g = {{0 0 0 0}}; % zero prescribed external forces everywhere
% PERFORM LINEARIZATION, DEGREE-OF-FREEDOM ASSIGNMENT, AND ASSEMBLE INITIAL CONDITION
fem = femdiff(fem); fem.xmesh = meshextend(fem);
 fem.sol=asseminit(fem);
% DEFINE STRUCTURE FOR COMPUTING RIGHT-HAND-SIDE IN ADJOINT PROBLEM [SEE EQUATION (29)]
femadj = fem;
femadj.equ.ga = {{{'diff(A,ux)' 'diff(A,uy)'} {'diff(A,vx)' 'diff(A,vy)'} ... 
                               {'diff(A,px)' 'diff(A,py)'} {'diff(A,gammax)' 'diff(A,gammay)'}}}; 
femadj.equ.f = {{'diff(A,u)' 'diff(A,v)' 'diff(A,p)' 'diff(A,gamma)'}}; 
femadj.bnd.g ={{'diff(B,u)' 'diff(B,v)' 'diff(B,p)' 'diff(B,gamma)'}};
femadj.xmesh = meshextend(femadj);
% GET INDICES OF DESIGN VARIABLE IN THE GLOBAL SOLUTION VECTOR (fem.sol.u)
i4 = find(asseminit(fem,'Init',{'gamma' 1},'Out','U'));
% COMPUTE VOLUME BELOW DESIGN VARIABLE BASIS FUNCTIONS
L = assemble(fem,'Out',{'L'});
Vgamma = L(i4); 
Vdomain =sum(Vgamma);
% GET INDICES OF VELOCITY-PRESSURE VARIABLES
i123 = find(asseminit(fem,'Init',{'u' 1 'v' 1 'p' 1},'Out','U'));
% DEFINE VARIABLES AND PARAMETERS FOR MMA OPTIMIZATION ALGORITHM [SEE REFERENCES [11,12,13]]
a0 = 1; a = 0; c = 20; d = 0; xmin = 0; xmax = 1;
 xold = fem.sol.u(i4); xolder = xold; low = 0; upp = 1;
% DESIGN LOOP FOR THE ACTUAL TOPOLOGY OPTIMIZATION
for iter = 1:100
% SOLVE NAVIER-STOKES FLOW PROBLEM TO UPDATE VELOCITY AND PRESSURE
fem.sol = femnlin(fem,'Solcomp',{'u' 'v' 'p'},'U',fem.sol.u);
% SOLVE ADJOINT PROBLEM FOR LAGRANGE MULTIPLIERS
[K N] = assemble(fem,'Out',{'K' 'N'},'U',fem.sol.u); 
[L M] = assemble(femadj,'Out',{'L' 'M'},'U',fem.sol.u);
femadj.sol = femlin('In',{'K' K(i123,i123)' 'L' L(i123) 'M' zeros(size(M)) 'N' N(:,i123)});
% SENSITIVITY ANALYSIS
gamma=fem.sol.u(i4); 
Phi = postint(fem,'A','Edim',2) + postint(fem,'B','Edim',1); 
dPhidgamma = L(i4) - K(i123,i4)'*femadj.sol.u;
% PERFORM MMA STEP TO UPDATE DESIGN FIELD
x = gamma; f = Phi/Phi0; 
g = gamma'*Vgamma/Vdomain - beta; 
dfdx = dPhidgamma/Phi0; dgdx = Vgamma'/Vdomain; 
d2fdx2 = zeros(size(gamma)); d2gdx2 = zeros(size(gamma'));
[xnew,y,z,lambda,ksi,eta,mu,zeta,s,low,upp]=mmasub(1,length(gamma),iter,x,xmin,xmax,xold,xolder,f,dfdx,d2fdx2,g,dgdx,d2gdx2,low,upp,a0,a,c,d);
xolder = xold; xold = x; gamma = xnew;
% TEST CONVERGENCE
if iter >= 100 | max(abs(gamma-xold)) < 0.01 
  break 
end
% UPDATE DESIGN VARIABLE
u0 = fem.sol.u; u0(i4) = gamma; fem.sol = femsol(u0);
% DISPLAY RESULTS FOR EACH ITERATION STEP
disp(sprintf('Iter.:%3d Obj.: %8.4f Vol.: %6.3f Change: %6.3f', ...
iter,f,xold'*Vgamma,max(abs(xnew-xold))))
postplot(fem,'arrowdata',{'u' 'v'},'tridata','gamma','trimap','gray') 
axis equal; 
shg; 
pause(0.1)
end