#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
This module is dedicated to hyperelasticity
Created on Thu April 2024 14:09:29 2024
@author: nchevaug
"""
from __future__ import annotations
import re
import time
from typing import Callable, Tuple, cast
import numpy as np
from numba import guvectorize, float32, float64, njit, prange
from scipy.special import lambertw as sp_lambertw
from pyxfem.core import smallsize_tensor
from pyxfem.core.config import FloatArray, FloatVars, npfloat
from pyxfem.core.smallsize_tensor import det_invt, cross_ijkl, df_num,\
random_lintransform, random_lintransform2dto3d,\
metric, det, trace, trace2, inv, invt, dinvt, check_relerr, \
Sy4, Sp4
from pyxfem.core.exec_tools import Profile
from pyxfem.solidmechanics.elasticity import MaterialHypothesis
#import pytest
#from sym_hyperelasticity import II
PLANESTRAIN = MaterialHypothesis.PLANESTRAIN
PLANESTRESS = MaterialHypothesis.PLANESTRESS
DIM3 = MaterialHypothesis.DIM3
HYP2DIM = {PLANESTRESS:2, PLANESTRAIN:2, DIM3:3}
import numpy as np
[docs]
def real_lambertw(z):
"""
Computes the principal branch of the Lambert W function for real, positive z
using vectorized Halley's method.
"""
# A robust initial guess for z >= 0
w = np.log(1.0 + z)
# Halley's method usually converges to machine precision in 3-4 iterations
for _ in range(5):
ew = np.exp(w)
w1 = w + 1.0
f = w * ew - z
# Halley step
step = f / (ew * w1 - ((w + 2.0) * f) / (2.0 * w1))
w -= step
# Optional: break early if max error is small enough,
# but for 5 fixed iterations, it's often faster to just loop.
if np.max(np.abs(step)) < 1e-12:
break
return w
#lambertw = lambda x : np.real(sp_labertw(x))
lambertw = real_lambertw
[docs]
def vecroots(f:Callable, fprime:Callable, x0:FloatVars,
atol:float=1.e-20, rtol=1.e-12, maxiter:int = 100, verbose = False)->FloatVars:
""" compute roots of an array of independant equations using newton's
method with line search so that all x stay > 0"""
r = f(x0)
if verbose:
print(0, x0, np.abs(r), 1.)
err0 = np.max(np.abs(r))
err = err0
if err0 < atol:
return x0
it = 0
while it < maxiter:
alpha=1.
j = fprime(x0)
s = -r/j
x = x0 + s
while np.any(x<=0.):
alpha = 0.5*alpha
x = x0 + alpha*s
r = f(x)
while np.max(np.abs(r)) > err:
alpha = 0.5*alpha
x = x0 + alpha*s
r = f(x)
err = np.max(np.abs(r))
if err < atol or err < rtol*err0 :
return x
x0 = x
it += 1
if verbose:
print(it, x0, np.abs(r), alpha)
msg = f"Failed to converge after {maxiter:d} iterations"
raise RuntimeError(msg)
[docs]
class HyperElasticMaterial():
""" Virtual base class for HyperElastic material """
# pylint: disable=C0103
# pylint: disable=W2301
def __init__(self, hyp:MaterialHypothesis=DIM3, rho0 = 1.):
self.hyp, self.rho = hyp, rho0
[docs]
def potentialF(self, F:FloatArray )-> FloatVars:
""" Compute the Elastic Potential as a Function of F"""
raise NotImplementedError()
[docs]
def potentialC(self, C:FloatArray )-> FloatVars:
""" Compute the Elastic Potential as a Function of C"""
raise NotImplementedError()
[docs]
def PKI(self, F:FloatArray)->FloatArray:
""" Compute Piola Kirchhoff I Stress tensor as a function of F """
raise NotImplementedError()
[docs]
def PKII(self, C:FloatArray)->FloatArray:
""" Compute Piola Kirchhoff I Stress tensor as a function of C = F^TF """
raise NotImplementedError()
[docs]
def Cauchy(self, b:FloatArray)->FloatArray:
""" Compute Cauchy Stress tensor as a function of b = FF^T """
raise NotImplementedError()
[docs]
def dPKIdF(self, F:FloatArray)->FloatArray:
""" Compute the deriviatives of the Piola Kirchhoff I Stress tensor
with regard to F
"""
raise NotImplementedError()
[docs]
def dPKIIdC(self, C:FloatArray)->FloatArray:
""" Compute the deriviatives of the Piola Kirchhoff II Stress tensor
with regard to C
"""
raise NotImplementedError()
[docs]
def dPKIIdE(self, C:FloatArray)->FloatArray:
""" Compute the deriviatives of the Piola Kirchhoff II Stress tensor
with regard to C
"""
return 2*self.dPKIIdC(C)
[docs]
def spacial_elasticty_tensor(self, b:FloatArray)->FloatArray:
""" Compute the spacial elasticity tensor (1/J FFFF dPKIIdE)
"""
raise NotImplementedError()
# class SimpleLaw0(HyperElasticMaterial):
# """ A very simple hyper elastic law. Mostly to check derivatives of I_C """
# # pylint: disable=C0103
# def __init__(self):
# super().__init__(PLANESTRAIN)
# def potential(self, F:FloatArray)-> FloatArray|float:
# m = metric(F)
# I_C = 1.+trace(m)
# return 0.5*I_C
# def PKI(self, F:FloatArray)->FloatArray:
# return F
# def dPKIdF(self, F:FloatArray)->FloatArray:
# tshape = F.shape[-2:]
# nbfirstaxis = F.ndim-2
# I40 = np.eye(tshape[0]*tshape[1]).reshape(tshape+tshape)
# return I40[(np.newaxis,)*nbfirstaxis]
# class SimpleLaw1(HyperElasticMaterial):
# """ A very simple hyper elastic law. Mostly to check derivatives of log(J) """
# # pylint: disable=C0103
# def __init__(self):
# super().__init__(PLANESTRAIN)
# def potential(self, F:FloatArray)-> FloatArray|float:
# return 0.5*np.log(det(metric(F)))
# def PKI(self, F:FloatArray)->FloatArray:
# return invt(F)
# def dPKIdF(self, F:FloatArray)->FloatArray:
# return dinvt(F)
# class SimpleLaw2(HyperElasticMaterial):
# """ A very simple hyper elastic law. Mostly to check derivatives of J """
# # pylint: disable=C0103
# def __init__(self):
# super().__init__(PLANESTRAIN)
# def potential(self, F:FloatArray)-> FloatArray|float:
# return np.sqrt(det(metric(F)))
# def PKI(self, F:FloatArray)->FloatArray:
# J, invFT = det_invt(F)
# return np.atleast_1d(J)[..., np.newaxis, np.newaxis] * invFT
# def dPKIdF(self, F:FloatArray)->FloatArray:
# J, invFT = det_invt(F)
# J_11 = np.atleast_1d(J)[..., np.newaxis, np.newaxis]
# PKI = J_11 * invFT
# J_1111 = J_11[..., np.newaxis, np.newaxis]
# return cross_ijkl(PKI, invFT) + J_1111*dinvt(F)
[docs]
@njit(inline='always', fastmath=True)
def fill_nh_sigma_2d(J, b00, b11, b01, mu, lamb, out_sigma):
""" Remplit la contrainte de Cauchy 2D en place """
S_mu = mu / J
S_lamb = (lamb * np.log(J)) / J
out_sigma[0, 0] = S_mu * (b00 - 1.0) + S_lamb
out_sigma[1, 1] = S_mu * (b11 - 1.0) + S_lamb
out_sigma[0, 1] = S_mu * b01
out_sigma[1, 0] = out_sigma[0, 1]
[docs]
@njit(inline='always', fastmath=True)
def fill_nh_c_2d(J, mu, lamb, out_c):
""" Remplit le tenseur d'élasticité spatial 2x2x2x2 en place """
S1 = lamb / J
S2 = 2.0 * (mu - lamb * np.log(J)) / J
S2_half = 0.5 * S2
out_c.fill(0.0)
out_c[0, 0, 0, 0] = S1 + S2
out_c[1, 1, 1, 1] = S1 + S2
out_c[0, 0, 1, 1] = S1
out_c[1, 1, 0, 0] = S1
out_c[0, 1, 0, 1] = S2_half
out_c[0, 1, 1, 0] = S2_half
out_c[1, 0, 0, 1] = S2_half
out_c[1, 0, 1, 0] = S2_half
@guvectorize([
(float64[:,:], float64, float64, float64[:,:,:,:]),
(float32[:,:], float32, float32, float32[:,:,:,:])
], '(m,m),(),()->(m,m,m,m)', nopython=True, cache=True)
def _NH_spacial_elasticity_plane_strain_guvec(b, lamb, mu, out):
det_b = b[0, 0] * b[1, 1] - b[0, 1] * b[1, 0]
J = np.sqrt(det_b)
fill_nh_c_2d(J, mu, lamb, out)
@guvectorize([
('float64[:,:], float64, float64, float64[:,:]'),
('float32[:,:], float32, float32, float32[:,:]')
], '(m,m),(),()->(m,m)', nopython=True, cache=True)
def _NH_cauchy_plane_strain_guvec(b, lamb, mu, out_sigma):
"""
Calcule la contrainte de Cauchy 2D (Déformation Plane) pour un
matériau Néo-Hookéen Compressible.
Inputs:
b : Tenseur de Cauchy-Green gauche (2x2)
lamb : Paramètre de Lamé lambda (scalaire)
mu : Module de cisaillement mu (scalaire)
Output:
out_sigma : Tenseur des contraintes de Cauchy (2x2)
"""
# 1. Calcul du déterminant et du Jacobien (J = sqrt(det(b)))
det_b = b[0, 0] * b[1, 1] - b[0, 1] * b[1, 0]
J = np.sqrt(det_b)
fill_nh_sigma_2d(J, b[0,0], b[1,1], b[0,1], mu, lamb, out_sigma)
[docs]
@njit(inline='always', fastmath=True, cache=True)
def neo_hookean_law_local(F, props, out_sigma, out_c):
"""
Calcule la contrainte de Cauchy et le tenseur d'élasticité spatial 2D.
Zéro allocation, remplissage en place.
"""
mu = props[0]
lamb = props[1]
# --- Cinématique ---
# J = det(F)
detF = F[0, 0] * F[1, 1] - F[0, 1] * F[1, 0]
J = detF
# b = F * F^T (Tenseur de Cauchy-Green gauche)
b00 = F[0, 0]*F[0, 0] + F[0, 1]*F[0, 1]
b01 = F[0, 0]*F[1, 0] + F[0, 1]*F[1, 1]
b11 = F[1, 0]*F[1, 0] + F[1, 1]*F[1, 1]
fill_nh_sigma_2d(J, b00, b11, b01, mu, lamb, out_sigma)
fill_nh_c_2d(J, mu, lamb, out_c)
[docs]
class CompressibleNeoHookeen(HyperElasticMaterial):
""" Implementation of the NeoHookeen Hyperelastic material """
# pylint: disable=C0103
def __init__(self,*, # pylint: disable=R0913
E:float|None=None,
nu:float|None=None,
lamb:float|None=None,
mu:float|None=None,
rho0:float|None = None,
hyp:MaterialHypothesis=DIM3):
super().__init__(hyp)
if isinstance(E, float) and isinstance(nu, float):
self.E = E
self.nu = nu
self.lamb = E*nu/(1.+nu)/(1.-2*nu)
self.mu = E/2./(1+nu)
elif isinstance(lamb, float) and isinstance(mu, float):
self.lamb = lamb
self.mu = mu
self.E = mu*(3.0*lamb + 2.0*mu)/(lamb + mu)
self.nu = 0.5*lamb/(lamb + mu)
else:
raise ValueError("E and nu or lambda and mu must be set (exclusive)")
self.rho0 = rho0
def _checktensorsize(self, F:FloatArray):
if F.shape[-2:] == (3,3) and self.hyp == DIM3:
return
if F.shape[-2:] in [(2,2), (3,2)]:
if self.hyp in [PLANESTRAIN, PLANESTRESS]:
return
raise ValueError("invalid shape for F. F.shape must be (...2,2) or (...,3,3)or (..., 3,2)")
[docs]
def l3_planestress(self, J2d:FloatVars)->FloatVars:
""" return the principal elongation lambda_3 in the direction orthogonal to the plane """
lamb, mu = (self.lamb, self.mu)
a = mu
b = lamb/2
c = lamb*np.log(J2d) -mu
#return np.sqrt(b*np.real(lambertw( a* np.exp(-c/b)/b)) /a)
return np.sqrt(b*lambertw( a* np.exp(-c/b)/b) /a)
# def resl3(l3):
# return lamb*np.log(J2d*l3) + mu*(l3*l3 - 1)
# def jacl3(l3):
# return lamb/l3 + 2*mu*l3
# if isinstance(J2d, float) or isinstance(J2d, int):
# x0 = 1.
# else:
# x0 = np.ones(J2d.shape)
# l3 = vecroots(resl3, jacl3, x0, atol=1.e-20, rtol=1.e-12, maxiter=100)
#return l3
# def _dl3dF2d_planestress(self, l3, invFT2d):
# lamb = self.lamb
# mu = self.mu
# return -l3/(1+2*mu*l3**2/lamb)*invFT2d
# def dl3F2d_planestress(self, F2d):
# """ Compute the derivative of the third principal strain as a fonction of
# the plane part of F (Assuming Plan Stress)"""
# J2d, invFT2d = det_invt(F2d)
# l3 = self.l3_planestress(J2d)
# return self._dl3dF2d_planestress(l3, invFT2d)
def _potential(self, J:FloatVars, I_C:FloatVars)-> FloatVars:
logJ = np.log(J)
return self.mu*(0.5*(I_C-3)- logJ) + self.lamb*logJ**2/2
def _potential_3d(self, F:FloatArray)-> FloatVars:
J = det(F)
I_C = trace2(F)
return self._potential(J, I_C)
def _potential_planestrain(self, F2d:FloatArray)-> FloatVars:
if F2d.shape[-2:] == (3,2):
m = metric(F2d)
J = np.sqrt(det(m))
I_C = 1.+trace(m)
elif F2d.shape[-2:] == (2,2):
J = det(F2d)
I_C = trace2(F2d) +1
else :
raise ValueError('_potential_planestrain, F2d must be of shape (...,3,2) or (...,2,2)')
return self._potential(J, I_C)
def _potential_planestress(self, F2d:FloatArray)-> FloatVars:
if F2d.shape[-2:] == (3,2):
m = metric(F2d)
J2d = np.sqrt(det(m))
I_C2d = trace(m)
elif F2d.shape[-2:] == (2,2):
J2d = det(F2d)
I_C2d = trace2(F2d)
else :
raise ValueError('_potential_planestress, F2d must be of shape (...,3,2) or (...,2,2)')
l3 = self.l3_planestress(J2d)
J = l3*J2d
I_C = l3**2+I_C2d
return self._potential(J, I_C)
def _PKI_3d(self, F:FloatArray)->FloatArray:
'''Compute Piola Kirchhoff I Stress tensor as a function of F: 3D case (and Plane Strain)'''
if F.shape[-2:] not in [(3,3), (3,2), (2,2)]:
raise ValueError("F must be a nd.array of shape[-2:] in [(3,3), (3,2), (2,2)")
J, invFT = det_invt(F)
#print("### ", J, F, invFT)
return self._PKI_exp(J, F, invFT)
def _PKI_exp_v0(self, J, F, invFT)->FloatArray:
''' Compute Piola Kirchhoff I Stress tensor as a function of J, F, invFT'''
return self.mu*(F-invFT) + self.lamb*np.log(J)[..., np.newaxis, np.newaxis]*invFT
def _PKI_exp(self, J, F, invFT) -> FloatArray:
''' Compute Piola Kirchhoff I Stress tensor as a function of J, F, invFT'''
# 1. Factor out the scalar coefficients to eliminate array subtractions
coef = (self.lamb * np.log(J) - self.mu)[..., np.newaxis, np.newaxis]
# 2. Allocate exact output memory once
ret = np.empty_like(F)
# 3. Compute in-place to bypass temporary array copies
np.multiply(F, self.mu, out=ret)
ret += coef * invFT
return ret
def _PKI_planestress(self, F2d:FloatArray)->FloatArray:
''' Compute Piola-Kirchhoff I Stress tensor as a function of F : plane stress case)'''
if F2d.shape[-2:] not in[(2,2), (3,2)]:
raise ValueError("F must be a nd.array of shape (...,3,3)")
J2d, invFT2d = det_invt(F2d)
if np.any(J2d<=0.):
raise ValueError("det(F2d) must be > 0.")
l3 = self.l3_planestress(J2d)
J = l3*J2d
return self._PKI_exp(J, F2d, invFT2d)
def _dPKIdF_exp_v0(self, c0, c1, c2, invFT)->FloatArray:
c0 = np.atleast_1d(c0)[(...,) + (np.newaxis,)*4]
c1 = np.atleast_1d(c1)[(...,) + (np.newaxis,)*4]
c2 = np.atleast_1d(c2)[(...,) + (np.newaxis,)*4]
tshape = invFT.shape[-2:]
I4_ = np.eye(tshape[0]*tshape[1]).reshape(tshape+tshape)
I4_ = I4_.reshape((1,)*(invFT.ndim-2)+I4_.shape, copy=False)
c0_I4 = c0*I4_
invFTxinvFT = cross_ijkl(invFT, invFT) # iFT x d(log(J)) = iFTxiFT
c2_invFTxinvFT = c2*invFTxinvFT
c1_dinvFTdF = -c1*(invFTxinvFT.transpose(list(range(invFTxinvFT.ndim-4)) + [-4,-1,-2,-3]))
#c1_dinvFTdF = c1*dinvt(invFT, invFT)
ret = c0_I4 + c1_dinvFTdF +c2_invFTxinvFT
if invFT.shape[:-2] == ():
return ret.squeeze()
return ret
def _dPKIdF_exp_v1(self, c0, c1, c2, invFT) -> FloatArray:
# 1. Expand scalars/arrays for 4D tensor broadcasting
c0 = np.atleast_1d(c0)[..., None, None, None, None]
c1 = np.atleast_1d(c1)[..., None, None, None, None]
c2 = np.atleast_1d(c2)[..., None, None, None, None]
# 2. Build the 4th-order identity tensor I4
d = invFT.shape[-1]
I4 = np.eye(d * d).reshape((d, d, d, d))
# 3. Allocate full memory block to match the number of quadrature points
full_shape = invFT.shape[:-2] + (d, d, d, d)
ret = np.empty(full_shape, dtype=invFT.dtype)
ret[...] = c0 * I4
# 4. Compute and subtract the c1 * dinvFTdF term in-place
ret -= c1 * (invFT[..., :, None, None, :] * invFT.swapaxes(-1, -2)[..., None, :, :, None])
# 5. Compute and add the c2 * invFTxinvFT term in-place
ret += c2 * (invFT[..., :, :, None, None] * invFT[..., None, None, :, :])
if invFT.ndim == 2:
return ret.squeeze()
return ret
def _dPKIdF_exp_v2(self, F)->FloatArray:
mu, lamb = self.mu, self.lamb
tshape = F.shape[-2:]
I4_ = np.eye(tshape[0]*tshape[1]).reshape(tshape+tshape)
I4_.shape = (1,)*(F.ndim-2)+I4_.shape
J, iFT = det_invt(F)
J_1111 = np.atleast_1d(J)[(...,) + (np.newaxis,)*4]
return mu*I4_ + lamb*cross_ijkl(iFT, iFT) + (lamb*np.log(J_1111)-mu)*dinvt(F)
def _dPKIdF_exp_v3(self, c0, c1, c2, invFT) -> FloatArray:
d = invFT.shape[-1]
I4 = np.eye(d * d).reshape((d, d, d, d))
# Expand constants for broadcasting
c0_b = np.atleast_1d(c0)[..., None, None, None, None]
c1_b = np.atleast_1d(c1)[..., None, None, None, None]
c2_b = np.atleast_1d(c2)[..., None, None, None, None]
# 1. Base allocation (c0 * I4)
full_shape = invFT.shape[:-2] + (d, d, d, d)
ret = np.empty(full_shape, dtype=invFT.dtype)
ret[...] = c0_b * I4
# 2. Term: -c1 * invFT_il * invFT_jk
# Directly computes the transposed term in C without copying memory
ret -= c1_b * np.einsum('...il,...kj->...ijkl', invFT, invFT)
# 3. Term: c2 * invFT_ij * invFT_kl
# Replaces cross_ijkl with direct C-loop
ret += c2_b * np.einsum('...ij,...kl->...ijkl', invFT, invFT)
if invFT.ndim == 2:
return ret.squeeze()
return ret
@staticmethod
@njit(parallel=True, fastmath=True, cache = True)
def _dPKIdF_numba_core(c0, c1, c2, invFT, ret):
N = invFT.shape[0]
d = invFT.shape[1]
for g in prange(N):
c0_g = c0[g]
c1_g = c1[g]
c2_g = c2[g]
for i in range(d):
for j in range(d):
for k in range(d):
for l in range(d):
val = 0.0
if i == k and j == l:
val += c0_g
val -= c1_g * invFT[g, i, l] * invFT[g, k, j]
val += c2_g * invFT[g, i, j] * invFT[g, k, l]
ret[g, i, j, k, l] = val
def _dPKIdF_exp(self, c0, c1, c2, invFT):
orig_shape = invFT.shape
d = orig_shape[-1]
if invFT.ndim > 2:
N = np.prod(orig_shape[:-2], dtype=int)
else:
N = 1
invFT_flat = invFT.reshape(N, d, d)
c0_flat = np.broadcast_to(c0, (N,)).astype(invFT.dtype)
c1_flat = np.broadcast_to(c1, (N,)).astype(invFT.dtype)
c2_flat = np.broadcast_to(c2, (N,)).astype(invFT.dtype)
ret_flat = np.empty((N, d, d, d, d), dtype=invFT.dtype)
# Call the static method via self
self._dPKIdF_numba_core(c0_flat, c1_flat, c2_flat, invFT_flat, ret_flat)
if invFT.ndim == 2:
return ret_flat.squeeze()
return ret_flat.reshape(orig_shape[:-2] + (d, d, d, d))
[docs]
def potentialF(self, F:FloatArray)-> FloatVars:
""" return the Hyperelastic potential as a function of F"""
self._checktensorsize(F)
if self.hyp == DIM3:
return self._potential_3d(F)
if self.hyp == PLANESTRAIN:
return self._potential_planestrain(F)
if self.hyp == PLANESTRESS:
return self._potential_planestress(F)
raise ValueError("MaterialHypothesis" + str(self.hyp) +"Not implemented")
[docs]
def PKI(self, F:FloatArray)->FloatArray:
''' Compute Piola Kirchhoff I Stress tensor as a function of F'''
if self.hyp in [PLANESTRAIN, DIM3]:
return self._PKI_3d(F)
return self._PKI_planestress(F)
[docs]
def dPKIdF(self, F:FloatArray)->FloatArray:
''' Compute the deriviatives of the Piola Kirchhoff I Stress tensor
with regard to F at
'''
self._checktensorsize(F)
J, invFT = det_invt(F)
#print('invFT shap', invFT.shape)
mu =self.mu
lamb = self.lamb
c0 = mu
if self.hyp in [DIM3, PLANESTRAIN]:
c1 = lamb*np.log(J)-mu
c2 = lamb
elif self.hyp == PLANESTRESS :
J2d = J
l3 = self.l3_planestress(J2d)
J = l3*J2d
c1 = lamb*np.log(J)-mu
c2 = lamb*(1.-lamb/(lamb+2*mu*l3*l3))
else:
raise ValueError(f"dPKIdF not implemented for {str(self.hyp):}")
if F.shape[-2:] == (3,2):
I4_ = np.eye(6).reshape((1,)*(F.ndim-2)+(3,2,3,2), copy = False)
c1 = np.atleast_1d(c1)[(...,)+ (np.newaxis,)*4]
c2 = np.atleast_1d(c2)[(...,)+ (np.newaxis,)*4]
return c0*I4_ + c2*cross_ijkl(invFT, invFT) + c1*dinvt(F)
return self._dPKIdF_exp(c0, c1, c2, invFT)
[docs]
def C33_planestress(self, II_C2D:FloatVars) -> FloatVars:
""" Return the square of the thichkess dilatation in plane stress Case,
knowing the secons invariant of the planar componant of C
mu(C33 -1) + kappa ln(J) = sig33
"""
mu, lamb= (self.mu, self.lamb)
a = 0.5*lamb/mu
ia = 1/a
#b = np.log(II_C2D)
# r0 = a*np.real(lambertw( np.exp(ia-b)*ia ))
# r1 = a*np.real(lambertw( ia/II_C2D*np.exp(ia))
# print("ga", r0-r1)
return a*np.real(lambertw( ia/II_C2D*np.exp(ia)))
[docs]
def dC33dII_planestress(self, II_C2D:FloatVars, C33:FloatVars|None = None) -> FloatVars:
""" Return the square of the thichkess dilatation in plane stress Case,
knowing the secons invariant of the planar componant of C """
mu, lamb= (self.mu, self.lamb)
II = II_C2D
if C33 is None :
C33 = self.C33_planestress(II)
a = cast(npfloat, 0.5*lamb/mu)
dC33dII = -a*C33/(II*(a+C33))
return dC33dII
[docs]
def potentialC(self, C:FloatArray) -> FloatArray:
""" Compute The elastic potential as a function of C"""
mu, lamb= (self.mu, self.lamb)
if C.shape[-2:] == (3,3):
I = trace(C)
logIII = np.log(np.linalg.det(C))
elif C.shape[-2:] == (2,2):
I_C2D = trace(C)
II_C2D = det(C)
if self.hyp == MaterialHypothesis.PLANESTRESS :
C33 = self.C33_planestress(II_C2D)
elif self.hyp == MaterialHypothesis.PLANESTRAIN:
C33 = 1.
else:
raise ValueError("self.hyp must be either PLANESTRESS or PLANESTRAIN")
I = I_C2D + C33
logIII = np.log(II_C2D*C33)
else :
raise ValueError("C must be of shape (..., 2,2) or (..., 3, 3)")
return mu*(I - 3 - logIII)/2. + lamb*logIII**2/8
[docs]
def PKII(self, C:FloatArray) -> FloatArray:
""" Return second Piola Kirchoff Stres as a function of C """
lamb, mu = (self.lamb, self.mu)
if C.shape[-2:] == (3,3) or \
(C.shape[-2:] == (2,2) and self.hyp == MaterialHypothesis.PLANESTRAIN) :
dim = C.shape[-1]
III = np.linalg.det(C).reshape( C.shape[:-2] + (1,1))
logIII = np.log(III)
Id = np.eye(dim).reshape( (1,)*(C.ndim-2) + (dim, dim) )
iCT = invt(C, III)
kappa = 0.5*lamb*logIII - mu
# print("fff", mu*Id[0,2,2], kappa[0], iCT[0,2,2])
elif C.shape[-2:] == (2,2) and self.hyp == MaterialHypothesis.PLANESTRESS :
Id = np.eye(2).reshape ((1,)*(C.ndim-2) + (2,2))
II_2D = det(C).reshape(C.shape[:-2]+(1,1))
iCT = invt(C, II_2D)
C33 = self.C33_planestress(II_2D)
logIII = np.log(C33*II_2D)
kappa = mu*C33*(lamb*(logIII-1) - 2*mu )/(2*mu*C33+lamb)
else :
raise NotImplementedError()
return mu*Id + kappa*iCT
[docs]
def dPKIIdC(self, C:FloatArray) -> FloatArray:
lamb, mu = (self.lamb, self.mu)
iC = inv(C)
#iCxiC = cross_ijkl(iC, iC)
#diCdC = dinv(C, iA= iC)
diCTdC = dinvt(C)
iCT = smallsize_tensor.trans(iC)
iCTxiCT = cross_ijkl(iCT,iCT)
#diCdC = dinv_sym(C, iA= iC)
#diCTdC = dinvt(C, iCT)
if C.shape[-2:] == (3,3) or \
(C.shape[-2:] == (2,2) and self.hyp == MaterialHypothesis.PLANESTRAIN) :
III = det(C).reshape( C.shape[:-2] + (1,1,1,1))
logIII = np.log(III)
# kappa = 0.5*lamb*logIII- mu
# IIIdkappadIII = 0.5*lamb
# dkappa = IIIdkappadIII
return (0.5*lamb)*iCTxiCT + (0.5*lamb*logIII- mu)*diCTdC
elif C.shape[-2:] == (2,2) and self.hyp == MaterialHypothesis.PLANESTRESS:
II_2D = det(C).reshape(C.shape[:-2]+(1,1,1,1))
C33 = self.C33_planestress(II_2D)
III = C33*II_2D
logIII = np.log(III)
a = 2*mu*C33+lamb
kappa = mu*C33*(lamb*(logIII-1) - 2*mu)/a
II_2DxdkappadII_2D = mu*lamb*C33*( a**2 - 2*mu*lamb*(C33-1) -lamb**2*(logIII) ) /a**3
dkappa = II_2DxdkappadII_2D
#II_2DxdkappadII_2D = mu*lamb*C33*( a**2 - 2*mu*lamb*(C33-1) -lamb**2*(logIII) ) /a**3
return dkappa*iCTxiCT + kappa*diCTdC
#return dkappa*iCxiC + kappa*diCdC
else :
raise ValueError("incorrect size for C (..., 2,2) or (..., (3,3))")
[docs]
def Cauchy(self, b:FloatArray, out_sigma:FloatArray|None = None)->FloatArray:
lamb, mu = (self.lamb, self.mu)
if b.shape[-2:] == (2,2) and self.hyp == MaterialHypothesis.PLANESTRAIN:
if out_sigma is None :
out_sigma = np.empty(b.shape[:-2] + (2,2))
_NH_cauchy_plane_strain_guvec(b, lamb, mu, out_sigma)
return out_sigma
if b.shape[-2:] == (3,3):
dim = b.shape[-1]
J = np.sqrt(det(b).reshape( b.shape[:-2] + (1,1)))
Id = np.eye(dim).reshape( (1,)*(b.ndim-2) + (dim, dim) )
elif b.shape[-2:] == (2,2) and self.hyp == MaterialHypothesis.PLANESTRESS :
II_2D = det(b).reshape(b.shape[:-2]+(1,1))
Id = np.eye(2).reshape ((1,)*(b.ndim-2) + (2,2))
b33 = self.C33_planestress(II_2D)
J = np.sqrt(b33*II_2D)
else :
raise NotImplementedError(f"Cauchy not implemented for \
{str(b.shape[-2:])} and {self.hyp.name}" )
return mu/J*(b-Id) + lamb*np.log(J)/J*Id
[docs]
def spacial_elasticty_tensor(self, b:FloatArray, out_c: FloatArray|None = None)->FloatArray:
"""
Compute the spatial (Eulerian) elasticity tensor for a Compressible Neo-Hookean material.
This tensor (c) relates the Lie derivative of the Kirchhoff stress to the
rate of deformation tensor. It is derived from the push-forward of the
material elasticity tensor to the current configuration.
The formula used is:
c = (lamb / J) * (I ⊗ I) + 2 * (mu - lamb * ln(J)) / J * Is
where J = det(F) = sqrt(det(b)), I ⊗ I is the identity for the volumetric
part, and Is is the symmetric fourth-order identity tensor.
Args:
b (ndarray): Left Cauchy-Green deformation tensor (Finger tensor).
Shape: (dim, dim) or (batch_dims, dim, dim).
Supported dims: 3D or 2D (under Plane Strain hypothesis).
Returns:
ndarray: Fourth-order spatial elasticity tensor.
Shape: (dim, dim, dim, dim) or (batch_dims, dim, dim, dim, dim).
Raises:
NotImplementedError: If the dimensions are not (3,3) or if (2,2) is
provided without Plane Strain hypothesis.
"""
lamb, mu = (self.lamb, self.mu)
if (b.shape[-2:] == (2,2) and self.hyp == MaterialHypothesis.PLANESTRAIN) :
if out_c is None:
batch_shape = b.shape[:-2]
out_c = np.empty(batch_shape + (2,2,2,2), dtype=b.dtype)
_NH_spacial_elasticity_plane_strain_guvec(b, self.lamb, self.mu, \
out=out_c) # type: ignore
return out_c
if b.shape[-2:] == (3,3):
nbtens = b.shape[:-2]
bcshape = (1,)*len(nbtens)
J = np.sqrt(det(b)).reshape((-1,1,1,1,1))
IxI = Sp4[3].reshape(bcshape+(3,3,3,3))
Is = Sy4[3].reshape(bcshape+(3,3,3,3))
c = lamb/J*IxI + 2*(mu-lamb*np.log(J))/J*Is
if b.ndim == 2 :
return c.squeeze()
return c
else :
raise NotImplementedError()
[docs]
def get_properties(self):
""" return material properies mu and lambda """
return np.array([self.mu, self.lamb], dtype = npfloat)
[docs]
@staticmethod
@njit(inline='always', fastmath=True, cache=True)
def sigma_c_kernel(F, props, out_sigma, out_c):
"""
Calcule la contrainte de Cauchy et le tenseur d'élasticité spatial 2D.
Zéro allocation, remplissage en place.
"""
mu = props[0]
lamb = props[1]
# --- Cinématique ---
# J = det(F)
J = F[0, 0] * F[1, 1] - F[0, 1] * F[1, 0]
# b = F * F^T (Tenseur de Cauchy-Green gauche)
b00 = F[0, 0]*F[0, 0] + F[0, 1]*F[0, 1]
b01 = F[0, 0]*F[1, 0] + F[0, 1]*F[1, 1]
b11 = F[1, 0]*F[1, 0] + F[1, 1]*F[1, 1]
fill_nh_sigma_2d(J, b00, b11, b01, mu, lamb, out_sigma)
fill_nh_c_2d(J, mu, lamb, out_c)
[docs]
def check_material_law_F(law:HyperElasticMaterial, F:FloatArray, verb = True)->Tuple[bool, bool]:
''' check PKI(F) and dPKIdF(F) return PKI_ok and dPKIDF_ok
PKI_ok and dPKIdF_ok, respectively True if PKI (dPKIdF) is close to
the finite difference approximation, False otherwise.
'''
# pylint: disable=C0103
#law = law
PKI = law.PKI(F)
PKI_num = df_num(law.potentialF, F, F.shape[-2:], ())
PKI_ok = check_relerr(PKI, PKI_num, PKI.shape[-2:] )
dPKIdF = law.dPKIdF(F)
dPKIdF_num = df_num(law.PKI, F, F.shape[-2:], F.shape[-2:])
dPKIdF_ok = check_relerr(dPKIdF, dPKIdF_num, PKI.shape[-2:] + F.shape[-2:], verbose= verb)
return PKI_ok, dPKIdF_ok
[docs]
def check_material_law_C(law:HyperElasticMaterial, C:FloatArray, verb = True)->Tuple[bool, bool]:
''' check PKII(C) and dPKIIdC(C) return PKII_ok and dPKIIDC_ok
PKII_ok and dPKIIdC_ok, respectively True if PKII (dPKIIdC) is close to
the finite difference approximation, False otherwise.
'''
# pylint: disable=C0103
#verb = True
PKII = law.PKII(C)
shape = C.shape[-2:]
PKII_num = 2.*df_num(law.potentialC, C, shape, ())
PKII_ok = check_relerr(PKII, PKII_num, PKII.shape[-2:] , verbose= verb)
dPKIIdC = law.dPKIIdC(C)
dPKIIdC_num = df_num(law.PKII, C, shape, shape)
#dPKIIdC_num2 = df_num(lambda C : law.PKII( C.swapaxes(-1,-2)), C, shape, shape)
#dPKIIdC_num = 0.5*(dPKIIdC_num + dPKIIdC_num2)
dPKIIdC_ok = check_relerr(dPKIIdC, dPKIIdC_num, PKII.shape[-2:] + C.shape[-2:], verbose = verb)
return PKII_ok, dPKIIdC_ok
[docs]
def benchmark(law:HyperElasticMaterial, n:int, mem = False, prof = True):
""" benchmark : compute and profile law.PKI and law.dPKIDF """
# pylint: disable=C0103
dim = HYP2DIM[law.hyp]
if law.hyp in [PLANESTRESS, PLANESTRAIN] and mem:
F = random_lintransform2dto3d(n)
print(f'computing PKI and dPKIDF for {n:.3e} F32, in in {law.hyp}')
else :
gen = random_lintransform
F = gen(n, dim)
print(f'computing PKI and dPKIDF for {n:.3e} F, in {law.hyp}')
Profile.doprofile = prof
Profile.enable()
start = time.time()
__ = law.PKI(F)
__ = law.dPKIdF(F)
end = time.time()
Profile.disable()
Profile.print_stats()
print('time : ', end - start)
if __name__ == "__main__":
lawb = CompressibleNeoHookeen(E = 1., nu = 0.3, hyp = PLANESTRESS)
benchmark(lawb, 1000000, prof=True)
# computing PKI and dPKIDF for 1.000e+06 F, in MaterialHypothesis.DIM3
# 550 function calls (532 primitive calls) in 1.540 seconds
# lawPSTRAIN = CompressibleNeoHookeen(E=1.2, nu = 0.3, hyp=PLANESTRAIN)
# lawPSTRESS = CompressibleNeoHookeen(E=1.2, nu = 0.3, hyp=PLANESTRESS)
# # benchmark(lawPSTRESS, 10000, mem =True, prof =True)
# law1 = SimpleLaw0()
# law2 = SimpleLaw1()
# law3 = SimpleLaw2()
# F2d0 = np.array([[2.,0.1],[1.,1.3]])
# F2d1 = np.array([[[2.,0.1],[1.,1.3]]])
# F2d2 = np.array([[[2.,0.1],[1.,1.3]],
# [[2.,0.5],[3.,1.1]]
# ])
# Fm = np.array( [
# [[1.,0.],[0.,1.],[0,0]],
# [[1.,0.],[0.,0.],[0,1]]
# ])
# #Fm = np.array([[1.,0.],[0.,1.],[0,0]])
# Fm2d = np.array([[2.,0.1],[1.,1.3],[0.,0.]])
# Fm = np.array([[[2.,0.], [0.,1.], [0,0]],
# [[2.,0.1],[0.1,1.],[0,0]]])
# Fm0= np.array([[[3.,0.1], [1., 0.1], [1.,1.]],
# [[2.,0.1],[0.1,1.],[0.1,0.2]]])
# Fm1= np.array([[[3.,0.1], [1., 0.1], [1.,1.]],
# [[2.,0.1],[0.1,1.],[0.1,0.2]]])
# lawu = lawPSTRESS
# Fu = Fm1
# print(test_material_law(lawu, Fu))
# computing PKI and dPKIDF for 1.000e+06 F, in MaterialHypothesis.DIM3
# 62 function calls in 1.243 seconds
# Ordered by: cumulative time
# ncalls tottime percall cumtime percall filename:lineno(function)
# 1 0.013 0.013 1.031 1.031 hyperelasticity.py:221(dPKIdF)
# 1 0.682 0.682 0.865 0.865 hyperelasticity.py:186(_dPKIdF_exp)
# 2 0.035 0.017 0.310 0.155 smallsize_tensor.py:69(det_invt)
# 2 0.262 0.131 0.262 0.131 smallsize_tensor.py:25(_cofac)
# 1 0.000 0.000 0.211 0.211 hyperelasticity.py:214(PKI)
# 1 0.003 0.003 0.211 0.211 hyperelasticity.py:164(_PKI_3d)
# 1 0.182 0.182 0.182 0.182 smallsize_tensor.py:79(cross_ijkl)
# 1 0.052 0.052 0.052 0.052 hyperelasticity.py:171(_PKI_exp)
# 2 0.014 0.007 0.014 0.007 smallsize_tensor.py:41(_det)