Constraint interface preconditioning for topology optimization problems
Abstract
The discretization of constrained nonlinear optimization problems arising in the field of topology optimization yields algebraic systems which are challenging to solve in practice, due to pathological illconditioning, strong nonlinearity and size. In this work we propose a methodology which brings together existing fast algorithms, namely, interiorpoint for the optimization problem and a novel substructuring domain decomposition method for the ensuing largescale linear systems. The main contribution is the choice of interface preconditioner which allows for the acceleration of the domain decomposition method, leading to performance independent of problem size.
opology optimization, domain decomposition, NewtonKrylov, preconditioning, interior point
65K10, 65N55, 65F10, 90C51, 49N90
1 Introduction
The aim of topology optimization is to determine the optimal distribution of a certain amount of material within a prescribed value in order to minimise the strain energy (or compliance) of a given structure. The distinguishing feature separating this approach from shape optimization involves the introduction of new boundaries, allowing for the consideration of a broader range of feasible solutions. The main problem when pursuing such an approach is that a large number of design variables are required in the discrete formulation in order to maintain the quality of the contours in the final design. Solutions are typically obtained through the use of iterative optimization techniques, requiring repeated discretizations via the finite element method, corresponding to a sequence of linearized problems. As a result, even problems resulting from using relatively coarse discretization parameters can be computationally demanding. Our aim in this paper is to introduce solution methods adapted to the complex nature of this class of problems.
Attempts to alleviate such difficulties can involve the application of a faster finite element solver, or the use of efficient discretisation techniques [9]. Standard approaches based around Picard iterations target the illconditioned equilibrium equations, from which the bulk of computational effort resides. In [31], MINRES coupled with recycling is explored based on the observation that the densities are only expected to undergo minor changes after a relatively small number of iterative steps. The illconditioning is dealt with through a preconditioning strategy involving both rescaling and an incomplete Cholesky decomposition.
In terms of parallel computing, the application of the preconditioned conjugate gradient (PCG) method coupled with Jacobi preconditioning has been considered in a number of references, including [7, 11, 18, 22, 30]. Two additional approaches are considered in [30], namely preconditioning based on an ILU factorisation, as well as condensation through substructuring coupled with a diagonal preconditioner for the resulting interface problem.
Alternatively, primaldual Newton methods can be considered, with particular focus on interior point approaches. Examples illustrating the application of such approaches for solving large scale topology optimization problems can be found in [3, 14, 21]. The KKT conditions from the resulting nonlinear equality constrained optimization problem are then solved using Newton’s method. Evidently, for large scale problems, obtaining solutions to the resulting system of equations will become expensive and even prohibitive in certain cases. In [21], Maar and Schultz applied multigrid to the resulting system, and from their results were able to witness an approximately linear overall complexity with respect to the number of unknowns used in the problem.
In this paper, we propose to apply domain decomposition to the resulting Newton system described above, which will then be solved using GMRES coupled with an appropriate preconditioning strategy. An important component within our preconditioner is based on targeting the resulting interface problem, which will be achieved through the consideration of an appropriate fractional Sobolev norm. Our paper will conclude by illustrating that results can be obtained without dependence on the chosen mesh parameter.
2 Problem description
Consider a material occupying an open and connected domain with Lipschitz boundary . We assume that the elasticity tensor for the material can be modeled as , where is the density, and is a prescribed constant tensor. This expression describes a common approximation of the material properties known as the Solid Isotropic Material with Penalisation (SIMP) model. In this paper we will be interested in the Variable Thickness Sheet (VTS) problem, corresponding to the choice , for which existence of solutions holds (see, e.g., [4, p.272–274]). We assume that the material is clamped (i.e., we assume Dirichlet conditions) along the subset of the boundary and that both body forces and boundary tractions act on the material, resulting in a displacement which satisfies the following equilibrium equations in :
where, by Hooke’s law, the stress tensor is .
We employ next a weak formulation in order to define our topology optimization problem. The natural spaces arising in this context are the standard spaces and
For the description of the preconditioners associated with a domain decomposition approach with corresponding interface between subdomains, we will require the fractional Sobolev space
which is an interpolation space of index 1/2 (see [19] for details).
2.1 Weak formulation
Let and let and be defined below
where is the strain tensor corresponding to a displacement and
We seek such that for all
In order to formulate our topology optimization problem, we define the following admissible set for our design variable :
where denotes the amount of material available, and and denote upper and lower limits on the density function, respectively.
Consider now the following nonlinear minimization problem:
(1)  
subject to:  (2) 
The state variable corresponds to the displacement of the material , while the design variable is the density function .
Let now denote finitedimensional subspaces of and , respectively and consider the following discrete weak formulation of our minimization problem:
(3)  
subject to:  (4) 
Using finite element bases for and , the minimization problem (3–4) yields the following discrete constrained minimization problem:
(5)  
subject to:  (6)  
(7)  
(8) 
where . Working with a piecewise constant approximation , we can write the assembly process of the stiffness matrix as follows:
where is global representation of the elemental matrix corresponding to simplex . This will allow for a certain simplified expression for the Jacobian matrix when considering the firstorder conditions for our minimization problem.
3 Interior point method
Traditional and popular approaches for solving optimization problems such as (5–8) involve separate treatment of both the design objective and the equilibrium equations. Typically, for an initial given design the stiffness matrix is assembled and used to solve the equilibrium equations for the displacement . This is then used to obtain an appropriate update to the design variables, which is then checked for suitability based on previous values. If an appropriate solution has yet to be found, the process is repeated. Typical approaches used to obtain an update to the design variables include both the Optimality Criteria (OC) method [4, p.308] and the Method of Moving Asymptotes (MMA) [26], usually followed by both sensitivity and filtering analysis to cater for the general SIMP setting [25].
More recently, fullycoupled approaches have been receiving considerable attention within the PDE constrained optimization community. In these approaches, all constraints are included and no subproblems are being solved separately. An important feature within these methods is that the equilibrium equations are embedded within the optimization routine, allowing for the simultaneous treatment of all constraints within the problem. Examples highlighting the benefits, and in particular the savings in computational time for such methods, can be found in a number of sources, including [5, 6, 14]. We describe below an interior point method as applied to the topology optimization problem introduced in (58).
Interior point methods are used to solve both convex linear and nonlinear optimization problems iteratively by considering updates confined to the feasible region (cf. [8, 10, 33]). As well as obtaining solutions in a polynomial time, these methods have been used to determine solutions to previously intractable problems, meaning that they are useful from both a theoretical as well as a practical view point. We begin by rewriting the formulation (5) slightly to incorporate the inequality constraints within the objective function. This will be achieved through the use of logarithmic barrier terms as illustrated below:
(9)  
subject to:  (10)  
(11) 
with and where .
The stationary points are defined by setting to zero the relevant partial derivatives of the Lagrangian (12):
(13)  
(14)  
(15)  
(16) 
In the above, and
It is important to note that the condition number of the Hessian of the Lagrangian may pose an issue for densities close to either or as both and tend to zero. To alleviate this issue, auxiliary nonnegative variables and are introduced in the following way
(17) 
Using (17), we see that
where and . Through this substitution, and the elimination of the Lagrange multiplier (which, by (13) and (15), is equal to ), the first order optimality conditions can be written as
By setting , Newton’s method applied to the above nonlinear optimality conditions has the following form
(18) 
where the Jacobian matrix is given below
(19) 
Remark 3.1
In the following, we will assume that problems of the form (9–11) yield Jacobian matrices of the form (19), which are nonsingular for in a neighbourhood of the solution. This property will be assumed to hold for both Dirichlet and mixed boundary conditions, due to the wellposedness of problem (1–2) for both mixed DirichletNeumann and Dirichletonly boundary conditions.
Despite being both nonsymmetric and indefinite, its condition number is expected to be bounded under reduction of the barrier parameters and . Therefore, it is important to consider appropriate strategies for obtaining an accurate update through (18). Work by Forsgren, Gill and Shinnerl [12] uses the diagonal structure of both and to transform into a symmetric matrix. Another possibility is to consider appropriate techniques to condense the matrix via block elimination or using a Schur complement approach. This approach is known to lead to illconditioning; however, the effect on the accuracy of the resulting solution can be benign, as discussed by Wright in [32]. Therefore, block elimination techniques remain a practical option, particularly for the situation where the original matrix is large. The drawback in this case is the loss of sparsity and the lack of obvious (and efficient) preconditioners. For this reason, we look to exploit the sparsity present in the original unreduced system by using an iterative method coupled with an appropriate preconditioning strategy, with the original matrix shown to exhibit favourable spectral properties as shown in recent work by Greif et. al. [13]. The preconditioner employed for this work will be based on a decomposition of the domain into subdomains, described in detail from the next section onwards. In particular, we provide a description of an interface preconditioning technique which was first introduced in [2] for solving scalar elliptic problems, and which is now adapted to the case of our constrained PDE problem.
4 Domain decomposition
In order to formulate a decomposition of our problem that is suitable for a parallel environment, we need to ensure that aside from the physical decomposition of the domain, the resulting subproblems are wellposed. It turns out that we can achieve this through a simple reformulation of our minimization problem, which targets the mass constraint (11).
4.1 Standard definitions and notation
Consider a subdivision of into nonoverlapping subdomains with boundaries such that
We denote the resulting interface by :
For each , we define the nodal index set to be the set of nodes strictly contained in and the simplex index set to be the set of indices of all simplices contained in . We further define .
4.2 Problem reformulation
One of the obstacles in decomposing the VTS problem (5–8), or equivalently the interior point formulation (9–11), is the global mass constraint (7), or (11), respectively. A suitable way to ’decompose’ this constraint is provided by the following equivalent formulation:
where and represent the respective subvectors of and corresponding to the index set . We note here that this approach introduces additional unknowns , together with additional constraints. Note also that given the definition of , represents the mass of subdomain .
The modified formulation corresponding to (9–11) is included below.
(20)  
subject to:  (21)  
(22)  
(23) 
The above problem is equivalent to minimization problem (58) and hence is wellposed.
Following the same procedure of differentiating the corresponding Lagrangian function, one can derive the firstorder optimality conditions and set up Newton’s iteration in the form (18). The resulting Jacobian matrix, also denoted by , now has the form
(24) 
where is described below
Using a node ordering comprising nodes interior to the subdomains followed by nodes on the interface , the Jacobian will have the following permuted block form (with subscripts and indicating this ordering)
(25) 
Under a further permutation which lists all the unknowns corresponding to subdomains , for each , it can be seen that it has the blockdiagonal structure . Each block represents the instance of the Jacobian for a minimization problem posed over the subdomain with Dirichlet boundary conditions on and which has the familiar algebraic form (19):
Here, are the counterparts of assembled on the interior of and with Dirichlet boundary conditions applied. By the wellposedness of problems of the form (9–11) for the case , we conclude that the associated firstorder optimality conditions are wellposed and yield a Jacobian matrix which is nonsingular (cf. Remark 3.1).
The directsum structure of was achieved through the reformulation of the mass constraint and clearly allows for a parallel implementation of the inverse of , which will be a key building block in the solution algorithm described in the next section. It is also an interesting feature in itself to be able to decompose the global Jacobian into local Jacobians associated with similar local minimization problems. One consequence is that the local problems inherit the wellposedness associated with the global problem; the resulting decomposition is hence useful, with the local problems invertible independently of each other.
4.3 A DirichletDirichlet approach
At each step of the Newton iteration we need to solve linear systems of the form
where is the Jacobian matrix (25) corresponding to the first order conditions for the reformulated problem (20–23). Due to the size and structure of , we seek to solve such systems using an iterative solver with a suitable preconditioner. A proven candidate is the block uppertriangular matrix below employed as a right preconditioner
(26) 
where is an approximation to the Schur complement matrix
Due to the directsum property of , this block approach can be classified in domain decomposition terminology as a nonoverlapping DirichletDirichlet procedure, where each subdomain block represents a Jacobian matrix arising in some topology optimisation subproblem posed on where the material is clamped on . Moreover, at each Newton iteration, the Schur complement can be seen as the finite element discretization of a generalised SteklovPoincaré operator corresponding to the interface problem generated by the decomposition. We remark here again that this approach is only available via the reformulation described above
The preconditioner inverse can be written as the product of three matrices:
(27) 
The application of would initially involve the action of on the skeleton problem corresponding to the interface . Next, a boundarytodomain update would be applied through , before applying the inversion in parallel of on subdomains. With the exception of , the potential for parallelism in (26), or equivalently in (27), is evident. Therefore, the task is to seek an appropriate representation to so that the preconditioner can be assembled, stored and applied in an efficient manner. This is discussed in detail in the next section.
5 Constraint interface preconditioners
The Schur complement has the following block structure
(28) 
where the matrices , and . The matrix is negative definite and has a nearlydiagonal structure with reducing entries as the Newton iteration progresses, while can be computed cheaply in parallel. The main focus will therefore be the matrix which is associated with the interface displacement nodes. This block dominates the Schur complement for the Jacobian, and so the aim is to design a preconditioning procedure based on the structure of above and in a manner that provides a suitable approximation of . We note here that one can view the blockstructure (28) as the discretization of an operator corresponding to and constrained by conditions incorporated in the remaining blocks.
Remark 5.1
The Schur complement is symmetric and indefinite, with a negative inertia equal to , where is the number of subdomains. Standard optimization approaches, such as projected CG, cannot be used in this case. Similarly, we expect positive definite preconditioners to be less effective. Consequently, we devised a novel approach which incorporates the indefiniteness (and the constraints) into our preconditioner. The resulting method can be loosely described as being of constrained type (cf. [17]).
Given the structure of we propose a preconditioner with a similar block form but with replaced by a suitable approximation. A direct calculation using the block form (25) of the Jacobian matrix yields
(29) 
where
(30) 
is the Schur complement arising in a DirichletDirichlet nonoverlapping domain decomposition method for the elasticity equations. The above splitting of suggests the following candidate for a preconditioner:
(31) 
This choice is computationally expensive due to the dense matrix and can only be seen as an ideal preconditioner, very much like the Schur complement itself. However, the (11)block can afford approximations which make the application of efficient. In the following, we restrict our attention to . This matrix is the finite element representation of the SteklovPoincaré operator corresponding to the interface problem for the elasticity equations, which is known to be continuous and coercive on . The above properties can be shown using the continuity and coercivity of the bilinear form on ; for details, see [24]. The restriction to , in the context of the finite element discretization of our problem, preserves these properties; in turn, can be shown to be spectrally equivalent to the matrix representation of a norm on . We describe this discrete norm below.
5.1 Discrete fractional Sobolev norms
In order to describe the relevant norms further, we first describe the relevant function spaces defined on the interface . Let represent the tangential gradient of a scalar function such that
corresponding to the projection of the gradient of onto the plane tangent to at the point . We now define
Let denote the set of points on the Dirichlet boundary of the domain which are also on the interface ; when this set is nonempty, we define the space
We recall here that the space introduced above is the fractional Sobolev space of index 1/2. More generally, we define the scale of spaces as the interpolation spaces of index corresponding to the pair of spaces in the sense of Lions and Magenes [19]:
Using a formulation involving discrete interpolation spaces with application to finite element spaces, one can define a finite dimensional space , together with the following matrix representation of a discrete interpolation norm [1]
(32) 
where and are the mass and Laplacian matrices assembled on using the restriction of the finite element basis for to the interface. Using this norm representation, the continuity and coercivity properties of the elasticity Schur complement translate into the following spectral equivalence [27, p. 129]
where . This equivalence is the basis of our candidate for the approximation of the Schur complement .
The matrix is in general full and expensive to compute. However, an implementation of the action of the inverse of on a given vector can be achieved cheaply using a Krylov subspace approximation constructed via a generalised inverse Lanczos iteration which only requires the application of the inverse of the sparse matrix for a small number of steps [1]. In the following section we indicate how to extend this procedure to the case of our proposed preconditioner.
5.2 Constraint preconditioners
Given the constraint form (28) of the interface Schur complement, we propose the following choice of preconditioner which preserves the constraining blocks while replacing the (1,1)block by a spectrally equivalent matrix:
(33) 
The notion of constraint preconditioning has been extensively analysed in the context of saddlepoint problems [17, 20]. In our case, the matrix structure is of a different nature; however, given the spectral equivalence between and , we hope to achieve a useful similarity between and . Our motivation is the following result. {proposition} Consider the generalized eigenvalue problem
where
with nonsingular. Let be a basis for the nullspace of . Then

with multiplicity ;

the remaining eigenvalues satisfy the eigenvalue problem
where .
Remark 5.2
The above result indicates that the increase in size by due to our problem reformulation (20–23) is automatically taken care of by our constraint interface preconditioner, which maps eigenvalues to 1. The remaining eigenvalues will depend on the closeness of our preconditioner to . As described above, we tried to incorporate this property by a choice of that is spectrally equivalent to the dominant part of , namely .
The main issue with the definition of is the full matrix arising in the (1,1)block. This matrix needs to be computed via an expensive matrix squareroot calculation; moreover, it needs to be employed in order to implement the action of the inverse of on a given vector. A practical alternative can be derived from the following constrained Lanczos decomposition of .
Let
denote the generalized Lanczos decomposition of the pencil in exact arithmetic, where is an orthogonal matrix and is a tridiagonal, symmetric and positivedefinite matrix. Define . Then (see also [1]). Consider the QRfactorization
where is orthogonal and define
We obtain the following orthogonal factorization of
We will refer to the above representation of , seen as a
twobytwo block matrix, as the constrained
Lanczos factorization of .
Given our preconditioning task, consider now the product
An approximation to can be constructed using a partial factorization:
where we define
The preconditioning operator implicit in the definition of will be denoted by :
Its implementation requires:

a partial Lanczos factorization of the (1,1)block of ; this yields a matrix of Lanczos vectors and a block tridiagonal matrix ;

a QRfactorization of the (1,2)block of multiplied by ; this yields the factors arising in and , respectively.
The number of Lanczos vectors in is expected to be small; correspondingly, the sizes of will also be small and the resulting matrix will be easy to invert. The overall complexity for the above procedure is of order . Thus, the application of the interface preconditioner will not dominate the cost of a subdomain solve provided we work with subdivisions for which , where we assumed that the cost of inverting a subdomain Jacobian matrix is the cost that a direct method requires to invert a banded matrix of size with bandwidth .
6 Numerical experiments
In order to illustrate the performance of the solution method described in Section 6.1, we consider a model test problem based on compliance design. The problem involves a cantilever beam posed on a rectangular design domain as illustrated in Figure 1, with clamping applied along the left edge and a force applied in the middle of the right edge. The density contour plot corresponding to the optimal design is displayed in Figure 2. Symmetry is a feature which may be exploited computationally; however in this paper we chose to retain the original design domain in order to test the performance of our solution method on the full problem.
6.1 Implementation details
Finite Element Method
We use a subdivsion of into square elements of size . We seek approximations
for all , where denotes the space of degree polynomials defined on . We employ uniform refinement, in order to exhibit the behaviour of our preconditioning technique with respect to the mesh parameter . We note here that the space is not constructed explicitly, and that it was only introduced in order to provide the mathematical description of our problem.
Interior Point Algorithm
We use a textbook interiorpoint algorithm; the key details are included below; for full details see [23, Ch. 19].
We repeat the following steps until convergence:

Solve a sequence of Newton systems (18) to get .

Find the step length (see below).

Update the solution

Update the penalty parameters via
We start with and continue until they are brought below a tolerance of . A more sophisticated version of the code, with an adaptive choice of the penalty parameters and can be found in [15]. However, we have never observed any difficulties with this simple update of and .
Step length. In the interiorpoint method we cannot take a full Newton step resulting from the solution of (18), as this would lead to an infeasible point. We propose a technique that will effectively involve finding such that for such that and such that for such that . Therefore, we consider obtaining and as follows
The constant represents an appropriate shortening of the Newton step to the interior of the feasible region. In the event that both and are greater than , the step will be shortened appropriately:
This procedure is relatively simple and straightforward to both implement. A more sophisticated line search procedure is described in [15] and could potentially be used here. However, it will be illustrated later that the current technique was able to yield desirable results.
Initial approximation. We start the interiorpoint algorithm with a uniform distribution of the design variable chosen with respect to the mass constraint (11)
The initial displacements are then computed from the equilibrium equation
Domain decomposition
We used only regular subdivisions into rectangular subdomains. The corresponding sizes of the resulting interfaces are illustrated in Table 1. It is evident (and wellknown) that in order to balance the complexities of the interface and subdomain problems, the increase in the number of subdomains should be paralleled by a decrease in . Aside from regular decompositions, one could also decompose the domain in an adaptive fashion based on the changing nature of the design. This could be carried out either at each outer iteration, or alternatively once after a fixed number of outer iterations based on the (previously mentioned) observation that the density will only be subject to minor changes after a relatively small number of iterative steps. In terms of a nonregular subdivision, the graph partitioning tool METIS [16] may be used in order to partition the finite element mesh into nonregular subdomains.
41,355  384  1,140  2,604  
164,619  768  2,292  5,292  
657,027  1,536  4,584  10,668 
Gmres
We used an inexact NewtonGMRES method preconditioned by either ; flexible GMRES was used for the case of . The GMRES stopping criterion was the reduction of the norm of the initial residual by a factor of .
Constrained Lanczos factorization
The implementation of the action of requires first to generate the matrices and . This is achieved as an additional preprocessing step involving one set of subdomain solves. The number of Lanczos vectors is taken to be , so that the overall complexity of using the preconditioner is of order .
6.2 Numerical results
Table 2 displays the results for our test case for a range of mesh and subdomain sizes. The results were obtained using a Linux machine with an Intel® Core™ i7 CPU 870 2.93 GHz with 8 cores. The upper and lower limits on the density were set at and respectively, with the permissible volume in each test case defined to be .
The first observation arising from our numerical experiments is that has performance independent of meshsize and almost independent of the number of subdomains; moreover, the number of iterations is low (averages between 5–8 iterations on the largest problem), with only a small departure from the optimal count of 2 iterations corresponding to the case where the exact Schur complement is used. This confirms the earlier assumption that the matrix in (29) is negligible and that the preconditioner is optimal in the same sense as the Schur complement.
Avg GMRES (Newton)  Total GMRES Its.  

No. Subdomains:  4  16  64  4  16  64  
Preconditioner  \backslashbox  0.5  0.6  0.7  0.5  0.6  0.7 
1/64  7.29 (14)  10.57 (14)  13.86 (14)  102  148  194  
1/128  5.93 (15)  8.13 (15)  8.93 (15)  89  122  134  
1/256  5.25 (16)  6.81 (16)  7.50 (16)  84  109  120  
1/64  11.36 (14)  25.64 (14)  35.00 (13)  159  359  455  
1/128  10.60 (15)  21.80 (15)  33.12 (16)  159  327  530  
1/256  10.25 (16)  19.94 (16)  29.82 (17)  164  319  507  
1/64  12.21 (14)  31.36 (14)  50.46 (13)  171  439  656  
1/128  12.00 (15)  27.47 (15)  43.47 (15)  180  412  652  
1/256  11.00(16)  26.41 (17)  40.06 (17)  176  449  681 
With regard to parameter dependence, we note that the number of iterations decreases with in all experiments. This is somewhat expected, given that is essentially an approximate implementation of , which incorporates in the (11) block the finite element discretization of the continuous SteklovPoincaré operator for the elasticity equations. An interesting fact is that the preconditioning technique appears to outperform occasionally the preconditioner ; however, different preconditioners lead to different Newton convergence histories given the adaptive stopping criterion employed, which may result in overall complexities more favourable for the constrained preconditioning approach.
The table also indicates what appears to be a logarithmic dependence on the number of subdomains. We found this difficult to analyze, but this behaviour is not unlike that exhibited by similar substructuring preconditioning techniques in [1]. This suggests that the preconditioner has the ability to match the properties of the discrete fractional Sobolev norm in a constrained setting. This represents a novel approach which could be useful in other constrained optimization settings and for other PDE models.
The values of listed in Table 2 were chosen by experimentation based on similar observations for the linear elasticity problem reported in [28], where it is noted that different values of may be able to provide a closer approximation to the decay of the associated SteklovPoincaré operator. Similar findings were found numerically in this work also, with the best values of used in order to produce the results reported in the table.
Despite the fact that a logarithmic dependence is noted for an increasing number of subdomains, the computational benefits gained as a result of distributing calculations amongst an increasing number of processors can lead to a significant speedup when compared to solving the original problem on the global domain . Exact calculations displaying this behaviour are not included here, however the interested reader is referred to [28], where notable speedup was observed under the preconditioner described in (32) for the Optimality Criteria method for topology optimization, which, unlike the interior point approach in this paper, requires the solution of a linear elasticity problem at each step.
7 Conclusions
We described a novel domain decomposition approach coupled with an interior point method for compliance minimisation problems arising in the field of topology optimization. The problem was reformulated in order to allow for wellposed subdomain problems which could be viewed as local Jacobian solves. This was an important step which allowed for the domain decomposition method to be welldefined. The resulting interface Schur complement problem yielded an indefinite matrix which included a global volume constraint; consequently, an indefinite preconditioner was devised in order to incorporate the properties and the structure of the interface Schur complement matrix. This resulted in a technique requiring the inversion of a socalled constrained discrete fractional Sobolev norm, the application of which was performed via a certain constrained Lanczos procedure. We tested the resulting method on a standard test problem of topology optimization, with experiments indicating independence on the mesh parameter, although a depedence on the number of subdomain was noticed.
Some of our current investigations include nonregular decompositions, as well as adaptive decompositions based on current iterates. The constrained Lanczos factorization, which allowed the sparse implementation of our interface preconditioner, will be the subject of further study. The precise role played by the parameter also requires more analysis and experimentation in order to quantify the most appropriate value of for a decomposition into any given number of subdomains. Finally, the reformulation of the problem resulting in Jacobian subproblems on subdomains is also worthy of further investigation as it points to nonlinear domain decomposition approaches, such as that described in [29].
Footnotes
 School of Mathematics, University of Birmingham, Edgbaston, Birmingham B15 2TT, UK, and Institute of Information Theory and Automation, Czech Academy of Sciences, Pod vodárenskou věží 4, 18208 Praha 8, Czech Republic. The work of this author has been partly supported by the EU FP7 project AMAZE and by grant A100750802 of the Czech Academy of Sciences
 School of Mathematics, University of Birmingham, Edgbaston, Birmingham B15 2TT, UK
 School of Mathematics, University of Birmingham, Edgbaston, Birmingham B15 2TT, UK
References
 M. Arioli, D. Kourounis, and D. Loghin, Discrete fractional Sobolev norms for domain decomposition preconditioning, IMA J. Numer. Anal., 33 (2013), pp. 318–342.
 M. Arioli and D. Loghin, Discrete Interpolation Norms with Applications, SIAM J. Numer. Anal., 47 (2009), pp. 2924–2951.
 A. BenTal, M. Kočvara, A. Nemirovski, and J. Zowe, Free Material Design via Semidefinite Programming: The Multiload Case with Contact Conditions, SIAM review, 42 (2000), pp. 695–715.
 M. P. Bendsøe and O. Sigmund, Topology Optimization: Theory, Methods, and Applications, Springer Verlag, Providence, Rhode Island, 2003.
 G. Biros and O. Ghattas, Parallel Lagrange–Newton–Krylov–Schur Methods for PDEConstrained Optimization. Part I: The Krylov–Schur Solver, SIAM Journal on Scientific Computing, 27 (2005), pp. 687–713.
 , Parallel Lagrange–Newton–Krylov–Schur Methods for PDEConstrained Optimization. Part II: The Lagrange–Newton Solver and Its Application to Optimal Control of Steady Viscous Flows, SIAM Journal on Scientific Computing, 27 (2005), pp. 714–739.
 T. Borrvall and J. Petersson, Largescale topology optimization in 3D using parallel computing, Computer Methods in Applied Mechanics and Engineering, 190 (2001), pp. 6201–6229.
 R. H. Byrd, M. E. Hribar, and J. Nocedal, An interior point algorithm for largescale nonlinear programming, SIAM Journal on Optimization, 9 (1999), pp. 877–900.
 G. C. A. DeRose Jr and A. R. Diaz, Solving threedimensional layout optimization problems using fixed scale wavelets, Computational mechanics, 25 (2000), pp. 274–285.
 A. S. ElBakry, R. A. Tapia, T. Tsuchiya, and Y. Zhang, On the formulation and theory of the Newton interiorpoint method for nonlinear programming, Journal of Optimization Theory and Applications, 89 (1996), pp. 507–541.
 A. Evgrafov, C. J. Rupp, K. Maute, and M. L. Dunn, Largescale parallel topology optimization using a dualprimal substructuring solver, Structural and Multidisciplinary Optimization, 36 (2008), pp. 329–345.
 A. Forsgren, P. E. Gill, and J. R. Shinnerl, Stability of Symmetric IllConditioned Systems Arising in Interior Methods for Constrained Optimization, SIAM J. Matrix Anal. Appl., 17 (1996), pp. 187–211.
 C. Greif, E. Moulding, and D. Orban, Bounds on Eigenvalues of Matrices Arising from InteriorPoint Methods, SIAM J. Optim., 24 (2014), pp. 49–83.
 R. H. W. Hoppe and S. I. Petrova, Primal–Dual Newton Interior Point Methods in Shape and Topology Optimization, Numerical Linear Algebra with Applications, 11 (2004), pp. 413–429.
 F. Jarre, M. Kočvara, and J. Zowe, Interior point methods for mechanical design problems, SIAM Journal on Optimization, 8 (1998), pp. 1084–1107.
 G. Karypis and V. Kumar, A fast and high quality multilevel scheme for partitioning irregular graphs, SIAM J. Sci. Comput., 20 (1998), pp. 359–392 (electronic).
 C. Keller, N. I. M. Gould, and A. J. Wathen, Constraint preconditioning for indefinite linear systems, SIAM J. Matrix Anal. Appl., 21 (2000), pp. 1300–1317.
 T. S. Kim, J. E. Kim, and Y. Y. Kim, Parallelized structural topology optimization for eigenvalue problems, International journal of solids and structures, 41 (2004), pp. 2623–2641.
 J. L. Lions and E. Magenes, Problèmes aux limites non homogènes et applications. I, vol. 3, Dunod, Paris, 1968.
 L. Lukšan and J. Vlček, Indefinitely preconditioned inexact Newton method for large sparse equality constrained nonlinear programming problems, Numer. Linear Algebra Appl., 5 (1998), pp. 219–247.
 B. Maar and V. Schulz, Interior point multigrid methods for topology optimization, Structural and Multidisciplinary Optimization, 19 (2000), pp. 214–224.
 A. Mahdavi, R. Balaji, M. Frecker, and E. M. Mockensturm, Topology optimization of 2D continua for minimum compliance using parallel computing, Structural and Multidisciplinary Optimization, 32 (2006), pp. 121–132.
 J. Nocedal and S. Wright, Numerical Optimization, Springer, 2nd edition ed., 2006.
 A. Quarteroni and A. Valli, Domain Decomposition Methods for Partial Differential Equations, Numerical Mathematics and Scientific Computation, The Clarendon Press Oxford University Press, New York, 1999. Oxford Science Publications.
 O. Sigmund, A 99 line topology optimization code written in Matlab, Structural and Multidisciplinary Optimization, 21 (2001), pp. 120–127.
 K. Svanberg, A class of globally convergent optimization methods based on conservative convex separable approximations, SIAM Journal on Optimization, 12 (2002), pp. 555–573.
 J. Turner, Application of Domain Decomposition to problems in Topology Optimization, PhD thesis, University of Birmingham, http://etheses.bham.ac.uk/5842/, 2014.
 J. Turner, M. Kočvara, and D. Loghin, Parallel Solution of the Linear Elasticity problem with applications in Topology Optimization, in Proceedings of the 4th Annual BEAR PGR Conference, University of Birmingham, UK, arXiv:1501.06211v2 [math.OC], 2013.
 , A nonlinear domain decomposition technique for scalar elliptic PDEs, in Domain Decomposition Methods in Science and Engineering XXI, J. Erhel, M. J. Gander, L. Halpern, G. Pichot, T. Sassi, and O. Widlund, eds., vol. 98 of Lecture Notes in Computational Science and Engineering, Springer, 2014.
 K. Vemaganti and W. E. Lawrence, Parallel methods for optimality criteriabased topology optimization, Computer Methods in Applied Mechanics and Engineering, 194 (2005), pp. 3637–3667.
 S. Wang, E. de Sturler, and G. H. Paulino, Largescale topology optimization using preconditioned Krylov subspace methods with recycling, Internat. J. Numer. Methods Engrg., 69 (2007), pp. 2441–2468.
 M. H. Wright, IllConditioning and Computational Error in Interior Methods for Nonlinear Programming, SIAM Journal on Optimization, 9 (1998), pp. 84–111.
 S. J. Wright, PrimalDual InteriorPoint Methods, vol. 54, SIAM, 1997.