Computation of Von Mises stresses

federicobetti · February 3, 2024, 9:25am

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): 
       return ufl.sym(ufl.nabla_grad(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.

dokken · February 6, 2024, 8:13pm

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

Please also add the full error trace of your code.

federicobetti · February 7, 2024, 8:10am

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'

dokken · February 7, 2024, 10:26am

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.

federicobetti · February 7, 2024, 3:19pm

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):
return ufl.sym(ufl.nabla_grad(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)

dokken · February 7, 2024, 4:04pm

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.

federicobetti · February 8, 2024, 9:07am

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?

dokken · February 8, 2024, 4:14pm

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

August · February 8, 2024, 9:03pm

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,
				    const uint8_t* quadrature_permutation,
				    int num_quadrature_points,
				    const double* quadrature_points,
				    const double* quadrature_weights,
				    const double* quadrature_normals)

This means that for every integration over a domain \Omega not fitting the mesh, one must pass the appropriate quadrature. This includes L2 projections.

dokken · February 8, 2024, 9:30pm

This kind of interface aligns with what I had in mind a a second signature for tabulate_tensor.
What are quadrature normals in your code?

federicobetti · February 10, 2024, 11:13am

Hi @dokken,

I would like to thank August for clarifying. Indeed, for the development of my applications, I used his library extensively over the last few months.

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.

Topic		Replies	Views
Issue with Von Mises stress calculation General	5	287	November 15, 2023
Operator on an integral	1	283	September 29, 2021
Von mises plasticity displacement control problem General	0	105	January 9, 2024
Convert Function to UFL type dolfinx	18	1139	May 8, 2023
Stress contour in frequency domain	0	318	March 4, 2022

Computation of Von Mises stresses

Related Topics