Topology Optimization: A Compact Review and Comparation on Level Set Technique Development in Single and Multi-material

Problem of Unified Theoretical Framework

In this study, papers using the level set method are selected based on criteria such as convergence speed, computation cost, optimizers that accelerate convergence or enable void creation, and volume fraction, to form a unified theoretical framework for the entire analysis. This approach distinguishes this review from others (Feng et al., 2018; Gibou et al., 2018; J. Liu et al., 2018) in its synthesis of studies, as it presents analyses based on results according to the four criteria mentioned above. These criteria were selected for the analysis for the following reasons:

1. Convergence speed is often one of the factors reflected in the results of topology optimization research. Although not explicitly stated, some studies use the number of convergence iterations to demonstrate that a new method is more optimal than a previous one . However, the total time required for convergence has not been analyzed as an important parameter, which must be related to the total number of iterations. There are two reasons that necessitate this correlation: firstly, optimization studies for the smoothness of the optimized shape or for multi-material may require up to 200, 300, 500, or even 1500 iterations (O. Sigmund, 2007; Sugimoto et al., 2024; F. Wang et al., 2011; Zhao & Zhang, 2021; R. Zheng et al., 2024) . Secondly, when comparisons are needed, the time per iteration is a better comparative metric than the number of iterations only, because a small number of iterations can result in a very long time.

1. Table 1. Visualization of the review purposes on the level set method for single material.

1. One challenge in topology optimization research is the hardware requirements of computer configurations (Wegert et al., 2025; S. Xu et al., 2022; Yago et al., 2022) . This research focuses on studies that can work efficiently on personal computers, which are suitable for most new researchers. This efficiency not only achieves optimal performance on the current problem, but it also assesses the ability to develop complex structures, such as optimizing 3D structures or multi-material problems, while maintaining optimal performance on personal computers.

2. The presence of voids increases the complexity of the optimized geometry during topology optimization. In other words, creating many voids is considered optimal for the method. Some level set approaches, including conventional level set method (Wei & Paulino, 2020), find it difficult to create new voids, so optimizers such as volume integration schemes are needed to generate them in the initial step, thereby speeding up convergence (Wei et al., 2018; Wei & Paulino, 2020) .

3. The volume fraction is one of the final results of the optimization process. Most studies reported this value only as numerical values that converged, such as 0.5 or 0.45, without deeper discussion. Only a few studies show that this result depends quite heavily on intermediate elements and the Lagrangian augmentation method (Y. Cui et al., 2023), which is widely used in research. For single-material studies, the constrained volume is usually large, for example, 0.5, so the volume fraction result is usually very close to the constrained volume. However, when optimizing across multiple materials, each material has a very small volume, leading to a significant discrepancy between the constrained and optimized volumes. This summary will demonstrate this discrepancy in multi-material optimization studies, where authors also avoid providing these values.

From the optimization perspective, studies that focus on the level set method for single-material to improve convergence speed, reduce computational cost, enhance the optimizer, and improve volume fraction are considered. Each purpose is presented in Table 1, and the studies corresponding to each purpose are detailed below.

Radial Basis Function

The work of Wei et al. (2018) for Radial Basis Function (RBF) is chosen as the reference benchmark in this study. The first reason is that Wei et al. (2018) introduced 88 lines of MATLAB code for educational purposes. Second, this work used a volume integration scheme to increase optimization speed by starting with a lower volume fraction (see Figure 2).

Wei and Paulino (2020) improved the RBF method by using unstructured polygonal finite element models, enabling optimization of more complex shapes, higher accuracy, and clearer material-void boundaries without the need for additional smoothing methods. Wei and Paulino (2020) also replaced the Hamilton-Jacobi equation with a system of ordinary differential equations (ODE) to update the level set function. However, the results showed no significant improvement in convergence speed.

Andrade and Santos (2022) used the Hilbertian extension to update the sensitivities along with the RBF method. In this method, a velocity field is created, reducing the dependence on the mesh of the previous method (which only used the natural normal velocity field). Andrade and Santos’s (2022) research also used the Ersatz material approach (Laurain, 2018) to achieve smoother density element values. Specifically, the Young's modulus of the voids is assigned a very small value (i.e., 1e-4) rather than 0, and the material's Young's modulus is set to approximately 1, creating a continuous cover zone around the boundary and resulting in a more seamless boundary. While this offers more advantages than the previous approach, using the Hilbertian extension increases computation costs.

In 2024, Sobczak and Sokół used RBFs and Hilbertian descent direction to approximate the level set function and performed their analysis using the Wolfram Mathematica code platform. The novel contribution of Sobczak and Sokół (2024) was the application of the analytical exponential formula in level set problems to reduce the number of iterations required for convergence. However, the study did not directly compare computation time with other studies due to differences in coding environments. In addition, this study focused on capturing the optimal shape, ignoring the presence of some free solid elements in the void domain.

In the same year, J. Zheng et al. (2024) published an improvement to the traditional RBF method that significantly improved computational performance. J. Zheng et al.’s (2024) method combined RBF with finite difference to create a meshless method based on local differentiation construction, namely RBF-FD-PLSM, thereby controlling the computational matrix.

Improvement of the Optimizer (RBF Approach)

The Importance of Multi-material Optimization

The advantages of topology optimization methods have been described and demonstrated in studies of single materials. To apply these advantages to practice, problems with multi-materials are necessary. With a large number of outstanding studies, multi-material optimization using the SIMP technique has achieved significant results. Using optimization for multi-materials can achieve better strength-to-weight ratios, optimize costs when combining different materials for the same structure (Bohrer & Kim, 2021; Doan et al., 2020) , or be applied to layered-material structures (Tajs-Zielińska & Bochenek, 2021). In practical applications, such as in the automotive field, Shah et al. (2023) developed a prototype design for a lightweight electric all-terrain vehicle (eATV) chassis by combining titanium, aluminum, and carbon fiber sheet molding compound. In the aerospace field, critical aircraft structures such as the aircraft pylon, engine bracket, and front fuselage have also been studied for topology optimization applications (Zhu et al., 2016b).

In optimization involving graded porous materials, studies by Bruggi et al., (2022), Wu et al. (2018), and Schmidt et al. (2019) show that the distributions of the solid phase of the material, a graded porous phase, and a void are a form of multi-material optimization, where the graded porous phase value lies between 0 (void) and 1 (solid phase). In practice, graded porous materials are used for bone screws or bone-like structures (Wu et al., 2018; Zhao & Zhang, 2021) to improve the mechanical and biological properties of porous osseointegrated screws (Yang et al., 2024).

The development of level set methods for multi-material optimization emerged nearly a decade after SIMP, in 2004 (M. Y. Wang & Wang, 2004) rather than 1997 (J. Liu et al., 2018; 0 Sigmund & Torquato, 1997; R. Zheng et al., 2024) , respectively. Nevertheless, research and applications of level set methods have made important progress, particularly MM-LS techniques, which offer high customization and the potential to develop new approaches.

The section on multi-materials will discuss methods for distributing different material types within the design domain by controlling level set functions. There are many different ways to control level set functions, and this study will classify and analyze these control methods.

“Color” Level Set Model Method (CLSM)

As shown in Figure 3, a method called CLSM (M. Y. Wang & Wang, 2004), defines up to four phases for the distribution of two materials defined by two level set functions. The intersection of 2 materials is defined to exist when ${\varphi }^{\langle 1\rangle }>0$ and ${\varphi }^{\langle 2\rangle }>0$, material 1 appears when ${\varphi }^{\langle 1\rangle }>0$ and ${\varphi }^{\langle 2\rangle }<0$, material 2 appears when ${\varphi }^{\langle 1\rangle }<0$ and ${\varphi }^{\langle 2\rangle }>0$, and voids appear when both ${\varphi }^{\langle 1\rangle }$ and ${\varphi }^{\langle 2\rangle }$ negative. Then, material 1 is distributed in the intersection area. Therefore, for each $M$ of material, the number of phases is $2^M$.

The fact that the number of phases is $2^M$, where $M$ is the number of materials to be optimized in the computational domain, leads to a significant increase in the number of phases as the number of materials $M$ increases. This makes it difficult to control the phases, or adds pressure on the computational equation.

Due to the disadvantage of increasing the number of materials to be optimized, the CLS method has not been widely adopted, and most studies stop at 4 phases, including voids (corresponding to 2 level set functions). Nevertheless, some studies have taken advantage of the CLS method's such as ability to generate 3 material phases from 2 level set functions, as well as its smooth phase transitions at boundaries and prevent the material overlap, such as Zhou et al. (2022) researched on multiphase auxetic metamaterials, Ghasemi et al. (2018) researched on flexoelectric composites, and P. Liu et al. (2016) used the union phase to define the bonding interface.

Multi-material Level Set (MM-LS) Method

To overcome the shortcomings of the CLSM method, Wang et al. (2015) introduced a completely new approach in which the total number of materials and void in the computational domain to be optimized equals the total number of phases, M+1, where M is the total number of materials. Thus, the number of phases for the material (excluding voids) is also equal to the number of level set functions. This approach makes controlling the level set functions easier and makes increasing the number of materials to be optimized more feasible. Furthermore, the fact that the number of phases for the material (excluding voids) equals the number of level set equations also facilitates the development of new methods based on MM-LS, as analyzed below.

First and original, the expression for cases with up to three materials in Wang et al.’s (2015) research is represented in Figure 4. The ${\varphi }^{\langle 1\rangle }$ is used to define the solid or void phase; therefore, in any area where ${\mathrm{\varphi }}^{\langle 1\rangle }<0$, the void has appeared. The ${\varphi }^{\langle 2\rangle }>0$ is used to present material 2 in the solid phase when ${\varphi }^{\langle 1\rangle }>0$. When ${\varphi }^{\langle 2\rangle }<0$, the rest of the solid phase is defined by material 1. Material 3 appears in the material 2 zone when ${\varphi }^{\langle 3\rangle }>0$. In case ${\varphi }^{\langle 3\rangle }<0$, the material 2 zone is not replaced.

According to Y. Wang et al. (2015), material occurrence was governed by level set functions. The main difference from the CLSM method was that a single level-set function (${\varphi }^{\langle 1\rangle }$) controlled the occurrence of all materials. Under the conditions for material occurrence, this level set function should always be positive; when it is negative, voids will appear. For consistency with other methods, Table 6 presents the conditions applied for three materials: code A indicates material occurrence, while code B indicates voids.

1. Table 2. Summary table of convergence rates of several older studies using the Radial Basis Function based level set method.

2. Table 3. Computational cost of Radial Basis Function based level set method.

3. Table 4. Computational cost of Reaction-Diffusion Equation based level set method.

4. Table 5. Computational cost of the Velocity Field Level Set method.

5. Table 6. The condition of the materials and the void phases in 3-material case.

Material 1, material 2, and material 3 denote materials of increasing stiffness: softest, medium, and hardest, respectively. Figure 5 shows the order and appearance conditions of these materials in Y. Wang et al.’s (2015) study, where the hardest material (material 3) appeared first.

Nguyen-Tran et al.’s (2025) research applies the exact volume constraint (Y. Cui et al., 2023) with MM-LS. In the 3-material problem, material 3 appears under condition $\left\{{\varphi }^{\langle 1\rangle }>0,{\varphi }^{\langle 2\rangle }>0,{\varphi }^{\langle 3\rangle }>0\right\}$, and material 2 appears under condition $\left\{{\varphi }^{\langle 1\rangle }>0,{\varphi }^{\langle 2\rangle }>0,{\varphi }^{\langle 3\rangle }<0\right\}$. Material 1 occurs for all values of ${\varphi }_3$ $\left\{{\varphi }^{\langle 1\rangle }>0,{\varphi }^{\langle 2\rangle }<0,\forall {\varphi }^{\langle 3\rangle }\right\}$, while the void state occurs for all combinations of ${\varphi }^{\langle 2\rangle }$ and ${\varphi }^{\langle 3\rangle }$, but only when ${\varphi }^{\langle 1\rangle }$ is less than 0 $\left\{{\varphi }^{\langle 1\rangle }<0,\forall {\varphi }^{\langle 2\rangle },\forall {\varphi }^{\langle 3\rangle }\right\}$. The study by Nguyen-Tran et al. (2025) used the same material reconstruction conditions as in Y. Wang et al.’s (2015) original study. The difference is that the softest material (material 1) appeared first under the same A-0 and A-1 conditions, as shown in Figure 6.

Next, M. Cui et al. (2016) applied the RDE technique to the improved MM-LS method, using all level set functions simultaneously to control the occurrence of each material. In M. Cui et al.’s (2016) study, the application of level set functions for material reconstruction was modified. Unlike Wang's original study, in which all conditions in code A were used to reconstruct the materials, M. Cui et al.’s (2016) study adjusted the application to use only one condition for each material, as shown in Figure 7.

In this case, the condition A-2 in code A is also presented as void (see Figure 8). However, since condition A-2 still has ${\varphi }^{\langle 1\rangle }>0$, a material currently under condition A-2 can be transformed into conditions A-1 and A-3 if any level set function, ${\varphi }^{\langle 2\rangle }$ or ${\varphi }^{\langle 3\rangle }$, changes sign. This phase shift is weak, but it is still possible if the conditions defined by the sign of the level set function are met. The advantage of M. Cui et al.’s (2016) method is that the material reconstruction rate for each material is mathematically identical, making it suitable for more than three materials when each is defined by only one condition.

However, the initial appearance of the weakest material can make it difficult to replace voids later in practice. To explain this disadvantage in detail, initially, the softest material appears throughout the entire domain. When hard materials replace the load-bearing areas, the softer materials must redistribute to support the structure. Nevertheless, these areas may already have been replaced by voids (condition code B). This makes it difficult for the softest material to replace them, because it increases compliance. There is evidence in Cui's results that the softest material appears as small dots around the harder material.

The approach in Wang et al.’s (2015) and Nguyen-Tran et al.’s (2025) method showed that all conditions with ${\varphi }_1>0$ are used to determine the material, and that there is always a preference for a material with more conditions to appear, with Wang et al. (2015) prioritizing the hardest material, while Nguyen-Tran et al. (2025) prioritizing the softest materials. This selection can lead to over prioritization when increasing the number of materials to be optimized, and also create differences in the allocation of conditions for material reconstruction. M. Cui et al.’s (2016) method can adjust a material's condition to address this problem. However, improvements are still needed in the order of material appearance and in the systematic selection of conditions.

Using an alternating active-phase formulation and a difference-set-based multi-material level set (DS-MMLS) design, Sha et al. (2021) presented a revolutionary multi-material level set topology optimization technique. Sha et al. (2021) developed a new approach by using the MM-LS method; however, by applying the alternating active-phase methodology, the challenge of multi-material topology optimization can be transformed into a succession of binary-phase sub-problems (Sha et al., 2021). The sub-problems help reduce the number of variables and the volume constraints (Sha et al., 2021). First, material 1 will appear in the design domain where ${\varphi }^{\langle 1\rangle }>0$. Next, material 2 will appear in the region where ${\varphi }^{\langle 1\rangle }<0\cap {\varphi }^{\langle 2\rangle }>0$, it must not overlap with the region where material 1 already exists. Then, according to the above principle, material 3 appears when ${\varphi }^{\langle 1\rangle }<0\cap {\varphi }^{\langle 2\rangle }<0\cap {\varphi }^{\langle 3\rangle }>0$, it must not overlap with the region of material 1 and material 2 before (see Figure 9).

Although Liu and Ma (2018) and Sha et al. (2021) also used the MM-LS method, their approaches shared some similarities with the CLSM method in defining the intersection region. Specifically, to optimize the level set using component-length-scale control on multi-material structures, Liu and Ma (2018) directly defined each level set function to represent each material phase, as shown in Figure 10, thereby freeing the dependent length-scale control. However, Liu and Ma’s (2018) research is limited to two materials, with materials 1 and 2 corresponding to level set functions 1 and 2, respectively, when they are positive. The intersection region between the two materials is defined as the intersection of two positive level set functions, assigned a stiffness value equal to 0.5 times the stiffness value of the softest material.

Recent studies on multi-material level-set-based topology optimization (MM-LS) show that the method has been widely adopted and is developing rapidly. Only Sha and Nguyen's study optimized more than 3 materials, specifically 5. Other studies, especially Kumar and Rakshit (2025), although introduced in the methods section as optimizing for up to 3 materials, reported results that optimized only for 2 materials (see Appendix B).

The common point of these methods is that they all use the exact Heaviside function to determine the elastic stiffness for the materials, where when the conditions are met, the exact Heaviside function will set the corresponding level set values ​​to 0 (void) and 1 (material), and assign them to the elastic stiffness values ​​through equations (2) to (7).

The Limitations of the Previous MM-LS and Future Works

The disadvantages are most evident when optimizing for multiple materials. First, the results do not show the total time to complete the iterations (except Kumar and Rakshit (2025)), so it is not possible to determine which approach is faster. The number of iterations varies across studies, depending on the author's purpose. Therefore, if the computation time per iteration were available, a comparison of convergence speeds between methods would be more appropriate.

The second concern is the volume fraction. Regarding volume fractions, materials within the multi-material computational domain often have low volume constraints, leading to larger discrepancies in volume fraction results when combined with intermediate elements and the augmented Lagrangian method. Specifically, previous studies did not provide precise results for the volume achieved after optimization; the results were only shown graphically (approximately). In particular, there were cases where the volume did not converge to the given constraint, such as Sha et al.’s (2021) study, which yielded a volume of approximately 0.28, while the constraint volume was 0.2. Similarly, in Kumar and Rakshit’s (2025) study, neither volume result converged to 0.25 as the constraint volume. Currently, only Nguyen-Tran et al. (2025) used the exact volume constraint to address this issue.

Analysis of the aforementioned multi-material optimization studies reveals that the results proposed in these previous studies have not been fully demonstrated by the previous authors, or may even be omitted or misjudged (Kumar & Rakshit, 2025b; Sha et al., 2021) . This creates a suitable research gap, where one or more new approaches are needed that fully represent computation speed (time per iteration), computational cost (to support the popular application of the method), and volume fraction (convergence criteria).

From these observations, it is clear that, in future research, the application of the MM-LS method in multi-material optimization still presents significant research gaps. First, the MM-LS method should be improved to use a single material per condition, and verified with more than three materials, as M. Cui et al.’s (2016) study failed to do. Next, the volume fraction value (accurate) should be improved compared to the constraint volume while still applying the Lagrangian augmentation method. Finally, the computation time value could be considered for presentation to develop the first comprehensive comparative study.

Comparative Analysis for Single Material Approaches

From the summary Tables 3 to 5 above, studies provide relatively complete results for the number of loops, but the total time required to complete them is not fully published. When comparing the time required to execute one loop, the results show that the study of Andrade & Santos (2022); J. Zheng et al. (2024) - RBF; C. Wang et al. (2022) – RDE; and Zhuang et al. (2024) are the most optimal in iteration number and convergence time. Note that the device's computing power affects the computation time per iteration, but, interestingly, the results of these research show that the personal computer can be used to achieve the optimal results.

The Reaction-Diffusion Equation (RDE) (Otomori et al., 2015), Velocity Field Level set (VFLS) (Y. Wang & Kang, 2021), and Radial Basis Function (RBF) (Wei et al., 2018) methods all provide useful, open research code and are all executable on personal computers, making them more accessible to a wider range of researchers. However, the studies selected as references for each approach do not fully demonstrate the values ​​that this review aims to achieve. Furthermore, the reason for selecting these studies as references for each approach is that they present the original methodology. From there, this review paper work will continue to calculate the time per iteration for studies Wei et al. (2018), Otomori et al. (2015), and Y. Wang & Kang (2021) in section Results to compare them on the same system.

The RBF and VFLS approaches both solve the Hamilton-Jacobi equation, creating symmetry holes in the initial design to improve convergence and reduce computation time. The RDE approaches, on the other hand, do not solve the Hamilton-Jacobi equation but instead solve the reaction-diffusion equation, which reduces computation time by eliminating the need for re-computation as in the Hamilton-Jacobi problem, and helps achieve the optimal result more naturally without creating holes in the initial computational domain.

In this comparison, to ensure a fair comparison of different approaches, the following criteria will be considered:

1. The AMD Ryzen 5 5600H (3.3GHz~4.2GHz) multicore CPU with 64 GB of RAM is used to run this program in MATLAB R2023b (License 40908710) for the numerical examples. 64 GB of RAM can provide sufficient memory for all numerical benchmarks.

2. Mesh parameters: all numerical computation cases are performed with the same finite element meshes. 5 mesh parameters are compared, including $60\times 30$, $80\times 64$, 12$0\times 60$, $160\times 128$, and 24$0\times 120$.

3. Volume fraction: Each approach optimizes the problem with a constrained volume of 0.5.

4. Parameters for each approach: published studies confirm that the parameters selected for each method have been optimized using the algorithm; therefore, the parameters are retained, including each approach's convergence criteria.

5. Objective function: In Wang et al.’s (2021) study, to reduce the influence of design domain size differences on compliance results, Wang et al.’s (2021) used normalized compliance values. In the current study, the design domain dimensions are similar across methods, so the comparison results are minimal compliance. Additionally, Otomori et al. (2015) code has been corrected by the current research for the incorrect placement of force points in the finite element model, resulting in accurate compliance.

6. The compliance criteria are kept by the original works, which may lead to differences in the total number of iterations required for convergence.

7. Therefore, computation speed will be compared using time per iteration. The stability of each method is determined by the stability of the iteration number results due to different mesh parameters.

The cantilever beam with the design domain as shown in Figure 2(a) is optimized using RBF, RDE, and VFLS for 2D compliance minimization problems. The 5 mesh parameters and results are presented in Table 9, whose volume fraction is set to 50%, a concentrated force $F=1$, the Poison's ratio is $\nu =0.3$.

Regarding convergence speed, for a $60\times 30$ mesh parameter, all three methods yield very fast per-iteration computation times. The computation time per iteration increases with increasing mesh size. However, the RDE method shows a significant improvement over the other two methods, requiring 108 seconds per iteration at a $240\times 120$ mesh size, whereas the RDE and VFLS methods require less than 1 second per iteration at this mesh size. Figure 11 illustrates the trend in computation time for the three approaches across 5 increasing levels of mesh parameters, all tested on the same personal computer.

The convergence iteration results indicate that RBF and RDE methods are stable, with little difference in iterations across mesh sizes. On the contrary, the VFLS method relies on parameters, leading to fast computation but unstable convergence. This conclusion is based on experiments using problems from Y. Wang & Kang’s (2021) study, which also showed highly parameter-dependent results.

1. Table 7. The comparison results of cantilever beam numerical example in 5 mesh parameters using 3 level set approaches.

The methods also showed different computational resource usage. Proportional to computation time, RAM usage increased with increasing mesh size. The RDE and VFLS methods only required between 50MB and 200MB of RAM for the smallest to largest mesh sizes in Table 7. The RBF method, on the other hand, used a lot of RAM, 400MB for the smallest mesh size, and increased rapidly to 42.6GB for the largest mesh size. Therefore, in the comparison presented in this review, the computer system used could not compute mesh sizes larger than $240\times 120$ due to exceeding the available RAM capacity.