FEM and Julia

code
simulation
Author

Lark

Published

December 17, 2025

FEM.jl

About

This repository is about implementing a \(n\)-Dimensional finite element method for solving differential equations in Julia.

Since I do not know much about this topic, this is mostly a learning project for myself.

And since this is a learning project, below there will be technical notes on implementation and theory in excruciating detail.

Theory

The basic problem we want to tackle is to solve a PDE on a domain \(\Omega\) with some boundary conditions. As an example I’ll use good old Poissons equation \[ -\nabla^2 u = f \] on the domain \(\Omega \in \mathbb{R}^n\) with \(u=0\) on its boundary. I may deviate from this later when neccesary.

Weak formulation

Trying to solve this equation directly on a computer has a major drawback, in that we require a high differentiability on our representation of the solution \(u\), depending on the order of the differential equation.

(As far as I understand,) there is another way to represent the problem which reduced (weakens) the differentiability requirement on the solution, which is called the weak formulation.

For the solver I want to program, you should probably already know the weak formulation of your problem, but generally the way you obtain it is by multiplying the PDE by a test function \(v\) and the integrating over the domain. In our example, this yields \[ -\int_{\Omega} (\nabla^2 u) v dx = \int_\Omega fvdx. \] Finally, integration by parts and assuming \(v=0\) on \(\partial \Omega\) yields the final weak formulation \[ \int_\Omega \nabla u \cdot \nabla v dx = \int_\Omega fvdx. \]

To abstract this a bit more, we can identify the left hand side to a bilinear form \(a(\cdot , \cdot ): V \times V \rightarrow \mathbb{R}\) and the right hand side with a linear form \(f(\cdot): V \rightarrow \mathbb{R}\). \(V\) is just any Hilbert space, but since we are looking for functions, our Hilbert space will be a function space.

Then the generic form of the weak formulation is: Find \(u\in V\) such that \[ a(u,v)=f(v) \qquad \forall v \in V. \]

This is called the Galerkin equation.

Galerkin dimension reduction

So far, we didn’t really gain anything. But now we are in a position to discretize the weak formulation. All we do is that we choose a subspace \(V_m \in V\) of dimension \(m\) and look at the projected problem:

Find \(u_m \in V_m\) such that for all \(v_m \in V_m\), \(a(u_m,v_m)=f(v_m)\).

Matrix form

Now we choose a basis \(\{e_i\}\) for \(V_m\). Since any \(v_m \in V_m\) can be constructed from just the basis, it is sufficient to fulfill the Galerkin equation for just the basis elements, so we have \[ a(u_m, e_i) = f(e_i) \qquad i=1,\dots,m. \] But clearly \(u_m \in V_m\) might just as well be expressed in the basis. Now we overload our notation a bit and say that \(u_m \in V_m\) is our solution function, but \(u_j \in \mathbb{R}\) are the coefficients for representing \(u_m\) in our basis: \[ u_m = \sum_{j=1}^m u_j e_j \]

If we insert that into the Galerkin equation as well, we obtain \[ a\left(\sum_{j=1}^mu_je_j,e_i\right) = \sum_{j=1}^m u_j a(e_j,e_i) =f(e_i) \qquad i=1,\dots,m, \] which is actually a linear system of equations \(Au=f\) with \(A_{ij}=a(e_j,e_i), \quad f_i=f(e_i)\).

We just need to plug that equation into a linear equation solver and we are done.

Basis functions

But hold on, what are these basis functions?

It seems that we often tile the domain with simplices or cubes, which define edges, faces, volumes, …, between the points.

The only choice of basis function I have seen so far is derived from the points we get from the tesselation of the domain. Each point brings with it a basis function that is \(1\) at that point, \(0\) at all other tesselation points, and linearly interpolates inbetween. This is explicitly NOT a delta function, since we require some differentiability in our function space.

For my own sanity, I will restrict the tesselation geoemtry to be the \(n\)-dimensional simplex, which should make interpolation manageable.

Implementation

Point

This is just a \(n\)-dimensional point.

Simplex

A collection of \(n\) \(n\)-dimensional points. Its main task is to provide a interpolate function.

Technical notes

Simplex math is quite interesting. Basically everything you need from a simplex can easily determined with the help of the so called barycentric coordinates.

Given the coordinates of the points in the simplex \(\{\vec{A}_i\}\) and a test point \(\vec{P}\), the barycentric coordinates \(\{a_i\}\) of \(\vec{P}\) are defined as the solution to \[ \begin{bmatrix} \vec{P}\\ 1 \end{bmatrix} = \begin{bmatrix} \vec{A}_0 \dots \vec{A}_n \\1\dots 1 \end{bmatrix}\begin{bmatrix} a_0\\\dots\\a_n \end{bmatrix} \]

These barycentric coordinates have some nice properties:

  1. If any barycentric coordinate is negative, the point is outside the simplex (=> Easy check if point is inside simplex)
  2. Linear interpolation within the simplex is performed by taking the dot product between the barycentric coordinates of a point and the values of the function at the points of the simplex.

Domain

I need to think about the implementation, but this is probably the cartesian product of intervals. This needs to somehow provide a tesselate function to obtain the simplices tesselating the domain and support for different tesselation strategies.

Interval

A closed interval on \(\mathbb{R}\).