## 3 Matlab Implementation
### Iterative Solution (Lines 54-137)
**(b)** For three materials, the compliance is 0.268, where blue, green, and yellow are the materials with \(\mathrm{E}=1\), \(\mathrm{E}=2\), and \(\mathrm{E}=5\), respectively. **(c)** For five materials, the compliance is 0.242, where blue, cyan, green, orange, and yellow are the materials with \(\mathrm{E}=1\), \(\mathrm{E}=2\), \(\mathrm{E}=3\), \(\mathrm{E}=4\), and \(\mathrm{E}=5\), respectively. variables equal to 1. Hence, it will drive the optimization model convergence to a clear result of material 1 and material 2 to be voids or solids. To further explore the effect of parameter \(p\) on the optimization results, the topology optimization of three materials with the half MMB beam using different P-norm values is also studied. The values of parameter \(p\) are set as follows:

Case 1: \(p=1\);

Case 2: \(p=6\);

Case 3: \(p=18\);

Case 4: \(p=32\);

The [PERSON]'s modulus of each material is set as par.E(1) = 1; par.E(2) = 2; par.E(3) = 5. The volume constraint is set to 0.5/3 for each material, and other parameters are the same as those in Section 4.1. The optimized structures are shown in Figure 5. It can be found that there are mixed materials at the boundary of the optimized structure when \(p=1\) because that p-norm function cannot approximate the maximum of the design variable of Equation (4). By increasing the value of parameter \(p\), the boundary of each material becomes clear for the optimized structure, and the compliance of the optimized structures is comparable. However, it will influence the stability of the optimizer with a very large value of parameter \(p\). Therefore, we choose \(p=6\) in the following examples to balance the structural stiffness and the stability of the optimizer [29]. Figure 4: The function diagram of interpolated Young’s modulus with respect to parameter \(p\). Figure 5: The optimized structure with different values of parameter \(p\) for topology optimization of three materials, where blue, green, and yellow are the materials with E = 1, E = 2, and E = 5, respectively. ## 5 Extension to 3D

The Matlab implementation described in Section 3 is remarkably easy to extend to 3D problems, which can be found in Appendix B. ### 3D Matlab Implementation

The 3D Matlab implementation is verified using a cantilever beam example here. The boundary conditions and loading conditions for the cantilever beam are shown in Figure 4, and the corresponding Matlab code is defined in lines 11-24. A notable modification compared to the 2D Matlab code is the change in the discrete unit element from four-node rectangle elements to eight-node cubic elements. The element stiffness matrix lk is organized in lines 25-40, and the matrix edofMate containing the node IDs for each element is obtained in lines 41-44. The iterative details for the three-material cantilever beam can be obtained by calling the code in lines 157-163. The final design is displayed based on the code top3d.m proposed by [PERSON] [20] in lines 164-200, and the optimized structures with multiple materials are shown in different colors. ### 3D Numerical Examples

#### 5.2.1 Example 1: The Number of Materials

The cantilever beam is taken as an example for the topology optimization of a 3D continuum structure with two, three, and five materials, respectively. The design domain and boundary conditions of the 3D cantilever beam are shown in Figure 6. The left part of the design domain is fixed, and a distributed load \(F=-1\) is applied at the lower right edge. The design domain is discretized into 80 \(\times\) 40 \(\times\) 20 cubic elements with Poisson's ratio \(\mu\) = 0.3. The initial density of each material is set to 0.5. All other parameters are the same as those in Section 4.1. The optimized structure and the converged history for topology optimization of different numbers of materials are shown in Figure 7. It can be found that the convergence histories of the compliance are very similar under different numbers of materials, and they converge rapidly in the first 50 iterations, which is similar to the 2D structures.