A stress-based formulation of the free material design problem with the trace constraint and multiple load conditions

The paper deals with minimization of the weighted sum of compliances related to the load cases applied non-simultaneously. The design variables are all components of the Hooke tensor, subject to the isoperimetric condition bounding the integral of the sum of the Kelvin moduli. This free material design problem is reduced to an equilibrium problem – in two formulations – of an effective body with locking. The stress-based formulation reduces to minimization of an integral of a certain norm of stress fields over the stress fields which equilibrate the given loads. The equivalent displacement-based formulation involves a locking locus defined by using a norm being dual to the previous one. The optimal Hooke tensor is determined by using the stress fields solving the auxiliary locking problem. To make the optimal Hooke tensor positive definite one should consider at least 3 load conditions in the 2D case and not less than 6 load conditions in the 3D case.


Introduction
The topology optimization of elastic structures provides the methods of the optimal layout: of two or several materials (1), of the bars of pin-jointed frameworks (2), of the elastic moduli within a given feasible domain (3), see Bendsøe (1995), Bendsøe and Sigmund (2003). In the case of a single load, a standard merit function is the compliance of the whole structure. The theory of this problem is the subject of the extensive studies, see Lipton (1994), Allaire (2002), Cherkaev (2000) and Haslinger et al. (2010) and the references therein. In the multiple load case the optimum design should be considered within the Pareto framework of the vector optimization. The Pareto solution can be found on the basis of the solution to the optimum design problem with the merit function being the weighted sum of the merit functions for subsequent load cases, cf. Hillermeier (2001) and Jahn (2004). To construct the Pareto solution one should solve the family of minimization problems over all possible weighting factors η i , the sum of which being 1. Therefore, the present paper is aimed at minimizing the weighted sum of the compliances corresponding to the given loads applied non-simultaneously, for fixed values of the weighting factors η i .
Although the topology optimization comprises other problems than compliance minimization, like optimal choice of the first eigenfrequency, the compliance minimization is the first on the list as the only problem being self-adjoint: the adjoint problem coincides with the initial one.
In problems named above by (1) the amounts of the materials are given thus providing the isoperimetric conditions expressed by the integrals of the mass densities. The problem is to place the materials to minimize the merit function (e.g. the compliance); the problem is rationally posed due to the isoperimetric conditions imposed, which prevent the solutions from being trivial.
In problems (2), usually the total mass (or weight) of the structure is bounded. The mass is determined by the areas of the cross sections and member lengths. Alternatively, the stresses in bars can be bounded. Such conditions reduce the problems (2) to trusses, since the bent bars are eliminated as stressed non-uniformly, and hence stressed non-optimally.
Yet it is much less obvious how to augment the problems (3) by isoperimetric conditions. The design variables are elastic moduli C ij kl . Consequently, the isoperimetric conditions should be expressed directly by these variables. Since there is no link between the moduli and the mass density, the conventional isoperimetric conditions concerning the mass density distribution cannot be applied. Probably, just this unclear question was the reason why the first papers on the problem of optimal choice of elastic characteristics have appeared only in 1990's, although the underlying theory of the compliance minimization has been elaborated in 1980's, while the theorems on invariants of the 4th order tensors were available in 1970's and found the mature state in 1980's, see Betten (1986Betten ( , 1987a. The first paper on the optimal choice of the elastic moduli in which the isoperimetric condition has been expressed in terms of invariants of the Hooke tensor is due to Ringertz (1993), yet only the paper by Bendsøe et al. (1994) put forward the crucial formulae and cast the problem in the clear mathematical setting. The extension of the previous paper to the multiload case is due to Bendsøe et al. (1995). This setting is called there Free Material Design (FMD) and this term will be used in the present paper. The merit function is chosen as a weighted sum of the compliances. The isoperimetric condition is expressed by the integral over the fixed design domain of a scalar function φ(C). Two cases have been discussed: Here trC is the trace of C while || · || is the Frobenius norm of C. The scalar function should involve the invariants of C. The functions mentioned above are admissible invariants, yet other choices are also worth considering. One can make use of the theory of invariants of 4th order tensors outlined in Betten (1986Betten ( , 1987a, Zheng (1994), Jemioło and Telega (1997). The invariance property of the scalar φ(C) is not sufficient to make the FMD problem well posed. Recently Barbarosie and Lopes (2008) and Haslinger et al. (2010) have formulated the sufficient conditions for φ for making the FMD problem solvable. The scalar functions φ chosen in Bendsøe et al. (1994Bendsøe et al. ( , 1995 satisfy these conditions. One can note that the main effort of the hitherto existing works on FMD has been put on the numerical procedures to solve the optimization problem. These algorithms are usually based upon the displacement formulation which leads to a saddle point problem, see e.g. (7) in Zowe et. al (1997) and (2.8) in Haslinger et al. (2010). The minimization operation over admissible Hooke tensors is preceded by the maximization operation over the kinematically admissible displacements. Further modification of this formulation is possible; the design variables can be eliminated to arrive at the dual formulation, see (13) in Kočvara et al. (2008) and Werner (2000). The difficulty of this dual formulation lies in the nonlinearity of the bounding conditions. Moreover, to find the optimal moduli C ij kl one should solve an auxiliary problem, see p. 87 in Kočvara et al. (2008). The contemporary codes to solve the non-linear semi-definite programming problems (SDP) make it possible to find effectively the solution to the FMD problem in its dual formulation. Yet this procedure is complex, needs the most advanced numerical tools and is by far not direct, see Section 2.4 ibidem. Within this framework one can take into account the state constraints concerning the displacements and stresses, see Stingl (2007, 2012) and Section 2.6 in Haslinger et al. (2010).
The papers referred to above do not enter into the algebraic structure of the Hooke tensors. The specific properties of symmetry of these tensors determine their spectral representation. The eigenvalues of the Hooke tensor of the known materials are positive; they are called Kelvin moduli and are denoted by λ 1 , λ 2 , . . . , λ m ; m = d(d + 1)/2, d being the dimension of the problem setting. Therefore, in the planar problem the number of the Kelvin moduli is 3, while in the 3D setting it is equal 6. The eigentensors are diads of the form ω ⊗ ω, where ω is a symmetric tensor of order 2. The theorems on this spectral representation are provided in the articles by Rychlewski (1984), Blinowski et al. (1996), andSutcliffe (1992), Norris (2005Norris ( , 2006, Moakher (2008), Mehrabadi and Cowin (1990), Jemioło and Telega (1997). This representation makes it possible to find a new interpretation of the results by Bendsøe et al. (1994), cf. Turteltaub and Washabaugh (1999). The spectral representation has also been used in the paper by Banichuk (1996) on the optimal local orientation of a given anisotropic material. In the FMD approach the condition λ K > 0 is relaxed to λ K 0 to achieve the best material properties.
The history of the development of the two-material layout theory shows that the most convenient approach to the compliance minimization is to express the compliance as the minimal value of the complementary energy over the statically admissible stress fields, see Lurie and Cherkaev (1986), cf. (2.4) in Allaire and Kohn (1993). The layout problem is thus reduced to a static problem of a body of nonlinear elastic characteristics given by a new, effective density of complementary energy, see (2.7) in Allaire and Kohn (1993). This final stress-based formulation of the minimum compliance problem does not involve the design variables. The optimal values of the design variables are expressed by the stress fields which make the effective complementary energy minimal.
The result mentioned above has been an inspiration for a similar reformulation of the FMD problem. Indeed, the FMD problem for the single load case can be reduced to a similar two-level scheme. The first step is to solve an auxiliary minimization problem over the statically admissible stress fields. Having found the minimizer one can determine all the optimal elastic moduli, see Czarnecki and Lewiński (2012).
The present paper is aimed at extension of the latter result to the case of arbitrary number of independent loads. The minimization over the tensors ω will be performed first thus leading to the explicit formulation serving as the point of departure for the subsequent step: the minimization operation over the Kelvin moduli obeying the isoperimetric condition involving the sum of these moduli as the integrand. This condition is equivalent to the isoperimetric condition with the integrand φ(C) = tr C, assumed in most papers on the FMD problems. The elimination of the Kelvin moduli reduces the FMD problem to the minimization problem of a functional depending on independent stress fields. Its integrand is the sum of singular values of the matrix composed of the components of these stress fields. The derivation of this result is one of the main objective of the present paper. Since this integrand has a linear growth with respect to its argument, the relevant minimization problem can be interpreted as an equilibrium problem of an effective body of locking properties, seeČyras (1972), Demengel and Suquet (1986), Telega and Jemioło (1998). To have a complete insight into the deformation of the body made of a locking material one should solve two mutually dual problems: the stress-based problem (minimization over statically admissible stresses, or rather stress rates) and the kinematic problem (with maximization over kinematically admissible displacements). In the considered problem of n loading conditions we have to deal with n stress fields and n strains, the latter lying within a locking locus at each point of the design domain. This result is crucial as linking the locking material theory with the FMD optimization problem. Yet this link is not a surprise. The known Michell truss problem is also expressed by two mutually dual formulations, both of them lying within the framework of mechanics of materials with locking, see Rozvany (1976), Strang and Kohn (1983), Lewiński (2004). Moreover, similar problems appear in the variable sheet problem, see Czarnecki and Lewiński (2013).
The present paper puts forward a new numerical method to solve the planar FMD problem for two load cases in its stress-based setting. The method is an extension of the numerical approach proposed in Czarnecki and Lewiński (2012) for the single load case. Although the meshing of the design domain is a starting point, the numerical method proposed does not belong to the finite element (FE) methods, since the integrand of the functional has a linear growth, while most FEM codes are applicable for the integrands of quadratic growth. Consequently, the stresses are not linked with strains by constitutive equations, but by the Kuhn-Tucker conditions. The kinematic formulation does not lead to the set of equations with a quadratic stiffness matrix. No aggregation procedure is applied to construct the global system of the governing equations. Instead, we have to solve the stress based problem by a direct minimization of the prescribed integral. The trial stress fields are interpolated element-wise by polynomials (which is the only common feature with FEM) and subject to the set of equilibrium equations implied by the variational equations of equilibrium. The set of equilibrium equations results in the algebraic under-determinate system of linear equations. The solution of this system is represented numerically with using the singular value decomposition (SVD). The redundant parameters are then determined by minimization of the prescribed functional. It is remarkable that these parameters are not subject to any constraints, hence the available, highly efficient optimization codes specially developed for the unconstrained optimization problems can be applied. The computational problem considered has never been discussed before. The numerical results can only be compared with the available numerical solutions to the FMD problem posed different way, as a saddle point problem. The final layouts found in the present paper compare favorably with those published previously by Werner (2000) and Kočvara et al. (2008), thus confirming the correctness of the stress-based method proposed.
The optimal layouts published in Czarnecki and Lewiński (2011) refer to the 2D case and to the arbitrary choice of tensors ω and given distribution of the Kelvin moduli. The layouts shown in the present paper refer to variation of both: tensors ω and Kelvin moduli. The results published in the present paper prove that admitting the Kelvin moduli to vary changes the final design substantially and contributes to a visible decrease of the compliance.
In the notation of the minimization problems the following convention is adopted if the structure of the condition x ∈ X is complicated. Moreover, the following standard notation is used. The domain of the body whose material properties are designed is denoted by Ω.  (v , . . . , v d ) in Ω determines the virtual strain field ε(v); ε ij (v) = (v i,j +v j,i )/2; (·) i = ∂(·)/∂x i . Locally ε ∈ E 2 s , E 2 s being the set of symmetric tensors of second rank. The scalar product of ε, κ ∈ E 2 s is defined by ε · κ = ε ij κ ij ; i, j run over 1, . . . , d. The Euclidean norm of ε ∈ E 2 s is defined by ||ε|| = (ε · ε) 1/2 . Tensors of E 2 s can be viewed as vectors in R m , m = d(d + 1)/2, according to the rules: This correspondence is explained in Moakher (2008). Thus elements of E 2 s will be treated as tensors or as vectors, depending upon the context. The p-th eigenvalue of a square matrix B is denoted by μ p (B). The p-th singular value of an arbitrary matrix A is defined by The indices i, j, k, l run over 1, . . . , d; d being the dimension of the problem. The indices K, L run over 1, . . . , m.
The indices α, β take values 1, . . . , n, n being the number of the load conditions. The set of matrices of m rows and n columns is denoted by M m×n .

The FMD problem in its stress-based formulation. Case of n loading conditions
The aim of the present section is to formulate the FMD problem for n independent loading conditions. For simplicity, the body forces will be omitted. The FMD means designing characteristics of an elastic body. We assume that body occupies a fixed Ω in R d . Thus its geometry is not subject to any change. Its segment Γ 2 of the boundary Γ where displacements are kept zero is also fixed. The remaining part Γ 1 of the boundary is subject to given tractions T α , α ∈ {1, . . . , n}. This load will not be subjected to any change.
Introduce the linear form representing the virtual work of T α on the trial displacement field v ∈ V (Ω); V (Ω) represents the space of kinematically admissible virtual displacements; v = (v 1 , . . . , v d ).
Consider virtual stress fields τ = (τ ij ) in Ω of class L(Ω, E 2 s ), which means that locally τ (x), x ∈ Ω is an element of the set E 2 s . The fields τ satisfy the regularity assumptions following from the equilibrium equation: Such virtual stress fields are called statically admissible stresses. They constitute a linear affine set Σ α (Ω). The (6) is called a virtual work equation. The design variables of the FMD problem are the elastic moduli C ij kl for a.e. x ∈ Ω. These components are referred to the basis e i ⊗ e j ⊗ e k ⊗ e l ; they represent a tensor C(x). The notation C ∈ E 4 s comprises the symmetry conditions: C ij kl = C klij ; C ij kl = C jikl . Any tensor C ∈ E 4 s admits the following spectral decomposition, see Rychlewski (1984), Sutcliffe (1992) where λ 1 λ 2 . . . λ m 0 are called Kelvin moduli and ω K (x) ∈ E 2 s are called eigenstates, see Rychlewski (1984). They satisfy the orthogonality conditions: If λ K = 0 for K = m 1 + 1, . . . , m, then C −1 (x) will be interpreted as where m 1 = rank(C). Note that then CC −1 will not be equal to a unit tensor in E 4 s . Assume the triplet (u α , σ α , ε(u α )) such that u α ∈ V (Ω), σ α and ε(u α ) are linked by σ α (x) = C(x)ε(u α (x)) and let σ α ∈ Σ α (Ω); α is fixed. The triplet constitutes a solution to the elasticity problem corresponding to the T α load. The main condition of uniqueness of the solution is λ m > 0. The existence problem is discussed in Duvaut and Lions (1976).
The quantitỹ is called the compliance. In case of λ m > 0 this quantity can be equivalently expressed as below The equalityΥ α = Υ α reflects the celebrated theorem of Castigliano, see Duvaut and Lions (1976). In the present paper the notion of compliance is extended to the case of rank(C) < m by using the formula (11) with C −1 given by (9). Let η α ∈ [0, 1] be fixed numbers such that η 1 + . . . + η n = 1. In the FMD problem considered all components of the tensor field C(x) for a.e. x in Ω are design variables. This will be expressed by writing Υ α = Υ α (C).
Introduce the functional to be minimized A detailed discussion of solutions to the problem of minimization of F η for η α linked by η 1 + . . . + η n = 1 makes it possible to find the whole Pareto front of the solutions. This question will be explained in Section 6.
The trace of C is defined by C ij ij or Since λ K 0 we can write tr C = ||λ|| 1 , where λ = (λ 1 , . . . , λ m ) and || · || p is a p-norm. Let · be averaging over Ω and let E 0 be a referential elastic modulus. Let H (Ω) be the set of appropriately regular tensor fields C on Ω such that for a.e. x in Ω, C(x) ∈ E 4 s and its Kelvin moduli are non-negative. We say that C ∈Ĥ (Ω) if C ∈ H (Ω) and For fixed η 1 , . . . , η n satisfying η 1 + . . . + η n = 1 the FMD problem considered is assumed in the form Let us substute (11), (12) into (16) and switch minimization operations. We find wherẽ It turns out that the functionalJ η, E 0 can be explicitly constructed by performing: first minimization over ω K and then over λ K . This construction is shown in the next section.
3 Reduction of the FMD problem to a locking material setting

Case of arbitrary, finite number of load conditions
Let us introduce the representation (9) of C −1 into (18). Let us express the integrand ψ of the functional minimized in (18) as follows The argument x is suppressed for simplicity. The matrix S η is constructed of the column vectors √ η α τ α using the representation (2) Let us define the mapping treating τ α as column vectors m×1, using the representation (2), (3). Hence The mapping S transforms the set of n column vectors in R m into a matrix of m × n dimensions. Definẽ Note that We know from algebra that rankS η = rankŜ η = rankS η . In case of n = m the matrixS η is a Gram matrix. This special case will be commented later.
Let us define Assume that and recall that Let us invoke the ,,rearrangement inequality", cf. Hardy et al. (1999) We see that to minimize the quantity ψ, see (19), over the fields ω K we should, for a.e. x ∈ Ω, arrange the quantities a K according to (25). The biggest value of a K , or a 1 , is attained by This maximum is attained by ω * 1 (x) which is the eigenvector corresponding to the first eigenvalue ofŜ η computed at x. Thus, for fixed x, the pair (a 1 , ω * 1 ) is the first solution to the eigenvalue problem: find (λ, y) such that We shall write a 1 = μ 1 (Ŝ η ), according to the notation convention assumed in Section 1. The next pairs (a K , ω * K ) are solutions to the same eigenvalue problem; a K , K = 1, . . . , m 1 are ordered as in (25). Note that for given x the number m 1 is fixed. In fact m 1 depends on x.
Concluding, minimization in (18) over the fields ω K satisfying the conditions: ω K · ω L = δ KL pointwise rearranges (17) to a new minimization problem where λ = (λ 1 , . . . , λ m 1 ) and The integrand of (31) reads where whereŜ η is given by (20). Let us substitute (32), (33) into (30). We find Minimization over λ K can be done analytically, by extending the results of the Appendix of Czarnecki and Lewiński (2013). The final result is Note that Moreover, note that the summation over K in (36) can be extended to K = m, since s K (S η ) = 0 for K = m 1 + 1, . . . , m. Instead of (36) it is more clear to write since now the upper limit of K is x independent. The optimal λ K = λ * K are expressed by minimizers τ * α of problem (38) We note that λ * . The optimal eigenstate ω * K is the K-th eigenvector of problem (29). For given x ∈ Ω the optimal Hooke tensor is expressed by Let us write and note that Thus the integral of (38) is a homogeneous function of degree 1. Instead of solving (38) it is sometimes easier to solve the problem dual to (38). In the following we shall show the passage from (38) to its dual. This passage is inspired by the work of Strang and Kohn (1983) concerning Michell trusses.
Let us disclose the equilibrium constraints in (38) and rearrange this problem to a saddle point problem: Now we interchange the min and max operations. The nested problem reads where ε α = ε(v α ). By analogy with Section 3.1. in Czarnecki and Lewiński (2013) we compute and The result (45) reduces the problem (43) to the form Let us now re-write (38) in the form The problems (48), (49) are dual to each other. They can be interpreted as two problems reflecting mechanical behavior of a certain n-components mixture of locking properties. The behavior of stresses is governed by (49), while strains are governed by (48). The relations between stresses and strains are not expressed by equations but by Kuhn-Tucker conditions. The noted link with the theory of materials with locking is not a surprise, since the theory of Michell trusses can be treated as a theory of a locking material, see Rozvany's (1976) comments and the footnote on Zenon Mróz suggestion, see page 64 in this book. The present paper does not deal with choosing spaces of admissible functions in the problems (48), (49). It seems that the stresses in (49) should satisfy the conditions mentioned in Strang and Kohn (1983) and Demengel and Suquet (1986), ; M 1 being the space of bounded measures on Ω. Yet the problem remains how to relax the boundary conditions on Γ 1 . The known methods of F. Demengel cannot be applied here directly.
Assume that both the problems (48), (49) have been solved and the tensor (40) has been constructed. The problem arises whether the elasticity problems corresponding to the T α loads and the anisotropy given by (40) are solvable, taking into account that, in general, C * is not positive definite. The affirmative answer is delivered by Th. 1 in Zowe et al. (1997), based on the displacement approach. Other existence theorems are formulated in Haslinger et al. (2010).
The results (50)-(53) have been for the first time published in Czarnecki and Lewiński (2012). These results are compatible with the results published in Bendsøe et al. (1994) and in the next papers on FMD. Let us stress that the FMD version developed in Bendsøe et al. (1994), Zowe et al. (1997), Hörnlein et al. (2001) is a displacement -based version and these papers neither include the formulation (50) nor the formulation dual to it, to be shown later.
According to Th. 6 in Zowe et al. (1997) there exist stresses σ ij and strains ε ij satisfying all equations of the elasticity problem for a body with anisotropy given by (53). It is worth noting that the stress σ is collinear with the auxiliary stress field τ * solving (50). One can derive: The problem dual to (50) is a special case of (48). It reads The locking locus is a unit ball: since ||ε|| * = ||ε||.
Step 3 Solve the eigenvalue problem for the matrixŜ η (x) to find the eigentensors ω * In 2D case a direct construction of ω * 1 , ω * 2 is as follows: Compute Find an angle α(x) such that α ∈ 0, π 2 and for selected x ∈ Ω. The optimal eigenstates are determined by Step 4 Compute the components of the optimal Hooke tensor by Step 5 The optimal compliance is given by Step 1 exceeds the framework of the conventional finite element method, since the integrand in (63) is of linear growth. Consequently, no constitutive equations are associated with this problem and in the discretized setting no stiffness matrix arises. To solve the problem (63) one should develop a new numerical scheme, which will be constructed in Section 4. Note, that for η = 0 and η = 1 the problem (63) corresponds to a single load case, see problem (50).
The problem dual to (63) is a special case of (48). It has the form where the locking locus is the ball: defined by the norm with |||(σ , τ )||| given by (62).

Three and more load conditions
In case of n = 3 and d = 2 the matrix S η , see (63), is a square 3 × 3 matrix; here n = m = 3. In case of n = 6 and d = 3 the matrix S η is a square 6 × 6 matrix. In these two cases all optimal Kelvin moduli attain positive values and the optimal tensor C * becomes non-singular. Let us conclude: a) d = 2. If tensor C * is nonsingular then n 3 b) d = 3. If tensor C * is nonsingular then n 6 Therefore, to make C * nonsingular, the number of the load conditions should be equal to (or bigger than) the number of stress (or strain) components.

Numerical treatment of the FMD problem in its stress-based setting. The case of two load condition
The stress-based FMD problem (63) for two load cases, with the trace constraint, or the problem: will be solved numerically for some selected 2D problems with using the newly developed scheme of construction of statically admissible stress fields defined element-wise, along with the optimizer solving the minimization problem. The numerical approach starts from division of the design domain Ω into 4-node, quadrilateral, isoparametric finite elements with bilinear shape functions interpolating the eight stress fields  (1 ± ξ)(1 ± ζ ). The conditions τ α ∈ Σ α (Ω) mean that the following variational equations are fulfilled for all the kinematically admissible test fields v = (v 1 , v 2 ). These test fields will be interpolated within an e-th element Ω e similarly as the stress field, or where p 2i+0 , p 2i+1 (i = 0, 1, 2, 3) are the unknown nodal virtual displacements v 1 , v 2 at the vertices 0, 1, 2, 3 respectively. Substitution of (74), (76) can be expressed as are the vectors that span the s-dimensional kernel of the matrix B u and are the arbitrary, fundamental solutions of the two sets of linear equations (78). The representations (79) can be performed with using e.g. the singular value decomposition algorithm (SVD) and just this method has been applied in the present paper. The short description of the SVD algorithm has been outlined in Czarnecki and Lewiński (2012). Let us mention that the QR decomposition of the matrix makes it possible to find the representation (79) more effectively than the SVD decomposition, due to a smaller computational complexity. Yet the Gaussian elimination as well as the LU decomposition methods fail to give satisfactory results for bigger systems of rectangular linear equations. Upon constructing the two solutions of the two linear, rectangular algebraic systems (78), we obtain in each e-th finite element Ω e the approximationsΣ α (Ω e ) of the two statically admissible sets of the stress fields Σ α (Ω e ), each determined by s global parameters  The numbering of indices j in α t * j and I j in b k I j (j = 0, 1, . . . , 11) is written in the local (element) notation and depends on the global allocation matrix.
Recall that the integrand of (73) is a norm in R 6 , hence is convex.
The problem (73) is now approximated by the following unconstrained minimization problem The partial derivatives of the integrand Π = |||( √ ητ 1 , √ 1 − ητ 2 ||| with respect to β i , i = 1, 2, . . . , S are computed by the rules The partial derivatives  (85) is now reformulated to the fully algebraic, finite dimensional and unconstrained problem of nonlinear mathematical programming where w = w(ς Q ) and det ∇x(ς Q ) are the weights and the determinants of the Jacobian matrix of the transformation (from the reference to the given finite element) calculated at the Gauss integration points ς Q = (ξ Q , ζ Q ) ∈ ω, respectively. We implement the well known assembly technique from the classical FEM: the formulae (87) for partial derivatives of the function = Π(ς , η, β 1 , β 2 , . . . , β S ) at arbitrary, but fixed point ς = (ξ, ζ ) ∈ ω make it possible to calculate all components of the gradient gradΠ = gradΠ(β 1 , . . . , β S ) of the objective function The algorithm Assume that coefficient η is given. The computational procedure consists of the following steps Step 1. Set the two rectangular systems of linear algebraic equations (77) or BΞ = α Θ, α = 1, 2, in accordance with FEM.
Step 2. Separate upper sub-matrix B u and two upper, subvectors α Θ u , α = 1, 2 corresponding to the unknown degrees of freedom of the FEM.

Case studies
The aim of this section is to show the optimal layouts of Kelvin moduli λ 1 , λ 2 and elastic moduli C ij kl within rectangular plates of various kinematic boundary conditions subject to two kinds of in-plane surface loads α T , α = 1, 2 that are non-simultaneously applied at Γ 1 and contribute to the functional (12) with weight factors η 1 = η and  Fig. 2). The finite element mesh is defined by n x × n y = 40 × 20 = 800 quadrilateral modules. The total number of nodes n o = (n x + 1)(n y + 1) = 861, which gives the total number of the columns (total number of the unknown nodal stress parameters) and rows (total number of the degrees of freedom) in matrix B equal to n s = 3n o = 2583 and m o = 2n o = 1722, respectively. The 9-and 3-points Hammer -Stroud rules of the Gauss integration for the two-dimensional element and one-dimensional segment are adopted, respectively.

The tractions
The values T max , w are assumed to be equal T max = 0.0376, w = 0.15, respectively, so all integrals parameter, see Figs. 4, 5, and 6. We consider three values of the weighting factor η: η = 0.9, η = 0.5 and η = 0.1. We note that the layouts depend heavily on the parameter η.
It is worth to note that the modulus λ 1 (with the layouts of "∨" or/and ,,∧" shapes) assumes extremal values in the zone near the applied forces and supports. The modulus λ 2 (of ,,|" shape) assumes extremal values only in the zone near the applied forces. The problem considered here has been previously discussed in Hörnlein et al. (2001). The result in Fig. 8 of Hörnlein et al. (2001) shows the layout of tr C (which equals to λ 1 + λ 2 ) for the case of η = 0.5. This result looks the same as the diagram of λ 1 shown in Fig. 5 because λ 2 is much smaller than λ 1 and that is why the plot of λ 1 + λ 2 has not been displayed. Let us note that neither the paper referred to above nor other hitherto published papers on FMD have presented the optimal distribution of the Kelvin moduli independently.
Example 2 Consider three plates supported at two non-sliding supports at their bottom corners. First plate (see Fig. 7a) is subject either to a horizontal load T 2 e 2 at its upper edge. Distribution of the values of the optimal Kelvin moduli λ 1 , λ 2 (λ 1 λ 2 ) and optimal components C ij kl referred to the (x 1 , x 2 ) coordinate system depend on the η parameter, see Figs. 8, 9 and 10. We consider three values of the weighting factor η: η = 0.9, η = 0.1 and η = 0.5. The third plate (Fig. 7c) is subject simultaneously to a load T e 1 + T e 2 at its bottom edge and to a load −T e 1 + T e 2 at its upper edge, where the function T = T (ζ ) has the form (92). For this single load condition, optimal layout of the one Kelvin modulusλ 1 and componentsC ij kl of the Hooke tensor are found on the basis of the algorithm presented in Czarnecki and Lewiński (2012) (see Section 3.2). If we e.g. assume that η = 0.5 for the two first plates, it may be interesting to compare the optimal results found for these three cases of loadings and notice that the result of the simultaneous vector minimization of the two independent compliances corresponding to the two independent loadings that was transformed into a problem of scalar-valued minimization of the weighted sum (12) of the two compliances (scalarization of the vectorvalued objective function), significantly differs from the minimization of the one, scalar-valued compliance corresponding to the weighted sum of the two loadings with the Distribution of the optimal λ 1 , λ 2 and C ij kl for η = 0.5, i, j, k, l = 1, 2, and optimal λ 1 andC ij kl i, j, k, l = 1, 2 (contour layouts in ChromaDepth scale) for the case a (left column) b (middle column) and c (right column) same value of the parameter η. The last approach could be interpreted as an alternative and a significantly simplified way of finding the optimal layouts in the case of two independent loadings. The optimal values of the objective function ( We see that the optimal value of the scalar-valued compliance corresponding to the resultant of the two loadings (J = 0.0215) is better than the both optimal values (J = 0.0636 and J = 0.1613) calculated on the basis of the vector optimization algorithm. However, the optimal layouts of the Hooke tensor calculated by scalarization of the vectorvalued objective function are a better compromise because they apply to the efficient points and Pareto optimal solutions, which is the best what mathematics can do at this stage of the setting of the problem.
Let us note that although the loads in Example 2a are either symmetric or skew-symmetric, the optimal distributions of the Kelvin moduli are symmetric or skewsymmetric in each case, see Figs. 8-10, the left columns. This property of the solution follows from the independent treatment of both the loads in the formulation (16): the trial stress fields τ α are taken from two independent sets. Contrary, the results of the Example 2c are asymmetric, because the only load applied is asymmetric.
Example 3 Consider the cantilever supported at two nonsliding supports at its left corners. The plate (see Fig. 11) is subject either to a horizontal load 1 T = 1 T 1 e 1 or to the vertical load 2 T = 2 T 2 e 2 at its right edge. Distribution of the values of the optimal Kelvin moduli λ 1 , λ 2 (λ 1 λ 2 ) and optimal components C ij kl referred to the (x 1 , x 2 ) coordinate are shown for the weighting factors η = 0.9 and η = 0.1, see Fig. 12.
The optimal values of the objective function occur to be J η=0.9 = 0.0514 and J η=0.1 = 0.1921.  In this example we have also tested numerically mesh independence of the presented algorithm. A rectangular plate was additionally divided into n x × n y = 20 × 10 = 200 and n x × n y = 60 × 30 = 1800, 4-node, quadrilateral, isoparametric finite elements. The optimal values of the objective function are set up in Table 1.
Similar comparison of the layouts of the optimal Kelvin moduli λ i is shown in Figs. 13 and 14.
The layouts shown in Figs. 13 and 14 and the values of the objective function set up in Table 1 confirm stabilization of the results as the mesh is finer and finer. Similar tendency was noted in all examples discussed in the present paper.
Remarks 1 Minimization of the functional (73) was performed by the numerical routine frprmn(...) that implemented Fletcher-Reeves-Polak-Ribiere algorithm or by the numerical routine dfpmin(...) based upon the Broyden-Fletcher-Goldfarb-Shanno variant of the Davidon-Fletcher-Powell algorithm. The C-codes of both the routines were implemented into the program originally published in the book by Press et al. (1992). The procedure was treated as the ,,black box" delivering the correct results, hence no iteration histories of the values of the functional (73) was checked. The initial value of the functional (73) depended on the randomly assumed values of the design parameters (86). Let it be stressed here that in all cases of randomly assumed values of the design parameters (86) we have always got the same numerical value of the objective function (73) and the same layouts of the design parameters. The number of iterations in the main loop of the above routine frprmn(...) (or dfpmin(...)) was observed to be dependent on the assumed loads cases, boundary conditions and numerical tolerances, but it was not bigger than 100. The time of executing the program depended heavily on the number of finite elements.

Comparison of the optimum free material designs with the optimum solutions of the fixed values of the Kelvin moduli
The designs reported in Section 5 correspond to the optimum simultaneous choice of all the parameters which determine at each point of Ω the tensor C of elastic moduli according to the spectral decomposition (7). On the other hand, the optimum designs presented in the previous paper (Czarnecki and Lewiński 2011) were found under the condition of the Kelvin moduli being fixed viz. the Kelvin moduli were not subject to optimization. A comparison analysis of both optimum design methods is pending. The aim of the present section is to assess the influence of releasing the Kelvin moduli on the final optimization result. Let us note that the tensors ω K determine non-dimensional characteristics of the underyling microstructure, like angles between the fibres, while the Kelvin moduli measure elastic stiffnesses. Thus one should expect that the optimization over all characteristics of the Hooke tensor, including the Kelvin moduli, shall lead to much stiffer designs, better formed to carry the given load or the system of loads.
Note that in the considered vector optimization for two load cases not one but a family of Pareto optimal solutions is constructed, indexed by the parameter η. A comparison of solutions in the sense of the work Czarnecki and Lewiński (2011) with the FMD solutions discussed here cannot be done directly, since we face the problem of comparing two sets of solutions.
Let us discuss this question with a greater detail. Let us formulate the problem of Section 3.3 in terms of Pareto optimization. Introduce a vector valued mapping whereĤ (Ω) is defined in Section 2 while Υ α is the compliance corresponding to the α-th load, see (11). The Pareto problem reads: among the Hooke tensor fields of clasŝ H (Ω) find all efficient points along with the pertaining Pareto tensor fields C * ∈Ĥ (Ω). Let us recall that Υ * ∈ Υ (Ĥ (Ω)) is an efficient point with respect to the order relation ≤ in R 2 if and only if one cannot find Υ = Υ * in Υ (Ĥ (Ω)) such that A point C * ∈Ĥ (Ω), Υ * = Υ (C * ) is called Pareto optimal if Υ * is efficient. Hence, in our case, the aim of the vector optimization is to find all efficient fields Υ * ∈ Υ (Ĥ (Ω)) ⊂ R 2 along with the Pareto optimal fields C * ∈Ĥ (Ω) pertaining to them. Thus the task is to find not one optimal field but the set of optimal fields. From this Pareto set of fields we can choose the particular solution using additional criteria not yet taken into account in the formulation of the problem. Among the many methods of finding Pareto optimal minimizers, the most popular and perhaps the simplest one is the weighting method that transforms the vector optimization problem into a problem of a scalar-valued optimization.
This family constitutes the efficient points being the images of the Pareto optimal solutions of the problem (96). Each efficient point Υ * of the contour of Υ (Ĥ (Ω)) is associated with a vector η perpendicular to the tangent line to the contour of this set, see Fig. 15a. The pre-images of the contour points are the Pareto solutions. The nine efficient points corresponding to the results of Example 3 (for the 40 × 20 mesh density) are found for the nine subsequent values: η = 0.1, 0.2, . . . , 0.9, see Fig. 15b. These efficient points are the images Υ (C * i ) of the nine various Pareto tensor fields C * i ∈Ĥ (Ω), i = 1, . . . , 9. The two of them for η = 0.1 and η = 0.9 are shown as Pareto optimal C * ij kl components of the Hooke tensor at the right and left column of Fig. 12, respectively. The two Cartesian coordinates (Υ 1 , Υ 2 ) of the nine efficient points are the two compliances of the nine optimally designed cantilevers corresponding to the nine Pareto optimal tensor fields C * i , i = 1, . . . , 9 and calculated in the cases of either the horizontal load T 1 (T 2 = 0) or the vertical load T 2 (T 1 = 0), respectively.
The two optimization problems: over ω K (considered in Czarnecki and Lewiński (2011)) and over (ω K , λ K ) put forward in the present paper determine two sets of Pareto efficient solutions. A comparison of both the sets seems impossible, because a uniform choice of the weighting coefficients η i from the interval [0, 1], i.e. the choice which divides the interval into equal segments, results in an irregular set of Pareto efficient points, see e.g. Hillermeier (2001), p. 21.
The comparison performed below will be confined to the case of η = 0.5. We consider the problem of Fig. 16. The plate is subject to two kinds of loads: symmetric and skewsymmetric. In the first step we solve the FMD problem (16) by solving the problem (63) for η = 0.5 and E o = 1 Pa. The optimal values of λ * 3 become zero. Since the loads are independent the optimal layouts of the moduli λ * 1 , λ * 2 are symmetric, see Fig. 17. The optimal layouts of the moduli C * 1111 , C * 1122 , C * 1221 , C * 2222 are symmetric while the layouts of C * 2212 , C * 1112 are skew-symmetric with respect to the vertical symmetry axis of the plate, see Fig. 18, the results at the left column.
Note that the effective domains of λ * 1 , λ * 2 (or the domains where they are not zero) form very thin strips.
Consider now the same plate of fixed and uniformly distributed Kelvin moduli of the values λ 1 = 0.9 Pa, λ 2 = 0.09 Pa, λ 3 = 0.01 Pa (97) within the domain Ω. The moduli thus fixed satisfy the isoperimetric condition (15) with the same value of E o as assumed in the results in Figs. 17 and 18. Optimization with respect to ω K leads to the Pareto optimal layouts of the moduli C # ij kl , see the middle column of Fig. 18. These layouts differ considerably from the layouts of C * ij kl found by the FMD optimization.
Let us check now whether the FMD layouts of C * ij kl can be repeated by performing optimization over ω K of the plate with the Kelvin moduli found by the FMD method, or Fig. 16 The plate subject to α T : the symmetric load (α = 1) or to the skew-symmetric load (α = 2)

Fig. 17
The solution to the problem in Fig. 16. Distribution of the Pareto optimal Kelvin modulus λ * 1 (first row) and λ * 2 (second row) for η = 0.5 (Scatter, Contours and Height Field plots in Rainbow Color map, respectively) with the Kelvin moduli λ * K reported at the left column in Fig. 17. This check is made difficult due to λ * 3 = 0, which is not allowed in the optimization over ω K as described in Czarnecki and Lewiński (2011). The condition λ 3 = 0 is replaced by λ 3 = 0.001. Moreover, minimizers of the problem (63) vanish, or almost vanish on some sub-domains of Ω and, consequently, the Kelvin moduli given by (39) assume there zero or very small values. Since the algorithm of Czarnecki and Lewiński (2011) of optimization over ω K requires positive values of the fixed Kelvin moduli, it was necessary to assign small values to λ 1 and λ 2 in those sub-domains. The optimal distributions of the Hooke tensor components Cî j kl , being the results of this optimization is shown in the right column of Fig. 18. We note that these layouts are similar to the optimal layouts found by the FMD algorithm, displayed in the left column of Fig. 18.
Let us discuss the question to what extent the optimization over the Kelvin moduli contributes to the decrease of the compliance. To this end consider the plate of Fig. 16, now subjected to one force of components (1 − η)T and ηT . We assume η = 0.5. If this plate is optimally designed by the FMD method, where both ω K and λ K are design variables, its normalized compliance equals J * = 0.0082. If this plate is made from the material of Kelvin moduli (97) with optimal eigenstates ω K found by the method of Czarnecki and Lewiński (2011) its normalized compliance equals J # = 0.0849, or about 10 times more. A comparison of this type has been reported in the paper Czarnecki and Lewiński (2012) concerning the single load case. The compliance of an optimal non-homogeneous and anisotropic plate was reported to be 7.2 times smaller than the compliance of an equivalent isotropic plate. As was to be expected the optimization over all characteristics of the Hooke tensor leads to much stiffer structures than optimization over the eigenstates ω K .

Remarks 2
The scalar optimization problem (16), aimed at minimizing a weighted sum of the compliances, is the only one from among many possible replacements of the primary vector optimization problem (94) by a suitable scalar optimization problem that allows to find all Pareto (called also Edgeworth-Pareto) optimal points. From among many scalarization methods used in the modern vector optimization (see the survey in Chapter 3.2 in Hillermeier (2001)) besides the most widely used (and used here) the weighted sum approach, the another method -weighted Chebyshev norm approach is frequently recommended -especially for the nonconvex multiobjective optimization problems, cf. Section 11.2.2, pp. 304-312 in Jahn (2004).

Final remarks
In virtue of expressing the compliance of the elastic body by the minimum of the complementary energy the optimum design problem of minimization of the weighted sum of n compliances over parameters characterizing the anisotropic properties of the body has been rearranged to a minimization problem involving: the behavioral variables (the components of n stress fields) and design variables (the Fig. 18 Solution of the problem of Fig. 16. Distribution of the optimal C * ij kl , C # ij kl and Cî j kl for η = 0.5, (scatter plot layouts in Rainbow Color map) parameters of the spectral decomposition of the Hooke tensor). This problem has been reduced to a new effective problem: to the minimization of an effective complementary energy of a mixture of materials with locking properties. This energy has the integrand of linear growth, a specific feature leading to non-smoothness of the solutions, experienced by the mathematicians engaged in the celebrated minimal surface problem, also characterized by the functional of the integrand of linear growth.
The locking material problems are expressed by two mutually dual problems: static and kinematic. In the present paper both the problems have been formulated, yet the static problem only has been the subject of the numerical treatment. This choice has been dictated by a direct link between the minimizers of the static problem and the optimum values of the design variables: the method developed delivers explicit rules for the latter formulae.
The effective stress-based problem involves minimization over statically admissible stress fields, corresponding to the independent loads applied. Any numerical approach necessitates interpolation of such fields. This interpolation problem has been solved by using the SVD decomposition method. Then the minimization operation has been performed by applying carefully chosen gradient oriented non-constrained optimization solvers. Indeed, the SVD decomposition introduces free parameters, not subject to any constraints, which paves the way for the unconstrained optimization.
A characteristic feature of the FMD solutions is that the stress fields which solve the locking problem (38) can vanish on some sub-domains of the design domain. There the material is not necessary. Such behavior of the solutions is admissible within the mathematical setting due to the linear growth of the integrand of the minimized functional in (38). Therefore, the method put forward discloses the domains where the material is not necessary, which makes it possible to circumvent other optimization methods aimed at locating the holes in the design domain as well as the shape optimization methods. The optimal shape of the domain is directly determined by the effective domain of the stress fields being minimizers of the stress-based auxiliary locking problem (38).
The method applies to the three dimensional case. To make it numerically efficient one should develop fast codes for solving the SVD problem. In the 3D case six stress components are assigned to each node of the mesh. Consequently, the number of the unknowns is much bigger than in 2D case, where three stress components at a node should be found. Moreover, the number of Kelvin moduli is six, and six tensors ω K should be determined. A huge number of unknowns will be an obstacle difficult to overcome, unless efficient parallelized SVD codes are developed.