How does one write general tensor contractions in the Python-based finite element package Nutils? For example, how does one write the contraction of a fourth-order elasticity tensor $\boldsymbol{C}$ with a second-order strain tensor $\boldsymbol{\varepsilon}$, to yield the second-order stress tensor $\boldsymbol{\sigma}=\boldsymbol{C}:\boldsymbol{\varepsilon}$?


For high dimensional tensor manipulations you should avoid using numpy.dot, which is very limited in scope, and instead use a combination of multiplication and summation. Although this would cost performance when used directly on Numpy arrays this is not the case in Nutils which uses a delayed execution model, and in fact uses dot or einsum under the hood.

In the following I use _ = numpy.newaxis as is customary in Nutils:

Inner product of two vectors

For $\mathbf{a} \in \mathbb{R}^n$ and $\mathbf{b} \in \mathbb{R}^n, c = \mathbf{a} \cdot \mathbf{b} \in \mathbb{R}$ (Einstein notation: $c = a_i b_i$)

c = numpy.dot(a,b) # numpy only
c = (a*b).sum(-1) # numpy/nutils
c = a['i']*b['i'] # nutils only

Outer product of two vectors

For $\mathbf{a} \in \mathbb{R}^n$ and $\mathbf{b} \in \mathbb{R}^m$, $\mathbf{C} = \mathbf{a} \otimes \mathbf{b} \in \mathbb{R}^{n \times m}$ (Einstein notation: $C_{ij} = a_i b_j$)

C = a[:,_]*b[_,:] # numpy/nutils
C = a['i']*b['j'] # nutils only

Matrix-vector product

For $\mathbf{A} \in \mathbb{R}^{m \times n}$ and $\mathbf{b} \in \mathbb{R}^n$, $\mathbf{c} = \mathbf{A}\mathbf{b} \in \mathbb{R}^{m}$ (Einstein notation: $c_{i} = A_{ij} b_j$)

c = numpy.dot(A,b) # numpy only
c = A[:,:]*b[_,:]).sum(-1) # numpy/nutils
c = a['ij']*b['j'] # nutils only

Double contraction of 4th and 2nd dimensional tensors

For $\mathbf{A} \in \mathbb{R}^{q \times r \times s \times t}$ and $\mathbf{B} \in \mathbb{R}^{s \times t}$, $\mathbf{C} = \mathbf{A}:\mathbf{B} \in \mathbb{R}^{q \times r}$ (Einstein notation: $C_{ij} = A_{ijkl} B_{kl}$)

C = A[:,:,:,:]*B[_,_,:,:]).sum([-1,-2]) # numpy/nutils
C = A['ijkl']*B['kl'] # nutils only

Not the answer you're looking for? Browse other questions tagged or ask your own question.