"""
Continuum Mechanics : Tensor Transformations and Acoustic Analysis Module.
This module provides a suite of functions to perform pull-back and push-forward
operations between different stress measures (Cauchy, PK1, PK2) and their
associated elasticity tangents. It also includes tools for computing acoustic
tensors and wave speeds in finite strain frameworks.
Supported measures:
- Deformation gradients (F) and Cauchy-Green tensors (C, b).
- Stress tensors: Cauchy (s), First Piola-Kirchhoff (P), Second Piola-Kirchhoff (S).
- Elasticity tensors: Material (C), Spatial (c), and Mixed (dPdF).
Compatibility:
All functions (but the acoustic tensors) are vectorized.
They support single tensors of shape (dim, dim)
or batches of tensors of shape (..., dim, dim) using NumPy broadcasting
via the `smallsize_tensor` backend.
"""
# pylint: disable = C0103
import numpy as np
from typing import Optional
from numba import guvectorize, float64, float32, objmode
import numba as nb
from pyxfem.core.config import FASTMATH, FloatArray
from pyxfem.core.smallsize_tensor import metric, sqrt22, trans, inv, det, \
det_invt, cofac, cross_ikjl, Id2, eigval22sym, AAT22, sqrt22_kernel
[docs]
def rightCauchyGreen(F):
"""
Compute the Right Cauchy-Green deformation tensor C = F^T @ F.
Args:
F (ndarray): Deformation gradient tensor(s) of shape (..., dim, dim).
Returns:
ndarray: The material deformation tensor(s) C.
"""
return metric(F)
[docs]
def leftCauchyGreen(F:FloatArray, out:FloatArray|None = None) -> FloatArray:
"""
Compute the Left Cauchy-Green deformation tensor b = F @ F^T.
Args:
F (ndarray): Deformation gradient tensor(s) of shape (..., dim, dim).
Returns:
ndarray: The spatial deformation tensor(s) b.
"""
if out is None :
out = np.empty_like(F)
if F.shape[-2:] == (2,2):
return AAT22(F, out)
out[:] = F@trans(F)
return out
[docs]
def PKII_to_PKI(PKII, F):
"""
Transform Second Piola-Kirchhoff stress (S) to First Piola-Kirchhoff stress (P).
Formula: P = F @ S
Args:
PKII (ndarray): PK2 stress tensor(s) of shape (..., dim, dim).
F (ndarray): Deformation gradient tensor(s) of shape (..., dim, dim).
Returns:
ndarray: First Piola-Kirchhoff stress tensor(s).
"""
return F@PKII
[docs]
def PKI_to_PKII(PKI, F):
"""
Transform First Piola-Kirchhoff stress (P) to Second Piola-Kirchhoff stress (S).
Formula: S = inv(F) @ P
Args:
PKI (ndarray): PK1 stress tensor(s) of shape (..., dim, dim).
F (ndarray): Deformation gradient tensor(s) of shape (..., dim, dim).
Returns:
ndarray: Second Piola-Kirchhoff stress tensor(s).
"""
return inv(F)@PKI
[docs]
def PKI_to_Cauchy(PKI, F):
"""
Push-forward First Piola-Kirchhoff stress (P) to Cauchy stress (sigma).
Formula: sigma = (1/J) * P @ F^T
Args:
PKI (ndarray): PK1 stress tensor(s) of shape (..., dim, dim).
F (ndarray): Deformation gradient tensor(s) of shape (..., dim, dim).
Returns:
ndarray: Cauchy stress tensor(s).
"""
J = det(F)[..., np.newaxis, np.newaxis]
return PKI@trans(F)/J
[docs]
def Cauchy_to_PKI(Cauchy, F):
"""
Pull-back Cauchy stress (sigma) to First Piola-Kirchhoff stress (P).
Formula: P = J * sigma @ F^-T = sigma @ cofac(F)
Args:
Cauchy (ndarray): Cauchy stress tensor(s) of shape (..., dim, dim).
F (ndarray): Deformation gradient tensor(s) of shape (..., dim, dim).
Returns:
ndarray: First Piola-Kirchhoff stress tensor(s).
"""
#return det(F)*Cauchy@invt(F)
print("## ", Cauchy.shape, F.shape, cofac(F).shape)
return Cauchy@cofac(F)
[docs]
def PKII_to_Cauchy(PKII, F):
"""
Push-forward Second Piola-Kirchhoff stress (S) to Cauchy stress (sigma).
Formula: sigma = (1/J) * F @ S @ F^T
Args:
PKII (ndarray): PK2 stress tensor(s) of shape (..., dim, dim).
F (ndarray): Deformation gradient tensor(s) of shape (..., dim, dim).
Returns:
ndarray: Cauchy stress tensor(s).
"""
J = det(F)[..., np.newaxis, np.newaxis]
return F@PKII@trans(F)/J
[docs]
def Cauchy_to_PKII(Cauchy, F):
"""
Pull-back Cauchy stress (sigma) to Second Piola-Kirchhoff stress (S).
Formula: S = J * F^-1 @ sigma @ F^-T
Args:
Cauchy (ndarray): Cauchy stress tensor(s) of shape (..., dim, dim).
F (ndarray): Deformation gradient tensor(s) of shape (..., dim, dim).
Returns:
ndarray: Second Piola-Kirchhoff stress tensor(s).
"""
J, FmT = det_invt(F)
J_view = np.reshape(J, (-1, 1, 1)) if np.ndim(J) > 0 else J
return J_view*trans(FmT)@Cauchy@FmT
[docs]
def dPKIIdE_to_dPKIdF(dPKIIdE, PKII, F):
"""
Transform the material tangent (dS/dE) to the mixed tangent (dP/dF).
Args:
dPKIIdE (ndarray): Fourth-order material elasticity tensor(s).
PKII (ndarray): PK2 stress tensor(s).
F (ndarray): Deformation gradient tensor(s).
Returns:
ndarray: Mixed elasticity tensor(s) dP/dF.
"""
dPKIdF = np.einsum("...iI, ...IJKL, ...kK->...iJkL", F, dPKIIdE, F)
dim = F.shape[-1]
Id = Id2[dim].reshape( (1,)*(F.ndim-2) + (dim, dim))
dPKIdF += cross_ikjl(Id, PKII)
return dPKIdF
[docs]
def dPKIdF_to_dPKIIdE(dPKIdF, PKI, F):
"""
Transform the mixed tangent (dP/dF) to the material tangent (dS/dE).
Args:
dPKIdF (ndarray): Mixed elasticity tensor(s) dP/dF.
PKI (ndarray): PK1 stress tensor(s).
F (ndarray): Deformation gradient tensor(s).
Returns:
ndarray: Material elasticity tensor(s) dS/dE.
"""
iF = inv(F)
PkII = iF@PKI
dim = F.shape[-1]
Id = Id2[dim].reshape( (1,)*(F.ndim-2) + (dim, dim))
tmp = dPKIdF - cross_ikjl(Id, PkII)
dPKIIdE = np.einsum("...Ii, ...iJkL, ...Kk->...IJKL", iF, tmp, iF)
return dPKIIdE
[docs]
def dPKIIdE_to_c(dPKIIdE, F):
"""
Push-forward the material tangent (dS/dE) to the spatial elasticity tensor (c).
Args:
dPKIIdE (ndarray): Material elasticity tensor(s) (dPKII/dE).
F (ndarray): Deformation gradient tensor(s).
Returns:
ndarray: Spatial elasticity tensor(s) c.
"""
J = det(F).reshape((-1,1,1,1,1))
c = np.einsum("...iI, ...jJ, ...kK, ...lL, ...IJKL -> ...ijkl",
F,F,F,F,dPKIIdE)/J
return c
[docs]
def c_to_dPKIIdE(c, F):
"""
Pull-back the spatial elasticity tensor (c) to the material tangent (dS/dE).
Args:
c (ndarray): Spatial elasticity tensor(s).
F (ndarray): Deformation gradient tensor(s).
Returns:
ndarray: Material elasticity tensor(s) dS/dE.
"""
J = det(F).reshape((-1,1,1,1,1))
iF = inv(F)
C = J*np.einsum("...Ii, ...Jj, ...Kk, ...Ll, ...ijkl -> ...IJKL",
iF,iF,iF,iF,c)
return C
[docs]
def dPKIdF_to_c(dPKIdF, PKI, F):
"""
Convert the mixed tangent (dP/dF) to the spatial elasticity tensor (c).
Args:
dPKIdF (ndarray): Mixed elasticity tensor(s).
PKI (ndarray): PK1 stress tensor(s).
F (ndarray): Deformation gradient tensor(s).
Returns:
ndarray: Spatial elasticity tensor(s) c.
"""
J = det(F)
dim = F.shape[-1]
s = PKI@trans(F)/J.reshape((-1,1,1))
Id = Id2[dim].reshape( (1,)*(F.ndim-2) + (dim, dim))
c = np.einsum("...jJ, ...iJkL, ...lL->...ijkl", F, dPKIdF , F)/J.reshape((-1,1,1,1,1))
c -= cross_ikjl(Id, s)
return c
[docs]
def c_dPKIdF(c, s, F):
"""
Convert the spatial elasticity tensor (c) to the mixed tangent (dP/dF).
Args:
c (ndarray): Spatial elasticity tensor(s).
s (ndarray): Cauchy stress tensor(s).
F (ndarray): Deformation gradient tensor(s).
Returns:
ndarray: Mixed elasticity tensor(s) dP/dF.
"""
J = det(F)
dim = F.shape[-1]
Id = Id2[dim].reshape( (1,)*(F.ndim-2) + (dim, dim))
tmp = J.reshape((-1,1,1,1,1))*(c + cross_ikjl(Id, s))
iF = inv(F)
dPdF = np.einsum("...Jj, ...Ll, ...ijkl -> ...iJkL", iF,iF, tmp)
return dPdF
[docs]
def acoustic_tensor_mixte(dPdF, N, J):
"""
Compute the acoustic tensor using the mixed tangent dP/dF.
Args:
dPdF (ndarray): Mixed elasticity tensor(s).
N (ndarray): normal vector(s) in reference config (..., dim).
N is the piola pull back of a unit normal vector in current configuration
N = (1/J) F^T n with |n| = 1
J (float/ndarray): Determinant of F.
Returns:
ndarray: The acoustic tensor(s) of shape (..., dim, dim).
Qik = J AiJkL NJ N L
"""
return np.einsum( "I,iIjJ,J->ij", N, dPdF, N)*J
[docs]
def acoustic_tensor_space_v0(c, s, n):
"""
Compute the acoustic tensor using spatial variables (c, sigma).
Args:
c (ndarray): Spatial elasticity tensor(s).
s (ndarray): Cauchy stress tensor(s).
n (ndarray): Direction vector(s) in current config (..., dim).
Returns:
ndarray: The spatial acoustic tensor(s).
"""
l = len(n.shape)-1
dim = n.shape[-1]
return np.einsum( "...j,...ijkl,...l-> ...ik", n, c, n) + np.einsum("...j,...jl,...l,...ik->...ik", n,s,n, Id2[dim].reshape(l*(1,)+(2,2) )) #n@(s@n)*Id2[dim]
[docs]
def acoustic_tensor_space_v1(c, s, n, out=None):
""" Version déroulée spécifiquement pour la 2D """
if n.shape[-1] != 2:
raise NotImplementedError()
if out is None:
out = np.empty(n.shape+(2,))
#out = np.empty_like(s)
n0, n1 = n[..., 0], n[..., 1]
# --- Calcul de la contrainte normale (scalaire) ---
# sigma_n = n0*(s00*n0 + s01*n1) + n1*(s10*n0 + s11*n1)
# Note : On utilise des vues temporaires, Python gère ça très vite
sigma_n = n0 * (s[..., 0, 0]*n0 + s[..., 0, 1]*n1) + \
n1 * (s[..., 1, 0]*n0 + s[..., 1, 1]*n1)
# --- Calcul de A00 ---
out[..., 0, 0] = (c[..., 0, 0, 0, 0] * n0 * n0 +
c[..., 0, 1, 0, 0] * n1 * n0 +
c[..., 0, 0, 0, 1] * n0 * n1 +
c[..., 0, 1, 0, 1] * n1 * n1) + sigma_n
# --- Calcul de A01 ---
out[..., 0, 1] = (c[..., 0, 0, 1, 0] * n0 * n0 +
c[..., 0, 1, 1, 0] * n1 * n0 +
c[..., 0, 0, 1, 1] * n0 * n1 +
c[..., 0, 1, 1, 1] * n1 * n1)
# --- Calcul de A10 ---
out[..., 1, 0] = (c[..., 1, 0, 0, 0] * n0 * n0 +
c[..., 1, 1, 0, 0] * n1 * n0 +
c[..., 1, 0, 0, 1] * n0 * n1 +
c[..., 1, 1, 0, 1] * n1 * n1)
# --- Calcul de A11 ---
out[..., 1, 1] = (c[..., 1, 0, 1, 0] * n0 * n0 +
c[..., 1, 1, 1, 0] * n1 * n0 +
c[..., 1, 0, 1, 1] * n0 * n1 +
c[..., 1, 1, 1, 1] * n1 * n1) + sigma_n
return out
[docs]
def acoustic_tensor_space(c: np.ndarray, s: np.ndarray, n: np.ndarray, out=None) -> np.ndarray:
# ...
if out is None:
batch_out = np.broadcast_shapes(c.shape[:-4], s.shape[:-2], n.shape[:-1])
out_dtype = n.dtype
out = np.empty(batch_out + (2, 2), dtype=out_dtype)
_acoustic_tensor_2d_guvec(c, s, n, out=out)
return out
#acoustic_tensor_space = acoustic_tensor_space
[docs]
@nb.njit(cache=True, fastmath = FASTMATH)
def acoustic_tensor_2d_kernel_njit(c, s, n, out):
n0 = n[0]
n1 = n[1]
sigma_n = n0 * (s[0, 0]*n0 + s[0, 1]*n1) + \
n1 * (s[1, 0]*n0 + s[1, 1]*n1)
out[0, 0] = (c[0, 0, 0, 0] * n0 * n0 + c[0, 1, 0, 0] * n1 * n0 +
c[0, 0, 0, 1] * n0 * n1 + c[0, 1, 0, 1] * n1 * n1) + sigma_n
out[0, 1] = (c[0, 0, 1, 0] * n0 * n0 + c[0, 1, 1, 0] * n1 * n0 +
c[0, 0, 1, 1] * n0 * n1 + c[0, 1, 1, 1] * n1 * n1)
out[1, 0] = (c[1, 0, 0, 0] * n0 * n0 + c[1, 1, 0, 0] * n1 * n0 +
c[1, 0, 0, 1] * n0 * n1 + c[1, 1, 0, 1] * n1 * n1)
out[1, 1] = (c[1, 0, 1, 0] * n0 * n0 + c[1, 1, 1, 0] * n1 * n0 +
c[1, 0, 1, 1] * n0 * n1 + c[1, 1, 1, 1] * n1 * n1) + sigma_n
@guvectorize([
(float64[:,:,:,:], float64[:,:], float64[:], float64[:,:]),
(float32[:,:,:,:], float32[:,:], float32[:], float32[:,:])
], '(m,m,m,m),(m,m),(m)->(m,m)', nopython=True, cache=True)
def _acoustic_tensor_2d_guvec(c, s, n, out):
acoustic_tensor_2d_kernel_njit(c, s, n, out)
[docs]
def acoustic_impedance_tensor(rho, c, s, n, outZ = None):
""" compute the acoustic tensor impedance in direction n
c is the 4 order tangent spacial material tensor
s is the cauchy stress
n is a direction (norm 1 ..)
rho is the density
return OutZ : acoustic tensor in direction n
"""
#print(rho.shape, c.shape, s.shape ,n.shape, outZ.shape)
#print(c.shape, s.shape, n.shape)
if outZ is None :
batch_out = np.broadcast_shapes(c.shape[:-4], s.shape[:-2], n.shape[:-1])
out_dtype = n.dtype
outZ = np.empty(batch_out + (2, 2), dtype=out_dtype)
acoustic_tensor_space(c,s,n, outZ)
#print("G", G.shape)
outZ *= rho[..., None, None]
return sqrt22(outZ, outZ)
#Z = sqrt22(rho.reshape(rho.shape+(1,1))*A)
#return outZ
[docs]
@nb.njit(parallel = True, cache = True, fastmath = FASTMATH)
def acoustic_impedance_tensor_pnjit(rho, c, s, n, out):
#print(rho.shape, c.shape, s.shape ,n.shape, out.shape)
N = rho.shape[0]
assert(N == c.shape[0] and s.shape[0] == N and n.shape[0] == N and out.shape[0] == N)
assert( len(n.shape) == 3) # typically we have more than one normal for each
for j in nb.prange(N):
jrho = rho[j, 0]
jc = c[j, 0]
js = s[j, 0]
jvn = n[j]
jvout = out[j]
for r in range(n.shape[1]) :
jrout = jvout[r]
acoustic_tensor_2d_kernel_njit(jc, js, jvn[r], jrout )
jrout[0,0]*=jrho
jrout[0,1]*=jrho
jrout[1,0]*=jrho
jrout[1,1]*=jrho
sqrt22_kernel(jrout, jrout)
return out
[docs]
def acoustic_tensor_material(C, S, N, F):
"""
Compute the acoustic tensor using material variables (C, S).
Args:
C (ndarray): Material elasticity tensor(s).
S (ndarray): PK2 stress tensor(s).
N (ndarray): normal vector(s) in reference config (..., dim).
N is the piola pull back of a unit normal vector in current configuration
N = (1/J) F^T n with |n| = 1
F (ndarray): Deformation gradient tensor(s).
Returns:
ndarray: The material acoustic tensor(s).
"""
dim = len(N)
J =det(F)
A = J*(np.einsum( "iI,J,IJKL,kK,L -> ik", F, N, C, F, N ) + N.dot(S@N)*Id2[dim])
return A
[docs]
def wavespeed(A, rho):
"""
Calculate wave speeds from the acoustic tensor and density.
For now only work for 2x2 acoustic tensor ....
Args:
A (ndarray): Acoustic tensor(s).
rho (float/ndarray): Material density.
Returns:
ndarray: Wave speeds (sqrt of eigenvalues of A/rho).
"""
if A.shape[-1] != 2 :
raise NotImplementedError("Wavespeed only works for acoustic tensor in 2d space")
eigs = eigval22sym(A)
np.sqrt(eigs/rho)
return eigs
[docs]
def test_acoustic_tensor():
""" Testing Acoustic tensor and wave seed computation """
# pylint: disable=import-outside-toplevel
from pyxfem.solidmechanics.hyperelasticity import CompressibleNeoHookeen
from pyxfem.solidmechanics.elasticity import MaterialHypothesis
law = CompressibleNeoHookeen(lamb = 3., mu= 7., hyp=MaterialHypothesis.PLANESTRAIN)
rho0 = 1.
#F = np.array([[1.,0.], [0., 1.]])
F = np.array([[3.,1.], [0.8, 2]])
n = np.array([1.,0.8])
n /= np.linalg.norm(n)
dPdF = law.dPKIdF(F)
J = det(F)
N = trans(F)@n/J
A1 = acoustic_tensor_mixte(dPdF, N, J)
#print("A1", A1)
b = leftCauchyGreen(F)
s = law.Cauchy(b)
c = law.spacial_elasticty_tensor(b)
A2 = acoustic_tensor_space(c, s, n)
S = law.PKII(rightCauchyGreen(F))
C = law.dPKIIdE(rightCauchyGreen(F))
A3 = acoustic_tensor_material(C, S, N, F)
print("A1", A1)
print("A2", A2)
print("A3", A3)
J = det(F)
rho = rho0/J
print('wavespeed', wavespeed(A1, rho))
print('wavespeed', wavespeed(A2, rho))
print('wavespeed', wavespeed(A3, rho))
print("test Acoustic tensor for multiple direction")
ns = np.array([[1.,0.],[0,1], [np.sqrt(2)/2, np.sqrt(2)/2], [n[0], n[1]]])
As = acoustic_tensor_space(c[np.newaxis], s[np.newaxis], ns)
print(As.shape)
print(As)
c = np.stack((c, 2*c), axis=0)
s = np.stack((s, 2*s), axis=0)
As = acoustic_tensor_space(c[:,np.newaxis], s[:,np.newaxis], ns[np.newaxis])
print("2 states, 4 vectors ")
print(As.shape)
print(As)
rho = np.array([rho, rho])
print('wavespeed', wavespeed(As, rho[np.newaxis,:]))
if __name__ == "__main__":
test_transport_formula()
test_acoustic_tensor()