Magnetostrain tensor product matrix

Hi,

When solving micromagnetics including elasticity, it is necessary to build a quadratic form of the kind
m1m1 -1/3 m1m2 m1m3
m2
m1 m2m2 -1/3 m2m3
m3m1 m3m2 m3*m3 -1/3

from a given magnetisation with components (m1,m2,m3). I had originally planned on building this using a GridFunction, but I am using a first order FE space, so it is inappropriate for this quadratic function when used in quadrature (building the linear system).

The documentation isn’t massively clear, but I am thinking I should build it using a CoefficientFunction as follows.
Code:

[code]fes_mag = VectorH1(mesh, order=1)
mag_grid_func = GridFunction(fes_mag)

give_random_magnetisation(mag_grid_func) # simply a function that fills each DoF with a random value such that the magnitude at each node is 1.
m1, m2, m3 = mag_grid_func.components
mymatrix = CoefficientFunction((m1m1 - 1/3,m1m2, m1m3, m2m1, m2m2 - 1/3, m2m3, m3m1, m3m2, m3*m3 - 1/3), dims=(3, 3))[/code]

If there is a better way, or built-in method for doing this I would appreciate it.

EDIT: I have attempted to upload a file of a basic code, but the add files button seems non-functional. It reports “Error Unable to get properties for the image.” despite me uploading either .py and .ipynb or a .zip file containing them.

mymatrix = OuterProduct (mag_grid_func, mag_grid_func) - 1/3*Id(3)
print (mymatrix)

best, Joachim