Implementation for isotropic, isochoric tensor multiplication.

Hi,

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


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.

Can you write down the expression you want to build using indices ? 

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 Z(m⊗m) = 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.

You can reshape tensors to matrices or vectors. To define

Zsigma_ij = Z_ijkl sigma_kl

you can do

(Z.Reshape( (9,9) ) * sigma.Reshape( (9,) ) . Reshape( (3,3) )

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.

Joachim

 

Can you clarify what you mean by "Internally, tensors are stored as vectors, where the rightmost index is always the innermost index."?

Thanks.