How to make Von Mises Stress look better

[MattFerraro]

I followed this excellent tutorial on the equations of linear elasticity and added the ability to import and mesh a .step file with 2nd order tetrahedrons. I can apply boundary conditions and solve for displacement successfully:

But when I compute the von mises stress exactly as in the tutorial, the result looks terribly blocky:

I understand that this line in the tutorial: V_von_mises = fem.functionspace(domain, ("DG", 0)) means the output will be solved to have one value per tet and will therefore be discontinuous between tets. So I expect to see the kind of discontinuity highlighted with the red line, but what are the diamond artifacts marked with the yellow line? Furthermore, why does von mises stress have to be solved in a discontinuous function space instead of a continuous one like Lagrange?

When I crank up the mesh resolution I can get something that looks pretty good:

So I don’t think my math is wrong, I think I’m missing something conceptually. I’m trying to achieve results that look more like commercial solvers, like this:

What am I missing that makes the commercial results look so good and mine so poor? The commercial solver used a 2nd order mesh with 2050 tets. My implementation used a 2nd order mesh with 2258 tets and is significantly worse looking. Why?

[dokken]

The Von-Mises stresses are in the tutorial in a Dg-0 space as they contain the gradient of u through sigma in:

s = sigma(uh) - 1. / 3 * ufl.tr(sigma(uh)) * ufl.Identity(len(uh))
von_Mises = ufl.sqrt(3. / 2 * ufl.inner(s, s))

As uh in the tutorial is a piecewise linear function, the gradient is a cell-wise constant. Thus, the stresses are not in a piecewise continuous space.

Often this is overlooked in commerical software, or they use some “averaging” over at the vertices for all neighbouring cells.

You can try to put the stresses in a higher order DG space.
Especially if your geometry is second order, and you use a second order lagrange space for the displacement, DG-0 will not be sufficient to capture the von-Mises stresses accurately.

Alternatively, you could visualize the stresses in quadrature points, using scifem:
https://scientificcomputing.github.io/scifem/examples/xdmf_point_cloud.html
Then you can choose how many points to visualize by changing the quadrature degree.

[Stein]

@MattFerraro

What am I missing that makes the commercial results look so good and mine so poor?

That this is precisely the strength of using a tool like FEniCS over commercial software. If the commercial software produces a smooth-looking von-Mises stress while using C^0 (e.g. piecewise linear) basis functions, then that is nonphysical manipulation in the post-processing phase. It is simply “wrong”. Yet, customers want to see a nice looking graph, so that is why those packages do this.

FEniCS makes you think about what you do, and enables you to do this right.

(If the reference commercial simulation actually uses a 2nd order basis for the displacement field, then @dokken’s comment about using DG1 space for the von-Mises stress should give you a much nicer visualization, and I retract most of my statement.)

[MattFerraro]

Thank you for the help! Here is the code I’m running now, straight from the tutorial but modified so that u is a 2nd order Lagrange space:

# Scaled variable
import pyvista
from dolfinx import mesh, fem, plot, io, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
from mpi4py import MPI
import ufl
import numpy as np

L = 1
W = 0.2
mu = 1
rho = 1
delta = W / L
gamma = 0.4 * delta**2
beta = 1.25
lambda_ = beta
g = gamma

domain = mesh.create_box(MPI.COMM_WORLD, [np.array([0, 0, 0]), np.array([L, W, W])],
                         [20, 6, 6], cell_type=mesh.CellType.hexahedron)
# The 2 here is important!
V = fem.functionspace(domain, ("Lagrange", 2, (domain.geometry.dim, )))

# Boundary Conditions


def clamped_boundary(x):
    return np.isclose(x[0], 0)


fdim = domain.topology.dim - 1
boundary_facets = mesh.locate_entities_boundary(domain, fdim, clamped_boundary)
u_D = np.array([0, 0, 0], dtype=default_scalar_type)
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(
    V, fdim, boundary_facets), V)
T = fem.Constant(domain, default_scalar_type((0, 0, 0)))
ds = ufl.Measure("ds", domain=domain)

# Variational formulation


def epsilon(u):
    # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)
    return ufl.sym(ufl.grad(u))


def sigma(u):
    return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2 * mu * epsilon(u)


u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(domain, default_scalar_type((0, 0, -rho * g)))
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx + ufl.dot(T, v) * ds

# Solve
problem = LinearProblem(a, L, bcs=[bc], petsc_options={
                        "ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()


# Visualize the deformation
# pyvista.start_xvfb() Doesn't work when running from terminal on mac
p = pyvista.Plotter()
topology, cell_types, geometry = plot.vtk_mesh(V)
grid = pyvista.UnstructuredGrid(topology, cell_types, geometry)

# Attach vector values to grid and warp grid by vector
grid["u"] = uh.x.array.reshape((geometry.shape[0], 3))
actor_0 = p.add_mesh(grid, style="wireframe", color="k")
warped = grid.warp_by_vector("u", factor=1.5)
actor_1 = p.add_mesh(warped, show_edges=True)
p.show_axes()
if not pyvista.OFF_SCREEN:
    p.show()
else:
    figure_as_array = p.screenshot("deflection.png")

# Stress Computation
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)


# Visualize the stresses
warped.cell_data["VonMises"] = stresses.x.petsc_vec.array
warped.set_active_scalars("VonMises")
p = pyvista.Plotter()
p.add_mesh(warped, show_edges=True)
p.show_axes()
if not pyvista.OFF_SCREEN:
    p.show()
else:
    stress_figure = p.screenshot(f"stresses.png")

The von mises visualization then looks like this:

Where I added a red square around one of the 2nd order hex elements.

I have tried using V_von_mises = fem.functionspace(domain, ("DG", 1)) but that yields:

Traceback (most recent call last):
  File "/Users/matthewferraro/Documents/git/project/backend/linear_elasticity.py", line 91, in <module>
    warped.cell_data["VonMises"] = stresses.x.petsc_vec.array
  File "/Users/matthewferraro/miniconda3/envs/project/lib/python3.10/site-packages/pyvista/core/datasetattributes.py", line 247, in __setitem__
    self.set_array(value, name=key)
  File "/Users/matthewferraro/miniconda3/envs/project/lib/python3.10/site-packages/pyvista/core/datasetattributes.py", line 609, in set_array
    vtk_arr = self._prepare_array(data, name, deep_copy)
  File "/Users/matthewferraro/miniconda3/envs/project/lib/python3.10/site-packages/pyvista/core/datasetattributes.py", line 773, in _prepare_array
    raise ValueError(
ValueError: data length of (5760) != required length (720)

And if I try with higher orders in the “DG” family I get the same error but with larger numbers on the left side of the error.

If I try with V_von_mises = fem.functionspace(domain, ("Lagrange", 1)):

Traceback (most recent call last):
  File "/Users/matthewferraro/Documents/git/project/backend/linear_elasticity.py", line 93, in <module>
    warped.cell_data["VonMises"] = stresses.x.petsc_vec.array
  File "/Users/matthewferraro/miniconda3/envs/project/lib/python3.10/site-packages/pyvista/core/datasetattributes.py", line 247, in __setitem__
    self.set_array(value, name=key)
  File "/Users/matthewferraro/miniconda3/envs/project/lib/python3.10/site-packages/pyvista/core/datasetattributes.py", line 609, in set_array
    vtk_arr = self._prepare_array(data, name, deep_copy)
  File "/Users/matthewferraro/miniconda3/envs/project/lib/python3.10/site-packages/pyvista/core/datasetattributes.py", line 773, in _prepare_array
    raise ValueError(
ValueError: data length of (1029) != required length (720)

You mention that the gradient of u which is computed in epsilon may be to blame. How should I modify epsilon (and/or sigma) to make it possible to have von mises live in a continuous function space?

[dokken]

@MattFerraro

Traceback (most recent call last):
  File "/Users/matthewferraro/Documents/git/project/backend/linear_elasticity.py", line 91, in <module>
    warped.cell_data["VonMises"] = stresses.x.petsc_vec.array
  File "/Users/matthewferraro/miniconda3/envs/project/lib/python3.10/site-packages/pyvista/core/datasetattributes.py", line 247, in __setitem__
    self.set_array(value, name=key)
  File "/Users/matthewferraro/miniconda3/envs/project/lib/python3.10/site-packages/pyvista/core/datasetattributes.py", line 609, in set_array
    vtk_arr = self._prepare_array(data, name, deep_copy)
  File "/Users/matthewferraro/miniconda3/envs/project/lib/python3.10/site-packages/pyvista/core/datasetattributes.py", line 773, in _prepare_array
    raise ValueError(
ValueError: data length of (5760) != required length (720)

You get this error as it is no longer “cell_data”.
What I would do if you want to have this within pyvista is to:

1. Interpolate your displacement into a Vector DG 1 space.
2. Get a component of that space to place the stresses in.
   Pseudocode:

uh = dolfinx.fem.Function(V) # V being your second order vector space
Q = dolfinx.fem.functionspace(mesh, ("DG", 1, (3, ))
Q0 = Q.sub(0).collapse()[0]
stresses = dolfinx.fem.Function(Q0)
# Place your von mises stresses in the stress variable
# ....

# Create pyvista grid
u_grid = pyvista.UnstructuredGrid(dolfinx.plot.vtk_mesh(Q))
# Then add point data (vector for displacement, scalar for von mises stresses)

[dokken]

@MattFerraro

I modify epsilon (and/or sigma) to make it possible to have von mises live in a continuous function space?

You can’t with a finite element method, without using finite elements with much higher regularity (H^2). This is not commonly done in FEM.

[MattFerraro]

Thanks for the help!

For posterity, all I had to do was make sure the von mises function space was ("Lagrange", 1)

    # Works and makes the solution piecewise linear and continuous
    V_von_mises = fem.functionspace(domain, ("Lagrange", 1))

    stress_expr = fem.Expression(
        von_Mises, V_von_mises.element.interpolation_points())
    stresses = fem.Function(V_von_mises)
    stresses.interpolate(stress_expr)

Then to save the result to an XDMF file I encountered the problem that the mesh and the function must be of the same order. I solved that by interpolating the von mises stress to a 2nd order Lagrange space. Note that this adds no new information, it just makes it so there is a stress value for every single node in the 2nd order mesh.

    # Interpolate stresses onto a higher order function space (from lagrange 1 to lagrange 2)
    V_high = fem.functionspace(
        domain, ("Lagrange", 2))
    identity = ufl.Identity(1) * stresses
    identity_expr = fem.Expression(
        identity, V_high.element.interpolation_points())
    stresses_2 = fem.Function(V_high)
    stresses_2.interpolate(identity_expr)

Then one last detail is that when writing the file I had to use two different timestamps for the two different function outputs:

with io.XDMFFile(MPI.COMM_WORLD, solution_folder + "/u.xdmf", "w") as file:
    file.write_mesh(domain)
    file.write_function(uh, 0)
    file.write_function(stresses, 1)

I don’t think this is “right” but it got me a file that works. Here it is pulled into Paraview, displaying the von mises stress scalar using the default options and colormap:

This is exactly what I was looking for, and I think it is much better!

[dokken]

@MattFerraro

This is exactly what I was looking for, and I think it is much better!

This is not better from a mathematical point of view, as I tried to stress above.
By moving something that is discontinuous into a continuous space, you are performing an operation that is not compatible.
Imagine that you have the heaviside function, i.e.
H(x)=\begin{cases} 0 \quad \text{ if } x< 0 \\ 1 \quad \text{otherwise} \end{cases}
If you have a mesh node at x=0, what value should H(x) take?
in the case of von Mises stresses, you get into this situation everywhere the stress is not continuous. What you have chosen is an arbitrary smoothing of the problem, which will not be consistent in parallel.

[MattFerraro]

Gotcha, I did not understand that von mises stress really is discontinuous between cells. Given that is true, I see why enforcing continuity is in some sense “lying” about the data, which I do not want to do.

So instead of:

V_von_mises = fem.functionspace(domain, ("Lagrange", 1))

I switched to:

V_von_mises = fem.functionspace(domain, ("DG", 1))

as you suggested above. This preserves the ability for stresses to be discontinuous at cell boundaries, but allows the stresses to vary linearly over each cell which should look much better than cellwise-constant rendering.

This creates a problem saving the stresses to an XDMF file:

RuntimeError: Degree of output Function must be same as mesh degree. Maybe the Function needs to be interpolated?

I can interpolate to a DG 2 space but that leads to:

RuntimeError: Function and Mesh dof layouts do not match. Maybe the Function needs to be interpolated?

Which I think makes sense because if my stresses are in a DG 1 space, then I need to specify multiple values per mesh node, which seems incompatible with how an XDMF file is laid out.

Is there some way to achieve writing these stresses to an XDMF file without having to perform the arbitrary smoothing that collapses each node down to a single stress value? Would it be possible with a VTX file instead?

[dokken]

VTXWriter and VTKFile supports arbitrary order continuous and discontinuous Lagrange spaces, and should be used for outputting these variables.

If you want to show the von Mises stresses on the deformed geometry, I would interpolate the displacement into the same DG-1 space, i.e.

V_DG = dolfinx.fem.functionspace(mesh, ("DG", degree, (mesh.geometry.dim,)))
stresses = dolfinx.fem.Function(V_DG, name="VonMises")
u_dg = dolfinx.fem.Function(V_DG, name="u")
bp_vonmises = dolfinx.io.VTXWriter(mesh.comm, output / "von_mises.bp", [stresses, u_dg])
s = sigma(u, mu, lmbda) - 1.0 / 3 * ufl.tr(sigma(u, mu, lmbda)) * ufl.Identity(len(u))
von_Mises = ufl.sqrt(3.0 / 2 * ufl.inner(s, s))
stress_expr = dolfinx.fem.Expression(von_Mises, V_DG.element.interpolation_points())
stresses.sub(0).interpolate(stress_expr)
u_dg.interpolate(u)
bp_vonmises.write(0.0)
bp_vonmises.close()