Plotting 3D piecewise constant stresses in paraview

Newbie to fenicsx here,
I’m trying to plot the Von Mises stresses found in the ‘The equations of linear elasticity’ tutorial in Paraview. The pyplot works fine in JupyterNotebooks, and when I save the Deformations into the xdmf file like this:

with io.XDMFFile(domain.comm, "out_test/deformation.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    uh.name = "Deformation"
    xdmf.write_function(uh)

Everything works fine in Paraview:

However when I try to save the stresses as well like this:

s = sigma(uh) -1./3*ufl.tr(sigma(uh))*ufl.Identity(len(uh))
von_Mises = ufl.sqrt(3./2*ufl.inner(s, s))
V_von_mises = fem.FunctionSpace(domain, ("DG", 0))
stress_expr = fem.Expression(von_Mises, V_von_mises.element.interpolation_points())
stresses = fem.Function(V_von_mises)
stresses.interpolate(stress_expr)

with io.XDMFFile(domain.comm, "out_test/deformation.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    uh.name = "Deformation"
    xdmf.write_function(uh)
    stresses.name = "stresses"
    xdmf.write_function(stresses)

(Bearing in mind this is the only change to the ipynb I made) Suddenly Paraview freaks out and doesn’t let me view either Deformation or stresses properly, they only appear as (partial):

Is it something to do with the stresses in a different FunctionSpace to the Deformation?

NB: correction, I also changed

L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds

to

L = ufl.inner(f, v) * ufl.dx + ufl.inner(T, v) * ds

as Jupyternotebooks doesn’t seem to like the complex numbers when using dot. This is why I get real_ and imag_ parts to both variable.

Use the extract block paraview filter to access each of the functions you have saved to the file.

Note that the stresses (if they are not in DG 0) are interpolated into (Lagrange, 1).

Thank you for the quick response!
Was struggling to get ExtractBlock to work but it turns out I was using Xdmf3ReaderS rather than Xdmf3ReaderT, switching to the correct one fixed it, thank you.

1 Like