An FE-DMN method for the multiscale analysis of fiber reinforced plastic components

In this work, we propose a fully coupled multiscale strategy for components made from short fiber reinforced composites, where each Gauss point of the macroscopic finite element model is equipped with a deep material network (DMN) which covers the different fiber orientation states varying within the component. These DMNs need to be identified by linear elastic precomputations on representative volume elements, and serve as high-fidelity surrogates for full-field simulations on microstructures with inelastic constituents. We discuss how to extend direct DMNs to account for varying fiber orientation, and propose a simplified sampling strategy which significantly speeds up the training process. To enable concurrent multiscale simulations, evaluating the DMNs efficiently is crucial. We discuss dedicated techniques for exploiting sparsity and high-performance linear algebra modules, and demonstrate the power of the proposed approach on an industrial-scale three-dimensional component. Indeed, the DMN is capable of accelerating two-scale simulations significantly, providing possible speed-ups of several magnitudes.


Problem setting
Injection molded short fiber reinforced components are frequently used for industrial applications as they combine favorable mechanical properties, free formability and short cycle times. As a result of the injection molding process, the fiber orientation and the fiber volume fraction may vary continuously within the component. Characterizing all possible orientation states for such materials is an arduous and expensive task, both experimentally and by simulative means. To illustrate the complexity to be handled routinely, Fig. 1(a) shows a quadcopter frame arm made from short fiber reinforced polyamide with a local fiber orientation determined by an injection molding simulation, see Section 5 below for details. The color scale encodes the local fiber orientation tensor. Magenta corresponds to a unidirectional, cyan to an isotropic and yellow to a planar isotropic fiber orientation state, see Fig. 3. We observe that, for the majority of the quadcopter arm, the fibers are almost aligned. Moreover, we encounter isotropic and planar isotropic fiber orientations in areas where   Fig. 1(b) represents an example for a component with a spatially varying complex microstructure. A monolithic finite element (FE) simulation of such a structure which resolves the heterogeneities is not feasible with the current computational power. If the microscopic heterogeneities fluctuate on a scale that is much smaller than the size of the component, homogenization methods may be used for obtaining so-called effective material models which account for the physical mechanisms of complex and highly heterogeneous microstructures. Effective models emerge naturally by solving a partial differential equation, the cell problem, on a suitable microstructure. For linear material models, the effective properties may be precomputed and cached. In this way, it is possible to compute the mechanical response of microstructured components.

State of the art
Treating inelastic materials is more difficult, as the internal variables live naturally on the microstructure, and cannot be "homogenized" to the macroscopic scale, in general. FE 2 methods, introduced by Renard-Marmonier [2] and subsequently refined [3][4][5][6], offer a solution by furnishing each Gauss point of the macroscopic finite element model with a microstructure on which the cell problem is solved, accounting for the evolution of the internal variables on the microscopic scale. To speed up the solution process on the microscale, FFT-based approaches [7,8] may be used, giving rise to the FE-FFT [9][10][11] method in the concurrent multiscale context. Despite recent progress in computational efficiency and parallel computing strategies, concurrent multiscale methods with full-field models on the microscopic scale are typically too computationally demanding for industrial use.
To mitigate the computational burden of concurrent multiscale methods, the microscale problem is considered as a parametric partial differential equation which needs to be solved repeatedly (for slightly different input parameters), and model order reduction techniques may be utilized. Motivated by classical meanfield methods [12,13], Dvorak and co-workers [14][15][16] introduced the transformation field analysis (TFA).
The TFA applies to small strain (visco-)plastic material models and assumes the inelastic strains to be piece-wise uniform on specific subdomains, accounting for the resulting elastic deformations via strain localization tensors. In this way, effective models with a finite number of internal variables arise, see also Chaboche et al. [17]. Furthermore, Liu and co-workers [18][19][20] introduced the self-consistent clustering analysis (SCA), which is similar in spirit to the TFA, but exploits the Hashin-Shtrikman variational principle [21][22][23] and is not per se limited to small strain (visco-)plasticity. Nevertheless, when considered as discretization methods, both TFA and SCA show a slow convergence rate in terms of the number of clusters [24] which is rooted in the weak approximation capabilities of piecewise uniform functions [25]. Inspired by recent variational estimates for nonlinear materials [26], the non-uniform transformation field analysis [27] (NTFA) relies upon non-uniform inelastic basis functions, permitting the approximation errors of the fields to be made as small as desired. However, the difficulty is transferred to prescribing suitable evolution equations [28,29] for the reduced inelastic strains, which should be independent of the basis, cheap to evaluate and consistent upon refinement. Possible remedies are Taylor series expansions of the force potential [30][31][32], mixed variational principles [33,34] or dedicated "reducible" models [35], giving rise to the FE 2R (R as reduced) method [36].
As an alternative to methods which approximate the solution of micromechanical problems on unit cells, it is possible to approximate the effective properties directly. Typically, this comes at the cost of operating on a high-dimensional domain of interest. Data driven methods, especially artificial neural networks (ANNs), are predestined for such approximation tasks. There is a number of works considering training ANNs to approximate the effective elastic energy of a medium, see Yvonnet and coworkers [37][38][39]. Furthermore, the regularity of the effective stress facilitates the direct approximation of the stress-strain relationship of inelastic problems, see the works of Jadid [40], Penumadu-Zhao [41] or Srinivasu et al. [42] for different approaches. Motivated by natural language processing, recurrent neural networks (RNN) may provide a framework for incorporating history dependence into the approximation of the stress-strain relationship. For instance, Gorji and coworkers [43] demonstrated that RNN are able to capture effects such as the Bauschinger effect, permanent softening or latent hardening in the context of elastoplasticity.
Using artificial neural networks accompanied by an on-the-fly switching to a reduced order model in a two-scale simulation was investigated by Fritzen et al. [44]. A shortcoming of such data driven methods appears to be that their predictive quality appears low far away from the training domain, and accounting Recently, Liu et al. [49] proposed a transfer learning approach to treat fiber reinforced composites by DMNs. More precisely, they proposed to interpolate the parameters of already identified DMNs, corresponding to different fiber orientation states. In this work, we go beyond this a posteriori approach, and seek DMNs covering the entire spectrum of fiber orientations arising in such a fiber reinforced component.

Contribution and outline
In this article, we investigate a multiscale methodology for direct deep material networks, which covers all possible variations of the second-order fiber orientation tensors and permits concurrent multiscale simulations with DMNs at the Gauss point level. Similar to the FE 2 , FE-FFT and the FE 2R methods, we call the ensuing method the FE-DMN method.
To account for a spatially varying fiber orientation, see Section 2, we augment direct DMNs by the fiber orientation interpolation concept introduced by Köbler et al. [50]. In contrast to Liu et al. [49], we do not consider a transfer learning strategy, i.e., to interpolate already identified models, but propose an a priori interpolation strategy. More precisely, we investigate deep material networks which explicitly depend on the fiber orientation, and identify the optimal model parameters jointly.
For this purpose, we utilize high-fidelity microstructures of fiber reinforced composites [51] which permit us to routinely cover the possible fiber orientations at industrial filler fraction and fiber aspect ratio. To this end, we sample the linear elastic training data from up to 31 microstructure realizations. We also improve upon the previous sampling strategy [45,46] for the constituents' stiffness tensors used in the offline training. We show by example that it is sufficient to cover those stiffnesses which arise as possible material tangents on the microscopic scale. For component scale simulations of industrial problems, it is necessary to optimize the user-defined material models based upon the identified DMNs, see Section 3. We take special care to ensure that the interpolated DMN generalizes accurately to the inelastic regime.
To this end, 78 additional microstructure realizations were generated, exclusively for the inelastic validations. We compute the stress response of each of the 109 generated microstructure representations by an FFT-based computational homogenization [7,8] code and compare the former to the predicted stress response of the interpolated direct DMN. The validation results show that with a maximum error of 5.5%, the DMN is capable of predicting the effective stress of all investigated 109 discrete fiber orientation states sufficiently. We refer to Section 4 for details.
To show that our approach is applicable to state of the art engineering computations, we consider the entire process chain of a quadcopter frame, see Fig. 1(b) and Section 5. We conduct an injection molding simulation, map the spatially varying fiber orientations upon a finite element mesh and conduct a twoscale simulation using the identified DMN surrogate model. We implement the DMN as a user-material (UMAT) subroutine in ABAQUS only relying upon the provided software interfaces. The quadcopter consists of four arms made from injection molded short fiber reinforced polyamide, two base plates which are made from aluminum and four legs which are made from pure polyamide. The simulation model of the quadcopter consist of more than two million elements resulting in almost ten million degrees of freedom. In more than 1.9 million elements, a deep material network is integrated implicitly at every Gauss points, accounting for the local microstructure information in the simulation. Last but not least, we discuss the computational costs accompanied by our approach in Section 6, and demonstrate that the educated guess of the potential computational power of DMNs made in the conclusion Liu et al. [49] was too pessimistic.

Direct deep material networks for variable fiber orientation
In this section, we extend direct deep material networks to short fiber reinforced composite microstructures parameterized by the second order fiber orientation tensor. The main technical tool is the fiber orientation interpolation technique, introduced by Köbler et al. [50], which we apply on the node level of the deep material network. First, we recall the basics of direct deep material networks and the fiber orientation triangle. Subsequently, we combine both concepts.

Direct deep material networks
Let GSM denote the set of all generalized standard materials (GSM) [52]. Then, any two-phase periodic microstructure Y ⊆ R d in d spatial dimensions gives rise to the (nonlinear) homogenization function which maps two input GSMs to the effective GSM that emerges by solving the cell problem of first order homogenization, see Gajek et al. [47] for details. Homogenization functions may be regarded as the basic objects of studying micromechanics at small strains. In general, evaluating homogenization functions requires significant computational resources. Only for special microstructures, the evaluation can be performed with minimal effort. An example for such a microstructure is given by a two-phase laminate, uniquely characterized by a direction n of lamination and the volume fractions c 1 and c 2 of the two phases. A direct deep material network is defined as a hierarchy of such two-phase laminates, see Fig. 2. The concept was originally introduced by Liu and co-workers [45,46] for laminates with fixed direction of lamination, but an additional rotation layer, and simplified by Gajek et al. [47]. By combining laminates in a hierarchical manner, the resulting homogenization function may be rather complex and, by a judicious choice of the involved laminates, may be used as an approximation of the homogenization function (2.1) corresponding to the original microstructure Y , which is significantly less demanding to evaluate. On a more formal level, a direct DMN is a perfect, ordered, rooted binary tree of depth K, where a two-phase laminate B i k is assigned to each node of the tree. We reserve the letter k for labeling the depth of a node, whereas the horizontal index is consistently indexed by the letter i. Our layer count only comprises the laminate layers, and the input is counted separately. Thus, the DMN comprises 2 K − 1 laminate nodes. For a two-phase DMN of depth K, the homogenization function is defined recursively by traversing the binary tree from the leaves, at level K, to the root Input materials are assigned in an alternating fashion, i.e., holds. We refer to Fig. 2(a) for a schematic. Liu et al. [45,46] noticed that parameterizing the involved laminates by the volume fractions c 1 and c 2 is not optimal. Indeed, if one of the volume fractions is zero, the entire corresponding sub-tree will have no further influence. It is more convenient to parameterize the laminates' volume fractions by assigning pairs of weights w 2i−1 K+1 and w 2i K+1 to each laminate on the input level K. These weights should be non-negative and sum to unity. Then, by traversing the binary tree from the leaves to the root, the weights on the k-th level are computed by the sum of weights of the respective laminates on the previous level, i.e., holds, see Fig. 2(b). The volume fractions c 1 and c 2 of each laminate B i k are then computed by normalization For fixed tree topology, a direct DMN is uniquely determined by the directions of lamination, one for each laminate, and the weights of the input layer. These free parameters are identified based on linear elastic precomputations, the so-called offline training, see Section 3.1. Once the free parameters are identified, the DMN can be applied to nonlinear and inelastic materials. This online evaluation is described in Section 3.2.

The fiber orientation triangle
Suppose a fiber orientation state is given in terms of a fiber orientation distribution (FOD) function ρ, which specifies the probability to find fibers in direction p. Advani-Tucker [53] introduced the second order fiber orientation tensor [53] as a compact measure for the current fiber orientation state. Despite its limited information content, its compact form makes it the typical quantity of interest for commercial injection molding simulations [54].
Higher-order moments of the FOD are then estimated by closure approximations, see Montgomery-Smith et al. [55].
The tensor A 2 is symmetric and positive definite with unit trace. Consequently, only five independent parameters are involved. In terms of an eigenvalue decomposition where the matrix Q ∈ SO(3) is orthogonal and the eigenvalues λ 1 ≥ λ 2 ≥ λ 3 are sorted in a descending order, the fiber orientation tensor A 2 may be described by two parameters λ 1 and λ 2 which satisfy the inequalities Thus, up to an orthogonal transformation, every tensor A 2 corresponds to a point in the triangle described by the inequalities (2.10), see Fig. 3. In this article, we follow Köbler et al. [50] and use the CMYK coloring scheme for encoding different fiber orientations as shown in Fig. 3. The manufacturing process induces local variations of the fiber orientation of short fiber reinforced thermoplastic components. Thus, for component-scale simulations, these variations need to be accounted for by the material models. If the fiber orientation state is described in terms of the second-order fiber orientation tensor A 2 , a family of effective material models, one for each such tensor A 2 , needs to be supplemented. By general covariance considerations, two fiber orientation states which differ only by an orthogonal transformation should give rise to effective material responses which differ only by this orthogonal transformation. Consequently, if the considered fiber orientation states are parameterized by the second order fiber orientation tensor, the essentially different fiber orientation states will be parameterized by the points inside the fiber orientation triangle (2.10). Thus, the material models to be identified are parameterized by a two-dimensional continuum. Furthermore, a certain continuity of the effective response of the material model, considered as a function of the fiber orientation tensor, is expected. Indeed, changing the fiber orientation only slightly is expected to change the effective mechanical response only slightly, as well, at least for non-critical loading. Unfortunately, the typical multiscale approach based on representative volume elements is unable to leverage this continuity. Indeed, although the effective material response depends continuously on the fiber orientation tensor, the representative volume element does not. Indeed, due to the stochastic nature of such fiber-filled volume elements, infinitely many different representative volume elements may be used to give rise to the same effective response. In particular, this reasoning has the following implication. Suppose that we furnish each point (λ 1 , λ 2 ) in the fiber orientation triangle (2.10) with a corresponding representative volume element Y λ1λ2 . Even if all these elements have the same size, the function (λ 1 , λ 2 ) → Y λ1λ2 will not be continuous in any useful way.
As an alternative, Köbler et al. [50] proposed a fiber orientation interpolation procedure on the level of effective stresses. This idea avoids the difficulty of interpolating internal variables which live in dif-ferent locations for different microstructures. However, this approach comes at a price: the number of stress evaluations is tripled by this approach. Indeed, for any fiber orientation state, the stress response associated to the three corners of the interpolating triangle needs to be evaluated. More precisely, we consider DMNs which are parameterized by points (λ 1 , λ 2 ) inside the fiber orientation triangle. As the DMN's weights are directly linked to the constituent volume fractions of the underlying microstructure [45][46][47], we will seek weights which are even independent of the fiber orientation. This appears reasonable, and we refer to Milton [56] and Torquato [57,Sec. 20.2.2] for background material.

Fiber orientation interpolation of deep material networks
In order to interpolate the lamination directions n i k on the fiber orientation triangle, we parameterize each normal n i k ∈ S 2 by spherical coordinates Then, we interpolate the angles α i k and β i k on the fiber orientation triangle (2.10) by a global (finite element) shape function. Please note the difference to Köbler et al. [50], who rely upon linear interpolation on a subtriangulation of the fiber orientation triangle. Particularly compact expressions for the finite element shape functions are obtained by transforming the parameters λ 1 and λ 2 to barycentric coordinates ϕ 1 , ϕ 2 and ϕ 3 , i.e., via solving the linear system see e.g., Vince [58]. We collect the parameters of the polynomial shape functions in a vector φ = [φ 1 , . . . , φ M ], where M denotes the number of shape functions. Then, the interpolated angles may be expressed as in terms of the parameter vectors p = [p 1 , . . . , p M ] ∈ R M and q = [q 1 , . . . , q M ] ∈ R M . In this article, we investigate linear, tri-linear and quadratic shape functions, see Tab. 1.
To sum up, extending DMNs to account for varying fiber orientation reduces to increasing the number of unknown parameters. Indeed, instead of identifying the angles α i k and β i k , the parameter vectors p i k and q i k are sought, in addition to the unknown weights w i K+1 .

Implementation
In this section, we explain how to evaluate linear and nonlinear homogenization functions of an interpolated direct DMN efficiently.

Offline training
The goal of the offline training is to identify the free parameters of the DMN, namely the weights and the vectors p i k and q i k used for interpolating the angles (2.13), which we collect in "long" vectors We insert the parameters of laminates on level K first and add the parameters of laminates for decreasing level index in their corresponding order. We enforce the non-negativity constraint on the weights represent the DMN's linear elastic homogenization function in the form DMN L Λ maps the input stiffnesses C 1 and C 2 , the fiber orientation parameters λ 1 and λ 2 and the unknown fitting parameters p, q and v to the DMN's effective stiffness. The specific binary tree structures of the DMN can be exploited to evaluate DMN L Λ efficiently. To this end, we assign the input stiffnesses C 1 and C 2 to the laminates of level K in an alternating fashion. The directions of lamination are interpolated on the fiber orientation triangle based on the given parameters λ 1 , λ 2 , p and q, see Section 2.3. The input stiffnesses are homogenized at level K for each laminate independently. In the next step, the homogenized stiffnesses serve as the input for level K − 1 and so forth, until the root is reached, giving rise to the DMN's effective stiffness C. The process of propagating stiffnesses from the K-th level to the root is illustrated in Fig. 4.  Figure 4: Schematic illustration of the stiffness propagation in a two-phase DMN of depth four Thus, computing the effective stiffness of a DMN reduces to computing a sequence of effective stiffnesses of two-phase laminates. According to Section 9.5 in Milton's book [59], the linear elastic homogenization function of a laminate may be determined by solving the equation , the set of symmetric d × d matrices, and P : Sym(d) → Sym(d) stands for a projection operator, which depends on the direction of lamination n i k , and reads (P(n)) mnop = 1 2 (n m δ no n p + n n δ mo n p + n m δ np n o + n n δ mp n o ) − n m n n n o n p (3.7) in Cartesian coordinates. Here, δ denotes the Kronecker symbol and α is a parameter which needs to be chosen either sufficiently large or suitably small, see Kabel et al. [60]. Keeping the former in mind, we turn our attention to the offline training. We sample N s quadruples of input stiffnesses and fiber orientations (C s 1 , C s 2 , λ s 1 , λ s 2 ), generate the respective microstructures and compute the effective stiffnesses C s . We denote the generated training data by as sequence of quintuples where s enumerates the sample index. The actual sampling process will be discussed in Section 4.4. For the moment, we assume the training data to be given and fixed. We express the offline training as an optimization problem The quadratic penalty term encodes the mixing constraint We implemented the offline training in PyTorch [61], see Gajek et al. [47], making use of the framework's automatic differentiation capabilities to solve the regression problem (3.8) by means of accelerated stochastic gradient descent methods using mini batches of size N b . An epoch j consists of evaluating 3.4 for all samples of the respective mini batch, evaluating the loss function (3.9), computing the gradients ∂J/∂ p ( p j , q j , v j ), ∂J/∂ q ( p j , q j , v j ) and ∂J/∂ v ( p j , q j , v j ) by means of automatic differentiation and updating the fitting parameters During offline training, it may happen that a portion of weights w i K+1 becomes equal to zero, and remains zero due to the vanishing gradient. Liu-Wu [46] removed such sub-trees from the binary tree by deleting nodes and merging the respective subtrees. In this work, we follow Liu-Wu and compress the binary tree to speed-up the training, and eventually, the online evaluation. Fig. 5 shows a schematic of how to remove laminates from the binary tree. The former illustrates a weighted tree with edge weights corresponding to the propagated weights w i K+1 . Remember that the volume fractions of the individual laminates are computed from these weights by normalization.
(a) Perfect binary tree Binary tree with removed laminates

Online evaluation
Eventually, we seek to employ DMNs to speed up a two-scale simulation, i.e., for every Newton iteration and at every Gauss point of a finite element model, a deep material network needs to be integrated implicitly. At the same time, we want to account for arbitrary eigenvalues λ 1 and λ 2 of the fiber orientation tensor at every Gauss point of the macro simulation. In this work, we extend the fast and flexible solution technique, introduced in Gajek et al. [47], to compute the effective stress of a two-phase DMN. We restrict to the two-potential framework of small-strain isothermal generalized standard materials (GSM) [52]. In d spatial dimensions, a GSM is a quadruple (Z, ψ, φ, z 0 ) ∈ GSM comprising a (Banach) vector space Z of internal variables, a free energy density ψ : Sym(d) × Z → R, a dissipation potential φ : Z → R ∪ {+∞} and z 0 ∈ Z serves as the initial condition. We assume that the dissipation potential φ is proper, convex, lower semi-continuous and satisfies φ(0) = 0 as well as 0 ∈ ∂φ(0), where ∂φ denotes the sub-differential of the convex function φ. For the complete documentation of the GSM model structure in a continuous setting, see, e.g., Gajek et al. [47]. Suppose that the two GSMs G 1 = (Z 1 , ψ 1 , φ 1 , z 0,1 ) and G 2 = (Z 2 , ψ 2 , φ 2 , z 0,2 ) are given. A time discretization of both phases i ∈ {1, 2} by the implicit Euler method gives rise to the formulae for discretized stress and Biot's equation, respectively, Here, t = t n+1 −t n denotes the time increment and the superscript n refers to the n-th time step at time t n . Due to the time discretization and freezing of the internal variables z n i at time t n , each GSM reduces to a nonlinear elastic material, see Lahellec-Suquet [62]. The condensed free energy is solely dependent on the input strain ε n+1 i and the internal variables z n i of the last (converged) time step. Then, the stress response reads For the sake of readability, we omit explicit reference to time step n+1. First, let us collect the lamination directions of all laminates in a single vector n ∈ (R d ) 2 K −1 with the same ordering that we used for the vectors p and q, i.e., We consider the displacement jump vector a ∈ (R d ) 2 K −1 , which inherits its ordering from n and the vector of strains ε = [ε 1 , ε 2 , . . . , ε 2 K ] ∈ (Sym(d)) 2 K . By introducing the gradient operator A λ1λ2 : w.r.t. the macro strain E and the unknown displacement jumps a. Here, the shorthand notation E = [E, E, . . . , E] ∈ (Sym(d)) 2 K is used. Indeed, for the work at hand, the gradient operator, which encodes the DMN's topology and lamination directions, depends on the fiber orientation parameters λ 1 and λ 2 .
To illustrate this concept, consider the following example. For a two-phase DMN of depth three, A λ1λ2 takes the following form with the symmetrization operators w.r.t. the lamination direction n i k,λ1λ2 = n i k (λ 1 , λ 2 ) as building blocks. Since we defined n i k (λ 1 λ 2 ) to depend on the parameters λ 1 and λ 2 explicitly, see Section 2.3, the gradient operator depends on the fiber orientation, as well. We account for this situation in our notation, i.e., we write A λ1λ2 and N i k,λ1λ2 . For the application at hand, i.e., integrating a DMN at every Gauss point during a two-scale simulation, the fiber orientation parameters λ 1 and λ 2 are held fixed. Indeed, we assume that the microstructure does not evolve under the applied load. Let us define the vector of internal variables of the last converged time step z n = z n 1 , z n 2 , z n 3 , . . . , z n 2 K ∈ Z := Z 1 ⊕ Z 2 ⊕ · · · ⊕ Z 1 ⊕ Z 2 and let Ψ : (Sym(d)) 2 K × Z → R be the averaged condensed free energy of the flattened laminate (3.19) alternating between the two given condensed free energies Ψ 1 and Ψ 2 . Then, we wish to solve the Euler-Lagrange equation of the DMN for the unknown displacement jumps a, where is the vector of phase stresses. The strain-wise "mass" matrix W : associates the weight w i K+1 , i = 1 . . . 2 K , to the corresponding phase strain ε i . We solve the Euler-Lagrange equation (3.20) by Newton's method. For an initial guess a 0 ∈ (R d ) N −1 , the displacement jump vector a is iteratively updated a j+1 = a j + s j a j , where the increment a j ∈ R d 2 K −1 solves the linear system

23)
A step size s j ∈ (0, 1] strictly less than unity may arise by backtracking. The Jacobian ∂ σ/∂ ε( E + A λ1λ2 a j , z n ) is a block-diagonal matrix containing the algorithmic tangents of the DMN's input materials, i.e., Upon convergence, the phase strains ε = E + A λ1λ2 a and, subsequently, the effective stress are computed by averaging. To determine the algorithmic tangent of the deep material network, for a start, the linear system is solved for ∂ a/∂E. Then, the algorithmic tangent may be represented in the form where ∂ σ/∂E denotes the vector of algorithmic tangents and [I s , I s , . . . , , is a vector of the identity operators on Sym(d).
By comparing equation (3.26) to equation (3.23), we observe that both problems share the same linear operator, but with different right hand sides. When using a direct solver, e.g., a Cholesky decomposition, it is recommended to reuse the matrix decomposition for reasons of efficiency.
To reduce the number of degrees of freedom and, thus, to speed up the solution process, we exploit that some weights become zero during training as explained in the previous section. We learned that in the offline training, we can dynamically build a binary tree with simplified topology but identical effective behavior. This is also true for the offline evaluation of the DMN. The DMN's topology is encoded by the gradient operator A λ1λ2 . Deleting laminate blocks from the binary tree is equivalent to deleting the associated rows and columns of A λ1λ2 . For the example shown in Fig. 5, we obtain a (reduced) gradient operator of the form with a reduced size, where we dropped the subscript (·) λ1λ2 of N i k,λ1λ2 for readability.

Identifying a DMN surrogate model
This section is dedicated to the identification of the DMN surrogate model. We discuss the pre-processing steps, i.e., finding the necessary resolution and the appropriate size of the volume elements, and investigate the discretization of the fiber orientation triangle. Subsequently, we explain the sampling of the training data, the offline training and the validation of the DMN surrogate model on the fiber orientation triangle. All computations were performed on a workstation equipped with two AMD EPYC 7642 with 48 physical cores each, enabled SMT and 1024 GB of DRAM.

Short glass fiber reinforced polyamide
For the work at hand, we focus on a short glass fiber reinforced polyamide. We consider E-glass fibers with a length of L f = 200 µm and a diameter of D f = 10 µm. The glass fibers are assumed to be isotropic, linear elastic. The fiber volume fraction is set to c f = 16% corresponding to a fiber mass fraction of approx. 30%. The matrix is assumed to be governed by J 2 -elastoplasticity, see Chapter 3 in Simo-Hughes [63] with an exponential-linear hardening The mechanical properties used in the simulation, taken from Doghri et al. [64], are summarized in Tab. 2. We rely upon the Sequential Addition and Migration (SAM) method [51] for generating periodic volume elements with prescribed volume fraction and second order fiber orientation tensor. The fiber length L f ,  Table 2: Material parameters for the short glass fiber reinforced polyamide [64] fiber diameter D f , fiber volume fraction c f and the axis aligned fiber orientation tensor A 2 , i.e., λ 1 and λ 2 , serve as input parameters for the SAM method. Please note that we only consider fiber orientation states with λ 3 ≥ 0.01, as purely planar fiber orientation states cannot be generated at high filler content essentially for geometric reasons, see Schneider [51] for a discussion.

On the necessary resolution and the size of the RVE
For a start, we study the resolution necessary to obtain accurate effective properties in the purely elastic case. For this purpose, we consider cubic microstructures with an edge length of L = 384 µm, i.e., roughly twice the fiber length of L f = 200 µm. We compute the effective stiffness with the help of an FFT-based computational micromechanics code [7,8] as described in Schneider [65], using the staggered grid discretization [66,67] and the conjugate gradient solver [68,69]. We consider the extreme orientations individually, i.e., unidirectional, isotropic and planar isotropic fiber orientation as shown in Fig. 3, and vary the resolution from 1.7 to 13.3 voxels per fiber diameter in equidistant steps. This corresponds to volume element discretizations with 64 3 to 512 3 voxels. We measure the error relative to the effective stiffness C and choose a resolution of 20 voxels per fiber diameter, i.e., a discretized by 768 3 voxels, as the reference.   Fig. 6(a) shows the relative error of the effective stiffness computed by the Frobenius norm of the corresponding Voigt matrices. For the crudest resolution, i.e., five voxels per fiber diameter, the relative error is well below 2% for all three considered fiber orientations. Notice that the relative error for the volume element with planar isotropic fiber orientations is consistently larger than the error for the unidirectional and isotropic orientations. As expected, the relative error decreases with increasing resolution. At a resolution of 6.7 voxels per fiber diameter, the relative errors of the unidirectional and isotropic fiber orientation fall below 1%. For 8.3 voxels per fiber diameter, the relative error of the planar isotropic fiber orientation is below 1%, as well. For this article, we consider a resolution of 6.7 voxels per fiber diameter as sufficient, i.e., relative errors below 1% for isotropic and unidirectional fiber orientation and an error slightly above 1% for the planar isotropic case. We fix this resolution and focus on finding the size of a representative volume element. For a resolution of 6.7 voxels per fiber diameter, we investigate volume elements with edge lengths L ranging from 1.44 up to 3.84 fiber lengths corresponding to volume element discretizations with 192 3 up to 512 3 voxels. To obtain our reference, we generated a volume element with edge lengths of 7.68 fiber lengths discretized with 1024 3 voxels. As for studying the necessary resolution, we again consider the relative error in the effective stiffness as our error measure. For the volume elements with edge lengths of 1.9 fiber length and above, the relative error is well below 0.5% and does not further decrease significantly for increasing volume element size. For this reason, we consider volume elements with an edge length of L = 384 µm as sufficient. To sum up, we finally choose a resolution of 6.7 voxel per fiber length, i.e., a voxel size of 1.5 µm and a discretization with 256 3 voxels for the article at hand.

Discretization of the fiber orientation triangle
To generate the linear elastic training data, we seek to sample the space of input stiffnesses and fiber orientations uniformly. Apart from sampling the input stiffnesses, it is possible to sample λ 1 and λ 2 as well, e.g., via a low-discrepancy sequence such as the Sobol sequence [70] or via Latin hypercube sampling [71]. Due to the high dimension of the input space, we follow the former approach for generating tuples of input stiffnesses. The considered fiber orientations are parameterized by a two-dimensional space, and traditional methods are more efficient. We discretize the fiber orientation triangle by partitioning it into four self-similar triangles, which may be subsequently partitioned, as well. We select the three points on the vertices of each triangle plus the centers of the triangles as sampling points for the parameters λ 1 and λ 2 . We start with the full orientation triangle, see Fig. 7(a). The four sampling points, illustrated by four hollow circles, are the three corners and the center of the orientation triangle. After the first splitting, the orientation triangle comprises four triangles and ten points, see Fig. 7(b). After dividing the triangles one more time, we arrive at 31 sampling points. For each of these points in fiber orientation space, we generate a single volume element using the SAM method [51] with edge lengths L = 384 µm and a discretization with 256 3 voxels, see Section 4.2. In this work, we consider the discretizations shown in Fig. 7, i.e., we discretize the orientation triangle with four, ten and 31 sampling points which we call D4, D10 and D31, respectively. Choosing the sampling points in a hierarchical manner permits us to re-use already generated volume elements for the next finer discretization. For instance, going from D10 to D31 only requires generating 21 new volume elements. This restriction is not imposed by the SAM [51] algorithm since generating volume elements with the given resolution and edge length is just a matter of milliseconds to seconds. Being able to reuse results from coarser discretizations comes in handy for the inelastic validations we perform in Section 4.6. In Section 4.6, we generate 78 additional volume elements to validate the DMN outside of its training regime and to check if the DMN generalizes sufficiently on the orientation space. For this purpose, we test several different load paths for each volume element, so being able to reuse already computed full field solutions comes in handy to keep the validation effort manageable.
To sum up, we consider the following eight degrees of freedom (a, e 1 , e 2 , e 3 , β, θ, ψ, φ) (4.10) with their respective domains specified above. To sample the input space evenly, we sample the parameters (4.10) by the Sobol sequence and, subsequently, construct the stiffness tensors (C 1 , C 2 ).

Offline training
For the offline training, we generate N s pairs of stiffnesses (C s 1 , C s 2 ) by the protocol described in Section 4.4. Then, we assign each stiffness tuple to one of the previously generated volume elements in a cyclic fashion. For instance, for the orientation discretization D4, we assign C 1 1 , C 1 2 to the volume element with fiber orientation λ 1 1 , λ 1 2 , C 2 1 , C 2 2 to the volume element with λ 2 1 , λ 2 2 and C 5 1 , C 5 2 to the volume element with λ 1 1 , λ 1 2 and so forth. For every quadruple (C s 1 , C s 2 , λ s 1 , λ s 2 ), we compute the associated effective stiffness C s with the help of an FFT-based computational micromechanics code [7,8]. The serves as training data for identifying the DMN. The number of samples depends on the discretization of the orientation triangle and is summarized in Table 3. For D4, we generate 800 samples in total which corresponds to 200 samples per volume element. To keep the sampling and training effort manageable, we reduce the number of generated samples to 100 and 50 per microstructure when increasing the number of discrete orientations to ten and 31, respectively. For D4, D10 and D31, we randomly split the pre-computed samples into a training and validation set, comprising 90% and 10% of the samples. We train the deep material network on mini-batches with a batch size of N b = 32 samples. More precisely, we draw the batches randomly from the training set and drop the last batch, should the remaining batch size be smaller than 32.

D4
D10 D31 Total 800 1 000 1 550 Per microstructure 200 100 50 Training set 720 900 1 395 Validation set 80 100 155 We consider deep material networks with eight layers. In general, eight layers are necessary to achieve a sufficient approximation quality, in particular for inelastic computations [45][46][47]. We train for 3 000 epochs and rely upon the AMSGrad method [73,74] combined with the warm restart technique suggested by Loshchilov-Hutter [75]. Furthermore, we make use of a modulation of the learning rate between a minimum learning rate α min and a maximum learning rate α max , i.e., α : N → R, m → γ m α min + 1 2 (α max − α min ) 1 + cos π m M , (4.11) where 2M corresponds to the period and M = 50 is chosen. Additionally, we decay the learning rate at a geometric rate with γ = 0.999. Since gradient descent is sensitive w.r.t. the proper choice of the step size, we determine the learning rates α p , α q and α v by a learning rate sweep as introduced by Smith-Topin [76]. The resulting learning rates are almost identical for all three parameter groups, and we set α max = 1.5 · 10 −2 . The minimum learning rate is chosen to be an order of magnitude smaller than the maximum learning rate, i.e., α min = 1.5·10 −3 .
We sample the initial weights v from a uniform distribution on [0, 1] and rescale the weights to sum to unity. The entries of the parameter vectors p and q are sampled from a uniform distribution on [0, 2π].
The penalty parameter of the objective function (3.9) is set to λ = 10 3 . Additionally, we set the exponents to p = 1 and q = 10, see Gajek et al. [47], i.e., we enforce the maximum of the component-wise mean error to be minimized. To assess the accuracy of the fit, we define the sample-wise error where · 1 refers to the Frobenius-1 norm defined by the 1 -norm of the stiffness components in Voigt notation. Additionally, we define the maximum and mean errors of all samples where N s denotes the number of elements in the training or validation set, depending on the considered scenario.
In Fig. 8, the training progress for the D31 orientation discretization and the investigated linear, tri-linear and quadratic orientation interpolations is shown. In the first 500 epochs, the effect of the learning rate modulation becomes apparent. The loss function and mean error fluctuate noticeably. In the last 500 epochs, the decay of the learning rate ensures convergence of the trained parameters. Conforming to intuition, increasing the degrees of freedom, i.e., choosing a tri-linear or quadratic orientation interpolation over a linear orientation interpolation, decreases the loss function at convergence. This trend carries over to the mean training errors, as well. For the mean validation error, however, the linear orientation interpolation provides the best mean validation error. Such overfitting phenomena are not uncommon for training deep neural networks, where increasing the degrees of freedom not necessarily yields better generalization and validation results. As training progresses, no increasing validation errors can be observed for linear, tri-linear and quadratic orientation interpolation. Thus, no significant model over-fitting is observed during training. Tri-linear 1.175 · 10 −2 0.700% 1.928% 0.879% 6.018% 65% Quadratic 9.757 · 10 −3 0.620% 1.361% 0.841% 7.051% 68% Table 4: Training results of the short glass fiber reinforced polyamide  In Table 4, we summarized the training results for the investigated orientation discretization and interpolation. Additionally, we listed the number of non-zero weights w i K+1 at the end of training. In general, we observe the following trends: Mean and maximum training errors will decrease if a tri-linear or quadratic orientation interpolation is chosen instead of a linear interpolation. The same observation holds for the mean validation error. The maximum validation error, on the other hand, does not necessarily decrease by introducing additional fitting parameters. Furthermore, the maximum validation error shows significant fluctuations. Indeed, for D4 and the linear interpolation, it reaches almost 12%. Since the maximum validation error is dominated by a single sample, see Liu et al. [45,46] or Gajek et al. [47], we consider the mean validation error to be a more appropriate indicator of the quality of the training results. If we go from D4 to D10 and D31, the loss function as well as the training and validation errors will increase, in general. Indeed, the DMN has to predict the effective behavior of significantly more volume elements with different fiber orientations and this result does not come unexpected.

Online evaluation
We implemented Newton's method, as described in Section 3.2, as a user-material subroutine in ABAQUS. The pseudo-code for the implementation can be found in Gajek et al. [47]. In terms of implementation, the major difference compared to Gajek et al. [47] is that, for the work at hand, the gradient operator depends on the fiber orientation parameters λ 1 and λ 2 , see Section 3.2. This does not infer any additional challenges, since the microstructure characteristics do not change during computation. Both parameters λ 1 and λ 2 are fixed during the online evaluation, and after assembling A λ1λ2 , the effective stress and algorithmic tangent are computed by Newton's method in the proposed manner. Our goal in Section 5 is to employ deep material networks for a two-scale simulation using ABAQUS. For this purpose, we seek to speed up the inelastic computations as much as possible. As a first step, the elimination procedure described in Section 3 proves effective. The deleting is performed in an upstream pre-processing step after the offline training to avoid unnecessary computational overhead. Secondly, we exploit the sparsity pattern of all involved linear operators, both, the sparsity pattern of the gradient operator and the Jacobian ∂ σ/∂ ε, containing the algorithmic tangents of the phases. To this end, we rely upon the library Eigen3 [77] for all linear algebra operations and use sparse matrices whenever possible. For Newton's method, we use the following convergence criterion where we set the tolerance tol to 10 −12 and · F refers to the Frobenius norm defined by the 2 -norm of the involved matrices in Voigt notation. We solve the linear system by means of a sparse Cholesky decomposition. Before employing the DMN in a two-scale simulation, we turn our attention to validating the predicted effective stress of the deep material network. To evaluate the approximation error in a quantitative way, we introduce the following error measures. For fixed orientation parameters λ 1 and λ 2 , we define the relative error in the stress component (i, j) as where T = [0, T ] denotes the considered time interval. Furthermore, the mean and the maximum error are defined by and respectively. To investigate whether deep material networks are capable of accurately interpolating the effective stress outside of the training regime, additional volume elements were generated. We use D4 as our point of departure and subdivide the orientation triangle three more times yielding 105 additional sampling points on the orientation triangle. For D10 and D31, we subdivide the orientation triangle two and one more times, giving rise to 99 and 78 additional sampling points, respectively. In Fig. 9, for all three discretizations, the sampling points used for obtaining the training data (hollow circles) and the sampling points exclusively used for the inelastic validations (black filled circles) are shown. For every additional sampling point, we generated a volume element using SAM [51]. We keep using the volume elements already generated for offline training such that we have 109 generated volume elements in total for every D4, D10 and D31.  (c) Quadratic interpolation Figure 12: Minimum, mean and maximum errors w.r.t. mean η mean λ1λ2 and maximum η max λ1λ2 error. Fig. 12, which compares 5 886 computed DMN load paths with 654 full field simulations, permits us to draw the following conclusions. For fixed type of orientation interpolation, a finer orientation discretization reduces the mean and maximum errors. The dependence of the mean and maximum error on the polynomial degree is more complicated. Compared to the linear orientation interpolation, both, the mean and the maximum error increase significantly for tri-linear and quadratic interpolation. The former holds true for D4, D10 and D31. Similar to the offline training, the loss and training errors decrease for higher order interpolations. Only the validation errors shows a slight increase. This suggests an overfitting of the online phase, so that the linear approaches to the orientation interpolation is recommended, in the end.

A computational example
After validating that the DMN is able to to provide sufficiently accurate results, we turn our attention to a component of industrial complexity. We choose a quadcopter frame, see Fig. 1(b), whose CAD geometry is publicly available [1]. We assume that the arms of the quadcopter are manufactured by injection molding. To obtain realistic fiber orientation data, we conducted a mold filling simulation for a single quadcopter arm (and utilize these data for all four identical quadcopter arms). We use the publicly available software InjectionMoldingFoam [80], which is based upon the two-phase, incompressible flow solver of OpenFOAM [81]. We choose identical settings and material parameters as Köbler et al. [50], i.e., we assume a homogeneous fiber volume fraction and select the following Carreau-WLF equation [54] η(θ,γ) = η 0 e −A2(θ−θ ref ) Here, θ denotes the absolute temperature andγ refers to the norm of the strain-rate tensor. The parameters for the injection molding simulation are summarized in Table 5 and originally stem from Bhat et al. [82]. The results of the injection molding simulation are shown in Fig. 13, both for a top and a Density: 1 410 kg/m 3 Folger-Tucker coefficient:  Table 5: Parameters used in the injection molding simulation [82] side view, and at three distinct instances of time, corresponding to a volume coverage of 50%, 75% and 100%, respectively. The computed fiber orientations are represented by the color scale shown in Fig. 3.
We observe that, at the flow fronts, planar and isotropic fiber orientations dominate. For the rest of the drone arm, fiber orientations close to the unidirectional case are prevalent. Fig. 13(d) illustrates the principal fiber orientations, i.e., the eigenvector corresponding to the largest eigenvalue of A 2 , after the mold is filled.  First, we perform a structural simulation on a single drone arm using the FE software ABAQUS. We mesh the drone arm by quadratic tetrahedron elements and investigate five different mesh densities ranging from 63 580 up to 1 005 862 elements in order to analyze convergence, wall-clock times and memory requirements, see Table 9. The computed fiber orientation tensors serve as the input for the simulation, i.e., the eigenvectors of A 2 are mapped onto the ABAQUS mesh and determine the material orientation. The eigenvalues λ 1 and λ 2 are provided to the DMN subroutine via pre-defined fields. We apply a loading of F = 80 N on the motor mount via a surface force and fix the left side of the drone arm, see Fig. 14.  The loading is applied in ten equidistantly spaced time increments. Fig. 14 shows the results for the finest discretization of about one million tetrahedron elements. The former took about 165 min to complete on 96 threads and required 133 GB of DRAM. On the left hand side, the von Mises stress distribution for the last time increment and for an assumed homogeneous and isotropic fiber orientation is shown. The right hand side of Fig. 14 shows the computed stress for the mapped anisotropic, inhomogenous fiber orientation. For the mapped anisotropic fiber orientation, stress fluctuations, especially in the vicinity of weld lines are clearly visible. In contrast, stress and strain concentrations at weld lines cannot be predicted for a homogeneous fiber orientation. Accounting for the entire process chain appears imperative in order to exploit the full lightweight potential of injection-molded fiber-reinforced components, as becomes evident when comparing the predicted total deflections. Indeed, for the assumed isotropic fiber orientation, the macro simulation predicts a deflection of 6.13 mm. A deflection of 3.99 mm is predicted for the anisotropic fiber distribution. Thus, the isotropic variant underestimates the actual stiffness of the component by a factor of two.
To demonstrate the capabilities of the introduced multiscale method, we investigate the entire drone frame in a mechanical simulation, see Fig. 15. The four drone arms are manufactured from injection molded, short fiber reinforced polyamide with mapped anistotropic fiber orientation. A deep material network is integrated at every Gauss point. Both, the upper and lower plates, which the drone arms are attached to, are made of aluminum. For this material, we use a J 2 -elastoplasticity model with power law hardening The material parameters are taken from Segurado et al. [83] and summarized in Table 6. We assume the Aluminium E = 75 GPa ν = 0.3 σ Y = 75 MPa k = 416 MPa m = 0.3895 Table 6: Material parameters of the aluminum plate [83] drone legs to be made of pure polyamide. For the latter, we assign a linear elastic material behavior, see Table 2. The simulation model consists of about two million elements with almost ten million degrees of freedom, see Table 10. The drone arms are loaded as shown in Fig. 15 Figure 15: Simulated quadcopter drone frame

Computational cost
Last but not least, we discuss the computational cost of deep material networks accounting for the offline training and online evaluation separately. The material sampling is performed in parallel, i.e., we compute six load steps in parallel using 16 threads for each individual simulation. All deep material networks are trained in parallel on 4 threads each. The wall-clock times of the material sampling and the offline training are summarized in Table 7. Apparently, the sampling and offline training effort increases linearly with the number of samples. Incidentally, the number of fitting parameters, which varies depending on the type of the orientation interpolation, has no significant influence on the runtime of the offline training.  Turning our attention to the online evaluation, we focus on the computational costs of the DMN evaluated at a single Gauss point. Integrating a deep material network at a single Gauss point for a prescribed macro strain increment takes about 2 ms on a single thread. This is about 120 000 times faster than a full-field simulation of the micro problem using an FFT-based computational mechanics solver on a microstructure discretized by 256 3 voxels, also running on a single thread. For the work at hand, we exclusively consider DMNs with eight layers. For applications which permit using DMNs with a smaller number of layers, even higher speed-up factors can be reached.
FFT (1 thread Table 8: Wall-clock times and speed-up (compared to an FFT-base computational micromechanics solver) for a single time step of the inelastic micro simulation Next, we focus on the component scale simulation of the quadcopter arm. The wall-clock times of all five examined discretizations ranging from 63 580 up to 1 005 862 quadratic tetrahedron elements and, computed on 96 threads, are summarized in Table 9. We observe that the DRAM footprint is roughly proportional to the number of elements. The wall-clock times, however, increase super-linearly. We attribute this effect to the complexity of the direct solver used by ABAQUS. Apparently, the applicability of the method is more restricted by the memory requirements, and the computational effort plays a minor role.  Table 9: Wall-clock time and memory consumption of the single drone arm for different mesh sizes jumps of the last converged time step are stored as well, i.e., 768 + 255 × 3 = 1 533 scalars need to be kept in memory. Since we rely upon the thinned binary tree as introduced in Section 3.2, 1 533 serves as an upper bound. For the application at hand, the actual number of internal variables of the DMN surrogate model is 1 297.
To analyze the drone arm, using less than one million elements for the discretization would be sufficient. Rather, by choosing such a fine discretization, we demonstrate that deep material networks easily scale to component scale simulations with higher complexity. The hardware requirements implicated by Table 9 can be provided by any state-of-the-art workstation.
Computing all ten load steps on 96 threads for the entire quadcopter took 267 min and required 252 GB or DRAM, see Table 10. This corresponds to over two million elements and about ten million degrees of freedom. In our opinion, this clearly shows that DMNs are a promising technique to enable two-scale simulations of industrial complexity with manageable resource expenditure.

Conclusion
In this work, we investigated the capabilities of deep material networks to provide a digital twin for short fiber reinforced plastic microstructures, which can be used in concurrent multiscale simulations. To realize the full lightweight potential of short fiber reinforced components, it is imperative to account for the locally varying fiber orientation in mechanical simulations on component scale. Building upon the work of Köbler et al. [50], we proposed a robust and computationally efficient approach to utilize direct deep material networks for variable fiber orientations. Instead of identifying multiple deep material networks and interpolating the effective stress, we interpolated the DMN's microstructure characteristics on the fiber orientation triangle. Assuming that the local fiber volume fractions of the individual laminates in the hierarchy are independent of the local fiber orientation, it suffices to fix the fiber volume fraction and to interpolate the lamination directions only. This procedure gives rise to a single DMN surrogate model covering all fiber orientations. Presumably, the scheme easily extends to incorporating local variations in the fiber volume fraction by interpolating the DMN's volume fractions as well. By sampling the training data from up to 31 microstructure realizations with different fiber orientation, we fitted the DMN to multiple fiber orientations simultaneously. Subsequently, we showed that the DMN generalizes to the entire fiber orientation triangle with small error, also for the inelastic regime.
To evaluate the ensuing performance of our approach, we simulated the entire process chain of a quadcopter frame starting from an injection molding simulation. We mapped the computed fiber orientations upon a finite element mesh of the complete quadcopter frame and conducted a two-scale simulation of the full component. Our results indicate that deep material networks enable two-scale simulations of structures with industrial complexity with moderate hardware requirements. In this way, the FE-DMN method finally realizes the promise of concurrent multiscale simulations and constitutes a powerful piece of technology which promises to become a standard tool for industrial applications. It should be interesting to extend the FE-DMN method to problems involving damage or fracture [84][85][86], finite strains [87,88] or thermomechanical coupling [89,90] and to extend the range of applicability, for instance in the context of polycrystalline materials [91,92] or sheet molding compound composites [93,94]. Despite the apparent success in practice, there is still a need for theoretical results which shed light on the approximation capabilities (and the limitations) of DMNs. Indeed, whether every fixed two-phase microstructure has a microstructure twin of hierarchical laminate with identical effective properties, appears to be unresolved, see Problem 4 in Milton [95]. Interestingly, there exist counter examples for five-phase composites where the former is false, see Milton [56].