Integrating vector fields

I would like to integrate vector fields over a given domain. The docs of the assemble
method seem to suggest the following:

from dolfin import *

mesh = UnitSquareMesh(32, 32)
V = VectorFunctionSpace(mesh, "Lagrange", 1)
f = Function(V)

F = PETScVector()                         
assemble(f*dx, tensor = F)

However, this returns an error message:
UFLException: Can only integrate scalar expressions. The integrand is a tensor expression with value shape (2,) and free indices with labels ().

What is the problem here?

In this example you can not directly multiply f by dx as f has been defined on a vectorial function space. You may want to incorporate an inner product by defining a test function. Then you will have a scalar variable as an integrand:

from dolfin import *
mesh = UnitSquareMesh(32, 32)
V = VectorFunctionSpace(mesh, "Lagrange", 1)
f = Function(V)
v = TestFunction(V)
F = PETScVector()
assemble(inner(f,v)*dx, tensor = F)

Thank you! However, this does not seem to be doing what I want to achieve. The result should contain two real numbers, namely the integrals of each component of the vector field. That is,

from dolfin import *
mesh = UnitSquareMesh(32, 32)
V = VectorFunctionSpace(mesh, "Lagrange", 1)
f = Constant((0,1))
f = project(f, V)
v = TestFunction(V)
F = PETScVector()
assemble(inner(f,v)*dx, tensor = F)

should return the vector (0,1), but this is not the case.

Does it work for your purposes to obtain the components as a Python tuple? That can be done by calling assemble for each component of f as follows:

from dolfin import *
mesh = UnitSquareMesh(32, 32)
V = VectorFunctionSpace(mesh, "Lagrange", 1)
f = Function(V)
f.interpolate(Constant((0,1)))
v = (assemble(f[0]*dx),assemble(f[1]*dx))
print(v)

Yes, this would be a workaround. However, this can become quite tedious for higher order tensors. I thought that the tensor keyword of the assemble function was meant for this. Is this not correct?

The term “tensor” in the context of the assemble() function refers to tensors defined on finite element function spaces (e.g., the right-hand-side vector or left-hand-side matrix from the linear system that you would solve) rather than the typical usage in mechanics, referring to tensors defined on the tangent space to your PDE domain (e.g., velocity, stress, etc.). The rank of tensor you get depends on how many trial/test functions are in the form you are assembling (so you get a scalar if there are no test/trial functions, a vector if there is test function, and a matrix if there is a test and a trial function).

2 Likes

Thank you for clarifying this!

You nee to use a “Real” element (i.e. Lagrange multiplier, defined over the whole mesh) to build your test functions:

from dolfin import *
mesh = UnitSquareMesh(32, 32)
V = VectorFunctionSpace(mesh, "Lagrange", 1)
VR = VectorFunctionSpace(mesh, "R", 0)
f = Constant((0,1))
f = project(f, V)
v = TestFunction(VR)
F = PETScVector()
assemble(inner(f,v)*dx, tensor = F)
print(F[:])