Bug in symmetric TensorFunctionSpace?

Question

I am trying to initilialise a symmetric TensorFunctionSpace,

>>> from dolfin import *
>>> mesh = RectangleMesh(Point(0.,0.),Point(2.,1.),20,10)
>>> T = interpolate(Constant(((3.,1.),(1.,4.))), TensorFunctionSpace(mesh,'CG',1,symmetry=True))
>>> print (T.vector().array()==4).sum()
0
>>> T = interpolate(Constant(((3.,1.),(4.,1.))), TensorFunctionSpace(mesh,'CG',1,symmetry=True))
>>> print (T.vector().array()==4).sum()
231

Is this behaviour expected? Should FEniCS not read the last diagonal component as the fourth element of the Constant?

EDIT: To follow up, plotting looks OK, but evaluation worries me:

>>> T(0.5,0.5)
array([ 3.,  1.,  4.,  0.])
>>> plot(T[0,0])
>>> plot(T[0,1])
>>> plot(T[1,0])
>>> plot(T[1,1])
>>> interactive()

Answer 1

Hi, this is most likely a bug. Below is an example showing that a symmetric TensorFunctionSpace cannot exactly represent a constant symmetric tensor field. The last expression also shows why your plot looks OK.

from dolfin import *
import numpy as np

A = np.array([[10, 2],
              [1, -3]])
Sym = np.array([[1, 0.5],
                [0.5, 2]])

mesh = UnitTriangleMesh()
V = TensorFunctionSpace(mesh, 'CG', 1)
S = TensorFunctionSpace(mesh, 'CG', 1, symmetry=True)

# See if T can be represented exactly in S
test = lambda T, S: np.allclose(interpolate(Constant(T), S)(0.5, 0.5).reshape((2, 2))-T,
                                1E-13)
# Sanity
assert test(A, V)
assert test(Sym, V)
# Correcly comes out as False - can't represent not symm tensor. 
print test(A, S)           
print test(Sym, S)         # BUG: False despite symmetry. True only if Sym[1, 1]
Sym[1, 1] = 0
print test(Sym, S)      

Is it fine if I report this issue on bitbucket?

Comments

Yes, it is fine, if you report it - thank you.

I do not understand the last expression. If the Function looks correct in the plot, will it behave correctly, if I use it in a bilinear form?

---

The example shows that symmetric 2x2 tensor fields are not interpolated correctly in geneneral. It appears that there is an additional requirement of having zero in (1, 1) position. If you can guarantee that these constraints are met in your application you should be fine. Anyway, I would not recommend this.

---

Reported as Issue660.