# How to write stresses to a xdmf file?

**URL:** <https://fenicsproject.discourse.group/t/how-to-write-stresses-to-a-xdmf-file/11306>\
**Category:** General\
**Created:** [May 26, 2023, 4:46pm UTC](https://fenicsproject.discourse.group/t/how-to-write-stresses-to-a-xdmf-file/11306 "2023-05-26T16:46:56Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![aaquib](https://avatars.discourse-cdn.com/v4/letter/a/d2c977/32.png) [@aaquib](https://fenicsproject.discourse.group/u/aaquib)\
**Post date:** [May 26, 2023, 4:46pm UTC](https://fenicsproject.discourse.group/t/how-to-write-stresses-to-a-xdmf-file/11306/1 "2023-05-26T16:46:56Z")

</div>

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

```auto
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

```auto
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?

---

<div class="post-metadata">

**Author:** ![nav](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nav/32/5225_2.png) [@nav](https://fenicsproject.discourse.group/u/nav)\
**Post date:** [May 26, 2023, 5:17pm UTC](https://fenicsproject.discourse.group/t/how-to-write-stresses-to-a-xdmf-file/11306/2 "2023-05-26T17:17:42Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![aaquib](https://avatars.discourse-cdn.com/v4/letter/a/d2c977/32.png) [@aaquib](https://fenicsproject.discourse.group/u/aaquib)\
**Post date:** [May 26, 2023, 5:34pm UTC](https://fenicsproject.discourse.group/t/how-to-write-stresses-to-a-xdmf-file/11306/3 "2023-05-26T17:34:54Z")

</div>

Sure,

```auto
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!

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [May 26, 2023, 6:18pm UTC](https://fenicsproject.discourse.group/t/how-to-write-stresses-to-a-xdmf-file/11306/4 "2023-05-26T18:18:13Z")

</div>

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](https://jsdokken.com/dolfinx-tutorial/chapter3/em.html)  
E.g.

```python
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](https://jsdokken.com/fenics22-tutorial/heat_eq.html#post-processing-without-projections)  
or  
[https://jsdokken.com/dolfinx-tutorial/chapter2/linearelasticity\_code.html#stress-computation](https://jsdokken.com/dolfinx-tutorial/chapter2/linearelasticity_code.html#stress-computation)

---

<div class="post-metadata">

**Author:** ![aaquib](https://avatars.discourse-cdn.com/v4/letter/a/d2c977/32.png) [@aaquib](https://fenicsproject.discourse.group/u/aaquib)\
**Post date:** [May 27, 2023, 3:04pm UTC](https://fenicsproject.discourse.group/t/how-to-write-stresses-to-a-xdmf-file/11306/5 "2023-05-27T15:04:25Z")

</div>

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

```auto
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

```auto
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

```

---

<div class="post-metadata">

**Author:** ![nav](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nav/32/5225_2.png) [@nav](https://fenicsproject.discourse.group/u/nav)\
**Post date:** [May 27, 2023, 3:20pm UTC](https://fenicsproject.discourse.group/t/how-to-write-stresses-to-a-xdmf-file/11306/6 "2023-05-27T15:20:58Z")

</div>

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

Change

> [@aaquib](#):
>
> `W=fem.VectorFunctionSpace(domain,("DG",0))`

to

```python3
W=fem.FunctionSpace(domain,("DG",0))

```

---

<div class="post-metadata">

**Author:** ![aaquib](https://avatars.discourse-cdn.com/v4/letter/a/d2c977/32.png) [@aaquib](https://fenicsproject.discourse.group/u/aaquib)\
**Post date:** [May 27, 2023, 3:37pm UTC](https://fenicsproject.discourse.group/t/how-to-write-stresses-to-a-xdmf-file/11306/7 "2023-05-27T15:37:30Z")

</div>

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

```auto
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)

```