## 3 Results and Discussion
### Modifications of the Topology Optimization Algorithm
The rule introduced above has an exception. Consider the boundary conditions for the case shown in Figure 9(a), left. The shown body has three special points. There are two points where the body touches rollers and the point where force \(F\) is applied. Let us suppose that another body produces this force. Therefore, near these three points, we have a region where we cannot easily separate optimized body and external bodies. These points always contain a portion of material. Thus, Equation (19) does not hold for them. One implementation issue should be noted (see the MATLAB code in the Appendix A). Consider the initialization of the bilinear spline defining the body shown in Figure 9(a), right. It is a meaningful option of the body design, but it has holes at the boundary of \(D\). Small negative values are used in Equation (19) instead of zero to allow the algorithm to go over similar shapes. This modification of Equation (19) is also useful in one more situation. During the optimization process, at some iteration \(i\), we can obtain the defining function (7) with all positive \(\alpha_{j}\). It is true for all points except ones mentioned in (19). Replacement of zero with some small negative value prevents the algorithm error at the normalization step (10). The above-mentioned modifications make the algorithm more robust. We compared the proposed modified topology optimization algorithm with the original one [20]. The laptop with the following characteristics was used for tests: Windows 10 Enterprise 64-bit operating system, x64-based processor Intel(R) Core(TM) i5-8250U CPU @ 1.60 GHz 1.80 GHz, and 8.00 GB (7.86 GB usable) of installed memory (RAM). Both algorithms, the original one and the modified one, were run 10 times with following initial datasets. The domain grid is as shown in Figure 9(c). The initial geometry is as shown in Figure 9(a), right. Boundary conditions are taken from Figure 9(a), left (see details in Appendix A). The desired volume, \(V_{0}\), was chosen as 45% of the volume of the entire \(D\). The [PERSON]'s modulus of stainless steel, \(E_{0}\), is assumed to be \(205 GPa\), and the [PERSON]'s ratio, \(\
u\), is assumed to be 0.27. Stainless steel is a material obtained of the powder used for 3D printing. The properties of material were studied in [25]. In the original algorithm all parameters of the updating scheme for \(f_{i}(x,y)\) were preserved, except the time step. The time step was changed to 0.2. These parameters were also used for the modified algorithm. Both algorithms exploit an isotropic material model at the FEM step. However, the research [25] showed that 3D printed parts have anisotropy. Nevertheless, the validation results (see Section 3.3) showed that the isotropic model of material is suitable for the considered cases. Results of the computational test are shown in Table 1. We can see that the modified algorithm is approximately twice as fast as the original one. It can be easily explained using different spline basis functions. When the bilinear spine is used instead of the MQ-spline,  
Figure 10: Optimization cases and initial assumptions: (**a**) the three-point bending beam case; (**b**) the beam is supported horizontally in two bottom points, and equal loads applied vertically in upper corners; (**c**) the cantilever beam case; (**d**) the beam is fixed along the bottom side, and the horizontal load linearly distributed along the left side of the beam; (**e**) the beam is fixed horizontally in two bottom points, and equal loads applied horizontally in upper corners  
the time-consuming operation of matrix inversion is not needed. Two obtained optimized bodies are shown in Figure 7a. Four more optimization cases shown in Figure 10 were considered to check the observed tendency. Figure 7 shows the results of optimization algorithms testing. Efficiency of the algorithms for these cases is compared in Table 2. The results show that the modified algorithm is faster. The difference in efficiency is much higher for finer grids of basis functions in the domain.