Assumed strain methods in micromechanics, laminate composite voxels and level sets

This work deals with the composite voxel method, which—in its original form—furnishes voxels containing more than one material with a surrogate material law accounting for the heterogeneity in the voxel. We show that the laminate composite voxel technique naturally arises as an assumed strain method, that is, the general framework introduced by Simo‐Rifai, for a specific choice of enhanced strain field. As a consequence, laminate composite voxels may be regarded as a kinematic assumption within a discretization scheme rather than a part of material modeling, as suggested originally. We discuss how to seamlessly integrate composite voxels into the framework of a level‐set description of the microstructure, in particular the accurate and efficient computation of normals and cut‐volume fractions. In contrast to more traditional strategies based on subvoxelizations, the introduced method avoids systematic errors when computing composite voxel properties. We demonstrate the applicability of the developed technology for a number of relevant computational examples.

The advent of high-fidelity digital images of microstructures, for example, serial sectioning, 6,7 optical microscopy, 8,9 FIB-SEM, [10][11][12] electron back-scatter diffraction, [13][14][15] x-ray diffraction microscopy, 16,17 and micro-computed tomography, 18,19 has changed the field completely, offering new possibilities and posing corresponding challenges.For instance, digital images are given on voxels (volume pixels), typically at a high resolution.In particular, a high number of such image points need to be processed.Moreover, the rectangular nature of voxels does not permit to determine the interfaces between the materials in a natural way.Rather, the classical segmentation process leads to a jagged interface.This circumstance is made worse by the inherent complexity of industrial-scale microstructured materials.In particular, when it comes to computationally resolving such microstructure models, it turns out to be rather challenging to create interface-conforming meshes with high quality.4][25] An efficient alternative operates on the regular grid provided by the voxel mesh directly. 26,279][30] Particularly well-performing are approaches based on the fast Fourier transform (FFT), which were introduced by Moulinec-Suquet 31,32 in the 1990s.Originally, their approach was motivated by the Lippmann-Schwinger integral equation for the strain field [33][34][35] and exploited a truncated set of Fourier frequencies to represent the sought periodic strain field.Designed to account for inelastic constitutive behavior from the beginning, 31,32 the available efficient implementations of the FFT 36 led to an ever-growing community of adherents exploiting the computational advantages of the approach, see the recent review articles [37][38][39] for an overview.
In its original form, Moulinec and Suquet presented their scheme in an integrated fashion, not distinguishing discretization scheme and solution method.In the succeeding years, the original solution method, the basic scheme, 31,32 was replaced by other solution methods with a superior convergence behavior.These improvements come in two flavors.6][47] For the second category, the kinematic compatibility is only satisfied upon convergence.4][55] Of course, it is also possible to swap kinematic compatibility with equilibrium and work with dual schemes. 56,57t the beginning of the 21st century, there was a debate within the community whether FFT-based methods converge for problems with infinite contrast, most prominently porous materials.Due to the regular grid structure and the use of FFT, also these void voxels needed to be included into the evaluation of the constitutive law, that is, the classical strategy of simply eliminating such regions from the grid is not feasible.Originally, this problem was thought of as being the concern of the solution scheme and appropriate improvements were sought.However, depending on the microstructure, constitutive law and convergence criterion, contradictory results were reported.Subsequently, it turned out that the root of the problem lies in the discretization scheme based on trigonometric polynomials whose global support may lead to an ill-conditioned micromechanical problem in the presence of pores.This defect was demonstrated by a computational experiment of a high-porosity foam 58( §4.2) where the predicted effective Young's modulus converges to zero eventually with increasing iteration count, contradicting physical intuition, that is, that the foam structure has a finite stiffness.As a remedy, different discretization schemes were exploited, building upon earlier work of Willot et al. 59,60 Indeed, it was recognized that on a regular grid, Lippmann-Schwinger type equations and solvers could be designed and used, as long as the discretization scheme respects the regular grid structure and the periodic boundary conditions.Such alternative discretization schemes include Fourier-Galerkin methods, [61][62][63] piece-wise constant strains, 64,65 finite differences, 58,60,[66][67][68][69] finite elements, [70][71][72] finite volumes, 73,74 splines, 75 and wavelets. 76espite the computational prowess of FFT-based methods, performance improvements are still valuable, in particular in case of expensive constitutive laws [77][78][79] or whenever a multitude of microstructures [80][81][82] or load scenarios [83][84][85] needs to be considered.A key technology for accelerating FFT-based methods is the composite voxel technique, originally introduced by two groups 86,87 independently at the same time.Interestingly, similar ideas were introduced for phase-field problems. 88,89The key insight is that in case the microstructure is given either analytically or in terms of a high-fidelity subvoxel image, there are voxels which comprise more than a single phase.To account for this heterogeneity within the voxel, it appears natural to furnish the voxel with an appropriately chosen composite constitutive law.Different possibilities were tested and for finitely contrasted materials it turned out to be favorable to equip the voxel with the constitutive law of a laminate 90( §9) with the proper volume fractions and a lamination direction which is close to the normal between the phases within the voxel.Interestingly, such a strategy does also work at finite strains 91 and for inelastic constitutive laws with internal variables. 92Subsequent works also exploited thin interphases 93,94 and composite boxels 95 on non-equi-axed elements.

Contributions
This work revisits the composite voxel method, mainly for the following two reasons.For a start, there appears to be no thorough explanation for the unreasonable effectiveness of laminate composite voxels.Indeed, in its original form, 87 the composite voxel methodology operates on the level of the constitutive law, touching neither kinematic compatibility nor the balance equation.In particular, pretty much any other "mixing rule" would be admissible.Originally, different such rules were tested, but not deemed performing.Subsequent works 91,92,94 highlighted the advantages of the laminate strategies also for general inelastic and possibly finite strain mechanics, leading to the impression that there may be a deeper reason for the prowess of the composite voxel technology.
The second motivation was a more practical one.Originally, composite voxels were coupled with a subvoxel description and an analytic description of the interfaces given via computer-aided design (CAD) data.For both strategies, accurate interface data is either not available at all or requires serious processing.An alternative which is naturally integrated with voxel data is the level-set method, [96][97][98] and we wanted to study how to design efficient composite voxel methods in this setting.
More to the point, it is necessary to obtain both the volume fractions within a composite voxel and an accurate interface normal.When testing our proposed techniques, see Section 3 for details, we uncovered that the original subvoxel-based normal-estimation technique introduced by Merkert et al. 99 is not multi-grid convergent.As this fact is of independent interest, we will discuss it in this central and prominent part of the manuscript.The strategy introduced by Merkert et al. 99 proceeds as follows.Suppose there is a (possibly non-normalized) level-set function L e ∶ [0, 1] 3 → R on the voxel in question which is sampled on a subvoxel grid with m 3 voxels, that is, the signs of the values are essential.Then, based on the centroids of the discretized phases and the centroid x c = (1∕2, 1∕2, 1∕2) of the voxel, the estimated normal is computed in either of the equivalent forms 95 Restricting to a two-dimensional setting, an example with a linear interface is shown in Figure 1.The strategy (4) provides the correct normal n e ≡ (1, 2)∕ √ 5 for an interface which cuts the left upper corner in Figure 1A.However, when shifting the interface to the right, the center of mass of the blue sublevel set shifts, as well.At the same time, the interface normal remains the same.In this way, a difference between the correct and the estimated normal is introduced, see Figure 1B.For the shown scenario with n ∞ e ≡ (11, 4)∕ √ 137, the difference between the two normals is approximately 6.58 • .Of course, using a finite number m 3 of integration points does not help, either, and the discussed scenario prevails as m → ∞.
Returning to the principal scientific contributions of the article, we provide a derivation of laminate composite voxels from first principles.More precisely, we derive the latter within the framework of assumed strain methods, see Section 2. To provide some background, consider, as an example, a finite element model of a beam-like component under bending.Then, linear elements have serious trouble representing the linearly varying strain field as the chosen shape functions lead to piece-wise constant strain fields only.Clearly, utilizing higher-order ansatz functions would mitigate the problem.Yet, such higher-order shape functions are characterized by a higher implementation and computational effort.Therefore, different strategies were exploited [100][101][102][103][104][105][106][107] which enrich the first-order shape functions by appropriate modes that permit to represent a linearly varying strain.These additional degrees of freedom were typically chosen in such a way that they are local to each element, that is, independent of the degrees of freedom associated to nearby elements.In particular, it was possible to eliminate these additional degrees of freedom in preprocessing, leading to a modified strain-displacement matrix per element.The resulting approaches share the same number of degrees of freedom as conforming first-order methods but are able to approximate bending-dominated loadings much more accurately.These methods go under different names, for example, the method of incompatible modes, [100][101][102] as the element-local additional shape functions lead to a violation of displacement continuity if interpreted as a displacement enrichment.Alternative denominations include B-bar methods 103,104 due to the altered strain-displacement matrix which arises from the static condensation or enhanced strain methods [105][106][107] which consider the additional degrees of freedom as augmenting the strain approximation.Another name is assumed strain methods, where the focus is laid on the additive splitting of the strain into a kinematically compatible and an enhanced strain part.Irrespective of the denomination, Simo-Rifai 108 provided a unified description and analysis of such methods in a conforming three-field variational framework.They provided three simple conditions which ensure well-posedness of the finite-element problem and convergence under grid refinement.A particular consequence is that the enhanced part of the strain must have a vanishing mean on each element.Moreover, assumed strain methods typically do not improve the convergence rate of the finite element discretization under grid refinement.However, the improvement of the constant in front of the convergence rate is significant and-in combination with the comparatively low inherent computational effort-responsible for the success of such methods.
In the manuscript at hand, we interpret the laminate composite voxels as such an assumed strain method.Instead of bending of homogeneous components, which the assumed strain methods were originally designed for, we are interested in homogenization or localization problems on heterogeneous microstructures.In case the regular mesh does not conform to the interfaces between the constituents, elements cut by the interfaces are unable to provide the appropriate mechanical response at the interface.Using appropriate enhanced strains in these voxels appears to be the natural way to proceed.More precisely, within the composite voxels, the difference in normal strains across the interface cannot be accounted for by first-order elements that do not resolve the interface.Thus, it is natural to provide such additional degrees of freedom for the composite voxel.Accounting for the vanishing mean constraint required by Simo-Rai's analysis 108 directly leads to laminate composite voxels, see Section 2.3.In particular, we interpret the composite voxel method as an augmentation of kinematic compatibility, more precisely the approximation thereof in a finite element context, independent of the constitutive modeling, which appears to be much more natural than the previous interpretation.
We provide a self-contained derivation and analysis of assumed strain methods in Sections 2.1 and 2.2, culminating in sharp convergence estimates for the method.
Our second contribution concerns the seamless integration of level-set methods and the composite voxel technology.We propose a number of strategies to estimate the normal and volume fractions in composite voxels, see Section 3. Here, the problem is the ambiguity of the normal in case the interface is not flat.We introduce ideas based on averaging, regression planes for the roots of the level-set function and a combinatorial strategy that was reported in literature. 109After finding a suitable linear interface intersecting the composite voxel, the corresponding volume fractions may be computed by calculating the volume of the emerging polyhedra.We discuss the efficient and numerically robust algorithmic resolution proposed by Mirtich 110 which is based on a repeated application of integral theorems to reduce the dimension of the domain to be integrated to the edges of the polyhedron where Gaussian quadrature is known to be optimal, see Hughes 111 and the references therein.
The article at hand concludes with a number of computational experiments discussed in Section 4, where the performance of the introduced strategies is assessed based on a number of microstructures with analytically available reference values or based on a high-fidelity computed reference.

The homogenization problem
We are concerned with a rectangular cell in d spatial dimensions (d = 2, 3) and suppose a heterogeneous strain-energy function w ∶ Y × Sym(d) → R, (x, )  → w(x, ), (6)   to be given, where Sym(d) refers to the space of symmetric second-order tensors in d dimensions.We assume that the function w is differentiable in the second argument and denote the associated stress mapping as follows: A convenient framework for nonlinear solid mechanics at small strains without softening involves strongly monotone and Lipschitz-continuous stress operators, 39 that is, there are positive constants  ± , s.t. the strong monotonicity and the Lipschitz continuity condition are satisfied.Here, the norm of strains and stresses is measured in the Frobenius norm, that is, Moreover, it is assumed that the unstrained stress state has finite elastic energy, that is, the condition holds.The symbol ⨍ refers to the mean integral, that is, the integral weighted in such a way that integrating the constant function 1 yields unity.Under the conditions (8), (9), and (11), it can be shown that the stress operator  gives rise to a well-defined (nonlinear) operator on the space of square-integrable Sym(d)-valued fields which is furthermore  − -strongly monotone and  + -Lipschitz continuous, see Schneider 112( §2.1) for details.For a given prescribed macroscopic strain  ∈ Sym(d), we seek an associated periodic displacement fluctuation field u ∈ H 1 # (Y ; R d ) with vanishing mean, s.t. the quasistatic balance of linear momentum div (⋅,  + ∇ s u) = 0 (12)   is satisfied.Under the strong monotonicity and Lipschitz condition on the stress operator, the cell problem ( 12) is well-posed, and the corresponding effective stress arises via volume averaging Traditionally, the problem ( 12) is interpreted as the Euler-Lagrange equation of the variational problem This variational problem is the typical starting point for numerical discretizations in micromechanics. 58,60,70,71,113,114imo-Rifai 108 discuss an alternative mixed variational principle which we adapt to the micromechanics framework.Assumed strain methods consider a decomposition of the local strain field into a kinematically compatible part  + ∇ s u and an additional, enhanced strain ε.At the continuum scale, kinematic compatibility forces this latter part to vanish-yet it may be non-trivial and useful after discretization.Retaining the continuous setting, the variational problem with the spaces is obviously equivalent to the original variational problem (14).Indeed, the constraint ε = 0 enforces the enhanced strain to vanish and we are left with the problem (14).
The constrained optimization problem (17) may be transformed into a saddle point problem by introducing a suitable Lagrangian multiplier  in the space Here, the minus sign is a convention, s.t. the KKT conditions associated to the saddle-point problem attain the form div (⋅,  + ∇ s u + ε) = 0, and the Lagrange multiplier  equals the local stress field.

Assumed strain methods in micromechanics
Following Simo-Rifai, 108 the starting point of assumed strain methods is the three-field variational principle (18) for the triple and the spaces Assumed strain methods are obtained by choosing suitable finite-dimensional subspaces corresponding to a discretization of the domain Y into finite-elements {Y e } n elm e=1 .More precisely, Simo-Rifai 108 seek a triple In their analysis of assumed strain methods, Simo-Rifai 108 proposed three conditions which cover a variety of previously defined assumed strain, enhanced strain and B-bar methods, yet still admit a unified discussion.

Condition (i)
The spaces  h and h have trivial intersection in the sense Condition (ii) The spaces h and  h are L 2 -orthogonal, that is, the condition Condition (iii) The space  h contains element-wise homogeneous stresses, that is, the inclusion holds for arbitrary  e ∈ Sym(d) and the characteristic function  e of the e-th element The first condition (26) encodes the intuitive fact that the enhanced strains εh ∈ h have vanishing compatible part, that is, the equation implies u h = 0 and εh = 0.
Combined with the third condition (28), the second condition (27) ensures that the enhanced strains have vanishing mean on each finite element Y e , This is readily seen by choosing a stress field (28) which is non-zero only in a single element when evaluating the condition (27).Thus, as the discretization parameter h goes to zero, the Equation ( 31) encodes, in a sense, that the enhanced strain εh goes to zero, as well.
Due to the second condition ( 27), the three-field saddle-point problem (25) simplifies to a two-field variational problem The corresponding Euler-Lagrange equation reads Under the conditions (i)-(iii) and the assumptions ( 8) and ( 9) on the constitutive behavior, existence, and uniqueness of solutions to the problem (33) can be shown.Moreover, the convergence estimate holds for a constant C which depends only on the material constants  ± .Here, we introduced the stress field  = (⋅,  + ∇ s u h + εh ) and its element-wise average For details on deriving the inequality (34), we refer to Appendix B. Suppose the subspace Then, the first term of the right-hand side of the Strang estimate (34) becomes infinitesimal.The second term goes to zero, as well.This is a direct consequence of Lebesgue's dominated convergence theorem provided the stress  is square-integrable (which is true by construction of the stress operator ).Thus, qualitative convergence of assumed strain methods in micromechanics is established.
To provide quantitative converge results, we follow the discussion in the literature 114,115 and assume conditions on the microstructure and the material behavior, s.t.natural higher regularity results hold for the displacement field u.More precisely, let us assume that the cell Y is decomposed into a finite number of non-overlapping subdomains Ω a (a = 1, … , K), s.t. the stress operator  is homogeneous as well as sufficiently smooth on each domain and the interface I is sufficiently regular.Then, we assume that the solution u ∈ H 1 # (Y ; R d ) has additional regularity in the sense that the associated strain field ∇ s u is essentially bounded, that is, the finiteness of norm condition holds true, and that the restriction of the displacement field u to each of the homogeneous domains Ω a has square-integrable second-order (weak) derivatives, The higher regularity conditions ( 37) and ( 38) hold for heterogeneous linear elastic solids with sufficiently smooth interfaces, 116 see Schneider 114 for a discussion.For nonlinear solids, we are not aware of a suitable generalization, and we will treat the conditions ( 37) and (38) as assumptions.
Suppose that the finite-element space  h has the first-order approximation property, that is, there is a constant C s.t. for every w ∈ H 2 (Y ; R d ), the estimate inf holds.Then, the estimate inf holds for any field u ∈ H 1 # (Y ; R d ) satisfying the additional regularity assumptions ( 37) and ( 38) for a generic constant C.This is readily seen by a smoothing argument involving a function  h which is zero on an h-neighborhood of the interface I, see Schneider 114 for details.
With this discussion at hand, we conclude for the first term in the estimate (34)   inf for a generic constant C independent of h.
To handle the second term, we distinguish two cases.If the finite element Y e is fully contained in a homogeneous domain Ω a , the Poincaré inequality implies with a constant independent of h in case of a quasi-uniform triangulation.Otherwise, the (trivial) estimate applies.Thus, we obtain the estimate As the interface I is a codimension one surface in the cell Y , we have the count estimate which implies the bound All in all, we are led to the convergence behavior Thus, the assumed strain discretization (25) converges with the same rate as the traditional Galerkin discretization 70,71,115 with first-order accurate ansatz functions or the Moulinec-Suquet discretization. 31,32,114ast convergence for the effective properties may be established as usual 114,117,118 and under the assumptions stated in these works, where holds for the solution (u, εh ) to the problem (33).

Laminate composite voxels as an assumed strain method
Up to this point, general assumed strain methods in computational micromechanics were analyzed.This section is devoted to establishing the connection to the composite voxel method, 86,87,92 more precisely laminate composite voxels.Laminate materials 90( §9) are special microstructures which are layered in a direction n, a unit vector in R d .Due to the translation invariance of such a laminate in directions orthogonal to the lamination direction n, both the kinematic compatibility and the equilibrium condition, which constitute partial differential equations, in general, reduce to ordinary differential equations.§9) Also, for nonlinear and inelastic constitutive laws, the effective properties of laminate materials may be computed with relative ease, 91,92 again due to the intrinsic single dimensionality.This fact has been used by Kabel and co-workers to study laminate composite voxels in the context of voxel-based micromechanics.As illustrated in Figure 2 for a two-phase microstructure, whenever the analytic description of the microstructure phases is available, there are non-homogeneous voxel elements, that is, elements which contain more than a single phase of the microstructure.Kabel et al. 87 proposed to furnish these composite voxels with the effective properties of an equivalent laminate, that is, a microstructure layered in a direction which approximates the normal to the interface between the materials and the appropriate cut volume fractions.
In its original form, 87 the composite voxel methodology is part of material modeling, as the discrete versions of both the kinematic compatibility and the equilibrium equation remain untouched.Rather, the constitutive law of composite voxels is altered compared to a homogeneous voxel on either side of the interface.However, there is no inherent reason for selecting the effective properties of an appropriate laminate as the surrogate material model for the composite voxel.Indeed, there is no layering or periodicity of the "microstructure" within the voxel.

F I G U R E 2
Illustration of a two-phase material with smooth interface I on a background regular grid and magnification of a single element.
In this section, we discuss an alternative point of view leading to "laminate composite voxels."Let Y e be a composite voxel, as illustrated on the right hand side of Figure 2. Denote by the volume fractions of the composite domains Ω ± ∩ Y e ⊆ Y e and let n e ∈ S d−1 be a unit vector representing an approximation of normal to the interface between the phases Ω ± inside the voxel Y e .Denote by G e the tensor fields defined by where ⊗ s denotes the symmetrized dyadic product.Before discussing the inherent properties let us mention that the fields G e are to be used as enhanced strains in the assumed strain method (33), together with element-wise constant Lagrange parameters see Equation (25).Please notice that we will tacitly assume G e ≡ 0 and a e = 0 in homogeneous, that is, non-composite, voxels.
The tensor field G e (x) defined in Equation ( 52) is concentrated in the element Y e , that is, vanishes outside of this voxel.Moreover, for a non-vanishing vector a e , the strain field G e (x) ⋅ a is discontinuous across the interface I, mimicking the strain discontinuity across a material discontinuity of the continuous solution to the field Equation (12).This discontinuity also shows that Simo-Rifai's first condition (26) is satisfied as long as a conforming, in particular continuous, ansatz for the displacement field u is used.Indeed, the Equation ( 30) can only be satisfied in the trivial case, as restricted to the element Y e , the Equation ( 55) involves a continuous field on the left-hand side and a discontinuous field on the right-hand side, unless a e = 0. Simo-Rifai's third condition ( 28) is satisfied by construction (54), whereas the second condition ( 27) is implied by the voxel-wise vanishing mean value of the fields G e , that is, as a consequence of the identity The latter argument also reveals the reason for choosing the prefactors  ± e with opposite signs in the definition (52).After defining the enhanced strain space (53), we evaluate the equation governing the mechanical behavior of the assumed strain method.The equilibrium Equation ( 33) specializes for the case (53) at hand with w h = 0 as well as γh = G e ⋅ b e for arbitrary b e ∈ R d to the condition where we used that the field G e is supported on the element Y e .Decomposing the latter equation into the Ω ± -parts, we obtain the equation Eliminating the common prefactor we are thus led to the condition equivalent to the assumed part of the equilibrium equation (33) and which encodes the continuity of the (averaged) normal stresses inside the voxel Y e .Thus, the kinematic assumption (53) of a strain-field jump along the normal direction n e , accounting for the zero-average condition (31), directly leads to the laminate-type problem (62).
A convenient computational resolution of the problem (33) proceeds via a Schur complement approach, as suggested by Simo-Rifai. 108First, the Equation ( 62) is solved for the vector a e in each element, treating the compatible strain ∇ s u h as a parameter.Subsequently, the identified vector a e is regarded as a function of the compatible strain ∇ s u h and inserted into the original problem (33), leading to a statically condensed problem for the displacement field u h ∈  h .§3)

Use of level-set and signed distance functions
Classically, the composite voxel method 86,87 uses a much finer voxel background mesh to compute the quantities required for the method, that is, the cut volume fractions and the normal.The original works 87,99 also considered the option where the geometry is given by computer-aided design (CAD) data and established an appropriate workflow.However, the CAD strategy turned out to be unfavorable.Indeed, the subvoxel-based strategy has a natural link with microstructures given as (high-fidelity) voxel data, for example, as a result of micro-computed tomography scans, 18,19 where extracting explicit CAD interfaces is non-trivial.However, we saw in Section 1.2 that the strategy introduced by Merkert et al. 99 to estimate the interface normal based on connecting the center of mass of one of the phases within the voxel and the center of the voxel does not converge upon mesh refinement.Put differently, this method comes with a systematic error that is not easily eliminated, compare Keshav et al. 95 The work at hand exploits an alternative approach to extract the data necessary for the composite voxel method which, at the same time, yields accurate volume fraction as well as normal data and is naturally compatible to digital images.We consider the case that the microstructure Y is decomposed into two phases Y ± .More phases can be handled similarly, and we restrict to two-phase media essentially for the purpose of exposition.We suppose that a (sufficiently regular) level-set function 96 is given, s.t. the two phases Y − and Y + arise as the sub-and superlevel of the level-set function L at level zero whereas the interface I between the phases is encoded as the zero iso-contour of the level-set function L Many interesting properties of the phases Y ± and the interface I can be readily computed from the level-set function (63).However, the description ( 64) is highly non-unique and the computational treatment may suffer from an inauspicious choice of the level-set function L. Therefore, it is convenient to assume that the level-set function is actually normalized appropriately.More precisely, we consider level-set functions which arise as signed distance functions of the interface 96 where dist Y denotes the periodic distance between two points on the rectangular cell Y (5).
If the interface I is piecewise smooth, the signed distance function is differentiable almost everywhere, and its gradient satisfies the eikonal equation If the interface is twice continuously differentiable, the distance function dist Y is continuously differentiable close to the interface I 119(Lemma 14.16) and satisfies the equation where n ∶ I → R d denotes the out-ward pointing unit normal.In particular, the signed distance function (66) may be regarded as a function whose gradient is normalized and serves as a differentiable extension of the unit-normal vector field on the interface.With using the assumed strain method discussed in Sections 2.2 and 2.3 in mind, we are actually interested in an approximation of the continuous level-set function (66) in terms of Lagrangian finite-element trilinear shape functions {N A } n e nodes A=1 which are also used to represent the displacement field u h , that is, a level-set function of the form Such a function is completely determined by the values {L A } at the nodes {x A }. Therefore, in case the interface I in Equation ( 65) is known analytically, one may define the discretized level-set function ( 69) by point-wise evaluation As these values need to be computed only once in preprocessing, the cost of evaluating the level-set function ( 66) is less important.For the microstructures studied in this article, see Section 4, the interpolation procedure (70) is feasible.More precisely, we consider rotated laminates and inclusion-matrix composites with spherical and cylindrical inclusions.In these cases, the signed distance to the interfaces may be computed readily.
][122] Last but not least let us discuss how to store composite voxel data.There is a simple method to detect a composite voxel Y e based on the eight level-set values L A at the eight adjacent nodes.Indeed, in case the signs of these eight values are either all positive or all negative, the voxel belongs entirely to phase Y + or Y − , respectively.Otherwise, the voxel is a composite voxel (where we include the tautological cases where one level-set value is identically zero, as well).As the composite voxel method (62) enhances composite voxels only, data associated to pure voxels can be discarded.
There are two practical ways to proceed.Either the eight level-set values of the adjacent nodes are stored for each composite voxel, for example, in an array, or the appropriate processed data like the volume fraction and the normal is stored for each voxel.Both strategies have their merits and shortcomings, whereas the latter is less memory-demanding, it is also less flexible, as only derived quantities are stored.
The remainder of this section is devoted to describing techniques for extracting a normal vector n e and the volume fractions of the cut volumes Y ± e for each composite voxel Y e .

Extracting normal data and the best-approximating linear interface
In this section, we report several numerical strategies to calculate a unit normal for a fixed composite voxel to be used in the composite voxel method (62).These strategies come with different computational efforts, and their merits and limitations will be closer investigated in Section 4.

Averaging
The first and most simple strategy uses the fact that the gradient of the signed distance function L is normalized to unity see Equation (67), and extends the inward-pointing unit normal see Equation ( 68), provided the interface I is sufficiently smooth.As the discrete level-set function ( 69) serves as an approximation to the actual signed-distance function L, a simple approximation reads with the average As we deal with regular voxel grids and trilinear shape functions, the average gradient g e may also be computed by evaluating the gradient ∇L h at the center x e c of the voxel Y e .

Regression
Suppose that, for a fixed composite voxel, the level-set description L h is given in terms of trilinear ansatz functions and the eight nodal values {L A } of the adjacent nodes.Then, the interface I within the voxel is nonlinear, in general.In particular, there is no unique and "natural" normal vector associated to the interface.Despite the nonlinearity, when restricted to the 12 edges of the voxel, the level-set function L h is linear.In particular, computing the roots of the level-set function L h on the 12 edges is performed with ease.Excluding pathological cases, there can be between three and six intersection points x 1 , x 2 , … , x n , in general, see Figure 3 for an illustration.Three intersection points determine the values of the linear level-set function with x 0 ∈ Y e and a unit normal n uniquely.For more intersection points, we wish to determine the best-fitting linear level-set function ( 75) by regression, that is, by minimizing the quadratic least-squares fit objective function Differentiating the function F w.r.t.
that is, we may choose the centroid of the point cloud {x i } as the point x 0 .With this insight at hand, we are led to minimizing the objective function ( 76) to be minimized over the set of vectors n with unit norm.Rewriting the function (79) in tensorial form with the covariance matrix enables us to solve the problem (79) with ease.Indeed, by the expression (80), the minimizing direction n corresponds to the smallest eigenvalue of the symmetric and positive (semi-)definite second-order tensor M.
To sum up, to identify the coefficients of the linear level-set function (75), we first compute the intersection points of the given trilinear level-set function L h on the edges, calculate the centroid (78) as well as the covariance matrix (81), then compute an eigenvalue decomposition of the latter.The unit-normal vector n e of the composite voxel then arises as the direction u i corresponding to the smallest eigenvalue  i .It might be necessary to change the sign of the normal to keep the sub-and super-levelsets of the function L lin consistent with the trilinear function L h .

Minimax
The previous method was based on the identified intersection points x 1 , x 2 , … , x n with the voxel edges and used a plane regression for identifying an appropriate linear interface.For general intersection points, none of these intersection points will actually lie exactly on the identified plane.In particular, denoting the intersection points of the linear interface with the edges by {x lin i }, the sets are disjoint, in general.In their analysis of an immersed finite element method, Guo and Lin 109 faced a similar problem and introduced a different strategy.More precisely, they noticed that three points in general position determine a plane uniquely.Thus, they considered selecting three of the intersection points {x 1 , x 2 , … , x n }, and choosing the corresponding linear interface I lin e .Based on the three chosen points, a linear signed-distance function with appropriate unit-normal n e is readily identified.
Guo and Lin 109 were interested in mathematical analysis and sought to select the three intersection points judiciously to avoid mesh intersections with bad interior angles.Following ideas by Chen, 123 they selected those three intersection points whose maximum internal angle in the determined triangle is the smallest.We call this strategy minimax.
The differences and similarities of the introduced strategies are illustrated for a cylindrical inclusion oriented along the space diagonal (1, 1, 1)∕ √ 3 in Figure 4. Notice that both regression and minimax are able to retain the planar interface at the faces of the cylinder, in contrast to averaging.Both averaging and regression lead to more or less symmetric configurations w.r.t. the symmetry axes.In contrast, the results of minimax look more like "scales" which protrude from the surface.

Computing cut volume fractions
For the composite voxel method (62) to work it is necessary to compute the cut volume fractions ( 50) Restricting to spatial dimension d = 3 we assume the composite voxel to coincide with the unit cube s.t.we need not worry about the prefactor in their definition (84).As we assume that the voxel Y e is fully occupied by the materials Y ± e , that is, the identity holds, it actually suffices to compute either  + e or  − e .For a given level-set function L h , it is possible to approximate the volume fraction via quadrature on a subvoxel grid m × m × m via the counting of elements symbol #, see Merkert et al. 99 The strategy ( 88) is straightforward to implement and robust.Moreover, it converges to the correct volume fraction with the rate 1∕M.In particular, reaching high fidelity in the calculated volume fractions is not computationally feasible with this approach.Indeed, to reach an accuracy of 10 −5 , on the order of 10 15 operations need to be carried out.
We consider an alternative based on the linear approximation of the interface introduced in the previous Section 3.2.In this case, the cut volumes Y ± e form convex polyhedra, that is, the finite intersection of half-spaces, see Figure 5A.As the volume fractions  ± e sum to unity, it suffices to compute the volume of one of the two polyhedra Y ± e .We choose the one with fewer corners, see Figure 5A, where the selected polyhedron is framed in red.
There is a number of algorithms for computing the volume of general polyhedra, which is generally considered to be a hard problem. 124However, for the problem at hand, the number of possible polyhedron topologies is rather limited, see Figure 3, and an efficient treatment is possible.Of special concern is keeping numerical errors small. 125,126For this purpose, we use the numerically accurate and efficient strategy introduced by Mirtich 110 which concerns the integration of arbitrary polynomials over general polyhedra.The basic idea is to reduce volume integrals to surface integrals by the divergence theorem.Subsequently, each surface integral may be reduced to a series of line integrals by Stokes' theorem.Evaluating the resulting line integrals with polynomial integrands is performed analytically.
To add a little more detail, we consider the polygon Ω P inside the unit cube and introduce the vector field in the standard Euclidean frame with coordinates (x, y, z) which enables us to invoke the divergence theorem to write where n denotes the outward-pointing unit normal on the polyhedron boundary Ω P .The latter boundary is decomposed of planar faces  with constant normals n, that is, we have the identity The outlined strategy reduces the task of computing the volume of the polyhedron to calculating the integral of the x-coordinate over the polygons  constituting the boundary of the polyhedron Ω P .Mirtich proposed to actually use the dimension-reduction strategy once again.To do so, one first needs to rewrite the surface integral as a two-dimensional integral in parameterized form on a suitable parametric domain Π and a parametrization  ∶ Π → R 3 .In case the projection along the z-axis has positive area, one may choose the projection for the linear interface described by the equation and the parametrization see Figure 5B.In this case, the surface integral (92) becomes where we used the transformation F I G U R E 5 Geometrical considerations for Mirtich's formulae. 110(A) Schematic of the polyhedron whose volume is to be determined.
(B) Illustration of the projection process required to represent the interface as a graph.
a consequence of the unit length of the face normal-vector.Mirtich proposed to rewrite the integral ( 96) by Green's theorem which may be further decomposed into individual linear segments.Some care has to be taken when using the representation (96) for numerical computations, as the z-component of the normal vector may become arbitrarily small.In these cases, Mirtich 110 suggested to use projections along either the y-or the z-axis to retain numerical stability, leading to corresponding equivalents of the formula (98).More precisely, the axis to project onto is selected based on the norm of the normal component which is largest, see Figure 5B.

Setup
We implemented the level-set based approach to composite voxel methods, as discussed in Section 3, into an existing in-house FFT-based micromechanics solver.The solver was written in Python with Cython extensions and an OpenMP parallelization.Cython is used for performance reasons. 127e used Willot's discretization, 60 which corresponds to trilinear finite elements with one-point quadrature, 70 and the linear conjugate gradient method 42 to resolve the balance of linear momentum to a tolerance of 10 −5 with the natural convergence criterion discussed in Schneider. 39We implemented a Newton method to determine the optimal level-set cutoff value to reach to analytically known targeted volume fractions of each phase to an accuracy of 10 −5 .The simulations were performed on a desktop computer with 64 GB RAM and 20 cores with 3.6 GHz each except for the reference fiber computation where a workstation featuring two AMD EPYC 7642 with 48 physical cores each was put to use.The employed isotropic and linear elastic material parameters are listed in Table 1.

Laminate structure
§9) More precisely, we consider a laminate material which is exactly represented on a periodic cell Y but whose lamination direction does not coincide with the axes of periodicity of the cell Y .We use the lamination direction n ≡ (1∕ √ 10, −3∕ √ 10, 0) and a period 16 Then, any such laminate may be exactly represented on a cubic volume Y = [0, 16] 3 (m) 3 .For the study at hand, we chose equal volume fractions.
As the level-set function of a laminate is linear in each voxel, all the three strategies presented in Section 3.2 recover the laminate interface exactly, as shown in Figure 6A.Also, the volume fractions are computed exactly by the Mirtich formulae presented in Section 3.3, up to numerical precision.
We start by investigating the quality of the normal estimated by the subvoxel-based formula (4), that is, where the normal n m e is given by the direction connecting the centroid of one phase with the centroid of the voxel.For each composite voxel Y e , we consider the following error measure which quantifies the deviation between the exact normal vector n e and the subvoxel approximation n m e with m 3 subvoxels.

TA B L E 1
Material parameters considered for the computational experiments. 128bers The quantity (99) measures the Frobenius norm of the difference of the orthogonal projectors onto the one-dimensional spaces defined by the respective normal vectors.Elementary algebraic transformations lead to the identity in terms of the angle  m between the unit vectors n e and n m e , implicitly defined via Recent work 95 introduced an improved subvoxel technique.More precisely, to reduce the inaccuracies inherent to the original formulation (4), Keshav et al. 95 introduced a weighted least-squares regression at a subvoxel scale.In detail, they compute the Laplacian of the discrete subvoxel image to determine the interface subvoxels.Then, for each interface subvoxel center, a difference vector is computed to the voxel center, similar to the subvoxel-based formula (4).The normal of the voxel is then obtained by a least-squares regression of these interface subvoxels with the absolute values of the Laplacian as weights.We implemented their approach 129 and refer to this strategy as composite boxels (ComBo) throughout this subsection.For a resolution with 16 3 voxels, the distribution of the error is shown in Figure 7A.For 4 3 subvoxels, the normals are systematically false with a relative error around 40%.This inaccuracy is also reflected when inspecting the linear interfaces in each composite voxel, see Figure 6B for an illustration.To create the images, we set up a linear level-set function in each composite voxel whose gradient matches the calculated normal n m e and whose cut volume fractions equal the quadrature-based estimates (88).Notice that such a voxel-based level-set function may lead to discontinuities of the interface.
Increasing the subvoxel count to 8 3 reduces the error to slightly more than 30% on average and decreases the standard deviation significantly.This improvement is also reflected in the geometry, see Figure 6B.Higher subvoxel counts further decrease the error.However, the error does not go to zero as the number of subvoxels goes to infinity.In view of the example discussed in Section 1.2, this observation does not come as a surprise.For ComBo, a similar behavior is observed.The mean absolute error of the estimated normals and ComBo is approximately 15% lower for all considered subvoxel counts.The standard deviation remains at a similar level to that of the subvoxel approach.On the visual side, there is not much difference between 4 3 and 8 3 subvoxels, see Figure 6C,D, respectively.Notice, however that the lack of convergence of the estimated normals is also reflected in the visualizations, that is, via the lack of continuity of the interface.
The statistics of the estimated normals is only computed based on the composite voxels, that is, the plain voxels do not enter into the consideration.Moreover, the fraction of composite voxels among all voxels gets smaller upon grid refinement, that is, their influence on the computed effective properties decreases, as well.Therefore, we would like to assess the influence of these inaccurate normals.For this purpose, we consider the effective elastic properties to be a reasonable quantity of interest.More precisely, we furnish the two phases of the laminate with the material parameters of glass and polyamide, see Table 1, and record the error in the normal stress in x-direction in response to a strain loading  = e x ⊗ e x in Figure 7B.Here, the analytic solution 90( §9) serves as the reference.
We notice that the effective stress converges linearly for all considered scenarios, supporting the result (48).Using no composite voxels at all-referred to as none in the legend-leads to an error above 2% even in the highest resolution.For exact normal and volume-fraction data, shown in red, the relative error is slightly larger than 2% for N = 16 voxels per edge.Less than 1% error is reached for N = 64 voxels.The subvoxel-based composite voxels lead to a consistently larger relative error.For 4 3 subvoxels, the relative error is about 9% on the coarsest resolution and below 0.6% on the finest resolution.Increasing the subvoxel count consistently leads to a more accurate solution.However, the difference between 8 3 and 16 3 subvoxels is rather small.Notice that an error below 1% is reached for N = 128 voxels per edge.An increase in the subvoxel count to 32 3 doesn't lead to a significantly better solution.For this study, the composite voxels based on exact normal data are about as accurate as subvoxel-based composite voxels on the double resolution per edge, that is, on the eightfold total voxel count.The relative error of the ComBo approach is generally in between the subvoxel method and the exact normal method.For coarser resolutions, the improvement to the subvoxel method is rather small.However, ComBo achieves the same level of accuracy as the exact normal method for finer resolution.At the highest resolution, the relative errors are indistinguishable.A smaller improvement is also noticed for increasing subvoxel count compared to using the subvoxel method.

Hashin's coated sphere
We continue our investigations on the performance of the approximation techniques for both the normal vector and the volume fractions within composite voxels introduced in Section 3. We again consider a scenario with an available analytic solution for the effective properties.However, in case of Hashin's neutral inclusion, 130 the interfaces are non-planar.Hashin's coated sphere consists of three phases, a spherical core an annular region and the remaining "matrix" with the geometric center x c of the microstructure Y = [0, L] 3 and two radii r 1 < r 2 < L. Here, dist Y denotes the periodic distance on the cell Y .The microstructure is illustrated in Figure 8. Endowing each of these three phases with isotropic elastic properties, Hashin 130 demonstrated that under macroscopic compression it is possible to select the elastic parameters of the core Ω 1 and the shell Ω 2 in such a way that the effective compression modulus of the microstructure coincides with the compression modulus of the matrix phase Ω 3 .The shear moduli  1 and  3 of both the core and the matrix do not enter into these considerations.For the investigation at hand, we fix the compression modulus K 3 ≡ K eff of the matrix and the ratio  2 = 3 K 2 ∕5 of the intermediate region.There is an interconnection 130 between the compression moduli K 1 and K 2 , that is, where We subsequently assume that the compression modulus K 1 exceeds K 3 ≡ K eff .Then, the parameter K 2 is smaller than K 3 ≡ K eff .For increasing K 1 , K 2 decreases.As K 1 → ∞, K 2 remains bounded away from zero.Thus, in the limit K 1 → ∞, we are faced with a rigid inclusion with special coating inside the matrix.Then, we may consider different material contrasts in the compression modulus.
For the study at hand, We use the radii where e = ∑ ∞ k=0 1 k! is Euler's number.The resulting volume fractions  1 and  2 are irrational, for example,  1 = (4 ∕(7 e)) 3 , avoiding pathological scenarios where the interfaces intersect the nodes frequently.
F I G U R E 8 Schematic cut of Hashin' coated sphere.
We investigated the performance of different precomputing techniques, as introduced in Section 3, for the composite voxel method for Hashin's coated sphere.
The results for varying resolution and material contrast (106) are recorded in Figure 9.For a comparatively small material contrast  = 10, shown in Figure 9A, we observe that using no composite voxels at all-referred to as none in the legend-leads to a linear convergence of the effective compression modulus.Moreover, even for the coarsest resolution, the relative error is well below 1%.Having established this setup devoid of composite voxels as our point of reference, we move on to discuss the different preprocessing techniques available for composite voxels.We start with the classical subvoxel-based strategy for estimating the normal (4) and the volume fraction (88).We use m = 4 subvoxels in each of the three coordinate directions as it represents a convenient trade-off between accuracy and computational speed in preprocessing.We notice that using this type of composite voxels consistently improves upon using no composite voxels at all.For the finer resolutions starting at 64 3 , the gain in accuracy is between a quarter and half an order of accuracy compared to using no composite voxels at all.
Taking a look at the three strategies introduced in Section 3.2, we observe that they all lead to a comparable performance, except for the coarsest resolution 16 3 considered.At this coarsest resolution, the averaging strategy (74) performs worst, that is, the resulting error exceeds its subvoxel-based pendant.The two other strategies regression (79) and minimax, see Figure 4C, lead to a slightly lower error on the coarsest resolution.For higher resolution, using each of the three strategies averaging, regression or minimax consistently improves upon the subvoxel-based composite voxels.For some resolutions, like 64 3 voxels, the improvement is small, but for others, like 128 3 voxels, half an order of magnitude can be gained in terms of accuracy.All in all, compared to using no composite voxels at all, these composite voxel strategies lead to an improvement between a half and a full order of magnitude.Thus, the novel composite voxel strategies are able to improve upon the subvoxel-based composite voxels.Indeed, the latter strategy may lead to a seemingly unreasonable lack of improvement for some resolutions, like at 32 3 voxels, which are probably caused by inaccurately computed interface normals.
Taking a look at the higher contrast  = 100, shown in Figure 9B, these observations are essentially confirmed.However, there are small differences.For a start the total error level is about half an order of magnitude larger than for a contrast of 10.Moreover, at the coarsest resolution, we observe that using no composite voxels at all (none) is better than using subvoxel-based composite voxels.Using the regression and minimax strategies turns out to be favorable and the biggest improvement is provided by the averaging strategy.
Perhaps more interesting is the case with a contrast of 1000, shown in Figure 9C.We observe that using the subvoxel-based estimates for the volume fractions and the normal within the composite voxels does not lead to an improvement of the error compared to using no composite voxels at all.Similar observations were reported previously. 87n these works, the lack of performance of laminate-based composite voxels were attributed to the possible degeneracy of laminates for infinite material contrast.In the work at hand, we observe that the failure to improve upon using only plain voxels is more likely rooted in the use of incorrect normal data.Indeed, any of the three introduced and more advanced normal-estimation techniques, see Section 3, leads to an improvement in the relative compression modulus compared to using no composite voxels.
To sum up, the investigations on Hashin's coated sphere, see Figure 8, revealed that using laminate-based composite voxels with either of the novel normal-estimation techniques provided in Section 3 improves the quality of the computed effective properties throughout the entire considered range of contrasts and resolutions.Moreover, some unreasonable performance lows of the subvoxel-based composite voxels could be removed which hints at the influence of incorrectly estimated normal data.

Fiber reinforced composite
Our final computational example is concerned with a short-fiber reinforced composite microstructure.More precisely, we consider fibers with a diameter of 10m and a fiber length of 200m, that is, an aspect ratio of 20, reinforcing a matrix at 20% volume fraction.We generated a periodic fiber microstructure with the SAM algorithm 131 in a cubic unit cell with 256m edge length and a fiber orientation described by the second-order fiber-orientation tensor w.r.t. the coordinate system of the unit cell, employing the exact closure. 132,133Thus, the fiber orientation is transversely isotropic w.r.t. the e 1 -axis, and the majority of the reinforcing fibers is pointing in this direction.We set the minimum inter-fiber distance to 4.5 m to avoid triple composite voxels, that is, more than two fibers intersecting a voxel, on the coarsest considered resolution of 64 3 voxels.The generated microstructure comprises 213 fibers and is shown in Figure 10 in all three coordinate planes.The polyamide matrix and the E-glass fibers were furnished with the material parameters listed in Table 1.
The results of the different composite voxel strategies are shown in Figure 11 and compared to using only plain, that is, no composite, voxels.As no analytical solution for the effective properties is available for this microstructure, we use a high-fidelity computation on 1024 3 voxels with regressed composite voxels as our reference when computing relative errors.Taking a look at the relative errors for the full stiffness tensor, see Figure 11A, where we use the Frobenius norm of the Mandel representation, we observe that all considered discretizations converge linearly, that is, support the result (48).There is a clear ordering between the composite voxel strategies.Averaging performs worst, and the subvoxel strategy (with 4 3 subvoxels) is only slightly better.The two other strategies for estimating the normal, regression and minimax, are close to each other and turn out to be consistently better than their competitors.The difference is largest at the coarsest resolution, 64 3 , at about half an order of magnitude, and shrinks to about a quarter order of magnitude for the finest considered discretization at 512 3 .Interestingly, the averaging and subvoxel strategies perform worse than the strategy avoiding composite voxels at all for the two coarsest discretizations considered.Supposedly, this is a consequence of the coarse resolution of the fibers.Still, the question arises why averaging and subvoxel suffer from such a bad performance for the coarse resolution in this scenario, whereas the other two strategies, regression and minimax, appear to work well, that is, they consistently outperform using no composite voxels even for a coarse resolution.In contrast to the microstructures considered in this computational examples section, the discontinuous fibers considered for the example at hand are cylindrical and naturally come with a signed distance function which has a kink at the rims of the cylinders' ends.
In particular, the signed distance function is not continuously differentiable, and the interpolated level-set function may not serve as an appropriate surrogate for the continuous level set function when quadrature-based estimates for the normal vector (4) are employed.The strategies regression and minimax, on the other hand, are based on linear interpolation along the voxel edges of exact values of the signed distance function evaluated on the voxel corners and appear not to suffer from the same handicaps.These differences mostly manifest for the two coarsest considered discretizations.At 256 3 voxels, all considered discretization scenarios match closely.Such a resolution is also required to ensure a relative error below 1%.For the 512 3 example, all composite voxel strategies outperform using no composite voxels.
We investigated the full stiffness matrix first because it serves as the primary quantity of interest to engineering practice when considering the effective elastic properties of such short-fiber composites.However, as the stiffness in different directions may differ strongly, such a combined measure does not permit to distinguish the approximation quality in the various directions.For instance, one expects the stiffness in the e 1 -direction to be higher than in the transverse directions.Depending on the context, it may be of interest to assess the quality of the approximation in these transverse directions, as well.
We start by investigating the relative error in the e 1 -direction, see Figure 11B, where the majority of fibers is oriented.We measure the error of the entire effective stress tensor, again in the Euclidean norm of the Mandel representation.The performance of averaging remains suboptimal.More precisely, the errors for both averaging and none almost exactly match the full-stiffness case in Figure 11A, except for an outlier of none at the 512 3 resolution, where an unusually small relative error is encountered.This apparent advantage of the subvoxel strategies probably comes from the increased importance of the volume fractions for these scenarios.Indeed, for a unidirectional continuous fiber composite, the Voigt average 134 of the Young's moduli of fiber and matrix serves as a good estimate for the Young's modulus in this direction.The subvoxel-based estimate for the volume fraction, on the other hand, is actually a rather accurate one.
For this scenario, the minimax strategy turns out to be both highly accurate and reliable.Even for the coarsest considered resolution, a relative error of 1% is reached.The regression strategy performs similarly, but comes with a higher variance, for example, a larger error for 64 3 voxels and a much smaller error at 128 3 voxels.
The results for the transverse direction, shown in Figure 11C, actually match the full-stiffness case in Figure 11A rather well, and the same conclusions apply, at least for the different composite voxel strategies.Using no composite voxels turns out to be a bit better for this case and a coarse resolution, where it turned out to be a little worse in fiber direction.This effect may be caused by the higher stresses encountered when loading in fiber directions.
To sum up the findings in this example, the averaging strategy is not recommended for composites with non-smooth interfaces, whereas both regression and minimax lead to a consistent and reliable improvement over using no composite voxels.The advantage is about half an order of magnitude for the coarse resolution which roughly amounts to using only one eighth of the resolution required for a specific accuracy and avoiding the use of composite voxels.

CONCLUSION
The work at hand was devoted to understanding the origins of the laminate-type composite voxel strategy 108 which helps to conceptually understand composite voxels.Indeed, in the same way that first-order finite elements fail to represent a gradient in the strain accurately, interface-nonconforming discretizations have difficulties providing accurate fields in the vicinity of the interface.In this analogy, laminate composite voxels play the same role in enhancing accuracy in micromechanics as linear assumed strain fields for bending-dominated problems.
We provided a straightforward and thorough analysis of the assumed strain method for the heterogeneous case, providing optimal convergence estimates for a rather large class of assumed strain methods.In fact, it might be possible to design more powerful assumed strain methods which enhance the bulk accuracy in micromechanical finite-element models in addition to providing accuracy close to the interface as supplemented by composite voxels.Moreover, it would be desirable to account for constitutive laws at the interface in the framework of composite voxels, complementing the available approaches. 135,136e provided a number of strategies for computing approximations to both the cut volume fractions and the normals to the interfaces inside a composite voxel.As these calculations are performed in preprocessing and in an independent fashion for each composite voxel, ensuring numerical stability of the algorithms is actually more important than pure performance.Indeed, due to the complete decoupling of the calculations, it is possible to parallelize them in a more or less perfect manner.We identified Mirtich's approach to computing polyhedral volume as a suitable candidate for satisfying this criterion.
Our numerical findings revealed the limitations of subvoxel-based composite voxels.Of course, there are scenarios where only a subvoxel grid is available, and there is no alternative.However, for analytically given interfaces, level-set based composite voxel represent an alternative with superior accuracy properties.For instance, we saw that there are cases where subvoxel-based composite voxels do not improve upon using no composite voxels at all, in contrast to some level-set based composite voxels.This circumstance also supports the importance of accurate normal data when using composite voxels.
We found that the averaging scheme did not work that well for cylindrical inclusions, and we recommend either the regression or the minimax strategy.Interestingly, using these level-set based composite voxels did turn out to improve the quality of the computed effective properties consistently, and its use is thus firmly recommended.measures the consistency error both of the space W h and the operator g h .Indeed, the element u solves the equation which is different from  * h g h (u), unless a conforming Galerkin framework is used.Let us proceed with deriving the estimate (A8).For any element v h ∈ W h and the solution u h ∈ W h of the Equation (A7), we notice where we used the condition (A7) for w h = u h − v h .Thus, adding a zero, we obtain the identity By definition of the inclusion mapping (A9), we observe that is, the identity holds.In view of the monotonicity condition, the triangle inequality leads to the estimate Dispensing with the term ||u h − v h || common to both sides implies the inequality where we used the Lipschitz condition (A5) in the second line.The triangle inequality applied to the additive splitting implies the estimate where we made use of the previous result (A19).As the element v h is arbitrary, minimizing over the right hand side of the estimate (A21) yields the desired inequality (A8).

APPENDIX B. ANALYSIS OF ASSUMED STRAIN METHODS FOR HETEROGENEOUS MATERIALS
The purpose of this appendix is to show the estimate (34) for the solution (u h , εh ) ∈  h × h of the Equation ( 33) This operator is well-defined thanks to the conditions ( 9) and (11).Moreover, it is  − -strongly monotone and  + -Lipschitz continuous as a consequence of the conditions (8) and (9).Then, the problem find  ∈ W 0 , s.t.g()[] = 0 for all  ∈ W 0 (B7) has a unique solution , which is furthermore of the form  = ∇ s u for the unique displacement solution u of the original variational problem (14).
Similarly, we may encode the problem (B2) in the form find  h ∈ W h , s.t.g( h )[ h ] = 0 for all  h ∈ W h (B8) with the connection  h = ∇ s u h + εh for u h ∈  h and εh ∈ h .Due to the strong monotonicity and Lipschitz continuity, existence and uniqueness of solutions to the problem (B8) follow directly.Moreover, we are in the position to apply the theory developed in Appendix A with g = g h to get an estimate for the discretization error.More precisely, the inequality (A8) implies the estimate To treat the second term in the inequality (B9), we first notice the identity due to the decomposition and the validity of the momentum balance (12).Furthermore, the element-wise averaged stress field  h defined in Equation ( B3) is an element of the space  h by condition (iii), see the inclusion (28).Then, we have the identity of the second term in the Strang inequality (B9).In particular, the desired estimate (B1) emerges.

1
Explicit example for a systematic error when computing the normal n e with center of mass method (4).(A) Example 1. (B) Example 2 (Figure 1A normal in red).
Possible cases of a linear interface intersection a cube up to rotation and reflection.(A) Case 1. (B) Case 2. (C) Case 3. (D) Case 4. (E) Case 5.

F I G U R E 4
Illustration of the normal-estimation methods introduced in Section 3.2 for a cylinder.(A) Averaging.(B) Regression.(C) Minimax.

F I U R E 7
Error in the normal vector and the normal stress in x-direction for the laminate shown in Figure 6A.(A) Statistics of the error (99) in the estimated normal versus number of subvoxels for 16 3 voxels.(B) Relative error in effective normal stress in x-direction versus resolution.

F 11
I G U R E 10 Different views on the generated short-fiber microstructure resolved by 512 3 voxels.(A) x-direction.(B) y-direction.(C) z-direction.Relative errors of computed stress tensor and stresses in principal fiber direction and transverse to it for the short-fiber microstructure shown in Figure 10.(A) Effective stiffness tensor.(B) Stress in fiber direction.(C) Stress in transverse direction.
Y (⋅,  + ∇s u h + εh ) ∶ (∇ s w h + γh ) dx = 0 for all (w h ,  h ) ∈  h × h (B2)and the stresses = ( + ∇ s u) as well as  h = n elm ∑ e=1  e  e with  e = ⨍ Y e  dx (B3)associated to the continuum solution u ∈ H 1 # (Y ; R d ) of the micromechanical balance Equation(12).To do so, we will re-write the problem (B2) in the language of Appendix A, which contains a version of Strang's second lemma137(Theorem 4.2.2)   suitable for monotone operators.We are concerned with the spaceW = L 2 (Y ; Sym(d)),(B4)and the subspacesW 0 = ∇ s H 1 # (Y ; R d ) as well as W h = ∇ s  h ⊕ h , (B5) together with the operator g ∶ W → W ′ , g()[] = ⨍ Y ( + ) ∶  dx, ,  ∈ W. (B6) ⨍ Y  ∶ γh dx = ⨍ Y ( −  h ) ∶ γh dx,  h ∈ h , holds, where the mapping  * h ∈ L(W ′ , W ′ h ) is the adjoint of the canonical inclusion  h ∶ W h → W, that is, we have( * h )[w h ] = [w h ] for  ∈ W ′ and w h ∈ W h .(A9)Thus, the last term in the estimate (A8) may also be written in the form|| * h g h (u)|| W ′ h = sup w h ∈W h , ||w h || W h ≤1 g h (u)[w h ],(A10)which is typically preferred when stating the linear version of Strang's second lemma.137(Theorem4.2.2)The interpretation of Strang's estimate (A8) is the following.The first terminf v h ∈W h ||u − v h ||,(A11)quantifies the approximation error of the solution u with elements of the space W h .The second term