The Weak Form of Linear Elasticity

This document provides an overview of the physics implemented in Scikit-Topt, focusing on linear elasticity and its formulation in the finite element method (FEM). We explain the governing equations, weak formulation, boundary conditions (Neumann and Dirichlet), and the resulting linear system solved inside the optimizer.

Linear Elasticity

In linear elasticity, the displacement field \(\mathbf{u} : \Omega \rightarrow \mathbb{R}^d\) describes how each point in the domain \(\Omega\) moves under external forces.

The stress–strain relationship follows Hooke’s law:

\[\boldsymbol{\sigma} = \mathbb{C} : \boldsymbol{\varepsilon},\]

where

  • \(\boldsymbol{\sigma}\) = Cauchy stress tensor

  • \(\boldsymbol{\varepsilon}\) = infinitesimal strain tensor

  • \(\mathbb{C}\) = 4th-order elasticity tensor (depends on Young’s modulus \(E\) and Poisson ratio \(\nu\)).

The strain tensor is defined as:

\[\boldsymbol{\varepsilon}(\mathbf{u}) = \frac{1}{2}(\nabla \mathbf{u} + \nabla\mathbf{u}^T).\]

Governing Equation (Strong Form)

Mechanical equilibrium in absence of body forces is written as:

\[-\nabla \cdot \boldsymbol{\sigma} = \mathbf{0} \quad \text{in } \Omega.\]

Boundary conditions are applied on the domain boundary \(\partial\Omega\):

  • Dirichlet BC (fixed boundary): \(\mathbf{u} = \bar{\mathbf{u}}\) on \(\Gamma_D\)

  • Neumann BC (traction/force): \(\boldsymbol{\sigma}\mathbf{n} = \mathbf{t}\) on \(\Gamma_N\)

Here \(\mathbf{n}\) is the outward normal vector.

Weak Formulation

We multiply the equilibrium equation by a test function \(\mathbf{v}\) and integrate over the domain:

\[\int_\Omega \boldsymbol{\sigma}(\mathbf{u}) : \boldsymbol{\varepsilon}(\mathbf{v}) \, d\Omega = \int_{\Gamma_N} \mathbf{t} \cdot \mathbf{v} \, d\Gamma.\]

This is the weak form of linear elasticity.

Substituting Hooke’s law gives:

\[\int_\Omega (\mathbb{C}\boldsymbol{\varepsilon}(\mathbf{u})) : \boldsymbol{\varepsilon}(\mathbf{v}) \, d\Omega = \int_{\Gamma_N} \mathbf{t} \cdot \mathbf{v} \, d\Gamma.\]

Finite Element Discretization

We approximate the displacement field as:

\[\mathbf{u}_h(\mathbf{x}) = \sum_{i=1}^{N} u_i \, \phi_i(\mathbf{x}),\]

where \(\phi_i\) are FEM shape functions and \(u_i\) the unknown nodal degrees of freedom.

Inserting this into the weak form yields the discrete system:

\[\mathbf{K} \mathbf{u} = \mathbf{f},\]

where:

  • \(\mathbf{K}\) is the stiffness matrix

    \[\mathbf{K}_{ij} = \int_\Omega \boldsymbol{\varepsilon}(\phi_j)^T \, \mathbb{C} \, \boldsymbol{\varepsilon}(\phi_i) \, d\Omega.\]

  • \(\mathbf{f}\) is the load vector

    \[\mathbf{f}_i = \int_{\Gamma_N} \mathbf{t} \cdot \phi_i \, d\Gamma.\]

Boundary Conditions

Dirichlet BC (fixed boundaries)

Enforced by modifying the stiffness matrix system:

  • Rows corresponding to fixed DOFs are zeroed.

  • Columns corresponding to fixed DOFs are also zeroed to preserve symmetry.

  • Diagonal entries are set to 1

  • The right-hand side is corrected as \(\mathbf{f} \leftarrow \mathbf{f} - \mathbf{K}_{:,i}\,\bar{u}_i\) so that the prescribed displacement is enforced consistently.

This approach follows the enforce() function in scikit-fem, which modifies the full matrix without reducing its size. An important advantage is that the resulting stiffness matrix remains identical across all load cases, allowing a single LU factorization to be reused efficiently. For this reason, Scikit-Topt uses enforcement rather than matrix condensation.

Neumann BC (surface forces)

Traction loads add contributions to \(\mathbf{f}\):

Solving the Linear System

After applying boundary conditions, we solve:

\[\mathbf{K}_e \mathbf{u} = \mathbf{f}_e.\]

In Scikit-Topt:

  • Default: sparse direct solver (scipy.spsolve or splu)

  • Optional: iterative solvers or PyAMG

Typical code:

K_e, f_e = skfem.enforce(K, f, D=dirichlet_dofs)
u = scipy.sparse.linalg.spsolve(K_e, f_e)