u/perebal — reddlx

Space-time FEM for elastic wave propagation — no time stepping

I was experimenting with a space-time finite element formulation for a simple elastodynamics problem and thought the result might be interesting here.

The example is a 1D elastic bar immediately after impact with a rigid wall. I do not model the impact itself; the calculation starts from the post-impact initial state and follows the subsequent stress-wave propagation, reflection and release.

Instead of discretizing space first and then advancing the solution in time, I introduce

y = ct

and treat (x,y) as an ordinary 2D finite-element domain.

Using particle velocity v and normalized stress

s = σ/(ρc)

the first-order system becomes

∂v/∂y − ∂s/∂x = 0
∂s/∂y − ∂v/∂x = 0

I used a least-squares formulation, which leads to four bilinear forms:

Kvv = ∫(Grad(V) ⋅ I ⋅ Grad(V))
Kvs = ∫(Grad(V) ⋅ C ⋅ Grad(S))

Ksv = ∫(Grad(S) ⋅ C ⋅ Grad(V))
Kss = ∫(Grad(S) ⋅ I ⋅ Grad(S))

The complete coupled system is then just

K = SystemMatrix([
    Kvv  Kvs
    Ksv  Kss
])

followed by one solve.

There is no time-stepping loop. The complete evolution over the chosen time interval is solved as one space-time finite-element problem.

What I especially like about the result is that the two wave fronts appear directly as characteristic lines in the (x,ct) domain.

The solution reproduces the classical 1D result: after impact, a compressive wave travels from the constrained end toward the free end. It reflects there as a release wave and travels back toward the wall.

During the compressed phase,

σ = -ρ c v₀,

and the release wave returns to the wall at

t = 2L/c.

The implementation uses LowLevelFEM.jl, but the notebook also contains the derivation of the formulation.

I'd be interested in thoughts from people who have worked with space-time FEM or least-squares formulations for hyperbolic problems.

The notebook is available via the link in the first comment.

u/perebal — 1 day ago

▲ 40 r/Julia

An executable FEM weak form in Julia is now faster than my original problem-specific implementation

One of the things I wanted to achieve with LowLevelFEM.jl was to keep finite element code reasonably close to the mathematical formulation.

For example, the stiffness matrix for a 3D linear elasticity problem can be written as

K = ∫(SymGrad(Pu) ⋅ D ⋅ SymGrad(Pu))

and the surface load as

f = ∫(Pu ⋅ [1.0, 0.0, 0.0], Γ="right")

rather than calling a dedicated elasticity assembly routine.

Originally, I considered this mainly an abstraction/readability feature. I expected the more general operator-based formulation to come with some performance cost.

After working on the assembly implementation — particularly direct assembly into a precomputed CSC sparsity pattern and multithreading — that is no longer necessarily the case.

In a small 3D elasticity example on my machine:

problem-specific high-level solve: 323 ms, 299 MiB
operator/weak-form solve: 120 ms, 52 MiB

Even the stress recovery can be written directly as field algebra:

ε = (u ∘ ∇ + ∇ ∘ u) / 2

σ = E / (1 + ν) * (ε + ν / (1 - 2ν) * trace(ε) * I)

For this particular example, that version is also slightly faster and uses considerably less memory than the older dedicated stress routine.

I don't mean these numbers as a general benchmark — they are just one mesh and one machine. What I find interesting is that the more general formulation no longer seems to require choosing between readable mathematical notation and reasonable performance.

For me, that was an important milestone in the development of the package.

I'd be interested in what people working on FEM/PDE software think about this kind of operator-level interface, especially where you would draw the line between mathematical expressiveness and implementation transparency.

reddit.com

u/perebal — 5 days ago

▲ 30 r/Julia

Writing finite element weak forms almost exactly as in textbooks (LowLevelFEM.jl)

One thing that has always bothered me when implementing finite elements is how quickly the code diverges from the mathematical formulation.

Over the last few months I've been experimenting with a Julia DSL where the weak form itself becomes executable code instead of something that first has to be translated into loops and element matrices.

For example, a standard bilinear form can be written as

K = ∫(Grad(Pu) ⋅ C ⋅ Grad(Pu))

while more complicated formulations are written in essentially the same style:

K = ∫((A⋅Grad(Pu) + G⋅Pu)' ⋅ C ⋅ (A⋅Grad(Pu) + G⋅Pu) * (2π*r))

B = A⋅Grad(Pu) + G⋅Pu
K = ∫(B' ⋅ C ⋅ B * (2π*r))

These are not symbolic expressions or macros generating another language. They are actual Julia expressions assembled directly into finite element matrices.

The same mechanism currently supports

scalar, vector and tensor fields,
multifield formulations,
user-defined operators,
variable coefficients,
full and reduced integration,
arbitrary operator compositions.

The interesting part for me was not only making the syntax close to the mathematics, but also keeping it reasonably efficient. After some recent refactoring, the compound operator assembly became significantly faster while keeping exactly the same high-level notation.

The package (LowLevelFEM.jl) has recently been published in JOSS.

GitHub: github.com/perebalazs/LowLevelFEM.jl

More examples can be found in the documentation:

perebalazs.github.io/LowLevelFEM.jl/stable/tutorials

I'd be interested to hear how others approach this.

Do you prefer writing PDEs as executable operator expressions like these, or do you find more explicit element-level assembly easier to understand, debug and maintain?

reddit.com

u/perebal — 19 days ago