Computation of Von Mises stresses

Hello,

I am trying to compute the Von Mises stresses from the displacement field. This is my script, mainly taking inspiration from

Implementation — FEniCSx tutorial (jsdokken.com)

``````import dolfinx, ufl, json
from petsc4py import PETSc
from mpi4py import MPI
import numpy as np

E = 1.0
nu = 0.3
mu = E / (2.0 * (1.0 + nu))
lmbda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))

def epsilon(v):

def sigma(v):
return 2.0 * mu * epsilon(v) + lmbda * ufl.tr(epsilon(v)) * ufl.Identity(len(v))

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10, dolfinx.mesh.CellType.quadrilateral)

V = dolfinx.fem.VectorFunctionSpace(mesh, ('CG', 1))
u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
a = dolfinx.fem.form(ufl.inner(sigma(u), ufl.nabla_grad(v)) * ufl.dx)

left_edge = dolfinx.mesh.locate_entities_boundary(mesh, 0, marker=lambda x: np.isclose(x[0], 0))
right_edge = dolfinx.mesh.locate_entities_boundary(mesh, 1, marker=lambda x: np.isclose(x[0], 1))

bc = dolfinx.fem.dirichletbc(np.array((0, ) * mesh.topology.dim, dtype=np.double), left_edge, V)

A = dolfinx.fem.petsc.assemble_matrix(a, bcs=[bc])
A.assemble()

facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
facet_indices.append(right_edge)
facet_markers.append(np.full(len(right_edge), 1))
facet_indices, pos = np.unique(np.array(np.hstack(facet_indices), dtype=np.int32), return_index=True)
facet_markers = np.array(np.hstack(facet_markers)[pos], dtype=np.int32)
facet_tag = dolfinx.mesh.meshtags(mesh, fdim, facet_indices, facet_markers)

ds = ufl.ds(domain=mesh, subdomain_data=facet_tag)
t = ufl.as_vector([1.0, 0.0])
L = dolfinx.fem.form(ufl.inner(t, v) * ds(1))
b = dolfinx.fem.petsc.assemble_vector(L)

uh = dolfinx.fem.Function(V)
ksp = PETSc.KSP().create(A.comm)
ksp.setOperators(A)
ksp.setType("preonly")
ksp.getPC().setType("lu")
ksp.getPC().setFactorSolverType("mumps")
vec = b.copy()
ksp.solve(b, vec)
uh.vector.setArray(vec.array)
uh.name = "Displacement"

sigma_dev = sigma(uh) - (1/3) * ufl.tr(sigma(uh)) * ufl.Identity(len(uh))
sigma_vm = ufl.sqrt((3/2) * ufl.inner(sigma_dev, sigma_dev))
mises_space = dolfinx.fem.FunctionSpace(mesh, ('DG', 0))
stress_expr = dolfinx.fem.Expression(sigma_vm, mises_space.element.interpolation_points())
stress = dolfinx.fem.Function(mises_space)
stress.interpolate(stress_expr)
``````

I am getting the following weird error:

`AttributeError: 'ExpressionIR' object has no attribute 'has_runtime_qr'`

when using `dolfinx.fem.Expression`. Any help? Thanks in advance.

What version of DOLFINx are you using?
I cannot reproduce this error with v0.7.3.

Hi @dokken,

Thanks for your answer. I am currently using a docker image based on v0.6.0.

The full error trace is below. Indeed, I just realized that the line of code which is currently giving an error is a modification to the standard FFCX backend, which allows to have runtime quadratures in FEniCSx. This is why I get this error whereas you don’t, probably.

Therefore, maybe a better question would be an alternative way to compute stresses? In legacy dolfin this was done with the `project` function, is there a similar feature in dolfinx?

``````Traceback (most recent call last):
File "/root/stress.py", line 63, in <module>
stress_expr = dolfinx.fem.Expression(von_Mises, V_von_mises.element.interpolation_points())
File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/fem/function.py", line 134, in __init__
self._ufcx_expression, module, self._code = jit.ffcx_jit(mesh.comm, (ufl_expression, _X),
File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/jit.py", line 56, in mpi_jit
return local_jit(*args, **kwargs)
File "/usr/local/dolfinx-real/lib/python3.10/dist-packages/dolfinx/jit.py", line 211, in ffcx_jit
r = ffcx.codegeneration.jit.compile_expressions([ufl_object], options=p_ffcx, **p_jit)
File "/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/jit.py", line 240, in compile_expressions
impl = _compile_objects(decl, expressions, expr_names, module_name, p, cache_dir,
File "/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/jit.py", line 259, in _compile_objects
_, code_body = ffcx.compiler.compile_ufl_objects(ufl_objects, prefix=module_name, options=options)
File "/usr/local/lib/python3.10/dist-packages/ffcx/compiler.py", line 106, in compile_ufl_objects
code = generate_code(ir, options)
File "/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/codegeneration.py", line 52, in generate_code
code_expressions = [expression_generator(expression_ir, options) for expression_ir in ir.expressions]
File "/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/codegeneration.py", line 52, in <listcomp>
code_expressions = [expression_generator(expression_ir, options) for expression_ir in ir.expressions]
File "/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/expressions.py", line 35, in generator
backend = FFCXBackend(ir, options)
File "/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/backend.py", line 36, in __init__
self.access = FFCXBackendAccess(ir, self.language, self.symbols,
File "/usr/local/lib/python3.10/dist-packages/ffcx/codegeneration/access.py", line 29, in __init__
self.has_runtime_qr = ir.has_runtime_qr
AttributeError: 'ExpressionIR' object has no attribute 'has_runtime_qr'
``````

Could you please print the following:

``````import ffcx
print(ffcx.__version__)
import dolfinx
print(dolfinx.__version__)
``````

could you also be so kind to supply a reference to what docker image you are using?

I would like to note that `has_runtime_qr` is not part of ffcx at the time of writing (or in 0.6.0 has that line:
ffcx/ffcx/codegeneration/access.py at v0.6.0 · FEniCS/ffcx · GitHub)
and I cannot find it in any version of ffcx.

Hi @dokken,

Indeed, as I wrote in the previous message (sorry, maybe it is not highlighted enough) I had to make some changes to the standard FFCX in order to have custom quadratures in FEniCSx (for CutFEM methods). Thus, it makes total sense that you do not see that line, because I added it myself.

That being said, would you be so kind to guide me in another approach to compute stresses without calling` dolfinx.fem.Expression`, which seems to be creating problems. I am doing this for the moment, but perhaps there is a better way:

``````import dolfinx, ufl, json
from petsc4py import PETSc
from mpi4py import MPI
import numpy as np

E = 1.0
nu = 0.3
mu = E / (2.0 * (1.0 + nu))
lmbda = E * nu / ((1.0 + nu) * (1.0 - 2.0 * nu))

def epsilon(v):

def sigma(v):
return 2.0 * mu * epsilon(v) + lmbda * ufl.tr(epsilon(v)) * ufl.Identity(
len(v)
)

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 100, 100, dolfinx.mesh.CellType.quadrilateral)
V = dolfinx.fem.VectorFunctionSpace(mesh, ('CG', 1))

u, v = ufl.TrialFunction(V), ufl.TestFunction(V)
a = dolfinx.fem.form(ufl.inner(sigma(u), ufl.nabla_grad(v)) * ufl.dx)
left_edge = dolfinx.mesh.locate_entities_boundary(mesh, 0, marker=lambda x: np.isclose(x[0], 0))
right_edge = dolfinx.mesh.locate_entities_boundary(mesh, 1, marker=lambda x: np.isclose(x[0], 1))
bc = dolfinx.fem.dirichletbc(np.array((0, ) * mesh.topology.dim, dtype=np.double), left_edge, V)

A = dolfinx.fem.petsc.assemble_matrix(a, bcs=[bc])
A.assemble()

facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
facet_indices.append(right_edge)
facet_markers.append(np.full(len(right_edge), 1))
facet_indices, pos = np.unique(np.array(np.hstack(facet_indices), dtype=np.int32), return_index=True)
facet_markers = np.array(np.hstack(facet_markers)[pos], dtype=np.int32)
facet_tag = dolfinx.mesh.meshtags(mesh, fdim, facet_indices, facet_markers)
ds = ufl.ds(domain=mesh, subdomain_data=facet_tag)

t = ufl.as_vector([1.0, 0.0])
L = dolfinx.fem.form(ufl.inner(t, v) * ds(1))
b = dolfinx.fem.petsc.assemble_vector(L)

uh = dolfinx.fem.Function(V)
ksp = PETSc.KSP().create(A.comm)
ksp.setOperators(A)
ksp.setType("preonly")
ksp.getPC().setType("lu")
ksp.getPC().setFactorSolverType("mumps")
vec = b.copy()
t = dolfinx.common.Timer()
ksp.solve(b, vec)
uh.vector.setArray(vec.array)
uh.name = "Displacement"

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 = dolfinx.fem.FunctionSpace(mesh, ("DG", 0))
stress = dolfinx.fem.Function(V_von_mises)
stress_trial, stress_test = ufl.TrialFunction(V_von_mises), ufl.TestFunction(V_von_mises)
a = ufl.inner(stress_trial, stress_test) * ufl.dx
L = ufl.inner(von_Mises, stress_test) * ufl.dx
projection_problem = dolfinx.fem.petsc.LinearProblem(a, L, [])
stress = projection_problem.solve()

file = dolfinx.io.XDMFFile(MPI.COMM_SELF, "example.xdmf", "w")
file.write_mesh(mesh)
file.write_function(stress)
``````

As you have modified FFCx, it is hard for me to give you any guidance as to how to bypass your modifications.

When you write `has_runtime_qr`, do you mean having different quadrature rules on different elements?

You can always project the stresses into an appropriate space, as you showed in your code above.

Hi @dokken,

Indeed, the problem stays in the modifications I did to FFCx. Thus I guess projecting as I did in the code above is a good enough solution for the moment. Thanks as usual for your help.

Yes, when I write `has_runtime_qr`, I mean having different quadrature rules on different elements. Actually, do you plan on integrating this feature in `dolfinx` sooner or later?

It is a feature that I find desirable.
@IgorBaratta has been doing some rewriting of FFCx that would make the integration of these tidier.

Hi everyone!

To clarify, the `has_runtime_qr` is part of my ffcx hack in GitHub - augustjohansson/ffcx-custom: ffcx custom quadrature to allow for passing cell-wise quadrature at runtime using the “custom” tabulate tensor function interface:

``````void tabulate_tensor_integral_3edb7c068402923a697d72e1b03e0957554f29c3(double* A,
const double* w,
const double* c,
const double* coordinate_dofs,
const int* entity_local_index,
Concerning your question, I might answer: `quadrature_normals` simply refers to the unit vectors normal to the surface, i.e. the normalized gradients of the level set function for an implicitly defined geometry. Hope this clarifies.