Divergence of the resulting stresses is always zero

mathematics
1

Hello everyone,

I’m running the tutorial on linear elasticity (Implementation — FEniCSx tutorial).
I added a function to compute the divergence of the stresses given a displacement field, i.e.

def div_sigma(u):
    return ufl.nabla_div(sigma(u))

which I call and plot with

div_sigs = div_sigma(uh)
divS = Function(V)
divSigs_expr = fem.Expression(div_sigs, V.element.interpolation_points())
divS.interpolate(divSigs_expr)

with XDMFFile(MPI.COMM_WORLD, "divSig.xdmf", "w") as file:
    file.write_mesh(domain)
    file.write_function(divS)

where $$uh$$ is the solution to the linear problem.

This might be a naive question, but given

\begin{split} -\nabla \cdot \sigma (u) &= f && \text{in } \Omega\\ \sigma(u)&= \lambda \mathrm{tr}(\epsilon(u))I + 2 \mu \epsilon(u)\\ \epsilon(u) &= \frac{1}{2}\left(\nabla u + (\nabla u )^T\right)\end{split}

I’d expect the divergence of the stresses to equal the negative body forces $$-f$$. However, the result is always zero in every coordinate.

Am I missing something obvious?

//Lucas

2

V is a P^1 space, hence div(sigma(u)), which contains second order derivatives, is zero.

3

oooh sure, that makes perfect sense, thank you. Do you think it makes sense to project the solution uh into some higher order space and then compute div(sigma(u))? If so, how would that work?

4

No. If it’s a simple projection or interpolation on the same spatial discretisation, note that P^1 is contained in P^2. You’ll just be doing the same computation.

5

makes sense, too. This is solved, thank you!