Integration of a gradient of function

I/O
1

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

2

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
3

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

4

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

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

8

Please provide a minimal code reproducing your error message.

9

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

10

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
11 
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?

12

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.

2 Likes
13

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?

14

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.