Implicit Gyroid structure

import pyvista
import torch

from torchfem import Solid
from torchfem.materials import IsotropicElasticity3D
from torchfem.mesh import cube_hexa

pyvista.set_plot_theme("document")
torch.set_default_dtype(torch.float64)

# Material parameters
E = 1000.0
nu = 0.3
material = IsotropicElasticity3D(E=E, nu=nu)

# Gyroid parameters
scale = 0.5 * torch.ones(3)
thickness = 0.2

# Mesh parameters
N = 51

# Homogenization parameters
disp = 0.1
rho_min = 0.01

Create voxel mesh

# Create a mesh
nodes, elements = cube_hexa(51, 51, 51)

# Create model
model = Solid(nodes, elements, material)

Evaluate signed distance function

def gyroid_sdf(points: torch.Tensor, scale: torch.Tensor) -> torch.Tensor:
    """Approximate signed distance to a gyroid surface with unit cell size `scale`."""
    x, y, z = (2 * torch.pi * points / scale).unbind(-1)
    f = (
        torch.sin(x) * torch.cos(y)
        + torch.sin(y) * torch.cos(z)
        + torch.sin(z) * torch.cos(x)
    )
    grad = torch.stack(
        [
            torch.cos(x) * torch.cos(y) - torch.sin(x) * torch.sin(z),
            torch.cos(y) * torch.cos(z) - torch.sin(x) * torch.sin(y),
            torch.cos(z) * torch.cos(x) - torch.sin(y) * torch.sin(z),
        ],
        dim=-1,
    )
    return f / (torch.norm(grad, dim=-1) + 1e-10)


sdf = gyroid_sdf(nodes, scale)

model.plot(node_property={"SDF": sdf}, cmap="coolwarm", clim=[-1.0, 1.0])

Plot shell gyroid

# Positive inside the wall
wall = thickness - sdf.abs()

model.plot(
    node_property={"Wall": wall},
    clip=("Wall", 0.0),
    color="lightblue",
    show_edges=False,
    show_outline=True,
)

# Set constraints
model.constraints[nodes[:, 0] == 0.0, 0] = True
model.constraints[nodes[:, 1] == 0.5, 1] = True
model.constraints[nodes[:, 2] == 0.5, 2] = True
model.constraints[nodes[:, 0] == 1.0, 0] = True
model.displacements[nodes[:, 0] == 1.0, 0] = 0.1

# Create nodal density field
mask = torch.abs(sdf) > thickness
rho_nodes = torch.ones_like(sdf)
rho_nodes[mask] = rho_min

# Integrate element density field
rho_elems = model.integrate_field(rho_nodes)
vol_elems = model.integrate_field(torch.ones_like(rho_nodes))
rho_elems /= vol_elems

model.plot(element_property={"Density": rho_elems}, cmap="coolwarm")

Compute homogenized properties

# Reduce stiffness with density field
model.material.C *= rho_elems[:, None, None, None, None]

# Solve
u, f, σ, F, α = model.solve()

ε = 0.5 * (F.transpose(-1, -2) + F) - torch.eye(3)
E_xx = torch.mean(σ[:, 0, 0] / ε[:, 0, 0])
vol_frac = (rho_elems * vol_elems).sum()
print(f"Effective E_xx ({100 * vol_frac:.2f}%): {E_xx:.2f}")

Effective E_xx (20.22%): 194.00

model.plot(
    node_property={"Disp": u, "Wall": wall},
    clip=("Wall", 0.0),
    scalars="Disp",
    show_edges=False,
    show_outline=True,
)