Acceleration of a Multi-objective Topology Optimisation in 2D Electro-Magnetic Field Based on the Level-Set Method and the Boundary Element Method by the $\mathscr{H}$-Matrix Method

Topology optimisation is the most flexible configuration optimisation method based on computer simulation and has extensively been researched so far. In topology optimisations, configuration of design objects is expressed with a function defined in a design domain and the function is updated in a way to minimise an objective function, which expresses the performance of design objects. Many varieties of methods for the shape representation are proposed, such as a homogenisation-based method [1] and solid isotropic material with penalization (SIMP) [2]. In this study, we employ the level-set method, which can express the boundaries of design objects clearly [3].

We update the configuration of design objects iteratively based on a topological derivative, which is the variation of the objective function when an infinitesimal disc, for 2D, is allocated in the design domain. To evaluate the topological derivative, we need to calculate the electro-magnetic responses for two kinds of incident waves; forward and adjoint. The finite element method (FEM) is usually employed to solve the forward and adjoint problems. In the FEM, however, the infinite domain should be approximated by a finite domain with an appropriate boundary condition. Hence, the FEM needs a large analysis domain in order to assure enough accuracy, which may cause large meshing costs. We apply the BEM to calculate the magnetic responses since the BEM requires mesh generation only on the boundary and can treat the infinite domain strictly.

The level-set-based topology optimisation with the BEM has been applied to many varieties of problems such as acoustics [4], heat conduction [5] and electromagnetics [6]. In many cases, algebraic equations that appear in a sensitivity analysis have the same coeffi-cient matrix. Hence, we can solve these algebraic equations in one time with direct solvers such as the LU decomposition. The BEM, however, leads fully populated coefficient matrix after discretisation, which causes large numerical costs to solve the algebraic equation with direct solvers. To realise an efficient direct solver, we apply the $\mathscr{H}$-matrix [7] method to accelerate the LU decomposition. The H-matrix method is an acceleration method for matrix operations based on hierarchical blocking of the matrix and low rank approximation. We can easily apply this method to direct solvers. Calculation cost for the LU decomposition can be reduced to O (N log N) by the use of the H-matrix method. Furthermore, in some multi-objective topology optimisation problems, the required number of electro-magnetic field analyses increases in proportion to the number of the objective functions. Hence, for this kind of problems, further acceleration is expected by the use of the $\mathscr{H}$-matrix method and the LU decomposition. We demonstrate the effectiveness of our method by an example of multi-objective optimization problems for 2D electro-magnetic field.

In this section, we formulate the electro-magnetic field analysis in a transverse electric (TE) polarised field with the BEM. We consider a vacuum domain in which dielectric materials are allocated (Fig. 1). We denote the vacuum domain and dielectric materials as $\Omega_1$ and $\Omega_2$, respectively. Also, we define $\Gamma$ as the boundary of $\Omega_2$. The magnetic field $u$ satisfies the following boundary value problem in a TE polarised field:

where $k_i$ and $\varepsilon_i$ denote the wave number and permittivity in the domain $\Omega_i$, respectively. Also, $\partial / \partial n$ denotes the inward normal derivative on the boundary $\Gamma$. Superscripted variables $u^i$ and $(\partial u=\partial n)^i$ in boundary conditions (3) and (4) express the limit from the domain $\Omega_i$ to the boundary $\Gamma$. In this study, we employ the PMCHWT formulation [8] to avoid the irregular eigenfrequency problem and obtain the following boundary integral equations:

where $u^{\text {inc }}$ denotes an incident wave. Also, $w$ and $w^{\text {inc }}$ are defined as follows:

$S i, D i, D_i^*, N i$ denote, respectively, the following operators:

By discretising the boundary integral equations (6), we obtain the algebraic equation to solve for the unknown quantities.

Topology optimisation is a method to determine the optimal structure of engineering devices to maximise its performance by iteratively updating the configuration in a design domain D based on computer simulations. In order to formulate a topology optimisation problem, we need to define an objective function, which expresses the effectiveness of a device. In this study, we define the objective functions in terms of the functions of the magnetic responses at some observation points $x_i^{\text {obs }}$, as follows:

In the following subsections, we present a level-set-based topology optimisation, which uses a level-set function to express the configuration of the design objects, and the topological derivative to update the configuration.

In this study, we express the configuration of design objects with the level-set method. In the level-set method, we define the level-set function $\phi(x)(-1 \leq \phi(x) \leq 1)$ in a design domain $D$. Each domain is expressed based on the value of the level-set function as follows:

We update the level-set function by the following equation:

where $t$ is a fictitious time corresponding to an optimisation step. The first term in the right hand side (RHS) denotes the direction and scale of the update. $C$ is a constant which is used to adjust the scale. Further, $\mathscr{F}(x)$ is the topological derivative, a variation of the objective function caused by allocating an infinitesimal disc in a design domain $D$. The second term of RHS is used to adjust the complexity of the obtained configuration. $\tau$ is a positive constant and $l$ denotes a characteristic length of the design domain $D$. The design domain $D$ is divided into a structured grid. We update $\phi$ at each grid point by solving the equation (17). This problem does not require regeneration of mesh. Further, the analysis domain is finite. We, therefore, employ the FEM to solve equation (17).

Topological derivative is defined as the first coefficient of the asymptotic expansion of the objective function by area of an infinitesimal disc s(ε):

where $\delta J$ is a variation of the objective function. In a TE polarised field, the topological derivative when an infinitesimal circle of dielectric element is allocated in the vacuum domain $\Omega_1$ is derived as follows:

where $\tilde{u}(x)$ is the solution of the following adjoint problem:

When a dielectric element $\Omega_2$ is pierced with an infinitesimal circle of vacuum, the topological derivative is derived by swapping $\varepsilon_1$ and $\varepsilon_2$ of (19).

In the case of multi-objective optimisation problems, we employ the KS function [9] to define the objective function as follows:

in which each $J_i$ is a function of magnetic response and $w_i$ denotes its corresponding weight. The topological derivative for the objective function (25) is calculated by the chain rule as follows:

where $\mathscr{F}_i(x)$ denotes the topological derivative when we define $J_i$ as the objective function.

The topology optimisation algorithm is illustrated in Fig. 2. Firstly, we determine an initial configuration of a design object and initialise the level-set function value in accordance with the initial configuration of the material. The boundary mesh is generated from the level-set function. Then, we conduct the forward analysis with the BEM. As a result of the forward analysis, we obtain the value of the objective function and the electro-magnetic response in the design domain. When the objective function does not satisfy a convergence criteria, we conduct the adjoint analysis. By using the forward and adjoint electro-magnetic responses, we calculate the topological derivative. Then, we update the level-set function by the equation (17).We repeat the above process until the objective function satisfies the convergence criteria.

In this study, we accelerate the calculation of the topological derivative by the $\mathscr{H}$ - matrix method. $\mathscr{H}$-matrix method is an acceleration method of matrix operations based on hierarchical blocking of a matrix and low rank approximation. With the help of the $\mathscr{H}$-matrix method, we can reduce calculation costs for matrix operations. By the use of such efficient operations, we can reduce the calculation cost for LU decomposition to $O(N \log N)$ for the number of boundary element $N$. We call the accelerated LU decomposition with the $\mathscr{H}$-matrix method as HLU. We can solve the forward problem consisting of equations (1)-(5) and the adjoint problem consisting of equations (20) - (24) in a fast manner simultaneously by the HLU. In this section, we show the method for blocking a coefficient matrix and low rank approximation.

Each row and column of a coefficient matrix which is derived by the BEM corresponds to each boundary element. Hence, we divide the coefficient matrix into some block matrices based on clustering of the boundary. We generate a boundary cluster by the following procedure:

•  We make a rectangle enclosing a boundary, which we are going to divide. We define a set of nodes in the rectangle as a ‘cluster’. Also, we denote the number of divisions that is required to obtain the cluster as ‘level’.

• We divide the longer side of the rectangle into two parts and define each obtained rectangle as a new cluster.

• We repeat the above process until the number of nodes in each cluster is less than a preset parameter $n_{\min }$. We define a cluster in which the number of nodes is less than $n_{\min }$ as a ‘leaf cluster’.

By the above process, we obtain a binary tree structure of the boundary clusters. Figure 3 shows an example of boundary clustering. ‘Level’ denotes a depth of each cluster in the binary tree structure and corresponds to the required number of clustering to obtain each cluster. Combination of two clusters in the same level denotes a block matrix of the coefficient matrix. When the two clusters satisfy the following admissibility condition, we define the block matrix as a far-field effect block:

where $\eta$ is a constant that arranges strictness of the admissibility condition. diam $C_i$ denotes the maximum distance between two nodes included in the cluster $C_i$. Also, $\operatorname{dist}\left\{C_i, C_j\right\}$ denotes the minimum distance between two nodes included in each $C_i$ and $C_j$ (Fig. 4). If a block matrix $C_s \times C_t$ is not a far-field effect block and at least one of $C_s$ and $C_t$ is a leaf cluster, we define the block matrix as a near-field effect block.

We approximate far-field effect blocks by low rank matrices with the adaptive cross approximation (ACA) [7]. The ACA is a method to approximate a matrix by the algebraic methods and can reduce the calculation cost and required memory for the coefficient matrix.

In this section, we apply the efficient sensitivity analysis with the H matrix method to a multi-objective topology optimisation problem. We show that the proposed method is advantageous to a conventional method.

We consider a problem to find the optimal configuration of dielectric element in a design domain $D=[0,90] \otimes[0,90]$, which minimises the following objective function for multiangle incidence:

in which $J_i$ is a sum of electro-magnetic fields for an incident wave of angle $\theta_i$ and is expressed as follows:

We consider the case where the angles of incidence are $0^{\circ}, 3^{\circ}, 6^{\circ}, 9^{\circ}$ and $12^{\circ}$. The wave length $\lambda$ is 40 . We allocate 60 observation points in a domain $[100,120] \otimes[35,55]$. The permittivity of the dielectric element $\varepsilon_2=5$. As an initial configuration, we put four circles of dielectric element whose radii are all 9.75 (Fig. 5). The weight function of the KS function $w i$ is fixed to 1.0 in this example. As is shown in Section 3, we have to calculate the topological derivatives $\mathscr{F}_i$ for $J_i$ to calculate the topological derivative $\mathscr{F}$ for the objective function $J$. It causes us to solve algebraic equations for $n$ incident waves, which have the same coefficient matrix. We solve these equations with $\mathscr{H}$-matrix method and the LU decomposition (HLU) in one time. In this study, parameters for $\mathscr{H}$-matrix method $n_{\text {min }}, \eta$ and tolerance for addition of two $\mathscr{H}$-matrices and the ACA are fixed to $32,1.4$ and $10^{-6}$, respectively.

Figures 6 and 7 show histories of the objective functions and the optimal configuration in the case where HLU and GMRES is employed as a solver. Sum of the magnetic fields of each angle of incidence is reduced to the same degree. Also, the history of the objective function is almost the same as each other. This means that the approximation error of ACA is relatively small compared to the discretisation error. Calculation time of each step is shown in Fig. 8, in which the horizontal axis shows the number of boundary elements for each optimisation step. By using HLU, the total calculation time from step 1 to 80 is reduced approximately by 75%.

We developed an efficient sensitivity analysis approach for topology optimisations using the $\mathscr{H}$-matrix method and the LU decomposition. We confirmed that computation cost spent to obtain the optimal configuration for a multi-objective optimisation problem is reduced by approximately $75 \%$ compared to the conventional method with GMRES. We are going to extend the proposed method to apply for optimisation problems in which we have to treat more complicated structures such as design problems of cloaking devices. For such problems, the structure of the $\mathscr{H}$-matrix becomes further more complicated, increasing the numerical cost for matrix operations.