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 maximum thickness $d_{max}=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.

Four bar truss

In [1]:
import matplotlib.pyplot as plt
import torch
from torchfem import Planar
from torchfem.materials import IsotropicPlaneStress
from tqdm import tqdm

torch.set_default_dtype(torch.double)

This is the planar FEM model from last exercise with an increased number of elements per direction $N=30$.

In [2]:
# Define Material
material = IsotropicPlaneStress(E=1000.0, nu=0.3)

# Create nodes
N = 30
n1 = torch.linspace(0.0, 1, N + 1)
n2 = torch.linspace(0.0, 1, N + 1)
n1, n2 = torch.stack(torch.meshgrid(n1, n2, indexing="xy"))
nodes = torch.stack([n1.ravel(), n2.ravel()], dim=1)

# Create elements connecting nodes
elements = []
for j in range(N):
    for i in range(N):
        n0 = i + j * (N + 1)
        elements.append([n0, n0 + 1, n0 + N + 2, n0 + N + 1])
elements = torch.tensor(elements)

# Create and plot the domain
square = Planar(nodes, elements, material)

# Define masks for boundary conditions
top = nodes[:, 1] == 1.0
left = nodes[:, 0] == 0.0
right = nodes[:, 0] == 1.0
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 = square.solve()

# Compute von Mises stress
mises = square.compute_stress(u, mises=True)

# Plot the result
square.plot(u=u, element_property=mises, cmap="inferno")
No description has been provided for this image

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.

In [3]:
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
In [4]:
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, xi=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.

In [5]:
def optimize(fem, rho_0, rho_min, rho_max, V_0, iter=50, xi=0.5, m=0.2, p=1.0, r=0.0):
    k0 = torch.einsum("i,ijk->ijk", 1.0 / fem.thickness, fem.k())
    rho = [rho_0]
    vols = d * fem.areas()

    # a) Check if there is a feasible solution before starting iteration

    # b) Precompute filter weights
    if r > 0.0:
        pass

    # Iterate solutions
    for k in tqdm(range(iter)):
        # c) Adjust thickness variables and compute FEM solution

        # d) Compute sensitivities

        # e) Filter sensitivities (if r provided)
        if r > 0.0:
            pass

        # f) For a certain value of mu, apply the iteration scheme
        def make_step(mu):
            pass

        # g) Constraint function
        def g(mu):
            pass

        # h) Find the root of g(mu)

        # Append variable to solution list

    return rho

i) Set up the initial design variables to $\rho_0=0.5, \rho_{min}=0.01, \rho_{max}=1.0$ for all elements and a volume constraint $V_0= 0.4 V_{max}$ with the maximum design volume $V_{max}$.

In [6]:
# Initial thickness, minimum thickness, maximum thickness

# Initial volume (40% of maximum design volume)

j) Perform the optimization with 50 iterations and the following parameters: $$p=3$$ $$r=0$$

In [ ]:

k) Plot the evolution of design variables vs. iterations. What does the graph tell you?

In [ ]:

l) Perform the optimization with 100 iterations and the following parameters $$p=3$$ $$r=0.05$$

In [ ]:

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.