Topology optimization of a cantilever beam in 3D¶

Open In Colab Binder

Topology optimization of a cantilever beam in 3D.

In [1]:
import torch
from scipy.optimize import bisect
from tqdm import tqdm

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

torch.set_default_dtype(torch.float64)

Model setup¶

We start by defining the base problem without considering the optimization yet.

In [2]:
# Material model
material = IsotropicElasticity3D(E=100.0, nu=0.3)

Nx = 20
Ny = 10
Nz = 15

# Create mesh
nodes, elements = cube_hexa(Nx + 1, Ny + 1, Nz + 1, Nx, Ny, Nz)

model = Solid(nodes, elements, material)

# Load at tip
tip = nodes[:, 0] == Nx
bottom = nodes[:, 2] == 0
model.forces[tip & bottom, 2] = -1.0
model.forces[tip & bottom & (nodes[:, 1] == 0), 2] = -0.5
model.forces[tip & bottom & (nodes[:, 1] == Ny), 2] = -0.5

# Constrained displacement at left end
model.constraints[nodes[:, 0] == 0.0, :] = True

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

# Plot
model.plot(u=u, node_property={"U": u}, bcs=True)

EmbeddableWidget(value='<iframe srcdoc="<!doctype html>\n<html lang=&quot;en&quot;>\n  <head>\n    <meta chars…

Optimization parameters¶

We define the optimization parameters, i.e. the volume fraction, the penalization factor, the move limit, the filter radius, and the number of iterations.

In [3]:
volfrac = 0.5
p = 3
move = 0.2
R = 1.5
TORCH_SENS = False
In [4]:
# Initial, minimum, and maximum values of design variables
rho_0 = volfrac * torch.ones(len(elements))
rho_min = 0.05 * torch.ones_like(rho_0)
rho_max = torch.ones_like(rho_0)

# Volume fraction
V_0 = volfrac * Nx * Ny * Nz

# Analytical gradient of the stiffness matrix
k0 = model.k0()
C0 = model.material.C.clone()

# Precompute filter weights
ecenters = nodes[elements].mean(dim=1)
dist = torch.cdist(ecenters, ecenters)
H = R - dist
H[dist > R] = 0.0

Optimization with optimality constraints.¶

This may take a minute to run.

In [5]:
rho = [rho_0]

# Iterate solutions
for k in tqdm(range(20)):
    # Adjust thickness variables
    model.material.C = torch.einsum("n,nijkl->nijkl", rho[k] ** p, C0)

    # Compute solution
    u_k, f_k, _, _, _ = model.solve()

    # Evaluation of compliance
    compliance = torch.inner(f_k.ravel(), u_k.ravel())

    # Compute analytical sensitivities
    u_j = u_k[elements].reshape(model.n_elem, -1)
    w_k = torch.einsum("...i, ...ij, ...j", u_j, k0, u_j)
    sensitivity = -p * rho[k] ** (p - 1.0) * w_k

    # Filter sensitivities (if r provided)
    sensitivity = H @ (rho[k] * sensitivity) / H.sum(dim=0) / rho[k]

    # For a certain value of mu, apply the iteration scheme
    def make_step(mu):
        G_k = -sensitivity / mu
        upper = torch.min(rho_max, (1 + move) * rho[k])
        lower = torch.max(rho_min, (1 - move) * rho[k])
        rho_trial = G_k**0.5 * rho[k]
        return torch.maximum(torch.minimum(rho_trial, upper), lower)

    # Constraint function
    def g(mu):
        rho_k = make_step(mu)
        return rho_k.sum() - V_0

    # Find the root of g(mu)
    with torch.no_grad():
        mu = bisect(g, 1e-10, 100.0)

    rho.append(make_step(mu))

100%|██████████| 20/20 [00:08<00:00,  2.27it/s]
In [6]:
model.plot(
    element_property={"rho": rho[-1]},
    clip=("rho", 0.5),
    show_edges=False,
    show_outline=True,
    color="lightblue",
    bcs=True,
)

EmbeddableWidget(value='<iframe srcdoc="<!doctype html>\n<html lang=&quot;en&quot;>\n  <head>\n    <meta chars…