Design of cellular materials for multiscale topology optimization: application to patient-specific orthopedic devices

A flexible problem-specific multiscale topology optimization is introduced to associate different areas of the design domain with diverse microstructures extracted from a dictionary of optimized unit cells. The generation of the dictionary is carried out by exploiting micro-SIMP with AnisoTropic mesh adaptivitY (microSIMPATY) algorithm, which promotes the design of free-form layouts. The proposed methodology is particularized in a proof-of-concept setting for the design of orthotic devices for the treatment of foot diseases. Different patient-specific settings drive the prototyping of customized insoles, which are numerically verified and successively validated in terms of mechanical performances and manufacturability.


Introduction
Cellular materials, characterized by a porous microstructure which properly alternates solid and void, have been engineered in the last years to artificially reproduce the lightweight and the strength properties exhibited by some biological systems, such as bones, honeycombs, sponges, and wood. In parallel, the rise of innovative manufacturing techniques, such as additive processes (e.g., 3D printing), promoted the production of cellular materials in very diverse fields of application, for instance, in healthcare, aerospace, defense, construction, and automotive (Xu et al. 2018;Ivarsson et al. 2020).
The strong interest in cellular materials has favoured the proposal of a wide range of analytical, numerical, and experimental methods for an efficient design of new materials. In this context, topology optimization represents the reference mathematical methodology. Actually, this technique has been extensively adopted not only to create optimized structures (Rozvany 2012;Bendsøe and Kikuchi 1988;Bendsøe and Sigmund 2003;Sigmund and Maute 2013;Allaire et al. 2004), but also to address the design of unit cells optimized with respect to a specific goal (see, e.g., Sigmund 1994). Starting from the seminal work by Sigmund (1994), inverse and direct homogenization techniques have been massively employed to engineer and validate periodic microstructures (Andreassen and Andreasen 2014;Allaire et al. 2019).
However, it is well known that the topology optimization carried out at a monoscale (either at a macroscale or at a microscale level) may limit structural performances. Multiscale topology optimization has been proposed as a remedy to overcome such limitations (see, e.g., Rodrigues et al. 2002;Sanders et al. 2021). This strategy consists in identifying the optimal alternation of void, solid, and microstructures inside the design domain (Auricchio et al. 2020;Arabnejad Khanoki and Pasini 2012). The distribution of void, solid, and microstructures can be carried out by the user through a trial-and-error approach, or automatically 1 3 79 Page 2 of 26 driven by a topology optimization tool. A single microstructure or a multiplicity of different cells can be used to handle the transition areas between solid and void. The first option confers the same effective properties on the whole design domain, while the coexistence of several topologies at the microscale allows us to locally diversify the property of the macro domain. The layout of the employed microstructures can be selected a priori, starting from consolidated dictionaries of unit cells (Allaire et al. 2019;Vigliotti and Pasini 2013), or designed from scratch to match the requirement of interest via homogenization techniques (Coelho et al. 2008;Nakshatrala et al. 2013;Djourachkovitch et al. 2021;Ferrer et al. 2016;Ferro et al. 2020b;di Cristofaro et al. 2021).
The scenario which provides us the highest flexibility with a view to a problem-specific design is the one where the microstructures are ad hoc designed and change across the optimized macro domain. This option clearly represents the most challenging choice in terms of modeling (dedicated optimizations are employed both at the macro-and at the microscale and, possibly, follow different goal quantities, see, e.g., Gao et al. 2019), computational effort (especially, if the optimizations at the macro-and at the microscale are concurrent, see, e.g., Coelho et al. 2008;Xia and Breitkopf 2014) and manufacture (due to intrinsic limits in manufacturability, for instance, to tackle the transition among different cells, or possible defects or irregularities characterizing the optimized layouts, see, e.g., Du et al. 2018).
In this paper, we propose a highly flexible procedure for multiscale topology optimization characterized by a computationally affordable burden, based on a recent and already successfully validated structural design methodology, which provides some improvements in terms of manufacturability. This procedure consists of two phases. During the first one, we create a dictionary of unit cells characterized by the same or by a different topology, sharing a common mechanical objective. In the second phase, we exploit the multiscale topology optimization strategy proposed in Wang et al. (2017), combined with a suitable density material thresholding, to identify the areas of the macro domain to be associated with the different microstructures.
The computational complexity of the proposed procedure remains limited thanks to a sequential optimization of the unit cells and of the macrostructure. The independence between the macro-and the microscale allows us to increase the flexibility of the modeling, since different goal quantities can be adopted to drive the optimization of the cells and of the macro layout. The generation of the dictionary relies on the recent algorithm micro-SIMP with AnisoTropic mesh adaptivitY (microSIMPATY; Ferro et al. 2020b;di Cristofaro et al. 2021;Ferro et al. 2020cFerro et al. , 2019, which combines homogenization techniques together with the solid isotropic material with penalization (SIMP) method and a customized selection of the computational mesh, generalizing the SIMPATY algorithm proposed in Micheletti et al. (2019) to a microscale. SIMPATY algorithm has been successfully applied to the design of structures optimized with respect to the static compliance and to the mass of the final layout Ferro et al. 2020a). In Ferro et al. (2020c) SIMPATY has been enriched by shape optimization to enhance the lightweight and mechanical properties of the optimized configuration, towards an out-of-the-box, free form paradigm.
In Ferro et al. (2019) the authors propose an innovative way to circumvent the computational burden typical of the SIMP approach. The idea is to resort to a proper orthogonal decomposition on SIMP snapshots to predict a rough structure, which is eventually finalized via SIMPATY, in the spirit of an offline/online paradigm (Kunisch and Volkwein 2002).
The first main novelty of this paper is the employment of microSIMPATY algorithm into a multiscale topology optimization setting. This leads us to devise a flexible problemspecific methodology to yield products that match prescribed design and mechanical requirements by alternating diverse microstructures. As a second contribution, we particularize the proposed workflow to a medical context. In particular, we focus on the design of patient-specific insoles for the treatment of foot problems ). Indeed, among the different interventions and pain management actions, foot orthotics represent a very common conservative solution, for foot diseases such as musculoskeletal disorders, foot deformity and ulceration (Iseli et al. 2021). The reduction of plantar pressure represents a reference target to improve patient's condition (Van Netten et al. 2018), typically ensuring the attenuation of symptoms for different pathologies (Collings et al. 2021). For this reason, several methods have been proposed to optimize the design and the production of patient-specific foot orthotics, guaranteeing a therapeutic effect (Ahmed et al. 2020) The methodology proposed in this paper provides a proof-of-concept of a new versatile design pipeline for the manufacturing of patientspecific orthoses.
The paper is organized as follows. Section 2 introduces the SIMPATY algorithm after a brief recap on the SIMP method. Section 3 is devoted to the homogenization techniques, with a particular emphasis on the inverse methodology. In Sect. 4, we exploit inverse homogenization to generate two dictionaries of new unit cells, which differ for the number of diverse microcellular topologies. The adopted multiscale topology optimization is described in Sect. 5. In Sect. 6 we establish a simplified setting with a view to the design of optimized orthotic devices. Two different patientspecific configurations are analyzed for the prototyping of customized insoles. The two scenarios are numerically verified. The most meaningful setting from a medical viewpoint is successively validated in terms of mechanical properties and 3D printing manufacturability. Finally, in the last section, some conclusions are drawn and some hints are provided on possible future developments.

A generic setting for structural topology optimization
Topology optimization aims at allocating the material inside the design domain in order to match a goal quantity strictly related to the application at hand, while satisfying specific physical and design constraints. Originally, topology optimization has been settled to deal with structural problems (see, e.g., the landmark papers by Michell 1904;Rozvany 2012Rozvany , 2009Bendsøe and Kikuchi 1988). In this context, the stiffness, the stress, the vibration modes and the mass of the structure are classical examples of target quantities or constraints to the optimization, in combination with the state equation which models the physical law the structure is subject to. Successively, the crucial impact of topology optimization has been verified in several engineering fields different from structural mechanics (see, e.g., Jenkins and Maute 2016;Dapogny et al. 2017;Villanueva and Maute 2017 and the references therein). The employment of topology optimization in a wide range of different settings motivated the proposal in the literature of several mathematical and numerical methodologies.
Among the most popular approaches, we mention density-based schemes (Bendsøe and Sigmund 2003;Sigmund and Maute 2013), level-set methods (Allaire et al. 2004), topological derivative procedures (Sokolowski and Zochowski 1999;Amstutz and Andrä 2006), phase field techniques (Bourdin and Chambolle 2003), evolutionary approaches (Xie and Steven 1993), homogenization (Bendsøe and Kikuchi 1988), performance-based optimization (Liang 2005), and more recently, model reduction and machine learning tools (Alaimo et al. 2018;Watts et al. 2019;Ferro et al. 2019;Sibileau et al. 2018). The main feature distinguishing these methodologies is the modeling expedient adopted to track the solid/void alternation in the design domain. In this paper, we opt for a density-based scheme, i.e., the SIMP method.
SIMP introduces an auxiliary function ∈ L ∞ ( , [0, 1]) , referred to as density or design variable, which models the distribution of the material in the design domain ⊂ ℝ 2 , here identified with a two-dimensional setting. Function is assumed to take only the extreme values, 0 and 1, where = 1 identifies the material, whereas = 0 characterizes the void. In practice, takes all the values in [0,1]. This leads very often to an over-diffused material-void interface in correspondence with intermediate densities, whose physical interpretation is not straightforward. To overcome this issue, SIMP penalizes the intermediate densities by means of a suitable function of .
The generic formulation of the SIMP method can be formulated as where G(⋅, ⋅) is the target functional to be optimized, the first constraint coincides with the weak form of the equation modeling the physics, the second one enforces a control on the system, with c m and c M the corresponding lower and upper bounds, the last two-sided inequality guarantees the well-posedness of the weak form, being 0 < min < max ≤ 1 , and where U denotes a Sobolev space strictly dependent on the boundary data characterizing the physical model (Ern and Guermond 2004). We remark that the weak form includes the density function according to a law, which depends on the phenomenon under investigation. For instance, a suitable power of multiplies the stiffness tensor when optimizing structures in a linear regime (Bendsøe 1995); in a fluid-dynamic setting, it is standard to consider a density-weighted inverse permeability (Borrvall and Petersson 2003); in the design of microstructures, modifies the homogenization applied at the microscale (Sigmund 1994). Moreover, the system can be multi-constrained. In such a case, relation c m ≤ C(u( ), ) ≤ c M is replaced by a set of constraints (see Sect. 5).
Different algorithms are used in the literature to manage the constrained minimization in (1), such as MMA (Svanberg 1987), Interior Point OPTimizer (IPOPT; Wächter 2002), heuristic-based, genetic or machine learning routines (Chi et al. 2021).
When dealing with a structural topology optimization at the macroscale, it is common to identify the state equation in (1) with the linear elasticity problem in the small displacement regime. The bilinear form becomes where ∶ → ℝ 2 denotes the displacement of the structure; ( ) = (∇ + ∇ T )∕2 is the strain tensor; (1) min is the stress tensor penalized by a power law of the density, p , with p ≥ max{1∕(1 − ), 4∕(1 + )}, the Lamé coefficients, E Y and the Young's modulus and the Poisson's ratio, tr (⋅) the trace operator, and I the identity tensor. Different boundary conditions may characterize the configuration of the structure. The selected boundary data define the space U and the right-hand side, F ( ) , in (1).
Problem (1) combined with the state equation (2) may be characterized by multiple local minima, due to the nonconvexity of the functional G , so that the uniqueness of the solution is a priori not guaranteed. This is, for instance, the case of the minimum compliance benchmark problem (Bendsøe 1995;Sigmund and Petersson 1998).
We discretize the weak form in (1) by means of continuous finite elements, i.e., we look for the displacement, h , and the density, h , in U h = [V r h ] 2 ∩ U and V s h , respectively, with the space of the finite elements of degree q, T h = {K} a conforming triangular tessellation of the domain and ℙ q the set of polynomials of global degree q. Notice that the functions in V 0 h are discontinuous and piecewise constant on T h (Ern and Guermond 2004).
It is well-known that SIMP method suffers from several issues related to the discretization adopted in (1). Among the main drawbacks, we list the dependence of the final layout on the grid T h , the presence of checkerboard patterns, grayscale effects (Sigmund and Petersson 1998;Bendsøe 1995), jagged boundaries and too complex topologies which make the optimized design unpractical for manufacturing.
The mesh dependence can be ascribed to the non-uniqueness of the solution. It is generally tackled by constraining or filtering the design variable.
Checkerboard effects arise when solid and void elements alternate in an uneven way. This is a standard feature for a two-field formulation, unless an ad hoc combination of the finite element spaces for the displacement, h , and the density, h , is adopted (e.g., s = 0 , r = 2 for quadrilateral elements Bendsøe and Sigmund 2003 and s = 0 , r = 1 for triangular elements Bruggi and Verani 2011). Filtering techniques are an alternative viable remedy to get rid of checkerboard patterns, together with adapted computational meshes as detailed in the next section.
Jagged material/void interfaces can be traced back to excessively coarse grids, whereas thin struts often depend on the employment of too fine meshes. As a compromise between too coarse and too fine meshes, in Micheletti et al. (2019) and Ferro et al. (2020c) the authors propose a new algorithm, presented in the next section, which enriches the SIMP method with a customized choice of the computational mesh in the framework of structural optimization.

The SIMPATY algorithm
SIMPATY algorithm has been proposed in Micheletti et al. (2019) to assist the design of lightweight and stiff structures aimed at aerospace applications. This new algorithm consists of an iterative procedure, which sequentially alternates the topology optimization in (1), tackled by the IPOPT package ( Wächter 2002), with the generation of an anisotropic adapted mesh.
SIMPATY algorithm proved to have remarkable features. As an alternative to projection methods (see, e.g., Wang et al. 2011), SIMPATY algorithm allows us to sharply detect the density across the material/void interface thus limiting filtering, i.e., the post-processing required by standard design tools. In addition, the same discrete space can be adopted both for the displacement components and the density with a consequent reduction of the computational effort. These properties identify a very cost-effective design procedure, striking a balance between accuracy and computational demands.
In this section we sketch the SIMPATY procedure as a fundamental building block for the microSIMPATY algorithm employed in Sect. 4 for the design of new dictionaries of optimized unit cells. In particular, the focus will be on mesh adaptation and on the mathematical tool we adopt to drive the generation of a problem-specific computational grid.
Standard isotropic adapted meshes allow us to improve the solution accuracy (for a certain mesh cardinality) or to contain the number of elements (for a user-defined tolerance on the numerical approximation) by optimizing the element size. These improvements can be further enhanced by resorting to anisotropic grids, which tune the size, the shape, and the orientation of the mesh triangles in order to track the directional features of the modeled phenomena, such as steep boundary and internal layers, discontinuities, sharp fronts, shocks in compressible flows, more in general areas where the problem exhibits strong gradients (Formaggia and Perotto 2001;Dompierre et al. 2002;Formaggia et al. 2002;Belhamadia et al. 2014;Porta et al. 2012;Ferro et al. 2018). The generation of anisotropic grids deserves more technicalities when compared with an isotropic context. This justifies the limited availability of anisotropic mesh adaptation modules in current simulation software.
Anisotropic mesh adaptation can be driven by heuristic or theoretically sound quantifiers. The first class exploits information prompted by the application at hand in terms of numerical solution or associated variation (e.g., gradient or Hessian). The second class is represented by the a priori and the a posteriori error estimators, which can be classified according to the controlled quantity adjusting the allocation of the elements in the adapted mesh (Ainsworth and Oden 2000).
After discretizing both the displacement and the density with affine finite elements, SIMPATY algorithm resorts to an a posteriori anisotropic recovery-based error estimator, , associated with the density variable. This estimator turns out to be the ideal tool to identify the steep gradients of across the material/void interface, since it controls the H 1 -seminorm of the discretization error, e = − h . The idea formalized by a recovery-based error estimator is very straightforward (Zienkiewicz and Zhu 1987). It consists in replacing in the definition of |e | H 1 ( ) the exact gradient of the solution with the so-called recovered gradient, namely, Recovery-based error estimators have been exploited in several engineering contexts (see, e.g., Porta et al. 2012;Yan 2001;Mu and Jari 2013). The popularity of these estimators is justified by the good properties they have, among which we cite the independence of the specific problem and of the adopted discretization; the computational cheapness; the handy implementation; the high performance in very diverse fields of application.
Concerning the choice of the recovery operator P in (4), several proposals have been made in the literature, to target specif problem-or theoretical-driven requirements. In general, P(∇ h ) is expected to provide a better approximation to ∇ with respect to ∇ h . Standard rules identify P with a local projection or average of the discrete gradient. Here, we adopt the area-weighted average of ∇ h over the patch of elements, Despite the low polynomial degree, the recovered gradient in (5) proved to be a good approximation for ∇ (Farrell et al. 2011;Ferro et al. 2020a, b, c;Micheletti et al. 2019;Porta et al. 2012), due to the large number of elements contributing to the average.
In 2010, a generalization of the recovery-based error estimators to an anisotropic setting has been proposed in . This extension exploits the anisotropic setting proposed in Formaggia and Perotto (2001). In this paper, the generic element K of the mesh T h is characterized by the elementwise metric In more detail, the lengths i,K and the orthonormal vectors i,K identify the ellipse, E K , circumscribed to K, which coincides with the push-forward of the circle circumscribed to the reference triangle K via the standard affine map T K ∶K → K (Ern and Guermond 2004). The scalar quantities i,K measure the lengths of the semi-axes of E K , whereas the unit vectors i,K provide the corresponding directions. To quantify the deformation of the element K, it is standard to adopt the aspect ratio s K = 1,K ∕ 2,K ≥ 1 , with s K = 1 for equilateral triangles.
The anisotropic variant of the recovery-based error estimator in Zienkiewicz and Zhu (1987) relies on the correspondence, up to an area-dependent scaling factor, between the H 1 -seminorm of a function v ∈ H 1 ( K ) and the anisotropic counterpart proposed in Formaggia and Perotto (2001), where G K (⋅) is the symmetric positive semi-definite matrix with entries for i, j = 1, 2 , = (z 1 , z 2 ) ∈ [L 2 ( )] 2 . In (6), the first order derivatives defining the H 1 -seminorm are projected along directions i,K 's, to comply with the anisotropic characterization of the element K.
A cross-comparison between (4) and (6) yields the anisotropic error estimator where ∇ = P(∇ h ) − ∇ h denotes the recovered error associated with the material density.
, Following Micheletti and Perotto (2006), estimator is used to generate an anisotropic adapted mesh. This goal is pursued by resorting to an iterative procedure, whose generic -th iteration consists of the following steps: (i) computation of the density function h , solution to problem (1), on the mesh T h ; (ii) evaluation of the a posteriori error estimator (7) for We comment on steps (iii)-(v), separately, by neglecting the iteration index to simplify the notation.
Concerning item (iii), we look for the elementwise metric For the derivation of M we enforce a certain accuracy on h , while minimizing the mesh cardinality and equidistributing the error throughout the mesh. These critera lead to the constrained minimization problem to be solved on each elements K, with ij the Kronecker symbol, and Ĝ K ( ∇ ) = G K ( ∇ )∕| K | . Quantity J K follows from the local estimator in (7), being Notice that the area information is now confined into the term | K | . Due to the error equidistribution, it follows Since minimizing the mesh cardinality is equivalent to maximize the area of the mesh elements, we derive metric M by solving problem (8). In Micheletti and Perotto (2006, Proposition 26), the authors supply the explicit solution to this problem, namely characterizing metric M * K , whereas we have only the ratio of the two lengths * i,K . To compute * 1,K and * 2,K separately, we explicitly impose the error equidistribution, i.e., with a user-defined tolerance on the density h , #T h the mesh cardinality, and where relation | K | = 1,K 2,K |K| has been exploited, with As far as item (iv) is concerned, we commute the elementwise quantities in M into information associated with the mesh vertices, as it is required by standard mesh generators. For this purpose, we resort to an area-weighted average across the patch of elements associated with each vertex Farrell et al. 2011). This step identifies the new metric M .
Finally, metric M is provided as an input to a metricbased mesh generator (see step (v)). For the numerical verification in Sect. 4, we use the software BAMG (bidimensional anisotropic mesh generator), which is the default geometric discretization tool linked to the adopted solver, FreeFEM (Hecht 2012).

Homogenization: direct and inverse techniques
We tackle the design of cellular materials, obtained by the periodic repetition of a unit cell, in order to match a required property at the macroscale. A standard approach is the asymptotic homogenization theory, which is, in general, adopted both in a direct and in an inverse form (Andreassen and Andreasen 2014; Sigmund 1994). The basic idea behind direct homogenization is to inherit at the macroscale the effects associated with the microscale. In this case, the microscale is the known contribution, whereas the macroscale has to be characterized. Vice versa, when dealing with inverse homogenization, we start from a desired feature at the macroscale and we look for the microscopic layout ensuring the expected macroscopic property in the (10) homogenized scale. Some of the methodologies adopted to perform topology optimization at the macroscale have been employed to design the periodic internal structure of cellular materials (see, for instance, Sigmund 1994;Huang et al. 2011;Noël and Duysinx 2017). The high versatility of direct and inverse homogenization justifies the adoption of these techniques for diverse applications (see, e.g., Ivarsson et al. 2020;Bruggi and Corigliano 2019).
In this section, we resort to the asymptotic homogenization theory combined with the SIMP method to design new cellular materials. This leads us to solve an optimization problem which can be formalized as in (1), after a suitable redefinition of the quantities involved.
As a reference physical model, we adopt the linear elasticity framework. In particular, it is convenient to consider the componentwise stress-strain ( − ) relation in terms of the symmetric stiffness tensor E, given by according to the Voigt notation (Helnwein 2001).
Direct homogenization incorporates the contribution at the microscale into the macroscale for coherence model by modifying the stiffness tensor E (Andreassen and Andreasen 2014). To this aim, it is standard to employ an asymptotic representation of the displacement field, following the twostep procedure: (i) we compute the microscopic displacement, * ,ij , with ij ∈ I = {11, 22, 12} , by solving, in the periodic func- for ij, kl ∈ I. The stiffness tensor E H takes into account the effect of the microscale. Thus, when dealing with a homogenized context, the linear elasticity stress-strain relation in (11) is replaced by With reference to an inverse homogenization setting, following Sigmund (1994), we design the optimal unit cell ensuring the desired property at the macroscale by means of the SIMP method properly cast in the homogenization framework. To find the optimal distribution of material m (where the subscript m stands for microscale) in the unit cell, we solve problem (1) where we identify the design domain with Y, and the bilinear and linear forms with respectively, with ∈ U # and ij ∈ I . The design variable m is selected in V # = H 1 ↺ (Y) to extend the periodic conditions on * ,ij to the density function. This assumption simplifies the theoretical and numerical discussion below. Concerning the goal functional in (1) ijkl , for ij, kl ∈ I G , denotes the goal stiffness tensor component adopted by the user to control a physical quantity of interest. Finally, the first box inequality in (1) is reduced to the one-sided constraint on the structure mass with ∈ (0, 1] the maximum allowed volume fraction. To sum up, problem (1) turns into the optimization statement with ij ∈ I G .
As in Sect. 2.2, the minimization in (19) is carried out by using IPOPT. Among the input quantities to IPOPT, the gradient, ∇ m G( m ) , of the goal functional in (16) with respect to the density variable m requires a careful computation. To this aim, we exploit an adjoint-based Lagrangian formulation, by introducing the augmented functional for ij, kl ∈ I G . In particular, ∈ [U # ] g is the adjoint variable used to impose the state equations, with components mn .
By differentiating L with respect to * ,pq , we obtain the so-called dual problems with pq ∈ I G ; the derivative of L with respect to pq yields the g state equations in (19); the gradient of the goal functional for any ∈ V # , is obtained by deriving L with respect to m . Formulas (20) and (21) thanks to the use of SIMPATY algorithm in the design of the new cells, so that m,h ∈ V #,h and * ,ij h and with ij ∈ I G . The resulting algorithm, named microSIMPATY (Ferro et al. 2020b;di Cristofaro et al. 2021), is listed below.
The algorithm switches between an optimization phase (function in line 3) and a mesh adaptation step (function in line 4), based on the anisotropic adaptive strategy detailed in Sect. 2.2. In particular, three tolerances, together with a maximum number of global (optimization + adaptation) iterations, constrain the optimization.
In more detail, controls the stagnation of the number of elements between two consecutive mesh adaptations; sets the accuracy for the minimization problem; is the user-defined accuracy characterizing the equidistribution of the error throughout the mesh elements according to (10)

Generation of a cellular dictionaries
In this section, we provide two different strategies to generate a collection of different unit cells. We will refer to this ensemble as to a dictionary, D = {W 1 , … , W z } , consisting of z different words. The generic word W s = ( s , s , H ( s )) , with s = 1, … , z , gathers the layout s generated by micro-SIMPATY algorithm, the corresponding volume fraction, s , and the associated homogenized stiffness tensor, H ( s ) , defined as in (22). We consider words characterized by different volume fractions, ordered in an ascendant way, and optimized with respect to the same goal quantity G m , even though for different values E G .
We adopt a multi-objective framework, where topology optimization is contemporary driven by several goal quantities. This choice is typical of practical contexts such as the one analyzed in Sect. 6. A first instance of a multi-objective topology optimization carried out by the microSIMPATY algorithm is provided in di Cristofaro et al. (2021).
A multi-objective formulation of problem (19) implies a generalization of the definition for the goal functional G m , here identified by the convex combination of functionals as in (16), i.e., for ij t , kl t ∈ I G , with weights t ≥ 0 for t = 1, … , , such that ∑ t=1 t = 1 . The components of the goal stiffness tensor are collected in the vector G = [E ij t kl t ] t=1 in ℝ , which replaces the scalar input E G ijkl to microSIMPATY algorithm. By exploiting this optimization process, we generate two distinct dictionaries. The words of the first ensemble (referred to as single-cellular dictionary) share the same cell topology; vice versa, the second set (referred to as multi-cellular dictionary) consists of words with different topologies.
Algorithm 2 provides a unified scheme for the generation of the two dictionaries (with =' ' for the singlecellular dictionary, =' ' for the multi-cellular one). The details of the two procedures are addressed in the next sections. (23)

Generation of a single-cellular dictionary
A single-cellular dictionary is characterized by a unique topology. We choose the reference layout, R , yielded by microSIMPATY algorithm for the multi-objective goal functional identified by the reference vector G = G,R ∈ ℝ , and for an assigned volume fraction, v R f , both selected by the user (lines 3-5). In particular, the words in D are generated via a suitable post-processing of R , which essentially increases or reduces the thickness of the considered layout, according to the volume fractions in the input vector . For this purpose, we resort to a parametric design platform, which is based on the Rhino 3D CAD software 1 endowed with the embedded visual programming tool Grasshopper. 2 The definition of each word W s is completed by computing the homogenized stiffness tensor H ( s ) defined as in (22)  We remark that the R is the only topologically optimized layout in the dictionary. However, this is a standard practice in many application fields dealing with graded materials (Wang et al. 2017).
We exemplify this procedure by selecting the multi-objective goal functional Concerning the other input parameters to Algorithm 2, we set with f = [0.35, 0.4 ∶ 0.1 ∶ 0.9] , T h coinciding with a structured 50 × 50 mesh of the unit cell Y = (0, 10) 2 (mm 2 ), and where the flag is set to ' ' since interested in a single-cellular dictionary. The values chosen for the reference volume fraction and the homogenized stiffness tensor components are v R f = 0.5 and E G 1111 = E G,R 1111 = 2.50e−01, E G 2222 = E G,R 2222 = 1.25e−01, respectively. We remark that the components both of the reference and of the homogenized stiffness tensors are scaled with respect to the corresponding full material configuration (i.e., m = 1 in (17)), after setting E Y = 45.00 (N/mm 2 ) and = 0.49 to characterize the base material. Finally, the penalization exponent in (15) is taken equal to 4.
The reference layout, R , yielded by algorithm, is the one in the left panel of Fig. 1. According to the colorbar, red areas ( = 1 ) identify the full material, while the blue zones ( = 0 ) represent the void. The adapted mesh is superimposed to the layout and is highlighted in cyan, in order to emphasize the anisotropic features along the structure boundaries. Figure 1 (right) shows the convergence history characterizing the optimization of the unit cell, in terms of the volume fraction and the stiffness tensor components defining the goal functional. More iterations are required by E H 1111, R and E H 2222, R to stabilize when compared with the volume fraction. The vertical lines identify the adaptation steps. We remark that each mesh adaptation introduces a slight instability in the optimization process (to be ascribed to the variables projection onto the new mesh), which is dumped within few iterations.
Starting from the reference design R and the volume fractions in f , function returns the layouts in Fig. 2. The employment of a unique optimized topology modified via offset criteria may lead to layouts that are not so handy to be manufactured. This is the case of the layouts in panels (a) and (f) which exhibit thin struts and small holes, respectively. Moreover, we have verified that volume fractions less than 0.35 yield unprintable configurations.
Function computes the homogenized stiffness tensor characterizing s 's.

Generation of a multi-cellular dictionary
A multi-cellular dictionary assembles words which are associated with different topologies. As a consequence, we run microSIMPATY algorithm z times (line 16, after selecting a multi-objective goal functional and a certain volume fraction. This leads us to pick the s-th row, for s = 1, … , z , of the tensor G ∈ ℝ z× (line 14), collecting the goal values E G ij t kl t , with t = 1, … , , which define the functional in (23), and the s-th component of the vector f ∈ ℝ z gathering the chosen volume fractions. As far as the testing below, we select f = [0.1 ∶ 0.1 ∶ 0.9], (i.e., z = 9 ), and, for each s = 1, … , 9 , we choose the goal functional in (24) (i.e., = 2 ), with goal values specified in Table 2, panel 2. These targets prescribe the x-direction to be mechanically stiffer with respect to the y-direction, which turns out to be strategical with a view to the design of orthopedic devices tackled in Sect. 6. Also the selected volumes in f are strictly motivated by the application of interest. Moreover, since the volume fraction constraint and the homogenized goal values may be in conflict, inappropriate target inputs might lead to an empty solution space. To overcome this issue, we selected values for the components E G 1111 , E G 2222 so that the optimization turns out to be feasible. As far as the other input parameters to dictioSIMPATY are concerned, we preserve the same values as in (25), while adopting a structured N T 0 h × N T 0 h initial mesh, T h , the flag being ' ' since interested in a multi-cellular dictionary, D . Notice that, to design the nine layouts, we adopt an initial tessellation of the unit cell Y = (0, 10) 2 (mm 2 ) characterized by a different value for N T 0 h (see Table 2, panel 5). The   Figure 3 collects the microSIMPATY layouts s , for the diverse volume fractions, overlapped with the corresponding final anisotropic adapted mesh. The stretched elements of the grids allow to precisely detect the void/material interface, leading to very smooth boundaries, while regularizing the possible sharp angles of the layouts (as a clear example, see the unit cell in panel (d)). Moreover, following Ferro et al. (2020c), to improve the accuracy of the mechanical analysis, we restore an isotropic tessellation in the material portion of the cell, so that the stretched elements act only in sharply detecting the topology.
The obtained configurations respond to the different design prescriptions. In more detail, except for the unit cell in panel (a) which is far from being manufacturable, the x-direction turns out to be the preferential one in terms of mass allocation, at least until the volume fraction is sufficiently far from 1. This evidence is fully in agreement with the imposed mechanical requirement which reduces the stiffness in y-direction. More quantitative information about the final adapted mesh, here denoted by T h to simplify the notation, is provided in Table 2, panel 5, where the cardinality and the maximum aspect ratio of the grid are gathered. We observe that complex topologies require a larger number of elements, whereas the stretching of the mesh elements is strictly related to the size of the holes in the cell layout.
Concerning the homogenized stiffness tensor of the unit cells, we refer to Table 2, panel 3, for the components driving the topology optimization (and to Fig. 9 for the other components). We recognize a good matching between goal and homogenized components, as highligthed by the relative errors, %E 1111 and %E 2222 , in Table 2, panel 4.

Multiscale topology optimization
In this section, we address the design of a macroscale geometry in terms of optimal allocation of microscopic unit cells. This process, widely employed in practice, is generally referred to as multiscale topology optimization.
In accordance with Wang et al. (2017) and Allaire et al. (2019), we detail how to employ a dictionary, D , of unit cells, optimized to match a goal quantity G 1 at the homogenized scale, for the design of a structural configuration in order to satisfy a mechanical requirement, G 2 , at the macroscale.
In the literature, both pre-defined (Allaire et al. 2019;Vigliotti and Pasini 2013;Moussa et al. 2020) and ad hoc designed (Coelho et al. 2008;Nakshatrala et al. 2013;Ferro et al. 2020b;di Cristofaro et al. 2021) cells are employed as words of the dictionary D . We follow the second approach to exploit the dictionaries, D and D , generated in the previous section.
The reference framework is still represented by problem (1) which has to be properly customized according to a multiscale formulation. In particular, we solve the problem where G M = G 2 drives the optimization at the macroscale; the bilinear form of the state equation is  To explicitly define the components of the multiscale homogenized stiffness tensor, H , we follow an approach similar to the one adopted in Wang et al. (2017).
First of all, we build a possible trend for the six sets of values { s , E H ijkl, s } s=1,…,z by means of a polynomial approximation H ijkl ( ) , for ij, kl ∈ I . In particular, the dependence of E H ijkl, s on s is understood due to the implicit relation between the structure volume fraction and the layout.
As a successive step, we make matrix H ( M ) computable in to evaluate the bilinear form in (27). For this purpose, for any ∈ , we resort to the identification = M ( ) , and we assign the value H ijkl ( ) to the component H ijkl ( M ( )) (see Fig. 4 for a schematic representation of the implied correspondence between multiscale density and microscopic volume fraction).
From a computational viewpoint, consistently with the previous sections, we adopt linear finite elements to discretize both the multiscale density and the displacement, i.e., The whole multiscale topology optimization is sketched in Algorithm 3.
The algorithm combines, in a sequential way, the fitting of the data at the microscale (in line 2) with the topology optimization at the macroscale (in line 5). With a view to the optimization process, fitting functions H ijkl , depending on the microscopic volume fraction , are identified with H ijkl as functions of M . Successively, components H ijkl , together with their derivatives, are used by the optimize routine, according to (26), since instrumental to the setting of the state equation (27) and to the gradient-based optimization procedure, respectively. The output of the procedure is the discretized multiscale design variable.

Design of orthopedic devices
This section is devoted to the proposal of an innovative pipeline to design patient-specific insoles for the treatment of foot diseases. The impairments led by foot problems include a higher risk of falling, lower ability to execute standard activities of daily living, and a general lower level of quality of life (López-López et al. 2018).
The investigation carried out in this section focuses on patients affected by diabetes, who suffer from foot ulcers, one of the main causes of limb amputation. Several studies highlighted how mechanical stress can be a major indicator of ulcers progression, since they often appear in the highpressure regions of the foot (Armstrong and Lavery 1998;Ledoux et al. 2005). Due to the mechanical nature of the putative mechanism of ulcers worsening, several studies have investigated the effect of footwear design (Healy et al. 2013). In particular, the literature pointed out that wrong footwear causes over 40% of foot amputation, whereas properly designed diabetic insoles can locally reduce the mechanical stress in the areas of major pressure or in correspondence with ulcers location, thus preventing foot ulceration.
In this framework, our goal is to check whether the proposed multiscale topology optimization can be instrumental to the design of a patient-specific orthopedic device, starting Fig. 4 Multiscale topology optimization: correspondence between multiscale density and microscopic volume fraction from the well-established expertise gained by some of the authors in such an area (Mannisi et al. 2019). However, since the verification setting is here confined to a 2D framework, the results below have to be meant as a proof-of-concept investigation, preliminary to a 3D real design to be manufactured by a 3D printing technology.

Verification
In order to accomplish the design of the patient-specific insole, we resort to Algorithm 3, after introducing some simplifying hypotheses on the physical setting. The main simplification leads us to work on a transverse cut section of the 3D orthopedic device, thus reducing the design domain to a 2D slab in the xy Cartesian plane. Additionally, we assume that the patient-specific configuration is characterized only by three medical conditions for the sole of the foot, i.e., a healthy and an ulcerated zone, together with a possible transition region between these two areas.
With a view to the 3D printing of the optimized configuration, the 2D final layout is meant as to be extruded along the z direction.
From a mathematical viewpoint, the optimization of the insole can be cast into the generic framework in (26). In particular, the target functional G M ( ( M ), M ) coincides with the mass of the orthotic device The two-sided constraint in (26) is generalized to different controls to take into account that the stiffness varies according to the local clinical status of the tissue (healthy, ulcerated, or pre-ulcerated). In more detail, the three inequalities are imposed, where denotes the local static compliance computed on the boundary portion Γ Z ⊂ Γ N , with Z = H (healthy tissue), Z = U (ulcerated tissue) and Z = T (transition area between healthy and ulcerated tissue), scaled with respect to the corresponding compliance, in the full-material configuration. Thus, with reference to formulation (26), it follows that Finally, the PDE constraining the topology optimization is featured by the bilinear form in (27) t h e p r e s s u r e exerted by the heel, the pad, and the toe on the insole. The values are extracted from experimental data (Karia et al. 2016). We analyze two different clinical scenarios, i.e., the case of an abrupt transition between healthy and ulcerated areas (referred to as the H-U setting), and the complete configuration including the three possible status of the tissue (denoted by the H-U-T setting). The two choices essentially differ in the definition of portions Z , for Z = H, U, T (see the bottom panel in Fig. 5 for a sketch), being with The endpoints of the ulcerated and transition zones are derived from the intersection of the boundary portion N with two circular segment areas, in order to model the damaged and the neighboring tissue.
Concerning the adopted discretization, we use linear finite elements to approximate both the design variable, M , and the displacement, , on a 150 × 25 structured computational mesh T h .
The multiscale optimization of the insole will be now driven by the dictionaries, D and D , created in Sects. 4.1 and 4.2, respectively. In particular, the functional in (24) constraining the design of the words of the two dictionaries has been selected in order to minimize the stress component along the vertical direction on the patient's foot. In fact, active mechanical stress transmission to the ulcers, due to standing or walking, may lead to increase the ulcer size over time. As a consequence, it turns out to be beneficial for the patient to decrease deflection at the load application regions for reducing the peaks of pressure in ulceration zones.

Insole optimization based on the single-cellular dictionary
To design the orthopedic device by means of a graded single topology, we resort to Algorithm 3 with the input parameters m , M to be properly set according to the selected H-U or H-U-T configuration, 0 M,h = 0.95 , with the characteristic function associated with domain , and where D is the single-cellular dictionary in Fig. 2 being identically null, and the full material configuration being added to the ones in Fig. 2 for s = 8 . To this aim, we adopt a global polynomial least-square approximation of degree n. The value n = 5 is the lowest degree which provides us a sufficiently reliable data approximation. In particular, Fig. 6 shows the plot of the polynomials H ijkl ( ) , where the circle markers highlight the fitted values for the homogenized tensor components. We also provide the plot of the standard SIMP (i.e., 3 ), which exhibits a very similar trend when compared with the fitting polynomials, with a slightly more evident deviation for H 1111 . This can be justified by the promotion of the x-direction in the stiffness optimization process. Moreover, the four curves ensure tensor H ( M ) to be well-defined as a function of M . Finally, we have numerically checked that tensor H is definite positive for any M , which implies that the bilinear form a H M (⋅, ⋅) in (27) is coercive. The fitting in Fig. 6 becomes instrumental to function of Algorithm 3, since, via (28), it defines the PDE problem constraining the topology optimization in (26).
The lower and the upper bound for the two-sided constraint in (26) are set to and for the H-U and the H-U-T scenario, respectively corresponding to a requirement of high and low (and, if included, medium) stiffness for the insole.
The distribution of the multiscale density function M,h representing the output of Algorithm 3 is shown in Fig. 7 for the patient-specific H-U (top) and H-U-T (bottom) settings. The presence of the transition area between healthy and ulcerated zones induces a milder gradation of the material density in correspondence with the sick and the loaded portions of the boundary N . The total mass of the two insoles is given by 37.8% and 39.4% of the full material configuration for the H-U and the H-U-T scenario, respectively.
In order to grade the different unit cells across the optimized insole, we exploit again the relationship between M,h at the homogenized macroscale and at the microscale. In particular, we resort to a thresholding of M,h to identify the different subregions, s , of the design domain associated with a specific microscopic layout, s . The number of the  [8.00, 15.50, 9.00] T , subregions s will be, at most, equal to z + 1 , after assuming that the full material configuration corresponds to s = z + 1.
According to the thresholding criterion here adopted, we define the subregions for s = 1, … , z , with l,s , u,s ∈ ℝ properly selected by the user according to the density distribution range, and region being z+1 ∈ ℝ + set by the user. We have numerically investigated the effect of thresholding in terms of mass and stiffness of the structure. The percentage variations with and without thresholding are, in general, negligible (we refer to the next section for more quantitative details). By applying this thresholding to the configurations in Fig. 7, we obtain the four subregions i , with i = 1, … , 4 , shown in Fig. 8, after setting l,1 = 0 , u,1 = l,2 = −0.025 , u,2 = l,3 = u,4 = −0.05 , u,3 = l,4 = −0.065 . We observe that the inclusion of the pre-ulcerated portion in the insole design introduces an area of soft material just below the ulcer, and promotes the intermediate densities corresponding to the volume fractions = 0.4 and = 0.5 . These results support not only the reduction of pressure below the ulcer location but confirm the need of a pressure re-distribution around the damaged areas.

Insole optimization based on the multi-cellular dictionary
The topology optimization of the insole is now driven by exploiting the multi-cellular dictionary D generated in  Figure 9 gathers the polynomials H ijkl ( ) , the data to be approximated, together with the reference SIMP behaviour. A cross-comparison with Fig. 6 shows that the four curves exhibit a similar trend, except for polynomial H 1111 which increases slower to 1 when dealing with the multi-cellular dictionary. Moreover, the SIMP curve is now very close to the functions fitting the data for all the components.
The output of Algorithm 3 coincides with the distribution of M,h in Fig. 10, where the optimal design of the insole is distinguished between the H-U and the H-U-T scenarios (top and bottom, respectively). The distribution of the material density is very similar with a slight difference in the total mass of the insole for the H-U and H-U-T settings, equal to 29.8% and 32.4% of the full material domain, respectively. Moreover, the ulcer is taken into account in both the H-U and H-U-T designs in contrast to Fig. 7 (top panel) where the same density characterizes the ulcerated and the healthy portions. Concerning the optimization process, we provide in Fig. 11 the convergence history characterizing the minimization of the mass and the constraints on the local static compliance H , U , and T . The trend is very similar for all the quantities, with an initial transient phase followed by a gradual stagnation towards the optimal values.
With a view to the manufacturing of the optimized insoles in Fig. 10, we apply the thresholding procedure introduced in the previous section. This strategy leads us to identify three subregions, 2 , 3 , 4 , for both the scenarios, H-U and H-U-T, after selecting l,2 = 0 , l,3 = l,4 = u,s = −0.05 for s = 2, 3, 4 . It can be checked that the thresholding leads to slight changes with respect to the mass and the stiffness of the structure. In particular, the mass is essentially invariant (32.4% to be compared with 31.6% without and with thresholding); the compliance characterizing the ulcerated portion practically remains the same, with a 0.07% increment; the compliance in the healthy area, H , exhibits a 2.77% reduction, while T increases more significantly (10.4% with respect to the non-thresholded configuration), promoting stiffer and more flexible materials for the healthy and the pre-ulcerated tissue, respectively. Figure 12 shows the distribution of the three regions together with the associated unit cells. When compared with the single-cellular case (see Fig. 8), we recognize that the multi-cellular dictionary involves fewer unit cells in the prototype of the insole and turns out to be less sensitive to the inclusion of the pre-ulcerated zone in the design of the orthotic device. Finally, the two panels in Fig. 12 emphasize that the presence of the transition area between the healthy and the ulcerated tissue increases the extension of region 4 associated with the volume fraction = 0.4 . The transition areas favour the re-distribution of the pressure, thus avoiding undesired localization of the load. All these considerations lead us to select the setting in Fig. 12 for the validation and the comparison with standard 3D manufacturing choices.

Validation
The optimization procedure set for the design of orthopedic devices is finally tested in practice via a mechanical validation, in the case of the patient-specific H-U-T configuration of Fig. 12. In detail, we computationally explore the macroscopic mechanical behaviour of the prototype via a finite element analysis. Then, we investigate the 3D printing manufacturability of the design, via fused deposition modelling (FDM).

Mechanical performance evaluation
For the finite element analysis, we refer to the 2D geometry in Fig. 13, where the insole is regarded as a porous continuum body and includes a thin layer on the top, to simulate the foot tissue (i.e., fat and skin) in contact with the orthotic device. The considered configuration coincides with the explicit representation of the microcellular material shown in Fig. 12b.
To evaluate the mechanical performance of the optimized insole, we assume as the reference case the standard honeycomb material. In more detail, we mechanically test three insole settings, i.e., (1) honeycomb design with constant volume fraction of 30%; (2) honeycomb design with variable volume fraction, with the same distribution as in Fig. 12b; (3) the optimized H-U-T configuration in Fig. 12b.
The foot tissue layer is assumed linearly elastic, homogeneous and isotropic, as well as the thermoplastic polyether-polyurethane elastomer (TPE) material composing the insole. In particular, Young's modulus and Poisson's ratio are set equal to E TPE = 45 (MPa), TPE = 0.49 , and E foot = 1 (MPa), foot = 0.4 . The loading condition applied to the top tissue layer coincides with the patient-specific load distribution in Fig. 5 (top) while the bottom boundary of the insole is constrained, as shown in Fig. 13. The mechanical analysis for cases (1)-(3) is carried out by using the software COMSOL Multiphysics. 3 Figure 14 compares the deformation of the three designs under the imposed load, in terms of spatial distribution of the displacement magnitude. It is evident that the three insoles exhibit a different mechanical behaviour. In particular, the optimized configuration yielded by the procedure in Sect. 6.1.2 exhibits a maximum value for the displacement which is 2 and 3 times lower with respect to the one characterizing designs (1) and (2), respectively. This is confirmed by Fig. 15 (top) which shows the displacement profile along the cut line highlighted by the black dashed line in Fig. 13. It is evident that the new designed configuration localizes and attenuates the displacement in the ulcer zone, thus preventing the stress concentration in the surrounding area.
In order to characterize the loading-transfer mechanism between insole and foot tissue, in Fig. 15 (bottom) we investigate the mechanical behaviour of the three designs in terms of the von Mises stress distribution along the same cut line as in Fig. 13. The optimized layout in Fig. 12b dampens the stress in the foot tissue and localizes the effect of the loading pressure. In particular, when employing the two honeycomb insoles, the pressure affects also the plantar arch zone (i.e., the region between 0.1 and 0.15 m), though this area is not loaded. This is not the case for the new proposed configuration, thus introducing an innovative and relevant advantage in the design of orthopedic insoles.

Manufacturing feasibility evaluation
To assess the feasibility of the orthopedic device in Fig. 12b in terms of manufacturing, we apply an axial extrusion of 0.01 (m) to the optimized insole in Fig. 13, in order to obtain a section of the 3D model. To this aim, we use the dedicated software of slicing, Simplify3D, 4 to generate a g-code for the 3D printer. The flexible TPE filament Filaflex, with a shore A equal to 82 (Filaflex 82A), has been selected to manufacture the extruded section shown in Fig. 16. This material has proved to be suitable for the 3D printing of insoles, since it has an optimal resistance to elongation and abrasion as well as a high tensile strength (Wang et al. 2020). These properties make Filaflex 82A a perfect material for the production of flexible, comfortable and resistant insoles.
Different slicing parameters and profiles have been tested to obtain an optimal result in terms of reliability and quality of production. After several tests, the main printing parameters are set as: print speed 2200 mm/min, extruder  Fig. 13 with three enlarged views of the different unit cells temperature 235 • C, cooling 60%, top/bottom/outline perimeters equal to 3, no-retraction. The model took approximately 10 h of printing.
These preliminary results confirm the feasibility of the proposed method to create innovative and patient-specific solution for foot problems. The different areas are welldefined (see the enlarged views on areas A, B and C in Fig. 16) and the process is reasonable in terms of printing time and overall quality. This proof-of-concept represents a promising starting point for the realisation of a full model of patient-specific insoles that can be used for prevention of ulceration and can lead to beneficial medical outcomes.

Conclusions
The multiscale topology optimization proposed in this paper provides a fully general, cost-effective tool for the design of macrostructures optimized in terms of local distribution of void, solid, and cellular materials. Among the relevant properties, we list the possibility to decouple the optimization at the macro-and at the microscale, possibly associated with different goal functionals. This feature allows us to control target quantities which are simultaneously pivotal in the design of the final product, in the spirit of a multi-objective optimization. Moreover, the non concurrent optimizations at the micro-and at the macroscale lighten the computational burden of the whole procedure. Indeed, the generation of the cell dictionary could be, a priori, confined to a preliminary phase, mimicking an offline/online paradigm.
A cross-comparison between the single-cellular and the multi-cellular dictionary optimizations leads us to an additional interesting consideration. Actually, the possibility to involve different topologies in the design of the optimized macrostructure results in the adoption of a lower number of different unit cells, i.e., in a structure simpler to be manufactured. Also, the multi-cellular dictionary admits very small volume fractions (for instance, values = 0.2 , 0.3 excluded by D ). Moreover, dictionary D seems to be less sensitive with respect to the analyzed patient-specific configurations. Indeed, the presence of the pre-ulcerated area in Fig. 12 essentially preserves the allocation of the unit cells when compared with the H-U scenario, just slightly increasing the extension of the zones associated with the highest density. This is in contrast to the optimized layouts in Fig. 8 associated with the single-cellular dictionary, which strongly differ in the global unit cells distribution and in the presence of a soft region just below the damaged area in the H-U-T configuration.
The validation carried out in Sect. 6.2 confirms that the developed multiscale topology optimization can lead to some benefits when compared with consolidated medical practices. In particular, the finite element analysis in Fig. 15 shows that the insole design in Fig. 12 for the H-U-T scenario reduces the stress in correspondence with the ulcerated zone and the plantar arch, when compared with standard homogeneous and heterogeneous honeycomb infills. Finally, the manufacturing analysis in Sect. 6.2.2 confirms the feasibility of the optimized prototype with respect to a 3D printing technology.
These preliminary results support the proposed multiscale topology optimization as an effective computational and manufacturing workflow for the design of custom-made insoles for specific foot pathologies.
As a first possible development, we plan to extend the methodology here presented to a 3D setting to provide a more effective solution to different foot pathologies and to comply with patient-specific clinical prescriptions.
In addition, both the optimization phases can be enriched with constraints taking into account the current limits in additive manufacturing (Moussa et al. 2020), as well as the presence of possible failures or geometric defects (Liu et al. 2017) in the final layouts.
A handy management of the transition between different cell layouts represents another crucial issue with a view to manufacturability, in particular when dealing with heterogeneous dictionaries.
Further steps towards a more comprehensive description of complex phenomena include a multi-physics optimization, where different physical problems coexist; the possibility to properly orient the unit cells according to the direction of the applied load; the inclusion of null density in the multiscale design process to concurrently modify the macroscopic topology and the cell allocation.

Replication of results
All the results in the numerical assessment can be replicated in a straightforward way, starting from the algorithms and the data provided in the paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.