Topology Optimization is a computational method for finding the best way to distribute material within a design space to create the most efficient structure while using the least resources. You can read more about it here. In this tutorial, we’ll explore topology optimization using FEniCSx.

FEniCSx is a powerful open-source tool that makes solving partial differential equations (PDEs) easy. One of its best features is that FEniCSx lets us write equations in code just like in math, thanks to the Unified Form Language (UFL). It also uses efficient libraries and tools like PETSc to help your code run faster.

By the end of this guide, you’ll know how to set up a topology optimization problem, perform Finite Element Analysis (FEA) using FEniCSx, and find the best material distribution using the optimality criteria method. But first, let’s start with some math!

---

Theory Behind the Magic

Linear Elasticity Equations

In linear elasticity with small deformations, the strain tensor \(\epsilon\) is the symmetric part of the gradient of the displacement vector \( u \) .

$$\epsilon = \frac{1}{2} \left( \nabla u + (\nabla u)^T \right)$$

The stress tensor \(\sigma\) is connected to the strain through the material's elastic properties, as described by Hooke’s law.

$$\sigma = \rho \left( \lambda \operatorname{tr} (\epsilon) I + 2\mu \epsilon \right)$$

where \(\lambda\) and \(\mu\) are Lamé’s parameters, \(\text{tr}\) is the trace operator, and \( I \) is the identity tensor. \(\rho\) represents the density field and determines the material distribution. The value of \(\rho\) ranges from 0 (indicating the void phase) to 1 (indicating the solid phase). In topology optimization, our goal is to find the optimal \(\rho\) .

Using the equilibrium equation, we can formulate the problem as:

$$\nabla \cdot \sigma + f = 0$$

where \( f \) is the body force per unit volume.

The Weak Form

In FEA, we work with the variational form of the PDEs. The variational formulation of the PDE consists of forming the inner product of the PDE with a vector test function \( v \) , and integrating over the domain \(\Omega\) :

$$\int_{\Omega} (\nabla \cdot \sigma) \cdot v \, dx + \int_{\Omega} f \cdot v \, dx = 0$$

Applying the divergence theorem, we get:

$$\small \int_{\Omega} \sigma : \nabla v \, dx = \int_{\Gamma} (\sigma \cdot n) \cdot v \, ds + \int_{\Omega} f \cdot v \, dx$$

where the colon operator is the inner product between tensors (summed pairwise product of all elements), \( n \) is the outward unit normal at the boundary \(\Gamma\) , and \( x \) is the spatial coordinates. The quantity \(\sigma⋅n\) is known as the traction or stress vector at the boundary, and often prescribed as a boundary condition.

The inner product of a symmetric tensor \( A \) and an anti-symmetric tensor \(B\) is zero. Therefore, we can replace \(\nabla v\) with its symmetric part \(\epsilon (v)\) . The equation above can be rewritten as:

$$\text{Bilinear} (\rho, u, v) = \text{Linear} (v)$$

where the Bilinear form is defined as:

$$\text{Bilinear} (\rho, u, v) = \int_{\Omega} \sigma (\rho, u) : \epsilon (v) \, dx$$

and the Linear form is:

$$\text{Linear} (v) = \int_{\Gamma} T \cdot v \, ds + \int_{\Omega} f \cdot v \, dx$$

This formulation can be easily implemented in FEniCSx to solve for the displacement field \( u \) .

The Optimization Problem

In topology optimization, we aim to find the optimal density field \(\rho\) that minimizes (or maximizes) the given objective function while meeting certain constraints. A well-known topology optimization problem is strain energy minimization under a volume constraint. This problem can be expressed as:

$$\begin{aligned} & \underset{\rho}{\text{minimize}} \quad \psi = \int_{\Omega} \sigma(\rho, u) : \epsilon(u) \, dx \\ & \text{subject to} \quad \int_{\Omega} \rho \, dV \leq \bar{V} \end{aligned}$$

where \( dV \) is the volume measure and \(\bar{V}\) is the prescribed volume fraction. A common approach in topology optimization is using the Solid Isotropic Material with Penalization (SIMP) method. In this method, the design variable \(\hat{\rho}\) is penalized to compute the density field as follows:

$$\rho(\hat{\rho}) = \hat{\rho}^p$$

where \( p \) is the SIMP penalty parameter. This penalization discourages intermediate values of the design variable, helping it converge towards 0 and 1.

Sensitivity Analysis

Topology optimization is an iterative process where the design variable is gradually updated to reach the optimal design. This update is achieved by calculating the gradient of the objective function (and possibly the constraints) and guiding the design towards the optimal solution.

In the problem of minimizing strain energy described in the previous section, the gradient of the objective function with respect to the design variable is:

$$\frac{\partial \psi}{\partial \hat{\rho}} = -p \hat{\rho} ^ {(p-1)} \psi$$

This gradient is then used by an optimization algorithm, such as the Optimality Criteria (OC) method or the Method of Moving Asymptotes (MMA), to gradually update the design variable.

Filter

There are a couple of challenges in topology optimization, including the formation of checkerboard patterns and mesh dependency. A common way to address these challenges is by applying a filter to the design or the sensitivities.

In each step of the optimization process, we use a filter to adjust an element's property by taking the weighted average of the properties of nearby elements. The nearby elements are those within a distance shorter than a specified value, \(r\). This is achieved using a convolutional operator.

$$\phi_i = \frac{\sum{\omega_{ij} \bar{\phi}_j}}{\sum{\omega_{ij}}}$$

where \(\phi\) is the property (like the design variable or sensitivity) and \(\omega_{ij}\) is the weight associated with the nearby element. The value of \(\omega_{ij}\) decreases as the distance between the nearby element (marked by red crosses) and the base element (marked with the blue cross) increases.

As you can see, applying this filter can prevent the formation of checkerboard patterns and makes the results independent of the mesh size.

---

Code

Okay! Enough with the theory. Let’s start coding…

We'll use NumPy and SciPy for basic math operations, FEniCSx for FEA, and Matplotlib for visualization, so make sure you have them installed. You can find the complete code on my GitHub repository.

Defining the problem

We begin by defining the design domain and creating the mesh. While there are several ways to create the mesh, I recommend using the Gmsh mesh generator. It's easy to use and free!

The problem we're tackling in this tutorial is a cantilever beam with a rectangular design area measuring \(1 \times 0.5\), as shown below. We've set up the domain geometry in the geo file and created the mesh.

Then, we load the mesh into the code and use the built-in function to find the midpoint coordinates of each cell.

domain, _, _ = read_from_msh('domain.msh', MPI.COMM_WORLD)
dim, facets_dim = domain.geometry.dim, domain.topology.dim - 1
num_cells = domain.topology.index_map(domain.topology.dim).size_local
midpoints = compute_midpoints(domain, domain.topology.dim, np.arange(num_cells))

Next, let's set up the problem parameters. We'll include the SIMP penalty parameter, the volume fraction, and the Lamé parameters.

p = 3
volume_fraction = 0.4
lame_lambda, lame_mu = 1.25, 1

To define the boundary conditions, we need to specify their location and value. In the cantilever beam problem, the left edge is fixed. Therefore:

def dirichlet_boundary(x):
    return np.isclose(x[0], 0, atol=0.05)

dirichlet_value = dolfinx.default_scalar_type((0, 0, 0))

and a point load is applied at the center of the right edge.

def neumann_boundary(x):
    return np.logical_and(np.isclose(x[0], 1, atol=0.05), np.isclose(x[1], 0.25, atol=0.05))

neumann_value = dolfinx.default_scalar_type((0, -1, 0))

To keep things flexible, we also define the body force.

body_force_value = dolfinx.default_scalar_type((0, 0, 0))

Formulating the PDE

We're working with two kinds of quantities here. First, we have discrete quantities, like the density field \(\rho\) and the strain energy field \(\psi\) . These are measured at each cell in a discrete way. Then, we have continuous quantities, like the displacement field \( u \) . These are measured at the vertices and interpreted using a first-order function.

To handle these, we need two different function spaces. The trial and test functions are set up in the continuous function space.

CG1 = functionspace(domain, ('CG', 1, (dim,)))
u = TrialFunction(CG1)
v = TestFunction(CG1)

and the strain energy and density functions are set up in the discrete function space. The density field is initialized with a uniform distribution.

DG0 = functionspace(domain, ('DG', 0))
energy = Function(DG0, name='StrainEnergy')
rho = Function(DG0, name='Density')
rho.x.array[:] = volume_fraction

The Dirichlet boundary conditions are applied by first identifying the boundary facets, then finding the degrees of freedom associated with those facets, and finally assigning values to them.

dirichlet_facets = locate_entities_boundary(domain, facets_dim, dirichlet_boundary)
dirichlet_dofs = locate_dofs_topological(CG1, facets_dim, dirichlet_facets)
dirichlet_bc = dirichletbc(dirichlet_value, dirichlet_dofs, CG1)

Applying the Neumann boundary conditions is a bit more complex. First, we locate the facets, then mark them with a value (like 1) and specify the integration measure ds.

neumann_facets = locate_entities_boundary(domain, facets_dim, neumann_boundary)
facet_values = np.ones(len(neumann_facets), dtype=np.int32)
facet_tags = meshtags(domain, facets_dim, neumann_facets, facet_values)
ds = Measure('ds', domain=domain, subdomain_data=facet_tags)
neumann_bc = Constant(domain, neumann_value)

The body force is applied across the entire domain, so its function is simply defined as:

body_force = Constant(domain, body_force_value)

Now let's go ahead and define the functions for strain, stress, and strain energy.

def epsilon(u):
    return sym(grad(u))

def sigma(rho, u):
    return rho * (lame_lambda * tr(epsilon(u)) * Identity(len(u)) + 2 * lame_mu * epsilon(u))

def psi(rho, u, v):
    return inner(sigma(rho, u), epsilon(v))

How cool was that! We just wrote the math exactly as it is, using the same symbols and notations. It's that simple, and it WORKS! It's like magic :)

Similarly, the bilinear and linear forms, along with the PDE, are defined in the same way. Keep in mind that 1 in ds(1) represents the tag for the marked boundary facets.

bilinear_form = psi(rho, u, v) * dx
linear_form = dot(body_force, v) * dx + dot(neumann_bc, v) * ds(1)
problem = LinearProblem(bilinear_form, linear_form, bcs=[dirichlet_bc],
                        petsc_options={'ksp_type': 'preonly', 'pc_type': 'cholesky'})

That's it! Now you can simply run problem.solve(), and FEniCSx will compute the displacement field for you.

Don’t Forget the Filter!

Using SciPy KDTree, we find the nearest neighboring elements and calculate the weight matrix omega.

def prepare_filter(coords, r):
    tree = KDTree(coords)
    distance = tree.sparse_distance_matrix(tree, r).tocsr()
    distance.data = (r - distance.data) / r
    return distance, np.array(distance.sum(1)).flatten()

omega, omega_sum = prepare_filter(midpoints, 0.02)

where omega is a square, symmetric sparse matrix in which only the elements within the radius \(r\) have non-zero values.

Optimization

Finally, we wrap things up with an iterative optimization process. But before that, there's a small step: we need to calculate the volume of each cell.

volume = TestFunction(DG0)
volume_vector = assemble_vector(form(volume * dx))
volume_vector.scatter_forward()
cell_volumes = volume_vector.array
total_volume = cell_volumes.sum()

In each iteration of the optimization process, we update the value of rho and use FEniCSx to do FEA and find the displacement field. Then, we perform sensitivity analysis and update the design variable using the OC method. The optimization process is considered converged when the number of intermediate elements in the design drops below a certain percentage.

losses = []
x = np.full(num_cells, volume_fraction)
for epoch in range(100):
    rho.x.array[:] = (x ** p).flatten()  # SIMP Interpolation
    disp = problem.solve()  # Finite Element Analysis

    # Objective Function
    energy_exp = Expression(psi(rho, disp, disp), DG0.element.interpolation_points())
    energy.interpolate(energy_exp)
    losses.append(energy.x.array.sum())

    # Sensitivity Analysis
    sensitivity = - p * (x ** (p - 1)) * energy.x.array
    # Filtering the Sensitivity
    sensitivity = (omega @ sensitivity) / omega_sum

    # Optimality Criteria Method
    ocp = x * np.sqrt(-sensitivity / cell_volumes)
    l1, l2, move = 0, ocp.sum(), 0.2
    while (l2 - l1) > 1e-4:
        l_mid = 0.5 * (l2 + l1)
        x = (ocp / l_mid).clip(x - move, x + move).clip(0.01, 1)
        volume = (x * cell_volumes).sum() / total_volume
        l1, l2 = (l_mid, l2) if volume > volume_fraction else (l1, l_mid)

    # Checking for Convergence
    if np.logical_and(0.15 < x, x < 0.85).sum() / num_cells < 1e-2: break

That’s it. You made it!

Now, let’s run our code and see what happens.

---

Results

Hit the run button, sit back, and enjoy the show!

Looks like everything is working as expected. The objective function is minimized while maintaining the volume constraint. Now, let's try it with a different problem, such as the tip-loaded cantilever problem.

All we need to do is change the load location like this:

def neumann_boundary(x):
    return np.logical_and(np.isclose(x[0], 1, atol=0.05), np.isclose(x[1], 0, atol=0.05))

Also, let's set the volume fraction to 0.5 this time.

Looks great to me! How about you?

---

Final Thoughts

By following this tutorial, you've laid down a solid foundation of knowledge. Now, it's your turn to take the reins and explore your creativity with what you've learned. Dive in, experiment, and make it your own!

If you encounter any challenges or have questions, don't hesitate to reach out. I'm here to help, and I'd love to hear from you! Leave a comment below or connect with me on social media. Happy coding and good luck! 🌟