Topology optimization of an acoustic diode?

By using topology optimization, we consider the problem of designing a passive acoustic device that allows for one-way flow of sound waves; such a device is often colloquially referred to as an acoustic diode. The Helmholtz equation is used to model the time harmonic linear wave propagation together with a Dirichlet-to-Neumann (DtN) type boundary condition, and the finite element method is used for discretization. The objective of this study is to maximize the wave propagation in one direction (from left to right) and minimize the wave propagation in the reverse direction (from right to left) for planar incoming waves. The method of moving asymptotes (MMA) solves the optimization problem, and a continuation approach is used for the penalizing intermediate design variables. The results for the optimized waveguide show that more than 99.8% of the power of planar incoming waves get transmitted from left to right while less than 0.3% gets transmitted in the reverse direction for planar incoming waves in the specified frequency range. Since a true diode is a non-reciprocal device and here we used a linear acoustic wave model, which is basically reciprocal, we discuss details about how it appears to be possible to obtain a one-way waveguiding effect using this linear model.


Introduction
Like an electric diode, which only allows current to flow in one direction, an acoustic diode only allows acoustic waves to propagate in one direction. The magnificent impact of the invention of the electric diode on human life inspired a lot of research work trying to control other forms of energy in the same way. Since ultrasound is widely used in biomedical imaging and non-destructive diagnostics, recently, the concept of acoustic diodes has raised interest among researchers in physics and engineering. One way to get directional acoustic wave propagation in a waveguide, is to use nonlinear materials to break the reciprocity in the transmission (Liang et al. 2009(Liang et al. , 2010Popa and Cummer 2014;Grinberg et al. 2018;Darabi et al. 2019;Li et al. 2019Li et al. , 2020. This method, however, suffers from low efficiency due to considerable losses in the nonlinear conversion (Chen et al. 2017;Song et al. 2016). It has also other limitations, such as significant alteration and distortion of the frequency content of the transmitted acoustic signal (Grinberg et al. 2018), as well as bandwidth limitations of working frequencies (Zhu et al. 2015b). This motivates the efforts on the investigation of linear and passive acoustic one-way structures, trying to find the design of sound-hard material so that the problem setup works as an acoustic diode (He et al. 2011;Li et al. 2012a, b;Zhu et al. 2015a;Chen et al. 2017;He and Kang 2018). Here, we use material distribution-based topology optimization to optimize the design of a passive acoustic one-way waveguide.
Topology optimization is one of the most general and effective design optimization methods, where the aim is to find the best layout of material in a design domain such that the objective function is maximized (or minimized). For material distribution-based topology optimization, we divide the design domain into elements, and for each element, the optimization algorithm decides whether it will hold material or not. This method was first introduced for structural mechanics using the homogenization method (Bendsøe and Kikuchi 1988). Since then, researchers have adopted this method in many fields and applications, such as acoustics (Wadbro and Berggren 2006;Dühring et al. 2008;Lee and Kim 2009), electromagnetics (Hassan et al. 2014;Aage and Johansen 2017), optics (Cox and Dobson 1999;Jensen and Sigmund 2011;Mori et al. 2019), and fluid mechanics (Borrvall and Petersson 2002). A summary of early research works on material distribution-based topology optimization can be found in the monograph by Bendsøe and Sigmund (2003). Here, we introduce a novel optimized design of a directionally selective device that provides a one-way waveguiding effect, for the specified frequency range.
Besides the research content, this paper also has an educational message, which will be clear to the reader before the end of the article.
For now, in order not to spoil the intended message, we simply ask the reader to (just as usual) have an open and critical mindset. Moreover, the structure of the presentation in this article is non-standard; the presentation is structured as follows. Sections 2 and 3 provide a brief explanation of the problem description and also present the results for the optimized design. The purpose in these two sections is to give an understanding of the problem and the methods we used to solve it without getting into details. At the end of Section 3, we introduce and discuss conceptual issues regarding the problem setup. The remainder of the paper is spent on resolving these issues. More precisely, Section 4 provides important aspects of the resulting wave propagation and revisits the optimization problem formulation. Sections 5 and 6 contain results in more detail and a discussion on the essential issues mentioned above. Finally, details about the filtering algorithm used in the topology optimization are provided in the Appendix.

Brief problem description
Consider the axi-symmetric setup illustrated in Fig. 1. This setup comprises two semi-infinite waveguides and the gray meshed region Ω D , which we will use as our design domain. Some parts of Ω D may be filled with sound-hard material to make this setup a one-way waveguide. From here on, outgoing waves refer to all waves propagating away from the design domain Ω D and incoming waves refer to all waves propagating toward the design domain Ω D .
Consider the following two cases: Case 1 There is an incoming planar acoustic wave in the left waveguide, propagating towards Ω D . Fig. 1 a Axi-symmetric setup with cylindrical design domain Ω D in the middle and one cylindrical waveguide on both sides. b Case 1: incoming wave at the left waveguide propagates to right. c Case 2: incoming wave at the right waveguide gets reflected to the right Case 2 There is an incoming planar acoustic wave in the right waveguide propagating towards Ω D .
As illustrated in Fig. 1, we seek a design so that the incoming wave in Case 1 propagates to the other side of the setup, while the incoming wave in Case 2 gets reflected.

Mathematical model
The wave propagation inside the fluid (air) region of the setup in Fig. 2 is governed by the linear wave equation. By considering time harmonic wave propagation, we assume acoustic pressure P (x, t) = Re{p(x)e −iωt }. This gives us Helmholtz equation for the complex amplitude function p of the acoustic pressure. In cylindrical coordinates, we have: where the wave number k = ω/c, c is the speed of sound in the fluid, ω is the angular frequency, andΩ is the fluid (air) region of the setup. For the computational experiment, we consider following dimensions for the setup. The radius and length of the Fig. 2 The computational domain, Ω, consists of design domain Ω D and two truncated waveguides. The sound-hard material inside Ω D is defined as Ω s , and the rest of the setup filled with fluid (air) is defined asΩ. The boundaries consist of axi-symmetric Γ sym (dash-dotted line), the sound-hard walls Γ s (solid line), and artificial truncated boundaries Γ L and Γ R (dashed lines) design region Ω D is r D = 50 mm and l D = 50 mm, respectively. The radius and length of the waveguides are r W = 30 mm and l W = 20 mm, respectively.
To deal with the infinite waveguides numerically, we truncate them and use a DtN (Dirichlet-to-Neumann) type perfect absorbing (or non-reflecting) boundary condition on the artificial boundaries Γ L and Γ R , as illustrated in Fig. 2. For a general description of the Dirichlet-to-Neumann boundary condition, we refer the reader to the book by Ihlenburg (1998, p. 63). Hence, the state equation reads: where g L and g R are the incoming waves from the left and right waveguides, respectively. Here, we consider a planar incoming wave for both case 1 and case 2. The variational formulation for the states (2)-(5) is: Remark 1 Throughout this paper, we do not write symbols for measure such as dΩ, dΓ R , or dΓ L in integrals as long as their omission is not a source of confusion.
To model the distribution of sound-hard material in Ω D (that is, the shape of region Ω s ), we define a material indicator function, α. The function α ≡ 0 in Ω s and α ≡ 1 in Ω \ Ω s . By using this function to extend the integration domain of the first two integrals in problem (6), fromΩ to Ω, we obtain the following problem: where Ω =Ω ∪ Ω s is the extended domain that contains both the air and solid material. To allow gradient-based optimization, the standard approach in material distributionbased topology optimization is to allow α to take values in the range [0, 1].

Remark 2
The model of letting the parameter control the design act directly on the domain integral can at least be traced back to 2006 (Wadbro and Berggren 2006) and has been used in many contributions thereafter. For a pure solidfluid design, where the fluid is represented by α = 1 and the solid scatterer by α = , Kasolis et al. (2015) proved that the error in the complex-pressure field converges linearly in to the solution of the problem with an exactly modeled scatterer.

Discretization
We use the finite element method to discretize and solve problem (7). To discretize the domain Ω, we use a uniform mesh of square elements, and define a functional space V h that consists of functions that are continuous and biquadratic on each element. We denote the basis functions of V h by ϕ j , j = 1, 2, . . . , N N , where N N is the number of nodes or degrees of freedom in the finite element approximation. We approximate the complex amplitude function p and the test function q by p h ∈ V h and q h ∈ V h , respectively. In addition, we use an element-wise constant function α h to approximate material indicator function α. The discretized version of problem (7) reads: We define the N N × N N stiffness K and mass M matrices with entries and also N N × 1 vectors b L and b R with entries respectively. Furthermore, we define α = α 1 , α 2 , . . . , α N D T that holds the element values of α h , where N D is the number of elements in Ω D , and p = p 1 , p 2 , . . . , p N N T that holds the nodal values of complex pressure amplitude in Ω. By using the definitions above (8) reads: where the matrix B corresponds to the Dirichlet-to-Neumann boundary conditions (4) and (5). More details about how to compute matrix B can be found in an earlier article by Wadbro (2014, Appendix A).

Objective function and optimization problem
The transfer parameter T XY from X to Y is defined as: Power of the outgoing wave on Γ Y given an incoming planar wave with unit power at Γ X , where X ∈ {L, R} and Y ∈ {L, R}, in which L and R denote the left and right waveguides, respectively. Section 4 provides detailed information on how the transfer parameters are evaluated for a given design.
Our goal is to maximize the power transmitted to the right in case 1 and minimize the power transmitted to the left in case 2. Which is to minimize the power at Γ L for both cases, based on conservation of power, whether the incoming wave is from the left or right waveguide. So for case 1, we minimize reflection T LL , and for case 2, we minimize transmission T RL . The objective function can thus be written as: where N f is the number of frequencies subject to optimization. So, the optimization problem consists of finding the distribution of sound-hard material in Ω D that minimizes the objective function (13). Figure 3 shows the axi-symmetric optimized waveguide, and Fig. 4 is its 3D visualization. For planar incoming waves, this design allows acoustic power to be transmitted from left to right but not from right to left (to within 0.3% of the incident power) with frequencies in the range 8-9 kHz. Figure 5 shows the transfer parameters as functions of frequency for the optimized design in Fig. 3. The top graph corresponds to case 1, where there is an incoming planar wave at left waveguide, and shows the resulting transmission (T LR ) and reflection (T LL ) as functions of frequency. This graph shows that throughout the frequency band of interest, the incoming planar wave is transmitted from left to right waveguide and there is essentially no reflection at left waveguide. More precisely, for all integer frequencies in the range 8-9 kHz, T LL < 1.4 × 10 −3 , which means that at most 0.2% of the input power is reflected back and thus more than 99.8% of the power is transmitted to the right waveguide. Throughout this paper, by integer frequencies in the range 8-9 kHz, we refer to the frequencies 8000, 8001, . . . , 9000 Hz. Similarly the bottom graph in Fig. 5 shows the transmission (T RL ) and reflection (T RR ) as functions of frequency for case 2, when there is an incoming planar wave at right waveguide. In this case for all integer frequencies in the range 8-9 kHz, T RL < 2.8 × 10 −3 which means more than 99.7 % of power is reflected back to the right waveguide and at most 0.3 % of the input power is transmitted to the left waveguide.

Results
The design presented above appears to act as a perfect acoustic diode for the frequency band of interest. However, the observant and critical reader may already have thought about raising the following questions and remarks: -What about acoustic reciprocity? -The diode is a non-linear device; if so, how and why is it possible to create such a device using a linear model?
The above questions are indeed valid, and we have on purpose withheld some information hoping that the reader  would keep a critical mindset. We stress that we have not made any claims that the design in Fig. 3 is an acoustic diode; we stay with the weaker statement that it appears to act as one. We devote the remainder of this paper to describe and detail what we have done to arrive at this result.

Detailed problem description
Here, we present details about how our waveguide setup works by including the discussion of modes. In Section 4.1, we will elaborate how to compute these eigenvalues and eigenmodes for the waveguide setup. In Section 4.2, we provide the full optimization problem.

Wave propagation in a circular duct
Consider semi-infinite waveguides on the left and right sides of Ω D , as illustrated in Fig. 1. Assume that the complex pressure amplitude p satisfies Helmholtz (2) and boundary condition (3). Using separation of variable, we write the general solution in the left waveguide as: and solution in the right waveguide as: where z L and z R are the positions of z-axis on Γ L and Γ R , the functions f m (r) are the modes at left and right waveguides (in the continuous case, it is well known that the functions f m are so-called Bessel functions), e is the base of the natural logarithm, i is the imaginary unit, and A L m and B L m are complex constants that determine the amplitude of incoming and outgoing waves, respectively, at the left waveguide. Similarly, A R m and B R m are complex constants that determine the amplitude of incoming and outgoing waves, respectively, at the right waveguide, and the constants k m are the so-called reduced wave numbers.
To obtain a discrete exact radiating boundary condition, we will, instead of using the Bessel functions, use eigenfunctions to a finite element problem on the boundary mesh on Γ L and Γ R .
In short, we solve one-dimensional eigenproblems in the radial direction on Γ L and Γ R . Both these problems are of the form: (Note that since the radii of the left and right waveguide are the same, we only need to solve one problem. In the text below, Γ refers to either Γ L or Γ R .) To solve the one-dimensional eigenproblem on Γ , we use a restriction of the finite element mesh on Ω to Γ . We let V h Γ be the space of functions that are bi-quadratic on each element of Γ . We solve a discrete version of problem (16) to find the M eigenvalues λ m and the corresponding eigenfunctions or modes f h m ∈ V h Γ , where m = 0, 1, . . . , M. We assume that the eigenvalues are sorted and the lowest eigenvalue λ 0 = 0. Due to orthonormality of modes, we have: Since the complex amplitude function p satisfies Helmholtz (2) and boundary condition (3), we conclude that the reduced wavenumbers satisfy: Note that for the wave to propagate, the reduced wave number k m must be real, so λ m ≤ k 2 . That is, for a given frequency f , there is only a finite number M p f of socalled propagating modes. Waves of modes higher than M p f are geometrically evanescent, and any such outgoing mode decays exponentially in the waveguide.
The Dirichlet-to-Neumann type boundary conditions allow us to set the amplitudes A m of the propagating incoming modes and perfectly absorb all outgoing waves.
We set the amplitude A L m and A R m of the incoming wave as in case 1, and in case 2.
Having solved governing equation (11), we can compute the outgoing wave amplitudes B L m and B R m . For X ∈ {L, R}, it holds that: where the first equality follows from expression (14) or (15) and the second equality follows from relation (17). This expression holds for both Γ L and Γ R , so where the N N × 1 vectors c L m and c R m have entries

Objective function with modes
The outgoing waves can be written as a sum over propagating modes, and the power of each of these modes can be monitored, which enables us to achieve a greater control and understanding of the resulting wave propagation. Thus, for X ∈ {L, R}, Y ∈ {L, R}, and i, j ∈ {0, 1, 2, . . . , M p f }, we define the so-called scattering parameters as: Power of mode j of the outgoing wave on Γ Y given an incoming wave of mode i with unit power at Γ X .
The power of a propagating wave of mode m is proportional to the square of its amplitude and the corresponding reduced wave number k m . Provided that the input wave is of mode i at Γ X , then the corresponding scattering parameters can be evaluated as: The transfer parameter T XY defined in Section 2.3 is related to scattering parameters S X i Y j as Remark that transfer parameters are defined for the planar incoming wave. Hence, only X 0 is included in the definition above. Using the definition above, the objective function (13) can be re-written as: where we minimize reflection of power over M p f modes at left waveguide for case 1, and transmission of power over M p f modes at left waveguide for case 2. In the numerical experiments, we use a combination of filtering and a penalization method to suppress the intermediate values in N D × 1 vector d that defines the material distribution inside Ω D before filtering. To this end, we let α := F(d), where F(d) is a filtered version of d (details in the Appendix) and add a quadratic penalization term (Allaire and Kohn 1993) to the objective function. In summary, we solve the optimization problem: where γ is a penalty parameter and is a set of admissible designs in Ω D .

Results revisited
Here, we revisit the results from Section 3. In Fig. 5, the frequency response of the optimized design shows transmission and reflection of an incoming planar wave in terms of transfer parameters. In this section, we will present all the details to show transmission and reflection of incoming wave in terms of scattering parameters. That is, how the output power is distributed over the propagating modes. With a given waveguide setup, the propagation of modes depends upon the frequency. For our waveguide setup (recall Fig. 2), the first non-planar mode starts propagating at 7 kHz, and the second non-planar mode starts propagating at 13 kHz. For the numerical experiments, we consider the size of the device and select frequencies that only allow two modes to propagate in our waveguide setup. That is, only the planar mode (m = 0) and the first non-planar mode (m = 1) propagate.
To obtain the design in Fig. 3, we minimize objective function (28) for the N f = 11 frequencies 8. 0, 8.1, 8.2, . . . , 9.0 kHz. More precisely, we solve the optimization problem for 25,600 design variables in Ω D using the method of moving asymptotes (MMA) by Svanberg (1987). For penalization of design variables, we use a two-step approach. In the first step, we set γ = 1 and initial design d to the vector with all ones, and solve the optimization problem. The optimization terminates when the residual norm of first order optimality condition computed by the MMA algorithm is less than 10 −3 . In the second step, we increase the penalty parameter by setting γ = 10 2 and solve the updated optimization problem using optimal design from the first step as initial guess to obtain our final optimized design. We use the same termination condition for MMA in both steps. Figure 6 presents graphs of transmission and reflection distributed over modes with planar incoming wave at left Fig. 6 Normalized power of reflected planar (solid line) and nonplanar (solid line with diamond marker) mode, and transmitted planar (dotted line) and non-planar (dotted line with diamond marker) mode as a function of frequency (kHz) for case 1 and case 2 and right waveguide, respectively. For case 1, the graph shows that nearly all the power of the incoming planar mode at left waveguide is transmitted to the first non-planar mode at the right waveguide, throughout the frequency band of interest. More precisely, more than 99.8% of the power is transmitted to the non-planar mode at the right waveguide for all integer frequencies in the range 8-9 kHz. Due to conservation of power, there is essentially no reflection at left waveguide, as values of S L 0 L 0 < 5.4 × 10 −4 and S L 0 L 1 < 9.0 × 10 −4 , respectively, for all integer frequencies in the range 8-9 kHz. Similarly, there is approximately zero transmission at the planar mode at the right waveguide, as value of S L 0 R 0 < 2.2 × 10 −3 for all integer frequencies in the range 8-9 kHz.

Planar mode as incoming wave
For case 2, the graph shows that nearly all the power of incoming planar mode at right waveguide is reflected back throughout the frequency band of interest. More precisely, more than 99.5 % of power is reflected back to the planar mode at the right waveguide for all integer frequencies in the range 8-9 kHz. Moreover, there is essentially no transmission to the left waveguide, as values of S R 0 L 0 < 2.2 × 10 −3 and S R 0 L 1 < 4.2 × 10 −4 , respectively, for all integer frequencies in the range 8-9 kHz. Similarly, there is approximately zero reflection at the first non-planar mode of the right waveguide, as value of S R 0 R 1 < 2.3 × 10 −3 for all integer frequencies in the range 8-9 kHz.
For case 1, the wave propagation is essentially in one direction, that is, from left to right waveguide throughout the frequency band of interest. Similarly, for case 2, there is essentially no wave propagation in the left waveguide, and incoming planar mode is reflected back to the right waveguide throughout the frequency band of interest. By using mode conversion for both case 1 and case 2, the numerical experiment successfully minimizes the objective function (28), and the waveguide setup appears to act as an acoustic diode for an incoming planar wave at the left and right waveguides, respectively.

First non-planar mode as incoming wave
To test if the optimized design in Fig. 3 acts as an acoustic diode for the first non-planar mode as incoming wave, we define two more cases: Case 3 There is an incoming wave of the first non-planar mode in the left waveguide, propagating towards Ω D . Case 4 There is an incoming wave of the first non-planar mode in the right waveguide, propagating towards Ω D . Figure 7 presents graphs of transmission and reflection distributed over modes with the first non-planar mode as incoming wave at the left and right waveguides, respectively. For case 3, the graph shows that nearly all the power of incoming first non-planar mode is reflected back throughout the frequency band of interest. More precisely, more than 99.7% of the power is reflected back to the left waveguide for all integer frequencies in the range 8-9 kHz. Moreover, there is essentially no transmission towards the right waveguide, as values of S L 1 R 0 < 4.2 × 10 −4 and S L 1 R 1 < 9.9 × 10 −4 , respectively, for all integer frequencies in the range 8-9 kHz. Similarly, there is approximately zero reflection at planar mode of the left waveguide, as value of S L 1 L 0 < 9.0 × 10 −4 for all integer frequencies in the range 8-9 kHz.
For case 4, nearly all the power is transmitted to planar mode at left waveguide with an incoming first non-planar mode of the right waveguide throughout the frequency band of interest. More precisely, more than 99.6% of the power is transmitted to the left waveguide for all integer frequencies in the range 8-9 kHz. There is essentially no reflection at the right waveguide, as values of S R 1 R 0 < 2.3 × 10 −3 and Fig. 7 Normalized power of reflected planar (solid line) and nonplanar (solid line with diamond marker) modes, and transmitted planar (dotted line) and non-planar (dotted line with diamond marker) modes as a function of frequency (kHz) for case 3 and case 4 S R 1 R 1 < 5.4 × 10 −4 , respectively, for all integer frequencies in the range 8-9 kHz. Similarly, there is approximately zero transmission towards the first non-planar mode of the left waveguide, as value of S R 1 L 1 < 9.9 × 10 −4 for all integer frequencies in the range 8-9 kHz.
For both case 3 and case 4, the device does not appear to act as an acoustic diode (from left to right) and there is wave propagation from right to left. Similarly, like case 1 and case 2, here the device uses mode conversion to achieve one-way waveguiding effect in reverse (from right to left). Thus, the one-way waveguiding effect depends on the input mode. Fig. 8 Plots of Re(p) at 8 KHz for case 1, case 2, case 3, and case 4. We provide videos of plots of Re(p) for all four cases as supplementary material Fig. 9 Four-port network: In case 1, the incoming planar mode (L 0 ) input is converted to the first non-planar mode (R 1 ). Similarly in case 4, the incoming first non-planar mode (R 1 ) is converted to planar mode (L 0 ). On the other hand, incoming planar (R 0 ) and non-planar modes (L 1 ) are reflected back in case 2 and case 3, respectively

Concluding discussion
The results of the numerical experiments presented in Fig. 6 show that the waveguide setup appears to act as an acoustic diode allowing the incoming planar wave to propagate from left to right. However, Figs. 7 and 8 show that for the first non-planar mode as incoming wave, it appears to act as an acoustic diode in reverse, allowing propagation from right to left. Hence, the device designed in this paper is in fact a high-quality mode-converter or a four-port network, illustrated in Fig. 9, instead of a diode.
Recall that the diode is a non-linear device that allows one-way transmission of any incoming wave. This is not possible for a linear device due to reciprocity. The reciprocal property of this setup is due to Helmholtz equation, where the time-harmonic approximation of acoustic pressure for the incoming P (x, t) = Re e −iωt p(x) and the outgoing wave P (x, t) = Re e iωt p(x) are the time-reverse of each other. Thus, S X i Y j = S Y j X i . For example, S L 0 R 1 (case 1) is equal to S R 1 L 0 (case 4) where the optimized design inside the design domain converts planar mode to non-planar mode and vice versa. It is in fact not possible to design a true acoustic diode in reciprocal systems.
Generally, all so-called linear acoustic diodes in the literature are just mode converters which work as multi-port networks, connecting different modes to each other. Despite all possible applications for a linear directional waveguide, diode is just the wrong term for it, as a diode is a non-linear device and to achieve that, we need to break reciprocity, which is an innate property of linear wave propagation. That being said, it is indeed possible to design a highly efficient linear directional waveguide. The results show that by using topology optimization and selecting frequency band that allow two or more propagating modes, we can design a multi-port network that works as an efficient directional waveguide for an incoming planar wave.
Thus, an executive summary of the above reads: -The optimized design in this paper is a high-quality mode converter, but not an acoustic diode. -Diodes are non-linear devices-there is no such thing as a diode in a reciprocal system. -The so-called linear diodes presented in the literature are in fact mode converters. (They should thus be treated and addressed as such and not as diodes.)

Filtering
In this work, we use non-linear density filters that are based on morphological operators (Sigmund 2007). To use gradient-based optimization algorithm, these morphological operators are approximated by harmonic means-based filters (Svanberg and Svärd 2013). For optimization, we use the design domain Ω D ; however, for filtering, we extend the domain on three sides, illustrated in Fig. 10. The extended domain for filtering is Ω F = Ω D ∪ Ω W ∪ Ω M , where Ω W represents the left and right waveguides, and Ω M is the domain with solid material. The material distribution is fixed in Ω W and Ω M , that is, α h (x) = 1 (air) in Ω W and α h (x) = 0 (solid) in Ω M . The discretization in Ω F consists of N F square elements of the same size as those in Ω. The elements in Ω F are arranged such that elements of Ω D come first, Ω M second, and Ω W last.
Furthermore, we define neighborhood N r,i of radius r around element i with centroid x i , where N r,i consists of j elements with centroid x j . We assemble a matrix G r with Fig. 10 Filtering domain Ω F = Ω D ∪ Ω W ∪ Ω M binary entries, where g ij = 1 if j ∈ N r,i , else g ij = 0. If the distance between element i and j is less than r, that is, x i − x j ≤ r then element j belongs to N r,i . We define a weight matrix W r = D −1 r G r , where D r = G r 1 N and 1 N is a N × 1 matrix with all entries equal to 1.
We define functions f (x) E r = 1/(x + β) and f (x) D r = f E r (1 − x) to compute discrete harmonic mean, where β > 0, and f −1 E r and f −1 D r are inverse of these functions. Using notation of fW-mean filters (Wadbro and Hägg 2015), we define discrete erode and dilate operators of the N F × 1 vector ρ with entries in the range 0, 1 by where . . , f (ρ N f ) D β are element-wise functions.
Using the definition of harmonic erode and dilate in expression (30), we define harmonic close that acts on N D × 1 vector d F(d) where I N d is an identity matrix, 0 N d ×(N W +N M ) is a zero matrix. Similarly, 0 N M is a zero vector and 1 N W is a vector with all entries equal to 1. The F(d) gives us filtered version of design vector d that holds element values of α h in Ω D .
Funding Open Access funding provided by Umea University.