{'fraction_non_alphanumeric': 0.037847058497212434, 'fraction_numerical': 0.013674312817241284, 'mean_word_length': 4.908956813898222, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 9, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Modular-topology optimization of structures and mechanisms with free material design and clustering

Topology optimization of modular structures and mechanisms enables balancing the performance of automatically-generated individualized designs, as required by Industry 4.0, with enhanced sustainability by means of component reuse. For optimal modular design, two key questions must be answered: (i) what should the topology of individual modules be like and (ii) how should modules be arranged at the product scale? We address these challenges by proposing a bi-level sequential strategy that combines free material design, clustering techniques, and topology optimization. First, using free material optimization enhanced with post-processing for checkerboard suppression, we determine the distribution of elasticity tensors at the product scale. To extract the sought-after modular arrangement, we partition the obtained elasticity tensors with a novel deterministic clustering algorithm and interpret its outputs within Wang tiling formalism. Finally, we design interiors of individual modules by solving a single-scale topology optimization problem with the design space reduced by modular mapping, conveniently starting from an initial guess provided by free material optimization. We illustrate these developments with three benchmarks first, covering compliance minimization of modular structures, and, for the first time, the design of non-periodic compliant modular mechanisms. Furthermore, we design a set of modules reusable in an inverter and in gripper mechanisms, which ultimately pave the way towards the rational design of modular architectured (meta)materials.


Introduction
This contribution addresses challenges emerging in the design of structures, mechanisms, and architectured materials (archimats) related to two contemporary research topics: (i) automation and individualization of the design and manufacturing processes (aligned with Industry 4.0 precepts), and (ii) environmental sustainability. In relation to the former challenge, topology optimization coupled with additive manufacturing seems to provide the definite answer [1,2]. The sustainability aspect, however, requires not only efficient material usage, but also structural reusability after the end of the structure lifetime [3], repairability [4] and reconfigurability. In general, these aspects have seldom been considered in topology optimization.
Aligned with these goals, we believe that modularity is an ingredient that would allow topology optimization to become more sustainable by forging a pathway, via modules, towards more circular element design, because modular architectured materials, (self-)assembled [5] from a number of patterned components called modules, boost structural sustainability with possibly (re)configurable, multi-functional structures and mechanisms [6,7,8,9,10] while enabling efficient mass manufacturing through prefabrication, e.g., [11,12].
However, optimal design of these modular materials is still in its infancy. In principle, the design problem involves solving a bilevel modular-topology optimization problem incorporating (a) the topological design of a small number of modules, and (b) the spatial arrangement of the modules, see and versatility of the final structure. While the lower level (a) alone introduces only a slight modification to standard topology optimization (TO) procedures [13,9,14], the upper-level assembly plan design (b) adds another combinatorial dimension to the structure of the problem [13]. For the case of continuum topology optimization, the module design problem (a) lacks convexity and appears to be strongly dependent on initial guesses, which makes a rigorous assessment of assembly plan performance challenging. Following a brief introduction to microstructural optimization in Section 1.1, we focus on modular designs and divide the contemporary methods for solving the modular-topology optimization problem into four groups: intuitionbased, multi-material-based and meta-heuristic methods covered in Section 1.2, and clustering-based procedures in Section 1.3. In general, these paradigms enable either a concurrent design scheme, in which the two levels are solved simultaneously and converged designs are (locally) optimal with respect to both the assembly plan and the module topologies, or a sequential scheme, in which the assembly plan design precedes the design of modules, or vice versa.
In this contribution, we develop a novel sequential clustering method for designing modular structures and mechanisms based on Wang tiling formalism, Section 1.4. Finally, the main challenges and our contributions are summarized in Section 1.5.

Microstructural optimization
Microstructural topology optimization has been a very vivid research area [14], particularly due to the link to numerical homogenization [15]. However, conversely to the homogenization, which evaluates the effective properties of periodic media, the inverse homogenization approach proposed by Sigmund [16] has established a new framework for designing periodic materials with prescribed constitutive properties. The local periodicity requirement for homogenization-based methods relies on an infinite scale separation in order to maintain material connectedness of non-uniform microstructures and theoretical, superior performance. Because non-uniform distributions are optimal in general [17,13], but since only finite-scale separations can be practically considered, the disconnectedness issue of neighboring microstructural cells has challenged researchers for years [18,19].
Except for periodic-unit-cell (PUC) designs, which are compatible by definition but usually provide poor structural performance [20,21,10], several conceptual ideas have been proposed to resolve this problem for practical finite-scale separations during the last decade: prescribing connectors of neighboring, often parametric microstructures [22], design post-processing [23,24,18], and maintaining natural mechanical connectivity via single-scale optimization of the whole domain [13,9] or subdomains [19].
Among these methods, the de-homogenization approach [25,24,26,27] has proven to be a superior tool for single-load-case compliance optimization problems. With it, optimal microstructures are known to be rank-2 or rank-3 laminates whose stiffness is aligned with principal stress directions. Consequently, the optimal structures are found by projecting and smoothing optimal parametric microstructures. In more general cases such as multiple loading scenarios, modular structures, or compliant mechanism designs, the topology of optimal microstructures is not known. Although well-performing microstructures can still be obtained for multiple loading scenarios and compliant mechanisms using suitable microstructural parametrizations, such as fiber-reinforced composites [28], these approaches do not extend to modular structures.

Design of modular (meta)materials
Unlike the theoretically-optimal non-uniform distributions of microstructures, modularity limits itself to a finite small number of pattern types. As a consequence, it only approximates the optimal distribution, unless the optimal distribution is periodic. There are two extreme cases of modular designs: PUC microstructures, which repeat a single module in the design domain; and the optimal, infinite number of module types. As a result, modularity drastically reduces the number of design variables compared to the theoretically-optimal solution, but also provides a better structural response than PUC designs [13]. Furthermore, modularity may also improve the susceptibility of structures to local damage or imperfections [10].
As we have already outlined, the grand challenge in the design of modular structures and mechanisms remains a coupled design of module topologies and their connectable spatially-varying arrangements. Hence, early methods for modular design relied on a predefined distribution of modules [29,30,31,32]. Although the optimized designs performed better than PUCs in general, their non-optimality with respect to the macro-scale module-assembly design problem substantially influenced their performance.
An inverse approach has often been adopted in modular metamaterial design. The particular topologies of modules are then usually designed by intuition, combining a few design patterns, and the parametric module properties with their spatial distribution drives the exotic material response [12], such as textured deformation [33], extreme stiffnessto-density ratio [12], tunable Poisson ratio [8], tunable band-gap properites [34], or on-demand acoustic wave filter, guide, lens, or cloak devices [35]. While some of these modular metamaterials rely on lattices [36,12,37,38,39], other are based on origami or kirigami principles [40,41] possibly involving compatible mechanisms [42]. Benefiting from a library of metamaterials already presented in the literature, Wu et al. [34] proposed a structural-evolution-based method to generate meta-atoms of metamaterials with specified properties.
Substituting the intuition-based design with topology optimization, several of the early methods for modulartopology optimization relied on meta-heuristics coupled with numerical optimization. Among these, Tugilimana et al. [21] optimized first the spatially-varying rotations of a truss module, and, ultimately, both rotations and placements of multiple module types [9]. In the latter work, a simulated annealing was used for selecting an appropriate assembly plan, and a non-convex mathematical program for designing module topologies. Another approach was developed by Tyburec et al. [13], who optimized truss modules and their assemblies in the form of Wang tilings. In this method, a genetic algorithm was combined with convex mathematical programming, which provided a reliable way of assessing and comparing the performance of individual assembly plans.
Another class of approaches to modular design relies on concurrent multi-material-based methods. Among these, Zhang et al. [22] developed an approach based on ordered SIMP interpolation. Its main idea encompassed a repeated solution to two subproblems: solving a multi-material distribution problem for fixed material properties, and finding the material properties for a given material distribution. Notice that the former subproblem is not considered in the above meta-heuristic approaches. A similar setting to [22] can also be found in [43], where a rotation of a single parameterized unit cell was optimized for compliant mechanism problems.

Clustering-based methods
The third class of methods relies on clustering algorithms. Because these methods are particularly relevant to us, we review relevant developments in a dedicated section. For the readers' convenience, we categorize the clustering methods according to the quantities that are clustered: densities, the stress or strain tensors, and material stiffness tensors.
Density-based clustering. Density-based clustering determines candidate modular arrangements from clustered pseudodensities of the variable thickness sheet problem. Since these densities usually interpolate isotropic materials, the isotropic homogenized properties of the modules are implicitly assumed. Due to the continuity of the density field, an infinite number of isotropic microstructures is needed in order to describe a general macro-domain. To effectively limit this number, Zhang and Sun [44], Li et al. [45] and Gao et al. [46,47] partitioned the pseudo-density field based on predefined range thresholds or unique values. Each of these clustered densities was subsequently associated with a unique microstructure found, e.g., by an inverse-homogenization-based method.
A related concept was adopted in Yan et al.'s [48] concurrent method, in which the strain energy of thermoelastic lattice structures was minimized and the discreteness of the macro-scale structure assembly was secured by the kmeans algorithm which clustered similar strain energies of the cells without the need to specify fixed thresholds a priori. In another work, Nikravesh et al. [49] partitioned the strain energies by a data-driven approach, with the final designs recovered using a mapping of isotropic microstructures onto the partition elements.
Stress-and strain-based clustering. The previous density-based clustering relying on the assumption of isotropic microstructures is usually too restrictive, since it reduces all microstructural information into a scalar field. Consequently, a more general viewpoint incorporating the directional information, i.e., the data provided by the stress and strain tensors, has emerged.
In one of the earliest works leading towards stress-based clustering, Xu and Cheng [50] investigated the concurrent setting of minimum-compliant multiscale topology optimization. In their work, the requirement of multiple unique microstructures was enforced by using so-called misplaced material volume constraints, which penalized the assignments of microstructures out of their assumed range of principal stress directions. Individual microstructures were designed by the inverse homogenization approach, but lacked mutual connectivity. Another recent concurrent method from Qiu et al. [51] proposed using dynamic clustering to cluster vectors containing the direction of the principal stresses as well as their ratios. The clustering task was solved by the kmeans algorithm. To ensure similar local optima of the clustering algorithm among the iterations of topology optimization, the cluster centers of kmeans were initialized based on the previous run.
In contrast to density-based clustering, stress-or strain-based clustering usually operates without microstructural volume fraction information. This drawback has been circumvented by the combined clustering proposed by Kumar and Suresh [52], who clustered the vectors containing densities and the ratios of principal strains.
Material stiffness clustering. However, even the latter perspective lacks full generality: it designs each microstructure based on the principal directions only, i.e., rank-3 laminates in 3D. Such settings thus cannot handle cases with different optimum microstructures properly. These cases include multiple loading scenarios in minimum-compliance optimization or compliant mechanism design.
The most structurally-efficient designs can only be found when exploiting all available information, i.e., for microstructural design, information following from the fourth-order material stiffness tensor. Considering isotropic materials only, Yang and Li [53] optimized the Poisson ratios of isotropic materials by introducing free isotropic material optimization. Their work combined a hierarchical clustering algorithm with a second-order conic formulation in the transformed polyhedral P-Q space, optimizing the structural complementary strain energy. A generalization, developed by Hu et al. [54], explored the entire design freedom of positive definite material tensors by free material optimization (FMO), e.g. [55]. The optimal spatially-varying material tensors were clustered by a hierarchical clustering and individual micro-structures were obtained using the inverse homogenization method with embedded physics-based connectivity control.

Wang tilings
In this work, we adopt the concept of Wang tiles-a natural generalization of the PUC-as a convenient formalism for describing modular assemblies. Visualized as unit squares with colored edges and a fixed orientation, the concept of Wang tiles was developed by, and is named in honor of, Hao Wang to simulate the ∀∃∀ decision problem of predicate calculus in mathematical logic. In an equivalent domino problem, Wang [56] assumed an arbitrary but fixed finite set of unitary tiles or dominoes and asked whether one can determine that their copies can cover an infinite plane while satisfying the domino rule of matching colors at shared edges of adjacent pieces; see Fig. 2 for an illustration.
An answer to the domino problem flourished in the work of Wang's student Berger [57], who proved that, in addition to tilesets that can cover the infinite plane by repeating periodic patterns, aperiodic tilesets also exist which satisfy the tiling objective yet forbid any periodic arrangement. Hence, the domino problem is not decidable. Since then, numerous aperiodic tilesets have emerged; see [58,Chapter 11] for a thorough introduction.
Berger's original proof relied on a reduction of the Turing machine halting problem into the domino problem. This reduction implies that Wang tiles are Turing-complete, and hence their (self)-assembly can simulate desired Turing computations. This property has motivated applications, e.g., at the DNA level [59,60], and, recently, by Jílek et al. [5] on the centimeter-scale. Visually appealing, nonperiodic yet compressed arrangements have impelled the adaptation of Wang tiles in computer graphics to generate naturally-looking textures [61,62] and Poisson disk distributions [61,63]. Sharing objectives, these applications have inspired the use of Wang tiles in modeling materials microstructures, both for geometrical representation [64,65,66] and calculations [67,68,69,70]. Because the cardinality of the Wang tileset balances the computation efficiency of the PUC with precise (micro-)structural description, Wang tilings also appear to be useful in securing automatic connectivity in the topology optimization of modular truss structures [13].

Aims and novelty
As follows from the literature survey, the design of modular architectured materials has only recently become an active research domain in the topology optimization field. Although several significant results have already appeared, there still remain fundamental questions to be answered.
Because the overall response of modular structures follows from two constituents, with the smooth macro-domain response driven by the effective constitutive parameters that are impaired by the presence of heterogeneities [70] in the form of module topologies, finding an appropriate assembly plan is the major challenge that, when successfully addressed, would allow attaining the desired structural effectiveness.
In our recent work [13], we introduced a concurrent modular-topology optimization method for designing minimumcompliant modular truss structures encoded with the corner Wang tiling formalism. The method relied on a bilevel optimization involving convex conic optimization subproblems (leading to the optimal compliance for a fixed assembly plan), and a genetic algorithm to optimize the module arrangement. The convexity of the module design subproblem is indeed an essential assumption; otherwise, the assembly plan performance cannot be assessed reliably. Unfortunately, for the case of continuum topology optimization, the module design problem lacks convexity in general, so the method in [13] lacks applicability.
Thus, in this contribution, we rely on the sequential clustering approach, recall Section 1.3, operating over material stiffnesses, and a standard single-scale topology optimization [71] extended to structural modularity. Each of our modules is defined as a quadruple of edge labels, thus defining a Wang tile, recall Fig. 2. In this dual point of view, we ask how do we place at most m different interface types to make the resulting modular structure behave efficiently?
The answer is indeed in the form of a valid Wang tiling over at most m colors.
To make our setting fully general, our clustering dataset originates from a solution to an FMO problem, which we extend to handle both minimum-compliance and compliant mechanism design objectives, Section 2.1.1. To the best of our knowledge, the procedures for designing modular assembly plans in linear elasticity problems have mostly been limited to minimum-compliance optimization. The only exception we are aware of is [43], optimizing rotations of a single unit cell for compliant mechanisms via multi-material topology optimization. However, a more general setting has already been requested [10].
Surprisingly, as we will show in Section 2.1.2, the optimal elastic properties obtained from the FMO problem often contain checkerboard patterns. Nevertheless, this artifact has not yet been investigated with connection to the FMO problems, although it has appeared in published results, e.g., [54]. Here, we mitigate the issue by enhancing the kinematic response though a refined design-element mesh.
The objective of the clustering algorithms is to partition a given dataset of points into a specific number of clusters. However, standard clustering algorithms have several disadvantages: they are usually stochastic and, therefore, require multiple evaluations to generate quality outputs. Moreover, they also generally fail to maintain structural symmetries. To overcome these shortcomings, we have developed a novel greedy agglomerative clustering algorithm that is deterministic while preserving structural symmetry of the elasticity field, when present (Section 2.1.3). Using this clustering algorithm on a module interface-based discretization produces a valid Wang tiling.
For appropriate assembly plans obtained via the clustering procedure, the module topologies must also be optimized. Unfortunately, the emergent optimization problems appear to be strongly dependent on the provided starting point. To make the initial choice problem-specific yet well-defined, we propose to use a solution to a macro-scale FMO problem with modularity to construct initial guesses based on the traces of the optimal stiffness tensors obtained from FMO (Section 2.2.1). A single-scale topology optimization extended to structural modularity is finally used for a computation of the optimized modular topologies in Section 2.2.2.
We illustrate the developed theory on four selected examples, including minimum-compliance optimization of the MBB beam problem, Section 3.1, and three compliant mechanisms: inverter (Section 3.2), gripper (Section 3.3), and a module set reusable among both the inverter and gripper mechanisms (Section 3.4). We conclude that our approach to designing modular structures provides fairly efficient and well-behaving designs.
In this contribution, we adopt the following notation: scalar variables are written in plain, lower-case italic letters (a). First-and second-order tensors appear in lower-case, bold italic (a) and bold regular (a) letters, respectively, and the fourth-order tensors are denoted as bold capital regular letters (A). Finally, column vectors and matrices appear in sans-serif lower and upper-case bold characters (a, A).

Methods
This section is devoted to establishing a sequential method that provides a heuristic solution to the modulartopology optimization problem with the modular assembly plan parameterized as valid Wang tilings, recall Section 1.4. Our approach to designing modular structures and mechanisms encompasses four sequential steps: (i) optimum free material design performed on a module interface-based mesh, Section 2.1.1, including a proper handling of emerging checkerboard patterns (Section 2.1.2), (ii) and clustering of the optimal elasticity matrices by a newly-developed agglomerative clustering algorithm to output valid Wang tilings, Section 2.1.3. For the resulting assembly plan, we (iii) optimize the material properties of individual modules via an FMO extended to modularity, Section 2.2.1, which provides an appropriate initial guess for (iv) the standard continuum topology optimization extended to modularity in Section 2.2.2.

Assembly plan design
Based on the inherent bilevel nature of the modular-topology optimization problem, we split the above four steps into two parts-a heuristic procedure for assembly plan design, and then topology optimization of the modules. In this section, we aim to answer the question of how should we design efficient assembly plans? To this goal, we introduce free material optimization in Section 2.1.1, handling of checkerboard patterns in Section 2.1.2, and a novel clustering procedure in Section 2.1.3.

Free material optimization
Let us consider a two-dimensional continuum body Ω under appropriate static and kinematic boundary conditions. In this domain, the free material optimization (FMO) method seeks the spatially-varying material elasticity tensors E that make the structural response optimal with respect to a specified performance functional a(u) [55] of a displacement field u.
In this manuscript, we assume a linear strain-displacement relation where ε denotes a strain tensor. The sought elastic material properties are represented by the effective elasticity tensor E that satisfies linear Hooke's law written in indicial notation as in which σ is a stress tensor.
In what follows, we adopt the (engineering) Voigt notation that allows us to rewrite the tensors ε, σ, and E as column and square matrices, respectively. To this goal, let be the strain column matrix, with γ 12 being the engineering shear strain, and let denote the stress matrix. Then, the matrix form of the Hooke law (2) reads as with the material stiffness matrix constructed from the (to-be-optimized) components E i jkl of the forth-order elasticity tensor E. To simplify the further developments, we introduce a vector form of the fourth-order tensor E as such that where • 2 denotes the standard Euclidean norm of the vector •, and • F is the Frobenius norm of the forth-order tensor.
Our approach to a numerical solution of the FMO optimization problems relies on approximate solutions to the governing equation obtained with the finite element method: we search for the kinematically admissible displacement field u ∈ R n dof satisfying the equilibrium equations where f ∈ R n dof denotes the column matrix of equivalent nodal forces, n dof stands for the number of degrees of freedom, the super-script • ( ) expresses an assignment of • to the element , and K E (1) , . . . , E (n ) ∈ S n dof 0 is the symmetric, positive definite structural stiffness matrix-a function of the discretized field of elasticity matrices E ( ) ∈ S 3 0 of n finite elements. In this notation, S • refers to a space of symmetric square matrices of the size • × •, and denotes positive definiteness. In addition, we define X to express that the eigenvalues of X are greater than or equal to . The structural stiffness matrix K follows from an assembly of element stiffness matrices K ( ) based on the element gather matrices L ( ) and a design-independent term K 0 ∈ S n dof 0 , where the element contributions K ( ) E ( ) ∈ S 0 evaluate as In (11), B ( ) denotes the strain-displacement matrix arising from the adopted displacement approximation, and Ω constitutes the element domain. From now on, we assume E ( ) to be constant over each element, so that K ( ) is a linear function of E ( ) .
In the spirit of standard topology optimization (TO) [71], FMO seeks the most efficient distribution of a limited amount of material, constrained by Eqs. (12c)-(12e), that satisfies the (elastic) equilibrium, Eq. (12b). However, compared to TO, the FMO setting opens up a broader design freedom, because the full material stiffness matrix is the subject of optimization. Here, we formalize the FMO problems of our interest as In (12), s ∈ R is a slack variable and a ∈ R n dof stands for a fixed vector representing the discretized performance of the functional a. Furthermore, ∈ R >0 and ∈ R >0 are lower-and upper-bounds on the (sum of the) eigenvalues of E ( ) , ∀ ∈ {1, . . . , n }, V ∈ 0, 1 stands for the volume fraction, v ( ) ∈ R >0 represents the volume of the -th element, and Tr is the trace operator, amounting to the sum of its argument's eigenvalues. Notice that in contrast to established formulations [55,72,54], we use an element volume-based weighted average in (12c) to facilitate optimization with non-structured spatial discretizations.
When setting a = f, the slack variable s becomes an upper bound on the compliance functional, (12b), rendering (12) to be a compliance minimization problem. Such optimization problem is convex: the functions in (12a) and (12c)-(12e) are linear (matrix) functions, and (12b) possesses a positive definite Hessian, see Appendix B.2, and it is therefore also convex. For the case of a f, the second term in function (12b) amounts to a T u. Consequently, formulation (12) facilitates optimization of compliant mechanisms by minimizing weighted sums of the displacement vector entries.
Considering Eqs. (12c) and (12e), a natural question arises of how to choose and . Clearly, the lower bound for eigenvalues must be sufficiently small to mimic void materials but should be large enough to avoid problems with a possibly non-existing inverse in (12b). In our case, we fix = 10 −3 . To determine the upper-bound on the sum of eigenvalues, , we link the FMO problem (12) with our ultimate goal of designing microstructures made of a mixture of an isotropic material and void, optimized with standard continuum topology optimization. Assuming a family of isotropic materials under plane stress that are parameterized by the Young modulus E ( ) ∈ 0, 1 and a fixed Poisson ratio ν, the trace of the isotropic material stiffness matrix E ( ) iso amounts to If the entire design domain is filled with the solid material, V = 1, the trace constraint (12c) implies that With the theoretical bounds on the Poisson ratio ν ∈ (−1, 0.5 , Eq. (14) implies that for all solid isotropic materials under plane stress ∈ 5 4 + √ 6 2 , ∞). The optimization problem (12) can readily be solved by nonlinear semidefinite programming solvers such as Penlab [73] or Pennon [74]. However, for the special case of a = f, one can also adopt any of the more widely available (convex) linear semidefinite programming optimizers. In this contribution, we rely on the Pennon solver kindly provided by its authors. For compliance minimization problems, we use both the first-and second-order derivatives, see Appendix B. For compliant mechanism problems, it appeared to be more efficient to use an identity matrix instead of the exact Hessian term from Eq. (46). In addition, practical computations revealed that it is more beneficial to solve a scaled version of (12), we refer the reader to Appendix C for details.

Checkerboard pattern
The FMO generates (locally-)optimal elasticity matrices E ( ) that are, however, prone to the checkerboard pattern phenomenon [75], especially in the case of the bilinear quadrilateral finite elements, which are the most widely used elements in the TO community [76]. In Fig. 3, we show this emerging issue on traces of the optimal material stiffnesses in compliance optimization of a bar under uniaxial compression. Please note that while this issue is apparent at the element stiffness traces in this case, it may be difficult to recognize in general because FMO operates directly over stiffness components and, hence, it has no scalar unknown field as in standard density-based TO, and the checkerboard may be lost as a result of visualization.
There are probably multiple reasons why the checkerboard issue occurs. First, we optimize over element design variables, whereas the finite element system (9) operates over a nodal-based displacement field. Therefore, the output of (12) shall be interpreted as optimal discretized (nodal) stiffnesses in the form of the stiffness matrix K E (1) , . . . , E (n ) , which seems to be checkerboard-free. However, discretized nodal stiffnesses result from a (weighted) sum of the stiffnesses of connected elements, recall (10). Clearly, the element stiffnesses do not play any role unless they change a value in their assembled sum K E (1) , . . . , E (n ) . Hence, the element stiffnesses may indeed form a checkerboard pattern, which is often falsely recognized in finite element analyses as structurally more efficient [75]. We note here that associating the design variables with the nodes does not resolve the problem, because the element stiffnesses (11) are usually evaluated at the Gauss points inside the finite element interior, not at the nodal points. Furthermore, the bilinear quadrilateral elements used in this manuscript are known to be prone to locking, and, therefore, they tend to overestimate the shear effects [77,78]. Again, this may result in a load transfer that is nonphysical.
To our knowledge, the checkerboard issue has not been investigated with connection to the FMO problem (12) yet, probably because only regular orthogonal mesh discretizations have been adopted, and, in that case, the issue is often not visually significant but indeed occurs, see e.g., [54,Figs. 4b and 4e]. For another discretizations, such as the rotated one used in this manuscript, Fig. 3, the checkerboard is more apparent.
Traditionally, two approaches have been adopted in the density-based topology optimization to mitigate the checkerboard patterns: the sensitivity and density filters [71]. These filters suppress the checkerboard pattern by a convolution of the first-order derivatives, or of the pseudo-density field. They can also be applied to the FMO problem directly by filtering the entries of the elasticity matrices independently because the space of positive semi-definite matrices is convex.
Another option is to enhance the kinematic response by adopting higher-order finite elements, or by a refined design-element discretization, where the original mesh is used for the material design variables and a refined variant for the evaluation of structural response [79]. Here, we follow the latter approach and adopt the refined discretization shown in Fig. 3b. Figure 4: Illustration of the proposed heuristic approach to generate efficient valid Wang tilings: First, we (a) discretize a design domain into modules, and (b) associate a finite element with each edge: elements belonging to horizontal and vertical edges are shown in gray and white, respectively. For a refined interface-based discretization, we (c) solve the free material optimization problem (12), represented in this figure by the traces of material stiffnesses. Subsequently, the material stiffnesses visualized in the (d) horizontal and (e) vertical edges as envelopes of the stress components σ x , σ y and τ generated by arbitrary unit strain vectors, are clustered independently based on the Frobenius norm, establishing (f) a Wang tiling assembly plan. The assembly plan forms the input to (g) topology optimization of the structure consisting of (h) a set of modules.

Clustering for modular Wang tile design
Using the optimal, checkerboard-free material stiffness field, our next goal is to determine structurally-efficient modular assemblies in the form of valid Wang tilings. In this work, we follow the direction of clustering the material stiffnesses because it represents the most general option: it generalizes both the density clustering for isotropic materials with a fixed Poisson ratio and the strain-based clustering for problems possessing optimal design in the form of rank-3 laminates. Our approach thus extends to the remaining linear elasticity problems including multi-load designs and compliant mechanisms. To this goal, this section develops a greedy agglomerative clustering algorithm that maintains structural symmetries when present.
As discussed already in Section 1.4, valid Wang tilings are determined uniquely by the edge color codes. Therefore, instead of the modules, Fig. 4a, we must cluster material stiffnesses corresponding to the horizontal and vertical edges. These stiffnesses are obtained by solving the free material optimization problem (12) on a rotated, module interface-based finite element mesh, see Fig. 4b. When comparing with Fig. 4a, it clearly follows that each of the modules is defined by a quadruple of (edge-related) elements in Fig. 4b. Hence, after clustering, the edge labels automatically define Wang tiles and a valid tiling, see Figs. 4f-h. We note here that a similar approach relying on rotated mesh has been applied to generate Wang tile-based textures in computer graphics [61] and microstructure modeling [65].
The common aim of the clustering algorithms is to partition a set of points P e (1) , . . . , e (n ) into (at most) m groups, such that a specified distance-based criterion is minimized [80]. Among these algorithms, kmeans is probably the most widely used. It starts with a random assignment of points to clusters, computes a centroid of each cluster based on a specified distance function, and re-assigns each point to the cluster with the nearest centroid. This procedure repeats until convergence, but the resulting clustering strongly depends on the initial partitioning.
The second group of clustering algorithms, spectral clustering, exploits a representation of the points P in a graph form. For this representation, one computes the graph Laplacian matrix and finds a partitioning based on its eigenvectors. Based on our testing (not included), neither kmeans, nor spectral clustering led to convincing designs.
In this manuscript, we rely on the third class of algorithms, the hierarchical clustering, because this class may avoid randomization and can be very simply tailored to satisfy additional constraints imposed on the clustering data.
In particular, the latter feature appears to be very difficult to achieve with the two previous classes of algorithms. Hierarchical clustering algorithms can be categorized into divisive clustering-starting with all points belonging to a single cluster, partitioned recursively-and agglomerative clustering-where each point initially defines a cluster, and these are joined repetitively based on a linkage criterion.
In our implementation, we follow the agglomerative clustering approach and adopt the squared Frobenius norm to measure the distances of the material stiffness tensors, which is equivalent to the Euclidean norm of the difference in the stiffness vectors e ( ) , recall Section 2.1.1. To compute a distance between clusters, we represent each cluster C i by its centroid with components where w ∈ R n >0 are weights of the input data and |C i | denotes the cardinality of points in C i , and by the cluster cost To determine which clusters shall be merged, we compute for all pairs of clusters (i, j) the deterioration cost where c i, j stands for the cost of the cluster C i ∪ C j . We accept the merger that least deteriorates the objective function. Notice that this linkage criterion is not equivalent to merging the nearest clusters, as more distant clusters containing fewer points may be cheaper to merge than the nearer ones with more data points. In this sense, we believe that this clustering procedure provides a more balanced number of points inside clusters. Assume now that the material stiffness distribution exhibits symmetry with respect to a vertical or horizontal axis. For such a setting, we find it reasonable that the resulting clusters preserve this symmetry pattern. In our implementation, this can be guaranteed by restricting which clusters cannot join, i.e., the deterioration cost (17) is set to ∞ for some non-admissible pairs (i, j).
To this goal, we split the material stiffnesses into three groups α, α and β. β contains material stiffnesses of the elements that are either located on the axis of symmetry or do not have symmetric counterparts, and groups α and α consists of the symmetric complements. More precisely, let us assume a fixed i ∈ {1, . . . , |C|} and j ∈ {1, . . . , |C|} \ i such that the entries of these centroids satisfy t i,1 = t j,1 , t i,2 = t j,2 , t i,3 = t j,3 , t i,4 = −t j,4 , t i,5 = −t j,5 , and t i,6 = t j,6 , which certifies symmetry of the corresponding material stiffness centroids with respect to the vertical or horizontal axis, cf. Appendix A for details. If no such j exists for a fixed i or when t i,4 = t i,5 = 0, we set t i ∈ β. If the corresponding pair exists, we set t i ∈ α if t i,4 > 0 or t i,4 = 0 ∧ t i,5 > 0, and t i ∈ α otherwise.
Using these groups, a finite deterioration cost is associated only with (i) a merger of the corresponding pair i ∈ α and j ∈ α , producing a cluster belonging to β, (ii) joining clusters {i, j : i, j ∈ β, i j}, with the result also belonging to β, and (iii) a merger of {i, j : i, j ∈ α, i j}, which requires concurrently merging the corresponding pair Figure 5: Illustration of all allowed cluster aggregations to generate Wang tilings preserving the symmetric distribution of material stiffness matrices when present. The set β consists of material stiffnesses of the elements located on the axis of symmetry and those that do not have symmetric counterparts, and α with α containing the symmetric complements.
{i, j : i, j ∈ α , i j}. In this last case, the number of clusters is decreased by two, with the merged outputs belonging to α and α , respectively. All the possible aggregations are illustrated in Fig. 5.
In each iteration, the aggregation procedure reduces |C| by 1 at least, building an aggregation tree called dendrogram. When |C| ≤ m, the clustering algorithm terminates, see Algorithm 1 for the pseudo-code. As an output, we obtain an assignment of the input data (horizontal or vertical edges) into clusters, Fig. 4e-f, and this partitioning is uniquely representable by a valid Wang tiling over (at most) m colors, Fig. 4g α, α , β, C ← mergeCheapest(pair 1 , cost 1 , pair 2 , cost 2 , pairs 3 , cost 3 ) 10: end if 11: end while 12: return C

Design of topologies of modules
Once an assembly plan-in the form of valid Wang tiling-is obtained, the dimension of design space in TO usually reduces significantly because n m < n , with n m being the number of Wang tiles or modules and n denoting the number of module positions in Fig. 4a. However, it seems that-at least for compliant mechanism design problemsthe landscape of the objective function used in TO is very complex and optimized topologies depend strongly on a supplied initial guess (see Section 3.5 for an example). This issue becomes even more relevant when dealing with the design of modules reusable among multiple structures (Section 3.4), since an inappropriate initial guess may prevent the mechanisms from functioning at all.
In this section, we first recall a variable linking procedure [13] to describe modular structures mathematically, and apply it to the FMO problem (12). The modular FMO allows for an inexpensive computation of the initial guesses for TO based on the element traces, see Section 2.2.1, and also for a simple comparison of the efficiency of different assembly plans. Description of the adopted TO method is given in Section 2.2.2.

Modular free material optimization
Let n m denote the number of modules and {g ∈ N n , ∀ ∈ {1, . . . , n } : g ≤ n m } be a group vector assigning a module number m ∈ {1, . . . , n m } to each element ∈ {1, . . . , n } of the module scale, recall Fig. 4a. In addition, let us define a group matrix G ∈ B n ×n m as Using this definition, we can construct the full-domain elasticity properties e k , k ∈ {1, . . . , 6}, by mapping the grouped modular elastic properties e g,k = e 1 k , . . . , e n m k T using the group matrix G, e k = Ge g,k , ∀k ∈ {1, . . . , 6}, with the notation • m expressing an association of • with the module m. Notice that an association with an element remains denoted by • ( ) .
Then, the optimization problem (12) extended to structural modularity reads as with the stiffness matrix assembled as Because of n m < n , the optimization problem (20) is usually solved much faster than (12). After its solution, we obtain the optimal elasticity matrices of individual modules, and their traces are used to build an initial guess for TO as constant pseudo-density fields of individual modules according to

Topology optimization
The above-described procedures based on the free material optimization yield an efficient assembly plan and related optimal material stiffnesses that are constant over individual modules. Ultimately, though, we want to obtain explicit design of modules manufactured from a single material. As mentioned already in Section 1, one option is to perform inverse homogenization for each module separately, optionally with enforcing geometrical or mechanical continuity across their boundaries [19,54]. Here, we opted for a different approach in which the functional continuity is resolved automatically by performing the single-scale topology optimization with a modularity projection following from the optimized assembly plan.
In particular, we adopted the widely-used Solid Isotropic Material with Penalization (SIMP) approach [81,71]. The design is thus parameterized by a pseudo-density field 0 ≤ ρ(x) ≤ 1 that, in turn, dictates the distribution of the material stiffness according to where E 0 0 and E min 0 are the stiffness tensors of the bulk material and void space (which attain small but positive eigenvalues in order to avoid a singular Hessian), respectively, and the penalization parameter p = 3 is introduced to promote "black-white" designs [71]. Throughout all the simulations, we assumed an isotropic material under plane stress conditions with fixed Poisson ratio ν = 0.3; the pseudo-density ρ thus effectively scales only the Young modulus.
Similarly to the free material optimization problem (12), the formulation of the topology optimization builds on the finite element discretization of the governing equation of linear elasticity with unknown vector ρ containing elementwise constant densities ρ ( ) . For notation simplicity, we assume that each of n m modules is discretized with n m finite elements. The modularity is then introduced in the formulation via a boolean group matrix G TO ∈ {0, 1} n ×(n m ·n m ) which maps the module-related densities collected in vector ρ g = ρ 1 , . . . , ρ n m T onto full-domain densities ρ, recall Eq. (19). Keeping the same structure as in Eq. (12), the modular topology optimization is eventually formalized as min s TO , ρ=G TO ρ g s TO (25a) The structural stiffness matrix K (ρ) follows from Eq. (10) (including the design-independent part K 0 ) with E ( ) depending on ρ ( ) through Eq. (23). The slack variable s TO is introduced here only formally to denote the objective. Contrary to the formulation (12), the volume constraint (25c) is here formulated directly in element densities ρ ( ) and no analogue to the eigenvalue constraints (12d) and (12e) is needed because they are already satisfied due to the interpolation scheme (23). On the other hand, this comes at the price of a restricted design space compared to the FMO in which the whole material stiffness matrix is unknown. Having effectively only volume and bound constraints, we solve the problem (25) with the first-order Optimality Criteria (OC) method [71] due to its simplicity and practical performance. Moreover, we modify the standard OC scheme with the gray-scale suppression proposed by Groenwold and Etman [82] to further promote faster convergence towards designs with clearly distinguished phases.
The sensitivity of the objective function ∂s TO ∂ρ g under satisfied equality equation (25b) follows from the chain rule where ∂s TO ∂ρ is obtained similarly to FMO, see Appendix B.1, and ∂ρ ∂ρ g = G TO , recall Eq. (24). The same chain rule applies also to the sensitivity of the volume constraint (25b), denoted as ∂V ∂ρ g further on. In order to avoid meshdependency and checkerboard patterns discussed already in Section 2.1.2 for FMO, we opt for the sensitivity filter [83] that modifies the objective sensitivity by weight averaging. In the discrete setup, where each ρ ( ) belongs to the center of gravity x ( ) of the corresponding element , the filtered non-modular sensitivity ∂s TO ∂ρ is defined as which leads to a linear relation with the filtering matrix F w following from Eq. (27). Here, we assume a linear weight function with the radius r.
Denoting the sensitivity vectors g s TO ∂V ∂ρ for terseness, the update of the design variables ρ g within the modified OC scheme takes the form where η and δ are the damping factor and the move limit of the original OC method, respectively, q is the power coefficient pertinent to the power-law gray-scale suppression, and µ facilitates the inconsistent modification [71] such that µ is set to −∞ for the compliance minimization and 10 −10 for mechanisms. Finally, Λ denotes the Lagrange multiplier related to the volume constraint (25b), whose optimal value is found with a bisection algorithm; see e.g. [71] for more details.
The entire implementation and related dataset are available at [84]. While the free-material optimization and clustering procedures described above were implemented in Matlab, relying on the Pennon optimizer [74], the topology optimization part was implemented in an in-house, object-oriented C++ code, which was linked against Intel ® oneAPI Math Kernel Library and utilized Pardiso direct sparse solver.

Results
We illustrate the presented theoretical developments with three selected classical problems in topology optimization. In particular, we show the results for the minimum-compliance optimization of the Messerschmitt-Bölkow-Blohm (MBB) beam in Section 3.1, Sections 3.2 and 3.3 then cover the design of compliant mechanisms (an inverter and a gripper, respectively), and a novel reusable modular gripper-inverter design is demonstrated in Section 3.4.
For all these problems, we limited the relative volume with V = 0.4. Modules were discretized with a regular grid of 100 × 100 bilinear fully integrated quadrilateral finite elements. We assumed isotropic, linear-elastic materials under the plane stress with E = 10 −9 and ν = 0.3 for the voids and E = 1 and ν = 0.3 for the bulk material, which in turn lead to 0 ≤ Tr E ( ) ≤ with the upper-bound on the sum of eigenvalues = 235/91, recall Eq. (14). In the topology optimization part, radius r of the sensitivity filter was set to r = 3.5h with h being the length of an element edge. Optimality Criteria's parameters were as follows: move limit δ = 0.1, damping coefficient η = 0.5 for compliance and η = 0.3 for mechanisms. We used a continuation strategy for the power coefficient q of the gray-scale suppression; starting with the 20 th iteration, the initial value q = 1.0 was iteratively updated as q new = min(1.01q, 2.0).
Within all examples, the maximal number of iterations was set to 150. The convergence criteria were based on the change in the objective, |δs TO | < 10 −12 , and the change in the design variables, δρ g ∞ < 10 −2 . To avoid premature convergence with significant gray areas, quantified with the measure [83] we added a convergence requirement M nd < 10 −3 . All topology optimization problems were initialized using the module-wise constant guess provided by FMO, as specified in Eq. (22).

Messerschmitt-Bölkow-Blohm beam
First, we examine the performance of our scheme for the Messerschmitt-Bölkow-Blohm (MBB) beam problem. The beam is 0.375 in height and 1.0 in width, and discretized by 32×12 equisized square modules as shown in Fig. 6a. The beam is supported along 1/128-wide segments at the bottom corners of the beam and loaded by a distributed load with the intensity 64 acting at a 1/64-wide segment at the top. We start with discretizing the MBB beam with a rotated mesh, recall Section 2.1.3, and solve the FMO problem (12) that provides us with an optimal compliance of s FMO = 6.11, see Fig. 6b. Using the optimal material stiffnesses, we proceed with the clustering procedure explained in Section 2.1.3 for a hierarchical merging of the varying material stiffness tensors into a predefined number of groups, one to four in our case. This clustering procedure is applied twice: once for horizontal and once for vertical edges. In the two-dimensional setting, this process can easily be visualized as follows: the stiffness matrices E ( ) are a linear mapping of the engineering strain vector ε ( ) onto the stress vector σ ( ) , Eq. (2). Hence, a sphere of all unit engineering strains is mapped by the elasticity matrix onto a 3-dimensional ellipsoid whose surface corresponds to all admissible stresses induced in the material by arbitrary unit-sized strains, determined using the 2-norm. The clustering process partitions the ellipsoid of all elements in FMO into a given number of bins, trying to replace similar ellipsoids with an approximate one corresponding to the cluster centroid. Consequently, the lengths and orientations of the ellipsoid's axes characterize the material stiffness. Partitioning induced by the clustering is denoted by coloring of these ellipsoids, see first two columns in Fig. 7.
After clustering, each module position encompasses a quadruple of edge labels, defining therefore a Wang tile, and these tiles are placed in an assembly plan satisfying the Wang tiling formalism, recall Section 1.4 and the developments in [61]. For such modular assemblies, we solve a reduced-sized modular-FMO problem (20) to receive effective

Horizontal clustering
Vertical clustering Optimized structure Wang module set Figure 7: Modular designs of the MBB beam problem for the numbers of edge colors ranging from 1 to 4. The first two columns illustrate the partitioning produced by the hierarchical clustering algorithm, the third and fourth columns visualize optimized modular topologies and their assembly. stiffnesses of individual modules, and use their traces to construct initial guesses for modular TO, recall Section 2.2.1. Although the TO algorithm described in Section 2.2.2 usually converges to a neighborhood of a locally-optimal solution, we believe that it is reasonable to assume that the performance of all modular designs shall be bounded by a PUC solution from above, and by a non-modular design from below. Please note, however, that despite the fact that we have not observed any violations of this assumption, it is not guaranteed in general due to the inherent non-convexity of the optimization problems. Using the parameters of TO listed at the beginning of the section, we discretize the interior of each module type using a regular 100 × 100 mesh, implying that the MBB beam contains 3.84 millions bilinear quadrilateral finite elements. Using the filter radius of 3.5 elements, standard TO converges to the design in Fig. 6c of compliance s TO,∞ = 20.69. As expected, imposing modularity degrades structural performance: the worstcase single-module, i.e., PUC, design obtained for a single-color clustering possesses the compliance s TO,1 = 56.91, see the first row in Fig. 7. Using two colors, we obtain a design in the second row of Fig. 7 composed of 12 modules, of an improved compliance of s TO,2 = 49.15. The three-color-based clustering produces a design with 29 modules, exhibiting the compliance s TO,3 = 31.24, see the third row in Fig. 7. Finally, the four-color clustering provides us with the design in the fourth row of Fig. 7, composed of 46 modules and compliance s TO,4 = 29.09, which is 41% more compliant than the non-modular design. We wish to emphasize that, because our hierarchical clustering algorithm respects the symmetric stiffness distribution, the resulting module sets contain multiple module pairs that have symmetric topologies.
In addition to structural performance, we have also investigated the influence of modularity on stress distribution. For the MBB beam problem, the maximum von Mises stress in optimized modular designs exceed the maximum stress value in the non-modular design, see Tab. 1, which is consistent with our findings in optimization of modular truss structures [13].

Force inverter
Our second illustrative problem involves designing a modular force inverter mechanism. We consider the design domain to be 1.0 both in width and height and to consist of 10 × 10 square-sized modules, as shown in Fig. 8a

Horizontal clustering
Vertical clustering Optimized structure Wang module set Figure 9: Modular designs of the force inverter problem for the numbers of edge colors ranging from 1 to 4. The first two columns illustrate the partitioning produced by the hierarchical clustering algorithm, the third and fourth columns visualize optimized modular topologies and their assembly.
to this point of loading and another, output horizontal spring of the stiffness k out = 0.05 is placed at the mid-height of the right edge of the domain. For these boundary conditions, we search for a material distribution that makes the unit input force minimize the motion of the output point to the right, denoted by u out in Fig. 8a. When solving the FMO formulation (12) on a rotated mesh we receive s FMO = −0.892, see Fig. 8b.
The clustering procedure, illustrated within the first two columns in Fig. 9, produced assembly plans shown in the third column of the same figure. For these assembly plans, we solved the modular-FMO problem (20) to receive the initial guesses for TO. Similarly to the MBB beam design, we may provide performance bounds of topologicallyoptimized designs established by a non-modular design, s TO,∞ = −0.398 in Fig. 8c, and a single color clustering PUC, s TO,1 = 0.000, see the top row of Fig. 9 showing that the PUC design does not connect the unitary input force to the output properly. However, based on our tests, we believe that such a PUC inverter design does not exist for the selected combination of materials used in TO and for the given modular discretization. Please note that although a PUC force inverter was presented in [10] for a symmetric part of the design domain, its extension to the full domain would require a second module type in our formalism.
As shown in Fig. 9, the two-color-based clustering provides a design composed of 11 modules with the perfor- Table 2: Overview of the modular inverter performances, where u in denotes the horizontal displacement at the force input, s TO the mechanism performance, and σ v,max the maximum von Mises stress. Values marked by † express a non-functional mechanism. Further, GA abbreviates geometrical advantage, evaluated as s TO /u in . Design

Gripper
As the third problem, we consider the traditional gripper mechanism design [85]. Similarly to the force inverter, we adopt a square design domain of the size 1.0 × 1.0 and split it into 10 × 10 equisized square modules, see Fig. 10a. Similarly to the force inverter and [85], we fix the top and bottom 0.025-long segments of the left edge of the domain.
At the center of the left edge, we place a spring with the stiffness k in = 1 and an input force of a magnitude 1, both aligned with the horizontal direction. In the central part of the opposite, right edge, we prescribe an empty space of 2×4 modules to facilitate a placement of an item that shall be grasped by the mechanism. The grasping occurs through fixed jaws placed at the top and bottom of the empty space. These jaws are 0.025 in height and 0.2 in width and made of a solid isotropic material with E = 1 and ν = 0.3, matching so the material adopted in TO and the maximum trace constraint in FMO (12d).
For the mechanism functionality, the input force must trigger closing of the gripper's jaws, which we quantify by minimizing the difference of the displacements u o2 − u u1 of the rightmost points of the jaws, each equipped with a stiffness k out = 0.025 in the vertical direction.
Again, we first solve the FMO optimization problem (12)

Horizontal clustering
Vertical clustering Optimized structure Wang module set Figure 11: Modular designs of the gripper problem for the numbers of edge colors ranging from 1 to 4. The first two columns illustrate the partitioning produced by the hierarchical clustering algorithm, the third and fourth columns visualize optimized modular topologies and their assembly.
design. Finally, consistently with the findings for the inverter mechanism, modularity in this case also decreases the maximum von Mises stress in the structure, see Tab. 3.

Module reusability in gripper and inverter
To demonstrate the merits of modularity, we address a reusable simultaneous design of the previously described inverter and gripper mechanisms as our last example. Based on the multi-functionality of the resulting module designs, it becomes apparent that modularity paves a way from structure-based designs to meta-material based ones. Table 3: Overview of the modular gripper performances, where u in denotes the horizontal displacement at the force input, s grip the mechanism performance, and σ v,max the maximum von Mises stress. Further, GA abbreviates geometrical advantage, evaluated as s TO /u in . For the reusability design problem, we adopt the discretizations and boundary conditions described in Sections 3.2 and 3.3. However, the reusable designs add another level of complexity as they allow for prioritizing of individual subproblems, i.e., they generate a multi-objective optimization problem. In what follows, we reduce this setting to a single-objective one by introducing weights w inv and w grip , and write the resulting objective as with s inv and s grip denoting the objective values of the inverter and gripper problems, respectively. Based on a com-Horizontal clustering Vertical clustering Optimized structure Wang module set Figure 13: Modular designs of reusable inverter-gripper problem for the numbers of edge colors ranging from 1 to 4. While the first two columns illustrate the partitioning produced by the hierarchical clustering algorithm, the third and fourth columns visualize optimized modular topologies and their assembly. Table 4: Overview of the modular reusable mechanism performances, where u in,inv and u grip denote the horizontal displacement at the force input of the inverter and gripper mechanisms, respectively, and s TO,inv with s TO,grip the corresponding objective function of individual mechanisms, s TO the combined one, and σ v,max the maximum von Mises stress. Values preceded by † symbol denote a non-functional mechanism. Further, GA abbreviates geometrical advantage, evaluated as the sum of the mechanism performances s TO to the sum of mechanisms inputs. parison of the objectives of modular gripper and inverter designs and numerical tests (not included), we believe that the choice of w inv = 1.0 and w grip = 0.60 is reasonable, because it balances the terms w inv s inv and w grip s grip for the non-modular designs, recall Sections 3.2 and 3.3. The first step in all previous examples was a solution to the FMO problem on the rotated FE mesh; however, we have already solved these subproblems in Sections 3.2 and 3.3 and, hence, we can recycle them here. Thus, instead of another solution, we cluster the former material distributions concurrently, Figs. 8b and 10b. In the clustering procedure, recall Section 2.1.3, we weight the stiffness matrices of the gripper mechanism by w grip and of inverter by w inv . Similarly to TO, these weights maintain more balanced performance. Finally, we reach the assembly plans and module sets shown in Fig. 13.
For TO, we aggregate the independent volume constraints and adopt V = 0.4 for the union of the domains. Then, the optimized non-modular designs in Fig. 12

Dependence on initial guess
For the case of the 4-color-based reusable modular mechanisms, we further illustrate the strong dependence of the TO part on supplied starting points. First, because the agglomerative clustering algorithm is deterministic, an increase in the number of module edge colors always generates assembly plans for which the starting points obtained for a lower number of edge colors are also feasible. As a simple example, we have initiated the 4-color-based TO with the starting point obtained from the solution to the 3-color-based modular FMO. The resulting design shown in Fig. 16a exhibits a slightly improved overall objective function s TO,4 = −0.561 compared to the results in Section 3.4; with the gripper improved significantly (s TO,4,grip = −0.585), and the performance of the inverter mechanism degraded (s TO,4,inv = −0.210). Interestingly, the maximum von Mises stress is also reduced, recall Tab. 4.
Another straightforward option is to adopt the common starting point of a uniform mass distribution, yielding s TO,4 = −0.353 and the designs in Fig. 16b. Here, the gripper exhibits s TO,4,grip = −0.589, but the inverter mechanism fails to work, showing that an appropriate choice of the starting point is very important and that it is difficult to rigorously compare the performance of different fixed assembly plans.

Conclusions
In this contribution, we have developed a novel heuristic sequential framework for designing modular structures and mechanisms based on Wang tiling formalism. To the best of our knowledge, this manuscript presents the first study which allows modular compliant mechanisms with non-uniform modules to be designed, thus answering one of the challenges outlined in [13,10].
Our strategy relies on a solution to the Free Material Optimization (FMO) problem (12) on a module interfacerelated finite element mesh, and clustering of the resulting elasticity tensors. Topologies of individual modules in the assembly plans generated by the previous steps are eventually found by standard Topology Optimization (TO) extended to structural modularity, which we initialize using a well-defined yet problem-specific initial guess. Compared to our prior study on the modular-topology optimization of truss structures [13], this framework constitutes a computationally more attractive scheme that requires neither a repetitive solution to the module topology optimization problem nor its convexity. Consequently, the new approach also applies to non-convex problems arising, e.g., in continuum topology optimization.
As the first step of our method, we solved the FMO problem (12) on a module-interface-based mesh discretization: solving for the (locally-)optimal elasticity tensors in the discretized domain using the Pennon [74] optimizer. Interestingly, the tensors are often prone to the checkerboard phenomenon, at least when the most common bilinear quadrilateral finite elements are employed. Although the resulting degradation of structural responses is very significant and can be verified simply by projecting the optimal material properties onto a refined mesh, the checkerboard issue has surprisingly never been studied or even reported with connection to FMO problems. The reason behind this can be attributed to the fact that the checkerboard patterns in tensorial field are hard to identify in general-contrary to, e.g., the scalar density fields in TO. However, the checkerboard issue in FMO is often discernible when visualizing traces of the optimized material stiffnesses, recall Figs. 3 or cf. [54,Figs. 4b and 4e]. Here, we eliminate the checkerboard issue by adopting a refined design-element discretization, in a manner similar to [79].
The particular shape of the finite element discretization used in the problem (12) was adopted to provide a simple interface-based definition of individual modules, i.e., each module being defined by four edges and the corresponding quadruple of elasticity matrices. This definition, which is equivalent to Wang tiles, allows us to control the number of interface types among the modules. In this manuscript, we have forbidden rotations of modules, so the clustering of the vertical and horizontal edges was performed independently.
To this goal, we developed a deterministic agglomerative clustering algorithm tailored to satisfy structural symmetries when present. While in very specific instances of linear elasticity one may rely on partitioning simpler datasets than with stiffness tensors, recall Section 1.3, our clustering covers the most general class of problems, and therefore also applies to multiple-load case minimum-compliance optimization or compliant mechanism design problems.
The resulting clustered designs form modular assembly plans in the form of valid Wang tilings, and their performance can readily be assessed using FMO formulation extended to structural modularity (20). However, the modular FMO problem also appears appealing from another viewpoint: its outputs in the form of the traces of module stiffnesses can be advantageously supplied to TO as initial guesses. Subsequently, TO-with the Optimality Criteria method with sensitivity filtering and a continuation-was used to arrive at optimized topologies.
We demonstrated the applicability of the developed framework on three classical problems: the design of a minimum-compliant modular MBB beam, of a force inverter, and of a modular gripper mechanism. Ultimately, we introduced a new problem for designing topologies and assembly of modules that can be reused in both the gripper and inverter mechanisms.
For all these problems, we have shown that with an increasing number of colors, the number of module types not only grows, but overall structural performance also improves. In the ultimate limit, we would reach the performance of a non-modular design. A second interesting observation emerged from the von Mises stress distributions. For the case of compliance minimization problems, modularity seems to increase the maximum stress value compared to non-modular designs, but the converse is true for compliant mechanisms-few strong singularities in compliant hinges in non-modular designs are usually spread into a larger number of weaker ones, reducing the overall maximum stress level.
This work can be extended in several directions, starting with a straightforward extension to three-dimensional problems. Assuming that simpler manufacturability is one beneficial future output of modularity, we will try to develop constraints to secure a continuous material distribution inside individual modules and adopt a robust scheme to suppress the compliant hinges in mechanism designs [10]. Another promising extension constitutes the design of multi-functional mechanisms under finite strains. Such extension, however, would require extending the clustering procedure in Section 2.1.3 to deal with (pseudo)time-dependent datasets. In addition, modularity can also be exploited within optimization and in finite element analyses by efficient preconditioners or domain decomposition methods [86]. Finally, to address the challenge of manufacturability, we plan to manufacture the selected modular designs and validate experimentally the structural response of the prototypes against numerical model's predictions.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
From the first three test strain fields it follows that the diagonal entries of the elastic material stiffness matrix must match, i.e., e (i) 1111 = e ( j) 1111 , e (i) 2222 = e ( j) 2222 , and e (i) 1212 = e ( j) 1212 . Using this observation, the fourth test strain yields e (i) 1122 = e ( j) 1122 . The relation of the two remaining components of the elastic tensor are determined from the stress balance induced by the fifth field. Then, e (i) 1112 = −e ( j) 1112 and e (i) 2212 = −e ( j) 2212 .

B. Derivation of sensitivities
Let us assume an unknown field denoted by y, which encompasses the elements of E in the case of (12) and ρ in (25). In addition, we define a (y) a T K(y) −1 f − s, (34) so that the inequalities (12b) or (25b) are equivalent to a(y) ≤ 0.
In the following subsections, we derive analytical expressions for the first-and second-order derivatives of a(y) with respect to y.
Let us now define matrices U, Λ with the i-th columns computed as Λ i (y) = ∂K(y) ∂y i λ.

C. Scaling of the optimization problem
In the case of f a, numerical solution to (12) appears to be sensitive to scaling of the design variables. To make the solution process more robust, we substitute E ( ) with Ẽ ( ) , which implies that Tr(Ẽ ( ) ) ∈ 3 / , 1 . Consequently, we arrive at the modified formulation min s sc ,Ẽ (1) ,...,Ẽ (n ) s sc (47a)