Level-set based design of Wang tiles for modelling complex microstructures

Microstructural geometry plays a critical role in a response of heterogeneous materials. Consequently, methods for generating microstructural samples are becoming an integral part of advanced numerical analyses. Here, we extend the unified framework of Sonon and co-workers, developed originally for generating particulate and foam-like microstructural geometries of Periodic Unit Cells, to non-periodic microstructural representations based on the formalism of Wang tiles. The formalism has been recently proposed as a generalization of the Periodic Unit Cell approach, enabling a fast synthesis of arbitrarily large, stochastic microstructural samples from a handful of domains with predefined microstructural compatibility constraints. However, a robust procedure capable of designing complex, three-dimensional, foam-like and cellular morphologies of Wang tiles has been missing. This contribution thus significantly broadens the applicability of the tiling concept. Since the original framework builds on a random sequential addition of particles enhanced with a level-set description, we first devise an analysis based on a connectivity graph of a tile set, resolving the question where should a particle be copied when it intersects a tile boundary. Next, we introduce several modifications to the original algorithm that are necessary to ensure microstructural compatibility in the generalized periodicity setting of Wang tiles. Having a universal procedure for generating tile morphologies at hand, we compare strictly aperiodic and stochastic sets of the same cardinality in terms of reducing the artificial periodicity in reconstructed microstructural samples, and demonstrate the superiority of the vertex-defined tile sets for two-dimensional problems. Finally, we illustrate the capabilities of the algorithm with two- and three-dimensional examples.



There are no comments yet.



Computational analysis of transport in three-dimensional heterogeneous materials

Porous and heterogeneous materials are found in many applications from c...

Force Analysis for Interactions beyond the Closest Neighbor in a Periodic Structure

Periodic structures are a type of metamaterials in which their physical ...

Handling boundary constraints for numerical optimization by particle swarm flying in periodic search space

The periodic mode is analyzed together with two conventional boundary ha...

An integral equation method for the simulation of doubly-periodic suspensions of rigid bodies in a shearing viscous flow

With rheology applications in mind, we present a fast solver for the tim...

Periodic Fast Multipole Method

A new scheme is presented for imposing periodic boundary conditions on u...

Fast multipole boundary element method for the acoustic analysis of finite periodic structures

In this work, two fast multipole boundary element formulations for the l...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Geometrical details of a material composition drive many macroscopic phenomena such as crack initiation and propagation (Nguyen et al., 2011) or meta-behaviour (Barchiesi et al., 2019). Advanced computational strategies thus tend to incorporate knowledge of material microstructures111We use the term “microstructure” in a broader sense, without referring to a specific scale length. to enhance their predictive power. Matouš et al. (2017)

When the micro and macro scales are well separated, the microstructural response is typically up-scaled through homogenization methods. Albeit computationally intensive, numerical homogenization (Matouš et al., 2017) now supersedes analytical methods (Torquato, 2002) thanks to its ability to handle complex microstructural geometries and non-linear problems. With increasing computational power, numerical models that fully resolve the microstructural geometry in the whole macroscopic domain, e.g. (Nguyen et al., 2017; Lloberas-Valls et al., 2011), or in regions of interest (Akbari Rahimabadi et al., 2015; Lloberas-Valls et al., 2012), emerge, addressing problems without clear scale separation.

In both cases, however, the accuracy of the strategies critically depends on the representativeness of the provided microstructural geometry, accenting the crucial role of microstructure modelling in multi-scale approaches.

Compared to materials with regular microstructures, characterized entirely by Periodic Unit Cells (PUCs), modelling random heterogeneous materials is more intricate, because any finite-size representation automatically implies information loss. The optimal microstructure representation should capture intrinsic randomness and fluctuations in a microstructure while remaining computationally tractable (for numerical homogenization) or inexpensive to construct (for fully resolved simulations).

Figure 1: Illustration of a synthesized sample of a closed-foam microstructure obtained as a tiling assembled from a set of Wang Cubes. Individual face codes are shown in semi-transparent colours.

1.1 State-of-the-art in modelling random microstructures

One of the widely adopted approaches to modelling heterogeneous materials rests on an extension of PUC generated such that its spatial statistics match that of a reference microstructure. This procedure appears in the literature under various names such as Statistically Optimal Representative Unit Cell (Lee et al., 2009), Repeating Unit Cell (Yang et al., 2018), Statistically Similar Representative Volume Element (Balzani et al., 2014), or Statistically Equivalent Periodic Unit Cell (SEPUC) (Zeman and Šejnoha, 2007), among others. The spatial statistics involved range from Minkowski functionals (Scheunemann et al., 2015)

to multi-point probability functions 

(Torquato, 2002), out of which the two-point probability (Jiao et al., 2007; Zeman and Šejnoha, 2001; Rozman and Utz, 2001), two-point cluster (Jiao et al., 2009), and lineal path (Lu and Torquato, 1992; Kumar et al., 2006; Zeman and Šejnoha, 2007; Havelka et al., 2016) functions are the most frequently used.

Following Povirk’s seminal work (Povirk, 1995), the majority of cell representations are generated using optimization procedures, minimizing the discrepancy between the statistical characterization of the reference microstructure and its compressed, PUC-like representation. The particular choice of optimization algorithm varies, including the simulated annealing (Yeong and Torquato, 1998; Kumar et al., 2006; Zeman and Šejnoha, 2007), genetic (Basanta et al., 2005; Kumar et al., 2008; Lee et al., 2009; Yang et al., 2018) and gradient (Povirk, 1995; Fullwood et al., 2008a), or phase-recovery (Fullwood et al., 2008b) algorithms.

The second approach to microstructure generation utilizes reference samples of the microstructure. New realizations are then obtained with a Markovian process, taking individual voxels (Liu and Shapiro, 2015) or a patch of voxels (Tahmasebi and Sahimi, 2013)

from the provided reference samples according to the proximity of the spatial statistics computed for their surroundings. Alternatively, searching for the statistics proximity can be replaced with a classification tree-based supervised learning model 

(Bostanabad et al., 2016).

The previous two approaches suffer from high computational cost related either to optimization or to training the learning model. Besides, the applicability of their outputs is sensitive to the spatial statistics considered, attesting to the ill-conditioning of microstructure reconstruction problem itself. Achieving a good match in selected statistics does not automatically guarantee similar overall behaviour; for instance, Biswal et al. (1999) demonstrated that realizations with similar two-point probability functions could differ significantly in percolation characteristics that govern overall transport properties.

Complementary to the statistics-informed methods, the third approach to generating microstructural realizations relies on meta-modelling of microstructure genesis. These methods range in complexity from the Monte-Carlo Potts (Saito and Enomoto, 1992) and phase field models (Krill III and Chen, 2002) of grain growth, to sedimentation-and-compaction models (Biswal et al., 1999), to various particle packing algorithms, e.g. (Falco et al., 2017; Alsayednoor and Harrison, 2016, and references therein), based on either Random Sequential Adsorption (RSA) (Cooper, 1988; Segurado and Llorca, 2002) or molecular dynamics (Lubachevsky and Stillinger, 1990; Donev et al., 2005; Ghossein and Lévesque, 2013).

The relevance of packing algorithms extends beyond simple particle-matrix microstructures because the resulting packings often serve as initial seeds for tessellation-based models applicable to polycrystals (De Giorgi et al., 2010; Chen et al., 2015), foams (Alsayednoor and Harrison, 2016), or cell tissues (Mebatsion et al., 2006; Chakraborty et al., 2013). Due to its straightforward implementation, Voronoï tessellation is the most common choice; however, the resulting geometry is oversimplified for many materials. For instance, Voronoï-based models overestimate overall stiffness for high porosity foams (Doškář and Novák, 2016). The curvature of the cell walls (De Giorgi et al., 2010; Simone and Gibson, 1998) and heterogeneity in cell size (Alsayednoor and Harrison, 2016; Chen et al., 2015) and wall thickness (Chen et al., 2015) must be additionally introduced to obtain realistic geometries. Similar effects can be achieved by modifying the distance measure used during tessellation, e.g. models based on the Laguerre variant generate microstructures with multi-mode cell size distribution (Chen et al., 2015; Alsayednoor and Harrison, 2016; Falco et al., 2017). Improving on this idea, Chakraborty et al. (2013) proposed Adaptive Quadratic Voronoï Tessellation, attributing a distinct non-isotropic metric to each seed and thus allowing for additional control over the resulting geometry.

Recently, Sonon et al. (2012) presented a method based on a level-set description of particles that significantly accelerates the standard RSA procedure (see discussion in (Sonon et al., 2012)). In comparison to other accelerations of RSA, e.g. the hierarchical BBox-based algorithm (Yang et al., 2018), the Sonon’s method readily facilitates generating complex microstructures using linear combinations of the nearest neighbour distance functions and dedicated morphing operations (Sonon et al., 2012; Sonon and Massart, 2013; Sonon et al., 2015; Massart et al., 2018). In a sense, this approach introduces non-isotropic pseudo-metrics222Level-set description based on the signed distance to the th nearest particle boundary does not fulfil all metric criteria. of (Chakraborty et al., 2013) in a geometrically-motivated way by considering arbitrarily-shaped particles.

Albeit faster than the original RSA or optimization-based approaches, the latter method still starts anew every time an additional realization is required, imposing overhead on, e.g., investigations of the Representative Volume Element (RVE) size that require multiple microstructural samples to be generated, see (Kanit et al., 2003; Gitman et al., 2007; Dirrenberger et al., 2014). Alternatively, larger microstructural realizations can be assembled from (SE)PUC; however, such construction introduces non-physical, long-range, periodic artefacts in the microstructural geometry and in its local response.

1.2 Wang tiling in microstructure modelling

(a) (b) (c)
Figure 2: Illustration of three abstract sets of the equal cardinality of 16 tiles: (a) an edge-based set with two horizontal and vertical codes, (b) a vertex-based set over four horizontal and vertical codes, and (c) the aperiodic Ammann’s set over six horizontal and vertical codes. The corresponding vertex codes—original (for (b)) and identified using the method outlined in Section 2.3 (for (a) and (c))—are depicted in light and dark grey. Note that from the graph analysis in Section 2.3 follows that only one vertex code exist in (a) and (c).

Inspired by applications in computer graphics (Cohen et al., 2003), we have introduced the formalism of Wang tiles as a suitable generalization of the SEPUC representation of heterogeneous materials (Novák et al., 2012). The formalism provides a compromise between the SEPUC approach and the use of a microstructure generator for each new realization. The concept of Wang tiles decomposes microstructure generation into an offline phase, which can possibly be computationally intensive, and a nearly instantaneous online phase. In the offline phase, information on a heterogeneous microstructure is compressed to a set of smaller domains—Wang tiles—with predefined compatibility constraints on the compressed microstructural geometry. In the online phase, microstructural realizations are assembled from these domains with a fast, linear algorithm that produces stochastic realizations with suppressed periodicity. Merits of the tiling concept for RVE size analyses were demonstrated in (Doškář and Novák, 2016; Doškář et al., 2018).

Optimization-based approaches developed initially for the SEPUC design can be extended to incorporate generalized periodic boundary conditions and used to generate the morphology of tiles (Novák et al., 2012, 2013). However, the extension amplifies the major weakness of optimization approaches—their computational cost—making them prohibitively expensive for complex three-dimensional models. As a remedy, we have proposed a method motivated by Cohen et al. (2003) that combines a sample-based approach with quantitative spatial statistics (Doškář et al., 2014). While this method is by orders of magnitude faster than the optimization approach, it struggles when handling complex, percolated microstructures such as foams, producing corrupted ligaments in sample overlaps (Doškář and Novák, 2016).

In this paper, we extend Sonon et al.’s method to Wang tiles in order to produce tile-based representations of microstructures intractable by the former methods, see Fig. 1. To this goal, in Section 2, we review the fundamentals of the Wang tile concept and discuss, in details, connectivity and mapping between vertex- and edge/face-based definitions of Wang tiles. Next, we outline Sonon et al.’s method (Sonon et al., 2012, 2015) and describe modifications necessary to accommodate generalized periodicity, see Section 3. Finally, equipped with the adapted procedure, we illustrate some outputs of the procedure in two and three dimensions. Besides, we compare three tile sets—two stochastic and a strictly aperiodic one —of the same cardinalities in terms of periodicity artefacts, complementing our previous study (Doškář et al., 2014) that dealt only with the distribution of tile types.

2 Wang tiles

Wang tiles constitute the building blocks of the abstract concept used in this work. Albeit the shape of the tile domains can be any parallelogram (or a parallelepiped in 3D), for simplicity we assume only square (or cubic) domains in the sequel. All tiles from the tile set have compatibility codes attributed to their edges (faces), illustrated with colours in Fig. 2. These codes play the role of constraints during an assembly process of tile instances in a tiling—a portion of a plane (space) without any voids or overlaps; only tiles with the same codes on the abutting edges (faces) can be placed side by side. In addition, tiles can be neither rotated nor reflected during the assembly. Even though the last two requirements can be eliminated by modifying the definition of compatibility codes (Robinson, 1971), we retain them for practical purposes as they preserve orientation of the microstructure compressed within the tile set. The particular version of an assembly algorithm depends on the type of the tile set as discussed next.

2.1 History and applications

Originally, the concept of Wang tiles was introduced in the first-order predicate calculus as a visual surrogate to a decision problem of statements (Wang, 1961, 1965). The initial conjecture that the whole plane can be covered only if a periodic tiling exists (Wang, 1961) was disproved shortly after (Berger, 1966), triggering a pursuit of the smallest tile set allowing for strictly aperiodic tiling of the plane, see the classical book (Grünbaum and Shephard, 2016) and (Jeandel and Rao, 2015) for historical overviews. The quest seems to be over; supported by an extensive computer-aided search, Jeandel and Rao (2015) announced the set of 11 tiles over 4 codes for each edge orientation, stating that no smaller set exists.

Similarly to applications in biology (Rothemund et al., 2004), physics (Leuzzi and Parisi, 2000), and computer graphics (Kopf et al., 2006; Cohen et al., 2003), we use the concept only in its geometrical interpretation as a suitable formalism describing mutual compatibility of small domains. Except for the comparison study in Section 4.1, we also limit ourself to the stochastic tile sets introduced by Cohen et al. (2003). Besides the fact that deterministic tilings of the aperiodic sets often exhibit locally ordered patterns, the stochastic sets offer higher flexibility in design, i.e. choosing the number of tiles and codes and their distribution within the set333The limitations of aperiodic sets is especially critical in 3D, because only one aperiodic set of Wang cubes has been published so far (Kari and Culik, 1995)..

The assembly algorithm for stochastic sets works sequentially: an initially empty grid is traversed in a scanline (Lagae and Dutré, 2006) way; at each grid node possible candidates compatible with previously placed tiles are identified in the set and one candidate tile is randomly chosen and placed. Thus, the only requirement for the set design is that there exist at least one tile (but optimally two or more to preserve stochastic nature of the assembly) for every possible combination of codes on the upper and left-hand edges. In principle, the random choice from the subset can be replaced with an informed selection preferring e.g. different phase volume fractions in different regions of a domain according to a pre-generated Gaussian random field. This modification allows for correlations in the microstructure at length scales larger than the tile size. For the sake of brevity and without loss of generality, however, we use only the standard stochastic algorithm in this work.

The tiling concept’s ability to generate naturally looking patterns from a limited amount of samples—which proved highly appealing in computer graphics (Cohen et al., 2003; Kopf et al., 2006)—stems from the reduced periodicity in tiling assemblies. This feature complies well with our ambition to replace a SEPUC-based description of random heterogeneous materials with its generalization that would allow for fast synthesis of stochastic microstructural samples or microstructure geometries for entire macroscopic domains.

The design of the compressed microstructural representation consists of two steps. First, the cardinality of the tile set and a particular distribution of edge (face) codes is chosen, which controls the frequency of tile occurrence in a tiling. Next, tile interiors are designed such that (i) the generated microstructure is continuous across the compatible edges and (ii) assembled tilings resemble the desired microstructure. Note that the latter requirement is not posed on individual tiles. On the contrary, variability in the compressed representation is the main merit of the tiling concept. Tile interiors together with the edge compatibility carry local microstructural characteristics, while the fluctuations over distances larger than the tile size are facilitated via the tile assembly algorithm.

2.2 Edge- vs. vertex-based tile definitions

Complementary to the edge-based specification introduced above, Wang tiles can be defined through vertex codes. Unlike the original setting, the vertex-based definition allows for a direct control of tile states across vertices, preventing from pronounced repetitiveness when a visually distinctive microstructural feature falls to the vertex region (Cohen et al., 2003). To avoid this “corner problem”, Cohen et al. (2003) proposed marking the tile vertices with an additional set of codes, essentially overlaying two tile set definitions. Subsequently, Lagae and Dutré (2006) retained only vertex codes, reporting superior spectral properties of assembled patterns compared to the edge-based sets of the same cardinality.

Albeit the two definitions yield the same cardinality of tiles () in the complete set, i.e. set that contains all combinations of  codes, they differ when it comes to minimal stochastic sets, i.e. sets that contain at least two tiles for each admissible code combination of already placed neighbouring tiles ( for edge-based sets vs. in the case of vertex-based ones). The difference is further pronounced in three dimensions: a full face-defined set contains cubes compared to cubes in a vertex-defined set. Theoretically, Wang cubes can be also defined with edges codes in three dimensions; however, the cardinality exceeds both the face- and vertex-based definitions without any known benefits.

Mapping vertex-defined tiles444The same procedure holds also for Wang cubes, where four vertices define one face code. to the original definition is straightforward: two vertex codes define one edge code. Hence all the above-mentioned design and assembly procedures directly apply to the vertex-based definition as well. In fact, the mapping provides an effortless way to produce a minimal stochastic edge-defined set with equal occurrence probability of each tile, the set characteristics sought for in Novák et al. (2013). On the other hand, because the mapping is injective only, it is not possible to map an edge-defined set to a vertex-based one in general.

2.3 Vertex analysis of a tile set

Figure 3: Illustration of the generalized periodicity effect on the particle copy process. Three shown steps illustrate how a particle initially placed at the top left corner of the first tile (grey particle on the left-hand side of the figure) is sequentially copied to other loci, following the edge codes it intersects. Dashed particle outlines loci where the particle will be copied in the next step.
(a) (b) (c)
Figure 4: Comparison of graphs pertinent to the tile sets from Fig. 2, identifying a potential vertex-based definition of the sets. Each graph node corresponds to a code on a half of a tile edge (either left/right for horizontal codes, or top/bottom for vertical codes), and graph arcs correspond to tile vertices. As there are multiple vertices—and hence parallel arcs—connecting the same nodes, the graphs should be drawn as a multigraph. However, multiple arcs are collapsed into one for the sake of brevity; only in (a) the multiplicity is indicated by tile indices attached to each arc. Independent sub-graphs indicating the vertex character of the set are shown in distinct light and dark grey, which correspond to the vertex codes in Fig. 2.

Here, we outline a tile set analysis capable of revealing the underlying vertex definition, if exists. The motivation behind the analysis is the pragmatic question arising in implementing the generalized periodicity: If a particle intersects a tile boundary, what other tiles should it be copied to? While the answer is straightforward for a particle intersecting a tile edge, vertex overlap is more involved. A particle overlapping a vertex is carried to other tiles by both the horizontal and vertical edges, and because the particle overlapped a vertex in the original tile, its images overlap vertices of the other tiles as well. Consequently, these images are propagated further by the newly affected edge codes, see Fig. 3. Whether the particle is eventually copied to all vertices or appears only in a selected subset depends on the allocation of codes to individual tiles.

Assume an undirected graph555Instead of the standard vertex-edge nomenclature for a graph, we use node and arc terms in order to avoid confusion with similar geometrical notions related to Wang tiles. where each node represents the code on a particular half of an edge (either top/bottom for vertical or left/right for horizontal codes, resulting in two occurrences of each edge code) and each arc corresponds to a tile vertex. The graph is by definition bipartite, because each vertex connects horizontal and vertical codes, and represents a multigraph as there are usually more vertices with the same adjacent edge codes in the tile set. To answer the aforementioned question, we identify connected components of the graph using the Depth-first algorithm (Sedgewick, 2002, Section 18.2). Each independent sub-graph then corresponds to a distinct vertex code and the arcs pertaining to the sub-graph determine the vertices to which a vertex-overlapping particle will be propagated. See Fig. 4 for a comparison of three graphs pertinent to the tile sets depicted in Fig. 2. The subgraphs (if present) are plotted in distinct grey colours, which correspond to the colours of vertex codes shown in Fig. 2.

In three dimensions, a set of Wang cubes can be analysed analogously with only a minor modification: a graph node represents one corner of a cube face, thus, each face code appears four times in the graph. In addition, cube edges with the same direction vector must be also classified, addressing the situation when a particle intersects cube edge. The classification follows the exact same procedure as vertex identification in two dimensions, neglecting the codes on faces perpendicular to the analysed edges.

3 Tile design using level-set functions

As mentioned in the Introduction, Random Sequential Adsorption (RSA)666Sometimes, the term Random Sequential Addition, e.g. (Talbot et al., 1991), is used interchangeably. is one of the most frequent algorithms for generating particle packings and microstructural geometries in general. It follows a simple idea of throwing particles of an arbitrary shape in a given domain and keeping those that do not overlap with previously placed ones. In its original setting, however, the algorithm poses a critical drawback for higher volume fractions: the success rate of accepting the particle position rapidly decays in later stages, because the majority of the randomly selected positions collide with the already placed particles. Optimally, the remedy would be to sample the new particle position directly from a domain that is certified to result in a non-overlapping composition.

To provide such a remedy, Sonon and coworkers Sonon et al. (2012); Sonon and Massart (2013); Sonon et al. (2015); Massart et al. (2018) adopted the implicit, level-set-based description of a microstructural geometry and demonstrated that it enables—in addition to the desired sampling from a valid domain loci—generating complex microstructural geometries (such as open and closed foams) in a unified framework, using suitable morphing operations. Moreover, both features are independent: if a particle packing is the desired output, the morphing operations can be neglected; and, on the other hand, any distribution of points or particles obtained by different packing methods can serve as the input for the morphing operations.

In the rest of this section, we recall individual steps of their framework and introduce necessary modifications facilitating the compatibility constraints arising in the concept of Wang tiles. For the sake of clarity, all procedures are presented in the two-dimensional setting, with comments on three-dimensional problems when necessary. Section 4 then presents both two- and three-dimensional results.

3.1 Generating particle packings

The original approach of Sonon et al. (2012) rests on describing a geometry of particle with a level-set function , which gives the signed distance of point from the nearest particle boundary (with negative values inside the particle), such that


with the signed distance function given by


Consequently, the zero iso-line of represents the particle boundary .

Assume for now that we have domain that contains set of already placed particles and their geometry is encoded in a single level-set field ,


see Fig. 5 for an illustration.

(a) (b)
Figure 5: Two consecutive steps of the level-set based RSA algorithm for Periodic Unit Cell: (a) level-set field with previously placed particles (values of outside the admissible sampling domain are displayed as translucent) at the th step; (b) updated state of and after placing a new particle at step .

Knowing the radius of the smallest circumscribed circle as the only characterization of the new particle, overlaps with the existing particles can be readily prevented by sampling its centre from the domain


After placing the new particle , the level-set field is updated with ,


This procedure repeats until either a pre-defined stopping criterion (e.g. a desired volume fraction) is met or .

To increase the volume fraction that can be achieved with this procedure, two additional restrictions can be posed on the admissible subdomain . First, a maximal distance from an existing particle boundaries can be added to Eq. (4), preventing too large gaps between particles. (If requested, a minimal distance between two particles can be further enforced as well.) Second, to prevent locally jammed states with large interparticle gaps, an additional field storing the shortest distance to the surface of the second nearest particle can be included in Eq. (4) and limited from above with , see (Sonon et al., 2012) for additional details. With all these constraints in action, Eq. (4) takes the form


Enforcing periodicity of the packing is also straightforward: upon a placement of a new particle, domain level-set field is updated also with the particle’s periodic images (eight in 2D, recall the grid in Fig. 6, and 26 in 3D).

3.1.1 Extension to Wang tiles

Figure 6: Neighbourhood grid identified for each tile in a tile set (here shown for the tile V6 from Fig. 2b) containing codes of tiles that can occur at given positions. Light and dark grey regions depict individual entities in the book-keeping structure preventing duplicated calculations. (During the update procedure described in the text, the blue particle can be accounted for either in tile 6 at the centre, or from tiles 3, 6, 11, or 14 at the first column, second-row position. With the book-keeping structure considered, only one update is performed.) The blue square illustrates the use of the acceleration via a precomputed patch; values of are computed only once and used for all particle occurrences during the update. If requested, values outside the patch are computed in a regular way.

Compared to generating PUCs, the additional compatibility constraints among individual Wang tiles necessitate several modifications in the original algorithm to generate continuous microstructural morphology of a tile set. As we use a level-set field defined on a regular grid for each tile in the tile set .

For each tile in a pre-processing stage, we identify potential neighbouring tiles and store their indices in a grid ( or in two or three dimensions, respectively), see Fig. 6, which will be used later in the updating phase.

For further description, we define a copy-inducer as a geometrical entity that is responsible for inducing particle copies to other tiles. Hence, a copy-inducer can be either a vertex, edge, face, or formally the tile itself and is related to a code identified in Section 2.3. The particular type of the copy-inducer stems from particle’s intersections with a tile boundary (if the particle does not intersect the boundary, tile is the copy-inducer, indicating that no copy is necessary; if it intersects one edge, the edge is the copy-inducer; etc.).

Modification 1: Update

The main structure of the modified algorithm follows the original one. After sampling a particle location from , which yields tile and centre where the new particle will be placed, particle’s intersections with the boundary of are determined, and—based on the intersections—a copy-inducer is identified. Images of the particle are subsequently copied to the relevant tiles, following the matching codes on the corresponding copy-inducers. Level-set field of the whole set is then updated tile-wise as follows.

For each tile, we sequentially loop over all tile positions in the neighbours grid shown in Fig. 6 (including the centre position as well) and check each of the potential tiles at that position for a new particles added in the current step. If any, the central tile’s level-set field is updated according to Eq. (5). Clearly, the new particle can appear several times in the same place. To avoid duplicated calculations, we set up a book-keeping structure that records which copy-inducer in the neighbouring grid was already used during the update, and we run only the unperformed updates, see Fig. 6.

Two implementation approaches are possible at this stage: (i) is computed anew for each tile and position in the book-keeping structure, or (ii) provided that is represented with discrete values on a regular grid, is precomputed for an auxiliary domain/patch aligned with the particle centre before any update and this field is then reused for all particle instances. Moreover, we combine both approaches with the pre-screening acceleration proposed in Sonon et al. (2012) that resorts to computing potentially demanding surface-distance functions only when the value of a circumscribed circle/sphere—which is fast to calculate—is less than the actual value for given .

While the pre-screening acceleration is significant in all cases, the benefits of a precomputed patch depends mainly on the patch size relative to the particle circumscribed radius, field resolution, and cardinality of the tile set. Especially in the later stages of the algorithm when only relatively small parts of are updated, pre-computing a patch for a particle that is not copied to other tiles might cause an unnecessary overhead. Therefore, we use the patch approach only for particles with either vertex, edge, or face copy-inducers.

For the three-dimensional set of 16 cubes used in Section 4.3 with ellipsoidal particles of circumscribed radius the patch pre-calculation delivered 10 % saving in computational time on average. Note also that the multi-query problem of finding the shortest distance to a particle boundary from a set of points is trivially parallelizable, changing the trade-off between patch pre-computations and direct calculations with available threads.

When all tiles are updated, the algorithm proceeds with identifying for the next particle.

Modification 2: Artificial level-set field

Definition of the admissible domain must be also modified in the case of Wang tiles; particle centres cannot be sampled directly from obtained from . The reason behind it is that the near-boundary parts of each tile must be informed about the interiors of related tiles to prevent insertions similar to the one depicted in Fig. 7a, where a copied instance of a newly placed particle collides with interior particle(s) of other tiles.

(b) (c)
Figure 7: (a) Illustrations of particle intersections that can happen due to copying particle images without considering artificially updated level-set field in the definition of . The problematic overlaps occur only at inducer places (highlighted in blue) with the same orientation as the particle copy-inducer (dashed outlined grey particle on the left-hand side). Figures (b) and (c) depict the original and its modified version that communicates necessary data across boundary regions to prevent intersections from (a).

The problem arises only when a particle image is added to a copy-inducing entity with the same orientation as the actual inducer (left-hand side edges in the illustration); particles copied to the opposite copy-inducers automatically meet the non-overlap requirement thanks to the first update modification described above.

To avoid this problem, we replace in Eq. (4) or Eq. (6) with a modified field which is equal to except for boundary regions, whose width is dictated by the radius of the circle/sphere circumscribed to the current particle. Each boundary region is then constructed as a point-wise minimum over all regions related to the same copy-inducers with the same code and the same orientation; see Fig. 7c. This construction propagates the close-boundary state across relevant tiles and thus it prevents the overlapping insertions similar to the one depicted in Fig. 7a. Note that a similar modification is not required for , because all information necessary to prevent the collisions is already contained in .

Modification 3: Breaking regular grid

The level-set field is typically implemented on a regular grid, the particle centres thus end up aligned with the grid. To break this artificial ordering, after sampling a new particle centre, we check whether the surrounding points belong to . The belonging points then define quadrants/octants in which the particle can be freely moved. We generate a random shift within these quadrant/octants and update the particle centre accordingly. However, this modification is possible only if the grid spacing is sufficiently small compared to the particle size; otherwise the implicitly assumed linear approximation of is inaccurate and might result in particle intersections at intermediate locations.

The above described modifications are sufficient to extend the original level-set-based packing algorithm of Sonon et al. to the concept of Wang tiles. Note that all control variables, e.g. particle radius or minimal and maximal distances and , can change during the algorithm operation. For instance, to achieve denser packings for multi-modal particle size distributions, it is preferential to start placing large particles first and sequentially proceed to the smaller ones.

Nonetheless, our extension inherits the weakness of the original algorithm when it comes to particles whose shapes are significantly different than their circumscribed circle/sphere, e.g. prolonged ellipsoids. In particular, the more the particle volume deviates from the volume of the circumscribed sphere, the less dense the final packing can get. While the geometry of the already placed particles is exactly described with (up to inaccuracies due to the grid spacing in

), the newly placed particle is represented only approximately during the centre selection, giving the pessimistic estimate of mutual intersections. The potential remedy is to translate the particle centre right after the sampling according to some heuristic rule aimed at producing denser packings; however, we leave this issue unaddressed in this paper.

3.2 Morphing operations

The acceleration of RSA by adopting the level-set approach is attractive by itself; however, the main appeal of the level-set framework of Sonon et al. is the elegance with which complex microstructures can be generated. As the simplest example, adding a constant to enables (i) fine-tuning of particle volume fractions, (ii) smoothing of sharp corners of polyhedral particles, and (iii) particle coating, e.g. for defining the Interfacial Transition Zone.

Combinations of and permit rendering interparticle bridges and controlled Voronoï-like tessellations; see (Sonon et al., 2012), where these morphing operations were used to mimic microstructural geometry of clay/sand mixed soils and irregular masonry. Considering also the shortest distance to the boundary of the third nearest particle , Sonon et al. (2015) demonstrated capability of their framework to produce highly adjustable models of foam-like microstructures, including features such as smooth transition from open to closed foams, concavity of foam ligaments, their coating and hollow interior, and variable thickness of foam ligaments and cell walls.

Similarly to the previous subsection, in the following, we present only the essentials of generating foam-like microstructures and introduce necessary modifications. The reader is referred to (Sonon et al., 2012, Section 4) and (Sonon et al., 2015, Sections 3 and 4) for details on the above-mentioned local adjustments, as they remained unaltered.

Assume we have at our disposal all three fields , , and , containing the shortest distance to the boundary of three nearest particles (either a posteriori computed for a given particle assembly or already tracked during the packing algorithm).

Similarly to the classical Voronoï tessellation, where a domain is partitioned by boundaries that have the same distance to the two closest seeds, thresholding the modified difference


yields a closed-foam-like geometry with indirectly controlling the thickness of foam cell walls. In the same spirit, the centrelines of the open-foam ligaments can be viewed as loci with the same distances to the boundary of three nearest particles; hence thresholding


produces open-foam ligaments, for which governs their thickness. Finally, combination of Eqs. (7) and (8)


allows for more realistic geometries with material concentration at wall intersections, see Fig. 8.

(a) (b)
Figure 8: Influence of morphing operations on the resulting (a) two- and (b) three-dimensional geometry. From the original particle distribution (depicted with light grey), the closed-foam-like geometry (shown in blue) is obtained by Eq. (7) and the open-like features (plotted in red) follow from Eq. (8).
(a) (b) (c)
Figure 9: Assembled level-set fields composed of tiles without (the top row) and with (the bottom row) a virtual boundary inset taken into account during particle placement. Note that even the nearest-neighbour distance (a) is not continuous across vertices when the virtual inset is not considered albeit the particle boundaries, i.e. zero-value contours outlined in grey, are. Severe discontinuities then appear for the second-nearest (b) and the third-nearest (c) neighbour distance without the inset boundaries.

Critical in extending the original framework to Wang tiles, with respect to the morphing procedures, is ensuring that all three fields are continuous across the relevant tile edges/faces. As illustrated in Fig. 9, copying particle images introduced in the previous subsection is insufficient in this regard. To guarantee the required continuity, we define a wider domain margin as a subdomain related to individual copy-inducers. Instead of computing intersections with a tile boundary, we compute them for a virtual inset boundary (inset by is usually enough) and copy the particles according to the virtual copy-inducers. Consequently, wider portions of tile domains are restrained near boundaries, resulting in the restored continuity. This modification reflects also in the construction of auxiliary field .

4 Results

4.1 Comparison of 16-tile sets

(a) (b) (c)
Figure 10: cross-sections along coordinate for tilings for (a) the edge-defined tile set C16, (b) the vertex-defined set V16, and (c) the Amman set A16, considered in Section 2.3. Thin grey lines depicts the average over 100 realizations of each tile set morphology, solid blue line shows the overall average (i.e. over realizations and morphologies) and green area captures the range between minimal and maximal values. Dashed grey line marks the asymptotic value , where the averaged volume fraction is considered.

First, we compare the periodicity reduction in three planar sets depicted in Fig. 2, which have the same cardinality but different origin. We supplement our earlier observations Doškář et al. (2014) regarding suppressed artificial periodicity compared to PUC-based reconstructions and better fitness of stochastic tile sets over aperiodic counterparts. Namely, we consider: (C16) an edge-based stochastic set over four colours,777This set corresponds to the set reported in Cohen et al. (2003). (V16) a vertex-based stochastic set over eight colours, and (A16) the aperiodic Ammann’s set Ammann et al. (1992) over 12 colours888For the aperiodic assembly we used the substitution rule from Stam (1997)..

Assuming ergodicity of a microstructure, we quantify the artificial periodicity in the tiling assemblies by means of the secondary peaks in the two-point probability function  Torquato (2002). states the chance of finding two points separated by in the same phase. Hence, attains its maximum at , where it equals the volume fraction of a chosen phase. For microstructures without any internal ordering, with because two sufficiently distant points are uncorrelated. On the other hand, if a microstructure is composed of a repeating (SE)PUC, exhibits secondary peaks of the same magnitude as at nodes of a regular grid whose spacing corresponds to the (SE)PUC size. Thus, in the intermediate case of Wang tiling, the magnitude of secondary peaks indicates the remaining artificial ordering in the reconstructed microstructure, see Novák et al. (2012); Doškář et al. (2014) for details. For computing

, we used the Fast Fourier Transform (FFT) instead of random sampling, as the bias in the statistics caused by the implicit periodicity introduced with FFT is negligible for reasonably large microstructural realizations 

Gajdošík et al. (2006).

For each of the three tile sets we ran the presented algorithm 25 times and thus generated 25 different tile set morphologies. Aiming at capturing the artificial periodicity related to the limited number of edge codes and tiles, we resorted to the simplest set-up with circular inclusion of a fixed radius in order to minimize the influence of particle shapes. In addition, we prevented the particles from overlapping tile vertices in order not to favour the vertex-based tile set a priori999From Fig. 4 follows that a particle will be copied to all tiles if it overlaps any vertex in the case of A16 and C16 tile sets.. Level-set fields of each unit-size tile () was discretized using a regular grid comprising points. Finally, we posed two constraints and on , recall Eq. (6).

Figure 11: Box-and-whisker plot of the normalized maximal secondary peaks

over 100 realizations and 25 tile set morphologies for each tile set. The central lines mark the median; edges of each box denote the first and third quartiles; and whisker lines contain all data not considered as outliers (marked as dark blue dots).

We assembled 100 microstructural realizations, i.e. tilings, containing tiles for each tile set and tile morphology and computed statistics for the inclusion phase. Cross-sections of along coordinate are plotted in Fig. 10 for the considered tile sets. As expected, the secondary peaks appear at loci with integer coordinates. Excluding the primary peak, we picked the second highest extreme normalised against the averaged asymptotic value for each realisation and plotted the obtained data with box-and-whisker diagrams in Fig. 11.

While locally the Ammann’s tile set A16 seems promising—note the suppressed secondary peaks at and in Fig. 10c—its deterministic aperiodic structure leads to pronounced secondary extremes in other regions, e.g. . This observation corroborates our previous conclusions, based solely on the distribution of individual tiles within tilings, that the strictly aperodic sets such as Ammann’s or Culik’s set considered in Doškář et al. (2014) lead to higher secondary peaks than stochastic sets Cohen et al. (2003) of the same cardinality, which feature more uniform secondary peaks. Figs. 10 and 11 clearly show that, out of the two stochastic sets, the vertex-based definition performs better as expected, even despite the prevented vertex overlaps. This superiority of vertex-based tile sets in two dimensions is further supported by additional results comparing vertex- and edge-based sets comprising 81 tiles, see Fig. 11.

4.2 2D example

Second, we present an example of a two-dimensional microstructure based on polygonal particles. Geometry of each particle was derived from an originally regular, randomly-oriented polygon with the number of vertices being sampled from a normal distribution with the mean value

and the standard deviation

, , and rounded towards the nearest integer. The angle between two rays connecting the particle centre and neighbouring vertices was perturbed with a value randomly chosen from . Finally, each vertex was placed on its corresponding ray at a distance sampled from (and capped by ) relative to the particle’s circumscribed radius.

Following the outcomes of Section 4.1, we picked the vertex-based set with 16 tiles depicted in Fig. 2b. Similarly to the previous set-up, the size of each tile was and the corresponding level-set fields were discretized using a regular grid comprising points. The initial value of the circumscribed radius was set to . After 40 algorithm steps, the radius was increased to , and eventually we reduced it to 0.06 after next 60 steps. Only the first two constraints on from Eq. (6) were posed. We fixed the minimal particle distance at ; the maximal distance was set to in the first 50 steps and decreased to later. The width of the virtual boundary inset was kept constant at .

(a) (b)
Figure 12: Results of the algorithm set-up described in Section 4.2: (a) A composed view of the tile set with highlighted edge and vertex codes. The particulate assembly is drawn in grey, the corresponding foam-like microstructure is shown in blue (similarly to Fig. 8); (b) An example of a tiling.

Fig. 12 depicts a composed view of the resulting particle assembly and the related closed-foam-like microstructural geometry obtained by considering and in Eq. (9).

4.3 3D example

Finally, we demonstrate a three-dimensional output of the algorithm for a set of 16 Wang cubes with two codes for each face orientation uniformly distributed in the set.

Again, we used centred tile domains of a unit size discretized with a grid of points. For the packing part of the algorithm, we considered ellipsoidal particles with a random orientation and a fixed circumscribed radius . The ratio between the middle and the major semi-axes lenghts was sampled from a uniform distribution ; the minor to the major semi-axes length ratio was randomly generated from . The admissible domain at each algorithm step was dictated both by and with constraints , , and in Eq. (6). Since we aimed at a foam-like microstructure, the width of a virtual inset boundaries was set to .

The final microstructure was obtained by performing the morphing operations according to Eq. (9) with and . For the sake of brevity, we don’t show individual Wang cubes, and we plot only a microstructural sample comprising tiles in Fig. 13.

Figure 13: A microstructural realization assembled as a tiling from the Wang tiles generated in Section 4.3.

5 Summary

In this work, we have extended the level-set based framework of Sonon et al. (2012) to the microstructural models based on the formalism of Wang tiles. The extension enables generating compressed representations of complex microstructural geometries such as open and closed foams or cells, which have been nearly impossible Doškář et al. (2014) or very expensive Novák et al. (2012) to generate for Wang tiles so far.

Advancing from the standard Periodic Unit Cell to the generalized periodicity of Wang tiles necessitated several modifications of the original algorithm that addressed difficulties originally not encountered in the case of PUC. Driven by a geometrically motivated question of where a vertex-overlapping particle should be copied, we have come up with a simple procedure based on a graph analysis capable of revealing an underlying vertex definition of a tile set if present. We have also demonstrated that a straightforward copying boundary-intersecting particles according to the edge and vertex codes is insufficient to prevent from spurious particle overlaps, because the information about tile interior is not automatically communicated across individual tiles. As a remedy, we have introduced an artificially updated level-set field that facilitates the necessary communication. The width of the updated region in this field is dictated by the required final geometry. If only a particle distribution is to be generated, the width equal to the particle radius suffices. On the other hand, if a foam-like microstructures are desired, the virtual boundary inset by is necessary to preserve continuity of the level-set fields and across tile edges, resulting in the updated region of width .

Additionally, we have devised two minor modifications: one capable of breaking a regular underlying grid of possible particle centres by introducing their random shifts within an admissible region, the other accelerating the updates by pre-computing on a patch. Importance of the latter modification grows with higher resolution and larger cardinality of a tile set.

Having a universal framework for generating Wang tile morphologies at our disposal, we have supplemented our previous comparison of strictly aperiodic and stochastic tile sets Doškář et al. (2014) and confirmed the superiority of vertex-defined tile sets in suppressing artificial periodicity in assembled microstructural samples. Given the same cardinality of the edge- and vertex-based tile sets in two dimensions, we recommend using the latter. On the other hand, the same comparison in three dimensions is more subtle and the choice should be always a compromise between available computational resources and required periodicity reduction, following e.g. an approximate formula in Novák et al. (2012).

Admittedly, the presented modification inherits certain limitations of the original framework, namely it provides only an indirect control of the resulting microstructural geometry and no spatial statistics are involved in the particle placing process. Albeit addressing these limitations remained out of the scope of this work, ideas conceptually similar to Yang et al. (2018), i.e. optimizing the particle positions a posteriori to minimize the discrepancy between the target and computed spatial statistics, could be potentially applied to steer the generated microstructure to the desired statistics.

The generalized periodicity also brings slightly increased computational cost compared to PUC generation. Recall though that the tile-based compression is intended mainly for generating numerous stochastic realizations of arbitrary size. Thus more time can be spent on the initial microstructure preparation as the increased cost will be amortized later by the repeated use of the same compression.

Related to the last remark, our current focus is on developing a robust finite element discretization tool that would enable meshing both outputs of the framework, i.e. analytical geometries of particulate microstructures and complex geometries implicitly defined via level-set fields and processed with the marching-cube algorithm. In both cases, special attention is paid to ensuring the topological and geometrical compatibility of the resulting finite-element meshes across the corresponding faces/edges.


M. Doškář, J. Novák, and D. Rypl acknowledge the endowment of the Ministry of Industry and Trade of the Czech Republic under the project No. FV10202. This research has been performed in the Center of Advanced Applied Sciences (CAAS), financially supported by the European Regional Development Fund (project No. CZ.02.1.01/0.0/0.0/16_019/0000778). M. Doškář was also partially supported by the Grant Agency of the Czech Technical University in Prague through a student project No. SGS18/036/OHK1/1T/11, and J. Zeman received partial support from the Czech Science Foundation project No. 17-04301S.


  • Nguyen et al. (2011) V. P. Nguyen, O. Lloberas-Valls, M. Stroeven, L. J. Sluys, Homogenization-based multiscale crack modelling: From micro-diffusive damage to macro-cracks, Computer Methods in Applied Mechanics and Engineering 200 (9-12) (2011) 1220–1236, ISSN 00457825, doi:10.1016/j.cma.2010.10.013.
  • Barchiesi et al. (2019) E. Barchiesi, M. Spagnuolo, L. Placidi, Mechanical metamaterials: a state of the art, Mathematics and Mechanics of Solids 24 (1) (2019) 212–234, ISSN 1081-2865, 1741-3028, doi:10.1177/1081286517735695.
  • Matouš et al. (2017) K. Matouš, M. G. Geers, V. G. Kouznetsova, A. Gillman, A review of predictive nonlinear theories for multiscale modeling of heterogeneous materials, Journal of Computational Physics 330 (2017) 192–220, ISSN 00219991, doi:10.1016/j.jcp.2016.10.070.
  • Torquato (2002) S. Torquato, Random heterogeneous materials: microstructure and macroscopic properties, no. 16 in Interdisciplinary applied mathematics, Springer, New York, ISBN 0-387-95167-9, 2002.
  • Nguyen et al. (2017) T. Nguyen, J. Yvonnet, M. Bornert, C. Chateau, F. Bilteryst, E. Steib, Large-scale simulations of quasi-brittle microcracking in realistic highly heterogeneous microstructures obtained from micro CT imaging, Extreme Mechanics Letters 17 (2017) 50–55, ISSN 23524316, doi:10.1016/j.eml.2017.09.013.
  • Lloberas-Valls et al. (2011) O. Lloberas-Valls, D. Rixen, A. Simone, L. Sluys, Domain decomposition techniques for the efficient modeling of brittle heterogeneous materials, Computer Methods in Applied Mechanics and Engineering 200 (13-16) (2011) 1577–1590, ISSN 00457825, doi:10.1016/j.cma.2011.01.008.
  • Akbari Rahimabadi et al. (2015) A. Akbari Rahimabadi, P. Kerfriden, S. Bordas, Scale selection in nonlinear fracture mechanics of heterogeneous materials, Philosophical Magazine 95 (28-30) (2015) 3328–3347, ISSN 1478-6435, 1478-6443, doi:10.1080/14786435.2015.1061716.
  • Lloberas-Valls et al. (2012) O. Lloberas-Valls, D. Rixen, A. Simone, L. Sluys, Multiscale domain decomposition analysis of quasi-brittle heterogeneous materials, International Journal for Numerical Methods in Engineering 89 (11) (2012) 1337–1366, ISSN 00295981, doi:10.1002/nme.3286.
  • Lee et al. (2009) H. Lee, M. Brandyberry, A. Tudor, K. Matouš, Three-dimensional reconstruction of statistically optimal unit cells of polydisperse particulate composites from microtomography, Physical Review E 80 (6) (2009) 061301, doi:10.1103/PhysRevE.80.061301.
  • Yang et al. (2018) M. Yang, A. Nagarajan, B. Liang, S. Soghrati, New algorithms for virtual reconstruction of heterogeneous microstructures, Computer Methods in Applied Mechanics and Engineering 338 (2018) 275–298, ISSN 00457825, doi:10.1016/j.cma.2018.04.030.
  • Balzani et al. (2014) D. Balzani, L. Scheunemann, D. Brands, J. Schröder, Construction of two- and three-dimensional statistically similar RVEs for coupled micro-macro simulations, Computational Mechanics 54 (5) (2014) 1269–1284, ISSN 0178-7675, 1432-0924, doi:10.1007/s00466-014-1057-6.
  • Zeman and Šejnoha (2007) J. Zeman, M. Šejnoha, From random microstructures to representative volume elements, Modelling and Simulation in Materials Science and Engineering 15 (4) (2007) S325–S335, ISSN 0965-0393, 1361-651X, doi:10.1088/0965-0393/15/4/S01.
  • Scheunemann et al. (2015) L. Scheunemann, D. Balzani, D. Brands, J. Schröder, Design of 3D statistically similar Representative Volume Elements based on Minkowski functionals, Mechanics of Materials 90 (2015) 185–201, ISSN 01676636, doi:10.1016/j.mechmat.2015.03.005.
  • Jiao et al. (2007) Y. Jiao, F. Stillinger, S. Torquato, Modeling heterogeneous materials via two-point correlation functions: Basic principles, Physical Review E 76 (3) (2007) 031110, ISSN 1539-3755, 1550-2376, doi:10.1103/PhysRevE.76.031110.
  • Zeman and Šejnoha (2001) J. Zeman, M. Šejnoha, Numerical evaluation of effective elastic properties of graphite fiber tow impregnated by polymer matrix, Journal of the Mechanics and Physics of Solids 49 (1) (2001) 69–90, doi:10.1016/S0022-5096(00)00027-2.
  • Rozman and Utz (2001) M. G. Rozman, M. Utz, Efficient reconstruction of multiphase morphologies from correlation functions, Physical Review E 63 (6) (2001) 066701, doi:10.1103/PhysRevE.63.066701.
  • Jiao et al. (2009) Y. Jiao, F. H. Stillinger, S. Torquato, A superior descriptor of random textures and its predictive capacity, Proceedings of the National Academy of Sciences 106 (42) (2009) 17634–17639, ISSN 0027-8424, 1091-6490, doi:10.1073/pnas.0905919106.
  • Lu and Torquato (1992) B. Lu, S. Torquato, Lineal-path function for random heterogeneous materials, Physical Review A 45 (2) (1992) 922–929, ISSN 1050-2947, 1094-1622, doi:10.1103/PhysRevA.45.922.
  • Kumar et al. (2006) H. Kumar, C. Briant, W. Curtin, Using microstructure reconstruction to model mechanical behavior in complex microstructures, Mechanics of Materials 38 (8-10) (2006) 818–832, ISSN 01676636, doi:10.1016/j.mechmat.2005.06.030.
  • Havelka et al. (2016) J. Havelka, A. Kučerová, J. Sýkora, Compression and reconstruction of random microstructures using accelerated lineal path function, Computational Materials Science 122 (2016) 102–117, ISSN 09270256, doi:10.1016/j.commatsci.2016.04.044.
  • Povirk (1995) G. L. Povirk, Incorporation of microstructural information into models of two-phase materials, Acta Metallurgica et Materialia 43 (8) (1995) 3199 – 3206, ISSN 0956-7151, doi:http://dx.doi.org/10.1016/0956-7151(94)00487-3.
  • Yeong and Torquato (1998) C. Yeong, S. Torquato, Reconstructing random media, Physical Review E 57 (1) (1998) 495–506, ISSN 1063-651X, 1095-3787, doi:10.1103/PhysRevE.57.495.
  • Basanta et al. (2005)

    D. Basanta, M. A. Miodownik, E. A. Holm, P. J. Bentley, Using genetic algorithms to evolve three-dimensional microstructures from two-dimensional micrographs, Metallurgical and Materials Transactions A 36 (7) (2005) 1643–1652, ISSN 1073-5623, 1543-1940,

  • Kumar et al. (2008) N. C. Kumar, K. Matouš, P. H. Geubelle, Reconstruction of periodic unit cells of multimodal random particulate composites using genetic algorithms, Computational Materials Science 42 (2) (2008) 352–367, ISSN 09270256, doi:10.1016/j.commatsci.2007.07.043.
  • Fullwood et al. (2008a) D. Fullwood, S. Kalidindi, S. Niezgoda, A. Fast, N. Hampson, Gradient-based microstructure reconstructions from distributions using fast Fourier transforms, Materials Science and Engineering: A 494 (1-2) (2008a) 68–72, ISSN 09215093, doi:10.1016/j.msea.2007.10.087.
  • Fullwood et al. (2008b) D. T. Fullwood, S. R. Niezgoda, S. R. Kalidindi, Microstructure reconstructions from 2-point statistics using phase-recovery algorithms, Acta Materialia 56 (5) (2008b) 942–948, ISSN 13596454, doi:10.1016/j.actamat.2007.10.044.
  • Liu and Shapiro (2015) X. Liu, V. Shapiro, Random heterogeneous materials via texture synthesis, Computational Materials Science 99 (2015) 177–189, ISSN 09270256, doi:10.1016/j.commatsci.2014.12.017.
  • Tahmasebi and Sahimi (2013) P. Tahmasebi, M. Sahimi, Cross-Correlation Function for Accurate Reconstruction of Heterogeneous Media, Physical Review Letters 110 (7) (2013) 078002, ISSN 0031-9007, 1079-7114, doi:10.1103/PhysRevLett.110.078002.
  • Bostanabad et al. (2016) R. Bostanabad, A. T. Bui, W. Xie, D. W. Apley, W. Chen, Stochastic microstructure characterization and reconstruction via supervised learning, Acta Materialia 103 (2016) 89–102, ISSN 13596454, doi:10.1016/j.actamat.2015.09.044.
  • Biswal et al. (1999) B. Biswal, C. Manwart, R. Hilfer, S. Bakke, P. Øren, Quantitative analysis of experimental and synthetic microstructures for sedimentary rock, Physica A: Statistical Mechanics and its Applications 273 (3-4) (1999) 452–475, ISSN 03784371, doi:10.1016/S0378-4371(99)00248-4.
  • Saito and Enomoto (1992) Y. Saito, M. Enomoto, Monte Carlo Simulation of Grain Growth, ISIJ International 32 (3) (1992) 267–274, ISSN 0915-1559, doi:10.2355/isijinternational.32.267.
  • Krill III and Chen (2002) C. Krill III, L.-Q. Chen, Computer simulation of 3-D grain growth using a phase-field model, Acta Materialia 50 (12) (2002) 3059–3075, ISSN 13596454, doi:10.1016/S1359-6454(02)00084-8.
  • Falco et al. (2017) S. Falco, J. Jiang, F. De Cola, N. Petrinic, Generation of 3D polycrystalline microstructures with a conditioned Laguerre-Voronoi tessellation technique, Computational Materials Science 136 (2017) 20–28, ISSN 09270256, doi:10.1016/j.commatsci.2017.04.018.
  • Alsayednoor and Harrison (2016) J. Alsayednoor, P. Harrison, Evaluating the performance of microstructure generation algorithms for 2-d foam-like representative volume elements, Mechanics of Materials 98 (2016) 44–58, ISSN 01676636, doi:10.1016/j.mechmat.2016.04.001.
  • Cooper (1988) D. W. Cooper, Random-sequential-packing simulations in three dimensions for spheres, Physical Review A 38 (1) (1988) 522–524, ISSN 0556-2791, doi:10.1103/PhysRevA.38.522.
  • Segurado and Llorca (2002) J. Segurado, J. Llorca, A numerical approximation to the elastic properties of sphere-reinforced composites, Journal of the Mechanics and Physics of Solids 50 (10) (2002) 2107–2121, ISSN 00225096, doi:10.1016/S0022-5096(02)00021-2.
  • Lubachevsky and Stillinger (1990) B. D. Lubachevsky, F. H. Stillinger, Geometric properties of random disk packings, Journal of Statistical Physics 60 (5-6) (1990) 561–583, ISSN 0022-4715, 1572-9613, doi:10.1007/BF01025983.
  • Donev et al. (2005) A. Donev, S. Torquato, F. H. Stillinger, Neighbor list collision-driven molecular dynamics simulation for nonspherical hard particles. I. Algorithmic details, Journal of Computational Physics 202 (2) (2005) 737–764, ISSN 00219991, doi:10.1016/j.jcp.2004.08.014.
  • Ghossein and Lévesque (2013) E. Ghossein, M. Lévesque, Random generation of periodic hard ellipsoids based on molecular dynamics: A computationally-efficient algorithm, Journal of Computational Physics 253 (2013) 471–490, ISSN 00219991, doi:10.1016/j.jcp.2013.07.004.
  • De Giorgi et al. (2010) M. De Giorgi, A. Carofalo, V. Dattoma, R. Nobile, F. Palano, Aluminium foams structural modelling, Computers & Structures 88 (1-2) (2010) 25–35, ISSN 00457949, doi:10.1016/j.compstruc.2009.06.005.
  • Chen et al. (2015) Y. Chen, R. Das, M. Battley, Effects of cell size and cell wall thickness variations on the stiffness of closed-cell foams, International Journal of Solids and Structures 52 (2015) 150–164, ISSN 00207683, doi:10.1016/j.ijsolstr.2014.09.022.
  • Mebatsion et al. (2006) H. Mebatsion, P. Verboven, B. Verlinden, Q. Ho, T. Nguyen, B. Nicolaï, Microscale modelling of fruit tissue using Voronoi tessellations, Computers and Electronics in Agriculture 52 (1-2) (2006) 36–48, ISSN 01681699, doi:10.1016/j.compag.2006.01.002.
  • Chakraborty et al. (2013) A. Chakraborty, M. M. Perales, G. V. Reddy, A. K. Roy-Chowdhury, Adaptive Geometric Tessellation for 3D Reconstruction of Anisotropically Developing Cells in Multilayer Tissues from Sparse Volumetric Microscopy Images, PLoS ONE 8 (8) (2013) e67202, ISSN 1932-6203, doi:10.1371/journal.pone.0067202.
  • Doškář and Novák (2016) M. Doškář, J. Novák, A jigsaw puzzle framework for homogenization of high porosity foams, Computers & Structures 166 (2016) 33–41, ISSN 00457949, doi:10.1016/j.compstruc.2016.01.003.
  • Simone and Gibson (1998) A. Simone, L. Gibson, The effects of cell face curvature and corrugations on the stiffness and strength of metallic foams, Acta Materialia 46 (11) (1998) 3929–3935, ISSN 13596454, doi:10.1016/S1359-6454(98)00072-X.
  • Sonon et al. (2012) B. Sonon, B. Franois, T. Massart, A unified level set based methodology for fast generation of complex microstructural multi-phase RVEs, Computer Methods in Applied Mechanics and Engineering 223-224 (2012) 103–122, ISSN 00457825, doi:10.1016/j.cma.2012.02.018.
  • Sonon and Massart (2013) B. Sonon, T. Massart, A Level-Set Based Representative Volume Element Generator and XFEM Simulations for Textile and 3D-Reinforced Composites, Materials 6 (12) (2013) 5568–5592, ISSN 1996-1944, doi:10.3390/ma6125568.
  • Sonon et al. (2015) B. Sonon, B. François, T. J. Massart, An advanced approach for the generation of complex cellular material representative volume elements using distance fields and level sets, Computational Mechanics 56 (2) (2015) 221–242, ISSN 0178-7675, 1432-0924, doi:10.1007/s00466-015-1168-8.
  • Massart et al. (2018) T. J. Massart, B. Sonon, K. Ehab Moustafa Kamel, L. H. Poh, G. Sun, Level set-based generation of representative volume elements for the damage analysis of irregular masonry, Meccanica 53 (7) (2018) 1737–1755, ISSN 0025-6455, 1572-9648, doi:10.1007/s11012-017-0695-0.
  • Kanit et al. (2003) T. Kanit, S. Forest, I. Galliet, V. Mounoury, D. Jeulin, Determination of the size of the representative volume element for random composites: statistical and numerical approach, International Journal of Solids and Structures 40 (13-14) (2003) 3647–3679, ISSN 00207683, doi:10.1016/S0020-7683(03)00143-4.
  • Gitman et al. (2007) I. Gitman, H. Askes, L. Sluys, Representative volume: Existence and size determination, Engineering Fracture Mechanics 74 (16) (2007) 2518–2534, ISSN 00137944, doi:10.1016/j.engfracmech.2006.12.021.
  • Dirrenberger et al. (2014) J. Dirrenberger, S. Forest, D. Jeulin, Towards gigantic RVE sizes for 3D stochastic fibrous networks, International Journal of Solids and Structures 51 (2) (2014) 359–376, ISSN 00207683, doi:10.1016/j.ijsolstr.2013.10.011.
  • Cohen et al. (2003) M. F. Cohen, J. Shade, S. Hiller, O. Deussen, Wang Tiles for image and texture generation, ACM Transactions on Graphics 22 (3) (2003) 287–294, ISSN 07300301, doi:10.1145/882262.882265.
  • Novák et al. (2012) J. Novák, A. Kučerová, J. Zeman, Compressing random microstructures via stochastic Wang tilings, Physical Review E 86 (4) (2012) 040104, doi:10.1103/PhysRevE.86.040104.
  • Doškář et al. (2018) M. Doškář, J. Zeman, D. Jarušková, J. Novák, Wang tiling aided statistical determination of the Representative Volume Element size of random heterogeneous materials, European Journal of Mechanics - A/Solids 70 (2018) 280–295, ISSN 09977538, doi:10.1016/j.euromechsol.2017.12.002.
  • Novák et al. (2013) J. Novák, A. Kučerová, J. Zeman, Microstructural enrichment functions based on stochastic Wang tilings, Modelling and Simulation in Materials Science and Engineering 21 (2) (2013) 025014, ISSN 0965-0393, 1361-651X, doi:10.1088/0965-0393/21/2/025014.
  • Doškář et al. (2014) M. Doškář, J. Novák, J. Zeman, Aperiodic compression and reconstruction of real-world material systems based on Wang tiles, Physical Review E 90 (6) (2014) 062118, ISSN 1539-3755, 1550-2376, doi:10.1103/PhysRevE.90.062118.
  • Robinson (1971) R. M. Robinson, Undecidability and nonperiodicity for tilings of the plane, Inventiones Mathematicae 12 (3) (1971) 177–209, ISSN 0020-9910, 1432-1297, doi:10.1007/BF01418780.
  • Wang (1961)

    H. Wang, Proving Theorems by Pattern Recognition - II, Bell System Technical Journal 40 (1) (1961) 1–41, ISSN 00058580,

  • Wang (1965) H. Wang, Games, Logic and Computers, Scientific American 213 (5) (1965) 98–106.
  • Berger (1966) R. Berger, The Undecidability of the Domino Problem, American Mathematical Society, Providence, R.I., ISBN 978-1-4704-0013-2 1-4704-0013-8, 1966.
  • Grünbaum and Shephard (2016) B. Grünbaum, G. C. Shephard, Tilings and Patterns: Second Edition, Dover Publications, Inc, Mineola, New York, ISBN 978-0-486-46981-2, 2016.
  • Jeandel and Rao (2015) E. Jeandel, M. Rao, An aperiodic set of 11 Wang tiles, arXiv:1506.06492 .
  • Rothemund et al. (2004) P. W. K. Rothemund, N. Papadakis, E. Winfree, Algorithmic Self-Assembly of DNA Sierpinski Triangles, PLoS Biology 2 (12) (2004) e424, ISSN 1545-7885, doi:10.1371/journal.pbio.0020424.
  • Leuzzi and Parisi (2000) L. Leuzzi, G. Parisi, Thermodynamics of a tiling model, Journal of Physics A: Mathematical and General 33 (23) (2000) 4215–4225, ISSN 0305-4470, 1361-6447, doi:10.1088/0305-4470/33/23/301.
  • Kopf et al. (2006) J. Kopf, D. Cohen-Or, O. Deussen, D. Lischinski, Recursive Wang tiles for real-time blue noise, ACM Transactions on Graphics 25 (3) (2006) 509, ISSN 07300301, doi:10.1145/1141911.1141916.
  • Kari and Culik (1995) J. Kari, K. Culik, An Aperiodic Set of Wang Cubes, JUCS - Journal of Universal Computer Science 1 (10) (1995) 675–686, ISSN 0948695X, doi:10.3217/jucs-001-10-0675.
  • Lagae and Dutré (2006) A. Lagae, P. Dutré, An alternative for Wang tiles: colored edges versus colored corners, ACM Transactions on Graphics 25 (4) (2006) 1442–1459, ISSN 07300301, doi:10.1145/1183287.1183296.
  • Sedgewick (2002) R. Sedgewick, Algorithms in C++, Part 5: Graph Algorithms, Addison-Wesley, Boston, 3rd edn., ISBN 978-0-201-36118-6, 2002.
  • Talbot et al. (1991) J. Talbot, P. Schaaf, G. Tarjus, Random sequential addition of hard spheres, Molecular Physics 72 (6) (1991) 1397–1406, ISSN 0026-8976, 1362-3028, doi:10.1080/00268979100100981.
  • Ammann et al. (1992) R. Ammann, B. Grünbaum, G. C. Shephard, Aperiodic tiles, Discrete & Computational Geometry 8 (1) (1992) 1–25, ISSN 0179-5376, 1432-0444, doi:10.1007/BF02293033.
  • Stam (1997) J. Stam, Aperiodic Texture Mapping, Tech. Rep. R046, url: www.autodeskresearch.com/publications/aperiodictexture, 1997.
  • Gajdošík et al. (2006) J. Gajdošík, J. Zeman, M. Šejnoha, Qualitative analysis of fiber composite microstructure: Influence of boundary conditions, Probabilistic Engineering Mechanics 21 (4) (2006) 317–329, ISSN 02668920, doi:10.1016/j.probengmech.2005.11.006.