## 2 Software Structure and Methods
### Topology Optimization  
The developed system uses the modified topology optimization algorithm presented in [20]. For the 2D topology optimization case the compliance minimization problem used in the system can be expressed as follows:  
\[ \begin{array}{l}\text{minimize(over $f$) }J_{0}(u)=\int_{D}E_{ijkl}\epsilon_{ ij}(u)\epsilon_{kl}(u)H(f)d\Omega,\\ \text{subject to: }a(u,v,f)=\left\{v,f\right\},\quad\forall v\in U,\\ \left.u\right|_{\Gamma_{1}}=u_{0},\\ \left.E_{ijkl}\epsilon_{ij}(u)\right|_{\Gamma_{2}}=p,\\ \int_{D}H(f)d\Omega\leq V_{0},\end{array} \]  
where  
\[ \begin{array}{l}a(u,v,f)=\int_{D}E_{ijkl}\epsilon_{ij}(u)\epsilon_{kl}(v)H(f)d \Omega,\\ \left.l(v,f)=\int_{D}gv+\
abla(p\
u n)H(f)d\Omega,\\ H(f)=\begin{cases}1&\text{if }f\geq 0,\\ 0&f<0,\end{cases}\end{array}  \tag{1}\]  
\(J_{0}(u)\) is an elastic energy, the objective function; \(u\)--the displacement field; \(\epsilon_{ij}(u)\) and \(\epsilon_{kl}(u)\)--the linearized strain tensors; \(E_{ijkl}\)--the Hook elasticity tensor; \(u_{0}\), \(p\) and \(g\)--the constant values that represent the given displacement, traction and body force respectively; \(a(u,v,f)\) is the linear elastic equilibrium equation written in its weak, variational form and \(l(v,f)\) is the load linear form, with \(v\) denoting a virtual displacement field in the space \(U\) of kinematically admissible displacement fields; \(V_{0}\) is the maximum allowable volume in the design domain \(D\); \(H(f)\) is the Heaviside function that represents the presence of material at a point of the level-set function value \(f\):  
\[H(f)= \begin{cases}1&\text{if }f\geq 0\\ 0&\text{if }f<0\end{cases}  \tag{2}\]  
This problem statement uses the level-set function \(f\) which, at the same time, is an FRep defining function \(f(x_{1},x_{2},...,x_{n})\) (see Section 2.2):  
Figure 1: The principal components of the FRep system. \[f(x_{1},x_{2},...,x_{n}) \begin{cases}>0,&\text{if }(x_{1},x_{2},...,x_{n})\in \Omega\\ =0,&\text{if }(x_{1},x_{2},...,x_{n})\in\partial\Omega\\ <0,&\text{if }(x_{1},x_{2},...,x_{n})\
ot\in\overline{\Omega}\end{cases}  \tag{3}\]  
where \(f\) is a real continuous function defined on \(E^{n}\), \(n\) equals 2 for 2D, and 3 for 3D topology optimization, respectively. Let us consider the 2D optimization. The proposed topology optimization algorithm incorporates a free-form defining function with bilinear basis functions:  
\[\phi(x,y)= \begin{cases}(1-|x|)(1-|y|)&\text{if }|x|\leq 1\text{ and }|y|\leq 1\\ 0&\text{otherwise.}\end{cases}  \tag{4}\]  
This function defines a surface drawn in Figure 2a. Moreover, in general case, we need to use this equation with constants \(h_{1}\) and \(h_{2}\) to define a rectangular non-zero region of the spline basis function, as follows:  
\[\phi^{\prime}(x,y)= \begin{cases}(1-\frac{1}{h_{1}}|x|)(1-\frac{1}{h_{2}}|y|)& \text{if }|x|\leq h_{1}\text{ and }|y|\leq h_{2}\\ 0&\text{otherwise}\end{cases}  \tag{5}\]  
However, the form (4) can be used here without the loss of generality. This form of the spline basis function is well suited for the further discussion. Space mapping is used to locate spline basis function at some portion of the domain. The bilinear spline basis function with the center of the non-zero region at \((x_{i},y_{i})\) is:  
\[\phi_{i}(x,y)= \begin{cases}(1-|x-x_{i}|)(1-|y-y_{i}|)&\text{if }|x-x_{i}|\leq 1\\ &\text{and }|y-y_{i}|\leq 1\\ 0&\text{otherwise.}\end{cases}  \tag{6}\]  
Therefore, the free-form defining function (3) used here can be expressed as follows:  
Figure 2: Bilinear spline construction: (**a**) the surface defined by the bilinear spline basis function (4); (**b**) the surface produced by two spline basis functions and its zero-level set (boundary); (**c**) rectangular domain D with 90 \(\times\) 25 elements and 91 \(\times\) 26 knots. \[f(x,y)=\sum_{i=1}^{N}\alpha_{i}\phi_{i}(x,y) \tag{7}\]  
with constants \(\alpha_{i}\). The main steps of the proposed modified algorithm are shown in Algorithm 1. ```
Step 1.