### Density-Based Topology Optimization  
Earlier in this section, the general problem was formulated based on the boundary-following mesh or Lagrangian model. In contrast, the Eulerian model (fixed mesh) maps the topology through a density function, \(\rho\), on the discretized design domain. Finally, this density distribution, \(\rho\), is fed to the optimization process as design variables. Such reformulation is called density-based topology optimization. Compared to the shape-based TO, the nodes of the structures' boundaries are the optimization variables, and LSTO considers the parameterized LSFs to be the design variables for the optimization process. Thus, the structure's boundaries are at least retained in the optimization process, even though it has blurred out the discretized setting of the continuum design space when utilizing density-based mapping. However, when elements or nodes of the mesh are optimized, such as sizing variables, the sense of the exact illustration of the structure's boundaries is physically lost. Mathematically, it takes the following form:  
\[ \begin{split}\min_{\rho}&:\Phi(\mathbf{U},\mathbf{\rho} )\\ &:=\sum_{i}\int_{\Omega_{i}}f(\rho_{i},\mathbf{U}(\rho_{i}))d \mathbf{x}\\ \text{s.t.}&:G_{0}(\mathbf{\rho})=\sum_{i}\mathrm{v}_{i }\rho_{i}-|\Omega_{d}|\leq 0,\\ &:G_{j}(\mathbf{\rho},\mathbf{U}(\mathbf{\rho}))\leq G_{j}^{*},\quad j=1, \ldots,Q\\ &:\rho=\begin{cases}\epsilon,&(\text{ void })\\ 1,&(\text{ solid })\end{cases}^{\prime}i=1,\ldots,N\end{split}  \tag{6}\]  
The density distribution, \(\mathbf{\rho}\), denotes the design variable vector of length N. It is defined such that \(\rho=1\) if \(\mathbf{x}\in\omega\) and \(\rho=0\) otherwise. For regions where \(\rho=0\), solutions to the boundary value problem are not guaranteed, as the energy's bilinear form is not coercive. Thus, the density function is defined as \(\epsilon+(1-\epsilon)\mathbf{\rho}\), in which \(\epsilon\ll 1\) is the esatz parameter:  
\[ \begin{split}\mathbf{U}&=\inf_{\mathbf{U}}\Pi(\mathbf{ \rho},\mathbf{U}),\text{ s.t. }\\ \Pi(\mathbf{\rho},\mathbf{U})&=\int_{\Omega}[ \epsilon+(1-\epsilon)\mathbf{\rho}]W_{0}(\mathbf{U},\mathbf{x})\mathrm{d}\mathbf{ x}-\int_{\Gamma_{N}}\mathbf{\tau}\cdot\mathbf{U}dS\end{split}  \tag{7}\]  
Thus, the above formulation is a binary problem representing the void and solid regions of the structure, hence named discrete density-based topology optimization (DDTO). The well-known DDTO method is the Bi-directional Evolutionary Structural Optimization (BESO), as it presently defines and adopts techniques commonly used in continuous DTO approaches. Recently, [PERSON] et al. [58] combined the features of BESO and the sequential integer linear programming for discrete topology optimization. Interested readers can find comprehensive reviews on the BESO methods in [44; 59]. Another outlook on approaching the discretizing problem is using a genetic algorithm [60] that can find a \"global minimum\" and allow the handling of a discrete variable, but this always comes at a higher computational cost. Furthermore, [PERSON] [61] questions the usefulness of non-gradient approaches in TO. With that, we concluded the discussion on DDTO as, in practice, it is recommended, when formulating the TO problem, to assume the continuous density field together with (7), as portrayed below:  
\[ \begin{split}\min_{\mathbf{\rho}}&:\Phi(\mathbf{\mathrm{ U}},\mathbf{\rho})\\ &:=\sum_{i}\int_{\Omega_{i}}f(\rho_{i},\mathbf{\mathrm{U}}(\rho_{i}))d \mathbf{\mathrm{x}}\\ \text{s.t.}&:G_{0}(\mathbf{\rho})=\sum_{i}v_{i}\rho_{i} -|\Omega_{d}|\leq 0,\\ &:G_{j}(\mathbf{\rho},\mathbf{\mathrm{U}}(\mathbf{\rho}))\leq G_{j}^{*},\quad j =1,\ldots,Q\\ &:0\leq\rho_{min}\leq\rho\leq 1,\quad i=1,\ldots,N\end{split}  \tag{8}\]  
The above Formulation (8) is a broadly received idea in the TO of continuum structures that utilizes continuous density design variables instead of binary density variables, thus enabling the use of gradient-based information.