step. We choose a length scale $\varepsilon$ for this projection step and reconstruct a periodic shape with length scale $\varepsilon$, approximating the optimal one.

### 6.6.1 Projection in the simple case without varying orientation, $\alpha \equiv 0$

First, we consider the case where there is no varying orientation, that is, $\alpha \equiv 0$. The unit cells (rectangular hole in a square) are defined by

$$ Y(m) = \left\{ \begin{array}{ll} \cos(2\pi y_1) \ge \cos(\pi(1 - m_1)), \\ y \in [0, 1]^2 : & \text{or} \\ \cos(2\pi y_2) \ge \cos(\pi(1 - m_2)) \end{array} \right\}. $$

The domain $D$ is paved with cells $\varepsilon Y(m)$. Since the hole size $m(x)$ is varying in $D$, the periodicity cell is macroscopically modulated and we define a projected lattice shape $\Omega_{\varepsilon}(m)$ as

$$ \Omega_{\varepsilon}(m) := \left\{ \begin{array}{ll} \cos\left(\frac{2\pi x_1}{\varepsilon}\right) \ge \cos(\pi(1 - m_1(x))), \\ x \in D : & \text{or} \\ \cos\left(\frac{2\pi x_2}{\varepsilon}\right) \ge \cos(\pi(1 - m_2(x))) \end{array} \right\}, $$

where $m_1(x), m_2(x)$ are functions defined on $D$ with values in $[0, 1]$.

**Remark 6.1.** The values of $m_1, m_2$ are not necessarily constant in each cell of the structure. Hence, the holes in the cellular structure $\Omega_{\varepsilon}(m)$ are not exactly rectangles. But, when $\varepsilon$ goes to 0, the sequence of cellular structures converges to the composite with local Hooke's law equal to $A_0^*(m)$.

The cellular structures can be defined using level-sets. We introduce two functions $f_{\varepsilon,i}^m$, one for each direction

$$ f_{\varepsilon,i}^m(x) := -\cos\left(\frac{2\pi x_i}{\varepsilon}\right) + \cos(\pi(1 - m_i(x))), $$

and the level-set function

$$ F_{\varepsilon}^{m} := \min(f_{\varepsilon,1}^{m}, f_{\varepsilon,2}^{m}). $$
The final structure $\Omega_{\varepsilon}(m)$ is then defined by

$$ \Omega_{\varepsilon}(m) = \{x \in D : F_{\varepsilon}^{m}(x) \leq 0\}. $$

The construction of a minimizing sequence of shapes is immediate: we just have to update the size $\varepsilon$ in the previous level-set function.

### 6.6.2 Projection in the general case with orientation, $\alpha \neq 0$

This section is based on the paper of Pantz and Trabelsi [PT2008]. The main idea of Pantz and Trabelsi is to find a map $\varphi = (\varphi_1, \varphi_2)$ from $D$ into $\mathbb{R}^2$ which distorts the regular square grid in order to orientate each square at the optimal angle $\alpha$ (or $\alpha + \pi$). Once $\varphi$ is found, we can proceed as before.

The final shape, now denoted $\Omega_{\varepsilon}(\varphi, m)$, is still defined by a level-set function:

$$ \Omega_{\varepsilon}(\varphi, m) = \{x \in D : F_{\varepsilon}^{\varphi, m}(x) \leq 0\} $$

with $F_{\varepsilon}^{\varphi, m} = \min(f_{\varepsilon,1}^{\varphi, m}, f_{\varepsilon,2}^{\varphi, m})$ and

$$ f_{\varepsilon,i}^{\varphi,m}(x) = -\cos\left(\frac{2\pi\varphi_i(x)}{\varepsilon}\right) + \cos(\pi(1 - m_i(x))). $$

Geometrically (in 2-d), the Jacobian $\nabla\varphi$ should be proportional to the rotation matrix defined by

$$ Q(\alpha) = \begin{pmatrix} \cos \alpha & -\sin \alpha \\ \sin \alpha & \cos \alpha \end{pmatrix}. \qquad (6.6) $$

One possibility (Pantz and Trabelsi) is to find $\varphi$ which (roughly) minimizes

$$ \int_D |\nabla\varphi - Q(\alpha)|^2 dx. $$

However, it is not straightforward because $Q(\alpha)$ is not smooth (the orientation is not coherent a priori). Pantz and Trabelsi proposed a complicated trick to avoid this coherent orientation issue. We also mention that Groen and Sigmund [GS2018] suggested another trick from image processing to obtain a coherent orientation.
### 6.6.3 The new approach of Allaire-Geoffroy-Pantz

In this section, we propose a new approach, based on the paper of Allaire-Geoffroy-Pantz [AGP2018].

Let us recall that geometrically (in 2-d), at every point $x \in D$, the Jacobian $\nabla\varphi$ should be proportional to the rotation matrix $Q(\alpha)$ defined by (6.6). Moreover, the proportions of the cell have to be preserved in order to converge to a true square and not simply to a rectangle. For this purpose, we impose $|\nabla\varphi_1| = |\nabla\varphi_2| = e^r$, where $r \in H^1(D)$ is a scalar dilation field. Then the Jacobian $\nabla\varphi$ should be

$$ \nabla\varphi = e^r Q(\alpha) \quad \text{in } D. \qquad (6.7) $$

This equation can be satisfied only if $\alpha$ satisfies the following conformality condition.

**Lemma 6.2.** Let $\alpha$ be a regular orientation field and $D$ be a simply connected domain. There exists a mapping function $\varphi$ and a dilation field $r$ satisfying $\nabla\varphi = e^r Q(\alpha)$ if and only if

$$ \Delta\alpha = 0 \quad \text{in } D. \qquad (6.8) $$

We recall that, for a vector field $u = (u_1, u_2)$, its curl is defined as $\mathrm{curl} \, u = \nabla \wedge u = \frac{\partial u_2}{\partial x_1} - \frac{\partial u_1}{\partial x_2}$, where $\wedge$ is the 2-d cross product of vectors. Of course, $\mathrm{curl} \, \nabla\varphi = (\mathrm{curl} \, \nabla\varphi_1, \mathrm{curl} \, \nabla\varphi_2) = 0$.

*Proof of Lemma 6.2.* Since the domain $D \subset \mathbb{R}^2$ is simply connected, By Poincaré's lemma, the map $\varphi$ exists if and only if the right hand side of (6.7) is curl-free. That is,

$$ \mathrm{curl} \, (e^r Q(\alpha)) = 0. $$

Let $a_1, a_2$ be the columns of $Q(\alpha)$. Then

$$ \mathrm{curl} \, (e^r Q(\alpha)) = 0 \iff \nabla r \wedge a_i = -\nabla \wedge a_i, \quad i = 1, 2. \qquad (6.9) $$

Since, for fixed $\alpha$, $(a_1, a_2)$ is an orthonormal basis of $\mathbb{R}^2$, by (6.9) we see that the vector $\nabla r$ can be decomposed as

$$ \nabla r = (-\nabla \wedge a_2)a_1 + (\nabla \wedge a_1)a_2. $$
We compute

$$ \nabla \wedge a_1 = \frac{\partial \alpha}{\partial x_1} \cos(\alpha) + \frac{\partial \alpha}{\partial x_2} \sin(\alpha) \quad \text{and} \quad \nabla \wedge a_2 = -\frac{\partial \alpha}{\partial x_1} \sin(\alpha) + \frac{\partial \alpha}{\partial x_2} \cos(\alpha). $$

It leads to

$$ \nabla r = \left( -\frac{\partial \alpha}{\partial x_2}, \frac{\partial \alpha}{\partial x_1} \right)^T . $$

Thus, again by Poincaré's lemma, the dilation factor $r$ exists if and only if the right hand side in the above is curl-free, which leads to the harmonic condition on $\alpha$. $\square$

The following proposition is a very useful property of conformal orientations.

**Proposition 6.3.** If there exists a map $\varphi = (\varphi_1, \varphi_2)$ from $D \subset \mathbb{R}^2$ to $\mathbb{R}^2$ such that

$$ \nabla \varphi = e^r Q(\alpha) \quad \text{in } D, $$

then all angles are preserved by the map $\varphi$. In particular, small square cells are deformed into almost square cells locally.

*Sketch of the proof.* Let $x = (x_1, x_2)$ be the origin of a small square $S$ of side length $\varepsilon > 0$ and edges $\varepsilon e_1, \varepsilon e_2$ (in other words, the vertices of $S$ are $x, x + \varepsilon e_1, x + \varepsilon e_2$ and $x + \varepsilon e_1 + \varepsilon e_2$). The map $\varphi$ then transforms $S$ into an “almost” square $\varphi(S)$ of vertices $\varphi(x), \varphi(x + \varepsilon e_1), \varphi(x + \varepsilon e_2)$ and $\varphi(x + \varepsilon e_1 + \varepsilon e_2)$. By a Taylor expansion, we see that $\varphi(S)$ is, up to terms of order $\varepsilon^2$, equal to a parallelogram of origin $\phi(x)$ and edges $\varepsilon \frac{\partial \varphi}{\partial x_1}$ and $\varepsilon \frac{\partial \varphi}{\partial x_2}$. Since $\nabla \varphi = e^r Q(\alpha)$, the two edges are orthogonal, so $\varphi(S)$ is “almost” a square with side length $\varepsilon e^r(x)$ (possibly rotated with respect to the initial square $S$). $\square$

Nevertheless, in the applications we face a problem. Since $\alpha$ is a stress eigen-direction, it has no reason to be a harmonic function in general. Even worse, $\alpha$ might not even be smooth at some places (for example, at corners or at the junction point of different boundary conditions, but at other places as well). A more profound reason lies in the following observation: both $\alpha$ and $\alpha + \pi$ give rise to the same orientation. By Pedersen's result, we can show that the rotated Hooke's law $A^*(m_1, m_2, \alpha) = R(\alpha)^T A^*(m_1, m_2, 0) R(\alpha)$ depends only on the double angle $2\alpha$.

From now on, we are going to be working with the double angle $\beta = 2\alpha$, thus, removing the indeterminate additive constant $\pi$. In what follows, we shall regularize the double angle $\beta = 2\alpha$ and make it harmonic.
<figure><img src="image_49.png" /><figcaption>Figure 40: A small square S and the distorted square $\varphi(S)$.</figcaption></figure>

### 6.6.4 Regularization of the double angle $\beta = 2\alpha$

As we mentioned before, the orientation $\alpha$ given by the optimization does not necessarily satisfy the conformality condition $\Delta\alpha = 0$. Thus, at each iteration of the algorithm, instead of minimizing locally (by using Pedersen's result) the following quantity

$$A^*(m_1, m_2, \beta)^{-1}\sigma \cdot \sigma,$$

let us consider to minimize globally

$$\int_D (A^*(m_1, m_2, \beta)^{-1}\sigma \cdot \sigma + \eta^2 |\nabla\beta|^2) \, dx$$

for a small parameter $\eta > 0$, under the harmonic constraint

$$\int_D \nabla\beta \cdot \nabla q \, dx = 0 \quad \forall q \in H_0^1(D).$$

This is a non-linear (and non-quadratic) constrained optimization problem. It turns out that working with the angle $\beta$ is not so easy since the Hooke's law $A^*(m_1, m_2, \beta)$ is highly nonlinear in terms of $\beta$. It is however quadratic with respect to the vector $b$, defined by

$$b = (\cos \beta, \sin \beta) \quad \text{and} \quad A^*(m_1, m_2, \beta) = S(b)^T A(m_1, m_2, 0) S(b),$$

where the matrix $S(b)$ is defined by $S(b) = R(\alpha)$. It can be shown that $S(b)$ is affine with respect to $b$ by a careful computation based on Pedersen's result.
Therefore, from now on we shall work with $b$ as the main unknown. We remark that $\nabla\beta = b \wedge \nabla b$. Indeed,

$$b \wedge \nabla b = (\cos \beta, \sin \beta) \wedge (-\sin \beta \nabla \beta, \cos \beta \nabla \beta) = \nabla \beta.$$

Therefore, the objective function becomes

$$\int_D (A^*(m_1, m_2, b)^{-1}\sigma \cdot \sigma + \eta^2 |b \wedge \nabla b|^2) \, dx, \qquad (6.10)$$

which in turn has to be minimized under the harmonic constraint

$$\int_D (b \wedge \nabla b) \cdot \nabla q \, dx = 0 \quad \forall q \in H_0^1(D). \qquad (6.11)$$

We iteratively solve the non-linear problem (6.10) under the minimization constraint (6.11) by a Newton-type approximation with an increment $\delta b$ as follows.

Find a step $\delta b^n \in H^1(D; \mathbb{R}^2)$ and a Lagrange multiplier $p^{n+1} \in H_0^1(D)$ such that, for any $\delta c \in H^1(D; \mathbb{R}^2)$ and $q \in H_0^1(D)$,

$$\int_D A^*(m)^{-1} S(b^n + \delta b^n) \sigma \cdot S'(\delta c) \sigma \, dx + \eta^2 \int_D (b^n \wedge \nabla(b^n + \delta b^n)) \cdot (b^n \wedge \nabla \delta c) \, dx \\ + \int_D (b^n \wedge \nabla \delta c) \cdot \nabla p^{n+1} \, dx = 0 \qquad (6.12)$$

and

$$\int_D (b^n \wedge \nabla(b^n + \delta b^n)) \cdot \nabla q \, dx = 0, \qquad (6.13)$$

where $S'(\delta c)$ is the directional derivative of $S(b)$ in the direction $\delta c$. Notice that, since $S$ is an affine function, then $S'(\delta c) = S(c + \delta c) - S(c)$. At each iteration, we update the vector field $b$ as follows

$$b^{n+1} = \frac{b^n + \delta b^n}{|b^n + \delta b^n|}. \qquad (6.14)$$

We now apply this regularization process together with the alternate minimization algorithm.

The above mentioned algorithm is structured as follows:
<figure><img src="image_50.png" /></figure>

**Figure 41:** Regularized orientation for the bridge case (left) and angle difference between optimized and regularized orientation.

**Algorithm 9 Regularization algorithm**

1. Initialization of the design parameters $m_1^0, m_2^0, b$ with the results of the optimization without the harmonic constraint.

2. Iterations until convergence, for $n \ge 0$:

(a) Computation of the stress tensor $\sigma^n$ through a problem of linear elasticity.

(b) Updating of the hole parameters $m^n$ by using the projected gradient algorithm (6.5) with the orientation $b^n$.

(c) Computation of the increment $\delta b^n$ by solving (6.12) and (6.13).

(d) Updating of the orientation with (6.14).

**Remark 6.4.** In numerical practice, starting from the optimal (but not necessarily smooth) orientation, a few tens of iterations of this regularization process are enough.

**Remark 6.5.** One advantage of this b-formulation is that it is insensitive to the $2\pi$-modulo of $\beta$.

As we can see by Figure 41, the regularization occurs mainly in areas where the density is close to 0 or 1, i.e. where the homogenized material is almost isotropic and the orientation has no significant impact.
<figure><img src="image_51.png" /><figcaption>Figure 42: Positive singularity (left) and negative singularity (right).</figcaption></figure>

**Remark 6.6.** Unfortunately this process does not always work in general. In other words, it does not always yield a smooth harmonic angle $\beta$. This is because of “true” singularities, i.e. singularities of the orientation that remain and thus do not allow the angle to be harmonic. The vector field is not coherently orientable in these cases, as the vector rotates by an angle of $\pm\pi$ along circles which are enclosing the singularities (see Figure 42. Such cases might be handled by means of other regularization processes (e.g. minimizing a Ginzburg-Landau energy [GE2018]).

### 6.6.5 Computation of the map $\varphi$

Once a harmonic angle $\alpha = \beta/2$ has been found, one needs to compute $r$ and $\varphi$ such that

$$ \nabla\varphi = e^r Q(\alpha) \quad \text{in } D. $$

Since the dilation field $r$ satisfies

$$ \nabla r = (-\nabla \wedge a_2)a_1 + (\nabla \wedge a_1)a_2 \quad \text{with } (a_1, a_2) = Q(\alpha), $$

one computes $r$ as the minimizer in $H^1(D)$ of

$$ \int_D |\nabla r + (\nabla \wedge a_2)a_1 - (\nabla \wedge a_1)a_2|^2 dx. $$

Once $r$ has been computed, a naive idea would be to compute $\varphi$ as a minimizer in $H^1(D; \mathbb{R}^2)$ of

$$ \int_D |\nabla\varphi - e^r Q(\alpha)|^2 dx. $$
However, we know that, even if $\beta$ is smooth, $\alpha$ may have jumps of the type $\pm\pi$ and thus $Q(\alpha)$ may have jumps of its sign (recall that $Q(\alpha + \pi) = -Q(\alpha)$).

To compute $\varphi$ there are two possibilities.

1. Find a coherent orientation of $\alpha$ (i.e. choose between $\alpha$ and $\alpha + \pi$ at every point): this is possible only if there are no singularities (this is the approach of Groen and Sigmund [GS2018]).

2. Leave the angle $\alpha$ as it is and extend $\varphi$ to be defined in an abstract manifold. This is the approach of Allaire-Geoffroy-Pantz [AGP2018] and it works also in the presence of singularities.

### 6.6.6 An abstract manifold setting

Let us introduce the cover space of $D$.

**Definition 6.7.** Denote by $T$ a rotation matrix field which is a candidate for being $Q(\alpha)$. Then we define

$$ \mathcal{D} = \{(x, T) \in D \times \mathrm{SO}(2) \text{ such that } T^2 = Q(\beta)\}, $$

where $\mathrm{SO}(2)$ is the set of rotations in $\mathbb{R}^2$.

We note that at every point $x \in D$ the rotation satisfies $T(x)^2 = Q(\beta)(x)$. If the angle $\alpha$ is globally orientable, then $T(x) = Q(\alpha)(x)$ or $T(x) = -Q(\alpha)(x)$, and that $\mathcal{D}$ is simply the union of two copies of $D$, consisting of the two possible signs of $Q(\alpha)$. We assume the simple case where $\alpha$ could be globally oriented (no singularity) but extend it to the singular cases.

We change our working space from $D$ to $\mathcal{D}$. The map $\varphi(x, T)$ is now defined on the manifold $\mathcal{D}$ by

$$ \nabla\varphi = e^r T, \qquad (6.15) $$

and the gradient operator in (6.15) defined by

$$ \nabla\varphi(x, T) = \nabla\varphi_U(x), $$

where $U$ is an orientable open subset of $D$,

$$ \varphi_U(x) = \varphi \circ g_U(x), $$
<figure><img src="image_52.png" /><figcaption>**Figure 43:** Two possible orientations (left) and the manifold obtained by gluing them together.</figcaption></figure>

and $g_U$ is one of the charts

$$ g_U^{\pm}: U \to \mathcal{D} \\ x \mapsto (x, \pm T_U(x)) $$

with $T_U^2 = Q(\beta)$ and $T_U \in C(U, \text{SO}(2))$. Moreover, without loss of generality, we can assume the antisymmetric property

$$ \varphi(x, -T) = -\varphi(x, T). $$

Indeed, if $\varphi$ satisfies $\nabla\varphi = e^r T$, then the map $(\varphi(x, T) - \varphi(x, -T))/2$ still satisfies (6.15) and is antisymmetric. Thus, if the orientation $\alpha$ satisfies the conformality condition (6.8), the map $\varphi$ can be defined as a minimizer of

$$ \min_{\varphi \in \mathcal{V}} \int_{\mathcal{D}} |\nabla \varphi - e^r T|^2 dx, $$

over all maps $\varphi$ in

$$ \mathcal{V} := \{ \varphi \in H^1(D, \mathbb{R}^2) : \varphi(x, -T) = -\varphi(x, T) \text{ for all } (x, T) \in \mathcal{D} \}. $$

In practice we face the problem of making the actual computations on the abstract manifold $\mathcal{D}$. In order to solve this problem we use a new idea, namely, non-conformal finite elements on $\mathcal{D}$.
<figure><img src="image_53.png" /><figcaption>**Figure 44:** Left: orientation of $\beta$ (black arrows) and of $\alpha$ (red arrows). Right: two possible coherent orientations of $\alpha$.</figcaption></figure>

On each triangle $K$ of the mesh we compute one continuous orientation $T_K$ such that $T_K^2 = Q(\beta)$. Then, we glue together (with $P_1$ discontinuous finite elements) these orientations. We compute

$$ \int_D |\nabla \varphi - e^r T|^2 dx = \sum_K \int_{g_K^+(K) \cup g_K^-(K)} |\nabla \varphi - e^r T|^2 dx \\ = \sum_K \int_K |\nabla (\varphi \circ g_K^+) - e^r T_K(x)|^2 dx \\ + \sum_K \int_K |\nabla (\varphi \circ g_K^-) - e^r T_K(x)|^2 dx $$

with $g_K^\pm = \text{Id} \times (\pm T_K)$. By the antisymmetry of $\varphi$, we obtain

$$ \int_D |\nabla \varphi - e^r T|^2 dx = 2 \sum_K \int_K |\nabla (\varphi \circ g_K^+) - e^{r(x)} T_K(x)|^2 dx. $$

Then, we minimize with respect to $\varphi$ in the space of $P_1$ discontinuous finite elements.

## 6.7 A final post-processing/cleaning of the lattice reconstruction

The shapes we obtained in Section 6.6 are not straightforwardly manufacturable. Indeed, there are disconnected components of the lattice structure and/or too thin members that should be removed. A final post-processing is made to cure these defects. Note that there is room for improvement in the process.
<figure><img src="image_54.png" /><figcaption>**Figure 45:** $\Omega_{\varepsilon}(\varphi, m)$ for several $\varepsilon$ in the case of the bridge.</figcaption></figure>

<figure><img src="image_55.png" /><figcaption>**Figure 46:** The map $|\varphi_i|$ (isolines) and the vectors $a_i$ (arrows) for $i = 1$ (left) and $i = 2$ (middle). On the right we have the projection of a regular grid by the map $\varphi$.</figcaption></figure>
Let $h_{\min}$ be the minimal manufacturable lengthscale or feature size, meaning the smallest possible width of bars and diameter of holes which can be effectively built. Recall that $\varepsilon$ is our choice of a global size of cells. After deformation, the cell size is $h_c(x) = \varepsilon e^{-r(x)}$. Hence the local widths of the bars and holes are respectively given by $(1 - m_i(x))h_c(x)$ and $m_i(x)h_c(x)$.

In the following, we distinguish two regimes, depending of the local size of the cell $h_c(x)$. First, if the cell size is too small, a hole and a bar of minimal width cannot coexist and then we have to choose a completely full or void cell. Hence, if $h_c < 2h_{\min}$, a thresholding is applied separately to each field $m_i$: it is assigned the value 0 if $m_i < 0.5$ and 1 otherwise.

Second, when $h_c \ge 2h_{\min}$, our post-processing criterion is satisfied if

$$ \frac{h_{\min}}{h_c} \le m_i \le 1 - \frac{h_{\min}}{h_c}, \quad i = 1, 2. $$

Otherwise, we simply threshold the values of $m_1$ and $m_2$, according to Figure 47, in order to reach void or full materials. The thresholded $m$ is then denoted by $\tilde{m}$.

<figure><img src="image_56.png" /><figcaption>Figure 47: Thresholding criteria.</figcaption></figure>

Let $O_{\varepsilon}(\varphi, \tilde{m})$ be the shape obtained from $\Omega_{\varepsilon}(\varphi, \tilde{m})$ by filling its closed holes.
Numerically, the complement of $\Omega_{\varepsilon}(\varphi, \tilde{m})$ is computed step by step, by evaluating the sign of $F_{\varepsilon}^{\varphi, \tilde{m}}$. If it is positive, the current vertex belongs to the complement $\Omega_{\varepsilon}^{c}(\varphi, \tilde{m})$ and then its neighbors, which are not already visited, are added to the list of vertices which should be tested. Otherwise, the current vertex does not belong to $\Omega_{\varepsilon}^{c}(\varphi, \tilde{m})$.

We will regularize the subset $O_{\varepsilon}(\varphi, \tilde{m})$ in order to remove the disconnected bars or the bars that have one free endpoint. Numerically, we explore all the vertices of the complement as follows. For any given vertex, we check each other vertex not further away than a distance $h_{\min}$: if this vertex belongs to the complement too, all vertices between them are added to the complement. In this way, we suppress all disconnected bars and all bars of $O_{\varepsilon}(\varphi, \tilde{m})$ that have one free endpoint, which are not too wide. This new subset is denoted by $\tilde{O}_{\varepsilon}(\varphi, \tilde{m})$.

Finally, the post-processed structure is given by the intersection $\tilde{\Omega}_{\varepsilon}(\varphi, \tilde{m}) := \Omega_{\varepsilon}(\varphi, \tilde{m}) \cap \tilde{O}_{\varepsilon}(\varphi, \tilde{m})$. Several post-processed structures $\tilde{\Omega}_{\varepsilon}(\varphi, \tilde{m})$ for the bridge case are displayed in Figure 48.

## 6.8 Other numerical examples

In this section, we show some numerical examples of the application of the whole method of this chapter. Here, we mention that the method is not always applicable because some singularities cannot be eliminated, as we will see later.

We will show the numerical results for the optimization of a cantilever (Figure 49), an MBB beam (Figure 50) and an L-beam (Figure 51). For each case, we have represented:

(a) the optimal orientation of the periodicity cells before regularization,

(b) the optimal orientation of the periodicity cells after regularization,

(c) the underlying lattice on which the optimal composite is built, i.e., the projection by $\varphi$,

(d)-(f) the sequence of shapes after post-processing for the case of $\varepsilon = 0.2, 0.1$ and $0.05$ respectively.
<figure><img src="image_57.png" /><figcaption>**Figure 48:** Post-processed $\bar{\Omega}_{\varepsilon}(\varphi, \bar{m})$ for several $\varepsilon$ in the case of the bridge.</figcaption></figure>

As shown in the case of L-beam (Figure 51), the singularities which appear in Figure 51-(a) are removed during the regularization step (Figure 51-(b)). This is a necessary condition in order to apply our method. Indeed, as we mentioned, there is a case where the singularity cannot be removed, the so-called electrical mast (see Figure 52). Figure 52-(a) shows that two negative singularities, located inside the domain, cannot neither be removed nor pushed toward the boundary during the regularization step.

Hence, the computed grid shown in Figure 52-(b) is clearly not correctly aligned with the optimal orientation of the cells in the vicinity of the singularities. To overcome this problem, at least two different strategies can be considered. One consists in modifying the regularization process in a way that forces more effectively the singularities to be eliminated (see Figure 52-(c) and (d)). Another approach is to adapt the projection step so that it is able to take singularities into account, and
<figure><img src="image_58.png" /><figcaption>(a) Opt. orientation</figcaption></figure>

<figure><img src="image_59.png" /><figcaption>(b) Regul. orientation</figcaption></figure>

<figure><img src="image_60.png" /><figcaption>(c) Projection by $\varphi$</figcaption></figure>

<figure><img src="image_61.png" /><figcaption>(d) $\bar{\Omega}_{\varepsilon}(\varphi, \bar{m}), \varepsilon = 0.2$</figcaption></figure>

<figure><img src="image_62.png" /><figcaption>(e) $\bar{\Omega}_{\varepsilon}(\varphi, \bar{m}), \varepsilon = 0.1$</figcaption></figure>

<figure><img src="image_63.png" /><figcaption>(f) $\bar{\Omega}_{\varepsilon}(\varphi, \bar{m}), \varepsilon = 0.05$</figcaption></figure>

<caption>Figure 49: Cantilever case.</caption>

will be discussed in the next section.

## 6.9 Further issues

As we mentioned in Section 6.6, the problem of removing the singularity might persist. To overcome this problem, some preliminary remedies are proposed in the PhD thesis of Perle Geoffroy [GE2018] by either trying to eliminate them by a Ginzburg-Landau approach or compute a map $\varphi$ with the previous approach and an enriched discontinuous finite element space. In [GE2018], one may find others extensions to cases, such as different objective functions, multiple loads and 3-d problems.
<figure><img src="image_64.png" /><figcaption>(a) Opt. orientation</figcaption></figure>

<figure><img src="image_65.png" /><figcaption>(b) Regul. orientation</figcaption></figure>

<figure><img src="image_66.png" /><figcaption>(c) Projection by $\varphi$</figcaption></figure>

<figure><img src="image_67.png" /><figcaption>(d) $\bar{\Omega}_{\varepsilon}(\varphi, \bar{m}), \varepsilon = 0.2$</figcaption></figure>

<figure><img src="image_68.png" /><figcaption>(e) $\bar{\Omega}_{\varepsilon}(\varphi, \bar{m}), \varepsilon = 0.1$</figcaption></figure>

<figure><img src="image_69.png" /><figcaption>(f) $\bar{\Omega}_{\varepsilon}(\varphi, \bar{m}), \varepsilon = 0.05$</figcaption></figure>

<caption>Figure 50: MBB beam case.</caption>

## 6.10 Exercises

**Problem 6.10.1.** Consider a compliance minimization problem (for one test case like cantilever, bridge, MBB beam, etc.) for the homogenized formulation. Choose an homogenized tensor corresponding to a non-isotropic microstructure (for example a square cell with a rectangular hole with fixed size). Then minimize the compliance with respect to the sole orientation of the microstructure (using Pedersen result).

**Problem 6.10.2.** For the same compliance minimization problem as in the previous exercise, fix now the orientation and let the parameters of the microstructure (for example, the lengths $m_1$ and $m_2$ of the rectangular hole, see Figure 22) become the optimization variables. Then minimize the compliance with respect to these parameters (with fixed orientation). Compare the optimal designs and the attained minimal compliances with the previous exercise. In particular, notice that the absence of orientation optimization yields a self-penalizing effect, namely the obtained designs feature almost no intermediate densities (somehow similarly to the SIMP method).

**Problem 6.10.3.** Combine the two previous optimization (with respect to the orientation and the size parameters) and minimize the compliance for the test case of the previous exercises.
<figure><img src="image_70.png" /><figcaption>(a) Opt. orientation</figcaption></figure>

<figure><img src="image_71.png" /><figcaption>(b) Regul. orientation</figcaption></figure>

<figure><img src="image_72.png" /><figcaption>(c) Projection by $\varphi$</figcaption></figure>

<figure><img src="image_73.png" /><figcaption>(d) $\bar{\Omega}_{\varepsilon}(\varphi, \bar{m}), \varepsilon = 0.2$</figcaption></figure>

<figure><img src="image_74.png" /><figcaption>(e) $\bar{\Omega}_{\varepsilon}(\varphi, \bar{m}), \varepsilon = 0.1$</figcaption></figure>

<figure><img src="image_75.png" /><figcaption>(f) $\bar{\Omega}_{\varepsilon}(\varphi, \bar{m}), \varepsilon = 0.05$</figcaption></figure>

<caption>Figure 51: L-beam case.</caption>
<figure><img src="image_76.png" /><figcaption>(a) Singularities appear</figcaption></figure>

<figure><img src="image_77.png" /><figcaption>(c) Improved reg. orientation</figcaption></figure>

<figure><img src="image_78.png" /><figcaption>(b) Projection by $\varphi$</figcaption></figure>

<figure><img src="image_79.png" /><figcaption>(d) Reconstructed grid</figcaption></figure>

**Figure 52:** Improved regularized orientation in the case of the electrical mast.
Bibliography

[Al1992] G. ALLAIRE, *Homogenization and two-scale convergence*, SIAM J. Math. Anal. **23** (1992), 1482-1518.

[Al2002] G. ALLAIRE, Shape optimization by the homogenization method, vol **146** of Applied Mathematical Sciences, Springer-Verlag, New York, (2002).

[Al2007-1] G. ALLAIRE, Conception optimale de structures, vol **58** of Mathématiques et Applications, Springer, Heidelberg, (2007).

[Al2007-2] G. ALLAIRE, Numerical Analysis and Optimization. An introduction to mathematical modelling and numerical simulation. Translated from the French by Alan Craig, Numer. Math. Sci. Comput., Oxford University Press, Oxford, UK, (2007).

[AGP2018] G. ALLAIRE, P. GEOFFROY-DONDERS AND O. PANTZ, *Topology optimization of modulated and oriented periodic microstructures by the homogenization method*, Computers & Mathematics with Applications, special issue SimAM, (2019).

[AP2006] G. ALLAIRE AND O. PANTZ, *Structural optimization with FreeFem++*, Struct. Multidiscip. Optim., **32** (2006), 173-181.

[BK1988] M. BENDSØE AND N. KIKUCHI, *Generating optimal topologies in structural design using a homogenization method*, Comput. Methods Appl. Mech. Engrg. **71** (1988), no. 2, 197-224.

[BS2003] M. BENDSØE AND O. SIGMUND, Topology optimization, Springer-Verlag, Berlin, (2003).
[BLP1978] A. BENSOUSSAN, J.-L. LIONS AND G. PAPANICOLAOU, Asymptotic analysis for periodic structures, studies in mathematics and its applications, North-Holland, Amsterdam, (1978)

[BGLS2006] J.F. BONNANS, C. GILBERT, C. LEMARÉCHAL AND C. SAGAS-TIZÁBAL, Numerical optimization. Theoretical and practical aspects. Second edition. Universitext. Springer-Verlag, Berlin, 2006.

[CH2000] A. CHERKAEV, Variational Methods for Structural Optimization, Springer Verlag, New York (2000).

[CD1999] D. CIORANESCU AND P. DONATO, An introduction to homogenization, Oxford Lecture Series in Mathematics and Applications, **17**, Oxford (1999).

[ET1999] I. EKELAND AND R. TÉMAM, Convex analysis and variational problems, Classics in Applied Mathematics, **28**, Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA (1999).

[GLRS2013] B. GEIHE, M. LENZ, M. RUMPF AND R. SCHULTZ, *Risk averse elastic shape optimization with parametrized fine scale geometry*, Math. Program., Ser. A., **141** (2013), 383–403.

[GE2018] P. GEOFFROY-DONDERS, *Homogenization method for topology optimization of structures built with lattice materials*, PhD thesis, Ecole Polytechnique, Université Paris-Saclay (2018).

[GAP2018] P. GEOFFROY-DONDERS, G. ALLAIRE AND O. PANTZ, *3-d topology optimization of modulated and oriented periodic microstructures by the homogenization method*, submitted. HAL preprint: hal-01939201 (November 2018).

[GRS2015] I. GIBSON, D. ROSEN, B. STUCKER, *Additive Manufacturing Technologies*, Springer, New York (2015).

[GS2018] J.P. GROEN AND O. SIGMUND, *Homogenization based topology optimization for high resolution manufacturable microstructures*, International Journal for Numerical Methods in Engineering, **113** (2018), 1148–1163.
[HS1963] Z. HASHIN AND S. SHTRIKMAN, *A variational approach to the theory of the elastic behavior of multiphase materials*, J. Mech. Phys. Solids, **11** (1963), 127–140.

[HM2003] J. HASLINGER AND R. MÄKINEN, *Introduction to shape optimization. Theory, approximation, and computation*, SIAM, Philadelphia (2003).

[He2012] F. HECHT, *New development in freefem++*, J. Numer. Math., **20** (2012), 251–265.

[HP2018] A. HENROT AND M. PIERRE, Shape variation and optimization. A geometrical analysis. English version of the French publication with additions and updates. EMS Tracts in Mathematics, 28. European Mathematical Society (EMS), Zrich, 2018.

[H1996] U. HORNUNG, Editor, *Homogenization and porous media*, Springer Verlag (1996).

[JKO1995] V. JIKOV, S. KOZLOV AND O. OLEINIK, *Homogenization of differential operators and integral functionals*, Springer, Berlin, (1995).

[KPTZ2000] B. KAWOHL, O. PIRONNEAU, L. TARTAR AND J.-P. ZOLSIO, Optimal shape design. Lectures given at the Joint C.I.M./C.I.M.E. Summer School held in Tria, June 16, 1998. Edited by A. Cellina and A. Ornelas. Lecture Notes in Mathematics, 1740. Fondazione CIME/CIME Foundation Subseries. Springer-Verlag, Berlin; Centro Internazionale Matematico Estivo (C.I.M.E.), Florence, (2000).

[MI2001] G. MILTON, *The theory of composites*, Cambridge University Press (2001).

[Mu1977] F. MURAT, *Contre-exemples pour divers problèmes où le contrôle intervient dans les coefficients*, Annali Mat. Pura Appl., **112** (1977), 49–68.

[MT1997] F. MURAT AND L. TARTAR, *H-Convergence*, in Topics in the Mathematical Modeling of Compositic Materials, Progress in Nonlinear Differential
Equations and their Applications., ed. A. Cherkaev and R. Kohn, Birkhäuser,
Boston, **31** (1997), 21–43.

[Ng1989] G. NGUETSENG, *A general convergence result for a functional related to the theory of homogenization*, SIAM J.Math. Anal. **20** (1989), 608–623.

[NW1999] J. NOCEDAL AND S. WRIGHT, Numerical optimization, Springer Science
(1999).

[PT2008] O. PANTZ AND K. TRABELSI, *A post-treatment of the homogenization method for shape optimization*, SIAM J. Control Optim., **47** (2008), 1380–1398.

[Pe1989] P. PEDERSEN, *On optimal orientation of orthotropic materials*, Structural optimization **1** (2) (1989), 101–106.

[ROZ1989] G. ROZVANY, *Structural design via optimality criteria*, Kluwer Academic Publishers, Dordrecht (1989).

[SK1992] J. SOKOLOWSKI, J.-P. ZOLÉSIO, *Introduction to shape optimization. Shape sensitivity analysis*, Springer Series in Computational Mathematics, 16, Springer, Berlin (1992).

[TA2000] L. TARTAR, *An introduction to the homogenization method in optimal design*, in Optimal shape design (Tróia, 1998), A. Cellina and A. Ornelas eds., Lecture Notes in Mathematics 1740, pp.47-156, Springer, Berlin (2000).