# Computer-Aided Design of Lacunar Fractal Polygons

Boris Bordeaux, Christian Gentil, Iván Rasskin, Lionel Garnier

► To cite this version:

Boris Bordeaux, Christian Gentil, Iván Rasskin, Lionel Garnier. Computer-Aided Design of Lacunar Fractal Polygons. 2025. hal-05407434

HAL Id: hal-05407434

https://hal.science/hal-05407434v1

Preprint submitted on 9 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Computer-Aided Design of Lacunar Fractal Polygons

Boris Bordeaux¹, Christian Gentil¹, Iván Rasskin², and Lionel Garnier¹

¹Université Bourgogne Europe, LIB UR 7534, 21000, Dijon, France

²LIS, Aix-Marseille Université, CNRS, 13288, Marseille, France

## Abstract

Imagining and designing complicated topologies, such as multi-scale lacunar objects, is challenging. Generally, designers generate these structures from ad hoc algorithms with limited options for parameterization. Some models, like Boundary Controlled Iterated Function Systems (BC-IFS), provide a complete formalism for describing and controlling both their topology and geometry. However, encoding the topology involves abstract knowledge and can be tedious to do manually. On the other hand, polyhedral crystallographic circle packings provide a complete construction of fractal shapes arising from polyhedra. The geometry of such packing relies on circles, but the simple input polyhedron entirely determines their topology. This article shows how the BC-IFS model can encode the topology of fractals obtained from polyhedral crystallographic circle packings. Our method automatically deduces the topological BC-IFS constraints directly from the polyhedron's structure. As a result, polyhedra serve as an intuitive topological editor for creating multi-scale lacunar polygons. Finally, BC-IFS allows the end-user to modulate the global aspect using control points and refine local details with subdivision points.

<figure><img src="image_2.png" /></figure>

## 1 Introduction

Designing complex shapes such as lattices and scaffolds to produce lacunar structures involves defining automatic methods and the development of suitable algorithms. These methods have been notably implemented for building lighter structures thanks to additive manufacturing [ZL24]. Sometimes, these lattice structures (macro, micro, or multi-scale) are the result of a topological optimization algorithm [WSG21], which highlights their interest in applications such as [AALS17] where self-similar structures of tree shapes appear. A possible approach involves assembling geometric patterns to produce more complex structures. One of the most common is to tile the space using an elementary volume (cube, tetrahedron) and place a 3D periodic pattern revealing voids inside each elementary volume. Youn Hong [HE21] starts from a volume representation (V-Rep), a trivariate tensor product splines, which he partitions regularly from a 3D grid defined in the parameter space. Then, each voxel of this grid is "replaced" by a geometric pattern defined in the parameter space. For trimmed volume assemblies, matching work between patterns must be carried out. Antolin [ABC+19], following this
same approach by V-Rep, proposes to adapt the structures according to constraints by either sizing
the geometric patterns or modulating the grid's resolution, impacting the voxel size. In the latter case,
it must process the connections between different resolution levels. Garnier [GSR22] constructs micro-
lattice structures using an appropriate anisotropic reaction-diffusion system. A tensor field resulting
from a stress analysis guides the diffusion system. Sitharam [SYNP19] proposes to conceive a multi-
scale structure based on a Kagome graph. This kind of graph has good physical properties, hence its
structures take advantage of it. However, the structures produced are only triangles or tetrahedra.
Applying the same approach again on each triangle or tetrahedron of the resulting structure induces
multi-scale structures.

While effective in generating lattice and multi-scale structures, these methods present some limita-
tions. A fixed cell topology facilitates junctions between elements but significantly restricts structural
variability. As a result, the generated designs may lack adaptability to specific constraints or functional
requirements. Additionally, handling multi-scale transitions remains a challenge. Refining the struc-
ture at different levels introduces abrupt scale changes, making handling junctions between resolution
levels difficult, especially when ensuring mechanical continuity or aesthetic coherence.

This paper aims to automatically generate various multi-scale lacunar cells and consistently assem-
ble them to create a structure with different local topologies. These multi-scale geometric patterns
are often present in nature and have a universally recognized aesthetic. Fractal geometry is a way
to describe and produce such shapes. The geometry and topology of these structures are complex,
with potentially infinite details and ramifications. There is little work to represent such structures.
Shapes are generally constructed using more or less specific algorithms, with limited functionality
and control [ZLL<sup>+</sup>23]. The consistency of the objects is algorithmically verified during the iterative
construction process. One of the main difficulties is controlling the topology of these structures and
finding a way to describe them.

Formalisms such as Hutchinson's Iterated Function System allow us to represent these structures. An extension of IFS, the Controled-IFS model, adds rules to control the iterative process. Rules define which transformations of the IFS to apply depending on a control graph. There are various ways to implement it. One can use formal languages by using Language Restricted IFS [TT93] or L-Systems [Pru86, PHHM96]. Barnsley introduced Recurrent IFS in [BEH89]. He uses an automaton as a control graph.

The generalization of recurrent IFS to Boundary Controlled Iterated Function System (BC-IFS) [SGGM15, GGS21] allows encoding of the topology of fractal structures by generalizing the notion of cellular decomposition to non-standard topologies (edges can be Cantor sets, faces fractal faces, etc.). Although BC-IFS describes the topology of fractal structures rigorously, this description can be tedious and complex for a non-expert due to the many topological constraints that need to be explicit. In addition, imagining self-similar multi-scale structures is not an easy exercise. To facilitate this process, Bordeaux [BG24] proposes to define lacunar polygons by choosing the type of edges and automatically constructing the lacunar subdivision of the polygon. This approach facilitates the design and assembly of structures by imposing the edge type.

On the other hand, crystallographic packings, introduced by Kontorovich and Nakamura in [KN19],
are a class of highly symmetric infinite $N$-sphere packings that densely fill the space $\mathbb{R}^N$, exhibiting a
rich algebraic structure and diversity [CRCS20]. Furthermore, their limit set forms a lacunar fractal
structure closely related to previous models constructed using BC-IFS [BG24, GGS21]. In particular,
we focus on a two-dimensional subclass known as polyhedral crystallographic packings [RAR25], which
depend on a three-dimensional convex polyhedron that dictates the topology of the limit set. However,
the construction process is geometrically constrained, relying exclusively on circles.

This article aims to use the combinatorial structure of polyhedral crystallographic packings to gen-
erate lacunar fractal polygons and formalize their topology using the BC-IFS model. The polyhedron
fully defines all incidence and adjacency constraints required to encode the fractal's topology. Once
this topological structure is established, the geometric embedding can be selected, as illustrated in this
study.

With this approach, we use polyhedra as topological editors for lacunar fractal polygons. Given an
initial polygon, the end-user selects an associated polyhedron that defines the topology of the fractal
structure. Additionally, we examine the necessary conditions for assembling distinct lacunar fractal
polygons while preserving the continuity of their boundaries. By incorporating the geometric embed-
ding provided by the BC-IFS model, end-users can construct three-dimensional patches of lacunar
fractal polygons that can be locally and globally deformed.

We introduce the BC-IFS model in Section 2.1. Section 2.2 describes the polyhedral crystallographic circle packing method. In Section 3, we detail how to automatically encode the topology of a polyhedral crystallographic circle packing with the BC-IFS model based on the polyhedron topology. Section 4 demonstrates a possible integration of our method into CAD software to model fractal structures. We also present a design rule that ensures topological consistency when assembling fractal polygons with different topologies. Finally, Section 5 discusses the limitations of our approach and outlines directions for future work.

# 2 Background

## 2.1 BC-IFS model

The model is based on an Iterated Function System (IFS) [BD85], which consists of a set of $N$ contractive transformations $T_i$ applied recursively. Each iteration approximates the resulting fractal object, known as the *attractor*, which is obtained after an infinite number of iterations. The attractor $A$ is the unique compact set satisfying $A = \bigcup_{i=0}^{N-1} T_i(A)$; it is composed of smaller copies of itself.

By defining the attractor in the barycentric space $BI^n$, it can be projected into the modeling space $(\mathbb{R}^2$ or $\mathbb{R}^3)$ using a set of $n$ control points. It allows for a global deformation of the attractor in the modeling space while preserving its topology. In the barycentric space, the transformations are encoded by square matrices of size $n \times n$ such that any point in $BI^n$ (a hyperplane in $\mathbb{R}^n$) is still in $BI^n$ after being transformed. It implies that the sum of the elements in each column must be 1.

The BC-IFS model controls the iterative process using an automaton. Each state represents a topological cell of the structure. Transitions between states are of two kinds: *subdivision* and *boundary* operators. Subdivision operators (denoted by $\div$) specify how a cell subdivides into several smaller sub-cells, corresponding to the application of IFS transformations. Boundary operators (denoted by $\partial$) specify a barycentric subspace of topological cells (i.e., a subset of control points). The automaton encodes the complete iterative process, where a type of cell can be subdivided into different types of cells, each with its own subdivision rules.

The BC-IFS model extends the C-IFS model by incorporating *incidence* and *adjacency* constraints. Incidence constraints are used to specify the cellular decomposition. These constraints take the form of equations involving subdivision and boundary operators. Solving these equations imposes a relationship between subdivision operators, meaning that certain values in one transformation matrix must equal corresponding values in another. Identical values in the matrices correspond to the same position in the barycentric space and, consequently, the same position in the modeling space. Figure 1 shows an example of a fractal created using the BC-IFS model.

In summary, the conception of a fractal structure with the BC-IFS model requires the definition of:

* The automaton, to define how each cell subdivides.

* Incidence constraints, which specify the cellular decomposition.

* Adjacency constraints, which define the topology.

For more details, the reader can refer to [GGS21]. An example of a fractal structure is shown in Figure 1. It consists of a hexagon subdivided into six smaller copies of itself recursively. Figure 2 illustrates the cellular decomposition of the shape, including a face, its edges, and its vertices. The edges have specific topologies: *Bézier curves* (in blue) and *Cantor sets* (in red).
<figure><img src="image_3.png" /><figcaption>Figure 1: Example of fractal created with the BC-IFS model. The model allows the deformation of fractal structures without changing the topology. (Top left) Control points in red are planar, and subdivision points in blue induce circular gaps. (Top Right) Subdivision points control the local geometry. (Bottom left) Control points control the global geometry. (Bottom right) Moving control points and subdivision points changes the global and local geometry.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 2: Cellular decomposition of a fractal face. (Left) A face has edges as boundaries, which have vertices as boundaries. (Center) Subdivision of the face into six sub-faces respecting the boundary subdivision. (Right) After a few iterations, approximation of the attractor. In black, four subdivision steps of Cantor and Bézier edges. Cantor edges subdivide into two disconnected sub-edges. Bézier edges subdivide into two connected sub-edges.</figcaption></figure>

All automaton states are denoted as follows: *F* for the face, *B* and *C* for Bézier and Cantor edges, respectively, and *V* for vertices. The face *F* has six boundary operators, $\partial_0^F, \dots, \partial_5^F$, one for each edge. Edges have two boundary operators each: $\partial_0^B$ and $\partial_1^B$ for *B*, and $\partial_0^C$ and $\partial_1^C$ for *C*. The boundaries of edges are vertices. Vertices have no boundary operators, as they are the smallest topological cells.

The face subdivides into six copies of itself. Therefore, there are six subdivision operators, $\div_i^F$, one
for each sub-face $F_i$. Both Bézier and Cantor edges have two subdivision operators each: $\div_0^B$ and $\div_1^B$ for $B$, and $\div_0^C$ and $\div_1^C$ for $C$. Vertices subdivide into themselves with one subdivision operator, $\div_0^V$, and behave topologically as a point. The automaton for this example is shown on the left of Figure 3.

<figure><img src="image_5.png" /><figcaption>Figure 3: The automaton of the fractal hexagon in Figure 2 with its incidence (middle) and adjacency (right) constraints.</figcaption></figure>

The automaton alone does not control the topology of the fractal. To address this, we need incidence and adjacency constraints. In our example, there are twelve incidence constraints and six adjacency constraints. All of these are listed in Figure 3 (middle and right). Incidence constraints ensure that sub-faces align with the subdivision of their face edges. For instance, Equation 1 shows an example of an incidence constraint, highlighted in green in Figure 4. This constraint specifies that the sub-edge 0 of the edge 0 must match the edge 0 of the sub-face 2.

$$ \partial_0^F \div_0^C = \div_2^F \partial_0^F \quad (1) $$

Adjacency constraints, on the other hand, ensure that sub-faces connect properly by sharing an edge. Equation 2 provides an example of an adjacency constraint, represented in orange in Figure 4. This constraint states that the edge 2 of the sub-face 2 must correspond to the edge 4 of the sub-face 3.

$$ \div_2^F \partial_2^F = \div_3^F \partial_4^F \quad (2) $$

For simplicity, we define Bézier curves using three control points, Cantor sets using two control points, and vertices using one control point. The subdivision matrix of the vertex state $\div_0^V$ is the one by one identity matrix. We set the subdivision operators of $B$ and $C$ as follows:

$$ \div_0^B = \begin{pmatrix} 1 & 1/2 & 1/4 \\ 0 & 1/2 & 1/2 \\ 0 & 0 & 1/4 \end{pmatrix} \qquad \div_1^B = \begin{pmatrix} 1/4 & 0 & 0 \\ 1/2 & 1/2 & 0 \\ 1/4 & 1/2 & 1 \end{pmatrix} \qquad \div_0^C = \begin{pmatrix} 1 & 2/3 \\ 0 & 1/3 \end{pmatrix} \qquad \div_1^C = \begin{pmatrix} 1/3 & 0 \\ 2/3 & 1 \end{pmatrix} $$

When applying the constraints, we obtain all subdivision operators $\div_i^F$ for $F$, which include both fixed and free values (denoted by ●). For instance, the matrices for $\div_2^F$ and $\div_3^F$ are:

$$ \div_2^F = \begin{pmatrix} 1 & 2/3 & \bullet & \bullet & \bullet & 1/4 & 1/2 \\ 0 & 1/3 & \bullet & \bullet & \bullet & 0 & 0 \\ 0 & 0 & \bullet & \bullet & \bullet & 0 & 0 \\ 0 & 0 & \bullet & \bullet & \bullet & 0 & 0 \\ 0 & 0 & \bullet & \bullet & \bullet & 0 & 0 \\ 0 & 0 & \bullet & \bullet & \bullet & 0 & 0 \\ 0 & 0 & \bullet & \bullet & \bullet & 1/2 & 0 \\ 0 & 0 & \bullet & \bullet & \bullet & 1/4 & 1/2 \end{pmatrix}, \div_3^F = \begin{pmatrix} 1/3 & 0 & 0 & 0 & \bullet & \bullet & \bullet & \bullet \\ 2/3 & 1 & 1/2 & 1/4 & \bullet & \bullet & \bullet & \bullet \\ 0 & 0 & 1/2 & 1/2 & \bullet & \bullet & \bullet & \bullet \\ 0 & 0 & 0 & 1/4 & \bullet & \bullet & \bullet & \bullet \\ 0 & 0 & 0 & 0 & \bullet & \bullet & \bullet & \bullet \\ 0 & 0 & 0 & 0 & \bullet & \bullet & \bullet & \bullet \\ 0 & 0 & 0 & 0 & \bullet & \bullet & \bullet & \bullet \\ 0 & 0 & 0 & 0 & \bullet & \bullet & \bullet & \bullet \end{pmatrix} $$
<figure><img src="image_6.png" /><figcaption>Figure 4: Illustration of incidence and adjacency constraints and their impact on the fractal structure's topology. In green, the incidence constraint in Equation 1 ensures the sub-face $F_2$ follows the subdivision of the first edge of the face. In orange, the adjacency constraint in Equation 2 ensures the two sub-faces $F_2$ and $F_3$ are connected by sharing an edge.</figcaption></figure>

Due to the incidence constraints, the subdivision matrices of the edges appear as sub-matrices within the subdivision matrices of $F$. It guarantees that the attractor of the edges are sub-attractor of faces. Since the transformations must be contractive, the sum of the values in each column must equal 1. Each column in the matrix represents the position of a subdivision point in the barycentric space (shown in blue in Figure 1). These subdivision points are the images of control points transformed by the matrices during the first subdivision step.

Modifying the matrix values changes the position of the subdivision points in the barycentric space, which locally deforms the geometry. Due to the self-similarity of the structure, these changes propagate recursively. Adjacency constraints ensure consistency by enforcing that some columns of one matrix match specific columns of another. For example, the blue (respectively, gray) column of $\div_2^F$ is identical to the blue (respectively, gray) column of $\div_3^F$. These blue and gray columns represent the shared edge between the sub-faces $F_2$ and $F_3$.

We employ a mass-spring system to generate a default geometry in the barycentric space. This system provides an initial layout for the subdivision points, serving as a starting point for further adjustments.

## 2.2 Polyhedral crystallographic circle packings

This section provides the necessary background for understanding polyhedral crystallographic circle packings and the related structures used throughout this paper. For further details, see [RAR25, Ras21]. For any point $v \in \mathbb{R}^3$ outside the unit sphere $\mathbb{S}^2$ (i.e. with $\|v\| > 1$), the *stereographic circle* of $v$ is defined as the image, under the stereographic projection, of the spherical cap

$$C_v = \{\langle u, v \rangle \ge 1 \mid u \in \mathbb{S}^2\},$$

where $\langle \cdot, \cdot \rangle$ denotes the Euclidean inner product. Geometrically, $C_v$ can be interpreted as the illuminated region on the unit sphere produced by a light source located at the point $v$.

We refer to a (convex) *polyhedron* $\Pi$ as the convex hull of four or more points of $\mathbb{R}^3$ in general position. The 0, 1 and 2-cells of $\Pi$ are called vertices, edges, and faces of $\Pi$, respectively. We label each vertex by a cardinal number and denote each edge and face by the concatenation of its vertices. The *combinatorial structure* of $\Pi$ is the lattice of all its dimensional cells ordered by inclusion relation. This lattice includes the empty set and the entire polyhedron, and its elements can be stratified by level according to dimension (see Figure 5).
<figure><img src="image_7.png" /><figcaption>Figure 5: A tetrahedron (left) with the cell notation, displayed next to its combinatorial structure (right).</figcaption></figure>

The *arrangement projection* of a polyhedron $\Pi$, whose vertices all lie outside the unit sphere, is the collection of stereographic circles corresponding to its vertices. For any edge $uv$ of $\Pi$, the spherical caps $C_u$ and $C_v$ are disjoint, tangent, or intersect with overlapping interiors if and only if the edge $uv$ intersects $S^2$ transversely, is tangent to it, or avoids it entirely, respectively. When $\Pi$ is edge-scribed (all its edges are tangent to $S^2$), its arrangement projection forms a *polyhedral circle packing* $C_{\Pi}$ (see Figure 6).

<figure><img src="image_8.png" /><figcaption>Figure 6: Illuminated region of three cubes with a different geometric realization (top figures) and their arrangement projection (bottom figures). When the cube is edge-scribed, its arrangement projection is a polyhedral circle packing.</figcaption></figure>

A polyhedral circle packing inherits the combinatorial structure of the corresponding polyhedron. The vertices and edges of $\Pi$ correspond bijectively to the circles and tangency relations of $C_{\Pi}$. The intersections of the faces of $\Pi$ with the unit sphere define the *dual circles* of $C_{\Pi}$. By inverting $C_{\Pi}$ through each of its dual circles and iterating this process ad infinitum, one obtains a *polyhedral crystallographic circle packing* [KN19, CRCS20] (see Figure 7).

Every polyhedron admits an edge-scribed realization inducing a polyhedral crystallographic circle packing. However, this property does not extend to higher dimensions [Ras24]. In the case of the
tetrahedron, the resulting structure is the well-known fractal known as the *Apollonian packing* (see Figure 8).

<figure><img src="image_9.png" /><figcaption>Figure 7: Example of crystallographic circle packing. (Left) A cubic circle packing $\mathcal{C}_\Pi$ with its dual circles (in dashed). (Middle) The obtained packing by inverting $\mathcal{C}_\Pi$ through its dual circles. (Right) The cubic crystallographic circle packing when iterating the inversions.</figcaption></figure>

### 3 The BC-IFS model of polyhedral crystallographic circle packings

The primary objective of this paper is to construct fractal structures that are topologically equivalent to polyhedral crystallographic circle packings using the BC-IFS model. This approach enables controlled deformations of the structure's shape. To achieve this, we start with the complement of the packing on the sphere, which is then subdivided into smaller sub-regions (see Figure 8).

<figure><img src="image_10.png" /><figcaption>Figure 8: The complementary regions of the initial packing and the first four iterations on the sphere of a tetrahedral crystallographic circle packing, commonly known as the Apollonian packing.</figcaption></figure>

When projecting the complementary regions onto the plane, we can choose the projection center to lie inside one of the packing circles. This choice results in a bounded complement, which induces bounded regions (see Figure 9).

We consider each inversion as a subdivision process which creates sub-regions. The formation of sub-regions directly results from these inversions. Intuitively, dual circles invert each region, so every sub-region corresponds to an inverted version of an existing region, as illustrated in Figure 10. This property reflects the general behavior of the system, where each region subdivides into all other regions.

The subdivision rules rely solely on the combinatorial structure of the polyhedron. This dependence allows us to fully automate the design process of the BC-IFS model when conceiving the topology of polyhedral crystallographic circle packings. In the following sections, we provide a detailed explanation of how each step of this automated design process works.

#### 3.1 States of the automaton

For each polyhedron $\Pi$, the states of the automaton $A_\Pi$ in the BC-IFS model that generates the corresponding lacunary fractal structure are the cells of the *primitive* shape $P_\Pi$, which we define as the complement of the interior of a polyhedral circle packing $\mathcal{C}_\Pi$ (see Figure 8, left). The primitive is topologically equivalent to a polygonal complex composed of a collection of polygons where the intersections are either empty or consist of a single vertex. The 0, 1, 2-cells of $P_\Pi$ will also be called
<figure><img src="image_11.png" /><figcaption>Figure 9: We consider inversion as subdivision. (Left) The stereographic projection onto the plane of a tetrahedral circle packing and its dual circles (in dashed), with its complementary regions (in gray). (Middle) Same configuration after one iteration of inversion through each dual circle. (Right) The subdivision process of one complementary region. In this case, each region subdivides into three sub-regions.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure 10: A vertex-centered projection of primal and dual circles of an edge-scribed tetrahedron with colored areas we consider as material. After inversion, each sub-area is an inverted area.</figcaption></figure>

vertices, edges, and faces, respectively. These cells can be identified in terms of the cells of the corresponding polyhedron by the following bijections:

(i) Each vertex $V_e$ of the primitive $P_Π$ corresponds to a unique edge $e$ of the polyhedron $Π$.

(ii) Each edge $E_v^f$ of the primitive $P_Π$ corresponds to a unique vertex-face pair $(v, f)$ of the polyhedron $Π$, where $v \in f$.

(iii) Each face $F_f$ of the primitive $P_Π$ corresponds to a unique face $f$ of the polyhedron $Π$.

Therefore, the states of the automaton $A_Π$ consist of a *vertex state* for each edge of $Π$, an *edge state* for each vertex within each face of $Π$, and a *face state* for each face of $Π$. The left side of Figure 11 illustrates an example of encoding of the cells of the primitive $P_Π$ when $Π$ is a tetrahedron, using the concatenation of the vertex names as described before.

## 3.2 Boundary operators

After defining the states of the automaton $A_Π$, we now define the transitions between the states, beginning with the boundary operators. For each face $f$ of $Π$ and for each vertex $v$ of $f$, we define a boundary operator from the face states to the edge states of $A_Π$ as

$$ \partial_v^{F_f} : F_f \mapsto E_v^f \qquad (3) $$
<figure><img src="image_13.png" /><figcaption>Figure 11: Example of encoding of the cells of the primitive $P_{\Pi}$ when $\Pi$ is a tetrahedron. (Left) Association between polyhedral and fractal cells on one face of the tetrahedron. (Right) After one iteration, the polyhedral topology appears in sub-faces due to inversions.</figcaption></figure>

For each face $f$ of $\Pi$ and for each edge $uv$ of $f$, we define a boundary operator from the edge states to the vertex states of $A_{\Pi}$ as

$$ \partial_{u}^{E_{v}^{f}} : E_{v}^{f} \mapsto V_{uv} \qquad (4) $$

Figure 12 illustrates an example of the boundary operators on the states of the cells of a single face $F_{123}$ when $\Pi$ is a tetrahedron. The full automaton with all the boundary operators is depicted in Figure 15.

<figure><img src="image_14.png" /><figcaption>Figure 12: A part of the automaton $A_{\Pi}$ including the boundary operators of a face $F_{123}$, when $\Pi$ is a tetrahedron.</figcaption></figure>

## 3.3 Subdivisions operators

Every face $F_f$ of the primitive subdivides into all the other faces $F_{f'}$. Hence, for every pair $f, f'$ of distinct faces of the polyhedron $\Pi$, we define the subdivision operators on the face states of $A_{\Pi}$ as

$$ \div_{f'}^{F_f} : F_f \mapsto F_{f'} \qquad (5) $$

Every edge $E_v^f$ of the primitive subdivides into the other edges $E_v^{f'}$ corresponding to the same vertex $v$ of the polyhedron (see Figure 13). Hence, for every pair of distinct faces $f, f'$ of the polyhedron $\Pi$ sharing a vertex $v$, we define a subdivision operator on the edge states of $A_{\Pi}$ as

$$ \div_{f'}^{E_v^f} : E_v^f \mapsto E_v^{f'} \qquad (6) $$
<figure><img src="image_15.png" /><figcaption>Figure 13: Example of edge subdivision. (Left) Highlighted edges share the same associated polyhedral vertex. (Right) The edge $E_1^{123}$ subdivides into the edges $E_1^{124}$ and $E_1^{134}$ sharing the same associated polyhedral vertex (the vertex 1), preserving the adjacency.</figcaption></figure>

Figure 14 describes the subdivision operators of the face states and edge states when $\Pi$ is a tetrahedron. Figure 15 shows the full automaton with all the subdivision operators.

<figure><img src="image_16.png" /><figcaption>Figure 14: Parts of the automaton $A_{\Pi}$ with subdivision operators applied to all the face states (left) and to the edge states corresponding to the same vertex of $\Pi$ (right), when $\Pi$ is the tetrahedron.</figcaption></figure>

As mentioned in Section 2.1, vertices play a special role when dealing with the fractal structures produced by the BC-IFS model. As shown in Figure 11, all vertices of the fractal structure remain unchanged after one iteration, subdividing into themselves. Therefore, for each edge $e$ of the polyhedron $\Pi$, we define a subdivision operator on the vertex states as follows:

$$ \div_{0}^{V_e} : V_e \mapsto V_e \qquad (7) $$

The states, along with all boundary and subdivision operators, form the complete automaton $A_{\Pi}$ of the BC-IFS model (see Figure 15 for the case of the tetrahedron).
<figure><img src="image_17.png" /><figcaption>Figure 15: The automaton $\mathcal{A}_{\Pi}$ with all the states, boundary operators (black arrows), and subdivision operators (dashed arrows) when $\Pi$ is the tetrahedron.</figcaption></figure>

We recall that while the automaton governs the subdivision process, it does not determine the topology of the attractor. To ensure topological consistency, we must define incidence and adjacency constraints, which can be derived from the combinatorial structure of the polyhedron.

## 3.4 Incidence constraints

Each vertex of the polyhedron $\Pi$ corresponds to a set of edges of the primitive $P_{\Pi}$. After one iteration, each of these edges subdivides into others while preserving adjacency (see Figure 13).

An incidence constraint ensures that every sub-edge of every edge is equal to an edge of the corresponding sub-face of the primitive. By combining this condition with the bijections between the cells of the primitive and those of the polyhedron, we derive algorithm 1 which lists all the incidence constraints.

## 3.5 Adjacency constraints

Unlike the hexagon example in Section 2.1, where two sub-faces are connected by a sub-edge, in this case, sub-faces are connected by a vertex, which corresponds to the edge of the polyhedron joining the corresponding faces. That means that each constraint involves two successive boundary operators:
<table><caption>**Algorithm 1** Algorithm to define all incidence constraints.</caption><tr><td>for</td><td>each face</td><td>f</td><td>of</td><td>Π</td><td>do</td></tr><tr><td>for</td><td>each vertex</td><td>v</td><td>of</td><td>f</td><td>do</td></tr><tr><td>for</td><td>each face</td><td>f'</td><td>≠</td><td>f</td><td>of</td><td>Π</td><td>containing</td><td>v</td><td>do</td></tr><tr><td>Define incidence constraint</td><td>∂<sub>v</sub><sup>F<sub>f</sub></sup> ÷<sub>f'</sub><sup>E<sub>v</sub><sup>f'</sup></sup> = ÷<sub>f'</sub><sup>F<sub>f</sub></sup> ∂<sub>v</sub><sup>F<sub>f'</sub></sup></td></tr><tr><td>end for</td></tr><tr><td>end for</td></tr><tr><td>end for</td></tr></table>

one to obtain an edge of a sub-face and another to get a vertex of that edge. Combining this condition with the bijections between the primitive cells and those of the polyhedron allows us to derive the algorithm 2, which lists all the adjacency constraints.

<table><caption>**Algorithm 2** Algorithm to define all adjacency constraints.</caption><tr><td>for</td><td>each face</td><td>f</td><td>of</td><td>Π</td><td>do</td></tr><tr><td>for</td><td>each edge</td><td>uv</td><td>∈</td><td>f</td><td>of</td><td>Π</td><td>do</td></tr><tr><td>Let</td><td>f', f''</td><td>be the two adjacent faces of Π</td><td>sharing</td><td>uv</td></tr><tr><td>Define adjacency constraint</td><td>÷<sub>f'</sub><sup>F<sub>f</sub></sup> ∂<sub>u</sub><sup>F<sub>f'</sub></sup> ∂<sub>v</sub><sup>E<sub>u</sub><sup>f'</sup></sup> = ÷<sub>f''</sub><sup>F<sub>f</sub></sup> ∂<sub>u</sub><sup>F<sub>f''</sub></sup> ∂<sub>v</sub><sup>E<sub>u</sub><sup>f''</sup></sup></td></tr><tr><td>end for</td></tr><tr><td>end for</td></tr></table>

# 4 Results

## 4.1 Automatic extraction of crystallographic circle packings topology

Our method automatically extracts the polyhedral crystallographic circle packing topology and encodes it using the BC-IFS model. It deduces the automaton and all incidence and adjacency constraints from the combinatorial structure of the polyhedron. Figure 16 illustrates the fractal structure obtained from a tetrahedron, where the topology of each face mirrors the Sierpinski triangle (a triangle subdividing into three smaller triangles).

<figure><img src="image_18.png" /><figcaption>Figure 16: First four iterations of the tetrahedral crystallographic circle packing material area. There are four fractal faces, each with the same topology as the Sierpinski triangle.</figcaption></figure>

Our approach also applies seamlessly to irregular polyhedra, as the automaton states distinguish each face of the polyhedron. For instance, Figure 17 shows the resulting fractal shape when the polyhedron is a square-based pyramid. Each triangular face subdivides into four sub-faces (three triangles and one quadrilateral) while the quadrilateral face subdivides into four smaller triangles.

By focusing in one of the faces of the primitives, we obtain a *lacunar fractal polygon* (see Figure 18 in the case of the cube). The full process of topology extraction is illustrated in the video available in the supplemental materials.
<figure><img src="image_19.png" /><figcaption>Figure 17: First four iterations of a square pyramidal crystallographic circle packing. (Top line) The complementary of the circles is visible on the sphere. (Bottom line) The extracted topology of the material area.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Figure 18: First four iterations of the cubic crystallographic circle packing. (Top line) The complementary of the circles is visible on the sphere. (Bottom line) A lacunar fractal square corresponding to one face.</figcaption></figure>

## 4.2 Application for a CAD editor of lacunar fractal polygons

Lacunar fractal polygons enables a straightforward design of fractal structures assembling several of them, each with its unique topology. We propose the following workflow for creating such structure in a CAD editor:

* Define the domain boundary.

* Partition the domain into several polygons.

* For each polygon, assign a suitable polyhedron.

Each initial polygon is a face of a primitive, which corresponds to the dual of the faces of the corresponding polyhedron; see Figure 18. Since polygons are self-dual, a suitable polyhedron is a polyhedron containing the polygon as one of its faces.

On the other hand, maintaining topological consistency on the edges shared by two assembled polygons imposes a constraint on face selection. Any two connected polygons must share a fractal edge with the same topology. This implies that the corresponding polyhedra are connected by a vertex. Since the edge subdivision of the primitive depends on the valency of the corresponding vertex of the polyhedra, the valency of both corresponding polyhedra at this vertex must be equal.

The selection of polyhedra for each polygon can be done by drawing its graph in the dual of the polygon, as shown in Figure 19. The boundary face of the graph represents the face we use for automatic topology extraction.
End-users can meet the valency constraint by following a simple design rule: two polygons can share an edge only when both have the same number of graph edges along the shared boundary.

<figure><img src="image_21.png" /><figcaption>Figure 19: The selected face in red on the polyhedron, whose graph is drawn for a polygon, with its corresponding fractal shape we can deform.</figcaption></figure>

For instance, the structure in Figure 20, shaped like the letter F, assembles several polygons while adhering to this design rule and preserving topological consistency.

<figure><img src="image_22.png" /><figcaption>Figure 20: From left to right: a cutting of a structure in the shape of the letter F; a drawing of a polyhedron graph on each cell ensuring the respect of the constraint to have a topologically consistent structure; first three iterations of the fractal structure.</figcaption></figure>

# 5 Conclusion and future work

This work provides a method for computer-aided design of fractal structures with an independent control of geometry and topology: an aspect that has traditionally been difficult to achieve.

We proposed a method to construct lacunar fractal polygons with the topology of polyhedral crystallographic circle packings using the BC-IFS model. The topology depends on the chosen polyhedron. This methodology allows designers to easily construct complex fractal structures by assembling lacunar fractal polygons and selecting a polyhedron for each polygon.

Additionally, the BC-IFS model allows precise control over local and global deformations through subdivision and control points while preserving topology, as shown in Figure 21.

Future research could extend this approach to three dimensions. Analogues of polyhedral crystallographic circles packings exist in higher dimensions, but dual spheres intersect in dimensions higher than
<figure><img src="image_23.png" /><figcaption>Figure 21: (Left) Fractal structure in the shape of the letter F. (Middle) Local deformation using subdivision points. (Right) Global deformation using control points.</figcaption></figure>

two [Ras24]. This intersection eliminates distinct material regions (the sphere packing complementary) and complicates the identification of subdivision rules.

A possible approach involves treating inversions as polytope subdivisions. Figure 22 illustrates the subdivision of a tetrahedron through successive inversion steps. Each newly introduced circle generates a corresponding vertex on the polyhedron, and the polyhedral graph is identical to the tangency graph of the circle packing.

<figure><img src="image_24.png" /><figcaption>Figure 22: First four iterations of the crystallographic subdivision on a tetrahedron.</figcaption></figure>

Future work could focus on extracting the topology of iterative sphere packings derived from modifying the subdivision rules of the polytopes. One could define a polytope and its subdivision rules, perform a subdivision step, recompute an edge-scribed realization if it is possible, and generate the corresponding sphere packing. The topology of such packings depends exclusively on the subdivision rules, which could simplify the extension to three dimensions.

## References

[AALS17] Niels Aage, Erik Andreassen, Boyan S. Lazarov, and Ole Sigmund. Giga-voxel computational morphogenesis for structural design. *Nature*, 550(7674):84–86, October 2017.

[ABC$^{+}$19] Pablo Antolin, Annalisa Buffa, Elaine Cohen, John F. Dannenhoffer, Gershon Elber, Stefanie Elgeti, Robert Haimes, and Richard Riesenfeld. Optimizing Micro-Tiles in Micro-Structures as a Design Paradigm. *Computer-Aided Design*, 115:23–33, October 2019.

[BD85] Michael F Barnsley and Stephen Demko. Iterated function systems and the global construction of fractals. *Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences*, 399(1817):243–275, 1985.

[BEH89] Michael F Barnsley, John H Elton, and Douglas P Hardin. Recurrent iterated function systems. *Constructive approximation*, 5:3–31, 1989.

[BG24] Boris Bordeaux and Christian Gentil. Automatic construction of fractal structures with locally controlled lacunarity. *Journal of WSCG*, 32(1-2):1–12, 2024.
[CRCS20] Devora Chait-Roth, Alisa Cui, and Zachary Stier. A taxonomy of crystallographic sphere packings. *Journal of Number Theory*, 207:196-246, 2020.

[GGS21] Christian Gentil, Gilles Gouaty, and Dmitry Sokolov. *Geometric Modeling of Fractal Forms for CAD*. John Wiley & Sons, 2021.

[GSR22] David-Henri Garnier, Martin-Pierre Schmidt, and Damien Rohmer. Growth of oriented orthotropic structures with reaction/diffusion. *Structural and Multidisciplinary Optimization*, 65(11):327, November 2022.

[HE21] Q Youn Hong and Gershon Elber. Conformal Microstructure Synthesis in Trimmed Trivariate Based V-Reps. *Computer-Aided Design*, 140:103085, November 2021.

[KN19] Alex Kontorovich and Kei Nakamura. Geometry and arithmetic of crystallographic sphere packings. *Proceedings of the National Academy of Sciences*, 116(2):436-441, 2019.

[PHHM96] Przemysław Prusinkiewicz, Mark Hammel, Jim Hanan, and Radomir Mech. L-systems: from the theory to visual models of plants. In *Proceedings of the 2nd CSIRO Symposium on Computational Challenges in Life Sciences*, volume 3, pages 1-32. Citeseer, 1996.

[Pru86] Przemysław Prusinkiewicz. Graphical applications of l-systems. *Proceedings of Graphics Interface '86 / Vision Interface '86*, pages 247-253, 1986.

[RAR25] Jorge L. Ramírez Alfonsín and Iván Rasskin. A polytopal generalization of apollonian packings and descartes' theorem. *Journal of Number Theory*, 275:67-103, 2025.

[Ras21] Iván Rasskin. *A polytopal approach to Apollonian packings and discrete knotted structures*. PhD thesis, Université de Montpellier, 2021.

[Ras24] Iván Rasskin. Regular polytopes, sphere packings and apollonian sections. *Geometriae Dedicata*, 218(6):105, 2024.

[SGGM15] Dmitry Sokolov, Gilles Gouaty, Christian Gentil, and Anton Mishkinis. Boundary Controlled Iterated Function Systems. In Jean-Daniel Boissonnat, Albert Cohen, Olivier Gibaru, Christian Gout, Tom Lyche, Marie-Laurence Mazure, and Larry L. Schumaker, editors, *Curves and Surfaces*, volume 9213, pages 414-432. Springer International Publishing, Cham, 2015. Series Title: Lecture Notes in Computer Science.

[SYNP19] Meera Sitharam, Jeremy Youngquist, Maxwell Nolan, and Jörg Peters. Corner-sharing tetrahedra for modeling micro-structure. *Computer-Aided Design*, 114:164-178, 2019.

[TT93] Joëlle Thollot and Eric Tosan. Construction of fractals using formal languages and matrices of attractors. In *Conference on Computational Graphics and Visualization Techniques (Compugraphics)*, pages 74-81, 1993.

[WSG21] Jun Wu, Ole Sigmund, and Jeroen P. Groen. Topology optimization of multi-scale structures: a review. *Structural and Multidisciplinary Optimization*, 63(3):1455-1480, March 2021.

[ZL24] Qiang Zou and Guoyue Luo. Geometric modeling for microstructure design and manufacturing: A review of representations and modeling algorithms. *Computer-Aided Design*, page 103834, 2024.

[ZLL+23] Tianyu Zhou, Hengxu Li, Xinming Li, Carlos F. Lange, and Yongsheng Ma. Feature-based modeling for variable fractal geometry design integrated into CAD system. *Advanced Engineering Informatics*, 57:102006, August 2023.