Jon Wright, 26 April (-> 1st May -> +AH 19th May) 2020
The observed data from a 3DXRD experiment are the scattering vectors in the sample frame, the "g-vectors". If we succeed to index these using orientation matrices, (UB), we assign hkl indices to spots. The matrix inverse of (UB) gives the real space unit cell vectors, called $\mathbf{L}$ here (for lattice, it is the ubi from ImageD11). In the ImageD11 code there is often confusion between row and column vectors as the code usually applies transformations from the left and does not distinguish row and column vectors. If you substitute h,k,l as unit vectors in (1) you can identify the unit cell vectors in reciprocal space:
\begin{equation} \begin{pmatrix} g_x \\ g_y \\ g_z \end{pmatrix} = \mathbf{L^{-1}} \begin{pmatrix} h \\ k \\ l \end{pmatrix} = (\mathbf{UB}) \begin{pmatrix} h \\ k \\ l \end{pmatrix} = \begin{pmatrix} \begin{pmatrix} x_a^* \\ y_a^* \\ z_a^* \end{pmatrix} \begin{pmatrix} x_b^* \\ y_b^* \\ z_b^* \end{pmatrix} \begin{pmatrix} x_c^* \\ y_c^* \\ z_c^* \end{pmatrix} \end{pmatrix} \begin{pmatrix} h \\ k \\ l \end{pmatrix} \tag{1} \end{equation}\begin{equation} \begin{pmatrix} h \\ k \\ l \end{pmatrix} = \mathbf{L} \begin{pmatrix} g_x \\ g_y \\ g_z \end{pmatrix} = ( \mathbf{UB} )^{-1}\begin{pmatrix} g_x \\ g_y \\ g_z \end{pmatrix} = \begin{pmatrix} (x_a & y_a & z_a ) \\ (x_b & y_b & z_b ) \\ (x_c & y_c & z_c ) \end{pmatrix} \begin{pmatrix} g_x \\ g_y \\ g_z \end{pmatrix} \tag{2} \end{equation}The real space lattice as row vectors is not quite so obvious. Instead you need to multiply from the left with a row vector in crystal space $(c_x,c_y,c_z)$ as unit vector and compute a column vector in laboratory space, $(l_x,l_y,l_z)$ as in (3):
\begin{equation} (c_x, c_y, c_z) \begin{pmatrix} (x_a & y_a & z_a ) \\ (x_b & y_b & z_b ) \\ (x_c & y_c & z_c ) \end{pmatrix} = \begin{pmatrix} l_x \\ l_y \\ l_z \end{pmatrix} \tag{3} \end{equation}The difference of row versus column vectors and coming from the left or the right of a matrix has been the cause of many years of confusion. Many thanks to Axel Henningsson for pointing out some errors here. A quick note to self: when we multiply two matrices A, B with shapes (rows,cols):(M,K) and (K,N), then the K has to be the same and $ C = A.B $ with shapes $ C[m,n] = A[m,k].B[k,n] $.
If a crystal is strained then we look how the atoms move compared to a reference which was strain free. If we are sitting on some atom in the initial state, the positions of the atoms in the next unit cells are given by a matrix $\mathbf{L_0}$. After strain the neighbor positions will now be given by the 3 lattice vectors:
$ \mathbf{a} = \mathbf{F . a_0 } $
$ \mathbf{b} = \mathbf{F . b_0 } $
$ \mathbf{c} = \mathbf{F . c_0 } $
$ \mathbf{L^T} = \mathbf{F . L_0^T }$
[To be confirmed: second attempt: I think this is the deformation gradient tensor from finite strain theory, but neglecting any rigid body translations].
This is easier for most people to read in terms of column vectors with $\mathbf{A = L^T}$ and $\mathbf{A = F.A_0}$.
Assuming this is the $\mathbf{F}$ of finite strain theory then everything else falls beautifully into place. If we have a reference unit cell, in any orientation, then we can compute $ \mathbf{F} $ which includes both a change in orientation and strain:
$ \mathbf{L^T} = \mathbf{F.L_0^T } $
$ \mathbf{F} = \mathbf{ L^T . L_0^{-T}} = \mathbf{ (U.B)^{-T} (U_0.B_0)^{T}} $
$ \mathbf{F} = \mathbf{ U^{-T}.B^{-T}.B_0^{T}.U_0^{T}} $
We see that when we form the products $\mathbf{F^T.F}$ and $\mathbf{F.F^T}$ we are going to rotate the strain into the laboratory or reference axes and annihilate the central rotation matrices.
Note that we can use any axis choice that we like, provided that we are consistent for $\mathbf{L}$ and $\mathbf{L_0}$. There is no requirement here for a selection of a set of orthogonal axes. If we change our definition of $\mathbf{U_0}$ then we may change $\mathbf{F}$, but it will come out in the wash.
The code to compute $\mathbf{F}$ from an ImageD11 grain and reference pair of ubi matrices is in the next cell:
In [1]:
import numpy as np
np.random.seed(4567) # prime, try to get reproducible numbers
def deformation_gradient_tensor(ubi, ubi_0):
"""
Computes the Deformation Gradient Tensor F such that :
ubi = np.dot( F, ubi_0 )
"""
F = np.dot( ubi.T, np.linalg.inv( ubi_0 ).T )
return F
def test_F():
ubi1 = np.random.random((3,3))
ubi0 = np.random.random((3,3))
F = deformation_gradient_tensor( ubi1, ubi0 )
for i in range(3):
assert np.allclose( np.dot(F, ubi0[i]), ubi1[i] )
return "OK, all three axes transformed as they should"
test_F()
Out[1]:
Now that we have $\mathbf{F}$ then it is fairly straightforward to copy and paste various equations from wikipedia. To get the orientation matrix we need to decompose into a rotation R and a strain (V or S). There are two options, first rotate and then strain (RS), or first do the strain and then rotate (VR). In the language of finite strain theory these are the left and right polar decompositions:
$\mathbf{F} = \mathbf{R.S} = \mathbf{V.R} $
I have replaced the U which is on wikipedia for S to avoid collisions with the orientation above. This decomposition can be computed using a singular value decompositon of F, as in the following cell:
In [2]:
def polar_decomposition( F ):
"""
Uses an SVD to compute the polar decomposition of F
F = V.R = R.S
V and S are symmetric deformations in lab and sample respectively
R is a rotation (orthonormal) that relates them
Can be done quicker iteratively, apparently.
Returns V,R,S
"""
w,sing,vh = np.linalg.svd( F )
R = np.dot( w, vh )
S = np.dot( vh.T, np.dot( np.diag(sing), vh ) )
V = np.dot( w , np.dot( np.diag(sing), w.T ) )
return V, R, S
def test_polar():
F = np.random.random((3,3))
V,R,S = polar_decomposition(F)
assert np.allclose( np.dot(V, R), F )
assert np.allclose( np.dot(R, S), F )
return "Polar decomposition looks OK"
test_polar()
Out[2]:
The thing we wanted to compute is the strain. It seems there are an infinite number of different measures we could use in finite strain theory and they should all come down to the same thing for infinitesimal strains. When there are changes in orientation as well as strain, then V is constant and S depends on how/where someone defined the standard orientation. In the past I think we have been using S to get an infinitesimal strain in fable/FitAllB (I am not 100% sure). For the Lagrangian and Eulerian Finite Strains we have:
Lagrangian Finite Strain in the Reference Crystal Axes: \begin{equation} \mathbf{E} = \frac{1}{2}( \mathbf{F^T.F - I } ) = \frac{1}{2} \mathbf{ ( (L^T . L_0^{-T})^T .( L^T . L_0^{-T} ) - I )} = \frac{1}{2} \mathbf{ ( L_0^{-1}.(L.L^T) . L_0^{-T} - I )} \end{equation}
Eulerian Finite Strain, in the Sample x,y,z axes: \begin{equation} \mathbf{e} = \frac{1}{2}( \mathbf{F.F^T - I } ) = \frac{1}{2} \mathbf{ (( L^T . L_0^{-T} ).(L^T . L_0^{-T})^T - I )} = \frac{1}{2} \mathbf{ (L^T .(L_0^{-T}.L_0^{-1}).L - I) } \end{equation}
Note: these come from the infinitesimal strain theory wikipedia page. Both of these equations correspond to a formula where $\epsilon = (d-d_0)/d_0$. If you want other strains you will need to adjust accordingly. The finite strain page has the opposite sign for the Eulerian strain.
In [3]:
def strains_from_deformation(ubi, ubi_0):
"""
Gets the sample (Lagrangian) and lab (Eulerian) system strain tensors
ubi and ubi_0 are 3x3 numpy arrays where the rows are the real
space unit cell vectors.
Also gives back the Polar Decompostion of F = V.R = R.S
Returns E_grain, E_sample, V, R, S
"""
F = deformation_gradient_tensor(ubi, ubi_0)
I = np.eye(3)
E_ref = 0.5 * ( np.dot(F.T,F) - I )
e_lab = 0.5 * ( np.dot(F,F.T) - I )
return E_ref, e_lab
In both of the previous expressions I have substituted for $\mathbf{F}$ and we see the metric tensors have shown up. For the reference state the metric tensor is $( \mathbf{G} = \mathbf{L.L^T} = \mathbf{L^T.L} )$. For the current lab axis state it is the reciprocal metric tensor: $( \mathbf{g_0} = \mathbf{L_0^{-1}.L_0^{-T}} = \mathbf{L_0^{-T}.L_0^{-1}}) $. Finding out how the strains relate to the metric tensor has been a critical problem in Jon's mind. It appears to come down to:
$ \mathbf{ E } = \frac{1}{2} \mathbf{ L_0^{-1}.( G - G_0 ).L_0^{-T} } $
$ \mathbf{e} = \frac{1}{2} \mathbf{ L^{T} .(g_0 - g ). L } $
If I now take the difference in metric tensors of the deformed and undeformed states as being the main observable from the experiment then I can then read finally read some sense into all of these equations. Taking the difference in real or reciprocal space and then how to transform them was not obvious at first sight(!)
The change in metric tensor was always a symmetric matrix that looked like it was going to be strain. Also it is easy to compute in the laboratory axes independently of grains and reference system choices. The problem has been to relate this strain in the crystallographic coordinate system to the sample x/y/z and the crystal orthogonal axes. The strains we found above are just this tensor transformed into the sample or reference axis system.
At last, we attempt to compute these different things in code:
In [4]:
def strains_from_metric(ubi_1, ubi_0):
"""
Strains computed via the shifts in metric tensors
- no matrix decomposition is used
- no choice of reference axes
- no orientation matrix is found
input: ubi and ubi_0 are 3x3 numpy arrays with rows as real
space unit cell vectors.
Returns E_sample E_grain
"""
G0 = np.dot( ubi_0, ubi_0.T ) # reference metric tensor
G = np.dot( ubi_1, ubi_1.T ) # measured metric tensor
ub_0 = np.linalg.inv(ubi_0) # to go to reference axes
# F = np.dot( ubi.T, ubi_0.-T ) = ubi.T, ub0.T
# 2e = FT.F - I = ub_0 ubi ubi.T ub_0.T
# 2e = FT.F - I = ub_0 (G - G0) ub_0.T
E_ref = 0.5*np.dot( np.dot( ub_0, G - G0 ) , ub_0.T ) # into ref space
# 2e = F.FT - I = ubi.T, ub_0.T ub_0 ubi - I
# 2e = F.FT - I = ubi.T( g_0 - g )ubi
dg = np.linalg.inv(G0) - np.linalg.inv(G)
e_lab = 0.5 * np.dot( np.dot( ubi_1.T , dg), ubi_1)
return E_ref, e_lab
def test_strain():
ubi1 = np.random.random((3,3))+np.eye(3)
ubi0 = np.random.random((3,3))+np.eye(3)
Eref1, elab1 = strains_from_deformation(ubi1, ubi0)
Eref2, elab2 = strains_from_metric(ubi1, ubi0)
assert np.allclose(Eref1, Eref2)
assert np.allclose(elab1, elab2)
return "Seems the strains are OK"
test_strain()
Out[4]:
In [5]:
a0 = 4.064 # Aluminium
e12 = a0*0.25 # 25% strain
ubi_0 = np.diag([a0,a0,a0])
ubi = np.array( [[a0, e12, 0],
[e12, a0, 0],
[0, 0, a0]] )
Eref1, elab1 = strains_from_metric( ubi, ubi_0 )
Eref2, elab2 = strains_from_deformation(ubi, ubi_0)
assert np.allclose(Eref1, Eref2)
assert np.allclose(elab1, elab2)
print("Eref\n"+str(Eref1))
print("elab\n"+str(elab1))
In [6]:
# Apply a large random rotation to the grain
from scipy.spatial.transform import Rotation
rot = Rotation.from_euler("XYZ",(123,456,789)).as_matrix()
ubi_rotated = np.dot( ubi , rot )
np.set_printoptions(suppress=True)
Eref3, elab3 = ret = strains_from_deformation(ubi_rotated, ubi_0)
F = deformation_gradient_tensor( ubi_rotated, ubi_0 )
V,R,S = polar_decomposition(F)
Eref4, elab4 = strains_from_metric( ubi_rotated, ubi_0 )
print("Ref,grain:")
print(Eref3)
print(Eref4)
print("LabXYZ:")
print(elab3)
print(elab4)
# Ensure we get a compatible rotation back
assert np.allclose( np.dot(R, rot), np.eye(3) )
# Ensure the Eulerian Strain is the same as before rotation
assert np.allclose( Eref3, Eref4 )
assert np.allclose( elab3, elab4 )
In [7]:
# Apply a large random rotation to the reference too grain
from scipy.spatial.transform import Rotation
rot1 = Rotation.from_euler("XYZ",(123,456,789)).as_matrix()
ubi_rotated = np.dot( ubi , rot1 )
rot2 = Rotation.from_euler("XYZ",(789,5,123)).as_matrix()
ubi_ref = np.dot( ubi_0 , rot2 )
np.set_printoptions(suppress=True)
Eref6, elab6 = ret = strains_from_deformation(ubi_rotated, ubi_ref)
F = deformation_gradient_tensor( ubi_rotated, ubi_ref )
V,R,S = polar_decomposition(F)
Eref7, elab7 = strains_from_metric( ubi_rotated, ubi_ref )
print("Ref,grain:")
print(Eref6)
print(Eref7)
print("LabXYZ:")
print(elab6)
print(elab7)
# Ensure we get a compatible rotation back
assert np.allclose( np.dot(R, np.dot(rot2.T, rot1)), np.eye(3) )
# Ensure the Eulerian Strain is the same as before rotation
assert np.allclose( Eref6, Eref7 )
assert np.allclose( elab6, elab7 )
In [ ]: