Non-equal weighting within neighborhoods can be achieved by sequentially applying the same equally weighted fW-mean filter twice (or more). More precisely, when applying an equally-weighted fW-mean filter twice, we have that
$$ \boldsymbol{F}(\boldsymbol{F}(\boldsymbol{\rho})) = \boldsymbol{f}^{-1}\left( \boldsymbol{W}^{2} \boldsymbol{f}(\boldsymbol{\rho}) \right). $$
(37)
The weight matrix W 2 has, in general and for each neighborhood, weights that decay from the neighborhood center. Figure 1 illustrates the weights corresponding to the first three powers of W = D −1 G, from top to bottom W, W 2, and W 3, for an octagonal (left column) and a rectangular neighborhood (right column). If the neighborhood shape is a convex polytope \(\mathcal {P}\subset \mathbb {R}^{d}\), then it can be shown that a neighborhood shape corresponding to W 2 is \(2\mathcal {P}\).
Fig. 1Illustration of the first three powers (from top to bottom) of the weight matrix, W = D −1 G, W 2, W 3, for an octagonal and a rectangular neighborhood (left and right column, respectively)
Full size image
We remark that, in general, evaluating the sums Gf(ρ) accounts for a significant portion of the computational cost of the filter application. Moreover, the computational effort required to evaluate these sums grows with the complexity of the neighborhood polytope. Thus, if one wishes to apply a filter with weights that decay with the distance from the neighborhood center, particularly in three space dimensions, one could save a great portion of the computational time by selecting a simple neighborhood. For example, the computational complexity for filtering over a box shaped neighborhood is approximately 10 times lower than the complexity for filtering over the significantly more complex rhombicuboctahedron (a polytope with 26 faces—twelve rectangular, six square, and eight triangular faces) neighborhood (Wadbro and Hägg 2015).
On the other hand, if one wishes to use one of the nonlinear filters that are designed to mimic min or max operators over the neighborhood, then the neighborhood shape is important. In this latter case, one cannot use the weighted version with a square or box shaped neighborhood to approximate a circular or spherical neighborhood, since the nonlinearity of the filtering process will essentially pick out the maximum/minimum of the element values in each neighborhood. However, by using neighborhoods of different but still relatively simple shapes, such as two square shaped neighborhoods (with a relative rotation of 45 degrees), one can get a filtering procedure where the neighborhood around each element is an octagon. A similar approach can also be used for the three dimensional case, where one can use four cubic neighborhoods (the scaled unit cube plus three other cubes, rotated 45 degrees in the x 1 x 2, x 1 x 3, and x 2 x 3 planes, respectively) to approximate a sphere. The use of a sequence of filters with simple neighborhood shapes streamlines the implementation. Moreover, the resulting memory access pattern can be made very regular, which paves the way for future highly efficient parallel implementations.
If the neighborhood shape and the weighting (possibly nonuniform) within each neighborhood is the same for all neighborhoods, then Wf(ρ) corresponds to a convolution. In this case fW-mean filters may be efficiently applied using an FFT based algorithm (Lazarov et al. 2016). The main idea behind FFT based filtering is that the convolution between the density and the filter kernel corresponds to elementwise multiplication in frequency domain. The application of an fW-mean filter can be performed as
$$ \boldsymbol{F}(\boldsymbol{\rho}) = \boldsymbol{f}^{-1}\left( \boldsymbol{\mathcal{F}}^{-1} \big\{ \boldsymbol{\mathcal{F}}\{\boldsymbol{w}\} \odot \boldsymbol{\mathcal{F}}\{\boldsymbol{f}(\boldsymbol{\rho})\} \big\} \right), $$
(38)
where \(\boldsymbol {\mathcal {F}}\) and \(\boldsymbol {\mathcal {F}}^{-1}\) represent the d-dimensional FFT and inverse FFT transform, respectively, ⊙ denotes the element-wise product, and w is a vector that represents the filter kernel. The asymptotic computational complexity of FFT based filters is O(nlogn) independent of the complexity of the neighborhood shape. For filtering using any given polygonal neighborhood, the FFT-based filters are hence asymptotically slower than the O(n) algorithm based on recursively adding and subtracting moving sums. For large-scale problems, with 105– 109 elements, the FFT based algorithm is 1.3–6.5 times slower than the moving sums based algorithm on a standard desktop equipped with an Intel Xeon E5-1650 v3 CPU. The FFT algorithm is potentially sensitive to round-off errors, similarly as the O(n) moving sums based algorithm (Wadbro and Hägg 2015). In many fields the FFT method is standard and there exists many highly optimized and parallel FFT routines, thus the FFT-based filters are a competitive alternative.
The impact of the filter on the sensitivities is found by computing v T∇F(ρ) for some vector \(\boldsymbol {v}\in \mathbb {R}^{n}\). In practice, v is the gradient of the objective or a constraint function with respect to the filtered densities. By using expression (35) we find that
$$\begin{array}{@{}rcl@{}} \sum\limits_{i=1}^{n} v_{i}\frac{\partial F_{i}(\boldsymbol{\rho})}{\partial \rho_{j}} &=& \sum\limits_{i=1}^{n} v_{i} w_{ij}\frac{f^{\prime}(\rho_{j})}{f^{\prime}\left( F_{i}(\boldsymbol{\rho}) \right)}\\ &=& f^{\prime}(\rho_{j})\sum\limits_{i=1}^{n} w_{ij}\frac{v_{i}}{f^{\prime}\left( F_{i}(\boldsymbol{\rho}) \right)}, \end{array} $$
(39)
that is, to modify sensitivities we need to carry out matrix multiplication by W T. In the special case of equal weighting within neighborhoods, multiplication by W T translates to multiplication by G T; or expressed differently, to perform summation over the transposed neighborhoods \(\mathcal {N}_{i}^{T} = \{j:i\in \mathcal {N}_{j}\} = \{j:g_{ji}=1\}\). If the neighborhoods are symmetric, then G T = G and the same summation algorithm can be used for both filtering and sensitivity calculation, which facilitates the implementation.
Assume that the neighborhoods are defined by a neighborhood shape \(\mathcal {N}\) so that \(\mathcal {N}_{i}= \{j:x_{j}-x_{i}\in \mathcal {N}\}\), where x i and x j denote the centroid of elements i and j, respectively. For each element \(j\in \mathcal {N}_{i}^{T}\), we have that \(i\in \mathcal {N}_{j}\) and by definition \(x_{i}-x_{j}\in \mathcal {N}\); that is, there exists \(y\in \mathcal {N}\) so that x i −x j = y or equivalently x j −x i = −y. Since all steps above are bidirectional, we have that \(\mathcal {N}_{i}^{T} = \{j:x_{j}-x_{i}\in - \mathcal {N}\}\); hence a neighborhood shape which defines the transposed neighborhoods is found by inversion of \(\mathcal {N}\) in the origin. Since \(\mathcal {P}\) and \(-\mathcal {P}\) are essentially the same polytope, this means that implementing the fast summation algorithm over \(-\mathcal {P}\) requires the same amount of work as implementing it over \(\mathcal {P}\).
If we work with a Cartesian grid and the design variables are stored using a standard slice/fiberwise numbering (row- or column-wise in the two dimensional case), then we can use that G T = PGP, where P is the flip or exchange matrix. That is, P is the matrix with ones along the anti-diagonal and zeros elsewhere. Hence, the same summation algorithm may be used both for filtering and sensitivity evaluation.
It has already been established that sequential application of filters is a means to arrive at filters with desirable properties, see for instance the open and close filters introduced in topology optimization by Sigmund (2007).
Assume that we are given a family of N fW-mean filters,
$$ \boldsymbol{F}^{(K)}(\boldsymbol{\rho}) = \boldsymbol{f}_{K}^{-1}\left( \boldsymbol{W}^{(K)}\boldsymbol{f}_{K}(\boldsymbol{\rho}) \right), K\in\{1,\ldots,N\}. $$
(40)
For each K∈{1,…, N}, we define the cascaded filter function of order K, C (K):[0,1]n→[0,1]n, to be the composition
$$ \boldsymbol{C}^{(K)} = \boldsymbol{F}^{(K)} \circ \boldsymbol{F}^{(K-1)} \circ {\ldots} \circ \boldsymbol{F}^{(1)}, $$
(41)
and we define
$$ \boldsymbol{\rho}^{(K)} = \boldsymbol{C}^{(K)}(\boldsymbol{\rho}). $$
(42)
The cascaded filter is naturally applied sequentially
$$ \boldsymbol{\rho}^{(K)} = \boldsymbol{F}^{(K)}\left( \boldsymbol{C}^{(K-1)}(\boldsymbol{\rho})\right) = \boldsymbol{F}^{(K)}\left( \boldsymbol{\rho}^{(K-1)}\right), $$
(43)
where we set ρ (0) = ρ. If the prerequisites of the algorithm presented in Wadbro and Hägg (2015) are met for all N fW-mean filters in the cascade (41), then sequential application of the cascaded filter requires O(Nn) operations.
We proceed to show how to modify sensitivities in the case of a cascade of N fW-mean filters. That is, given a vector \(\boldsymbol {v}\in \mathbb {R}^{n}\), we want to compute v T∇C (N)(ρ). To this end, we assume that N≥2, let v (N) = v, combine expressions (42) and (43), and apply the chain rule
$$\begin{array}{@{}rcl@{}} && \sum\limits_{i=1}^{n} v_{i}^{(N)}\frac{\partial}{\partial \rho_{j}}C_{i}^{(N)}(\boldsymbol{\rho})\\ &&= \sum\limits_{i=1}^{n} v_{i}^{(N)}\frac{\partial}{\partial \rho_{j}}\left( F_{i}^{(N)}\left( \boldsymbol{C}^{(N-1)}(\boldsymbol{\rho})\right) \right) \\ &&=\sum\limits_{i=1}^{n} v_{i}^{(N)} \sum\limits_{k=1}^{n} \frac{\partial F_{i}^{(N)}}{\partial \rho_{k}}\bigg\rvert_{\boldsymbol{\rho}^{(N-1)}}\frac{\partial}{\partial \rho_{j}}C_{k}^{(N-1)}(\boldsymbol{\rho})\\ &&=\sum\limits_{k=1}^{n} \left( \sum\limits_{i=1}^{n} v_{i}^{(N)} \frac{\partial F_{i}^{(N)}}{\partial \rho_{k}}\bigg\rvert_{\boldsymbol{\rho}^{(N-1)}}\right) \frac{\partial}{\partial \rho_{j}}C_{k}^{(N-1)}(\boldsymbol{\rho})\\ &&=\sum\limits_{k=1}^{n} v_{k}^{(N-1)}\frac{\partial}{\partial \rho_{j}}C_{k}^{(N-1)}(\boldsymbol{\rho}), \end{array} $$
(44)
where we in the last step have defined
$$ v_{k}^{(N-1)} = \sum\limits_{i=1}^{n} v_{i}^{(N)}\frac{\partial F_{i}^{(N)}}{\partial\rho_{k}}\bigg\rvert_{\boldsymbol{\rho}^{(N-1)}}. $$
(45)
The first and last lines of expression (44) are similar in form. However, the order of the cascade in the last line is one less than that in the first line. The steps in expression (44) may be repeated until the order of the cascade in the last line is 1, that is, the cascade consists of just one filter. We conclude that evaluating v T∇C (N)(ρ) corresponds to sequentially computing the vectors v (K−1)
$$ v_{k}^{(K-1)} = \sum\limits_{i=1}^{n} v_{i}^{(K)}\frac{\partial F_{i}^{(K)}}{\partial\rho_{k}}\bigg\rvert_{\boldsymbol{\rho}^{(K-1)}}, $$
(46)
for K = N, N−1,…,1 and where as before v (N) = v. The sum on the right hand side of expression (46) is precisely of the form (39).
Hence, if we store the intermediate designs ρ (K), K∈{1,…, N} when sequentially applying the filter and if the prerequisites of the fast summation algorithm are fulfilled for all filters in the cascade, then modification of sensitivities requires O(Nn) operations. The low computational complexity for modifying sensitivities comes at the cost of O((N−1)n) memory used to store the intermediate designs ρ (K), K∈{1,…, N−1}.
Existence of solutions to the continuous penalized minimum compliance problem in the case when a cascade of fW-mean filters is applied can be proven by following the same reasoning as in the proof in Section 2.
In the following section, we present numerical experiments performed in Matlab. For the cantilever beam experiments, we use a modified version of the 2D multigrid-CG topology optimization code by Amir et al. (2014). Below, we describe the major changes done to the code. First, we introduce a filter struct filterParam to hold all information needed to perform filtering and sensitivity modification, see Table 1. In the code excerpts, we have suppressed the struct name filterParam to increase the readability. For instance, we simply write N instead of writing filterParam.N. Listings 1–3 contain the new parts of code that needs to be included in order to use a filtering procedure composed of a cascade of generalized fW-mean filters. The Matlab code in Listing 1 computes the neighborhood sizes. The Matlab code in Listing 2 filters the vector rho by using the filterParam struct and the procedure outlined in Section 4.3. The observant reader notices that in fact it is not ρ (K) that is saved in the filtering, but (D (K))−1 G (K) f K (ρ (K−1)) since this enables the use of generalized fW-mean filters where \(g_{K} \neq f_{K}^{-1}\). The Matlab code in Listing 3 computes v T∇C (N)(ρ), as outlined in Section 4.3. We remark that filtering of densities must be moved inside the optimality criteria update whenever a non volume-preserving filter is used.
Table 1 The fields in thefilterParam
struct. Note that fieldN
and indexk
corresponds to the variables N and K, respectivelyFull size table
Listing 1Matlab code that computes the neighborhood sizes
Full size image
Listing 2Matlab code that filters the design variables by using a cascade of generalized fW-mean filters
Full size image
Listing 3Matlab code that modifies the sensitivities with respect to the physical design following the description in Section 4.3
Full size image