---
jupyter:
  description: 'An auxetic material designed by level-set shape optimization with planex2d and NGSolve: a periodic unit
    cell whose homogenized Poisson''s ratio is driven below zero under a volume and a stiffness constraint, its shape
    derivative checked by finite differences, and the result turned into a clean planex2d cell, checked on its own
    mesh and tiled into a lattice.'
  nbsphinx:
    timeout: 1200
---

# Designing an auxetic material

Can a material get wider when it is pulled? Most materials get narrower, with a Poisson's ratio between 0 and 0.5, but a
lattice of the right shape does the opposite. Its cells unfold as they are stretched, and its Poisson's ratio is
negative. Such auxetic materials are studied for stents, protective padding and morphing structures. In two dimensions,
the Poisson's ratio of a stable isotropic material cannot fall below −1.

The example designs such a lattice. Its unit cell is filled with material where a level set function is positive, and
its effective stiffness comes from homogenization. The shape of the material is moved step by step to lower the
effective Poisson's ratio, while constraints hold the amount of material and the stiffness. The shape derivative is
checked against finite differences, and the result becomes a clean planex2d cell, whose Poisson's ratio is checked on
its own mesh.

This page is part of a collection of [field and simulation examples](fields-and-simulation.md) using planex2d.

![A lattice of auxetic cells evolving from square holes to rotating petals during the optimization](../_static/gallery/auxetic-cell.gif){.result}

```python
import matplotlib.pyplot as plt
import numpy as np
import planex as px
from IPython.display import HTML
from matplotlib.animation import FuncAnimation
from ngsolve import *
```

```python nbsphinx="hidden"
from docs_figures import FACE, INK, PLANEX, Chart, Row
```

## The periodic cell

The cell is the unit square with its opposite sides identified, so that whatever the material does at one side, it
continues at the other, and the cell stands for an infinite lattice. It is meshed with triangles of side 0.02:

```python
h = 0.02
cell = px.QuickSketch().Rectangle(1, 1).Face()
for edge in cell.edges:
    (x0, y0), (x1, y1) = edge.bounding_box
    edge.name = "left" if x1 < 1e-9 else "right" if x0 > 1 - 1e-9 else "bottom" if y1 < 1e-9 else "top"
cell.Identify("left", "right")
cell.Identify("bottom", "top")
mesh = Mesh(px.OCCGeometry(cell).GenerateMesh(maxh=h))
print(mesh.ne, "triangles")
```

## The effective stiffness

The material fills the cell where the level set function is positive. Each triangle gets the share $\chi$ of its area
that lies there, and elsewhere a weak ersatz material stands in for the void, with the stiffness factor
$m = \varepsilon + (1 - \varepsilon)\chi$. For the two unit strains $E^1 = e_x \otimes e_x$ and $E^2 = e_y \otimes e_y$,
the periodic fluctuations $w^\alpha$ solve

$$\int_Y m \, C\big(E^\alpha + \epsilon(w^\alpha)\big) : \epsilon(v) \, dy = 0,$$

and with $A^\alpha = E^\alpha + \epsilon(w^\alpha)$ the effective, homogenized tensor and Poisson's ratio are

$$C^H_{\alpha\beta} = \int_Y m \, C A^\alpha : A^\beta \, dy, \qquad \nu^H = \frac{C^H_{12}}{\sqrt{C^H_{11} C^H_{22}}}.$$

The base material has a Young's modulus of 1 and a Poisson's ratio of 0.3, in plane strain. One factorization serves
both unit strains, and a tiny mass term removes the free translations:

```python
E, nu0, eps = 1, 0.3, 1e-3
mu, lam = E / (2 * (1 + nu0)), E * nu0 / ((1 + nu0) * (1 - 2 * nu0))
V = Periodic(VectorH1(mesh, order=2))
u, v = V.TnT()
chi = GridFunction(L2(mesh, order=0))
m = eps + (1 - eps) * chi


def strain(w):
    return Sym(Grad(w))


def C(A):
    return 2 * mu * A + lam * Trace(A) * Id(2)


unit = [CF((1, 0, 0, 0), dims=(2, 2)), CF((0, 0, 0, 1), dims=(2, 2))]
a = BilinearForm(m * InnerProduct(C(strain(u)), strain(v)) * dx + 1e-8 * InnerProduct(u, v) * dx)
loads = [LinearForm(-m * InnerProduct(C(Ea), strain(v)) * dx) for Ea in unit]
fluctuations = [GridFunction(V) for Ea in unit]
inverse = a.Assemble().mat.Inverse(V.FreeDofs(), inverse="sparsecholesky")
A = [Ea + strain(w) for Ea, w in zip(unit, fluctuations)]


def solve_cells():
    a.Assemble()
    inverse.Update()
    for f, w in zip(loads, fluctuations):
        f.Assemble()
        w.vec.data = inverse * f.vec


def CH(i, j):
    return Integrate(m * InnerProduct(C(A[i]), A[j]), mesh)
```

## The level set

The level set function $\varphi$ lives in a discontinuous space of order 1, in which it is moved, and is copied into a
continuous one, from which the share of each triangle is computed. For a linear function on a triangle, the region where
it is positive is the whole triangle, nothing, a corner cut off by a straight line, or the triangle less such a corner.
The corner at the lone vertex $\varphi_0$ has the share
$\varphi_0^2 / \big((\varphi_0 - \varphi_1)(\varphi_0 - \varphi_2)\big)$. The start is a lattice of rounded square
holes:

```python
phi = GridFunction(Periodic(L2(mesh, order=1)))
phi1 = GridFunction(H1(mesh, order=1))
convert = ConvertOperator(L2(mesh, order=1), H1(mesh, order=1))
phi.Set(0.1 - cos(2 * pi * x) * cos(2 * pi * y) + 0.3 * (cos(2 * pi * x) + cos(2 * pi * y)))
triangles = np.array([[p.nr for p in el.vertices] for el in mesh.Elements(VOL)])
points = np.array([p.point for p in mesh.vertices])


def shares(values):
    """the share of each triangle where the linear function with these vertex values is positive"""
    t = values[triangles]
    positive = t > 0
    count = positive.sum(axis=1)
    out = np.where(count == 3, 1.0, 0.0)
    for k in (1, 2):
        rows = np.flatnonzero(count == k)
        lone = np.argmax(positive[rows] if k == 1 else ~positive[rows], axis=1)
        tip, b, c = (t[rows, (lone + j) % 3] for j in (0, 1, 2))
        corner = tip**2 / ((tip - b) * (tip - c))
        out[rows] = corner if k == 1 else 1 - corner
    return out


def material():
    phi1.vec.data = convert * phi.vec
    chi.vec.FV().NumPy()[:] = shares(phi1.vec.FV().NumPy())


material()
print("volume fraction at the start", Integrate(chi, mesh))
```

## The shape derivative

$C^H_{\alpha\beta}$ is stationary in the fluctuations, because the unit strains are part of its integrand, so moving
the boundary by a velocity field $\psi$ changes it without an adjoint problem:

$$dC_{\alpha\beta}(\psi) = \int_Y m \Big[ (C A^\alpha : A^\beta) \, \mathrm{div}\,\psi
- C A^\alpha : \mathrm{sym}(\nabla w^\beta \nabla\psi) - C A^\beta : \mathrm{sym}(\nabla w^\alpha \nabla\psi) \Big] dy,
\qquad dV(\psi) = \int_Y \chi \, \mathrm{div}\,\psi \, dy.$$

The derivative of $\nu^H$ follows by the quotient rule. The constraints are a volume fraction of 0.4 and an average
stiffness $(C^H_{11} + C^H_{22})/2$ of at least 0.15 of the base material's $\lambda + 2\mu$. Without the stiffness
constraint, the optimizer would thin every hinge into a mechanism. The descent field $X$ represents the combined
derivative in the $H^1$ inner product. It is periodic and scaled to a largest value of 1:

```python
VX = Periodic(VectorH1(mesh, order=2))
X, psi = VX.TnT()
full = lam + 2 * mu


def dC(i, j):
    f = LinearForm(VX)
    f += m * (InnerProduct(C(A[i]), A[j]) * div(psi) - InnerProduct(C(A[i]), Sym(Grad(fluctuations[j]) * Grad(psi)))
              - InnerProduct(C(A[j]), Sym(Grad(fluctuations[i]) * Grad(psi)))) * dx
    return f


dC11, dC22, dC12 = dC(0, 0), dC(1, 1), dC(0, 1)
dV = LinearForm(chi * div(psi) * dx)
riesz = BilinearForm((InnerProduct(Grad(X), Grad(psi)) + InnerProduct(X, psi)) * dx).Assemble()
riesz_inverse = riesz.mat.Inverse(VX.FreeDofs(), inverse="sparsecholesky")
field, rhs = GridFunction(VX), riesz.mat.CreateColVector()


def descent(c11, c22, c12, ell_volume, ell_stiffness):
    """the descent field for nu_H with the two multipliers, scaled to 1; returns its size before scaling"""
    for f in (dC11, dC22, dC12, dV):
        f.Assemble()
    s = sqrt(c11 * c22)
    w = ell_stiffness / (2 * full)
    rhs.data = (1 / s) * dC12.vec - (c12 / (2 * s * c11) + w) * dC11.vec - (c12 / (2 * s * c22) + w) * dC22.vec
    rhs.data += ell_volume * dV.vec
    field.vec.data = riesz_inverse * rhs
    size = max(abs(field.vec.FV().NumPy()))
    field.vec[:] *= 1 / size
    return size
```

## Moving the boundary

The level set is carried along the field by the transport equation $\partial_t \varphi + X \cdot \nabla\varphi = 0$,
discretized with upwind fluxes between the elements and three explicit Runge-Kutta stages per step. Every ten
iterations it is reset to the signed distance from its zero line, which keeps its gradient near 1. The distance is
taken to the zero line's segments in the cell and its eight neighbors, so that it is periodic too:

```python
p, q = phi.space.TnT()
n = specialcf.normal(2)
wind = GridFunction(VX)
transport = BilinearForm(phi.space, nonassemble=True)
transport += -wind * p * grad(q) * dx
transport += wind * n * IfPos(wind * n, p, p.Other()) * q * dx(element_boundary=True)
mass_inverse = L2(mesh, order=1).Mass(1).Inverse()
tau, steps = h / 6, 3


def move():
    """three steps of the transport equation along wind, each with the three stages of SSP-RK3"""
    k1, k2, k3, w, s = (phi.vec.CreateVector() for _ in range(5))
    for _ in range(steps):
        transport.Apply(phi.vec, k1)
        w.data = phi.vec - tau * mass_inverse * k1
        transport.Apply(w, k2)
        s.data = k1 + k2
        w.data = phi.vec - tau / 4 * mass_inverse * s
        transport.Apply(w, k3)
        s.data = 1 / 6 * k1 + 1 / 6 * k2 + 2 / 3 * k3
        phi.vec.data -= tau * mass_inverse * s


def zero_line():
    """the segment where the continuous level set is zero in each cut triangle, as its two ends"""
    t, corners = phi1.vec.FV().NumPy()[triangles], points[triangles]
    cut = (t.max(axis=1) > 0) & (t.min(axis=1) < 0)
    t, corners = t[cut], corners[cut]
    edges = ((0, 1), (1, 2), (2, 0))
    with np.errstate(all="ignore"):
        ends = np.stack([corners[:, i] + (t[:, i] / (t[:, i] - t[:, j]))[:, None] * (corners[:, j] - corners[:, i])
                         for i, j in edges], axis=1)
    crossed = np.stack([t[:, i] * t[:, j] < 0 for i, j in edges], axis=1)
    order, rows = np.argsort(~crossed, axis=1, kind="stable"), np.arange(len(t))
    return ends[rows, order[:, 0]], ends[rows, order[:, 1]]


def reinitialize():
    """the level set reset to the signed distance from its zero line, periodic"""
    heads, tails = zero_line()
    shifts = [(i, j) for i in (-1, 0, 1) for j in (-1, 0, 1)]
    heads, tails = np.concatenate([heads + s for s in shifts]), np.concatenate([tails + s for s in shifts])
    d = np.empty(len(points))
    for i in range(0, len(points), 500):
        pts, ab = points[i:i + 500, None, :], (tails - heads)[None]
        s = np.clip(np.sum((pts - heads[None]) * ab, 2) / np.maximum(np.sum(ab * ab, 2), 1e-30), 0, 1)
        d[i:i + 500] = np.sqrt(np.sum((pts - heads[None] - s[..., None] * ab) ** 2, 2)).min(axis=1)
    values = phi1.vec.FV().NumPy()
    values[:] = np.where(values >= 0, 1, -1) * d
    phi.Set(phi1)


reinitialize()
material()
print("volume fraction after the reset", Integrate(chi, mesh))
```

## The derivative check

Moving the material by $\pm t\psi$ for a smooth periodic test field and taking central differences of $C^H$ and the
volume checks the shape derivative. They agree to within a few percent, which is what the finite step of the
transport allows:

```python
solve_cells()
c11, c22, c12 = CH(0, 0), CH(1, 1), CH(0, 1)
for f in (dC11, dC22, dC12):
    f.Assemble()
dV_check = LinearForm(chi * div(psi) * dx).Assemble()
test = GridFunction(VX)
test.Set(0.5 * CF((sin(2 * pi * x) * cos(2 * pi * y), cos(2 * pi * x) * sin(2 * pi * y))))
start = phi.vec.CreateVector()
start.data = phi.vec
moved = []
for sign in (1, -1):
    phi.vec.data = start
    wind.vec.data = sign * test.vec
    move()
    material()
    solve_cells()
    moved.append((CH(0, 0), CH(1, 1), CH(0, 1), Integrate(chi, mesh)))
t = steps * tau
for name, f, k in (("C11", dC11, 0), ("C22", dC22, 1), ("C12", dC12, 2), ("V", dV_check, 3)):
    print(name, "| shape derivative", InnerProduct(f.vec, test.vec),
          "| finite difference", (moved[0][k] - moved[1][k]) / (2 * t))
phi.vec.data = start
material()
```

## The optimization

Each iteration solves the two cell problems, computes $\nu^H$, the volume and the stiffness, updates the two multipliers
like a PI controller, and moves the boundary along the descent field. Two hundred iterations take the Poisson's ratio
from about 0.44 to about −0.77, and the volume to its target. The stiffness constraint is not met within them, as the
cell ends at about 0.04 of $\lambda + 2\mu$, but its multiplier grows as long as the constraint is violated and holds
the hinges against thinning further (Figure 1):

```python
iterations, target, stiffness = 200, 0.4, 0.15
cP, cI, clip = 50, 0.3, 0.02
ell, ell_s = 0.0, 0.0
history, frames = {"nu": [], "volume": [], "stiffness": []}, []
with TaskManager():
    for it in range(iterations):
        material()
        if it % 10 == 0:
            reinitialize()
            material()
        solve_cells()
        c11, c22, c12, volume = CH(0, 0), CH(1, 1), CH(0, 1), Integrate(chi, mesh)
        nu_h = c12 / sqrt(c11 * c22)
        e, g = volume - target, (stiffness * full - (c11 + c22) / 2) / full
        history["nu"].append(nu_h)
        history["volume"].append(volume)
        history["stiffness"].append((c11 + c22) / 2 / full)
        if it % 2 == 0:
            frames.append(phi1.vec.FV().NumPy().copy())
        descent(c11, c22, c12, ell + cP * max(-clip, min(clip, e)), max(0, ell_s + cP * max(-clip, min(clip, g))))
        wind.vec.data = -field.vec
        move()
        ell += cI * e
        ell_s = max(0, ell_s + cI * g)
        if it % 50 == 0 or it == iterations - 1:
            print("iteration", it, "| nu_H", round(nu_h, 4), "| volume", round(volume, 4),
                  "| stiffness", round(history["stiffness"][-1], 4))
material()
frames.append(phi1.vec.FV().NumPy().copy())
solve_cells()
nu_h = CH(0, 1) / sqrt(CH(0, 0) * CH(1, 1))
print("C^H", np.round([[CH(0, 0), CH(0, 1)], [CH(0, 1), CH(1, 1)]], 4), "| nu_H", nu_h)
# plt.plot(history["nu"]); plt.show()
```

```python nbsphinx="hide-input"
Chart("Iteration", "Value", curves={"Poisson's ratio ν_H": (range(iterations), history["nu"])},
      lines={"Volume fraction": (range(iterations), history["volume"])},
      title="The cell's Poisson's ratio and volume over the iterations")
```

Six cells of the lattice show how the square holes turn into rotating petals, joined by thin hinges:

```python nbsphinx="hide-input"
fig, ax = plt.subplots(figsize=(6, 4))


def frame(k):
    ax.clear()
    for i in range(3):
        for j in range(2):
            ax.tricontourf(points[:, 0] + i, points[:, 1] + j, triangles, frames[min(k, len(frames) - 1)],
                           levels=[0, 1e9], colors=[PLANEX])
    ax.set_xlim(0, 3)
    ax.set_ylim(0, 2)
    ax.set_aspect("equal")
    ax.axis("off")


animation = FuncAnimation(fig, frame, frames=len(frames) + 20, interval=60)
plt.close(fig)
HTML(animation.to_jshtml())
```

## A clean cell

`px.FromLevelSet` turns the region where the continuous level set is positive into a planex2d shape. The mesh is
periodic and its sides are named, so the crossings on opposite sides match and the shape keeps the identifications. Its
splines are fitted within 1e-2, and edges shorter than that are removed, as on the page
[From level sets](../integrations/netgen/levelset.md). The pieces on opposite sides stay matched, so the cell can be
meshed with its identifications:

```python
traced = px.FromLevelSet(phi1, edge_length=1e-3)
fit = dict(max_edge_length=0.1, tol=1e-2, ignore_edges_smaller=1e-3)
clean = traced.Simplified(remove_edges_smaller=1e-2, fit_splines=fit)
print("area", clean.mass, "| volume fraction of the optimization", Integrate(chi, mesh), "| faces", len(clean.faces))
print("identifications", list(clean.identifications))
# px.Plot(clean)
```

On its own mesh, without ersatz material, the clean cell is softer than the optimization says and more auxetic, with a
Poisson's ratio of about −0.93 against −0.77. The difference lies in the thin hinges, which the ersatz material and the
cut triangles of the optimization mesh smear over whole triangles and make stiffer:

```python
fine = Mesh(px.OCCGeometry(clean).GenerateMesh(maxh=h / 2))
fine.Curve(3)
W = Periodic(VectorH1(fine, order=2))
u, v = W.TnT()
k = BilinearForm(InnerProduct(C(strain(u)), strain(v)) * dx + 1e-8 * InnerProduct(u, v) * dx).Assemble()
k_inverse = k.mat.Inverse(W.FreeDofs(), inverse="sparsecholesky")
fine_strains = []
for Ea in unit:
    w = GridFunction(W)
    w.vec.data = k_inverse * LinearForm(-InnerProduct(C(Ea), strain(v)) * dx).Assemble().vec
    fine_strains.append(Ea + strain(w))
fine_CH = [[Integrate(InnerProduct(C(fine_strains[i]), fine_strains[j]), fine) for j in (0, 1)] for i in (0, 1)]
print("clean cell: C^H", np.round(fine_CH, 4), "| nu_H", fine_CH[0][1] / sqrt(fine_CH[0][0] * fine_CH[1][1]),
      "| optimization", nu_h)
```

Tiled, the clean cells make the lattice, and a DXF of one cell, 10 mm across, is ready for a laser cutter (Figure 2):

```python
lattice = px.Unify([clean.Translated((i, j)) for i in range(4) for j in range(4)])
px.WriteDXF(clean.Scaled((0, 0), 10), "auxetic-cell.dxf", unit=px.mm, scale=1)  # scale=1: the numbers as they are
print("lattice: faces", len(lattice.faces), "| area", lattice.mass)
# px.Plot(lattice)
```

```python nbsphinx="hide-input"
Row((clean, "One cell"), (lattice, "Four by four cells"), caption="The clean cell and its lattice")
```

See [From level sets](../integrations/netgen/levelset.md) for `px.FromLevelSet` and for this cell extruded and pulled
in 3D, and [Identifications](../integrations/netgen/identifications.md) for periodic meshes.