Tensor Arithmetics edit page

MTEX offers powerful functionalities to calculate with tensors and lists of tensors without the need of many nested loops.

Basic algebra

First of all, all basic operations like *, .*, +, -, .^ known from linear algebra work also on lists of tensors.

T1 = tensor.rand('rank',2);
T2 = tensor.rand('rank',2);

% addition and multiplication
T = T1 + 2 * T2;

% point-wise product
T = T1 .* T2;

Tensor Products

Tensor product are the canonical way how tensors interact with each other. As an example consider a rank 4 stiffness tensor

C = stiffnessTensor.load(fullfile(mtexDataPath,'tensor','Olivine1997PC.GPa'))
C = stiffnessTensor (xyz)
  unit: GPa              
  rank: 4 (3 x 3 x 3 x 3)
 
  tensor in Voigt matrix representation:
 320.5  68.2  71.6     0     0     0
  68.2 196.5  76.8     0     0     0
  71.6  76.8 233.5     0     0     0
     0     0     0    64     0     0
     0     0     0     0    77     0
     0     0     0     0     0  78.7

Then by Hooks law the stiffness tensor acts on a strain tensor, e.g.,

eps = strainTensor(diag([1 0 -1]))
eps = strainTensor (xyz)
  type: Lagrange 
  rank: 2 (3 x 3)
 
  1  0  0
  0  0  0
  0  0 -1

according to the formula

\[\sigma_{ij} =\sum_{k,l} C_{ijkl} \epsilon_{kl}\]

and turns it into the stress tensor \(\sigma\). In MTEX such tensor products can be computed in its most general form by the command EinsteinSum.

sigma = EinsteinSum(C,[1 2 -1 -2],eps,[-1 -2])
sigma = tensor (xyz)
  rank: 2 (3 x 3)
 
  248.9      0      0
      0  -8.65      0
      0      0 -161.9

here the negative numbers indicate the indices which are summed up. Each pair of equal negative numbers corresponds to one sum. The positive numbers indicate the order of the dimensions of the resulting tensor. Accordingly we can compute the outer product

\[ (a \otimes b)_{ij} = a_i b_j \]

between two rank one tensors

a = tensor([1;2;3],'rank',1);
b = tensor([0;2;1],'rank',1);

by the command

EinsteinSum(a,1,b,2)
ans = tensor (xyz)
  rank: 2 (3 x 3)
 
 0 2 1
 0 4 2
 0 6 3

and the inner product

\[ a \cdot b = \sum_i a_i b_i \]

by

EinsteinSum(a,-1,b,-1)
ans =
     7

As a final example we consider the linear compressibility in a certain direction v which can be computed by the formula

\[ c = \sum_{i,j,k} S_{ijkk} v_i v_j \]

where \(C = S^{-1}\) is the inverse of the compliance tensor, i.e. the stiffness tensor

v = xvector;
S = inv(C)
c = EinsteinSum(C,[-1 -2 -3 -3],v,-1,v,-2)
S = complianceTensor (xyz)
  unit            : 1/GPa            
  rank            : 4 (3 x 3 x 3 x 3)
  doubleConvention: true             
 
  tensor in Voigt matrix representation: *10^-4
  34.85  -9.08   -7.7      0      0      0
  -9.08  60.76  -17.2      0      0      0
   -7.7  -17.2  50.85      0      0      0
      0      0      0 156.25      0      0
      0      0      0      0 129.87      0
      0      0      0      0      0 127.06
c =
  460.2500

Here we used the command inv to compute the inverse of any rank 2 or rank 4 tensor. There are shortcuts in MTEX for specific tensor products. E.g. the relation between stress and strain can be more compactly written as a double dot product 

C * eps
C : eps
ans = stressTensor (xyz)
  rank: 2 (3 x 3)
 
  248.9      0      0
      0  -8.65      0
      0      0 -161.9
 
ans = stressTensor (xyz)
  rank: 2 (3 x 3)
 
  248.9      0      0
      0  -8.65      0
      0      0 -161.9

The double dot product between two rank two tensors is essentially their inner product and can be equivalently computed from the trace of their matrix product

T1 : T2
trace(T1 * T2')
trace(T1' * T2)
ans =
    1.7334
ans =
    1.7334
ans =
    1.7334

Determinant

For rank two tensors we can compute the determinant of the tensor by the command det 

det(T1)
ans =
    0.4989

Rotating a tensor

Rotation a tensor is done by the command rotate. Let's define a rotation

r = rotation.byEuler(45*degree,0*degree,0*degree);

Then the rotated tensor is given by

Trot = rotate(T1,r);
plot(Trot)

Here is another example from Nye (Physical Properties of Crystals, p.120-121) for a third-rank tensor

P = [ 0 0 0 .17 0   0;
      0 0 0 0   .17 0;
      0 0 0 0   0   5.17]*10^-11;

T = tensor(P,'rank',3,'propertyname','piezoelectric modulus')

r = rotation.byAxisAngle(zvector,-45*degree);
T = rotate(T,r)
T = tensor (xyz)
  propertyname: piezoelectric modulus
  rank        : 3 (3 x 3 x 3)        
 
  tensor in compact matrix form: *10^-12
    0    0    0  1.7    0    0
    0    0    0    0  1.7    0
    0    0    0    0    0 51.7
 
T = tensor (xyz)
  propertyname: piezoelectric modulus
  rank        : 3 (3 x 3 x 3)        
 
  tensor in compact matrix form: *10^-12
     0     0     0     0   1.7     0
     0     0     0  -1.7     0     0
  51.7 -51.7     0     0     0     0