On a topology optimization problem governed by two-dimensional Helmholtz equation

The paper deals with a class of shape/topology optimization problems governed by the Helmholtz equation in 2D. To guarantee the existence of minimizers, the relaxation is necessary. Two numerical methods for solving such problems are proposed and theoretically justified: a direct discretization of the relaxed formulation and a level set parametrization of shapes by means of radial basis functions. Numerical experiments are given.


Introduction
Helmholtz equation is used in modelling e.g. polarized electromagnetic waves in two dimensions or acoustic waves in two and three dimensions. Applications involving shape and topology optimization governed by the Helmholtz or Maxwell equation include among others: stealth technology [15], noise reduction [1], enhancing performance of acoustic devices (horn-lens combinations) [24], anti-reflection surfaces [9], and optimal distribution of conductivity minimizing the dissipative energy [11]. Further, in combination with the homogenization approach it can be used to improve of metamaterial properties in different wave propagation problems (e.g. electromagnetic waveguides, band gap structures using photonic chrystals) [7,12,14,17,19,21], and reconstruction of electromagnetic inclusions by boundary measurements [4].
In this paper we restrict ourselves to topology optimization problems governed by the scalar time-harmonic wave equation, i.e. the Helmholtz equation. Theoretical studies of the topic are rare. Usually they are proposed by the inverse problems community. For example Mark Kac asked in article [13] "Can one hear the shape of a drum?". Recently, Bao et al. [3], have proved some useful results applicable also in studies of certain topology optimization problems.
We formulate a model topology optimization problem that is motivated by the above mentioned papers dealing with electromagnetic waveguides and band gap structures. We assume that there exists an obstacle in the free space that has different (but constant) relative permittivity compared to that of the free space. When the incoming wave travels between media of different relative permittivity it may both reflect and refract off a material interface. Consequently, in some areas the wave field is amplified or damped. This phenomenon can be utilized to prevent the incoming wave from entering some regions or propagating in certain directions. Our model problem, however, does not directly correspond to any concrete engineering design problem as we want to keep the formulation compatible with the one presented in [3].
As we have mentioned, a numerical realization will be based on an optimal control approach. The state problem will be given by the two-dimensional Helmholtz equation in which the relative permittivity ε r plays the role of the control variable and has the following form: ε r = ε r,0 + (ε r,1 − ε r,0 )χ Ω , where ε r,0 , ε r,1 is the constant relative permittivity of the free space, and the scatterer, respectively, and χ Ω stands for the characteristic function of a measurable set Ω which represents the scatterer. The choice of a cost functional J depends on particular optimization goals. Our aim is to find the distribution of ε r,0 and ε r,1 in such a way that J attains its minimum. It is well-know that this class of optimization problems has no solution, in general, i.e. no minimizer of J over the set the admissible ε r introduced above exists. For this reason an extension of this set is necessary. Since ε r appears in the lower order term of the Helmholtz equation, the respective extension is represented by all measurable functions ranging between ε r,0 and ε r,1 . A discretization of this new relaxed formulation can be directly used for solving the problem. Besides this, we propose another approach in which the parametrization of the shape of the scatter is done by a level set method involving radial basis functions. The advantage of this approach is the fact that shape/topology optimization can be transformed into a parametric optimization problem.
The paper is organized as follows. In Sect. 2 we introduce the physical setting of the state problem. In Sect. 3 we present a class of shape optimization problems governed by the Helmholtz equation and derive their relaxed form. Section 4 is devoted to a discretization and convergence analysis of both approaches mentioned above. In Sect. 5 implementation aspects are discussed. Finally, numerical results of two model examples are shown in Sect. 6.
Throughout of the paper we use the following notation: if Q ⊂ R s is a domain then H k (Q), k ≥ 0 integer, denotes the space of all functions which are together with their generalized derivative up to order k square integrable in Q (L 2 (Q) = H 0 (Q)). The norm, seminorm in H k (Q) will be denoted by k,Q , and | | k,Q , respectively. L ∞ (Q) stands for the space of all bounded measurable functions in Q. Finally, f , f , stands for the real, and imaginary part of a complex valued function f , respectively andf for its conjugate.

The state problem
Let us consider electromagnetic wave propagation in a dielectric, non-magnetic threedimensional infinite medium where no charges and currents exit. Classical Maxwell's equations can be reduced to the vector wave equation for the electric field E = (E 1 , E 2 , E 3 ) as follows: where ε r is the relative permittivity of the material and ε 0 and μ 0 are the permittivity and permeability of the free space, respectively. If we assume the time-harmonic transverse electric (TE) mode, the x 3 -component of the electric field is of the form where u tot is the complex amplitude of the wave and ω is the angular frequency. From the above assumptions it follows that (1) reduces to the scalar Helmholtz equation where the wave number Let ε r,0 be the relative permittivity of the inifinite medium. A scatterer having relative permittivity ε r,1 > ε r,0 is represented by a possibly multiply connected set Ω. Next we assume that the scatterer is contained in a domain Ω.
Let us introduce a computational domain Fig. 1. The relative permittivity can be now represented in where χ Ω is the characteristic function of Ω. Without loss of generality, we may assume that ε r,0 = 1 so that ε r = 1 + q 1 χ Ω =: 1 + q with q 1 = ε r,1 − ε r,0 > 0 and q ∈ L ∞ (Ω). For purely technical reasons and to simplify notations we will assume that q 1 = 1 in Sects. 3 and 4. An incident plane wave u 0 = e ik 0 d·x propagating into the direction d ∈ R 2 , d = 1 is scattered from the obstacle. The total electric field u tot can be split into the incident field u 0 and the reflected field u: The scattered field u should not be reflected back from the artificial boundary ∂ D. Therefore, as usual, we impose an absorbing boundary condition as an approximation to the (physically correct) Sommerfeld radiation condition In the theoretical analysis of the problem we employ the first-order approximation of (7) as the boundary condition on the artificial boundary ∂ D: To give a weak formulation of (6) and (8), we introduce the bilinear form a : Here The weak formulation of the state problem reads:

Setting of the optimization problem
The two dimensional Helmholtz equation introduced in the previous section will be used as the state problem and the relative permeability as the control variable in a class of optimization problems.
Recall the classial formulation of the state problem: Find a scattered field u : D → C solving the following boundary value problem: where χ Ω is the characteristic function of a measurable set Ω ⊂ Ω. The system of all such sets will be denoted by O and it will be identified with The weak formulation of (10) reads as follows: given q ∈ U , (P(q)) Let J : H 1 (D) → R be a functional and J : L ∞ (D) → R be defined by where u(q) is a solution to (P(q)), q ∈ U and α ≥ 0, γ > 0 are given positive parameters. J will play the role of the objective function in the following problem: Remark 1 If α > 0, then the last term in (12) can be interpreted as a penalty functional related to the constraint meas Ω ≤ γ .
As we shall see, problem (P) has no solution, in general. For this reason we shall introduce its relaxed form. The set U will be extended to the larger set U # , where and the state problem (P(q)) will be defined also for q ∈ U # . The existence and uniqueness of a solution to (P(q), q ∈ U # has been established in [3]. In addition, the estimate (14) holds, where the constant C > 0 depends only on k 0 and D. The relaxed form of (P) is given by Our aim will be to show that (P # ) has a solution, i.e. a minimizer of J on U # exists and, in addition, min q∈U # J (q) = inf q∈U J (q). To this end we shall need the following continuity assumption on J : We start with an auxiliary stability result.
Proof The sequence {u n } is bounded in H 1 (D) as follows from (14). Thus, one can find a subsequence {u n j } such that u n j j→∞ u in H 1 (D) (16) and, hence Passing to the limit with j → ∞ in (P(q n j )) and using (16), (17) we see that u solves (P(q)). Since (P(q)) has a unique solution, (16) and (17) hold for the whole sequence {u n }. To prove strong convergence we proceed as follows: making use of (16), (17) and the definition of (P(q n )), (P(q)).

Theorem 1 Problem (P # ) has a solution.
Proof Let {q n }, q n ∈ U # be a minimizing sequence: There exists a subsequence {q n j } and q * ∈ U # such that and by Proposition 1. From (15), (19), and (20) we see that l = J (q * ).
If q * ∈ U # is a solution of (P # ) and u * ∈ H 1 (D) the respective solution of (P(q * )) then (q * , u * ) will be termed an optimal pair of (P # ).
Now we introduce an auxiliary problem with a piecewise constant approximation With any T κ we associate the set We easily obtain the following existence result.

Theorem 2 Problem (P # κ ) has a solution for any κ > 0.
Proof is parallel to the one of Theorem 1.
Next we establish a relation between solutions to (P # ) and (P # κ ) when κ → 0+.
Proof The existence of a subsequence satisfying (21) follows from the definition of U # and Proposition 1. To verify that (q * , u * ) is an optimal pair of (P # ) we use the fact that the system {U # κ }, κ → 0+ is dense in U # in the weak * topology: for anȳ q ∈ U # given there exists a sequence {q κ },q κ ∈ U # κ such that and by Proposition 1. From the definition of (P # κ j ), where κ j is the filter of indices from (21) it follows that Letting j → ∞ and using (15), (21), (22), and (23) we get Any q κ ∈ U # κ , κ > 0 fixed, being piecewise constant on T κ can be attained as the weak * limit of locally periodic, rapidly oscillating functions. We briefly describe the construction of such functions.
Let Y =]0, 1[×]0, 1[ be the unit periodic cell and denote by O Y the system of all measurable sets ω ⊂ Y . The characteristic function χ ω of ω ∈ O Y will be periodically extended from Y on R 2 . For any ε > 0 we define the scaled characteristic function χ ε ω (x) := χ ω (x/ε), x ∈ R 2 and the set Let q κ ∈ U # κ be given and denote q i κ := q κ |Gi , i = 1, . . . , n(κ). It is well known (see [6] Consequently, Proposition 2 For any κ > 0 it holds: Proof Clearly min Let q * κ be a solution to (P # κ ). Using the approach described above, one can find a sequence {q ε κ },q ε κ ∈ U ε κ such that and This and (15) entails the equality in (25). Now we are ready to prove the main result of this section.

Theorem 4 It holds: min
Proof From (24) and Proposition 2 we have On the other hand Since the opposite inequality is automatically satisfied we proved (26).

Discretization of (P # ). Convergence analysis
This section deals with an approximation of (P # ). We present two ways: (i) a direct discretization of (P # ) (ii) discretization based on the level set approach parametrized by radial basis functions.

Direct discretization of (P # )
In Sect. 3 we introduced auxiliary problem (P # κ ) which is already a partial discretization of (P # ). The admissible set U # in (P # ) was replaced by U # κ which consists of piecewise constant functions on T κ . In what follows we shall study the full discretization of (P # ).
Let {T h }, h → 0+ be a regular system of triangulations ofD which satisfies the standard assumptions on the mutual position of triangles T ∈ T h . Although T h , T κ are two different partitions ofD and Ω, respectively, they will not be entirely independent. Next we shall suppose that each G i ∈ T κ is polygonal and T h |Gi defines a standard triangulation ofḠ i , i = 1, . . . , n(κ). Moreover we shall suppose that the norm κ of T κ is a decreasing function of h and κ(h) → 0+ ⇐⇒ h → 0+. With any T h we associate the space of piecewise affine complex valued functions V h : Let κ, h > 0 be fixed. The fully discrete problem (P # κh ) reads as follows: is a solution of the Galerkin approximation of (P(q κ )): Using classical compactness arguments, one can easily prove that (P # κh ) has a solution for any h, κ > 0.
In what follows, we shall study the relation between solutions to (P # ) and (P # κh ) when h → 0+. To this end, we shall need the following stability result on the solutions Comments on the satisfaction of (27) will be done later on. We start with the following auxiliary result which is parallel of Proposition 1.

Proposition 3 Let (27) be satisfied and {q
and and u := u(q) is a solution to (P(q)).
Proof In view of (27), one can find a subsequence {u h j } such that Letṽ ∈ H 1 (D) be arbitrary but fixed.
The definition of (P h j (q κ j )), where κ j = κ(h j ) with {h j } as in (29) yields: Letting j → ∞, using (28), (29), (30) and the fact thatṽ ∈ H 1 (D) is arbitrary we see that u solves (P(q)) and (29) holds for the whole sequence. Strong convergence of {u h } to u can be established exactly as in Proposition 1.
The next theorem is parallel of Theorem 3.
Moreover (q * , u * ) is an optimal pair of (P # ). Any accumulation point of {(q * κ , u * h )} in the sense of (31) has this property.
Proof The existence of a subsequence {(q * κ j , u * h j )} satisfying (31) is obvious. Let q ∈ U # be an arbitrary element. Then there exists a sequence {q κ },q κ ∈ U # κ such thatq κ κ→0+q weakly * in L ∞ (D) and at the same time by Proposition 3. From the definition of (P # κ j h j ) it follows: where h j , κ j is a filter of indices for which (31) holds. Passing to the limit with j → ∞ and using the continuity assumption (15) we obtain

Remark 2
In what follows, we justify the validity of (27). Let f ∈ L 2 (Q) and G ∈ H 1 (Q) be given. Denote g := G |∂ Q and suppose that f, g satisfy the balance condition. It has been shown in [10] that every solution of the Neumann problem belongs to H 2 (Q) provided that Q is a convex polygonal domain and, in addition, there exists a constant C > 0 which does not depend on f , G and such that From this and (6), one can easily show that the unique solution u(q) to (P(q)) belongs to H 2 (D) for any q ∈ U # and there exists a constant C > 0 which depends only on k 0 and D such that Since the bilinear form a(·, ·) appearing in (P(q)) satisfies the Gårding inequality, there exist: h 0 > 0 and a constantC > 0 which depend only on k 0 and D such that (P h (q κ )) has a unique solution for any h ≤ h 0 and any [18,20]). From this and (32) we obtain (27).

Discretization by the level set method parametrized by radial basis functions
In this section we describe and theoretically justify the partial discretization of (P # ) using a level set approximation of U # but still keeping the continuous setting of the state problem (P(q)), q ∈ U # . Level-set techniques have been gained attention during the past decade. See [23] for a review on the topic. We have chosen radial basis function (RBF) approach which was already used by the authors in the paper [22].
First we introduce the following notation. For any square Q we denote by H (Q) the set of all positive rationals δ defined as follows: δ ∈ H (Q) if and only if there exists a partition T δ (Q) of Q into squares of size δ whose interiors are mutually disjoint. In the rest of the paper we shall assume that Ω is a square.

Hence,
It is readily seen that the equality holds in (33). Indeed, if q * ∈ U # κ is a solution to (P # κ ) and {Δ n } is a sequence such that Δ n 0 then q * Δ n := q * −Δ n → n→∞ q * in L ∞ (D), Since q * Δ n ∈ U # κΔ n , the opposite inequality holds in (33). Consequently making use of (15) and (24). Next we shall approximate functions from U # κΔ , κ ∈ H ( Ω), Δ ∈]0, 1[ by sequences of scaled characteristic functions of a special class of sets which are periodically distributed in each G i ∈ T κ ( Ω).
Here and in what follows we shall suppose that κ ∈ H ( Ω) and ε ∈ H (G i ) ∀G i ∈ T κ ( Ω). For any m ∈ N, m ≥ 1 we construct a regular grid of points {C kl } m k,l=1 , C kl = ( k m+1 , l m+1 ), lying in the unit periodic cell Y =]0, 1[×]0, 1[. Denote by B kl,r the disc with the center at C kl and of radius r and set ω m,r = m k,l=1 B kl,r . It is easy to see that for any Δ ∈]0, 1[ one can findm ∈ N,r > 0 which depends solely on Δ and such that Clearly, for any ρ ∈ [0, 1 − Δ] there exists r ∈ [0,r ] such that meas ωm ,r = ρ and ωm ,r ⊂ Y . To emphasize the role of ωm ,r , we denote by Y r , εY r the periodic cell containing ωm ,r , and its scaled form, respectively (see Fig. 2). With any G i ∈ T κ ( Ω) we associate the parameter r i which may take any value between 0 andr . The scaled cell εY r i will be "stamped" onto all G i j ∈ T κ (G i ) (see Fig. 3). The characteristic function of all inclusions in G i generated by εωm ,r i ⊂ εY r i (black part of G i ) will be denoted by χ ε m,r i in what follows. For any κ, ε > 0 and Δ ∈]0, 1[ we introduce the set (37) We prove the following result: Proposition 4 Let κ > 0 and Δ ∈]0, 1[ be given. Then for any q ∈ U # κΔ there exists a sequence {q ε }, q ε ∈ U ε κΔ such that q ε ε→0+ q weakly * in L ∞ (D).
Problem (P ε κΔ ), κ, ε > 0 and Δ ∈]0, 1[ can be easily formulated as a parameter identification problem. Indeed, let Φ : [0, ∞[→ R be a function which is positive in where {C kl } m k,l=1 is the regular grid of points in Y constructed above. Let Δ ∈]0, 1[ be given andm ∈ N andr > 0 be such that (35) and (36) are satisfied. Since supp Ψ kl,r = B kl,r then ωm ,r = supp Ψm ,r , where Ψm ,r : Y → R, If H denotes the Heaviside function then H (Ψm ,r ) is the characteristic function of ωm ,r ⊂ Y , r ∈ [0,r ]. Next, Ψm ,r will be periodically extented from Y to the whole R 2 (keeping the same notation) and then scaled: Ψ ε m,r (x) = Ψm ,r ( x ε ), x ∈ R 2 , ε > 0. It is easy to see that after an appropriate translation of the grid ε Z , where Z = (k, l), k, l integer, the restriction H (Ψ ε m,r ) |Gi is the scaled characteristic function χ ε m,r from (37).
Let κ, ε > 0 and Δ ∈]0, 1[ be fixed. Denote by U ad the following subset of R n(κ) : Clearly, any function q ∈ U ε κΔ can be identified with a vector R ∈ U ad and optimization problem (P ε κΔ ) can be equivalently written in the following form: Remark 4 Instead of (43), one can use a more general definition of Ψm ,r , namely where r ∈ [0,r ], α kl ∈ [α min , α max ] ∀k, l, α min < α max given and α ∈ R m 2 is the vector of α kl , k, l = 1, . . . ,m. The appropriate choice of the coefficients α kl gives us more flexibility in building supp Ψm ,r,α ⊂ Y . Indeed, if α kl = 0 for some k, l then the respective supp Ψ k,l,r contributing to ωm ,r is missing. If for example Δ ∈]1 − π 4 , 1[ then instead of m 2 small discs of radius r at C kl the union of which satisfies (36) one can use only one disc of radius r ∈ [0, 1 2 ] and center ( 1 2 , 1 2 ) ∈ Y to fulfill (36). Analogously (44), let Q be a function of α and R: Fig. 4 Graph of the radial basis function Ψ i j,r i j defined by (47) where R ∈ U ad and α = (α 1 , . . . , α n(κ) ) ∈ (R m 2 ) n(κ) . In this case, (P ε κΔ ) can be formulated as a minimization problem for the control variables α and R.

Numerical realization
In this section we describe numerical realization of the radial basis function approach. As in computations one is not able to use a high number of radial basis functions, we use just one regular grid of points {C i j } uniformly distributed in Ω (see Fig. 5) and define the level set function as the linear (affine) combination of RBF basis functions as follows: where suppΨ i j,r i j is the disc with center at C i j and of radius r i j . Thus, the design variables of the parametrized optimization problem are represented by the vector β := (α, R) ∈ R N , N := 2m 2 . The constant shift Ψ 0 is optional and its role will be established later. The function Φ in (42) used for constructing Ψ i j,r i j is chosen as ( [25]): The respective radial basis function Ψ i j,r i j is depicted in Fig. 4. The function q in (11) can be now written with the aid of the level set function Ψ β := Ψ α,R and the Heaviside step function as follows: The weak formulation of the state problem (10) is discretized by picewise linear finite elements as discussed in Sect. 4.1. By q β we denote any of the finite dimen-sional parametrizations of q introduced in the previous section. The finite element approximation of (10) then reads as follows: where Remark 5 In oder to obtain physically realistic results, the simplest first-order boundary condition (8) on ∂ D is inappropriate unless ∂ D is very far from ∂ Ω. Therefore, the second-order approximation of (7) proposed by Engquist and Majda [8] and Bamberger et al. [2] was implemented. Instead of a(q β ; z, w) we use its modified form where s is the unit tangent vector on ∂ D and C P denotes the set of four corners of ∂ D. The finite element approximation of (10) used in computations reads as follows: A structured and uniform triangulation T h of D for constructing V h is shown in Fig. 5.
Let h > 0 and m ∈ N be fixed. The approximate solution u h ∈ V h is represented as a linear combination u h = i c i ϕ i , where c i ∈ C and {ϕ i } are the standard real-valued Courant basis functions. The matrix form of (50) reads: where the matrix A(β) and the right-hand side b(β) can be decomposed as follows: The entries of the above matrices and of the vector f (β) are given by Matrices K , C, D, and G are evaluated exactly. Their sum is a pentadiagonal sparse matrix. To compute the design dependent "mass" matrix M(β) and the excitation vector f (β), we utilize numerical integration. Let V (T ) be the set of nodal numbers of vertices associated with the triangle T . Moreover, let x T , x l denote the center of gravity and the l:th nodal point of T , respectively. The following integration formulae are used: Note that the former yields a diagonal approximation of M(β). Finally, we arrive at the following fully discrete nonlinear programming problem: where is the set of admissible design parameters, where α min , α max and r i j,max ∀i, j are given. The discretized cost function J in (P) is not continuously differentiable due to the discontinuous Heaviside function H . If one wishes to use descent-type methods, then smoothing of H is necessary. Therefore, H will be replaced by a smooth function H s defined by where s > 0 is given. As H s (0) = 1 2 , it is necessary to add a negative shift Ψ 0 to the level set function Ψ β to ensure that H s (Ψ β (x)) = 0 for x / ∈ Ω. Alternatively, one can omit the shift and modify the definition of H s .
To compute the gradient of J , the standard adjoint variable approach will be used: Here p solves the adjoint problem where stands for the conjugate transpose of a matrix and the complex gradient ∇ u is defined by The partial derivatives of M(β) and f (β) in (53) are easy to compute. From (48) we have:  1). In the relaxed formulation, the coefficient q ∈ U # is represented by where G i j is the system of m 2 small squares which define the partition of Ω and q 1 = ε r,1 − ε r,0 (see (4)).
The level set function used in computations by the radial basis approach has the following form: where Ψ i j,r i j is defined by (42). Unlike (46) where the coefficients of the linear combination α i j are independent of r i j we set α i j = r i j in (55). The motivation for this choice of α i j is the following: the maximum of Ψ i j,r i j is always equal to one for any r i j > 0. Therefore, when r i j → 0 the basis function Ψ i j,r i j becomes an increasingly sharp peak resulting in numerical difficulties. The maximal radius r i j,max of supp Ψ i j,r i j in the definition of A ad is 1.4/(m+1) for all C i j which are closest to ∂ Ω and 2.8/(m +1) for the rest. The smoothing parameter s in (52) is chosen by using the following heuristic rule: s = s 0 h, where s 0 = 4. We do not employ more sophisticated adaptive techniques to determine s used e.g. in [22] as they would increase the nonlinearity of the optimization problem and the heuristic choice works well enough. The shift in (55) is chosen as The state and adjoint solver as well as the cost function evaluation were implemented using MATLAB [16]. Minimization was carried out by fmincon with 'sqp' option from the MATLAB Optimization Toolbox. All computations were done on a desktop computer.
Example 1 This example deals with a simplified problem of designing the so-called bandgap structures. The aim is to design the scatterer Ω in such a way that the wave having the wave number k 0 is prevented to enter a specified subdomain D 0 ⊂ D. We use the simplest choice of the cost function, namely i.e. the absolute value of the total wave u tot (see (5)) is minimized in D 0 . The physical parameters are k 0 = 6π and q 1 = 0.75. We choose D 0 = ]−0.6, 0.6[×]0.7, 1[, and the mesh size h = 2/300. The problem was again solved by both the relaxed and the radial basis function formulations. We used grid of 20 × 20 design variables in both cases.
The initial guess for the relaxed formulation was α i j = 0.5∀i, j. After 121 optimization iterations (122 function evaluations) the cost function was reduced to the value 3.40 × 10 −4 . The optimized design and the corresponding total wave are shown in Fig. 6.
The optimized design has a significant amount of grey cells. In topology optimization during optimization process, usually some filtering or penalizing techniques are used to force the design to the upper or lower bound [5]. However, sophisticated filtering techniques are out of scope of this paper. Therefore, we only threshold the optimized design variables α * i j to either 0 or 1, i.e. we setα * i j = H (α * i j − 1 The initial guess for the radial basis function formulation was r i j = 0.05∀i, j. After 47 optimization iterations (50 function evaluations) the cost function decreased to the value 5.95 × 10 −5 . The optimized design and the corresponding total wave are shown in Fig. 8.
The optimized design contains many tiny discs which are obviously due to numerical noise and having no physical meaning. Therefore we repeated the optimization by taking as the initial guess the final solution of the previous run but seting those r i j smaller than 0.08 to zero. After 38 optimization iterations (44 function evaluations) the cost function was reduced to 6.07 × 10 −5 . The results are shown in Fig. 9. The cost is only slightly higher but the obtained design is now almost free from tiny discs.
Comparing the optimized designs obtained by the relaxed and radial basis function formulations, we see that they differ significantly on details. However, both designs essentially meet the design target, i.e., they prevent wave entering into D 0 . It is obvious that the optimization problem has many global optima and/or has many local minima with slightly higher cost function value than the global minimum. Finally, we recomputed the example (20 × 20 radial basis function discretization, initial guess r i j = 0.05∀i, j) but with the regularized cost functionJ (u) = J (u) + α D H s (Ψ ), α = 10 −2 . After 117 optimization iterations (123 function evaluations) the cost function was reduced to 3.84 × 10 −3 . The optimized design and the corresponding total wave are shown in Fig. 10.
The optimized design is almost free of isolated tiny disks. However, it is not symmetric any more.

Example 2
In this example, the basic setting of the problem is the same as in Example 1 except the area constraint meas(Ω) = γ added as a penalty term into the cost function: We choose the target area as γ = 3 10 and the penalty parameter as α = 10. The initial guess for the relaxed formulation (30 × 30 basis functions) was α i j = 0.25∀i, j. After 250 optimization iterations (253 function evaluations) the cost function was reduced to the value 4.87 × 10 −3 . The optimized design and the corresponding total wave are shown in Fig. 11. Convergence of the cost functional with respect to optimization iterations is shown in Fig. 12. The initial guess for the radial basis function formulation (20 × 20 basis functions) was r i j = 0.05∀i, j. After 162 optimization iterations (183 function evaluations) the cost function decreased to the value 2.14 × 10 −3 . The optimized design and the corresponding total wave are shown in Fig. 13. Convergence of the cost functional with respect to optimization iterations is shown in Fig. 14.
As can be noticed, the additional area constraint also serves as a good way to avoid oscillating designs.

Conclusions
This paper deals with a class of topology optimization problems governed by the two dimensional Helmholtz equation and the permittivity as the control variable. Its main contribution consists in the rigorous convergence analysis using both, the classical relaxed formulation and a level set/radial basis function approach.
The numerical tests confirmed that this type of problems is difficult from the implementation point of view, since many designs produce very similar solutions to the state problem. Also the discretization of the control variable using radial basis functions having variable radii seems to increase the nonlinearity of the optimization problem. However, the presence of a penalty or regularization functional helped to stabilize the computations.
Computations which are based on the relaxed formulation with a penalization or filtering techniques are currently dominant in practical applications. On the other hand it is worthwhile to study alternative techniques, among them level set methods which attained much interest in the past years. The advantage of level set approaches is the fact that they work on the "safe side", i.e. the designs are automatically required to have values ε r,0 and ε r,1 except in a very small region near the material boundaries due to the smoothing of the Heaviside function. Also smooth boundaries of simple geometrical shapes are better represented by the radial basis function approach.
Variants of the level set/radial basis functions approach presented in this paper, in particular the case when a single fixed radius defines the basis functions is the object of future studies.