On minimum length scale control in density based topology optimization

The archetypical topology optimization problem concerns designing the layout of material within a given region of space so that some performance measure is extremized. To improve manufacturability and reduce manufacturing costs, restrictions on the possible layouts may be imposed. Among such restrictions, constraining the minimum length scales of different regions of the design has a significant place. Within the density filter based topology optimization framework the most commonly used definition is that a region has a minimum length scale not less than D if any point within that region lies within a sphere with diameterD > 0 that is completely contained in the region. In this paper, we propose a variant of this minimum length scale definition for subsets of a convex (possibly bounded) domain. We show that sets with positive minimum length scale are characterized as being morphologically open. As a corollary, we find that sets where both the interior and the exterior have positive minimum length scales are characterized as being simultaneously morphologically open and (essentially) morphologically closed. For binary designs in the discretized setting, the latter translates to that the opening of the design should equal the closing of the design. To demonstrate the capability of the developed theory, we devise a method that heuristically promotes designs that are binary and have positive minimum length scales (possibly measured in different norms) on both phases for minimum compliance problems. The obtained designs are almost binary and possess minimum length scales on both phases.

are several reasons to restrict the set of feasible designssome of mathematical nature related to the well-posedness of the topology optimization problem, and some of more practical nature related to manufacturing. Constraints on the minimum feature size can be used to ensure manufacturability and help controlling the manufacturing cost. Moreover, as pointed out by Allaire et al. (2016), imposing a minimum feature size constraint can compensate for simplified physical modeling or be introduced for aesthetical reasons.
In density based topology optimization, M ⊂ is represented by its characteristic function 1 M : → {0, 1}-often referred to as the material indicator function. To utilize gradient based optimization algorithms, which are suitable for handling the large-scale problem obtained after discretization, the range of the material indicator function is relaxed to [0,1]. The relaxed material indicator function ρ : → [0, 1] is referred to as the density. By introducing a penalty on densities in (0, 1), binary designs are promoted. To guarantee existence of solutions, the density based topology optimization problem needs in general to be suitably regularized. Borrvall (2001) presents a systematic overview of several regularization schemes. Regularization by density filtering, where the density is restricted to the range of a regularizing filter operator, is frequently used (see, for instance, the works by Bruns and Tortorelli 2001;Bourdin 2001;Guest et al. 2004;Sigmund 2007;Lazarov and Sigmund 2011;Svanberg and Svärd 2013). We have previously introduced the class of generalized f W -mean filters that provide a common framework for the vast majority of density filters used in topology optimization (Wadbro and Hägg 2015). In a recent paper we extended Bourdin's (2001) result on existence of solutions to the linearly filtered minimum compliance problem to cascades of generalized f W -mean filters (Hägg and Wadbro 2017).

Overview on minimum length scale control in density based topology optimization
This section presents an overview of different methods that aim at imposing minimum length scale in topology optimization, with emphasis on the various definitions of minimum length scales involved. A comprehensive review on different techniques based on density filtering has recently been given by Lazarov et al. (2016).
Early attempts to gain control over the minimum length scales in density based topology optimization, such as slope constraints (Petersson and Sigmund 1998) and linear density filtering (Bruns and Tortorelli 2001;Bourdin 2001), focused on limiting the spatial variation of the density field. The linear filter replaces the density in one point with a weighted linear average over some neighborhood. The typical choice has been to use d-dimensional spherical neighborhoods with a given diameter D > 0. Necessarily, both methods (the addition of slope constraints and linear density filtering) lead to blurred transitions between material and void regions. If the slope constraint reads ∇ρ(x) ≤ D −1 or if the filter neighborhood is a d-dimensional sphere with diameter D > 0, then the transitional region is at least of size D. As indicated before (Sigmund 2007) and as will be emphasized in this paper, other neighborhood shapes can be employed. However to fix ideas, we will stick to d-dimensional spherical neighborhoods with diameter D in this section. Concentrating the weights used in the linear filter around the center of each neighborhood, allows for sharper transitions between material and void regions. Ultimately, when the weights far from the center of the neighborhood become negligible in comparison to those close to the center, we are effectively using a smaller neighborhood-In finite precision arithmetic, this is certainly the case. Nonlinear density filters (Guest et al. 2004;Sigmund 2007;Svanberg and Svärd 2013) have the advantage of allowing for almost sharp transitions between material and void regions, even for uniform weighting.
The underlying idea of the Heaviside projection, is to threshold a linearly filtered density at a fixed density level η ∈ [0, 1] (Guest et al. 2004;Sigmund 2007). When thresholding at level η = 0, the resulting density in a point is 0 if and only if the density within the neighborhood of that point is identically zero. On the other hand ρ(x) > 0 implies that the resulting density will be 1 at any point that has x in its neighborhood. Thus, any point within a material region lies within a d-dimensional sphere with diameter D > 0 that is completely contained in the material region. The text in italics may be adopted as the defining property of a material region possessing a minimum length scale not less than D. To distinguish this overarching definition we will denote it by NEighborhood based minimum Length scale (NEL). Although, we have not been able to determine the exact origins of this definition, it has been extensively used-both implicitly and explicitly-in the literature (see, for instance, the works by Bruns and Tortorelli (2001), Bourdin (2001), Guest et al. (2004), Sigmund (2007), Guest (2009), and Svanberg and Svärd (2013)). A more elaborate definition of the NEL has been given by Zhang et al. (2014). Given the structural skeleton S(M) of M = M, consisting of all points in M that have at least two nearest points on ∂M, Zhang et al. (2014) define the minimum length scale as whereB r (x) denotes the closed (Euclidean) ball centered at x with radius r. For binary densities, thresholding filtered densities at η = 0 gives the same result as computing maxima over neighborhoods (Sigmund 2007). Thus, for binary densities, thresholding filtered densities at η = 0 provides a realization of the dilate operator from mathematical morphology. Other realizations for binary densities of the dilate operator are given by the geometric and harmonic dilate filters with the parameter α = 0 (Svanberg and Svärd 2013). Analogously, considering thresholding at η = 1 leads to a definition of minimum length scale for void regions, considerations of the erode operator and its various realizations for binary densities.
In applications using gradient based optimization algorithms, thresholding is replaced by application of a differentiable function that approximates the Heaviside step function (Guest et al. 2004). Similarly, differentiable filters that approximate the dilate and erode operators are used (Sigmund 2007;Svanberg and Svärd 2013). In all cases, the degree of approximation is controlled by a scalar parameter. As pointed out by Sigmund (2007), it is important to keep in mind that although such filters give the same result on binary densities in the limit where the approximation is exact, they will in general produce different results for non-binary densities and finite approximation errors. The overall conclusion is that a lower bound on the NEL of either material or void regions can be imposed by using standard filtering methods. We also note that a common feature is that such filters approximate the morphological dilate or erode operators on binary densities. (This also applies to filters approximating morphological open or close operators since these are compositions of dilate and erode operators.) By introducing two design fields, one controlling the layout of material and one controlling the layout of void, Guest (2009) was able to impose independent minimum length scales on both phases using Heaviside projection. Strictly speaking, the definition of the NEL is only applicable to binary designs. However, it seems that the common understanding is that this is less of an issue when the designs are almost binary with almost sharp transitions between regions of different materials.
From a manufacturing point of view in two spatial dimensions, designs with NEL not less than D > 0 on the material phase could be manufactured using a circular deposition tool with diameter D, while designs with NEL not less than D on the void phase could be manufactured using a circular machining tool with diameter D (Sigmund 2007). The same reasoning applies to NELs defined using non-circular neighborhoods; for instance, if the void phase has NEL defined using some neighborhood, one may use a punch tool in the shape of the neighborhood in manufacturing.
Related to constraints on minimum length scale are so called robust formulations aiming at designs that are robust toward geometrical uncertainties (Lazarov et al. 2016). Robust formulations typically require multiple solves of the state equation at each design iteration and are therefore in general computationally demanding. Wang et al. (2011) noted that if thresholding the filtered densities at any level η ∈ [η d , η e ] yields designs with the same topology, then any design obtained by thresholding the filtered densities at an intermediate level η i ∈ (η e , η d ) possesses NELs on both phases. Later, Zhou et al. (2015) introduced geometrical conditions aiming to guarantee that thresholding the filtered densities within a range of levels leads to designs with the same topology. However, the functioñ ρ(x, y) = 0.5 + 0.5 0.3 cos(4πx) + tanh(20(y − 0.5)) 0.3 + tanh(10) , defined on the unit square provides a counter example.
Since ∇ρ = 0 in [0, 1] 2 , their conditions are vacuously satisfied. Nevertheless, thresholdingρ at η d = 0.2, η i = 0.5, and η e = 0.8, respectively, yields three different topologies, as illustrated in Fig. 1. We note that function (2) can be excluded by requiring ∂ρ ∂n = 0 at the design  (2) at η d = 0.2 (left),η i = 0.5 (middle), and η e = 0.8 (right), respectively, yields three different topologies domain boundary. Whether, appending ∂ρ ∂n = 0 to Zhou et al.'s (2015) conditions make them sufficient in general is an open question. It is interesting to note that the response of designs with NELs on both phases need not be robust against uniform dilations and erosions of the design (Zhou et al. 2015).
To simultaneously control the minimum length scale of both materials, Poulsen (2003) introduced the MOLEmethod (MOnotonicity based minimum LEngth scale). For simplicity, we state the definition of the MOLE in two spatial dimensions. The MOLE is at least D > 0 if at any point x ∈ the restriction of the density to any of the line segments x + t (cos θ, sin θ), t ∈ [−D/2, D/2], θ ∈ {0, π/4, π/2, 3π/4} is monotonic in t. Imposing that the MOLE is at least D > 0 implies that the boundary between material and void regions consists of piecewise straight lines meeting at angles no less than 3π/4, and in particular that the smallest feature is an octagon with diameter (1 + √ 2)D. To enforce a lower bound on the MOLE, one extra constraint was added to the formulation of the topology optimization problem. Allaire et al. (2016) work with a level set representation of the design and make the following definition. The thickness of M is larger than D > 0 if for each point x ∈ ∂M and any t ∈ [0, D], it holds that x − tn(x) ∈ M, where n(x) denotes the outward unit normal to M at point x.
To impose minimum thickness via constraints or penalties in the problem formulation, the definition is reformulated using the signed distance function. For a point in the interior of M the value of the signed distance function d M is the negative distance from the point to the boundary, for a point on the boundary of M the value is 0, and for a point in the exterior of M the value is the distance from the point to the boundary. In the level set context the NEL can be conveniently expressed by using the signed distance function and the skeleton of M (Guo et al. 2014b;Xia and Shi 2015): d M (x) ≤ −D/2 for all x ∈ S(M). As pointed out by Allaire et al. (2016), the behavior of their minimum thickness and the NEL can be rather different. Consider, for instance, the unit square with rounded corners with rounding radius 0 < r 1, then the minimum thickness is 1 while NEL = 2r 1.
Within the Moving Morphable Components (MMC) approach (Guo et al. 2014a) with trapezoidal components, a different notion of minimum length scale considering the individual sizes of the components and the sizes of their intersection regions was introduced by Zhang et al. (2016).

Paper outline
The remainder of the paper is organized into two main parts: Foundations of minimum size control, and Computational approach and experiments. The first part starts with Section 4.1 that gives an introduction to standard mathematical morphology for subsets of R d . To handle subsets of a bounded convex domain, we introduce a modification of standard morphology in Section 4.2. In Section 4.3, we propose definition (31) of the NEL for open subsets of a bounded convex domain and elaborate on its relation to the morphological operations of Section 4.2. Section 4.4 discusses morphological operators on Cartesian grids, their connection to the discretized density based topology optimization problem, and concludes with condition (53) that can be used to impose independent NELs on both phases. To quantitatively asses the NELs of an optimized design, we introduce some quality measures in Section 4.5. The first part concludes with a discussion in Section 4.6. The second part starts with Section 5.1 that introduces a method for compliance problems that heuristically imposes NELs on both phases via condition (53). Section 5.2 describes in detail the morphology-mimicking filters used in the numerical experiments, which are presented in Section 5.3. Section 5.4 provides a discussion on the proposed numerical method and the obtained results. A concluding summary related to both parts of the paper is given in Section 6.

Preliminaries on mathematical morphology
Mathematical morphology is a branch of image analysis originating from the works of Georges Matheron and Jean Serra in the early sixties. The main idea of mathematical morphology is to gain information about a set M ⊂ R d by probing it by a smaller set B ⊂ R d . In mathematical morphology B is referred to as a structuring element. In image analysis (topology optimization) the set M represents the binary image (design) given by the characteristic function 1 M (x). Dilation and erosion are perhaps the simplest morphological operations and they serve as building blocks for more complex operations.
Following Heijmans (1995), we define the dilation of M by B as (3) that is, D B (M) equals the Minkowski sum M ⊕ B. An immediate consequence of definition (3) is that dilation is commutative, that is, The erosion of M by B is defined as By taking the complement of definition (5) and using definition (3), we arrive at the following duality relation between The opening of M by B and the closing of M by B are defined by respectively. The open operator is anti-extensive while the close operator is extensive, that is, Moreover, these operators are both idempotent, that is, We note that for all x and all M we have that Moreover, D B (1 M (x)) = 1 holds if and only if there exists , which shows that definitions (3) and (13) are

Morphological operations on bounded domains
Note that the morphological operations presented in the previous section are defined for subsets of R d . Here, our main interest is however to consider morphological operations on subsets of a bounded convex domain ⊂ R d . We will also restrict our attention to structural elements containing 0; that is, 0 ∈ B. To this end, we introduce the operations where M ⊂ . (Note that by taking = R d the morphological operators of the previous section are retrieved, that is, and so on.) One may verify that operations (16) to (19) satisfy similar relations as the operations on R d presented in Section 4.1. For example, erosion (17) satisfies the analog of relation (6): Morphological operations (16) and (17) are dual in the sense that Open operator (18) is anti-extensive while close operator Moreover, these operators are idempotent The corresponding definitions for grayscale images for x ∈ .

Minimum length scale of sets
Let the structuring element B 0 be open, bounded, symmetric, and convex. Then the Minkowski functional We define the local length scale relative B at x ∈ M of a nonempty open set M ⊂ as the radius of the largest ball with center in M containing x such that its intersection with lies within M. That is, where B r (y) With definition (29) in mind, we define is by assumption open. Furthermore, the following lemma holds; a proof is given in Appendix A.
Having defined a local length scale, we define the NEL of M = ∅ relative B as the smallest local length scale of M, that is, For M = ∅, it is natural to define R B (∅) = 0, rather than using definition (31) and the convention that inf ∅ = ∞ that would lead to R B (∅) = ∞. We note that the case M = ∅, is of little practical importance in topology optimization. For bounded , the main difference between NEL (1) and NEL (31) lies in the treatment of regions that are close to the boundary ∂ . By letting = R d in definition (31), we essentially obtain NEL (1) in as much as the two definitions agree for sufficiently regular sets. For future reference, we define rB-regular sets are simply called r-regular (Serra 1982, § 5;Pavlidis 1982, § 7). It can be shown that an r-regular set must be of class C 1 (Duarte and Torres 2014).
We proceed by three theorems that link B-open sets and NEL (31). (16) and (18), as well as expression (4), we find that This shows that any x ∈ M belongs to (y + rB) We note that Theorem 1 makes it possible to quantitatively asses the NEL of a set using morphological operators. Moreover, Theorem 1 in combination with the idempotence (24) of the open operator provide a way to impose a lower bound on the NEL of the designs in topology optimization. More precisely, the NEL is guaranteed to be at least r by requiring that the design is the morphological opening by rB of some set. In particular, this is true if M = O rB (M).
Proof Since the open operator is anti-extensive (23) (20) we conclude that y ∈ E rB (M). Hence, there exists y ∈ E rB and b ∈ B such that m = y + rb, so by definitions (16) and (18) Theorem 2 shows that a set with positive NEL relative B is rB-open for any r > 0 less than the NEL. We also conclude that a set whose interior and exterior both possess positive NELs must have a regular boundary.
Finally, we combine Theorem 1 and Theorem 2 to show that NEL (31) could equivalently be defined using open operator (18).
Theorem 3 For M = ∅, using the convention that sup ∅ = 0, it holds that On the contrary, assume thatR B (M) > 0. Then for each > 0 there exists a strictly positive number r that satisfies Expression (34) shows that NEL (31) of a nonempty open set can in principle be estimated to any desired accuracy by computing the morphological opening for different values of r and applying a bisection method.
Allowing for the possibility that the interior is rB-open while the exterior is r B-open with r = r , defines a subclass of min{r, r }B-regular sets. That is, where V = \ M denotes the exterior of M relative . Using duality (22) of the morphological operators, condition (35) takes the equivalent form According to Theorem 1, the interior and exterior of a set satisfying (36) may possess different minimum length scales.
We could generalize further to allow for (note the prime on B in the close operator) Remark that, when using condition (37), the NELs of the interior and exterior of the set are measured in different norms. In topology optimization, the idea of imposing minimum length scales in different norms was suggested by Sigmund (2007). However, this idea has received little, if any, attention in the literature.

Minimum length scale in topology optimization using mathematical morphology
With the aim of introducing conditions that guarantee NELs on either or both phases of materials in density based topology optimization, we exploit the results of the previous section regarding minimum length scale (31) and its connection to morphological operations on bounded domains. We only discuss the discretized case, where the domain has been divided into n elements. To this end, assume that the design domain is a hyperrectangle discretized using a regular grid. Moreover, we assume that a piecewise constant function ρ h is used to approximate the design and let ρ ∈ [0, 1] n ⊂ R n be a vector that holds the values of the design at each element. For a generic neighborhood shape N ⊂ R d , we define the neighborhood N i of element i as all elements j with centroids x j in x i +N , that is, The neighborhoods are said to be symmetric if it holds that i ∈ N j if and only if j ∈ N i . For symmetric neighborhoods, the discrete counterparts of the grayscale morphological operators (26) and (27) are given by respectively. Analogously to the previous section, we assume that the neighborhood shape is the open ball with radius r > 0 in the norm · B introduced in the previous section; that is, N = rB. In this case which implies that the neighborhoods are symmetric and that i ∈ N i for all i. In this case the discrete morphological operations (39) satisfy where ρ ≤ η if and only if ρ i ≤ η i for all i. We note that operations (39) are increasing; that is, The symmetry of the neighborhoods implies that for any Hence, the discrete open operator is anti-extensive; that is, for any i Analogous arguments show that the discrete close operator is extensive By applying E to expression (44) and using that the discrete erode is increasing (42), we find that E(D(E(ρ))) ≤ E(ρ).

Moreover, extensivity (45) of C applied to E(ρ) shows that E(D(E(ρ))) ≥ E(ρ). Hence, E(D(E(ρ))) = E(ρ), which implies that the discrete open operator is idempotent
Analogous arguments show that the discrete close operator is idempotent We note that idempotence (47) of the discrete close operator also follows from the idempotence of the discrete open operator by applying expression (46) on 1 − ρ, and using the duality where 1 = (1, . . . , 1) T ∈ R n . As already indicated, the concept of minimum length scale is naturally defined for sets; that is, for binary designs. In the discrete case, any binary design is given by some ρ that satisfies the relation We note that, if ρ satisfies (49), then so does M(ρ), where M denotes any of the discrete morphological operators. Condition (37), which expresses that both the interior and exterior of the design possess NELs (possibly measured in different norms), translates to where we display the (possibly different) neighborhood shapes N and N as subscripts on the open and close operator, respectively. One way of guaranteeing fulfillment of conditions (50) and (51) in the optimization process would be to introduce two binary auxiliary design vectors η 1 , η 2 , imposing the constraint O N (η 1 ) = C N (η 2 ), and define the (physical) design as ρ = O N (η 1 ) (= C N (η 2 )).
The idempotences (46) and (47) of the discrete open and close operators imply that such ρ satisfies conditions (50) and (51). However, as we now demonstrate, introducing one extra design variable per element is superfluous as far as there is no gain in design freedom. By anti-extensivity (44) of the discrete open operator and extensivity (45) of the discrete close operator we find that holds for any design ρ. Expression (52) implies that conditions (50) and (51)

Quality measures
Based on the theoretical development presented, we propose four quality measures-two that measure the difference between the opening and the closing of the design; and two that measure how close the design is to being morphologically open and morphologically closed, respectively. The first of these is the measure of difference between open and close Thus, M DOC quantifies any residual in condition (53). A related quality measure is that is, the fraction of elements in the vector C N (ρ) − O N (ρ) that are greater than a cut off value of 0.5. In the 0-1 case, F DOC measures the number of elements that are 1 after the close operation but 0 after the open operation. We remark that M DOC = 0 implies that F DOC = 0 but the reverse is not necessarily true. To be specific, a purely 0-1 design satisfying M DOC = 0 has NELs on both phases. As will be seen in the numerical results, the attained NELs may be different for the two phases and not tight on the imposed lower bounds. To quantify such deviations, we introduce two more measures. The measure of difference between identity and open quantifies any residual in condition (50), while the measure of difference between identity and close quantifies any residual in condition (51). To be specific, a binary design satisfying M DIO = 0 (M DIC = 0) has NEL on the material (void) phase. In the numerical experiments we use M DIO and M DIC to estimate the attained NELs on the material and void phases, respectively. Finally, we note that quality measures (54), (56), and (57) are related: Here, we use the 1-norm for all quality measures. One could, of-course, also use another norm to define these measures. The choice of norm is dictated by how much one wants small areas where minimum size conditions (50), (51), or (53) are not satisfied to influence the total quality measure. The 1-norm is rather forgiving toward having only small areas where the minimum size conditions are not satisfied, while, for example the max-norm would only focus on the largest element-wise deviation. In the binary case with exact morphological operators on a uniform mesh, the quality measures based on the 1-norm give the relative area (or volume) where the corresponding minimum size condition is not satisfied.

Discussion
As pointed out before (Sigmund 2007;Guest 2009), there are some complications regarding minimum length scales for parts of the design that are adjacent to the design domain boundaries. In particular, if ⊂ R 2 is a rectangle and B is the Euclidean unit ball, then R B ( ) = 0. Hence, M ⊂ and V = \ M cannot both possess NELs defined by expression (32). Moreover, if both phases are required to possess NELs defined by expression (32), then only one phase can occupy the region adjacent to a flat boundary. (Remember that NEL (32) and NEL (1) are essentially the same.) NEL (31) handles these issues by allowing the structuring element to be cut by , implying that the definition of local length scale varies within . Another approach to handle such issues has been presented by Xia and Shi (2015) within the level set framework, where they propose to trim the skeleton close to the design domain boundary before computing NEL (1).
Padding the design with zeros, in the optimization process, corresponds to considering a larger convex domaiñ ⊃ ⊕ rB ⊃ . DefiningṼ =˜ \ M and requiring that O˜ rB (M) = M and O˜ rB (Ṽ ) =Ṽ will, according to Theorem 1, guarantee that R˜ B (M) ≥ r and R˜ B (Ṽ ) ≥ r, respectively. In fact, it will hold that O rB (M) = M, which implies that R B (M) ≥ r. However, as will be shown in the following example (see Fig. 4 The problem of one node hinges is often discussed in relation to minimum length scale control. As illustrated in Fig. 5, imposing NELs on both phases cannot in general prevent the formation of one node hinges. Nevertheless, since four non-overlapping Euclidean balls cannot tangent in one point, imposing NELs on both phases with a Euclidean ball as structuring element excludes designs with one node hinges (for sufficiently fine grids). We note that Poulsen (2003) in a similar way concludes that the MOLE needs check monotonicity not only along the coordinate axes but also along the diagonals in order to exclude one node hinges.

Heuristic method for compliance problems
Inspired by the success of the SIMP (Bendsøe and Sigmund 2003, § 1.1.2) and RAMP (Bendsøe and Sigmund 2003, § 1.5.4) approaches, we propose a heuristic method to promote designs satisfying (53) for minimum compliance problems. We modify the standard SIMP method as follows; for a given design vector ρ, we define the physical densitỹ where ρ > 0 is small (and ensures existence of solutions to the linear equations describing the behavior of the studied system), p ≥ 1 is a SIMP penalty parameter, and O h (ρ) is a differentiable approximation of the discrete morphological opening. However, instead of using O h (ρ) in the volume constraint, we use a differentiable approximation of the discrete morphological closing C h (ρ).
(To increase readability, we have suppressed the index showing the dependence of the approximate morphological operators on the neighborhood shape.) More precisely, we define the set of admissible designs A as where v ∈ R n is a vector that holds the fractional volume (|E|/| |) of the elements, and V * is the maximum fractional volume of allowed to be occupied by material. The motivation for this approach is that intuitively, the use of the SIMP penalty parameter p > 1 promotes binary designs while the use of the open operator to define the physical design but the close operator in the volume constraint promotes designs satisfying condition (53); that is, designs with NEL on both phases. More precisely, if we use exact morphological operators, thenby relation (52) yields a relatively small contribution to the stiffness compared to its contribution to the volume constraint. Thus, designs that satisfy C N (ρ) = O N (ρ) p are promoted. For such designs, relation (52) Hence the design satisfies relations (49) and (53); that is, the design is binary and possesses NELs on both material and void.

Filtering strategy
In the numerical experiments presented in this paper, we use harmonic mean based f W -filters to approximate morphological operators. To be explicit, in expressions (59) and (60) (65) and (66) found below. Filters based on the harmonic mean were first introduced in topology optimization by Svanberg and Svärd (2013). We use equal weights within each neighborhood and thus define the weight matrix W ∈ R n×n with nonnegative entries by

, C h (ρ) is replaced with C H α (ρ), and O h (ρ) with O H α (ρ); where the harmonic close C H α (ρ) and harmonic open O H α (ρ) are defined in expressions
where D = diag{|N 1 |, . . . , |N n |} and G is a binary neighborhood indicator matrix, with elements g ij = 1 if j ∈ N i , else g ij = 0. By using the function f E H α (x) = (x + α) −1 and letting f −1 E H α denote the inverse function of f E H α , the discrete harmonic erode operator (Svanberg and Svärd 2013) with parameter α > 0 is defined as where the function evaluations are taken elementwise, that is, , the harmonic dilate operator (Svanberg and Svärd 2013) with parameter α > 0 can be written as , we get that the harmonic dilate and erode operators with parameter α > 0 satisfy That is, the harmonic morphological operators satisfy duality (48). Technically, as briefly mentioned in Section 2, the harmonic approximations only converge to the morphological operators as α → 0 for binary designs. By using the harmonic erode and dilate operator, we define the harmonic close and open operator with parameter α > 0 as and respectively.

Numerical experiments
We employ the heuristic method proposed in Section 5.1 to solve two standard test problems: minimizing the compliance of a cantilever beam and the so-called heat compliance problem. Before we present the results from our tests, we digress on the quality measures that we will use to judge our results. A typical quality measure used in topology optimization is the measure of non-discreteness suggested by Sigmund (2007). This measure is zero if the design only consists of elements with a 0 or 1 density and is strictly positive otherwise. Thus, M ND quantifies how close the design is to begin binary. Here, we use C H α (ρ) in the definition of M ND since this is the design that enters the volume constraint and thus the one that naturally would be suggested as the design to be build. We remark that if we achieve our goal to have (This is indeed the case for the optimized designs presented later in this section.) In addition to M ND , we will use the quality measures introduced in Section 4.5 that quantify how close a binary design is to satisfy NEL conditions on either or both phases of material. When computing M DOC , F DOC , M DIO , and M DIC , we approximate the morphological operators by their harmonic counterparts O H α (ρ) and C H α (ρ) with α = 10 −8 . In particular, we have a purely 0-1 design with NELs on both materials if M ND = 0 and M DOC = 0. All quality measures below are presented in per percent (%) instead of dimensionless numbers; for example, if the right hand side of expression (67) equals 0.018, then we write that M ND = 1.8%.

Cantilever beam optimization
The first test problem we consider is to minimize the compliance of a cantilever beam. The beam occupies a part of the domain ⊂ R 2 , illustrated in Fig. 6; is held fast at its left side D ; and subject to a vertical surface traction density that is uniformly distributed over F , the middle 10% of the domain's right side.
We discretize by n = n x ×n y square elements, and use the finite element method with bilinear elements to solve for the nodal displacements of the structure u ∈ R N , where N = 2(n x + 1)(n y + 1). The problem of minimizing the compliance of the structure given a limit on the available volume of material can be written as where u solves In governing equation (69), the element densitiesρ(ρ) are defined as in expression (59) with ρ = 10 −9 , K (ρ(ρ)) is the N × N stiffness matrix corresponding to the element densitiesρ(ρ), and f ∈ R N is the load vector. More precisely, letting ϕ i , i ∈ {1, . . . , N} denote the finite element basis functions, the element values of K (ρ(ρ)) are given by where (u) = (∇u + ∇u T )/2 is the strain tensor (or the symmetrized gradient) of u, the colon ":" denotes the scalar product of the two matrices, E is a constant fourth-order elasticity tensor, andρ h is an elementwise constant function with element valuesρ(ρ). Similarly, the element values of f are where t ∈ L 2 ( F ) d represents the surface traction density on the boundary portion F . To minimize the compliance of the cantilever beam, we use the gradient based optimality criteria method (Bendsøe and Sigmund 2003, §1.2) with damping parameter η = 0.5 together with a continuation approach for p, the SIMP parameter used in definition (59), and α, the nonlinearity parameter used to define the nonlinear filters. The continuation approach starts using a fixed nonlinearity parameter α = 10 and approximately solves problem (68) for p = 1, 1.5, . . . , 3. After that, the continuation continues using a fixed p = 3 and approximately solves problem (68) with α = 10 1−m/2 for m = 1, 2, . . . , 18. We refer to this scheme as the cautious continuation approach. For the first problem we employ a uniform initial guess that actively satisfies the volume constraint. At each following subproblem, we use the approximate closing of the design obtained in the previous step as initial guess. The stopping criterion used at each step was that either the maximum elementwise update of the design vector ρ was less than 0.01 or that 50 iterations had passed on that continuation step. Figure 7 shows optimized designs obtained by numerically solving minimum compliance problem (68) with volume fraction V * = 0.5 on a design domain discretized by 768 × 512 elements. The (approximately) circular neighborhoods used in the filtering are illustrated in the top right corner of each sub-figure. The radii of these neighborhoods are 4 (top row left), 6 (top row right), 8 (middle row left), 10 (middle row right), 12 (bottom row left), and 14 (bottom row right) elements, respectively. The value of RelObj, stated below each optimized beam, shows the corresponding relative objective function value; that is, the optimized beam's compliance divided by the compliance of a beam optimized without filtering and without penalty (compare with the variable thickness sheet problem). For all these results the quality measures defined above satisfy the following inequalities: M ND < 1.2 · 10 −5 %, M DOC < 5.0 · 10 −6 %, and F DOC = 0%. The small values of the quality measures, as well as the presented results, demonstrate that we indeed have minimum size control of both material phases. A closer analysis reveals that the maximum elementwise change required to get a fully 0-1 design is less than 10 −5 for all designs in Fig. 7.  Fig. 7 Cantilever beams with volume fraction V * = 0.5 optimized for minimum compliance for different filter radii using the cautious continuation approach The observant reader notices that the size of the individual beams in the optimized designs is typically larger than the imposed minimum feature size set by the filter neighborhood. We remark that when using the optimality criteria method, the best we can hope for is that it finds a local minimum to the optimization problem. The cautious continuation approach has steered the optimality critera method toward designs with slightly wider beams than imposed by the filtering. By using an aggressive approach for the continuation, using a constant p = 3 and stepping α = 10 4−3m for m = 1, 2, 3, 4 and employing the same filtering procedure as for the results in Fig. 7, yields the designs in Fig. 8. These results satisfy similar inequalities for the quality measures (M ND < 9.7 · 10 −6 %, M DOC < 3.9 · 10 −6 %, and F DOC = 0%) as those obtained when using the cautious continuation approach. Thus, we also have minimum size control of both material phases in this case. Comparing the resulting designs in Figs. 7 and 8 reveals that, for large filter radii, the two continuation approaches produce similar results. However for small filter radii, the cautious continuation approach produced design with finer features than the aggressive continuation approach. Note however that both approaches produce designs that up to the given tolerance satisfy the first order optimality conditions (the KKT conditions) of the optimization problem.  (50) is satisfied (within a given tolerance) only before the jump in M DIO , and analogously NEL condition (51) is satisfied (within a given tolerance) only before the jump in M DIC . In accordance with Theorem 3, we estimate the NEL of the material region as the largest radius where M DIO is below a certain threshold; analogously, we estimate the NEL of the void region as the largest radius where M DIC is below a certain threshold. Based on the values in Fig. 9, we estimate that R B (M) = 8 and R B (V ) = 4 for the top left design in Fig. 7; and that R B (M) = 11 and R B (V ) = 4 for the top left design in Fig. 8. In either case, the imposed lower bound on the NELs of both phases were 4.
For some of the other designs in Figs. 7 and 8, the situation appears more complicated as the quality measures start out at the lower level for small radii and then jumps up and down until they settle at the higher level for larger radii (see, for instance, the bottom two diagrams of Fig. 9). The reason for the oscillatory behavior may be that in the discrete case, unlike in the continuous case, (the discrete approximation of) a large circle may not be constructible  Table 1. Thus, for these designs the NELs of the void regions are tight on the bound, while the NELs of the material regions are not. Moreover, the deviations from the lower bounds are larger for the beams obtained using the aggressive continuation strategy.
As a final test for the cantilever optimization, we check if the results are "mesh convergent", that is, whether or not we obtain the same or at least a very similar design if the mesh is refined. Figure 10 shows optimized designs obtained by using a constant physical filter radius and the cautious continuation approach on a domain discretized by 1536 × 1024 elements (top left), 3072 × 2048 elements (top right), as well as 6144 × 4096 elements (bottom). The filter radii are 4, 8, and 16 elements on the 1536 × 1024, 3072 × 2048, and 6144 × 4096 element discretization, respectively. All these results satisfy the following inequalities for the quality measures M ND < 3.9 · 10 −6 %, M DOC < 3.8 · 10 −6 %, and F DOC = 0%. In all these cases, the maximum elementwise change required to get a fully 0-1 design is less than 4·10 −5 .
Judging from Fig. 10, the obtained designs appears mesh convergent.

Minimum heat compliance
We aim to minimize the heat compliance of a square plate that occupies the computational domain , illustrated in Fig. 11, and is subject to uniform heating. The plate is held at zero temperature along the boundary portion D  and is thermally insulated along the rest of the boundary N . For the minimum heat compliance problem, using standard SIMP with p = 3 and standard filtering schemes to impose NEL on one of the phases yields designs with large M ND (Lazarov et al. 2016). Moreover, there are numerical evidence that, for this problem-unlike the cantilever beam example-, an open-close filtering strategy (combined with adequate penalization of intermediate conductivities) fails to provide designs with NELs on both phases (Hägg and Wadbro 2017). Thus, the minimum heat compliance problem, although its apparent simplicity, is Fig. 11 Geometry for the minimum heat compliance problem rather challenging for methods that impose NELs in density based topology optimization.
We discretize by using n = n 2 x square elements, and use the finite element method with bilinear elements to solve for the nodal thermal equilibrium temperature u ∈ R N , where N = (n x + 1) 2 . The problem of minimizing the heat compliance of the structure given a limit on the available volume of material (with high thermal conductivity) is given by problem (68) together with equation system (69). Here, the physical element densitiesρ(ρ) are defined as in expression (59) with ρ = 10 −3 . The entries of K (ρ(ρ)) and f are given by respectively. Hereρ h is an elementwise constant function with element valuesρ(ρ) and ϕ i , i ∈ {1, . . . , N} are the basis functions. We use V * = 0.5 and the same cautious continuation strategy for the SIMP parameter p and the nonlinearity parameter α as employed for the cantilever beam example. As for the cantilever beam examples, we stopped when either the maximum elementwise update of the design vector ρ was less than 0.01 or when 50 iterations had passed on that step. Figure 12 shows optimized designs from the minimum heat compliance problem obtained on a design domain discretized by 512 × 512 elements. These results are optimized by using the cautious continuation approach; that is, first raising p (p = 1, 1.5, . . . , 3) with α = 10 and then decreasing α (α = 10 1−m/2 for m = 1, 2, . . . , 18) keeping p = 3. The (approximately) circular neighborhoods used in the filtering are illustrated in the top right corner of each sub-figure. The radii of these neighborhoods are 4, 6, and 8 (top row); 10, 12, and 14, (second row from top); 16, 18, and 20 (third row from top); and 22, 24, and 26 (bottom row) elements, respectively. For all these results the quality measures satisfy the following inequalities: M ND < 0.14%, M DOC < 0.011%, and F DOC < 0.0016%, that is, at most 4 elements differ in material when using the discrete open and close operator. We also remark that F DOC = 0 for eight of the optimized designs in Fig. 12.
Experiments using the aggressive continuation approach as well as tests for mesh convergence have been performed. Since the behavior is similar to that observed for the cantilever beam example, these results are not presented. We have also estimated the obtained NELs for the material and void regions, similarly as illustrated in Fig. 9 and Table 1 for the cantilever beams. In contrast to the cantilever beam problem, the NELs for the designs optimized with respect to heat compliance are tight on the imposed lower bound for the material (high conductivity) regions, while the NELs of the void (low conductivity) regions are not. Moreover, in Fig. 12 Optimized designs from the minimum heat compliance problem with volume fraction V * = 0.5 for different filter radii this case and in particular for the void regions, it is harder to give an exact estimate of the obtained NEL, since the quality measures does not have a clear gap between low and high values in all cases.
As a final experiment, we use different shapes of the filter neighborhoods for the harmonic open and the harmonic close operators. The left image in Fig. 13 shows an optimized plate obtained by using maximum volume fraction V * = 0.5, a square neighborhood (illustrated at the top right corner of this image) to impose a length scale on the material, and an octagonal shaped neighborhood to impose a length scale on the void regions. These neighborhoods are illustrated to the right of the optimized plate. This result was obtained on a domain discretized by 2048 × 2048 elements, and the quality measures for this design satisfy: M ND < 0.012%, M DOC < 0.0030%, and F DOC < 5.7 · 10 −4 %. Thus, the design in Fig. 13 satisfies the imposed lower bounds on the NELs of both phases. As indicated in Section 2, this design could be manufactured using either a square deposition tool or an octagonal punch with the same sizes as the neighborhoods in Fig. 13. However, it could not be manufactured using an octagonal deposition tool or a square punch with the same sizes as the neighborhoods in Fig. 13.

Discussion
The foundation of the heuristic method comprises the following three main components: the theory developed in Section 4; the SIMP methodolgy (Bendsøe and Sigmund 2003, § 1.1.2) that promotes binary designs; and the harmonic mean based filters (Svanberg and Svärd 2013) that Fig. 13 Optimized design from the minimum heat compliance problem with volume fraction V * = 0.5 using 2048 × 2048 elements and neighborhoods with different shapes for the length scale definition used for material and void regions provide a family of regular approximations of morphological operators. By combining these three components, we obtain a method that promotes designs with small M ND and M DOC . In a similar way as the basic SIMP method cannot guarantee convergence toward binary designs, the heuristic method used in the numerical experiments does not guarantee convergence toward binary designs with NELs on both phases. Despite that, we obtain mesh convergent essentially binary designs with independent NELs of both phases for both the minimum compliant cantilever beam as well as for the heat compliance problem-without the need of any postprocessing or projection step. We remark that it is misleading to classify the harmonic f W -mean filters, as well as any other f W -mean filter, as projection based. The reason for this stand point is that all f W -mean filters are internal, that is, Expression (73) is in contrast to projection based filtering methods that do not respect the range of values of their input.
As we now demonstrate, condition (53) relates to socalled combination filters (Sigmund 2007) where as before O h (ρ) and C h (ρ) denote differentiable approximations to the morphological opening and closing, respectively. If filter (74) produces a binary result, then O h (ρ) = C h (ρ). Thus, in the limit of exact morphological operators and when the resulting design is binary, filter (74) imposes lower bounds on the NELs of both phases. The combination filter was introduced as a low-cost substitute of open-close and close-open filters (Sigmund 2007 (Schevenels and Sigmund 2016). The idea of Guest (2009), can be captured in filters of the form where ρ M and ρ V denote two independent design vectors used for controlling the layout of material and void, respectively; D h (ρ) and E h (ρ) denote differentiable approximations of the morphological dilation and erosion, respectively. If filter (75) produces a binary result, then Thus, in the limit of exact morphological operators and when the resulting design is binary, filter (75) imposes lower bounds on the NELs of both phases. Apart from not needing two design vectors, imposing condition (53) has the advantage that the design vector itself becomes regularized in the sense that Regarding M ND , the optimized designs presented in Section 5.3 outperforms those obtained by state of the art techniques for imposing minimum length scales presented in the review by Lazarov et al. (2016). In what way the proposed method benefits from extra penalization of intermediate densities has not been fully investigated. On one hand, for elements i where 0 < O hi (ρ) < C hi (ρ) < 1, we are effectively using a higher value of the SIMP parameter. On the other hand, we note that condition (53) can be satisfied for intermediate densities; for instance c1 = O N (c1) = C N (c1) for any c ∈ (0, 1), a relation which also holds for the harmonic approximations employed in the numerical experiments.
By using M ND together with the quality measures introduced in Section 4.5, we can quantitatively verify that the optimized designs possess the imposed NELs, and in principle quantitatively estimate their actual NELs. To the best of our knowledge, quantitative estimation of the minimum length scale has not been previously employed in density based topology optimization. As mentioned in Section 2, If the design is (almost) binary, then filters that approximate the morphological dilate (erode), will guarantee NEL on the material (void) region. In that case quantitative verification, that the imposed lower bound on the NEL of one of the phases is satisfied, is superfluous. Nevertheless, quantitative estimation of the actual NELs can still be of interest. The numerical experiments indicate that how tight the actual NELs are to their lower bounds depends on the continuation strategy used for the filter parameter α and SIMP parameter p. We believe that, quantitative estimation of the MOLE could be performed, in a similar way as done here, by using the discrete MOLE functional (Poulsen 2003) instead of M DIO and M DIC .
To obtain the beams in Fig. 7 using the cautious continuation approach, the OC algorithm required 550-750 iterations, the corresponding number for the beams in Fig. 8 is 90-120 iterations. The reason for the large number of iterations is that we converged at each step in the continuation. A much faster algorithm, that produce similar results can be obtained by relaxing the termination criteria for all except the last step in the continuation approach. We remark that in particular for small filter radii, using a continuation strategy with many stages produces better results than using a continuation approach with few stages. This is in line with the results by, for example, Alexandersen et al. (2016), who remark that using a continuation scheme typically produces better results than starting with the end values. However, without imposing severe design restrictions, one cannot provide any convergence guarantees for continuation schemes in general (Stolpe and Svanberg 2001). Apart from the cost related to the continuation strategy, the extra cost (compared to the standard SIMP method) associated with the heuristic method introduced in Section 5.1 is that related to the introduction of the approximate close operator. If, as done in the numerical experiments, a cascade of generalized f W -mean filters are used to approximate the discrete morphological operators the extra computational complexity is O(n) on a regular grid using polygonal shaped neighborhoods (Wadbro and Hägg 2015;Hägg and Wadbro 2017). For general neighborhood shapes the extra computational complexity is O(n log n) by employing an FFT-based filtering strategy (Lazarov and Wang 2017).

Concluding summary
In this paper, we propose definition (31) of the NEighborhood based minimum Length scale (NEL) that is similar to definition (1) of Zhang et al. (2014); the main difference lies in the treatment of regions that are close to the boundary of the design domain. The main contribution of this paper is to rigorously establish the connection between NEL (31) and (modified) morphological operators for subsets of a convex (possible bounded) domain. We show that subsets with positive NEL are characterized as being morphologically open. A direct consequence of that characterization is that subsets whose interior and exterior both possess positive NELs are characterized by being simultaneously morphologically open and (essentially) morphologically closed. We show that, in the discretized setting of density based topology optimization, the latter translates to condition (53). We modify the SIMP-method for minimum compliance problems so that it promotes binary designs satisfying condition (53). The obtained designs are almost binary and possess positive NELs (possibly measured in different norms) on both phases of material. We make no claim toward applicability of this numerical method for other topology optimization problems-It merely serves as a demonstrator of the developed theory, which is generally applicable. The proper implementation of condition (53), in the optimization formulation of more challenging problems, is open for future investigation. The theory presented in Section 4 naturally extends to multiphase topology optimization. By requiring that each phase M i ⊂ is r i B i -open, independent minimum length scales are imposed. However, how to implement such scheme in density based topology optimization is an open problem that requires further research.
Moreover, since x ∈ B r (x) ⊂ B r (y ) for any 0 < r < r λ=1 , we have that x ∈ B r (x) ∩ ⊂ B r (y ) ∩ ⊂ M for any 0 < r < r λ=1 ; implying that E rB (M; x) = ∅ for all 0 < r ≤ r λ=1 . We have thus shown that E rB (M; x) = ∅ for all 0 < r ≤ r ≤ R B (M; x). Since this holds for each > 0, we conclude that E rB (M; x) = ∅ for all 0 < r < R B (M; x).