I’ve been working on a small Julia FEM framework called LowLevelFEM.jl. One of the ideas behind the project is to assemble PDE systems directly from their weak form.
The core idea is simple:
If you know the weak form of a PDE, you can almost copy it from paper into executable code.
The framework provides a small set of variational operators (Grad, Div, SymGrad, etc.) that allow weak forms to be written in a way that closely resembles their mathematical formulation.
Example: steady Navier–Stokes
μ = 1.0
γ = 1e-1 # grad-div stabilization
δ = 1e-6 # pressure stabilization
A = ∫((ε(Pv) ⋅ ε(Pv)) * 2μ)
D = ∫(Div(Pv) ⋅ Div(Pv) * γ)
B = ∫(Div(Pv) ⋅ Pp)
C = ∫(Grad(Pp) ⋅ Grad(Pp) * δ)
K = SystemMatrix([
A + D B'
B -C ])
The snippet above is not pseudocode — it is the actual solver code.
The same formulation can be used for different PDE systems such as
diffusion problems
linear elasticity
nonlinear elasticity
fluid flow
multiphysics coupling
Here is a small example result from the Navier–Stokes test case:
Sounds interesting, reminds me of FEniCSx, which I use heavily. However, there’s already Gridap.jl - I wonder how your package compares with it?
I have also wondered whether it would make sense to build something like this on top of Ferrite.jl by leveraging Symbolics.jl to generate assembly code from the weak form. I’m curious to hear your thoughts on this? Did you consider such an approach?
Good question — there are indeed similarities to frameworks like FEniCS or Gridap.
One of the motivations behind LowLevelFEM was to build something transparent and open rather than a “black box” solver. I wanted a codebase that I can use both for teaching and research, where the numerical formulation is clearly visible and easy to modify.
The goal is not to hide the numerical method, but to expose it in a clear and programmable way.
In that sense the framework behaves more like a toolbox for PDE experimentation than a traditional solver package.
The framework therefore keeps the variational operators explicit, so the weak form is written almost directly in code, e.g.
A = ∫(2μ * (ε(Pu) ⋅ ε(Pu)))
B = ∫(Div(Pv) ⋅ Pp)
The idea is that you can work at different levels: use it as a straightforward engineering solver, or go deeper and modify the operators and assembly if needed.
Regarding Ferrite + Symbolics: that’s an interesting direction. At the moment, the operators are implemented directly in Julia rather than generated symbolically, but I could definitely imagine experimenting with that approach in the future.
