-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathlattice_strain.py
196 lines (155 loc) · 8.33 KB
/
lattice_strain.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
import numpy as np
import utils
from utils import import_diff_params
from utils import tensorcomps_fromaxis
from defdap.quat import Quat
from defdap.crystal import CrystalStructure, Phase
def VE_phase_mask(phase_name, ve_response):
"""
phase_idx:
"""
# get indicies for desired phase...
phase_idx = ve_response['field_data']['phase']['meta']['phase_names'].index(phase_name)
phase_mask = ve_response['field_data']['phase']['data'] == phase_idx
return phase_idx, phase_mask
def lattice_strain(phase, phase_name, oris, ori_incs, phase_idx, phase_mask, Ee, Ee_incs, axis, tol=5):
"""
Calculate mean for component of elastic strain tensor along the defined axis direction
for material points which satisfy Bragg's condition in defined axis direction.
Parameters
----------
workflow: workflow.hdf5 file imported as python dictionary
phases: dictionary containing crystallographic parameters for defined phase_labels
axis: axis direction as string eg. "X","Y","Z" to test Bragg's condition and choose tensor component.
"""
tensor_comp, unit_vector = utils.tensorcomps_fromaxis(axis)
# ensure increments are correct
assert Ee_incs == ori_incs
incs = ori_incs
assert oris['type'] == 'quat'
assert oris['quat_component_ordering'] == 'scalar-vector'
assert oris['unit_cell_alignment']['x'] == 'a' # SAME AS DAMASK
# get orientation data for phase...
quat_comps = oris['quaternions'][:, phase_mask, :]
# convert to P = +1 convention
if oris['P'] == -1:
quat_comps[:, :, 1:] *= -1
ori = np.empty(quat_comps.shape[:-1], dtype=object)
for i, row in enumerate(quat_comps):
for j, val in enumerate(row):
ori[i, j] = Quat(*val)
latticestrain_phase = {}
plane_intensity_phase = {}
for plane_label, crystal_plane in phase.diffraction_planes.items():
print(f"Processing plane {plane_label}...")
latticestrain_plane = []
plane_intensity_plane = []
for inc in incs:
try:
inc_idx = incs.index(inc)
except ValueError:
print(f"Increment {inc} does not exist.")
continue
# determine lit up material points
lit_up = find_illuminated_points(ori[inc_idx], crystal_plane, phase, incs, unit_vector, tol)
# filter data by illuminated points
unit_vector /= np.sqrt(np.dot(unit_vector, unit_vector))
# assert np.allclose(measure_dir, measure_dir), f'expect measure dir to be [1, 0, 0]'
latticestrain_inc = Ee[inc_idx, lit_up, tensor_comp-1, tensor_comp-1]
# append lattice strain for plane
latticestrain_plane.append(latticestrain_inc)
# account for multiplicity of planes here?
plane_intensity_plane.append(np.count_nonzero(lit_up)/6)
# append lattice strain for plane to phase
latticestrain_phase[plane_label] = latticestrain_plane
plane_intensity_phase[plane_label] = plane_intensity_plane
return latticestrain_phase, plane_intensity_phase
def find_illuminated_points(oris, crystal_plane, phase, incs, measure_dir, tol=5, batch_size=100000):
"""
Find the crystal orientations that will contribute to the diffraction peak of
a given crystal plane for a given load direction.
Parameters
----------
oris : numpy.ndarray(defdap.quat.Quat)
Orientation of each point
crystal_plane : numpy.ndarray
Crystal plane expressed in Miller-Bravais indicies
phase : defdap.crystal.Phase
Defdap phase
measure_dir : numpy.ndarray
Direction strain is measured in. Currently only working for z direction (0,0,1)
tolerance : float
Tolerance for meeting the Bragg condition, in degrees
batch_size : int, optional
Number of orientations to process at one time
"""
sym_group = phase.crystalStructure.name
if sym_group == 'hexagonal':
crystal_plane_ortho = transform_plane(crystal_plane, phase.cOverA) # only for HCP
else:
crystal_plane_ortho = crystal_plane.astype(float)
crystal_plane_ortho /= np.sqrt(np.dot(crystal_plane_ortho, crystal_plane_ortho))
measure_dir /= np.sqrt(np.dot(measure_dir, measure_dir))
tol = np.cos(tol * np.pi / 180.)
num_points = oris.shape[0]
num_batchs = num_points // batch_size + 1
lit_up = np.zeros(num_points, dtype=bool)
for i in range(num_batchs):
oris_batch = oris[i*batch_size:(i+1)*batch_size]
# transform load_dir to crystal coords for all orientations and symmetries
quat_comps_sym = Quat.calcSymEqvs(oris_batch, sym_group)
# temp variables to use below
quat_dot_vec = (quat_comps_sym[:, 1, :] * measure_dir[0] +
quat_comps_sym[:, 2, :] * measure_dir[1] +
quat_comps_sym[:, 3, :] * measure_dir[2])
temp = (np.square(quat_comps_sym[:, 0, :]) - np.square(quat_comps_sym[:, 1, :]) -
np.square(quat_comps_sym[:, 2, :]) - np.square(quat_comps_sym[:, 3, :]))
# array to store crytal directions for all orientations and symmetries
measure_dir_crystal = np.empty((3, quat_comps_sym.shape[0], quat_comps_sym.shape[2]))
# (quat_comps_sym * measure_dir) * quat_comps_sym.conjugate
measure_dir_crystal[0] = (2 * quat_dot_vec * quat_comps_sym[:, 1, :] +
temp * measure_dir[0] +
2 * quat_comps_sym[:, 0, :] * (quat_comps_sym[:, 2, :] * measure_dir[2] -
quat_comps_sym[:, 3, :] * measure_dir[1]))
measure_dir_crystal[1] = (2 * quat_dot_vec * quat_comps_sym[:, 2, :] +
temp * measure_dir[1] +
2 * quat_comps_sym[:, 0, :] * (quat_comps_sym[:, 3, :] * measure_dir[0] -
quat_comps_sym[:, 1, :] * measure_dir[2]))
measure_dir_crystal[2] = (2 * quat_dot_vec * quat_comps_sym[:, 3, :] +
temp * measure_dir[2] +
2 * quat_comps_sym[:, 0, :] * (quat_comps_sym[:, 1, :] * measure_dir[1] -
quat_comps_sym[:, 2, :] * measure_dir[0]))
# normalise vectors (just in case)
measure_dir_crystal /= np.sqrt(np.einsum('ijk,ijk->jk', measure_dir_crystal, measure_dir_crystal))
dot_prod = np.einsum('ijk,i->jk', measure_dir_crystal, crystal_plane_ortho)
lit_up[i*batch_size:(i+1)*batch_size] = np.any(np.abs(dot_prod) > tol, axis=0)
# print("\rDone. {:} points of {:} will contribute to diffraction peak.".format(np.count_nonzero(lit_up), num_points))
return lit_up
def transform_plane(crystal_plane, c_over_a):
"""Transform crystal plane expressed in Miller-Bravais indicies to a unit vector in orthonormal basis."""
# Convert plane and dir from Miller-Bravais to Miller indicies
crystal_plane_m = crystal_plane[[0, 1, 3]] # (M-B hkil -> M hkl)
# Create L matrix. Transformation from crystal to orthonormal coords
l_matrix = CrystalStructure.lMatrix(1, 1, c_over_a, np.pi / 2, np.pi / 2, np.pi * 2 / 3)
# Create Q matrix for transforming planes to orthonormal coords
q_matrix = CrystalStructure.qMatrix(l_matrix)
# Transform into orthonormal basis and then normalise
crystal_plane_ortho = np.matmul(q_matrix, crystal_plane_m)
crystal_plane_ortho /= np.sqrt(np.dot(crystal_plane_ortho, crystal_plane_ortho))
return crystal_plane_ortho
def process_strain(Ee, lit_up, measure_dir, comp):
"""
Filter the strain data by the illuminated points.
Parameters
----------
Ee : numpy.ndarray
Elastic strain tensor at each point, shape (num_points, 3, 3)
lit_up : numpy.ndarray(bool)
Boolean array indicating if each point is illuminated, shape (num_points,)
measure_dir : numpy.ndarray
Direction strain is to be measured in. Currently only working for z direction (0,0,1)
"""
measure_dir /= np.sqrt(np.dot(measure_dir, measure_dir))
# assert np.allclose(measure_dir, measure_dir), f'expect measure dir to be [1, 0, 0]'
measured_strain = Ee[lit_up, comp, comp]
return measured_strain