Suppose Z is a fully symmetric fourth-order tensor, which is both isotropic and isochoric (volume preserving). I need to build a coefficient function of the form

2 ZC[ε(u) - ε_m(M)] M

where Z is the aforementioned tensor, C is the isotropic fourth-order elastic tensor (this can be stolen directly from here), ε is the usual strain Sym(Grad(u)), ε_M is the magnetostrain mentioned in my previous post here which is given by Z(M⊗M), and |M|^2 = 1. This needs to be used in quadrature, so projecting into my first-order finite element space is not appropriate (although this would make it easier as I could move to numpy/scipy).

My question is how to best calculate Z(MatrixCoefficientFunction). I am thinking that I should be able to simply use the same representation that is used for the C tensor (2µ ε + λTr(ε)δ_ij), with the coefficients chosen so that when Z is applied to M⊗M you get ε_M. A suitable one would be ε-Tr(ε)δ_ij, after a suitable scaling A, implemented as say

A(stresslike - Trace(stresslike)*Id(3))

Otherwise, I will have to perform a multiplication between Z and C[ε(u) - ε_m(M)], and I am not sure how to even represent fourth order tensors within ngsolve, let alone multiply a matrix coefficient function by one.

Yes. In general, it will look like Z{ijkl} σ{kl} m_{j}, where the stress is constructed as
σ_{kl} = C_{klpq}[ε_pq(u) - ε^M _pq(m)]

where ε_pq is the symmetric part of the gradient of u, and ε^M pq(m) is [b]Z/b = Z{pqab} (m⊗m){ab}, a function of m. Together, this is a large mess Z{ijkl} C{klpq}[ε_pq(u) - Z{pqab}(m⊗m){ab}] m{j}.

Here {i,j,k,l,p,q,a,b} are indices between 1 and 3. Z is a fourth-order minorly symmetric tensor (Z{ijkl} = Z{jikl} = Z_{ijlk}), and C is the usual fourth-order elastic tensor.

There should be a free index “i” left over, for use in an inner product of the form <Z{ijkl} σ{kl} m_{j}, v_i> during assembly.

Both Z and C can either be stored as Voigt (6x6) arrays, or in their full 3x3x3x3 form, it shouldn’t matter as long as correct transformations are made.

Yes. In general, it will look like Z{ijkl} σ{kl} m_{j}, where the stress is constructed as
σ_{kl} = C_{klpq}[ε_pq(u) - ε^M _pq(m)]

where ε_pq is the symmetric part of the gradient of u, and ε^M pq(m) is [b]Z/b = Z{pqab} (m⊗m){ab}= Z{pqab} m_{a}m_{b}, a function of m. Together, this is a large mess Z{ijkl} C{klpq}[ε_pq(u) - Z{pqab}m{a}m_{b}] m_{j}.

Here {i,j,k,l,p,q,a,b} are indices between 1 and 3. Z is a fourth-order minorly symmetric tensor (Z{ijkl} = Z{jikl} = Z_{ijlk}), and C is the usual fourth-order elastic tensor.

There should be a free index “i” left over, for use in an inner product of the form <Z{ijkl} σ{kl} m_{j}, v_i> during assembly.

Both Z and C can either be stored as Voigt (6x6) arrays, or in their full 3x3x3x3 form, it shouldn’t matter as long as correct transformations are made. I have attached the Voigt notation for the isotropic, isochoric case of Z. Some more details of Z can be found in “Tensor representation of magnetostriction for all crystal classes, Federico 2018” if needed.

The rest you can do using the usual matrix/vector operations.

Internally, tensors are stored as vectors, where the rightmost index is always the innermost index. With Reshape, you use the same underlying vector, but access it with a different tensor shape.