Hello everyone, I performed a 3D simulation of the Stokes flow against an obstacle. The outputs are represented by the fields (u, p) in R^3 and R respectively, in accordance with mathematics. Now I would like to integrate the symmetric part of the gradient of the u field (called epsilon(u) in many examples) on the volumetric domain. I don’t really know how can I set up the problem. Many thanks

Do you mean that you want to compute \int_\Omega \frac{1}{2}(\nabla u+\nabla^T u) dV ? If so, the resulting object is then a 2nd rank tensor and you must compute each component individually with

```
assemble(epsilon(u)[i, j]*dx)
```

for given i and j indices

I want to compute this integral, where D stands for the sym(grad) operator.

and the indexes are i,j = 1,2,3.

How about the following

```
def symgrad(u, i):
return sym(grad(u[i]))
K = np.zeros((3,3))
for i in range(3):
for j in range(3):
K[i,j] = assemble(2*mu*ufl.inner(symgrad(u, i), symgrad(u, j))*dx)
```

thank you dokken, but i get the following error

ufl.log.UFLException: Symmetric part of tensor with rank != 2 is undefined.

Please provide a minimal code reproducing your error message.

I’ve post the meshed domain and the code in the following free-access folder.

https://drive.google.com/drive/folders/1_pOWo9-23KmrvvI6uviZc6Hwmsn-pCQT?usp=sharing

For the sake of reproducibility and the maintainability of the forum, I would appreciate if you could make the code a minimal code example (using a built-in mesh), that you present directly on the forum (rather than a google drive link that might change/disappear) in the future.

You can post code directly on the forum using markdown syntax with 3x` encapsulation, as shown below:

```
def test():
print("HERE")
test()
```

```
w = Function(W)
solve(a == L, w, bcs)
u, p = w.split()
def symgrad(u, i):
return sym(grad(u))
K = np.zeros((3,3))
for i in range(3):
for j in range(3):
K[i,j] = assemble(2*mu*inner(symgrad(u, i), symgrad(u, j))*dx)
print("K_ij =", K[i,j])
```

This is the minimal code I think, what i would to do is to use the output u from the solver and then integrate over dx the formula above. You think this way is appropriate for the forum?

A minimal code needs the definition of the function spaces, and should be executable by anyone. The code below is an example of this:

```
from dolfin import *
import numpy as np
mesh = UnitCubeMesh(10, 10, 10)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
B = FiniteElement("Bubble", mesh.ufl_cell(), mesh.topology().dim() + 1)
V = VectorElement(NodalEnrichedElement(P1, B))
Q = P1
W = FunctionSpace(mesh, V * Q)
class Everywhere(SubDomain):
def inside(self, x, on_boundary):
return True
w = Function(W)
x = SpatialCoordinate(mesh)
W0 = W.sub(0).collapse()
u_ = project(as_vector((x[0], 2*x[1], 3*x[2]**2)), W0)
fa = FunctionAssigner(W.sub(0), W0)
fa.assign(w.sub(0), u_)
u, p = w.split()
def symgrad(u, i):
return sym(grad(u[i]))
mu = 1.49
K = np.zeros((3, 3))
for i in range(3):
for j in range(3):
K[i, j] = assemble(2*mu*inner(symgrad(u, i), symgrad(u, j))*dx)
print("K_ij =", i, j, K[i, j])
```

which produces your error message.

This is because you are trying to take the symmetric gradient of `u[i]`

, which is a scalar value.

The gradient of `u[i]`

is a vector and taking the symmetric part of a vector does not make sense.

Thus you need to specify what the `i`

and `j`

index is supposed to work on.

Yes, I think this is the right way to formulate the integral:

```
def symgrad(u):
return sym(grad(u))
K = np.zeros((3,3))
for i in range(3):
for j in range(3):
K[i,j] = assemble(2*mu*inner(symgrad(u)[i,:], symgrad(u)[j,:])*dx)
print("K_ij =", K[i,j])
```

do you agree?

It depends on what your initial integral

means.

If D\mathbf{u}_i means the ith row of D\mathbf{u} you are correct in your interpretation above.