# %%
import warnings
from math import sqrt
import numpy as np
from typing import Union
try:
import torch
Tensor = torch.Tensor
except ImportError:
Tensor = np.ndarray
##########################
# Crystal symmetries
##########################
def enforce_trigonal(S0 : np.ndarray) -> np.ndarray:
S = S0.copy()
S[1,1]=S[0,0]
S[1,2]=S[0,2]=S[2,0]=S[2,1]
S[3,3]=S[4,4]
S[5,5]=2*(S[0,0]-S[0,1])
S[1,3]=-S[0,3]
S[3,1]=S[1,3]
S[3,0]=S[0,3]
S[4,5]=2*S[0,3]
S[5,4]=S[4,5]
S[1,4]=-S[0,4]
S[4,1]=S[1,4]
S[4,0]=S[0,4]
S[3,5]=2*S[1,4]
S[5,3]=S[3,5]
return S
def enforce_tetragonal(S0 : np.ndarray) -> np.ndarray:
S = S0.copy()
S = 0.5*(S+S.T)
S[1,1]=S[0,0]
S[1,2]=S[0,2]=S[2,0]=S[2,1]
S[3,3]=S[4,4]
S[0:3,3:6]=S[3:6,0:3]=0
S[3,4:6]=S[4:6,3]=0
S[4,5]=S[5,4]=0
return S
def enforce_hexagonal(S0 : np.ndarray) -> np.ndarray:
S = S0.copy()
S = 0.5*(S+S.T)
S[1,1]=S[0,0]
S[1,2]=S[0,2]=S[2,0]=S[2,1]
S[3,3]=S[4,4]
S[5,5]=2*(S[0,0]-S[0,1])
S[0:3,3:6]=S[3:6,0:3]=0
S[3,4:6]=S[4:6,3]=0
S[4,5]=S[5,4]=0
return S
def enforce_orthotropic(S0 : np.ndarray) -> np.ndarray:
S = S0.copy()
S = 0.5*(S+S.T)
S[0:3,3:6]=S[3:6,0:3]=0
S[3,4:6]=S[4:6,3]=0
S[4,5]=S[5,4]=0
return S
# Compliance tensors in Voigt notation for different crystal symmetries
def orthotropic_S(Ex, Ey, Ez, Gyz, Gxz, Gxy, nuxy, nuyz, nuxz) -> np.ndarray:
nuzy = nuyz * Ez/Ey
nuzx = nuxz * Ez/Ex
nuyx = nuxy * Ey/Ex
S=np.array([[1/Ex,-nuxy/Ex,-nuxz/Ex,0,0,0],
[-nuyx/Ey,1/Ey,-nuyz/Ey,0,0,0],
[-nuzx/Ez,-nuzy/Ez,1/Ez,0,0,0],
[0,0,0,1/Gyz,0,0],
[0,0,0,0,1/Gxz,0],
[0,0,0,0,0,1/Gxy]])
return S
def monoclinic_S(Ex, Ey, Ez, nuxy, nuxz, nuyz, nuxxyz, nuyyyz, nuzzyz, nuxzxy, Gyz, Gxz, Gxy) -> np.ndarray:
s1 = 1/Ex; s2 = 1/Ey; s3 = 1/Ez
s4 = 1/Gyz; s5=1/Gxz; s6=1/Gxy
s7 = nuxy/Ex; s8 = nuxz/Ex; s9 = nuyz/Ey
s10 = nuxxyz/Gyz; s11=nuyyyz/Gxz; s12=nuzzyz/Gxy; s13=nuxzxy/Gxz
S=np.array([[s1,-s7,-s8,s10,0,0],
[-s7,s2,-s9,s11,0,0],
[-s8,-s9,s3,s12,0,0],
[s10,s11,s12,s4,0,0],
[0,0,0,0,s5,s13],
[0,0,0,0,s13,s6]])
return S
def tetragonal_S(Ex, Ez, nuxy, nuxz, Gyz, Gxy) -> np.ndarray:
s1 = 1/Ex
s2 = 1/Ez
s3 = 1/Gyz
s4 = 1/Gxy
s5 = nuxy/Ex
s6 = nuxz/Ex
S=np.array([[s1,-s5,-s6,0,0,0],
[-s5,s1,-s6,0,0,0],
[-s6,-s6,s2,0,0,0],
[0,0,0,s3,0,0],
[0,0,0,0,s3,0],
[0,0,0,0,0,s4]])
return S
def trigonal_S(Ex, Ez, nuxy, nuxz, Gyz, nuxxyz) -> np.ndarray:
s1 = 1/Ex
s2 = 1/Ez
s3 = 1/Gyz
s5 = nuxy/Ex
s6 = nuxz/Ex
s4 = nuxxyz/Gyz
S=np.array([[s1,-s5,-s6,s4,0,0],
[-s5,s1,-s6,-s4,0,0],
[-s6,-s6,s2,0,0,0],
[s4,-s4,0,s3,0,0],
[0,0,0,0,s3,2*s4],
[0,0,0,0,2*s4,2*(s1+s5)]])
return S
def cubic_S(E, nu, G) -> np.ndarray:
s1 = 1/E
s2 = nu/E
s3 = 1/G
S=np.array([[s1,-s2,-s2,0,0,0],
[-s2,s1,-s2,0,0,0],
[-s2,-s2,s1,0,0,0],
[0,0,0,s3,0,0],
[0,0,0,0,s3,0],
[0,0,0,0,0,s3]])
return S
def hexagonal_S(Ex, Ez, nuxy, nuxz, Gyz) -> np.ndarray:
s1 = 1/Ex
s2 = nuxy/Ex
s3 = nuxz/Ex
s4 = 1/Ez
s5 = 1/Gyz
S=np.array([[s1,-s2,-s3,0,0,0],
[-s2,s1,-s3,0,0,0],
[-s3,-s3,s4,0,0,0],
[0,0,0,s5,0,0],
[0,0,0,0,s5,0],
[0,0,0,0,0,2*(s1+s2)]])
return S
def isotropic_S(E=1, nu=0.3) -> np.ndarray:
G = 0.5 * E / (1 + nu)
S=np.array([[1/E,-nu/E,-nu/E,0,0,0],
[-nu/E,1/E,-nu/E,0,0,0],
[-nu/E,-nu/E,1/E,0,0,0],
[0,0,0,1/G,0,0],
[0,0,0,0,1/G,0],
[0,0,0,0,0,1/G]])
return S
############################
# Rotation functions
############################
[docs]
def rotate_4th_order(T : np.ndarray, Q : np.ndarray) -> np.ndarray:
"""Rotate 4th order tensor to a new coordinate system
Args:
T (np.ndarray): 3x3 rotation matrix
Q (np.ndarray): 4th order tensor
Returns:
np.ndarray: Rotated 4th order tensor
"""
out = np.einsum('ia, jb, kc, ld, abcd->ijkl', Q,Q,Q,Q,T)
return out
[docs]
def rotate_Voigt_stiffness(C : np.ndarray, R : np.ndarray) -> np.ndarray:
"""Rotate 6x6 stiffness matrix in Voigt notation
Args:
C (np.ndarray): 6x6 Voigt stiffness matrix
R (np.ndarray): 3x3 rotation matrix
Returns:
np.ndarray: Rotated 6x6 Voigt stiffness matrix
"""
K1 = R**2
K2 = R[:,[1,2,0]]*R[:,[2,0,1]]
K3 = R[[1,2,0],:]*R[[2,0,1],:]
K4 = R[[1,2,0],:][:,[1,2,0]]*R[[2,0,1],:][:,[2,0,1]] + R[[1,2,0],:][:,[2,0,1]]*R[[2,0,1],:][:,[1,2,0]]
K = np.block([
[K1, 2*K2],
[K3, K4]
])
return K @ C @ (K.T)
[docs]
def rotate_Voigt_compliance(S : np.ndarray, R : np.ndarray) -> np.ndarray:
"""Rotate 6x6 compliance matrix in Voigt notation
Args:
S (np.ndarray): 6x6 Voigt compliance matrix
R (np.ndarray): 3x3 rotation matrix
Returns:
np.ndarray: Rotated 6x6 Voigt compliance matrix
"""
K1 = R**2
K2 = R[:,[1,2,0]]*R[:,[2,0,1]]
K3 = R[[1,2,0],:]*R[[2,0,1],:]
K4 = R[[1,2,0],:][:,[1,2,0]]*R[[2,0,1],:][:,[2,0,1]] + R[[1,2,0],:][:,[2,0,1]]*R[[2,0,1],:][:,[1,2,0]]
K = np.block([
[K1, 2*K2],
[K3, K4]
])
K_T = np.block([
[K1, K2],
[2*K3, K4]
])
return K_T @ S @ np.linalg.inv(K)
[docs]
def Voigt_rot_matrix_stress_inplace(Q: Union[np.ndarray, Tensor], R: Union[np.ndarray, Tensor]) -> None:
"""Fill the 6x6 rotation matrix for Voigt stress vector
Args:
Q (Union[np.ndarray, Tensor]): 3x3 rotation matrix
R (Union[np.ndarray, Tensor]): 6x6 rotation matrix to be filled
"""
assert Q.shape[-1] == 3 and Q.shape[-2] == 3
assert R.shape[-1] == 6 and R.shape[-2] == 6
A11 = Q**2
inds0 = [1,0,0]
inds1 = [2,2,1]
A12 = 2 * Q[...,inds0] * Q[...,inds1]
A21 = Q[...,inds0,:] * Q[...,inds1,:]
B1 = Q[...,inds0,:]
B2 = Q[...,inds1,:]
A22 = B1[...,inds1] * B2[...,inds0] + B1[...,inds0] * B2[...,inds1]
R[...,0:3,0:3] = A11
R[...,0:3,3:6] = A12
R[...,3:6,0:3] = A21
R[...,3:6,3:6] = A22
[docs]
def Voigt_rot_matrix_stress_numpy(Q: np.ndarray) -> np.ndarray:
"""Generate 6x6 rotation matrix for Voigt stress vector
Args:
Q (np.ndarray): 3x3 rotation matrix
Returns:
np.ndarray: 6x6 rotation matrix for Voigt stress vector
"""
if Q.ndim == 2:
R = np.zeros((6,6))
elif Q.ndim == 3:
R = np.zeros((Q.shape[0], 6,6))
Voigt_rot_matrix_stress_inplace(Q, R)
return R
[docs]
def Voigt_rot_matrix_strain_inplace(Q: Union[np.ndarray, Tensor], R: Union[np.ndarray, Tensor]) -> None:
"""Fill the 6x6 rotation matrix for Voigt strain vector
Args:
Q (Union[np.ndarray, Tensor]): 3x3 rotation matrix
R (Union[np.ndarray, Tensor]): 6x6 rotation matrix to be filled
"""
assert Q.shape[-1] == 3 and Q.shape[-2] == 3
assert R.shape[-1] == 6 and R.shape[-2] == 6
A11 = Q**2
inds0 = [1,0,0]
inds1 = [2,2,1]
A12 = Q[...,inds0] * Q[...,inds1]
A21 = 2*Q[...,inds0,:] * Q[...,inds1,:]
B1 = Q[...,inds0,:]
B2 = Q[...,inds1,:]
A22 = B1[...,inds1] * B2[...,inds0] + B1[...,inds0] * B2[...,inds1]
R[...,0:3,0:3] = A11
R[...,0:3,3:6] = A12
R[...,3:6,0:3] = A21
R[...,3:6,3:6] = A22
[docs]
def Voigt_rot_matrix_strain_numpy(Q: np.ndarray) -> np.ndarray:
"""Generate 6x6 rotation matrix for Voigt strain vector
Args:
Q (np.ndarray): 3x3 rotation matrix
Returns:
np.ndarray: 6x6 rotation matrix for Voigt strain vector
"""
if Q.ndim == 2:
R = np.zeros((6,6))
elif Q.ndim == 3:
R = np.zeros((Q.shape[0], 6,6))
Voigt_rot_matrix_strain_inplace(Q, R)
return R
[docs]
def Mandel_rot_matrix_inplace(Q: Union[np.ndarray, Tensor], R: Union[np.ndarray, Tensor]) -> None:
"""Fill the 6x6 rotation matrix for Mandel stress/strain vector
Args:
Q (Union[np.ndarray, Tensor]): 3x3 rotation matrix
R (Union[np.ndarray, Tensor]): 6x6 rotation matrix to be filled
"""
assert Q.shape[-1] == 3 and Q.shape[-2] == 3
assert R.shape[-1] == 6 and R.shape[-2] == 6
A11 = Q**2
A12 = sqrt(2) * Q[...,[1,0,0]] * Q[...,[2,2,1]]
B1 = Q[...,[1,0,0],:]
B2 = Q[...,[2,2,1],:]
A21 = sqrt(2) * B1 * B2
A22 = B1[...,[1,0,0]] * B2[...,[2,2,1]] + B1[...,[2,2,1]] * B2[...,[1,0,0]]
R[...,0:3,0:3] = A11
R[...,0:3,3:6] = A12
R[...,3:6,0:3] = A21
R[...,3:6,3:6] = A22
[docs]
def Mandel_rot_matrix_numpy(Q: np.ndarray) -> np.ndarray:
"""Generate 6x6 rotation matrix for Mandel stress/strain vector
Args:
Q (np.ndarray): 3x3 rotation matrix
Returns:
np.ndarray: 6x6 rotation matrix for Mandel stress/strain vector
"""
if Q.ndim == 2:
R = np.zeros((6,6))
elif Q.ndim == 3:
R = np.zeros((Q.shape[0], 6,6))
Mandel_rot_matrix_inplace(Q, R)
return R
# %%
[docs]
def tens_2d_to_Mandel_inplace(s: Union[np.ndarray, Tensor], x: Union[np.ndarray, Tensor]) -> None:
"""Fill Mandel vector from 2nd order tensor
Args:
s (Union[np.ndarray, Tensor]): 2nd order tensor
x (Union[np.ndarray, Tensor]): Mandel vector
"""
assert s.shape[-1] == 3 and s.shape[-2] == 3
assert x.shape[-1] == 6
x[...,0] = s[...,0,0]
x[...,1] = s[...,1,1]
x[...,2] = s[...,2,2]
s2 = sqrt(2)
x[...,3] = s2 * s[...,1,2]
x[...,4] = s2 * s[...,0,2]
x[...,5] = s2 * s[...,0,1]
[docs]
def tens_2d_to_Mandel_numpy(s: np.ndarray) -> np.ndarray:
"""Convert 2nd order tensor to Mandel vector
Args:
s (np.ndarray): 2nd order tensor
Returns:
np.ndarray: Mandel vector
"""
assert np.allclose(s, s.swapaxes(-1,-2)), "s is not symmetric"
if s.ndim == 2:
x = np.zeros((6))
elif s.ndim == 3:
x = np.zeros((s.shape[0], 6))
tens_2d_to_Mandel_inplace(s, x)
return x
[docs]
def stress_to_Voigt_inplace(s: Union[np.ndarray, Tensor], x: Union[np.ndarray, Tensor]) -> None:
"""Fill Voigt notation vector from 2nd order stress tensor
Args:
s (Union[np.ndarray, Tensor]): 2nd order stress tensor
x (Union[np.ndarray, Tensor]): Voigt stress vector
"""
assert s.shape[-1] == 3 and s.shape[-2] == 3
assert x.shape[-1] == 6
x[...,0] = s[...,0,0]
x[...,1] = s[...,1,1]
x[...,2] = s[...,2,2]
x[...,3] = s[...,1,2]
x[...,4] = s[...,0,2]
x[...,5] = s[...,0,1]
[docs]
def stress_to_Voigt_numpy(s: np.ndarray) -> np.ndarray:
"""Convert 2nd order stress tensor to Voigt notation
Args:
s (np.ndarray): 2nd order stress tensor
Returns:
np.ndarray: Voigt stress vector
"""
assert np.allclose(s, s.swapaxes(-1,-2)), "s is not symmetric"
if s.ndim == 2:
x = np.zeros((6))
elif s.ndim == 3:
x = np.zeros((s.shape[0], 6))
stress_to_Voigt_inplace(s, x)
return x
[docs]
def strain_to_Voigt_inplace(e: Union[np.ndarray, Tensor], x: Union[np.ndarray, Tensor]) -> None:
"""Fill Voigt notation vector from 2nd order strain tensor
Args:
e (Union[np.ndarray, Tensor]): 2nd order strain tensor
x (Union[np.ndarray, Tensor]): Voigt strain vector
"""
assert e.shape[-1] == 3 and e.shape[-2] == 3
assert x.shape[-1] == 6
x[...,0] = e[...,0,0]
x[...,1] = e[...,1,1]
x[...,2] = e[...,2,2]
x[...,3] = 2 * e[...,1,2]
x[...,4] = 2 * e[...,0,2]
x[...,5] = 2 * e[...,0,1]
[docs]
def strain_to_Voigt_numpy(e: np.ndarray) -> np.ndarray:
"""Convert 2nd order strain tensor to Voigt notation
Args:
e (np.ndarray): 2nd order strain tensor
Returns:
np.ndarray: Voigt strain vector
"""
assert np.allclose(e, e.swapaxes(-1,-2)), "e is not symmetric"
if e.ndim == 2:
x = np.zeros((6))
elif e.ndim == 3:
x = np.zeros((e.shape[0], 6))
strain_to_Voigt_inplace(e, x)
return x
[docs]
def compliance_Voigt_to_4th_order(S : np.ndarray) -> np.ndarray:
"""Convert Voigt notation compliance matrix to 4th order tensor
Args:
S (np.ndarray): 6x6 Voigt notation compliance matrix
Returns:
np.ndarray: 4th order compliance tensor
"""
if S.ndim==2:
_S = S.reshape((1,6,6))
elif S.ndim==3:
_S = S
# Convert S_ij to S_abcd
assert np.allclose(_S, _S.transpose((0,2,1))), "S is not symmetric"
S_4 = np.zeros((_S.shape[0],3,3,3,3))
for i in range(6):
if i<3:
a = b = i
elif i==3:
a = 1; b = 2
elif i==4:
a = 0; b = 2
else:
a = 0; b = 1
for j in range(i,6):
if j<3:
c = d = j
elif j==3:
c = 1; d = 2
elif j==4:
c = 0; d = 2
else:
c = 0; d = 1
Sij = _S[:,i,j]
Sij = Sij if i<3 else Sij/2
Sij = Sij if j<3 else Sij/2
S_4[:,a,b,c,d] = Sij
S_4[:,a,b,d,c] = Sij
S_4[:,b,a,c,d] = Sij
S_4[:,b,a,d,c] = Sij
S_4[:,c,d,a,b] = Sij
S_4[:,c,d,b,a] = Sij
S_4[:,d,c,a,b] = Sij
S_4[:,d,c,b,a] = Sij
if S.ndim==2:
return S_4[0]
else:
return S_4
# %%
[docs]
def compliance_4th_order_to_Voigt(S : np.ndarray) -> np.ndarray:
"""Convert 4th order compliance tensor to Voigt notation
Args:
S (np.ndarray): 4th order compliance tensor
Returns:
np.ndarray: 6x6 Voigt notation compliance matrix
"""
if S.ndim==4:
_S = S.reshape((1,3,3,3,3))
elif S.ndim==5:
_S = S
# Convert S_abcd to S_ij
assert np.allclose(_S, np.transpose(_S, (0,1,2,4,3)))
assert np.allclose(_S, np.transpose(_S, (0,2,1,3,4)))
assert np.allclose(_S, np.transpose(_S, (0,3,4,1,2)))
S_2 = np.zeros((_S.shape[0],6,6))
for i in range(6):
if i<3:
a = b = i
elif i==3:
a = 1; b = 2
elif i==4:
a = 0; b = 2
else:
a = 0; b = 1
for j in range(i,6):
if j<3:
c = d = j
elif j==3:
c = 1; d = 2
elif j==4:
c = 0; d = 2
else:
c = 0; d = 1
Sij = _S[:,a,b,c,d] if j<3 else 2*_S[:,a,b,c,d]
Sij = Sij if i<3 else 2*Sij
S_2[:,i,j] = Sij
S_2[:,j,i] = Sij
if S.ndim==4:
return S_2[0]
else:
return S_2
# %%
[docs]
def stiffness_4th_order_to_Voigt(C: np.ndarray) -> np.ndarray:
"""Convert 4th order stiffness tensor to Voigt notation
Args:
C (np.ndarray): 4th order stiffness tensor
Raises:
ValueError: Wrong shape of C
Returns:
np.ndarray: 6x6 Voigt notation stiffness matrix
"""
if C.ndim==5:
_C = C
elif C.ndim==4:
_C = C.reshape((1,3,3,3,3))
else:
raise ValueError(f'Wrong shape of C: C.shape={C.shape}')
C_2 = np.zeros((_C.shape[0], 6,6))
for i in range(6):
if i<3:
a = b = i
elif i==3:
a = 1; b = 2
elif i==4:
a = 0; b = 2
else:
a = 0; b = 1
for j in range(i,6):
if j<3:
c = d = j
elif j==3:
c = 1; d = 2
elif j==4:
c = 0; d = 2
else:
c = 0; d = 1
Cij = _C[:,a,b,c,d]
C_2[:,i,j] = Cij
C_2[:,j,i] = Cij
if C.ndim==5:
return C_2
else:
return C_2[0,:,:]
# %%
[docs]
def stiffness_Voigt_to_4th_order(C: np.ndarray) -> np.ndarray:
"""Convert Voigt notation stiffness matrix to 4th order tensor
Args:
C (np.ndarray): 6x6 Voigt notation stiffness matrix
Returns:
np.ndarray: 4th order stiffness tensor
"""
# convert C_ij to C_abcd
if C.ndim==3:
_C = C
elif C.ndim==2:
_C = C.reshape((1,6,6))
assert np.allclose(_C, _C.transpose((0,2,1)))
C4 = np.zeros((_C.shape[0], 3,3,3,3))
# 1-based continuum mechanics notation
for i in range(1,7):
if i<=3:
a = b = i
elif i==4:
a = 2; b = 3
elif i==5:
a = 1; b = 3
elif i==6:
a = 1; b = 2
for j in range(i,7):
if j<=3:
c = d = j
elif j==4:
c = 2; d = 3
elif j==5:
c = 1; d = 3
elif j==6:
c = 1; d = 2
C4[:, a-1, b-1, c-1, d-1] = _C[:, i-1, j-1]
C4[:, b-1, a-1, c-1, d-1] = _C[:, i-1, j-1]
C4[:, a-1, b-1, d-1, c-1] = _C[:, i-1, j-1]
C4[:, b-1, a-1, d-1, c-1] = _C[:, i-1, j-1]
C4[:, c-1, d-1, a-1, b-1] = _C[:, i-1, j-1]
C4[:, c-1, d-1, b-1, a-1] = _C[:, i-1, j-1]
C4[:, d-1, c-1, a-1, b-1] = _C[:, i-1, j-1]
C4[:, d-1, c-1, b-1, a-1] = _C[:, i-1, j-1]
if C.ndim==3:
return C4
else:
return C4[0,:,:,:,:]
[docs]
def stiffness_Voigt_to_Mandel(C: np.ndarray) -> np.ndarray:
"""Convert Voigt notation stiffness matrix to Mandel notation
Args:
C (np.ndarray): 6x6 Voigt notation stiffness matrix
Returns:
np.ndarray: 6x6 Mandel notation stiffness matrix
"""
s2 = np.sqrt(2)
mask = np.block([[np.ones((3,3)), s2*np.ones((3,3))],
[s2*np.ones((3,3)), 2*np.ones((3,3))]])
C_2 = C[...,:,:] * mask
return C_2
[docs]
def stiffness_Mandel_to_Voigt(C: np.ndarray) -> np.ndarray:
"""Convert Mandel notation stiffness matrix to Voigt notation
Args:
C (np.ndarray): 6x6 Mandel notation stiffness matrix
Returns:
np.ndarray: 6x6 Voigt notation stiffness matrix
"""
s2 = np.sqrt(2)
mask = np.block([[np.ones((3,3)), s2*np.ones((3,3))],
[s2*np.ones((3,3)), 2*np.ones((3,3))]])
C_2 = C[...,:,:] / mask
return C_2
[docs]
def compliance_Voigt_to_Mandel(S: np.ndarray) -> np.ndarray:
"""Convert Voigt notation compliance matrix to Mandel notation
Args:
S (np.ndarray): 6x6 Voigt notation compliance matrix
Returns:
np.ndarray: 6x6 Mandel notation compliance matrix
"""
s2 = np.sqrt(2)
mask = np.block([[np.ones((3,3)), s2*np.ones((3,3))],
[s2*np.ones((3,3)), 2*np.ones((3,3))]])
S_2 = S[...,:,:] / mask
return S_2
[docs]
def compliance_Mandel_to_Voigt(S: np.ndarray) -> np.ndarray:
"""Convert Mandel notation compliance matrix to Voigt notation
Args:
S (np.ndarray): 6x6 Mandel notation compliance matrix
Returns:
np.ndarray: 6x6 Voigt notation compliance matrix
"""
s2 = np.sqrt(2)
mask = np.block([[np.ones((3,3)), s2*np.ones((3,3))],
[s2*np.ones((3,3)), 2*np.ones((3,3))]])
S_2 = S[...,:,:] * mask
return S_2
[docs]
def numpy_cart_4_to_Mandel(_C: np.ndarray) -> np.ndarray:
"""Convert 4th order tensor to 6x6 Mandel notation
Args:
_C (np.ndarray): 4th order tensor of shape (...,3,3,3,3)
Returns:
np.ndarray: 6x6 Mandel notation tensor
"""
s2 = np.sqrt(2)
mask = np.block([[np.ones((3,3)), s2*np.ones((3,3))],
[s2*np.ones((3,3)), 2*np.ones((3,3))]])
if _C.ndim == 4:
C = _C[None,...] # add batch dimension
else:
C = _C
C2_shape = C.shape[:-2] + (6,6)
C_2 = np.zeros(C2_shape)
_cart_4_tensor_to_Mandel_inplace(C, C_2, mask)
if _C.ndim == 4:
C_2 = C_2[0, ...] # remove batch dimension
return C_2
[docs]
def numpy_Mandel_to_cart_4(C: np.ndarray) -> np.ndarray:
"""Convert 6x6 Mandel notation tensor to 4th order tensor
Args:
C (np.ndarray): Mandel notation tensor (...,6,6)
Returns:
np.ndarray: 4th order tensor of shape (...,3,3,3,3)
"""
# convert C_ij to C_abcd
if C.ndim>2:
assert C.shape[-1]==6 and C.shape[-2]==6
C4_shape = C.shape[:-2] + (3,3,3,3)
_C = C
elif C.ndim==2:
C4_shape = (1,3,3,3,3)
_C = C[None,...] # add batch dimension
C4 = np.zeros(C4_shape)
_Mandel_to_cart_4_inplace(_C, C4)
if C.ndim==2:
C4 = C4[0, ...] # remove batch dimension
return C4
# %% Young's modulus in a specific direction
[docs]
def Youngs_modulus(S : np.ndarray, d : np.ndarray) -> np.ndarray:
"""Calculate Young's modulus in a specific direction
Args:
S (np.ndarray): Compliance tensor (3,3,3,3)
d (np.ndarray): Direction vector (...,3)
Returns:
np.ndarray: _description_
"""
assert S.shape==(3,3,3,3)
if d.ndim==1:
assert np.allclose(np.linalg.norm(d), 1.0)
Sc = np.einsum('i,j,k,l,ijkl',d,d,d,d,S)
return 1/Sc
elif d.ndim==2:
assert d.shape[1]==3
assert np.allclose(np.linalg.norm(d, axis=1), np.ones(d.shape[0]))
Sc = np.einsum('ai,aj,ak,al,ijkl->a',d,d,d,d,S)
return 1/Sc
# %% Scaling exponents
def scaling_law_fit(youngs_moduli: dict):
x = list(youngs_moduli.keys())
y = np.row_stack([youngs_moduli[rel_dens] for rel_dens in x])
x = np.log10(x)
y = np.log10(y)
fit = np.polyfit(x, y, deg=1)
exponents = fit[0,:]
constants = 10**(fit[1,:])
return exponents, constants
# %% FUNCTIONS SHARED BETWEEN NUMPY AND TORCH ###
def _cart_4_tensor_to_Mandel_inplace(C4: Union[np.ndarray, Tensor], C_2: Union[np.ndarray, Tensor], mask: Union[np.ndarray, Tensor]):
for i in range(6):
if i<3:
a = b = i
elif i==3:
a = 1; b = 2
elif i==4:
a = 0; b = 2
else:
a = 0; b = 1
for j in range(i,6):
if j<3:
c = d = j
elif j==3:
c = 1; d = 2
elif j==4:
c = 0; d = 2
else:
c = 0; d = 1
Cij = C4[...,a,b,c,d]
C_2[...,i,j] = Cij
C_2[...,j,i] = Cij
C_2[...,:,:] *= mask
def _Mandel_to_cart_4_inplace(C: Union[np.ndarray, Tensor], C4: Union[np.ndarray, Tensor]):
for i in range(1,7):
if i<=3:
a = b = i
elif i==4:
a = 2; b = 3
elif i==5:
a = 1; b = 3
elif i==6:
a = 1; b = 2
for j in range(i,7):
if j<=3:
c = d = j
elif j==4:
c = 2; d = 3
elif j==5:
c = 1; d = 3
elif j==6:
c = 1; d = 2
val = C[..., i-1, j-1]
if i>3:
val = val/np.sqrt(2)
if j>3:
val = val/np.sqrt(2)
C4[..., a-1, b-1, c-1, d-1] = val
C4[..., b-1, a-1, c-1, d-1] = val
C4[..., a-1, b-1, d-1, c-1] = val
C4[..., c-1, d-1, a-1, b-1] = val
C4[..., b-1, a-1, d-1, c-1] = val
C4[..., c-1, d-1, b-1, a-1] = val
C4[..., d-1, c-1, a-1, b-1] = val
C4[..., d-1, c-1, b-1, a-1] = val
# %% TORCH FUNCTIONS ###
#########################
try:
import torch
def stiffness_cart_4_to_Mandel(_C: torch.Tensor) -> torch.Tensor:
return tensor_cart_4_to_Mandel(_C)
def compliance_cart_4_to_Mandel(_S: torch.Tensor) -> torch.Tensor:
return tensor_cart_4_to_Mandel(_S)
def stiffness_Mandel_to_cart_4(_C: torch.Tensor) -> torch.Tensor:
return Mandel_to_cart_4_tensor(_C)
def compliance_Mandel_to_cart_4(_S: torch.Tensor) -> torch.Tensor:
return Mandel_to_cart_4_tensor(_S)
def tensor_cart_4_to_Mandel(_C: torch.Tensor) -> torch.Tensor:
s2 = np.sqrt(2)
mask = torch.tensor([[1,1,1,s2,s2,s2],
[1,1,1,s2,s2,s2],
[1,1,1,s2,s2,s2],
[s2,s2,s2,2,2,2],
[s2,s2,s2,2,2,2],
[s2,s2,s2,2,2,2]], device=_C.device, dtype=_C.dtype)
if _C.dim() == 4:
C = _C.unsqueeze(0)
else:
C = _C
C_2 = _C.new_zeros((C.size(0),6,6))
_cart_4_tensor_to_Mandel_inplace(C, C_2, mask)
if _C.dim() == 4:
C_2.squeeze_(0)
return C_2
def Mandel_to_cart_4_tensor(C: torch.Tensor) -> torch.Tensor:
# convert C_ij to C_abcd
if C.ndim==3:
_C = C
elif C.ndim==2:
_C = C.unsqueeze(0)
C4 = torch.zeros((_C.shape[0], 3,3,3,3), device=C.device, dtype=C.dtype)
_Mandel_to_cart_4_inplace(_C, C4)
if C.ndim==2:
C4.squeeze_(0)
return C4
def Mandel_rot_matrix_torch(Q: torch.Tensor) -> torch.Tensor:
if Q.ndim == 2:
R = torch.zeros((6,6), device=Q.device, dtype=Q.dtype)
elif Q.ndim == 3:
R = torch.zeros((Q.shape[0], 6,6), device=Q.device, dtype=Q.dtype)
Mandel_rot_matrix_inplace(Q, R)
return R
def tens_2d_to_Mandel_torch(s: torch.Tensor) -> torch.Tensor:
if s.ndim == 2:
x = torch.zeros((6), device=s.device, dtype=s.dtype)
elif s.ndim == 3:
x = torch.zeros((s.shape[0], 6), device=s.device, dtype=s.dtype)
tens_2d_to_Mandel_inplace(s, x)
return x
except ImportError:
warnings.warn('torch not imported. Some functions will not be available.')