Outer product evaluation

Linear Algebra
1

Hi
I was wondering how it is possible to extract the values of a tensor produced by outer product of two vectors. Let’s say we have 2 unit vectors and the outer product of these vectors creates a tensor. Now I want to extract the values of this tensor:

from fenics import *

mesh = BoxMesh(Point(0, 0, 0), Point(1, 1, 1), 1, 1, 1)

A = as_vector([[1], [0],[0]])
B = as_vector([[0], [1],[0]])

C = outer(A,B)

T = TensorFunctionSpace(mesh, 'CG', 1)
E = project(C, T)

#Tensor components extraction
print(E.vector().get_local())

It returns Shapes do not match error. However it works if I change outer to inner and replace the TensorFunctionSpace with FunctionSpace but I want to know why it does not work for outer product.

2

Hi,

The way in which you are defining the vectors A and B generates “vectors” of second order. Therefore the outer product of A and B gives a fourth order tensor that doesn’t match the order of the TensorFunctionSpace, which is a 3x3 second order tensor. The following example should work

from fenics import *

mesh = BoxMesh(Point(0, 0, 0), Point(1, 1, 1), 1, 1, 1)

A = as_vector([1, 0, 0])
B = as_vector([0, 1, 0])

C = outer(A,B)

# For sake of clarity, always is a good idea to define
# the shape of the tensor function space
T = TensorFunctionSpace(mesh, 'CG', 1, shape=(3,3))
E = project(C, T)

To get the i-th component component of the tensor as a dolfin Function object you can use E.sub(i). If you want the dof value you should consider using dofmaps (see here for instance).

3 Likes
3

Thanks for your clarification.

4

One more question:
By definition, the inner product of a vector in a rank 2 tensor should be a vector. For example using numpy:

import numpy as np

a = [[1, 5]]
b = [[1, 2], [3, 4]]

print(np.inner(a,b))

and it works fine. However same implementation in FEniCS:

from fenics import *

A = as_vector([1, 5])
B = as_tensor([[1, 2], [3, 4]])

C = inner(A,B)

Does not work and returns shape mismatch error. If I use np.inner insetead of inner then it works. It is a bit strange. Any idea why it is like that?

5

Hi, from UFL definition

The inner product is a contraction over all axes of a and b, that is the sum of all component-wise products. The operands must have exactly the same dimensions. For two vectors it is equivalent to the dot product.

whereas in Numpy inner performs a sum product over the last axes.

What you want is probably simply UFL dot

6

Thanks. So the np.inner is not predefined in FEniCS. I just want to calculate and compare the inner and dot product of two same rank tensors. Here is an example:

from fenics import *

n = 1
mesh = BoxMesh(Point(0, 0, 0), Point(1, 1, 1), n, n, n)

A = as_tensor([[2, 1],[4,3]])
B = as_tensor([[1, 2], [3, 4]])

DOT = dot(A,B)
INNER = inner(A,B)  # Up to this line there is no error

SPACE = TensorFunctionSpace(mesh, 'CG', 1)

DOT_PROJECT = project(DOT, SPACE)  #Error comes up here
INNER_PROJECT = project(INNER, SPACE)

print(DOT_PROJECT.vector().get_local()[0])
print(INNER_PROJECT.vector().get_local()[0])

What is the right way to do it?

7

Hello,
two points:

  • you use a 3D mesh and don’t specify the shape of the TensorFunctionSpace, by default it will deduce the dimension from that of the mesh. SPACE will therefore be 3x3 tensors but DOT is only 2x2 hence your first error.
  • INNER is a scalar since it is A_{ij}B_{ij}, you cannot project it on a TensorFunctionSpace
    I have trouble understanding why you are attempting to compare two different things
1 Like
8

Thanks. The difference in definition of the inner product in FEniCS and numpy got me confused. Now it is clear.