## Appendix B Multimaterial3d.m  
1 %% LOOP PARAMETERS
2 clear;
3 NumIter = 300;delta=1e-9;
4 nelx=60; nely=30; nelz = 15;
5 penal=3; rmin=5; NMaterial=3;
6 NDV=nelx*nely*nelz*NMaterial; NDV3=nelx*nely*nelz;
7 movelimit=0.1; beta =1;
8 %% MATERIAL PROPERTIES
9 par.E(1)=1; par.E(2)=2; par.E(3)=5;
10 nu = 0.3;
11 %% USER-DEFINED LOAD DOFs
12 il = nelx; jl = 0; kl = 0:nelz;
13 loadind = kl*(nelx+1)*(nely+1)+il*(nely+1)+(nely+1-jl);
14 loaddof = 3*loadind(:) - 1;
15 %% USER-DEFINED SUPPORT FIXED DOFs
16 [jf,kf] = meshgrid(1:nely+1,1:nelz+1);
17 fixedind = (kf-1)*(nely+1)*(nelx+1)+jf;
18 fixeddofs = [3*fixedind(:); 3*fixedind(:)-1; 3*fixedind(:)-2];
19 %% DEFINE LOADS AND SUPPORTS (HALF MBB-BEAM)
20 nele = nelx*nely*nelz;
21 ndof=3*(nelx+1)*(nely+1)*(nelz+1);
22 F=sparse(loaddof,1,-1,ndof,1);
23 U=zeros(ndof,1);
24 freedofs=setdiff(1:ndof,fixeddofs);
25 %% PREPARE FINITE ELEMENT ANALYSIS
26 ke=1/(1+nu)/(2*nu-1)/144 *([-32;-6;-6;8;6;6;10;6;3;-4;-6;-3;-4;-3;-6;10;3;6;8;3;3;4;-3;-3;-32;-6;-6;-4;-3;6;10;3;6;8;6;-3;-4;-6;-3;3;4;-3;3;8;3;...
27 3;10;6;-32;-6;-3;-4;-3;-3;4;-3;-6;-4;6;6;8;6;3;10;3;3;8;3;6;10;-32;6;6; -4;6;3;10;-6;-3;10; -3;-6;-4;3;6;4;3;3;8;-3;-3;-32;-6;-6;8;6;-6;10;3;3;4;...
28 -3;3;-4;-6;-3;10;6;-3;8;3;-32;3;-6;-4;3;-3;4;-6;3;10;-6;6;8;-3;6;10;-3;3;8;-32;-6;6;8;6;-6;8; 3;-3;4;-3;3;4;-3;6;10;3;-6;-3;2;6;-6;-4;3;3;8;-3;...
29 3;10;-6;-3;-4;6;-3;4;3;-32;6;3;-4;-3;-3;8;-3;-6;10;-6;-6;8;-6;-3;10;-3;2;6;-6;4;3;-3;8;-3;3; 10;-3;6;-3;6;-32;6;-3;10;-6;-3; 30 -32;3;-6;10;3;-3;8;6;-3;10;6;-6;8;-32;-6;6;8;6;-10;-6;3;-4;-6;3;-32;6;-6;-4;3;6;10;-3;6;8; -6;-32;6;3;-4;3;3;6;-4;-3;32;-6;-6;-4;6;3;10;-6;3; 31 -32;6;-6;8;-6;-6;10;-3;-32;-3;6;-4;-3;3;3;4;-32;-6;-6;8;6;6;-32;-6;-6;-4; -3;-32;-6;-3;-4;-32;-6;-3;-4;-32;
32 4;0;0;0;24;-12;-12;0;-12;0;0;-12;12;12;48;0;24;0;0;0;-12;-12;-24;0;-24;0;0;24;12;-12;12;0; -12;0;-12;-12;0;0;-12;12;12;48;0;24;0;0;0;-12;12;12;0;24;0;-24;-24;0;0;24;12;-12;12;0; -12;0;-12;-12;0;0;-12;12;12;0;24;0;-24;-24;0;...
33 0;-12;-12;0;0;-12;-12;0;-12;48;0;0;0;-24;0;-12;0;12;-12;12;0;0;0;-24; -12;-12;-12;0;0;48; 0;24;0;-24;0;-12;-12;-12;-12;12;0;0;24;12;-12;0;...
34 0;-12;0;48;0;24;0;-12;12;12;0;-12;12;24;-24;0;12;0;-12;0;-12;0;0;-12;48;0;0;0;-24;24;-12;0;0;-12;12;0;0;-12; 12;-12;0;0;-24;-12;-12;0;48;0;24;0;0;0;-12;0;-12;12;-12; 0;0;-12;12;0;0;24;48;0;12;-12;0;0;-12;12;0;-12;12;-12; 36 -12;0;0;-24;48;-12;0;-12;0;0;-12;0;12;-12;-24;24;0;48;0;0;0;-24;24;-12;0;12;0;24;0;48;0;24;0; 0;0;-12;12;-24;0;24;48;-24;0;0;0;-12;12;-12;0;-24; 37 0;48;0;0;0;-24;0;-12;0;-12;48;0;24;0;-12;12;48;0;-24;0;12;-12;-12;48;0;0;0;-24;-24;48;0; 24;0;0;48;0;0;48;0;48]);
38 lk(tril(ones(24))==1)=ke';
39 lk = reshape(lk,24,24 );
40 lk =lk +lk'- diag(diag(lk ) );
41 Num_node = (1+nely)*(1+nelx)*(1+nelz);42 nodenrs = reshape(1:Num_node,1+nely,1+nelx,1+nelz);
43 edofVec = reshape(3*nodenrs(1:end-1,1:end-1,1:end-1)+1,nelx*nely*nelz,1);
44 edofMat = repmat(edofVec,1,24)+ repmat([0 1 2 3*nely + [3 4 5 0 1 2] -3 -2 -1 3*(nely+1)*(nelx+1)+ [0 1 2 3*nely + [3 4 5 0 1 2] -3 -2 -1]],nele,1);
45 iK = reshape(kron(edofMat,ones(24,1))',24*24*nele,1);
46 jK = reshape(kron(edofMat,ones(1,24))',24*24*nele,1);
47 %% PREPARE FILTER
48 iH = ones(nele*(2*(ceil(rmin)-1)+1)^2,1);jH = ones(size(iH));sH = zeros(size(iH));
49 k = 0;
50 for k1 = 1:nelz
51 for i1 = 1:nelx
52 for j1 = 1:nely
53 e1 = (k1-1)*nelx*nely + (i1-1)*nely+j1;
54 for k2 = max(k1-(ceil(rmin)-1),1):min(k1+(ceil(rmin)-1),nelz)
55 for i2 = max(i1-(ceil(rmin)-1),1):min(i1+(ceil(rmin)-1),nelx)
56 for j2 = max(j1-(ceil(rmin)-1),1):min(j1+(ceil(rmin)-1),nely)
57 e2 = (k2-1)*nelx*nely + (i2-1)*nely+j2;
58 k = k+1;
59 iH(k) = e1;
60 jH(k) = e2;
61 sH(k) = max(0,rmin-sqrt((i1-i2)^2+(j1-j2)^2+(k1-k2)^2));
62 end
63 end
64 end
65 end
66 end
67 end
68 H = sparse(iH,jH,sH);
69 Hs = sum(H,2);
70 %% INITIALIZE ITERATION
71 loopbeta=0;
72 x=ones(NDV,1)*0.5;
73 for i=0:NMaterial-1
74 xTilde(1:nely,1:nelx,1:nelz,i+1)=reshape(x(1+i*NDV3:NDV3+i*NDV3),nely,nelx,nelz);
75 end
76 xPhys = (tanh(beta*0.5) + tanh(beta*(xTilde-0.5)))/(tanh(beta*0.5) + tanh(beta*(1-0.5)));
77 %% MMA PARAMETER INITIALIZE
78 a0=0; a=zeros(NMaterial,1);cc=1.06*ones(NMaterial,1); d=0*ones(NMaterial,1);
79 xmin=1e-3*ones(NDV,1); xmax=1*ones(NDV,1);
80 xold=x; xolder=xold;low=0; upp=1;
81 %% START ITERATION
82 for iter=1:NumIter
83 loopbeta=loopbeta+1;
84 %% FE-ANALYSIS
85 KE0=lk; p=6;
86 for i=1:NMaterial
87 KE(:,:,i)=par.E(i)*KE0;
88 end
89 KE_test=zeros(24*24,nely*nelx*nelz); g=zeros(NMaterial,1);
90 r1=xPhys(:,:,:,1); r2=xPhys(:,:,:,2); r3=xPhys(:,:,:,3);
91 for j=1:NMaterial
92 post_x(1:nely,1:nelx,1:nelz,j)=(xPhys(:,:,:,j).