Exercise 09 - Topology optimization for continua¶
Task 1 - Book shelf¶
Let us consider a bookshelf that needs a support structure. The design domain is given by a unit square $x \in [0, 1]^2$ and a thickness $d=0.1$. The left boundary of the domain $\partial \Omega_D$ is fixed to the wall and the top boundary $\partial \Omega_N$ is loaded with a uniform line load representing the weight of books.
import matplotlib.pyplot as plt
import torch
from torchfem import Planar
from torchfem.materials import IsotropicElasticityPlaneStress
from torchfem.mesh import rect_quad
from tqdm import tqdm
torch.set_default_dtype(torch.double)
# Create nodes and elements for a rectangular mesh
N = 30
L = 1.0
nodes, elements = rect_quad(Nx=N + 1, Ny=N + 1, Lx=L, Ly=L)
# Define Material
material = IsotropicElasticityPlaneStress(E=1000.0, nu=0.3)
# Create model
square = Planar(nodes, elements, material)
# Define masks for boundary conditions
top = nodes[:, 1] == L
left = nodes[:, 0] == 0.0
right = nodes[:, 0] == L
tip = top & right
# Load at top
square.forces[top, 1] = -1.0 / N
square.forces[tip, 1] = -0.5 / N
# Constrained displacement at left end
square.constraints[left, :] = True
# Thickness
d = 0.1
square.thickness[:] = d
# Solve the system
u, f, sigma, F, state = square.solve()
# Compute von Mises stress
mises = torch.sqrt(
sigma[:, 0, 0] ** 2
- sigma[:, 0, 0] * sigma[:, 1, 1]
+ sigma[:, 1, 1] ** 2
+ 3 * sigma[:, 1, 0] ** 2
)
# Plot the result
square.plot(
u=u,
element_property=mises,
cmap="viridis",
title="von Mises stress",
colorbar=True,
)
In addition, you are provided with a function that performs root finding with the bisection method and the computation of element surface areas from previous exercises.
def bisection(f, a, b, max_iter=50, tol=1e-12):
# Bisection method always finds a root, even with highly non-linear grad
i = 0
while (b - a) > tol:
c = (a + b) / 2.0
if i > max_iter:
raise Exception(f"Bisection did not converge in {max_iter} iterations.")
if f(a) * f(c) > 0:
a = c
else:
b = c
i += 1
return c
def compute_areas(truss):
areas = torch.zeros((truss.n_elem))
nodes = truss.nodes[truss.elements, :]
for w, q in zip(truss.etype.iweights, truss.etype.ipoints):
J = truss.etype.B(q) @ nodes
detJ = torch.linalg.det(J)
areas[:] += w * detJ
return areas
To save material, the bookshelf should use only 40% of the given design space, while being as stiff as possible to support many books without bending. We want to achieve this by topology optimization of the component.
To do so, we implement a topology optimization algorithm with optimality conditions in a function named optimize(fem, rho_0, rho_min, rho_max, V_0, iter=100, alpha=0.5, m=0.2, p=1.0, r=0.0) that takes the FEM model fem, the initial density distribution rho_0, the minimum and maximum thickness distributions rho_min, rho_max, the volume constraint, the maximum iteration count iter with a default value of 100, a SIMP penalty factor p with default 1, and a radius for sensitivity filtering r with a default 0.0.
a) Check if there is a feasible solution, i.e. if the design with minimum density has a volume smaller than the volume constraint. If not, raise an exception. Hint: You can compute the volume as the inner product of rho_min and the element volumes.
b) The filter weights can be precomputed before the optimization loop. Implement the computation of the filter weights if the radius is greater than 0.
Hints: Start by computing the center of each element and store it in a tensor of shape (Mx2) for M elements. Then, compute the distance between each element center using the function torch.cdist() and store it in a tensor of shape (MxM). Finally, compute the filter weights using the formula $H_{ij} = \max(0, r - d_{ij})$.
c) Add code that modifies the thickness according to the current design variables and solves the FEM problem in each iteration.
Hints: You can overwrite fem.thickness to set the thickness of the FEM object. You can use the thickness variable to inject the material density into the stiffness calculation and account for the penalty factor in the compliance calculation.
d) Compute the sensitivity of the compliance with respect to the design variables. Hints: This is equivalent to the compliance sensitivity in the previous exercise, but you need to account for the penalty factor.
e) Filter the sensitivity using the filter weights if the radius is greater than 0.
f) Define a function make_step(mu) that computes the design variable update for a given Lagrange parameter mu. Hints: Use the formula from the optimality conditions and clamp the result with the move limits.
g) Define a function g(mu) that evaluates the volume constraint. Hints: Use the make_step function to compute the design variable update and return the volume constraint violation.
h) Use the bisection method to find the Lagrange parameter that satisfies the volume constraint.
def optimize(
fem, rho_0, rho_min, rho_max, V_0, iter=100, alpha=0.5, m=0.2, p=1.0, r=0.0
):
k0 = torch.einsum("i,ijk->ijk", 1.0 / fem.thickness, fem.k0())
rho = [rho_0]
vols = d * compute_areas(fem)
# Check if there is a feasible solution before starting iteration
if torch.inner(rho_min, vols) > V_0:
raise Exception("rho_min is not compatible with V_0.")
# Precompute filter weights
if r > 0.0:
ecenters = fem.nodes[fem.elements].mean(dim=1)
dist = torch.cdist(ecenters, ecenters)
H = r - dist
H[dist > r] = 0.0
# Iterate solutions
for k in tqdm(range(iter), delay=0.1):
# Adjust thickness variables
fem.thickness = d * rho[k] ** p
# Compute solution
u, _, _, _, _ = fem.solve()
# Compute sensitivities
u_j = u[fem.elements, :].reshape(fem.n_elem, -1)
sens = (
-p * d * rho[k] ** (p - 1) * torch.einsum("...i,...ij,...j", u_j, k0, u_j)
)
# Filter sensitivities (if r provided)
if r > 0.0:
sens = (vols / rho[k]) * (H @ (rho[k] * sens / vols) / H.sum(dim=0))
# For a certain value of mu, apply the iteration scheme
def make_step(mu):
G_k = -sens / (mu * vols)
lower = torch.max(rho_min, (1 - m) * rho[k])
upper = torch.min(rho_max, (1 + m) * rho[k])
rho_trial = G_k**alpha * rho[k]
return torch.clamp(rho_trial, lower, upper)
# Constraint function
def g(mu):
rho_k = make_step(mu)
return torch.dot(rho_k, vols) - V_0
# Find the root of g(mu)
mu = bisection(g, 1e-10, 100.0)
rho.append(make_step(mu))
return rho
i) Set up the initial design variables $\rho_0=0.5, \rho_{min}=0.01, \rho_{max}=1.0$ and a volume constraint $V_0= 0.4 V_{max}$ with the maximum design volume $V_{max}$.
# Initial thickness, minimum thickness, maximum thickness
rho_0 = 0.5 * torch.ones(len(square.elements))
rho_min = 0.01 * torch.ones_like(rho_0)
rho_max = 1.0 * torch.ones_like(rho_0)
# Initial volume (40% of maximum design volume)
vols = d * compute_areas(square)
V0 = 0.4 * vols.sum()
j) Perform the optimization with the following parameters: $$ p=3 $$ $$ r=0 $$
# Optimize and visualize results
rho_opt = optimize(square, rho_0, rho_min, rho_max, V0, p=3.0, r=0.0)
square.plot(element_property=rho_opt[-1], cmap="gray_r")
100%|██████████| 100/100 [00:01<00:00, 56.55it/s]
k) Plot the evolution of design variables vs. iterations. What does the graph tell you?
plt.plot(torch.stack(rho_opt))
plt.xlabel("Iteration")
plt.ylabel("Values $ρ_j$")
plt.grid()
The graph shows that the design variables converge to 0 or 1 within the specified number of iterations, which is expected for the SIMP method.
l) Perform the optimization with 100 iterations and the following parameters $$ p=3 $$ $$ r=0.05 $$
# Optimize and visualize results
rho_opt = optimize(square, rho_0, rho_min, rho_max, V0, p=3.0, r=0.05)
square.plot(element_property=rho_opt[-1], cmap="gray_r")
100%|██████████| 100/100 [00:01<00:00, 53.16it/s]
m) How do you interpret the design? Decide which manufacturing process you would like to use and use a CAD software to create a design based on your optimization.
If we want to make the bookshelf from sheet metal, this could be an interpretation:
