Handle large (Green-Lagrange) strains in `xfab.tools.b_to_epsilon`

[AxelHenningsson]

The current implementation of xfab.tools.b_to_epsilon uses the small strain tensor definition

$\boldsymbol{\epsilon} = \dfrac{1}{2}(\boldsymbol{F}^T + \boldsymbol{F}) - \boldsymbol{I}$.

I suggest that the default behavior is changed to the large Green-Lagrange strain tensor

$\boldsymbol{E} = \dfrac{1}{2}(\boldsymbol{F}^T \boldsymbol{F} - \boldsymbol{I})$.

The back-transformation (xfab.tools.epsilon_to_b) is then achieved by a Cholesky factorization (when the strain tensor is feasible).

I am thinking about something along the lines of the below. Thoughts? 🙂

from scipy.linalg import cholesky

def strain_as_tensor( strain_vector ):
    e11, e12, e13, e22, e23, e33 = strain_vector
    return np.asarray([ [e11, e12, e13],
                        [e12, e22, e23],
                        [e13, e23, e33] ], float)

def strain_as_vector( strain_tensor ):
    return list(strain_tensor[0, :]) + list(strain_tensor[1, 1:]) + [strain_tensor[2, 2]]

def b_to_epsilon(B_matrix, unit_cell):
    B = np.asarray(B_matrix, float)
    B0 = form_b_mat(unit_cell)
    F = np.dot(B0, np.linalg.inv(B))
    strain_tensor = 0.5*(F.T.dot(F) - np.eye(3))  # large deformations
    return strain_as_vector( strain_tensor )

def epsilon_to_b(epsilon, unit_cell):
    strain_tensor = strain_as_tensor( epsilon )
    C = 2*strain_tensor + np.eye(3)
    eigen_vals = np.linalg.eigvalsh( C )
    if np.any( np.array(eigen_vals) < 0 ):
        raise ValueError("Unfeasible strain tensor with value: "+str(strain_as_vector(strain_tensor))+ \
            ", will invert the unit cell with negative deformation gradient tensor determinant" )
    F =  cholesky(C)
    B0 = form_b_mat(unit_cell)
    B = np.linalg.inv(F).dot(B0)
    return B

if __name__=='__main__':
    unit_cell = [4.926, 4.926, 5.4189, 90., 90., 120.]
    eps1 = 25 * 1e-4 * (np.random.rand(6,)-0.5)
    B = epsilon_to_b(eps1, unit_cell)
    eps2 = b_to_epsilon(B, unit_cell)
    print(eps1)
    print(eps2)
    assert np.allclose( eps1, eps2 )

[jonwright]

Hi Axel

I think that "epsilon" is the documented definition in the FitAllB papers, so changing it might cause some people to get different numbers, also for polyxsim. For small strains they are all supposed to be similar anyway. Maybe a different name for epsilon, like "strain" for example? Would "b_to_strain" and "strain_to_b" be better anyway?

There is a related code in ImageD11 to give a range of different strain tensors, it could be useful to debug that if it disagrees with this.

Best

Jon

[AxelHenningsson]

So you implemented like a set-Hill family of strain tensors in ImageD11, cool, I like that🙂 .

Perhaps we could put something like the below in the new xfab.laue then, not breaking anything and keeping the old functions for backwards comparability.

Will draft a PR on this with unit tests.

def strain_from_def_grad(F, strain_measure='small'):
    """Compute a strain tensor measure form a deformation gradient tensor F
    """
    if strain_measure=='green-lagrange':
        strain_tensor = 0.5*(F.T.dot(F) - np.eye(3))
    elif strain_measure=='small':
        strain_tensor = 0.5*(F.T + F) - np.eye(3)
    else:
        raise ValueError('No such strain_measure : '+str(strain_measure))
    return strain_as_vector( strain_tensor )

def b_to_strain(B_matrix, unit_cell, strain_measure='green-lagrange'):
    """Extract the strain tensor from the deformed crystal B matrix and a reference state unit cell.
    """
    B = np.asarray(B_matrix, float)
    B0 = form_b_mat(unit_cell)
    F = np.dot(B0, np.linalg.inv(B))
    return strain_from_def_grad(F, strain_measure)

def strain_to_b(strain, unit_cell, strain_measure='green-lagrange'):
    """Extract the deformed crystal B matrix from a strain tensor and a reference state unit cell.
    """
    strain_tensor = strain_as_tensor( strain )
    B0 = form_b_mat(unit_cell)

    if strain_measure=='green-lagrange':
        C = 2*strain_tensor + np.eye(3)
        eigen_vals = np.linalg.eigvalsh( C )
        if np.any( np.array(eigen_vals) < 0 ):
            raise ValueError("Unfeasible strain tensor with value: "+str(strain_as_vector(strain_tensor))+ \
                ", will invert the unit cell with negative deformation gradient tensor determinant" )
        F = cholesky(C)
    elif strain_measure=='small':
        L = 2 * (strain_tensor + np.eye(3))
        F = np.array([  [ L[0,0]/2., L[0,1],    L[0,2]    ],
                        [    0,      L[1,1]/2., L[1,2]    ],
                        [    0,        0,       L[2,2]/2. ]  ])
        if np.linalg.det( F ) <=0:
            raise ValueError("Unfeasible strain tensor with value: "+str(strain_as_vector(strain_tensor))+ \
                ", will invert the unit cell with negative deformation gradient tensor determinant" )
    else:
        raise ValueError('No such strain_measure : '+str(strain_measure))

    return np.linalg.inv(F).dot(B0)