Hi,
I’d like to compute a dot product of a tensor(matrix) field and a vector field as below
\mathbf v \left(\mathbf x\right) = \mathbf C \left(\mathbf x\right) \mathbf u \left(\mathbf x\right)
But for some reason, the matrix component numbering is way different from that I thought. I got,
import numpy as np
from dolfin import *
ndim = 2
mesh = UnitSquareMesh(3, 3)
VS = VectorFunctionSpace(mesh, 'CG', 1, dim=ndim)
TS = TensorFunctionSpace(mesh, 'CG', 1, shape=(ndim, ndim))
coo = VS.tabulate_dof_coordinates()[::ndim]
nvrt = coo.shape[0]
u = Function(VS)
C = Function(TS)
u0 = np.random.rand(nvrt, ndim)
C0 = np.random.rand(nvrt, ndim, ndim)
C0 = .5*(C0 + np.transpose(C0, axes=(0,2,1)))
u.vector().set_local(u0.flatten())
C.vector().set_local(C0.flatten())
v = project(dot(C,u), VS)
v0 = np.einsum('ijk,ik->ij', C0, u0)
I though v and v0 are identical, but
print(v(coo[0]))
print(v0[0])
it gives
[0.79744342 1.25202656]
[0.77438778 1.26000916]
In think set_local() might be wrong in here because the dof-component mapping was wrong.
For 2-dimensional vector field, I know that,
\mathbf u=\begin{bmatrix}u_0^{(0)}&u_1^{(0)}&\cdots&u_0^{(n)}&u_1^{(n)}\end{bmatrix}^\mathsf T
For a matrix field, I though its identical, like,
\mathbf C = \begin{bmatrix}C_{00}^{(0)}&C_{01}^{(0)}&C_{10}^{(0)}&C_{11}^{(0)}&\cdots&C_{00}^{(n)}&C_{01}^{(n)}&C_{10}^{(n)}&C_{11}^{(n)}\end{bmatrix}^\mathsf T
However, according to my test code, it is not.
How can I set local values of a tensor(matrix) field as I intended to?