How to write stresses to a xdmf file?

General
1

Newbie here. I tried to follow some tutorials and wrote the following code for a cantilever beam under a uniform load. However, I am not able to output the stresses to a xdmf file. I tried to write a syntax similar to what I found for displacement in the tutorials, but that does not help.
Here is the code snippet

problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
stress = sigma(uh)
with io.XDMFFile(domain.comm, "simple_beam.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    uh.name = "Deformation"
    xdmf.write_function(uh)
    sigma.name = 'Stresses'
    xdmf.write_function(sigma)

Here is the error message

Traceback (most recent call last):
  File "/Users/aaquib/Desktop/fenics/test1.py", line 60, in <module>
    xdmf.write_function(sigma)
  File "/Users/aaquib/opt/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/dolfinx/io/utils.py", line 150, in write_function
    super().write_function(getattr(u, "_cpp_object", u), t, mesh_xpath)
TypeError: write_function(): incompatible function arguments. The following argument types are supported:
    1. (self: dolfinx.cpp.io.XDMFFile, function: dolfinx.cpp.fem.Function_float64, t: float, mesh_xpath: str) -> None
    2. (self: dolfinx.cpp.io.XDMFFile, function: dolfinx.cpp.fem.Function_complex128, t: float, mesh_xpath: str) -> None

Invoked with: <dolfinx.io.utils.XDMFFile object at 0x166518900>, <function sigma at 0x1664e4430>, 0.0, "/Xdmf/Domain/Grid[@GridType='Uniform'][1]"

I saw some old questions based on legacy version of fenics, where TensorFunctionSpace was recommended, but apparently that does not work with fenicsx?

2

Could you provide a minimum working example? It’s tough to comment without knowing what stress or sigma are defined as, for example.

3

Sure,

L = 1
h = 0.1
mu = 1e10
#rho = 1e6
#delta = W/L
#gamma = 0.4*delta**2
#beta = 1.25
lambda_ = 1e10
#g = gamma
w= 1e6


import numpy as np
import ufl

from mpi4py import MPI
from petsc4py.PETSc import ScalarType

from dolfinx import mesh, fem, plot, io
model='plane_stress'
domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0,0]), np.array([L, h])],
                  [20,6], cell_type=mesh.CellType.quadrilateral)
V = fem.VectorFunctionSpace(domain, ("CG", 1))

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], dtype=ScalarType)
bc = fem.dirichletbc(u_D, fem.locate_dofs_topological(V, fdim, boundary_facets), V)

#T = fem.Constant(domain, ScalarType((0, 0)))
#ds = ufl.Measure("ds", domain=domain)

if model=='plane_stress':
    lambda_=2*mu*lambda_/(lambda_+2*mu)

def epsilon(u):
    return ufl.sym(ufl.grad(u)) # Equivalent to 0.5*(ufl.nabla_grad(u) + ufl.nabla_grad(u).T)
def sigma(u):
    return lambda_ * ufl.nabla_div(u) * ufl.Identity(len(u)) + 2*mu*epsilon(u)
print(fdim)
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = fem.Constant(domain, ScalarType((0, -w)))
a = ufl.inner(sigma(u), epsilon(v)) * ufl.dx
L = ufl.dot(f, v) * ufl.dx 

problem = fem.petsc.LinearProblem(a, L, bcs=[bc], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
stress = sigma(uh)
with io.XDMFFile(domain.comm, "simple_beam.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    uh.name = "Deformation"
    xdmf.write_function(uh)
    sigma.name = 'Stresses'
    xdmf.write_function(sigma)

Thanks for your help!

4

You should use dolfinx.fem.Expression to interpolate the stress into a suitable function space (DG of some order) and use dolfinx.io.VTXWriter to write it to file preserving the function space.
You can see how B is computed in Electromagnetics example — FEniCSx tutorial
E.g.

W = VectorFunctionSpace(mesh, ("DG", 0))
B = Function(W)
B_expr = Expression(as_vector((A_z.dx(1), -A_z.dx(0))), W.element.interpolation_points())
B.interpolate(B_expr)

https://jsdokken.com/fenics22-tutorial/heat_eq.html#post-processing-without-projections
or
https://jsdokken.com/dolfinx-tutorial/chapter2/linearelasticity_code.html#stress-computation

2 Likes
5

Thanks for the pointers!
Per your recommendation, I added the following bit to my code

W=fem.VectorFunctionSpace(domain,("DG",0))
sigma_xx=fem.Function(W)
sigma_xx_exp= fem.Expression(lambda_ * ufl.nabla_div(uh)  + 2*mu*ufl.sym(ufl.grad(uh))[0,0],W.element.interpolation_points())
sigma_xx.interpolate(sigma_xx_exp)

As you can see, I am trying to save each component of the stress as a scalar (here I just give the code for one component). The first three lines compile successfully, but the fourth gives the following error

Traceback (most recent call last):
  File "/Users/aaquib/Desktop/fenics/test1.py", line 61, in <module>
    sigma_xx.interpolate(sigma_xx_exp)
  File "/Users/aaquib/opt/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/dolfinx/fem/function.py", line 343, in interpolate
    _interpolate(u, cells)
  File "/Users/aaquib/opt/anaconda3/envs/fenicsx-env/lib/python3.10/functools.py", line 889, in wrapper
    return dispatch(args[0].__class__)(*args, **kw)
  File "/Users/aaquib/opt/anaconda3/envs/fenicsx-env/lib/python3.10/site-packages/dolfinx/fem/function.py", line 334, in _
    self._cpp_object.interpolate(expr._cpp_object, cells)
RuntimeError: Function value size not equal to Expression value size
6

You are currently trying to interpolate a scalar into a vector function space.

Change

to

W=fem.FunctionSpace(domain,("DG",0))
1 Like
7

Thank you both @nav and @dokken ! This solved my issue and its unfortunate that this site wont allow me to mark both of your comments as solution…
Anyhow, I am just attaching what worked for me so that people who come in later can benefit

W=fem.TensorFunctionSpace(domain,("DG",0))
sigma=fem.Function(W)
sigma_exp= fem.Expression(lambda_ * ufl.nabla_div(uh)* ufl.Identity(len(u))  + 2*mu*ufl.sym(ufl.grad(uh)),W.element.interpolation_points())
sigma.interpolate(sigma_exp)

with io.XDMFFile(domain.comm, "simple_beam_stress.xdmf", "w") as xdmf:
    xdmf.write_mesh(domain)
    sigma.name = "sigma"
    xdmf.write_function(sigma)
1 Like