Implicit Gyroid structure

import torch
from torchfem import Solid
from torchfem.materials import IsotropicElasticity3D
from torchfem.mesh import cube_hexa
from torchfem.sdfs import Gyroid

import pyvista

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.1

# 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

gyroid = Gyroid().scale(scale)

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

Plot shell gyroid

model.plot(
    node_property={"SDF": gyroid.sdf(nodes)}, contour=("SDF", [-thickness, thickness])
)

# 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(gyroid.sdf(nodes)) > thickness
rho_nodes = torch.ones_like(gyroid.sdf(nodes))
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]

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

E_xx = torch.mean(σ[:, 0] / ε[:, 0])
vol_frac = (rho_elems * vol_elems).sum()
print(f"Effective E_xx ({100*vol_frac:.2f}%): {E_xx:.2f}")

Effective E_xx (8.27%): 79.21

model.plot(
    node_property={"Disp": u, "SDF": gyroid.sdf(nodes)},
    contour=("SDF", [-thickness, thickness]),
    scalars="Disp",
)