## Appendix A MATLAB Implementation of the Optimization Algorithm  
```
1 nelx=90;nely=25;volfrac=0.45;
2%tLEVEL-SETFUNCTIONINITIALIZATION
3r=nely=0.1;rRadulusofINITIALHOLES
4 hx=nelx=[repmat([1/6,5/6],1,3),repmat([0,1/3,2/3,1],1,2),1/2];
5 hY=nely*[kron([0,1/2,1],ones(1,2)),kron([1/4,3/4],ones(1,4)),1/2];
6[X,y]=meshgrid(0.1:nelx,0.1:nely);
7 dx=bsxfun(@minus,repmat(X,[1,1,numel(hX)]),reshape(hX,1,1,numel(hX)));
8 dY=bsxfun(@minus,repmat([1,1,numel(hY)]),reshape(hY,1,numel(hY)));
9 Phi=max(-3,min(3,min(sqrt(dx,'2+dY,'2)-r,[,3]));
10%%RADIALBASIFICATIONINITIALIZATION
11 nNode=(nely+1)*(nelx+1);
12 PGpx=([Phi(i,'2,iend)zeros(nely+1,1)]-[zeros(nely+1,1)Phi(i,'1:end-1)]);
13 PGpx=([Phi(i:2:end),'zeros(1,nelx+1)]-[zeros(1,nelx+1)]-[;nbi(1:end-1,'i)]);
14 Boundary=setdiff(unique([1:1:nely+1,nNode-nely):1:nNode,...
15 nely+1:nely+1:nNode,1:nely+1:nNode-nely)],
16[(nely+1)*(ceil(nelx/2+1)+0,1:nelx*(nely+1)+1,0.9*nelx*(nely+1)+1]);
17 Phi(Boundary)==-1e-9;
18 Alpha=Phi(:);
19#FINITEELEMENT ANALYSISPREPARATION
20E0=1;Emin=1e-9;nu=0.3;MATERIAL PROPERTIES
21A1=[123-6-3;31230;-6312-3;-30-312];
22A12=[-6-303-3-6-3-6-3-6;0-3-63;3-63-6-3-6];
23B1=[-4-3-293-3-49-4;-2-9-4-3;94-3-4];
24B12=[2-3-4-9-3-3-29-2-4;4B2-3-2-3-2;1];
25K=-1(1-nu*2)/24/([A112A12A11A11A111B12B12B11)];
26 eleN1=repmat(([1:nely)*,1,nelx)+kron(0:nelx-1,(nely+1)*ones(nely,1));
27 elNode=repmat([e1N1,1,4)+repmat([0,nely+1,2],1,nelx*nely,1]);
28 edoMat=kron(deltaNode,[2,2])+repmat([-1,0],nelx*nely,4);
29 ix=reshape(kron(oddMat,ones(8,1))',64*nelx*nely,1);
30%t=reshape(kron(oddMat,ones(8,8))',64*nelx*nely,1);
31%BOUNDARYCONITIONDEFINITION- sparse(2*(((nely+1)*(ceil(nelx/2)+1)),1,-100,2+nNode,1);%NODAL LOADS
* fixeddofs = [2*0.1+nelx*(nely+1)+2,2*0.9*nelx*(nely+1)+2];

* SS% ITERATION OPTIMIIZATION
* 20; gamma = 0.05; comp = zeros(nLoop,1); vol = zeros(nLoop,1);

* % FINITLE ELEMENT ANALYSIS
* [s,t] = meshgrid(-1:0.1:1,-1:0.1:1); tmpPhl = (1-s(:)).*(1-t(:))/4*Phi(eleNode(:,1))'+(1+s(:)).*(1-t(:))/4*... * Phi(elecNode(:,2))'+(1+s(:)).*(1+t(:))/4*Phi(eleNode(:,3))'+... * [1-s(:)). *(1+t(:))/4*Phi(elecNode(:,4))';
* elevol = sum(tmphpla1)'/numel(s); vol(1)'+numel(s); vol(1)'-sum(eleVol)/(lenx*nely);
* sk = reshape(KE(:)*(Emin+eleVol'*(E0-Emin)),64*nelx*nelx*nely,1);
* sparse(ik,JK,SK); K
* U(freedofs,1) = K(freedofs,freedofs)F(freedofs,1);
* elecComp = sum((cup(dotMat)*KE).*(edofMat),2). *(Emin+eleVol*(E0-Emin));

* SS% DISPLAY RESULTS
* fprintf('No.%i, Obj:#f, Vol:#fUn', {iT,comp(iT),vol(iT)});
* figure(1); contourf(Phi,[0,0]);
* colormap([0,0,0]); set(qcf,'color','w'); axis equal; axis off;
* figure(2); sur(Phi); exals({-12,12});
* axis equal; axis({0,nelx,0,nely,-12,12}); view(3);
* figure(3); subject(2,1,1); plot(comp(1:iT),'-'); title('Compliance');
* subplot(2,1,2); plot(vol(1:iT),'-'); title('Volume fraction');
* SS% CONVERGENCE CHECK
* if iT=nRelax && abs(vol(iT)-volfrac)/volfrac<1e-3 &&...
* all(abs(comp(iT)-comp(iT-9:iT-1))/comp(iT)<1e-3) break;
* end
* SS% LAGRANGE MULTIPLIER
* if iTsnRelax
* lag = mu*(vol(iT)-vol(1)+(vol(1)-volfrac)*iT/nRelax);
* else
* lag = lag+gamma*(vol(iT)-volfrac);
* gamma = min(gamma+0.05,5);
* end
* SS% LEVEL-SET FUNCTION EVOLUTION
* [zeros(nely+1,1) Phi(:,1:end-1)]);
* [zeros(1,nelx+1); Phi(1:end-1,:)]);
* gradphi = sqrt((pQSPX).*2+(pQPY). *2);
* indexDelta = (abs(Phi(1))*beta);
* DeltaPhi = zeros(size(Phi));
* DeltaPhi(indexDelta) = 0.75/Delta*(1-Phi(indexDelta). *2/Delta*2);
* eleComp = reshape(eleComp,nely,nelx);
* eleCompLR = [eleComp(:,1),eleComp]+[eleComp,eleComp(:,end)];
* nodeComp = [(eleCompLR,eleCompLR(end,:)]+[eleCompLR(1,:),eleCompLR)]/4;
* B = (nodeComp(:)/median(nodeComp(:))-lag). *DeltaPhi(:)*Delta/0.75;
* Alpha = Alpha+dt*B;
* Alpha = Alpha/mean(gradPhi(unique(eleNode((eleVol<1 && eleVol>0),:))));
* Phi = reshape(Alpha,nely+1,nelx+1);
* Phi(Boundary)=min(Phi(Boundary),-1e-9);
* end

## Appendix B Modifications of the Optimization Algorithm

There are changes to achieve model Figure 7b.