## 2 Software Structure and Methods
### Topology Optimization
Define the number of grid elements in the rectangular domain \(nelx\times nely\); Step 2. Initialize coefficients of the bilinear spline; Step 3. Define boundary conditions for FEM analysis; Step 4. Initialize parameters of the optimization loop; Step 5. Perform FEM analysis of the domain; Step 6. If the algorithm converges, then quit, else update \(\alpha_{i}\) and go to Step 5. ```  
**Algorithm 1**The proposed topology optimization algorithm. The following equation is used to update the free-form function (3):  
\[ \begin{split}& V_{n}=V\cdot(-\frac{\
abla f}{|\
abla f|}),\\ &\frac{\partial f}{\partial t}-V_{n}|\
abla f|=0,\end{split}  \tag{8}\]  
where \(t\) is the pseudo time representing the evolution of the free-form function, and \(V_{n}\) is the outward normal velocity chosen based on the shape derivative of the optimization problem. The equation for updating \(\alpha_{i}\) is as follows:  
\[\alpha_{i}(t_{j+1})=\overline{\alpha}_{i}(t_{j})+\Delta t\Phi^{-1}B(\overline {\alpha}_{i}(t_{j}),t_{j}), \tag{9}\]  
where  
\[ \begin{split}&\overline{\alpha}_{i}(t_{j})=\frac{\alpha_{i}(t_{j})}{ mean(|f(x_{1}^{1},x_{2}^{1},...,x_{n}^{1},t_{j})|,|f(x_{1}^{2},x_{2}^{2},...,x_{n}^{2},t_ {j})|,...,|f(x_{1}^{r},x_{2}^{r},...,x_{n}^{r},t_{j})|},\\ &\Phi=\begin{pmatrix}\phi_{1}(x_{1}^{1},x_{2}^{1},...,x_{n}^{1} )&\cdots&\phi_{N}(x_{1}^{1},x_{2}^{1},...,x_{n}^{1})\\ \vdots&\ddots&\vdots\\ \phi_{1}(x_{1}^{N},x_{2}^{N},...,x_{n}^{N})&\cdots&\phi_{N}(x_{1}^{N},x_{2}^{ N},...,x_{n}^{N})\end{pmatrix},\end{split}  \tag{10}\]  
\((x_{1}^{1},x_{2}^{1},...,x_{n}^{1}),...,(x_{1}^{r},x_{2}^{r},...,x_{n}^{r})\) are the points near the solid boundary, \((x_{1}^{1},x_{2}^{1},...,x_{n}^{1}),...,\)\((x_{1}^{N},x_{2}^{N},...,x_{n}^{N})\) are the knot points of the bilinear spline (3), \(B(\overline{\alpha}_{i}(t_{j}),t_{j})\) is the velocity vector at the time step \(t_{j}\):  
\[B(\overline{\alpha}_{i}(t_{j}),t_{j})= \begin{pmatrix}V_{n}(x_{1}^{1},x_{2}^{1 },...,x_{n}^{1},t_{j})\delta(f(x_{1}^{1},x_{2}^{1},...,x_{n}^{1},\overline{ \alpha}_{i}(t_{j})))\\ &\vdots\\ V_{n}(x_{1}^{N},x_{2}^{N},...,x_{n}^{N},t_{j})\delta(f(x_{1}^{N},x_{2}^{N},...,x_{n}^{N},\overline{\alpha}_{i}(t_{j})))\end{pmatrix} . \tag{11}\]  
For the compliance minimization problem (1), the normal velocity \(V_{n}\) can be derived as follows:  
\[V_{n}=E_{ijkl}\epsilon_{ij}(u)\epsilon_{kl}(u)-\lambda \tag{12}\]  
where \(\lambda\) is the Lagrange multiplier to deal with the constraint of the volume fraction, which is calculated using the following augmented Lagrangian updating scheme:  
\[\lambda^{k+1}= \begin{cases}\mu G^{k}&k\leq n_{R}\\ \lambda^{k}+\gamma^{k}G&k>n_{R}\end{cases} , \tag{13}\]where \(\mu\) and \(\gamma^{k}\) are parameters in the \(k\)-th iteration of the optimization, and \(\gamma^{k}\) is updated using the following scheme:  
\[\gamma^{k+1}=min(\gamma^{k}+\Delta\gamma,\gamma_{max})\quad k>n_{R}. \tag{14}\]  
where \(\Delta\gamma\) is the increment and \(\gamma_{max}\) is the upper limit of the parameter \(\gamma\). Since the volume fraction of initial design usually does not meet the prescribed volume fraction, the volume constraint is relaxed in the first \(n_{R}\) iterations as below:  
\[G^{k}=\int_{D}H(f)d\Omega-\left(V_{0}-\left(V_{0}-V_{max}\right)\frac{k}{n_{R} }\right)\quad k\leq n_{R}. \tag{15}\]  
To avoid the unbounded growth of the free-form function, an approximate \(\delta(f)\) function is inserted into the velocity vector of (11), defined as follows:  
\[\delta(f)= \begin{cases}0&f<-\Delta\\ \frac{3}{4}(1-\frac{f^{2}}{\Delta^{2}})&-\Delta\leq f\leq\Delta\\ 0&f>\Delta\end{cases} , \tag{16}\]  
where the parameter \(\Delta\) is used to control the magnitude of the upper and lower bounds of the free-form function. The ersatz material model is used in the finite element analysis to obtain the naturally extended velocity field from the strain energy density in the whole design domain. The MATLAB implementation of the proposed optimization algorithm is given in the Appendix A.