## 2 Method
### Basic Definitions  
In topology optimization, the design domain is typically subdivided into elements by appropriate meshing. Figure 2 visualizes the elements (with one element hatched) and nodes. Figure 2: Design space overview with elements, nodes and dimension \(d\) (square case). In this work, we examined only square meshes with equal numbers of rows and columns. However, this method can be used for non-square and three-dimensional geometries. The total number of elements in the 2D case is as follows:  
\[d=d_{x}d_{y}, \tag{3}\]  
where \(d_{y}\) is the number of rows and \(d_{x}\) the number of columns (see Figure 2). In the square case, the number of rows and columns are equal: \(d_{x}=d_{y}=d\). The \(d^{2}\) design variables \(x_{i}\) {\(i=1,\ldots,d^{2}\)}, termed density values, scale the contributions of the single elements to the stiffness matrix. The density has a value of one when the stiffness contribution of the element is fully preserved and zero when it disappears. The density values are collected in a vector \(\mathbf{x}\). In general, the density values \(x_{i}\) are defined in the interval \([0,1]\). In order to prevent possible singularities of the stiffness matrix, a lower limit value \(x_{\text{min}}\) for the entries of \(\mathbf{x}\) is set as follows [2]:  
\[0<x_{\text{min}}\leq x_{i}\leq 1,\ \ \ i=1,2,\ldots,d^{2}. \tag{4}\]  
The vector of design variables \(\mathbf{x}\) can be transformed to a square matrix \(\mathbf{X}_{\text{M}}\) of order \(d\) by using the \(\mathcal{R}_{2d}\) operator:  
\[ \begin{array}{l}\mathcal{R}_{2d}(\mathbf{x})=\mathbf{X}_{\text{M}}=(x_{ \text{M}ij})\\ \text{where }x_{\text{M}ij}=x_{i+d(j-1)},\\ i=1,2,\ldots,d,\ \ \ j=1,2,\ldots,d.\end{array}  \tag{5}\]  
Although a binary selection of the density is desired (discrete TO, material present/not present), values between zero and one are permitted for algorithmic reasons (continuous TO). To get closer to the desired binary selection of densities, the so-called penalization can be used in the calculation of the compliance. The penalization is realized by an element-wise exponentiation of the densities by the penalization exponent \(p>1\)[26]. The arithmetic mean of all \(x_{i}\) defines the degree of filling of the geometry as follows:  
\[M_{\text{is}}=\frac{1}{d^{2}}\sum_{i=1}^{d^{2}}x_{i}. \tag{6}\]  
The target value \(M_{\text{tar}}\) is the degree of filling that is to be achieved by the predictor. Figure 3 shows an overview of the processing of boundary conditions. In this figure as well as in following figures, the symbols represent the kinematic boundary conditions (structural supports) and the vectors of the static boundary conditions (applied forces). Figure 3: Matrix representation of (**a**) kinematic boundary conditions (structural supports) and (**b**) static boundary conditions (applied forces). The kinematic boundary conditions are stored in two \((d_{\text{inp}}+1)\times(d_{\text{inp}}+1)\) Boolean matrices \(\mathbf{R}_{\text{k,x}}\) and \(\mathbf{R}_{\text{k,y}}\). The entries of \(\mathbf{R}_{\text{k,x}}\) are set to one if the x-component of the displacement in the corresponding node is fixed, and to zero otherwise. Analogously, the entries of \(\mathbf{R}_{\text{k,y}}\) are set according to the fixed y-components of the displacements. Both matrices can be transformed into vectors with the \(\mathcal{R}_{1d}\) operator, which is the inverse of the operator \(\mathcal{R}_{2d}\), and then arranged in sequence so that the vector \(\mathbf{r}_{\text{k}}\in\mathbb{R}^{\left(2d_{\text{inp}}+1\right)^{2}}\) is created. Analogous to the kinematic boundary conditions, two \((d_{\text{inp}}+1)\times(d_{\text{inp}}+1)\) matrices \(\mathbf{R}_{\text{s,x}}\) and \(\mathbf{R}_{\text{s,y}}\) are firstly built on the basis of the static boundary conditions. The x- and the y-components of the applied forces are placed, respectively, into the matrices in correspondence to their magnitude, while the remaining entries are set to zero.