Assembling and decompose SPD tensor

Question

Here's my problem:
I have to assemble a tensor (SPD) in TensorFunctionSpace(mesh, 'DG', 0) by calculating each entry.

card=Th.init(2)
idx=0
T=Function(TensorFunctionSpace(mesh, 'DG', 0))
v=numpy.empty((1, 4*card))

while idx<=card-1: 
      #calculation for every CELL and assembling in v
      v[0, idx] = ...
      v[0, card+idx] = ...
      v[0, 2*card+idx] = ...
      v[0, 3*card+idx] = ...


T.vector()[:] = #array previously calculated

Now, my questions are:

Is this correct or am I messing up with cells and dofs positions? How could I solve it?

The tensor is supposed to be SPD, but when I use numpy.linalg.eig I get negative eigenvalues. Any hint?

Answer 1

Hi, see it the following helps

from dolfin import *
import numpy as np

mesh = UnitSquareMesh(10, 10)
T = TensorFunctionSpace(mesh, 'DG', 0)
t = Function(T)
t_values = t.vector().array()

dofmap = T.dofmap()
# Build values of tensor
for cell in cells(mesh):
    # Get dof of each tensor component on the cell
    dofs = dofmap.cell_dofs(cell.index())
    # Or
    # dofs = [T.sub(i).dofmap().cell_dofs(cell.index())[0] for i in range(4)]
    dof0, dof1, dof2, dof3 = dofs

    # We make tensor values depend on position
    x, y = cell.midpoint().x(), cell.midpoint().y()
    t_values[dof0] = x**2+3  # [0, 0]
    t_values[dof1] = 1       # [0, 1]  
    t_values[dof2] = 1       # [1, 0]
    t_values[dof3] = y**2+2  # [1, 1]

t.vector().set_local(t_values)

# Local check for SPD
x = [0.22, 0.24]
mat = t(x).reshape((2, 2))
assert np.all(np.linalg.eigvals(mat) > 0), 'Not all eigs > 0'
assert near(mat[0, 1], mat[1, 0]), 'Not symmetric' 

Comments

Thank you, the code works but I still have problems with the eigenvalues. Probably I must check how I calculate the entries of the tensor.