class ElasticTensor:
"""
Class to represent an elastic tensor, along with methods to access it
"""
def __init__(self, mat):
mat = np.asarray(mat, dtype=np.float64)
if mat.shape != (6, 6):
# Is it upper triangular?
if list(map(len, mat)) == [6, 5, 4, 3, 2, 1]:
mat = [[0] * i + mat[i] for i in range(6)]
mat = np.array(mat)
# Is it lower triangular?
if list(map(len, mat)) == [1, 2, 3, 4, 5, 6]:
mat = [mat[i] + [0] * (5 - i) for i in range(6)]
mat = np.array(mat)
if mat.shape != (6, 6):
raise ValueError("should be a square matrix")
# Check that is is symmetric, or make it symmetric
if np.linalg.norm(np.tril(mat, -1)) == 0:
mat = mat + np.triu(mat, 1).transpose()
if np.linalg.norm(np.triu(mat, 1)) == 0:
mat = mat + np.tril(mat, -1).transpose()
if np.linalg.norm(mat - mat.transpose()) > 1e-3:
raise ValueError("should be symmetric, or triangular")
elif np.linalg.norm(mat - mat.transpose()) > 0:
mat = 0.5 * (mat + mat.transpose())
# Store it
self.c_voigt = np.array(mat)
# Put it in a more useful representation
try:
self.s_voigt = np.linalg.inv(self.c_voigt)
except np.linalg.LinalgError as e:
raise ValueError(f"Error inverting s_voigt: {e}") from e
vm = np.array(((0, 5, 4), (5, 1, 3), (4, 3, 2)))
def sv_coeff(p, q):
return 1.0 / ((1 + p // 3) * (1 + q // 3))
smat = [
[
[
[
sv_coeff(vm[i, j], vm[k, l]) * self.s_voigt[vm[i, j], vm[k, l]]
for i in range(3)
]
for j in range(3)
]
for k in range(3)
]
for l in range(3)
]
self.elasticity_tensor = np.array(smat)
@classmethod
def from_string(cls, s):
"""Initialize the elastic tensor from a string"""
if not s:
raise ValueError("no matrix was provided")
if isinstance(s, str):
# Remove braces and pipes
s = s.replace("|", " ").replace("(", " ").replace(")", " ")
# Remove empty lines
lines = [line for line in s.split("\n") if line.strip()]
if len(lines) != 6:
raise ValueError("should have six rows")
# Convert to float
try:
mat = [[float(x) for x in line.split()] for line in lines]
except ValueError as e:
raise ValueError(f"not all entries are numbers: {e}") from e
return cls(mat)
def youngs_modulus_angular(self, theta, phi):
a = angles_to_cartesian(theta, phi)
return self.youngs_modulus(a)
def youngs_modulus(self, a):
return 1 / np.einsum("ai,aj,ak,al,ijkl->a", a, a, a, a, self.elasticity_tensor)
def linear_compressibility_angular(self, theta, phi):
a = angles_to_cartesian(theta, phi)
return self.linear_compressibility(a)
def linear_compressibility(self, a):
return 1000 * np.einsum("ai,aj,ijkk->a", a, a, self.elasticity_tensor)
def shear_modulus(self, a, b):
return 0.25 / np.einsum(
"ai,aj,ak,al,ijkl->a", a, b, a, b, self.elasticity_tensor
)
def shear_modulus_angular(self, theta, phi, chi):
a = angles_to_cartesian(theta, phi)
b = angles_to_cartesian_2(theta, phi, chi)
return self.shear_modulus(a, b)
def poisson_ratio(self, a, b):
r1 = np.einsum("ai,aj,ak,al,ijkl->a", a, a, b, b, self.elasticity_tensor)
r2 = np.einsum("ai,aj,ak,al,ijkl->a", a, a, a, a, self.elasticity_tensor)
return -r1 / r2
def poisson_ratio_angular(self, theta, phi, chi):
a = angles_to_cartesian(theta, phi)
b = angles_to_cartesian_2(theta, phi, chi)
return self.poisson_ratio(a, b)
def averages(self):
A = (self.c_voigt[0, 0] + self.c_voigt[1, 1] + self.c_voigt[2, 2]) / 3
B = (self.c_voigt[1, 2] + self.c_voigt[0, 2] + self.c_voigt[0, 1]) / 3
C = (self.c_voigt[3, 3] + self.c_voigt[4, 4] + self.c_voigt[5, 5]) / 3
a = (self.s_voigt[0, 0] + self.s_voigt[1, 1] + self.s_voigt[2, 2]) / 3
b = (self.s_voigt[1, 2] + self.s_voigt[0, 2] + self.s_voigt[0, 1]) / 3
c = (self.s_voigt[3, 3] + self.s_voigt[4, 4] + self.s_voigt[5, 5]) / 3
KV = (A + 2 * B) / 3
GV = (A - B + 3 * C) / 5
KR = 1 / (3 * a + 6 * b)
GR = 5 / (4 * a - 4 * b + 3 * c)
KH = (KV + KR) / 2
GH = (GV + GR) / 2
return {
"bulk_modulus_avg": {
"voigt": KV,
"reuss": KR,
"hill": KH,
},
"shear_modulus_avg": {"voigt": GV, "reuss": GR, "hill": GH},
"youngs_modulus_avg": {
"voigt": 1 / (1 / (3 * GV) + 1 / (9 * KV)),
"reuss": 1 / (1 / (3 * GR) + 1 / (9 * KR)),
"hill": 1 / (1 / (3 * GH) + 1 / (9 * KH)),
"spackman": self.spackman_average(kind="youngs_modulus"),
},
"poissons_ratio_avg": {
"voigt": (1 - 3 * GV / (3 * KV + GV)) / 2,
"reuss": (1 - 3 * GR / (3 * KR + GR)) / 2,
"hill": (1 - 3 * GH / (3 * KH + GH)) / 2,
},
}
def plot2d(self, kind="youngs_modulus", axis="xy", npoints=100, **kwargs):
import matplotlib.pyplot as plt
import seaborn as sns
u = np.linspace(0, np.pi * 2, npoints)
v = np.zeros_like(u)
f = getattr(self, kind + "_angular")
fig, ax = plt.subplots()
fig.set_size_inches(kwargs.pop("figsize", (3, 3)))
font = "Arial"
xlims = kwargs.pop("xlim", None)
ylims = kwargs.pop("ylim", None)
# ax.set_title(f"${axis}$-plane", fontname=font, fontsize=12)
ax.set_xlabel(f"{axis[0]}", fontsize=12)
ax.set_ylabel(f"{axis[1]}", fontsize=12, rotation=0)
if axis == "xy":
v[:] = np.pi / 2
r = f(v, u)
x = r * np.cos(u)
y = r * np.sin(u)
elif axis == "xz":
r = f(u, v)
y = r * np.cos(u)
x = r * np.sin(u)
elif axis == "yz":
v[:] = np.pi / 2
r = f(u, v)
y = r * np.cos(u)
x = r * np.sin(u)
if xlims:
ax.set_xlim(*xlims)
if ylims:
ax.set_ylim(*ylims)
ax.xaxis.set_major_locator(plt.MaxNLocator(9))
ax.yaxis.set_major_locator(plt.MaxNLocator(9))
ax.plot(x, y, c="k", **kwargs)
sns.despine(ax=ax, offset=0)
ax.spines["bottom"].set_position("zero")
ax.spines["left"].set_position("zero")
if kwargs.get("grid", False):
ax.grid("minor", color="#BBBBBB", linewidth=0.5)
zero_tick = (len(ax.get_xticks()) - 1) // 2
for tick in ax.get_xticklabels():
tick.set_fontname(font)
tick.set_fontsize(4)
for tick in ax.get_yticklabels():
tick.set_fontname(font)
tick.set_fontsize(4)
ax.get_xticklabels()[zero_tick].set_visible(False)
ax.get_yticklabels()[zero_tick].set_visible(False)
ax.xaxis.set_label_coords(1.05, 0.53)
ax.yaxis.set_label_coords(0.49, 1.02)
return ax
def mesh(self, kind="youngs_modulus", subdivisions=3):
import trimesh
f = getattr(self, kind)
sphere = trimesh.creation.icosphere(subdivisions=subdivisions)
r = f(sphere.vertices)
sphere.vertices *= r[:, np.newaxis]
return sphere
def shape_descriptors(self, kind="youngs_modulus", l_max=5, **kwargs):
from chmpy.shape.shape_descriptors import make_invariants
from chmpy.shape.sht import SHT
sht = SHT(l_max=l_max)
f = getattr(self, kind)
points = sht.grid_cartesian
vals = f(points)
if kwargs.get("normalize", False):
vals = vals / self.spackman_average(kind=kind)
coeffs = sht.analyse(vals)
invariants = make_invariants(l_max, coeffs)
if kwargs.get("coefficients", False):
return coeffs, invariants
return invariants
def spackman_average(self, kind="youngs_modulus"):
mesh = self.mesh(kind=kind)
return np.mean(np.linalg.norm(mesh.vertices, axis=1), axis=0)
def voigt_rotation_matrix_6x6(self, rotation_matrix):
"""
Construct the 6x6 rotation matrix for transforming Voigt notation tensors.
Follows the standard approach for 4th-order tensor rotation in Voigt notation.
Args:
rotation_matrix (np.ndarray): 3x3 rotation matrix R
Returns:
np.ndarray: 6x6 rotation matrix for Voigt notation
"""
R = rotation_matrix
# Voigt index mapping
voigt_to_tensor = [(0, 0), (1, 1), (2, 2), (1, 2), (0, 2), (0, 1)]
T = np.zeros((6, 6))
for i in range(6):
for j in range(6):
p, q = voigt_to_tensor[i] # row indices
r, s = voigt_to_tensor[j] # column indices
# Apply the 4th-order tensor transformation rule:
# T_ij = sum over all valid combinations of R_pr * R_qs * R_rt * R_su
# where the indices match the tensor transformation pattern
# For symmetric tensors in Voigt notation, we use:
if i < 3: # diagonal terms (11, 22, 33)
if j < 3: # diagonal terms
T[i, j] = R[p, r] * R[q, s]
else: # off-diagonal terms (need factor of 2)
T[i, j] = 2 * R[p, r] * R[q, s]
else: # off-diagonal terms (23, 13, 12)
if j < 3: # diagonal terms
T[i, j] = R[p, r] * R[q, s]
else: # off-diagonal terms
T[i, j] = R[p, r] * R[q, s] + R[p, s] * R[q, r]
return T
def rotate_voigt_tensor(self, rotation_matrix):
"""
Rotate this elastic tensor using a 3x3 rotation matrix.
Uses a direct 6x6 transformation matrix construction.
Args:
rotation_matrix (np.ndarray): 3x3 rotation matrix R
Returns:
np.ndarray: Rotated 6x6 elastic tensor in Voigt notation
"""
T = self.voigt_rotation_matrix_6x6(rotation_matrix)
rotated_tensor = T @ self.c_voigt @ T.T
return rotated_tensor
def reoriented_into_standard_frame(self, vectors):
"""
Re-orient this elastic tensor from the provided cartesian frame
into a natural standard crystallographic frame for the crystal axes.
The standard frame uses the conventional crystallographic orientation:
- a along x-axis
- b in xy-plane
- c positioned to satisfy the lattice angles
Args:
vectors (np.ndarray): 3x3 matrix of lattice vectors (row major)
from the original coordinate system
Returns:
ElasticTensor: New elastic tensor in the standard frame
"""
from chmpy.crystal.unit_cell import UnitCell
from chmpy.util.num import kabsch_rotation_matrix
# Create unit cell from the provided vectors
original_uc = UnitCell(vectors)
# Create standard orientation unit cell with same lattice parameters
# This uses the standard crystallographic convention in chmpy
standard_uc = UnitCell(np.eye(3)) # temporary
standard_uc.set_lengths_and_angles(
[original_uc.a, original_uc.b, original_uc.c],
[original_uc.alpha, original_uc.beta, original_uc.gamma],
)
rotation_matrix = kabsch_rotation_matrix(original_uc.direct, standard_uc.direct)
rotated_c_voigt = self.rotate_voigt_tensor(rotation_matrix)
return ElasticTensor(rotated_c_voigt)
def __repr__(self):
s = np.array2string(
self.c_voigt, precision=4, suppress_small=True, separator=" "
)
s = s.replace("[", " ")
s = s.replace("]", " ")
return f"<ElasticTensor:\n{s}>"