 Integration of a gradient of function

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

1 Like

I want to compute this integral, where D stands for the sym(grad) operator. and the indexes are i,j = 1,2,3.

K = np.zeros((3,3))
for i in range(3):
for j in range(3):

thank you dokken, but i get the following error
ufl.log.UFLException: Symmetric part of tensor with rank != 2 is undefined.

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

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()
1 Like
w = Function(W)
solve(a == L, w, bcs)
u, p = w.split()

K = np.zeros((3,3))
for i in range(3):
for j in range(3):
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, 2*x, 3*x**2)), W0)
fa = FunctionAssigner(W.sub(0), W0)
fa.assign(w.sub(0), u_)

u, p = w.split()

mu = 1.49
K = np.zeros((3, 3))
for i in range(3):
for j in range(3):
print("K_ij =", i, j, K[i, j])

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.

2 Likes

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

K = np.zeros((3,3))
for i in range(3):
for j in range(3):