FastFEM

Sina Atalay

2024

FastFEM was a class project for Software Engineering for Scientific Computing at Princeton University, taught by Henry F. Schreiner and Romain Teyssier, and built with Sacha Escudier and Kentaro Hanson.

The finite element method (FEM) is a numerical way to approximate solutions of partial differential equations (PDEs). FastFEM uses it to solve scalar two-dimensional (2D) parabolic PDEs. Such a problem asks for an unknown field $f(x,y,t)$ , such as temperature, over a 2D domain $\Omega$ . Here, $\Omega$ is the geometry itself: it may be a square, an irregular shape, or a shape with holes. The PDE applies inside $\Omega$ , while boundary conditions apply along its edges, $\partial\Omega$ . FEM transforms the PDE into a solvable system of algebraic equations. Solving that system gives a continuous finite-element field $f_h(x,y,t)$ that approximates the exact solution $f(x,y,t)$ . FastFEM provides the complete workflow, from constructing the geometry ( $\Omega$ ) to visualizing the solution on that geometry.

FastFEM solves PDEs of the form

\nabla^2 f = h(f)\frac{\partial f}{\partial t} + g(x,y)

where $h(f)$ is the coefficient of the time derivative and $g(x,y)$ is a source term. To solve the PDE, we must also specify how $f$ behaves along the boundary. On each boundary segment, the value of $f$ can be held fixed (Dirichlet), or its derivative perpendicular to the boundary can be set to zero (Neumann), meaning $f$ does not change as we move across the boundary.

FEM rewrites the original differential equation (strong form) by multiplying it by any allowed function $v(x,y)$ , integrating over the entire geometry $\Omega$ , and integrating the Laplacian by parts. This gives

-\int_\Omega \nabla f\cdot\nabla v\,d\Omega = \int_\Omega h(f)\frac{\partial f}{\partial t}v\,d\Omega + \int_\Omega gv\,d\Omega.

If this integral equation (weak form) is satisfied for every allowed $v$ , it is equivalent to the original differential equation.

FEM makes the integral equation computable by dividing the geometry $\Omega$ into a mesh and defining one simple basis function $\phi_i$ around each mesh node. A basis function is just an ordinary function: it is $1$ at its own node, changes linearly to $0$ at neighboring nodes, and is $0$ outside the mesh elements touching that node. All finite sums of these building blocks form a limited space of available fields, and FEM searches only within this space. Instead of solving for an arbitrary function $f$ , it solves for the finite set of coefficients that multiply the basis functions. It chooses those coefficients by requiring the integral equation to hold for every basis function. Intuitively, among all fields that can be constructed by adding the basis functions with different coefficients, this selects the one that best satisfies the original PDE. The original field problem has now become one equation per coefficient and, after discretizing time, a finite system of algebraic equations.

Using these basis functions, FEM represents the approximate field and its time derivative as

f_h(x,y,t)=\sum_i F_i(t)\phi_i(x,y), \qquad \frac{\partial f_h}{\partial t}=\sum_i \dot F_i(t)\phi_i(x,y).

FEM sets $v=\phi_j$ for each basis function in turn, giving one equation for each mesh node:

-\sum_i K_{ji}F_i = \sum_i M^{(h)}_{ji}(F)\dot F_i+b_j.

Putting these equations together produces the finite-element system

M^{(h)}(F)\dot F = -KF - b

with

K_{ji}=\int_\Omega \nabla\phi_i\cdot\nabla\phi_j\,d\Omega, \qquad M^{(h)}_{ji}(F)=\int_\Omega h(f_h)\phi_i\phi_j\,d\Omega, \qquad b_j=\int_\Omega g\phi_j\,d\Omega.

FastFEM evaluates the integrals on each mesh element and assembles the local contributions into the global matrices. Discretizing time then turns the problem into algebraic equations for $F$ . For example, when $h$ is constant and $M$ is the resulting mass matrix, one backward-Euler step is

\left(\frac{M}{\Delta t}+K\right)F^{n+1} = \frac{M}{\Delta t}F^n-b.

Solving this system gives the field at the next time step.

To use FastFEM, users construct the 2D geometry $\Omega$ from points, lines, surfaces, and holes; label its domains and boundaries; generate a mesh with Gmsh; and define the initial field, $h(f)$ , $g(x,y)$ , and boundary conditions. FastFEM assembles the finite-element system and advances the solution in time. Its PyVista-based visualization tools display the mesh and map the solved nodal values back onto the geometry as plots or animations.

The example below solves the scalar wave equation on the unit square:

\frac{\partial^2 u}{\partial t^2} = c^2 \nabla^2 u, \qquad c = 0.2.

The field starts as a localized displacement with zero initial velocity and zero derivative perpendicular to the square’s boundary.

Transient scalar wave displacement $u$ on the unit square for $0 \le t \le 10$. The colorbar shows $u$. — Transient scalar wave displacement $u$ on the unit square for $0 \le t \le 10$ . The colorbar shows $u$ .

FastFEM’s field operations support NumPy and JAX. Using JAX enables automatic differentiation, just-in-time (JIT) compilation, and execution on GPUs.

Read the original report.

Check out the source code on GitHub: github.com/fastfem/fastfem