On the comparison of material interpolation schemes and optimal composite properties in plane shape optimization

This paper deals with the classical problem of material distribution for minimal compliance of two-dimensional structures loaded in-plane. The main objective of the research consists in investigating the properties of the exact solution to the minimal compliance problem and incorporating them into an approximate solid-void interpolation model. Consequently, a proposition of a constitutive relation for a porous material arise. The non-smoothness of stress energy functional known from the approach based on homogenization may be thus avoided which is beneficial for the numerical implementation of the scheme. It is next shown that the simplified variant of the proposed formula justifies and generalizes the RAMP (Rational Approximation of Material Properties) interpolation model of Stolpe and Svanberg (Struct Multidisc Optim 22:116–124). In the second part of the paper, explicit formulae for function θ: Ω → [0, 1] describing the distribution of basic isotropic material in the design space Ω ∈ ℝ2 are derived for various interpolation schemes by the requirement of optimality imposed at each x ∈ Ω. Theoretical considerations are illustrated by a code written in MATLAB for typical optimization problems of a cantilever and MBB beam.


Introduction
The classical optimization problem of minimal structural compliance was set in the late 1960s. Roughly speaking, the goal is to search the set of characteristic functions χ : → {0, 1}, where denotes the design area, for a minimizer χ opt such that the compliance J = J (χ) of a structure subjected to a given load achieves its minimal value J = J (χ opt ). At the early stage of research the abovementioned task had proven to be badly posed due to the general non-convergence of sequences {χ n } in the standard norm of L ∞ ( , {0, 1}) (see e.g. Kozłowski and Mróz 1969;Rozvany et al. 1982), hence χ opt could not be computed. This phenomenon, often referred to as "non-existence of classical solutions", revealed a need for regularizing the optimal design problem.
One of the regularization methods assumes replacing the classical design set L ∞ ( , {0, 1}) with its larger counterpart L ∞ ( , [0, 1]), i.e. the set of generalized designs θ whose main property is that these functions can take any value between 0 and 1. From the mechanical point of view, the extension of this type results in allowing the microstructural composites of basic material and void in the analysis of the problem. The mathematical foundation of such method, known as the homogenization theory, is being developed simultaneously with its mechanical applications from 1970s. Detailed exposition of this topic lays outside of scope of the paper, hence we refer the reader to Cherkaev (2000), Lewiński and Telega (2000), Milton (2002) and Tartar (2000) for further references.
It turns out that the shape optimization problem posed in the linear elasticity setting in two space dimensions has an explicit and accurate solution in terms of principal values N I (x), N II (x), x ∈ of the force resultant tensor N. Namely, one may analytically determine the material distribution function θ opt = θ opt (x) realizing J = J (θ opt ) as well as the local microstructural geometrical properties. Consequently, the macroscopic effective constitutive characteristics of a composite made of basic material and void locally mixed in proportions θ opt (x) and 1 − θ opt (x) respectively can also be computed. This feature is thoroughly dealt with in Allaire (2002). Nevertheless, the numerical implementation of a homogenization-based results is hampered by the non-smoothness of stress energy functional W (N, θ) for N I N II = 0. Namely, if this is the case, then the constitutive formula is not unique, hence it cannot be directly inverted into the stress-displacement form which serves as a natural base for a numerical implementation by FEM. As a possible remedy, one may either make use of the solution to the problem of optimal layout formulated for a two-material structure, see Lewiński (2004), or seek the smooth, thus invertible, approximation of W (N, θ). In the sequel, the latter option is investigated.
According to the homogenization approach, θ opt (x) ranges from 0 to 1 in the design area hence the generalized material distribution function may serve as a reference solution to certain suboptimal, though practical and manufacturable, "0-1" designs stemming from the simplified material-void interpolation models. The most popular SIMP (Solid Isotropic Material with Penalization) scheme emerged from the papers of Bendsøe (1989) or Zhou and Rozvany (1991) and it is now widely used in solving various engineering optimization problems, see (Bendsøe and Sigmund 2003) for their extensive review. The mathematical foundation of SIMP is also subject to continuous research, see e.g. Rietz (2001), Martínez (2005), Almeida et al. (2010), Amstutz (2011) and Azegami et al. (2011). Another power law-like model of the material-void composite constitutive behavior, called RAMP (Rational Approximation of Material Properties), was proposed by Stolpe and Svanberg (2001b) following Rietz (2001). Both SIMP and RAMP interpolation schemes incorporate a certain real parameter which can be adjusted during the optimization procedure to penalize the intermediate densities of the basic material in the effective composite thus the almost "0-1" designs can be created. Deep discussion of problems corresponding to penalization and numerous filtering techniques is a subject of e.g. Sigmund andPetersson (1998), Bourdin (2001), Bruns and Tortorelli (2001), Stolpe and Svanberg (2001a), Rietz (2007) and Sigmund (2007).
The constitutive properties of a composite and, consequently, the value of the energy accumulated in a particle of an effective medium depend on the material interpolation scheme. Hence, the range of the SIMP or RAMP parameters is restricted by the requirement W opt ≤ W app , where W opt denotes the stress energy related to the optimal solution based on the homogenization approach and W app stands for its approximation corresponding to the chosen interpolation scheme. Even more restrictive is the requirement of the macrostructural isotropy of the solidvoid composite. Bounds on the effective isotropic properties of a two-material mixture were introduced by Hashin and Shtrikman (1963) and improved by Cherkaev and Gibiansky (1993). It is worth pointing out, however, that both results coincide for the mixture of material and void hence in the sequel we will write W HS ≤ W app to denote the requirement of effective isotropy.
The paper is organized as follows: The notation necessary for the presentation of the topic is introduced in Section 2 together with a brief recall of optimal solution to the minimal compliance problem for material and void in a two-dimensional setting. Section 3 contains a derivation of an energy-based solid-void interpolation scheme in its general two-parameter form related to the homogenization approach. It is shown that a certain, simplified variant of the scheme allows for uniform estimation of the optimal stress energy functional and justifies the RAMP interpolation function of Stolpe and Svanberg (2001b). Hence, the proposed solid-void constitutive formula is nicknamed as generalized RAMP (GRAMP). Thus obtained scheme is next confronted with the Hashin-Shtrikman and SIMP estimates. Section 4 is devoted to the derivation of the explicit formulae for material distribution functions corresponding to various interpolation models. Numerical examples obtained with the help of GRAMP are presented in Section 5 and compared with the exact solutions and those based on SIMP and classical RAMP schemes.

Set
∈ R 2 with a Cartesian basis {e 1 , e 2 } for a middle plane of a thin plate. Let x ∈ and set t = 1 for the uniform thickness of a plate. Next, introduce the basis see e.g. Rychlewski (1995), allowing for representing the symmetric second-order tensors and Hooke's fourthorder tensors as vectors and symmetric matrices in R 3 respectively. From now on, Greek indices take values 1 or 2 while the Latin ones range from 1 to 3 and usual summation convention applies. Moreover, in the sequel we write f ,α = ∂ f /∂(x α ) for the differential operator; AB for the contraction of the Hooke tensor A (bold upright) and symmetric second-order tensor B (bold italic) and A · B for the scalar product of two symmetric second-order tensors.
Throughout the paper we frequently refer to the in-plane and isotropic compliance tensor expressed by Young's modulus E and Poisson's ratio ν. Formula ε = AN links the force resultant tensor with the deformation tensor ε(u) whose components are derived from the kinematically admissible displacement vector u by 2ε αβ = u α,β + u β,α .
Assume that the plate is loaded by the in-plane load p = p α e α . Tensor N is statically admissible and we write N ∈ S iff the equilibrium equation is satisfied for each v kinematically admissible. Let θ : → [0, 1] denote the function describing the layout of a given isotropic material A 0 in . Its total amount is limited by the isoperimetric condition The compliance of a plate related to given θ is thus expressed by or equivalently, by the Castigliano Theorem, as where W denotes the stress energy value calculated at given x ∈ , i.e. 2W (N, θ) = N · A(θ )N with N = N(x) and θ = θ(x) by abuse of notation. By substituting (4) in (3) and with N 11 = N I , N 22 = N II , N 12 = 0 standing for the principal values of N in (2) we may write Coefficients E(θ ), ν(θ) denote the constitutive properties of an isotropic and non-homogeneous plate locally composed of given material A 0 and void mixed in proportions θ, 1 − θ respectively. In what follows, the equivalent form of (9) is utilized, namely where ζ = N I /N II , N II = 0 and It is worth pointing out that all considerations in the sequel remain valid if one swaps N I with N II in (9) thus redefining ζ and the expressions in (10) accordingly. This operation is possible due to the assumed isotropy of A(θ ) and allows for the straightforward application of the subsequent results in the limiting cases of N I = 0 and N II = 0. We assume that N I = N II = 0 and W = 0 if θ = 0.
Optimal material distribution, or shape optimization problem, amounts to finding the minimal compliance of a structure through where θ ∈ L ∞ ( , [0, 1]) and denotes the Langrange multiplier for isoperimetric condition (6). Once minimizations in (12) are swapped over, it is possible that the "minimum over θ" operation is exchanged with integration due to the Rockafellar Theorem (Rockafellar 1976) which leads to As this Theorem is crucial for further considerations, let us discuss its conditions in optimal design setting. First, note that the existence of a minimizer θ min ∈ L ∞ ( , [0, 1]) in (12) is provided by the homogenization theory, see Allaire (2002, Th. 4.2.6). Moreover, structural compliance calculated for any material distribution is represented by a finite number, hence J < +∞ for given θ min . By this, one of the Rockafellar's Theorem assumptions is fulfilled.
In the original formulation of the Theorem it is expected that the minimizing function θ min is looked for in the decomposable linear space of measurable functions. Therefore, the set L ∞ ( , [0, 1]) should be extended to L ∞ ( , R) thus admitting θ min (x) < 0 and θ min (x) > 1 as possible solutions in (13). The extension, however, is formal and pose no additional problems. Indeed, by the inspection of results in the sequel one may check that the values θ min (x) and −θ min (x) saturate the extremum of the integrand in (13) at x ∈ , but the minimizer is always given by θ min (x) > 0. On the other hand, it is always possible to set θ min (x) = 1 if the upper bound of the [0, 1] interval is violated but such restriction has to be compensated by adjusting the Lagrange multiplier which is necessary to retain the isoperimetric condition (6).
The last two assumptions sufficient for (13) to hold are that the integrand is lower semicontinuous in θ for any x ∈ and that its epigraph is a measurable function. The former condition is fulfilled as a(θ ) and b(θ ) defining U (ζ, θ ), see (10), are continuous in θ for all cases considered in the remainder of this paper. Consequently, the latter condition is also satisfied.
Having disposed of this step, one may conclude that the optimal material distribution function is locally determined through the condition Functions satisfying (14) vary for different material-void interpolation schemes. Their explicit forms are derived in subsequent sections.
2.2 Properties of the explicit solution to the problem of material layout for minimal compliance From the considerations in Cherkaev (2000) and Allaire (2002) it turns out that for given (N, θ), the value of stress energy in (12) can be calculated in terms of a real parameter or explicitly in terms of principal values of tensor N as We refer the reader to the above-mentioned references for complete discussion of this result, here we only recall its certain properties that are necessary in further considerations.
Formula (15) for W opt is isotropic as it is a function of the invariants of N only. By substituting (4) with E = E 0 and ν = ν 0 in (17) and then in (16), W opt takes the equivalent form Indeed, for N I N II ≤ 0, formula (16) is given by Combining (19) and (20) in one expression leads to (18). Consequently, the optimal material distribution θ opt obtained through (14) is determined by Re-writing (18) similarly to (10) yields where, by making use of sgn(N I N II ) = sgn(ζ ), 3 Solid-void interpolation scheme based on the properties of optimal stress energy

Derivation of the interpolation function
Expressing the stress energy in the form given by (10) allows for analyzing the functional U instead of W in further considerations. Therefore, we continue in this fashion to obtain the approximation U app of U opt by taking (22) as the reference expression in determining a app (θ ) and b app (θ ) which are required to be independent of the sign of ζ and such that is satisfied for any pair (ζ, θ ). It turns out that by simplified variant of this approach one may also justify and generalize the RAMP interpolation scheme of Stolpe and Svanberg (2001b). For any θ, the sought approximation with the above required properties is realized at U app sharing the common tangent with U opt at two distinct points ζ 1 and ζ 2 . Setting ζ 1 = 1/q 1 and ζ 2 = −1/q 2 followed by leads to Requirement (24) restricts the range of q α to In this way, the two-parameter solid-void interpolation scheme is established with the approximate stress energy function W app defined through (10) and (26). It follows that U app shares the common tangent with U opt and U app = U opt at two distinct points iff q 1 = q 2 = 1.
The simplified, one-parameter variant of the above is obtained by setting q 1 = q 2 = q in (26). This choice is heuristic, however it seems quite legitimate as it leads to the uniform approximation of the functional U opt as it is shown in the sequel. As a result one may claim the generalized Rational Approximation of Material Properties (GRAMP) in the following form and (GRAMP) : where f q (θ ) denotes the common interpolation function for E(θ ) and ν(θ). Recall that the original RAMP scheme assumes ν(θ) = const. which gives The set of plots in Fig. 1 shows f q for different values of parameter q.
Assuming the solid-void interpolation scheme in a GRAMP form has a great impact on the quality of approximation U app . For given ζ 0 > 0, it can be measured as thus expressing a simple requirement of approximation uniformity for different signs of ζ . One may calculate that (U GRAMP ) = 1 which is not the case for neither of other approximation schemes discussed in this paper.

Relation to the Hashin-Shtrikman bounds and SIMP model
At this point of discussion it is worth following through the mutual relations among U GRAMP , U opt , U HS and U SIMP where the latter two stress energy functions are respectively related to the Hashin-Shtrikman and SIMP (Solid Isotropic Material with Penalization) interpolation models. Formula determining U HS can be expressed by (10) with while U SIMP admits where p stands for the penalization factor in SIMP scheme. Equation (11) links the coefficients of stress energy functionals and material properties E(θ ), ν(θ). Let us emphasize that both U HS and U opt establish certain lower limits on the stress energy accumulated in a particle of a structure. The difference in their formulae are, roughly speaking, due to the fact that the calculations of U HS incorporate the control over the isotropy of the effective material which is not the case for U opt .
Indeed, the derivation of Hashin-Shtrikman bounds consists of two independent steps. Application of a hydrostatic field N = N H E 1 allows for defining the minimal value of K (θ ), while obtaining the minimal value of L(θ ) requires subjecting a composite to the deviatoric field N = Conversely, U opt measures the energy of a composite subjected to an arbitrary field N. However, if we assume that e α stand for the eigenbasis of some N at given x ∈ , and if we set N I = N II = N (i.e. ζ = 1) then we obtain U HS = U opt , as N = ( √ 2N )E 1 is a hydrostatic field. On the contrary, setting N I = −N II = N (i.e. ζ = −1) does not lead to a similar conclusion on the equality of energies because applying the deviatoric field N = ( √ 2N )E 2 and controlling the response of a composite medium in the direction E 3 at the same time is impossible. Table 1 presents the lower bounds imposed on the parameters of GRAMP, RAMP and SIMP approximations. They stem from two requirements which has to be fulfilled by any approximate stress energy functional U app . Inequality U app − U opt ≥ 0, see (24), restricts the values of parameters to those for which the approximate energy takes the form of an isotropic function (10) while U app − U HS ≥ 0 imposes the effective isotropy of a composite material as an additional constraint. Partial results related to the Hashin-Shtrikman restrictions are reported in e.g. Bendsøe and Sigmund (1999) and Stolpe and Svanberg (2001b).
In case of the GRAMP model, the isotropy of an approximate stress function U GRAMP is satisfied for any q ≥ 1, see (27). By making use of (28) and (32) demanding in this way the macroscopic material isotropy.
which has to be fulfilled for any ζ . Consequently, one obtains q ≥ 3 as the function on the r.h.s. of (35) takes its maximum for ζ = −1. Similar calculations can be performed for SIMP and RAMP schemes. It is worth pointing out that in the GRAMP scheme parameters are independent of the basic material's Poisson ratio value ν 0 and density θ. On the contrary, in the course of calculations related to RAMP and SIMP schemes one has to set θ = 1 in order to assure the validity of expressions in Table 1 for any θ ∈ [0, 1]. The comparison of U opt , U HS and U GRAMP , U RAMP , U SIMP with parameters of lowest possible values are shown for ν 0 = 1/3 (Fig. 2) and ν 0 = −1/3 (Fig. 3). In the latter case of an auxetic material, values of parameters in the RAMP and SIMP Table 1 Requirements on the parameters in GRAMP, RAMP and SIMP interpolation schemes imposed by the conditions of optimality (U app ≥ U opt ) or material isotropy (U app ≥ U HS ) Constitutive properties predicted by Hashin-Shtrikman approach are realizable on certain isotropic microstructures like 3rd rank sequential laminates (Francfort and Murat 1986) or coated circles (Hashin 1962;Grabovsky and Kohn 1995a). On the other hand, composites of minimal compliance can be arranged as orthotropic 2nd rank orthogonal sequential laminates or, if det N > 0, in a form of the Vigdergauz microstructures (Vigdergauz 1989;Grabovsky and Kohn 1995b) or 4th rank sequential laminates (Allaire and Aubry 1999).
In light of the studies in this section it seems that one may neglect the requirement of the macrostructural isotropy if the best possible macroscale material distribution is considered as a primal goal of the optimization problem. Nevertheless, thus obtained solution may serve as a good starting point of the continuation method where the change of a parameter in the interpolation scheme (like e.g. q in RAMP or GRAMP) reflects in the tendency to the pure "0-1" distribution of a basic isotropic material in the domain of an analyzed structure. Problems related to the microscopic layout of material and void in optimized structures are detailed at length in e.g. Allaire (2002), Cherkaev (2000) and Bendsøe and Sigmund (2003), see also references therein.

Material distribution functions for minimal compliance
Our next objective is to derive the material distribution functions that correspond to various interpolation schemes. Repeated application of (14) where W (N, θ) is given by (10) and appropriate expressions for a(θ ) and b(θ ) enables us to write the sought formulae in their explicit forms similarly to (21). Indeed, for (28) we obtain the function whereas for (30), in the original RAMP scheme and the Hashin-Shtrikman framework with (32) provides θ HS (N) = min 1, Material distribution function predicted by the SIMP model with (33) is given by It follows that functions θ(N) in SIMP, RAMP and Hashin-Shtrikman interpolation schemes are sensitive to the sign of ζ and ν 0 which is not the case for the material distribution based on GRAMP or exact homogenization approach. Figure 4 illustrates the material distribution functions with the parameters set to the values for which the respective functionals U GRAMP , U RAMP and U SIMP do not violate the lower bounds on the stress energy given by U opt and U HS for ν 0 = 1/3. In this sense Fig. 4  Explicit formulae for material distribution can be directly applied in numerical codes for compliance minimization problem in this way allowing to avoid the heuristic density updating schemes.

General comments on numerical implementation
The GRAMP interpolation scheme is applied to find the suboptimal layout of basic isotropic material in classical problems of a cantilever and half of a MBB beam. Results obtained in this way are compared to the distribution of materials ensuing from the homogenization theory as well as those related to the original RAMP or SIMP approach. Designs associated with both exemplary problems are shown in Figs. 6, 7, 8 and 9. The dimensions of a rectangular design area are given by l × h = 2 × 1 and the material fraction of is set to m = 0.5. Young's modulus of an isotropic solid is denoted by E 0 = 1 and two cases of Poisson's ratio are considered in the course of calculations, i.e. ν 0 = ±1/3. The structure is due to a concentrated load P = 1. All calculations are coded in MATLAB and performed by FEM using the mesh of 200 × 100 square, four-node elements with two degrees of freedom per node. The idea of programming is to supplement the applications published at the TopOpt research group website (www.topopt.dtu.dk) hence the codes presented in the Appendices have a structure based on the program written by Andreassen et al. (2011) for topology optimization in 2D. Set then the optimization routine is stopped. Otherwise, the material distribution function is updated by applying (21) or one of the formulae (36)-(39) depending on the interpolation scheme. In this way vector x k is determined. Finally, the constitutive properties of a composite are modified and the (k + 1)-th step of the routine is executed. The l.h.s. in (40) can be understood as a counterpart of the L 2 norm of function N posed in a discrete setting.
Here we set = 10 −2 but this value can be adjusted to meet the required accuracy. It has to be pointed out that the convergence measured by (40) is fairly slow. Some authors prefer to rewrite this criterion in terms of the material density vector x, see the numerical codes in Bendsøe and Sigmund (2003) and Andreassen et al. (2011).
Minimal compliance problem tackled in this paper falls into the category of constrained optimization. The necessary conditions for optimum are thus obtained by applying the Karush-Kuhn-Tucker (KKT) Theorem, see e.g Haftka and Kamat (1985) for further details in context of structural optimization. Discussed formulation (12) is a double minimization with respect to N and θ. The KKT conditions of optimality admit the Lagrange multiplier for the isoperimetric constraint (6). The requirement of statical admissibility imposed on N by (5) can be accounted for in a similar way with the Lagrange multiplier assumed in a functional form. Both constraints are to be satisfied at the stationary point of an objective function. In this way, the stopping criterion of the numerical algorithm is provided.
In the corresponding numerical algorithm the unknowns N and θ are respectively referred to as vectors N and x. Iterative numerical implementation of (12) with alternating minimization in both variables guarantees that the objective function converges to a stationary point. Indeed, determining N k for fixed x k−1 reduces to solving the linear elasticity problem. Consequently, the stress energy is minimized hence the value of the objective function decreases. Upon updating the stress resultant vector, the objective function is minimized explicitly by x k obtained through (21) or one of (36)-(39) depending on the solid-void interpolation model. The Lagrange multiplier is adjusted in order to enforce the corresponding KKT condition. Formula (40) serves as the stopping condition of a numerical algorithm. Its fulfilment shows the convergence of N and the convergence of the objective function follows. Tables 2 and 3 give some details concerning the optimization procedure implemented in MATLAB for the half-MBB beam problem. The results were obtained using a DELL XPS M1330 laptop with Intel Core2 Duo CPU (model T9300), 4 GB RAM, Windows 7 (64-bit) and MATLAB 7.11.0 (R2010b).

Comments on the homogenization-based approach
Stress functional U opt obtained by the homogenization approach is continuous but not smooth at ζ = 0, see Figs. 2 and 3. This, in turn, poses a serious drawback of direct application of the stress-based results in numerical computations. However, the results reported in Lewiński (2004) for a two-material problem make it possible to reformulate the topology optimization problem in the displacement setting. Consequently, for ε I , ε II denoting the principal values of a deformation tensor ε and for k 0 = 1/K 0 , μ 0 = 1/L 0 one may set followed by thus defining the optimal stiffness tensor C opt represented by [C opt i j ] = diag 2k opt , 2μ opt , 2μ opt in basis (1). In order to reduce the checkerboard-like instability of solutions, all element-wise calculations in the homogenization-based approach are performed at four integration points. The stiffness matrix of e-th element is thus determined according to where b ji corresponds to i-th integration point and denotes the j-th row of a B i matrix containing the derivatives of element shape functions. Corresponding numerical code is presented in Appendix A, e.g. the composite solution in Fig. 8b is obtained by executing tophomog4(200, 100, 0.5), where the parameters respectively denote the number of elements in the horizontal and vertical directions and the solid fraction of the design area.

Comments on the GRAMP approach
Approximate function U GRAMP is smooth, see Figs. 2 and 3 hence it is invertible in the whole range of ζ . For the purpose of this paper, two MATLAB implementations of the GRAMP scheme were coded. The difference between them lies in the method of element stiffness matrix integration and in the number of points in which strains and stresses are computed. This, in turn, leads to the difference in obtained solutions. Application of the first code results in composite designs and it is similar to the code based on the homogenization approach. Slight modifications are detailed in Appendix A. This implementation of the GRAMP scheme, used for mimicking the composite solution of tophomog4, is illustrated in Fig. 8c and obtained by executing topgramp4(200, 100, 0.5, 1), where the last parameter stands for the GRAMP penalization factor q = 1. In the second code, presented in Appendix B, the element stiffness matrix K e is calculated as where and see (11) and (28). Strains and stresses are calculated in the middle of each element. The updated distribution of material is determined by making use of (36).

Comments on the filtered designs
For practical purposes it is desirable that the intermediate densities are filtered out from the final design. Application of such filter allows not only for obtaining nearly blackand-white material layouts but also helps in avoiding the checkerboard instability of numerical calculations. Roughly speaking, filtering the density θ(x) at given x ∈ works by averaging it over the neighborhood of given radius r . By this, the tendency to reproduce the fine-scale arrangement of solid and void is bypassed and composite regions are introduced in . In turn, these regions are penalized by proper parameter adjustment in the material interpolation formula. As a result, the topology of final design becomes nearly "0-1" and the width of transition between void and solid areas depends on the filtering radius r .
According to e.g. Bendsøe and Sigmund (2003), implementing a filtering technique in the formulation of optimal design problem at hand does not impose any additional limit on the space L ∞ ( , [0, 1]), i.e. the space of material distribution functions. Density filter can be applied at each x ∈ independently as a part of the material redistribution procedure.
Numerous filtering methods are reviewed in Sigmund (2007) and some issues related to the programming in MATLAB are tackled in Bendsøe and Sigmund (2003) and Andreassen et al. (2011). In the present paper we make use of the built-in MATLAB convolution function conv2, see MathWorks (2011), adjusted to the cone-shaped density filter of Bruns and Tortorelli (2001) and Bourdin (2001).
For simplicity, assume that the design area is mimicked by a k × l regular mesh of square elements with unit side length. Let the indices (i, j) respectively determine the horizontal and vertical location of the element and write θ (i, j) for the material density. Density θ mod (i, j) modified by making use of conv2 function is computed as the weighted average where and d (m,n) denotes the distance between the centers of given element (i, j) and the one located at (i + m, j + n) while r min stands for the filtering radius.
The exemplary solution in Fig. 8e is acquired by calling topgramp1 (200, 100, 0.5, 3, 1.5, 1), where the last three parameters respectively stand for the value of the GRAMP penalization factor, the filtering radius and a flag for filter application. It is worth noting that by setting r min = 1.5 the material density in given element is calculated as a weighted average over all adjacent ones. For the layouts related to the original RAMP and SIMP schemes, shown in Fig. 8d, f, the topgramp1 code requires slight modifications whose summary is given in Appendix B.

Conclusions
Considering the formula W opt (N, θ) for stress energy of a material-void mixture as a reference expression in a research for its approximation is motivated by the nonsmoothness of W opt for det N = 0. Consequently, a differentiable, two-parameter estimate W app (N, θ) is obtained. Additional requirement of approximation uniformity for different signs of det N and arbitrary θ ∈ [0, 1] leads to the simplified variant of a scheme. It is nicknamed GRAMP, as it justifies and generalizes the RAMP model of Stolpe and Svanberg (2001b) in a simple way, easy to implement in numerical subroutines. Moreover, besides the simplified composite solutions to the minimal compliance problem, GRAMP is able to produce the almost "0-1" layouts of basic material as a result of the continuation method possibly combined with certain design smoothing routines. However, application of these techniques are dealt with in this paper only in limited scope.