SIMP‐ALL: A generalized SIMP method based on the topological derivative concept

Topology optimization has emerged in the last years as a promising research field with a wide range of applications. One of the most successful approaches, the SIMP method, is based on regularizing the problem and proposing a penalization interpolation function. In this work, we propose an alternative interpolation function, the SIMP‐ALL method that is based on the topological derivative concept. First, we show the strong relation in plane linear elasticity between the Hashin‐Shtrikman (H‐S) bounds and the topological derivative, providing a new interpretation of the last one. Then, we show that the SIMP‐ALL interpolation remains always in between the H‐S bounds regardless the materials to be interpolated. This result allows us to interpret intermediate values as real microstructures. Finally, we verify numerically this result and we show the convenience of the proposed SIMP‐ALL interpolation for obtaining auto‐penalized optimal design in a wider range of cases. A MATLAB code of the SIMP‐ALL interpolation function is also provided.


INTRODUCTION
In the recent years, topology optimization has become a research field by itself, not only because it arouses interest on the scientific community but also because it provides successful solutions to industrial problems. 1 Traditionally, topology optimization seeks to design light-weight structural components without losing stiffness capabilities. In the last years, it addresses other fields of applications with equally success. For that propose, several approaches have emerged in the last decades.
Initially, topology optimization has been successfully addressed through regularization techniques. The SIMP method represents probably the most celebrated strategy. The characteristic function, usually used to describe the weak (white) and stiff (black) subdomains, is regularized allowing intermediate values, or in other words, allowing the presence of gray areas. Frequently, this regularized representation of the characteristic is commonly named density function. Usually, the challenge is to characterize the material properties of the intermediate values of the density function (gray areas). This characterization is commonly described by an interpolation function. As referred in the work of Bendsøe and Sigmund, 2 the SIMP method proposes a polynomial interpolation function that, in some cases, lies inside the Hashin-Shtrikman (H-S) bounds. This could lead to gray areas that sometimes cannot be interpreted as microstructures. Additionally, depending on the nature of the materials, the SIMP method could propose a weak interpolation penalization, and consequently, large gray areas could appear.

The exact topology optimization problem
Let Ω be a fixed domain split into Ω 1 and Ω 0 with enough regularity, fulfilling Ω 1 ∪ Ω 0 = Ω and Ω 1 ∩ Ω 0 = ∅. Let the boundary of Ω be defined as Γ ∶= Ω and split into Γ D and Γ N and let the spaces  , , and  be defined as Then, let us define the characteristic function ∈  and the fourth-order tensor C( ) as The constitutive tensors C 0 and C 1 in plane stress linear elasticity are defined as C 1 = 2 1 I + ( 1 − 1 )I ⊗ I and C 0 = 2 0 I + ( 0 − 0 )I ⊗ I, where I and I represents the fourth-and second-order identity tensors and 1 , 1 and 0 , 0 are the bulk and shear modulus of the material in Ω 1 and in Ω 0 , respectively. Introducing the bilinear and linear form of the linear elasticity problem as the exact topology optimization problem can be then defined as follows: find ∈  and u ∈  such that min.
,u J( , u) = l(u) where, in this case, the cost J( , u) stands for the compliance function.

The relaxed topology optimization problem
The difficulties of solving the exact problem are usually circumvented by proposing a relaxed version of it, this is by replacing the characteristic function with a density-like function (x) ∈ [0, 1], which, in contrast with the characteristic function, can take intermediate values. Thus, the exact topology optimization problem can be rewritten as the following regularized topology optimization problem: find ∈  = L ∞ (Ω, [0, 1]) and u ∈  such that min.
,u J( , u) = l(u) where the bilinear and linear form are now defined as In this case, the constitutive tensor takes the following form: are interpolation functions that must be proposed. Computation of the compliance gradient. Let the compliance function l(u) be defined as in (7) and let u ∈  be the solution of the elasticity problem defined in (6), then the gradient of the compliance takes the following form: where C ′ ( ) is the derivative of the constitutive tensor defined in (8). To see this result, let us define the implicit function and let us compute its differential. From ellipticity arguments, 15 we can state that ∀ v ∈  and ∀ ∈  that provides a coercive constitutive tensor, there exists a unique u ∈  that solves a( , u( ), v) = l(v). Thus, defining a density variatioñ ∈ , we can state that F( +̃) = F( ) = 0, and consequently, where D a( , u( ), v) is the Fréchet derivative with respect to the first argument of the bilinear form a ( , u, v) and Du( ) is the Fréchet derivative of u at . Then, the differential of the compliance can be computed as where we have used the symmetric behavior of the bilinear form and the result of (11). Note that, for neglecting higher-order terms, we have supposed that the reminder o(||̃|| 2 ) is small enough. The proof ends relating the gradient g( ) of the compliance with its Fréchet derivative Dl(u) by the Riesz representation theorem, ie, Dl( )̃= ∫ Ω g̃.

TOPOLOGICAL DERIVATIVE
On a given domain (unperturbed) Ω ⊂ R 2 , let us insert in a certain pointx a small circular inclusion B of radius > 0. We call the new (perturbed) domain Ω (see in Figure 1).

Definition 1.
Let J(Ω) and J (Ω) be the value of a function in the unperturbed and perturbed domain. For a certain positive function f( ) such that f( ) → 0 when → 0, typically, f( ) = 2 , we define the topological derivative of J(Ω) as the following limit: Remark 1. Note that to make expression (13) well defined, we need to ask that function J (Ω) admits the following asymptotic expansion: Definition 2. For a given Young's modulus and Poisson's ratio of the matrix Ω ⧵ B and the inclusion B , represented by E m , m and E i , i respectively, we define the fourth-order polarization tensor P as

FIGURE 1
Topological derivative concept where the parameters p 1 and p 2 take the following values: and the coefficients , , , 1 , 2 , and 3 are of the form Proposition 1. Let u ∈  be the solution of the elasticity problem defined in (6), let the compliance l(u) be defined as in (7), and let the polarization tensor P be expressed as in Definition 1, then the topological derivative at the pointx, defined as in Definition 1, is as follows: where ∇ s is the symmetric gradient operator and = C m ∶ ∇ s u the stresses with C m the constitutive tensor of the matrix domain.
Proof. This result is not trivial and is obtained after proposing an asymptotic expansion of the solution of the topologically perturbed problem and after analytically solving an associated exterior problem. See the work of Novotny and Sokołowski 16 for full details. For the same result but for any value of the Poisson ratio, see the work of Giusti et al. 17 For an easier comparison with the H-S bounds, we seek to rewrite the topological derivative only in terms of the shear and bulk modulus. .
Then, the topological derivative inx can be written only in terms of the shear and bulk modulus (of the matrix and inclusion) as follows: where the parameters d and d are To obtain this expression, let the shear and bulk polarization coefficients q and q be computed first as Inserting the lame parameters relations for plane stress E = 4 + , and = − + into the coefficients defined in (17) and (18) and replacing them in (16), we can write the shear and bulk polarization coefficients as follows: which are precisely q( m , I , m ) and q( m , I , m ). Then, replacing relation (23) in (15), we can directly rewrite the polarization tensor as P = p I + 1 2 (p − p )I ⊗ I. The proof ends by replacing this last relation and the constitutive tensor of the matrix C m = 2 m I + ( m − m )I ⊗ I in the topological derivative expression (19).
Remark 2. We can clearly identify two different cases.
(i) Stiffer inclusion insertion: We define this case by saying that, in a pointx ∈ Ω 0 with parameters m = 0 , m = 0 , we add an inclusion with parameters I = 1 and I = 1 . Then, the topological derivative g T 0 in this scenario takes in the pointx the following form: with (ii) Weaker inclusion insertion: Similarly, we define the opposite case by identifying in a pointx ∈ Ω 1 the matrix parameters as m = 1 , m = 1 and the parameters of the inclusion as I = 0 and I = 0 . Then, the topological derivative g T 1 takes in the pointx the following form: with where Remark 3. Following the seminal work of Amstutz, 13 one could think on proposing an interpolation scheme such that the gradient in Ω 1 and Ω 0 coincides precisely with the topological derivative. We distinguish again both scenarios.
(i) Stiffer inclusion insertion: Takẽ= 2 in B(x, ) and zero otherwise, where B(x, ) represents a small circular ball of centerx and radius . This means that, when replacing a small circular domain of the matrix by a small circular inclusion, we add a small amount of density. Then, if we want that Taylor's expansion (12) matches with the asymptotic expansion (14), we must define an interpolation such that, whenx ∈ Ω 0 ( = 0), (ii) Weaker inclusion insertion: In this case, we takẽ= − 2 in B(x, ) and zero otherwise. This means that, when replacing a small circular domain of the matrix by a small circular inclusion, we subtract a small amount of density. Similarly, to match both expansions, we should define an interpolation such that, whenx ∈ Ω 1 ( = 1),

CONNECTION BETWEEN THE TOPOLOGICAL DERIVATIVE AND HASHIN-SHTRIKMAN BOUNDS
H-S bounds. Let us suppose that a composite material with shear and bulk constitutive properties H and H is composed by two-phase materials with a fraction volume ∈ [0, 1]. Let us assume that the shear and bulk modulus of both constituents of the composite are 0 and 0 and 1 and 1 and satisfies 1 > 0 . Then, the shear and bulk modulus of the composite material H and H have an upper and lower bound regardless the arrangement of both phases. In addition, those bounds, commonly called isotropic H-S bounds, can be written, in the case of the shear modulus, as and similarly in the bulk modulus case as Those bounds were developed in the work of Hashin and Shtrikman. 18 See also reference books by Bendsøe and Sigmund 19 and Allaire. 20 Remark 4. Note that all bounds can be expressed in the following general form (generalized H-S bounds): where the shear and bulk H-S bounds are recovered by choosing the values described in the following table.
Remark 5. Note that the derivative of the H-S bounds in both components can be rewritten as follows: In fact, this way of rewriting the H-S bounds derivative can be obtained after simple calculations. From Equation (33), the generalized H-S bound can be rewritten as and consequently, the derivative of the generalized H-S bounds is expressed as Thus, its value in the components is just Inserting (34) into (37) and using definition of q in (20), the derivative of the H-S bounds for the shear modulus in Ω 0 ( = 0) and in Similarly, the derivative of the H-S bounds for the bulk modulus in Ω 0 ( = 0) and in Ω 1 ( = 1) can be written as Topological derivative as H-S bounds derivative. This is one of the two main results of this work. Let us define two material components with shear and bulk material properties 0 , 0 and 1 , 1 , respectively, and let us propose an interpolation function for modeling its mixture by means of its fraction volume = [0, 1] such that close to Ω 0 coincides with the H-S lower bound and close to Ω 1 with the H-S upper bound. Then, the gradient of the compliance coincides with the topological derivative in Ω 0 and Ω 1 .
To see this result, from Equation (9), we can write the gradient of the compliance (when using the isotropic lower and upper bound as interpolation function) as follows: Then, we can identify coefficients in (39) and (40) with precisely the ones established by the topological derivative in (29) and (30), and thus, Finally, using (42) in Equations (25), (27), and (9), we obtain Remark 6. This result confirms that the material in the neighborhood of the infinitesimal inclusion inserted in the perturbed domain used to compute the topological derivative can be understood as an homogenized microstructure with circular inclusions (H-S bounds microstructures) with a fraction volume value equivalent to the infinitesimal inclusion volume. This is to say that, in order to study a change on the compliance, the topological derivative concepts proposes in fact to insert H-S microstructures.

RATIONAL FUNCTIONS
Definition 3. Let us introduce the following particular family of rational functions: where A, B, C, and D ≠ −1 are scalar parameters to be determined by fixing the value and its derivative in the extremes, ie, Remark 7. Note that function ∈  is equivalent to the following expression: This result is easily obtained by computing the derivative of the rational function as and imposing conditions (45). Then, we can solve the system of four equations to find the value of the parameters A, B, C, and D. This is Remark 8. Note that, with the definitions above, the following property holds: Remark 9. The generalized H-S bounds and therefore the upper and lower bounds of the shear and bulk modulus are rational functions. This can be easily seen by imposing in the general rational Equation (46) and obtaining which is precisely the generalized H-S equation (33).

SIMP-ALL METHOD
The idea of the SIMP-ALL method is precisely to propose an interpolation such that the gradient of the compliance coincides with the topological derivative in both Ω 0 and Ω 1 . As already shown in Equation (43), in the case of isotropic materials, this is equivalent to impose the derivative of the H-S lower bound in Ω 0 and the derivative of the H-S upper bound in Ω 1 for the shear and bulk modulus interpolations. 1 . This is Remark 10. Note that the shear and bulk SIMP-ALL interpolation function SA ( ) and SA ( ) can be rewritten in the following compact form: where 0 = 0 0 2 0 + 0 , 1 = 1 1 2 1 + 1 , 0 = 0 , and 1 = 1 were previously defined in (34).
These expressions can be easily obtained by computing the rational function parameters d 01 , d 0 , and d 1 defined in Equation (46). In the case of the SIMP-ALL interpolation, they can be written as These parameters let us build the SIMP-ALL general interpolation function f SA ( ) as where the last expression is obtained multiplying the upper and lower parts of the fraction by Finally, we have to take the corresponding f 0 , f 1 , 0 , and 1 parameters values for the case of the shear and bulk modulus.
Difference of two rational functions. Let A ∈  and B ∈  be two rational functions such that ′ To obtain this result, we will show in four steps that the difference of f A and f B is positive ∀ ∈ [0, 1].
(i) Computing the difference function. Let us compute the function f C as the difference of the two rational functions f A and f B . This is where the numerator n C ( ) is written as and the denominator d C ( ) as (ii) Computing the difference function numerator. For sake of compactness, we write the numerator in matrix form as where c 0 = N 11 + N 22 + N 31 and c 1 = N 12 + N 22 + N 31 . Note that, from Remark 8, N 11 , N 12 , N 22 , and N 31 can be written as and the sum of N 22 and N 31 as N 22 . Thus, we can write coefficients c 0 and c 1 in a simplified form as (iii) Computing the difference function denominator. From Equation (58), this is simply (iv) Positiveness of the difference for all possible scenarios. The idea is now to analyze if the numerator and the denominator of the difference have always the same sign in all possible cases. First, note that, since A ∈  and B ∈ , their image should not be unbounded (should be between [f 0 , f 1 ]). Thus, from (62), we see that the following two scenarios are only possible: Then, the assumptions of ′ B (0) ≤ ′ A (0) and ′ A (1) ≤ ′ B (1) entail the following inequalities: In view of (64), we summarize all the possible scenarios in the following table.
In the table, Cases 1, 2, and 3 are the only possible cases that makes f A and f B satisfy condition (63), and consequently, the rest of cases are not included. Thus, we can guarantee positiveness of function f C reasoning as follows: Remark 11. Based on the results shown above, although it may be obvious, we could show that the H-S upper bound is above the H-S lower bound. Indeed, where we have used f 1 , f 0 ≥ 0 for coercivity reasons and 1 ≥ 0, 0 ≥ 0 is fulfilled by definition from table (34) for both the shear and the bulk modulus. Thus, as expected, the H-S bound for both the shear and bulk modulus is always above the H-S lower bound. This is

SIMP-ALL in between H-S bounds.
At this point, we can show that the SIMP-ALL interpolation for the bulk and shear modulus rest always between the H-S bounds, this is This result is certainly the second main contribution of this work. To show it, we only have to verify the conditions stated before Equation (55). In fact, the H-S upper bound and the SIMP-ALL interpolation function fulfills by construction Condition ′ UB (0) − ′ LB (0) ≥ 0 has been already shown in (67). We can proceed similarly for the H-S lower bound. By construction, and condition ′ SA (1) − ′ UB (1) = ′ LB (1) − ′ UB (1) ≥ 0 has been also already proven in (67).

COMPARISON BETWEEN SIMP AND SIMP-ALL
The celebrated SIMP interpolation function 19 proposes basically to combine the property of both materials as follows: where p is the penalization parameter, frequently taken as p = 3.

Interpolation function comparison
We seek in this section to compare both interpolations in the following two different situations.

(a) Material-void interpolation
In regularized topology optimization, the material of the domain (base material) is usually interpolated with a void material, modeled with an extremely weak material. As discussed in the work of Bendsøe and Sigmund, 2 standard SIMP proposes an isotropic interpolation, which its corresponding intermediate values (gray) cannot be interpreted in some situations as a composite material composed by the base material and the void material (the interpolation function remains outside the H-S bounds). To circumvent this limitation, an adaptive value of the penalization parameter (in terms of the Poisson ratio value of the base material 1 ) is proposed in the aforementioned work 2 for two-dimensional plane stress cases by changing the penalization parameter as follows: .
With this adaptation of the penalization parameter, the SIMP method remains in between the H-S bounds when 0 → 0 and 0 → 0. In the following, we will call this interpolation choice as adaptive SIMP. For comparing SIMP, SIMP-ALL, and adaptive SIMP, in this work, we have proposed a sequence of composite materials composed by an isotropic base material with properties (E 1 = 1, 1 ∈ {1∕3, 0, −0.5, −0.75, −0.9}) and a void material with properties (E 0 = 10 −3 , 0 = 1∕3). Note that the value of 0 when E 0 → 0 has no relevance. The SIMP, adaptive SIMP, and SIMP-ALL interpolation function have been depicted in Figure 2 in relation with the shear and bulk H-S bounds. In view of the results, when the Poisson ratio is = 1∕3, the three interpolations behave similarly. In fact, SIMP adaptive and SIMP coincides since p * (1∕3) = 3. Additionally, both the SIMP and the . The former coincidence is proved in the work of Bendsøe and Sigmund 19 while the latter is obtained by construction. In this case, both interpolation remain in between the H-S bounds. However, as soon as 1 decreases, SIMP, SIMP adaptive, and SIMP-ALL start differing. In fact, SIMP-ALL and SIMP adaptive interpolation function, in contrast with SIMP, remain in between H-S bounds as expected. Additionally, note that as soon as 1 decreases, eg, 1 = −0.9, p * = 40 (Case (E)), the adaptive SIMP interpolation, in comparison with the SIMP-ALL interpolation, increases extremely markedly resulting in possible limitations when solving topology optimization problems.

(b) Bimaterial interpolation
For a more general case, when the weak material is not necessary much less stiff than the base material, the topology optimization problem is in fact a bimaterial optimization problem. The aim is to decide in which part of Ω is appropriate to use one material or the other, or in case of "grays," a combination of both. Two different cases have been presented. The first one (Case (A)) seeks to model an Steel-Aluminum composite. In this sense, the material properties have been taken as E 0 = 1∕3, 0 = 0.35, E 1 = 1, and 1 = 0.3. The second example (Case (B)) intends to model an extreme case where one material has larger bulk modulus value and the other material has larger shear modulus value. To this aim, the material properties considered are E 0 = 0.9, 0 = −0.9, E 1 = 1, and 1 = 0.3. In the bimaterial problem, the value of the penalization parameter p for the SIMP interpolation function is unclear. In the following, we take a penalization parameter of value p = 3. The SIMP adaptive method is even not defined. The corresponding H-S bounds, SIMP, and SIMP-ALL method have been depicted in Figure 3. As expected, in both cases (A and B), SIMP function falls outside H-S bounds. However, although the space between H-S bounds is significantly narrow, the SIMP-ALL interpolation function comes to lie inside them, as stated in Equation (70). Thus, the use of SIMP-ALL interpolation seems as convenient as the SIMP interpolation for 1 ≃ 1∕3 and more convenient (from the physical point of view) for different values of 1 and when solving bimaterial problems.

Numerical examples comparison
In order to see the numerical effects of using different kind of interpolation functions, a cantilever beam example has been solved for the sequence of Cases described in Figure 2. The size of the domain is 2 × 1 with the left side clamped. A distributed vertical force of unit value is applied in a 0.2 centered part of the right side. The domain is equipped with a structured triangular mesh of 20 000 elements from a Cartesian grid by splitting each cell into four triangles. The density function is approximated with P 1 Lagrangian finite element functions 21 and the following filter is applied to the density function: where k is the density value in element k and S k is the set of nodes of element k, see the work of Amstutz et al 22 for further details. This filter is defined in a discrete sense, in terms of the shape functions. Other filters depending on the distance are also possible. We consider the compliance as the cost function and the volume (volume target of the stiff material V = 0.5) as a constraint of the topology optimization problem. Additionally, we consider the MMA algorithm, 23 widely used in the community, to solve the problem with a stopping criteria tolerance of 10 −4 and a maximum number of iterations of n max = 5000. The design variable is initialized with stiff material overall the domain ( = 1).
As depicted in Figure 4, we observe similar optimal topologies regardless the interpolation used when the Poisson ratio of the stiff material is closed to 1 = 1∕3. Additionally, the number of iterations and the finite-element solver evaluation are also similar. The examples are computed in MATLAB© with a standard PC (3.40 GHz processor in a 64-bit architecture). However, some differences appear when 1 becomes smaller. The optimal topologies when using the SIMP interpolation require large number of iterations and finite-element solver evaluations to converge and gray areas appear. This is probably due to a small penalization value appeared in that cases (see Figure 2). In the case of SIMP adaptive, in which the penalization value is rather large, although it seems to be almost free of numerical instabilities, nonstandard topologies appear. Likewise, a significant computational effort in terms of iterations is required to converge. In contrast, the SIMP-ALL interpolation behaves as adequately in terms of optimal topology and computational effort as  1∕3), all interpolations behave similarly. As 1 becomes more negative, gray areas appear when using SIMP and a larger number of convergence iterations are needed. Adaptive SIMP is presents no large gray areas but a noncompetitive cost function and a larger number of convergence iterations appear. SIMP-ALL interpolation presents no gray areas, the number of iterations remains small, and the objective function achieves considerable smaller values. In comparison with SIMP, the SIMP-ALL method evidences an auto-penalizing effect in a wider range of cases [Colour figure can be viewed at wileyonlinelibrary.com] when considering standard Poisson ratio values, typically = 1∕3. In Case (C), the SIMP method did not converge. Note that, comparing cost values between interpolation could be fuzzy, since in gray areas, each interpolation proposes different stiffness values. In some cases, SIMP interpolation may propose values stiffer than physically possible (out of H-S bounds). To have a stiffness comparison with physical sense, the cost of the optimal topologies obtained by the three interpolations have been evaluated using the SIMP-ALL interpolation (J SA ). The value of the cost (compliance) shown in Figure 4 is normalized with the value of the cost in the initial guess. Admittedly, using the SIMP-ALL interpolation, a decrease in the cost in Case (E) of 35% is obtained when comparing with the SIMP interpolation and of 30% when comparing with the SIMP adaptive interpolation. This result shows the convenience of using the SIMP-ALL interpolation when extreme isotropic materials (see, for example, the work of Greaves et al 24 ) or bimaterial problems are considered.
Finally, the gray areas appearing in the SIMP method is probably because intermediate values with this interpolation are competitive, as seen in Figure 2. Large stiffness values (more than physically possible) are provided for intermediate values. In this case, SIMP behaves not as an auto-penalizing interpolation function. In contrast, SIMP-ALL seems to be auto-penalizing in all the cases. This is not thoroughly surprising since, as shown in Equation (43), the topological derivative (used in the SIMP-ALL interpolation) proposes in fact to insert H-S microstructures near the extreme values. Thus, the possibility of adding a certain amount of density (or inserting certain stiffer inclusions) is not competitive since they provide small stiffness increment. Intuitively, microstructure made of weak material with a small number of circular inclusion transmits barely the load, making it, in general, not attractive from the optimization point of view. This may explain why the proposed SIMP-ALL interpolation, in contrast with SIMP, evidences in the results of Figure 4 an auto-penalizing behavior for all the studied cases.

CONCLUSIONS
The main contributions of this work can be summarized in the following four points.
First, we have showed in Equation (43) that the topological derivative of the compliance is precisely the linear approximation of the H-S bound (lower bound for the weak and upper bound for stiff material). In other words, using the topological derivative for the compliance is equivalent to inserting H-S bound microstructures (circular inclusions). This result strongly relates both concepts and allow us to give a new interpretation of the topological derivative.
Second, we have proposed the SIMP-ALL method, a new interpolation function for the shear and bulk modulus. It is defined as a rational function and build so that its derivative coincides with the H-S lower bound derivative in the weak material and the H-S upper bound derivative in the stiff material, making the gradient coincide, as shown in (43), with the topological derivative in the stiff and the weak material. The key point of the SIMP-ALL interpolation is that in Equation (70), we have shown that it stays always in between the H-S bound regardless the materials to be interpolated, allowing us to interpret its intermediate values always as real microstructures.
Third, we have compared the SIMP-ALL with the SIMP and SIMP adaptive interpolation function. We have observed similar behavior when 1 is closed to 1∕3 and quite different behavior as it becomes smaller. In contrast with SIMP, the SIMP-ALL interpolation has shown to stay in between the H-S bounds, numerically verifying the result of Equation (70), not only in the case of weak-stiff material but also in the more challenging case of interpolating two general isotropic materials.
Finally, concerning numerical examples, we observed also similar behavior of the optimal topologies around = 1∕3. However, as it becomes smaller, the SIMP-ALL method has shown to obtain up to 20% stiffer topologies, free of gray areas and smaller computational cost. In comparison with SIMP, this result shows that the SIMP-ALL method is auto-penalizing in a wider range of cases.
From the author's point of view, rather than an alternative, the SIMP-ALL interpolation method can be understood as a generalization of the SIMP method, adequate for all kind of interpolation materials, with physical meaning, free of heuristic parameters, and with no extra computational effort. In this sense, SIMP-ALL interpolation bridge also the gap between the topological derivative community and the SIMP community, allowing the first to use regularized algorithms and the second to use the topological derivative concept. A MATLAB implementation of the SIMP-ALL interpolation is available at https://github.com/FerrerFerreAlex/SimpALL and has been also provided in the Appendix.

SIMP-ALL MATLAB CODE
In this Appendix, we provide a simplified MATLAB code to facilitate the use of the SIMP-ALL interpolation to the reader. The code is also available at https://github.com/FerrerFerreAlex/SimpALL. The SIMP_all_interpolation.m function is devised to be used when starting the topology optimization code. The time of computing the interpolation function is negligible making it as usable as the well-established interpolation functions. The constitutive interpolated tensor and its derivative are provided since these are usually used when computing the compliance function and the gradient of the compliance. To obtain the interpolation function, we first compute the shear and bulk modulus and the additional parameters. Then, the coefficients and the rational function are computed. Finally, we differentiate with respect to the density variable to obtain the derivative.