Coverage for pyxfem/solidmechanics/hyperelasticity.py: 55%
480 statements
« prev ^ index » next coverage.py v7.17.0, created at 2026-10-24 07:17 +0000
1#!/usr/bin/env python3
2# -*- coding: utf-8 -*-
3"""
4This module is dedicated to hyperelasticity
5Created on Thu April 2024 14:09:29 2024
6@author: nchevaug
7"""
8from __future__ import annotations
9import re
10import time
11from typing import Callable, Tuple, cast
12import numpy as np
13from numba import guvectorize, float32, float64, njit, prange
14from scipy.special import lambertw as sp_lambertw
15from pyxfem.core import smallsize_tensor
16from pyxfem.core.config import FloatArray, FloatVars, npfloat
17from pyxfem.core.smallsize_tensor import det_invt, cross_ijkl, df_num,\
18 random_lintransform, random_lintransform2dto3d,\
19 metric, det, trace, trace2, inv, invt, dinvt, check_relerr, \
20 Sy4, Sp4
21from pyxfem.core.exec_tools import Profile
22from pyxfem.solidmechanics.elasticity import MaterialHypothesis
23#import pytest
25#from sym_hyperelasticity import II
27PLANESTRAIN = MaterialHypothesis.PLANESTRAIN
28PLANESTRESS = MaterialHypothesis.PLANESTRESS
29DIM3 = MaterialHypothesis.DIM3
30HYP2DIM = {PLANESTRESS:2, PLANESTRAIN:2, DIM3:3}
32import numpy as np
34def real_lambertw(z):
35 """
36 Computes the principal branch of the Lambert W function for real, positive z
37 using vectorized Halley's method.
38 """
39 # A robust initial guess for z >= 0
40 w = np.log(1.0 + z)
42 # Halley's method usually converges to machine precision in 3-4 iterations
43 for _ in range(5):
44 ew = np.exp(w)
45 w1 = w + 1.0
46 f = w * ew - z
48 # Halley step
49 step = f / (ew * w1 - ((w + 2.0) * f) / (2.0 * w1))
50 w -= step
52 # Optional: break early if max error is small enough,
53 # but for 5 fixed iterations, it's often faster to just loop.
54 if np.max(np.abs(step)) < 1e-12:
55 break
57 return w
59#lambertw = lambda x : np.real(sp_labertw(x))
60lambertw = real_lambertw
62def vecroots(f:Callable, fprime:Callable, x0:FloatVars,
63 atol:float=1.e-20, rtol=1.e-12, maxiter:int = 100, verbose = False)->FloatVars:
64 """ compute roots of an array of independant equations using newton's
65 method with line search so that all x stay > 0"""
66 r = f(x0)
67 if verbose:
68 print(0, x0, np.abs(r), 1.)
69 err0 = np.max(np.abs(r))
70 err = err0
71 if err0 < atol:
72 return x0
73 it = 0
74 while it < maxiter:
75 alpha=1.
76 j = fprime(x0)
77 s = -r/j
78 x = x0 + s
79 while np.any(x<=0.):
80 alpha = 0.5*alpha
81 x = x0 + alpha*s
82 r = f(x)
83 while np.max(np.abs(r)) > err:
84 alpha = 0.5*alpha
85 x = x0 + alpha*s
86 r = f(x)
87 err = np.max(np.abs(r))
88 if err < atol or err < rtol*err0 :
89 return x
90 x0 = x
91 it += 1
92 if verbose:
93 print(it, x0, np.abs(r), alpha)
94 msg = f"Failed to converge after {maxiter:d} iterations"
95 raise RuntimeError(msg)
97class HyperElasticMaterial():
98 """ Virtual base class for HyperElastic material """
99 # pylint: disable=C0103
100 # pylint: disable=W2301
101 def __init__(self, hyp:MaterialHypothesis=DIM3, rho0 = 1.):
102 self.hyp, self.rho = hyp, rho0
103 def potentialF(self, F:FloatArray )-> FloatVars:
104 """ Compute the Elastic Potential as a Function of F"""
105 raise NotImplementedError()
106 def potentialC(self, C:FloatArray )-> FloatVars:
107 """ Compute the Elastic Potential as a Function of C"""
108 raise NotImplementedError()
109 def PKI(self, F:FloatArray)->FloatArray:
110 """ Compute Piola Kirchhoff I Stress tensor as a function of F """
111 raise NotImplementedError()
112 def PKII(self, C:FloatArray)->FloatArray:
113 """ Compute Piola Kirchhoff I Stress tensor as a function of C = F^TF """
114 raise NotImplementedError()
115 def Cauchy(self, b:FloatArray)->FloatArray:
116 """ Compute Cauchy Stress tensor as a function of b = FF^T """
117 raise NotImplementedError()
118 def dPKIdF(self, F:FloatArray)->FloatArray:
119 """ Compute the deriviatives of the Piola Kirchhoff I Stress tensor
120 with regard to F
121 """
122 raise NotImplementedError()
123 def dPKIIdC(self, C:FloatArray)->FloatArray:
124 """ Compute the deriviatives of the Piola Kirchhoff II Stress tensor
125 with regard to C
126 """
127 raise NotImplementedError()
129 def dPKIIdE(self, C:FloatArray)->FloatArray:
130 """ Compute the deriviatives of the Piola Kirchhoff II Stress tensor
131 with regard to C
132 """
133 return 2*self.dPKIIdC(C)
135 def spacial_elasticty_tensor(self, b:FloatArray)->FloatArray:
136 """ Compute the spacial elasticity tensor (1/J FFFF dPKIIdE)
137 """
138 raise NotImplementedError()
141# class SimpleLaw0(HyperElasticMaterial):
142# """ A very simple hyper elastic law. Mostly to check derivatives of I_C """
143# # pylint: disable=C0103
144# def __init__(self):
145# super().__init__(PLANESTRAIN)
146# def potential(self, F:FloatArray)-> FloatArray|float:
147# m = metric(F)
148# I_C = 1.+trace(m)
149# return 0.5*I_C
150# def PKI(self, F:FloatArray)->FloatArray:
151# return F
152# def dPKIdF(self, F:FloatArray)->FloatArray:
153# tshape = F.shape[-2:]
154# nbfirstaxis = F.ndim-2
155# I40 = np.eye(tshape[0]*tshape[1]).reshape(tshape+tshape)
156# return I40[(np.newaxis,)*nbfirstaxis]
158# class SimpleLaw1(HyperElasticMaterial):
159# """ A very simple hyper elastic law. Mostly to check derivatives of log(J) """
160# # pylint: disable=C0103
161# def __init__(self):
162# super().__init__(PLANESTRAIN)
163# def potential(self, F:FloatArray)-> FloatArray|float:
164# return 0.5*np.log(det(metric(F)))
166# def PKI(self, F:FloatArray)->FloatArray:
167# return invt(F)
169# def dPKIdF(self, F:FloatArray)->FloatArray:
170# return dinvt(F)
172# class SimpleLaw2(HyperElasticMaterial):
173# """ A very simple hyper elastic law. Mostly to check derivatives of J """
174# # pylint: disable=C0103
175# def __init__(self):
176# super().__init__(PLANESTRAIN)
177# def potential(self, F:FloatArray)-> FloatArray|float:
178# return np.sqrt(det(metric(F)))
180# def PKI(self, F:FloatArray)->FloatArray:
181# J, invFT = det_invt(F)
182# return np.atleast_1d(J)[..., np.newaxis, np.newaxis] * invFT
184# def dPKIdF(self, F:FloatArray)->FloatArray:
185# J, invFT = det_invt(F)
186# J_11 = np.atleast_1d(J)[..., np.newaxis, np.newaxis]
187# PKI = J_11 * invFT
188# J_1111 = J_11[..., np.newaxis, np.newaxis]
189# return cross_ijkl(PKI, invFT) + J_1111*dinvt(F)
191@njit(inline='always', fastmath=True)
192def fill_nh_sigma_2d(J, b00, b11, b01, mu, lamb, out_sigma):
193 """ Remplit la contrainte de Cauchy 2D en place """
194 S_mu = mu / J
195 S_lamb = (lamb * np.log(J)) / J
197 out_sigma[0, 0] = S_mu * (b00 - 1.0) + S_lamb
198 out_sigma[1, 1] = S_mu * (b11 - 1.0) + S_lamb
199 out_sigma[0, 1] = S_mu * b01
200 out_sigma[1, 0] = out_sigma[0, 1]
202@njit(inline='always', fastmath=True)
203def fill_nh_c_2d(J, mu, lamb, out_c):
204 """ Remplit le tenseur d'élasticité spatial 2x2x2x2 en place """
205 S1 = lamb / J
206 S2 = 2.0 * (mu - lamb * np.log(J)) / J
207 S2_half = 0.5 * S2
209 out_c.fill(0.0)
210 out_c[0, 0, 0, 0] = S1 + S2
211 out_c[1, 1, 1, 1] = S1 + S2
212 out_c[0, 0, 1, 1] = S1
213 out_c[1, 1, 0, 0] = S1
214 out_c[0, 1, 0, 1] = S2_half
215 out_c[0, 1, 1, 0] = S2_half
216 out_c[1, 0, 0, 1] = S2_half
217 out_c[1, 0, 1, 0] = S2_half
220@guvectorize([
221 (float64[:,:], float64, float64, float64[:,:,:,:]),
222 (float32[:,:], float32, float32, float32[:,:,:,:])
223], '(m,m),(),()->(m,m,m,m)', nopython=True, cache=True)
224def _NH_spacial_elasticity_plane_strain_guvec(b, lamb, mu, out):
225 det_b = b[0, 0] * b[1, 1] - b[0, 1] * b[1, 0]
226 J = np.sqrt(det_b)
227 fill_nh_c_2d(J, mu, lamb, out)
229@guvectorize([
230 ('float64[:,:], float64, float64, float64[:,:]'),
231 ('float32[:,:], float32, float32, float32[:,:]')
232], '(m,m),(),()->(m,m)', nopython=True, cache=True)
233def _NH_cauchy_plane_strain_guvec(b, lamb, mu, out_sigma):
234 """
235 Calcule la contrainte de Cauchy 2D (Déformation Plane) pour un
236 matériau Néo-Hookéen Compressible.
238 Inputs:
239 b : Tenseur de Cauchy-Green gauche (2x2)
240 lamb : Paramètre de Lamé lambda (scalaire)
241 mu : Module de cisaillement mu (scalaire)
242 Output:
243 out_sigma : Tenseur des contraintes de Cauchy (2x2)
244 """
245 # 1. Calcul du déterminant et du Jacobien (J = sqrt(det(b)))
246 det_b = b[0, 0] * b[1, 1] - b[0, 1] * b[1, 0]
247 J = np.sqrt(det_b)
248 fill_nh_sigma_2d(J, b[0,0], b[1,1], b[0,1], mu, lamb, out_sigma)
251@njit(inline='always', fastmath=True, cache=True)
252def neo_hookean_law_local(F, props, out_sigma, out_c):
253 """
254 Calcule la contrainte de Cauchy et le tenseur d'élasticité spatial 2D.
255 Zéro allocation, remplissage en place.
256 """
257 mu = props[0]
258 lamb = props[1]
259 # --- Cinématique ---
260 # J = det(F)
261 detF = F[0, 0] * F[1, 1] - F[0, 1] * F[1, 0]
262 J = detF
263 # b = F * F^T (Tenseur de Cauchy-Green gauche)
264 b00 = F[0, 0]*F[0, 0] + F[0, 1]*F[0, 1]
265 b01 = F[0, 0]*F[1, 0] + F[0, 1]*F[1, 1]
266 b11 = F[1, 0]*F[1, 0] + F[1, 1]*F[1, 1]
267 fill_nh_sigma_2d(J, b00, b11, b01, mu, lamb, out_sigma)
268 fill_nh_c_2d(J, mu, lamb, out_c)
270class CompressibleNeoHookeen(HyperElasticMaterial):
271 """ Implementation of the NeoHookeen Hyperelastic material """
272 # pylint: disable=C0103
273 def __init__(self,*, # pylint: disable=R0913
274 E:float|None=None,
275 nu:float|None=None,
276 lamb:float|None=None,
277 mu:float|None=None,
278 rho0:float|None = None,
279 hyp:MaterialHypothesis=DIM3):
280 super().__init__(hyp)
281 if isinstance(E, float) and isinstance(nu, float):
282 self.E = E
283 self.nu = nu
284 self.lamb = E*nu/(1.+nu)/(1.-2*nu)
285 self.mu = E/2./(1+nu)
286 elif isinstance(lamb, float) and isinstance(mu, float):
287 self.lamb = lamb
288 self.mu = mu
289 self.E = mu*(3.0*lamb + 2.0*mu)/(lamb + mu)
290 self.nu = 0.5*lamb/(lamb + mu)
291 else:
292 raise ValueError("E and nu or lambda and mu must be set (exclusive)")
293 self.rho0 = rho0
295 def _checktensorsize(self, F:FloatArray):
296 if F.shape[-2:] == (3,3) and self.hyp == DIM3:
297 return
298 if F.shape[-2:] in [(2,2), (3,2)]:
299 if self.hyp in [PLANESTRAIN, PLANESTRESS]:
300 return
301 raise ValueError("invalid shape for F. F.shape must be (...2,2) or (...,3,3)or (..., 3,2)")
303 def l3_planestress(self, J2d:FloatVars)->FloatVars:
304 """ return the principal elongation lambda_3 in the direction orthogonal to the plane """
305 lamb, mu = (self.lamb, self.mu)
306 a = mu
307 b = lamb/2
308 c = lamb*np.log(J2d) -mu
309 #return np.sqrt(b*np.real(lambertw( a* np.exp(-c/b)/b)) /a)
310 return np.sqrt(b*lambertw( a* np.exp(-c/b)/b) /a)
312 # def resl3(l3):
313 # return lamb*np.log(J2d*l3) + mu*(l3*l3 - 1)
314 # def jacl3(l3):
315 # return lamb/l3 + 2*mu*l3
316 # if isinstance(J2d, float) or isinstance(J2d, int):
317 # x0 = 1.
318 # else:
319 # x0 = np.ones(J2d.shape)
320 # l3 = vecroots(resl3, jacl3, x0, atol=1.e-20, rtol=1.e-12, maxiter=100)
321 #return l3
323 # def _dl3dF2d_planestress(self, l3, invFT2d):
324 # lamb = self.lamb
325 # mu = self.mu
326 # return -l3/(1+2*mu*l3**2/lamb)*invFT2d
328 # def dl3F2d_planestress(self, F2d):
329 # """ Compute the derivative of the third principal strain as a fonction of
330 # the plane part of F (Assuming Plan Stress)"""
331 # J2d, invFT2d = det_invt(F2d)
332 # l3 = self.l3_planestress(J2d)
333 # return self._dl3dF2d_planestress(l3, invFT2d)
335 def _potential(self, J:FloatVars, I_C:FloatVars)-> FloatVars:
336 logJ = np.log(J)
337 return self.mu*(0.5*(I_C-3)- logJ) + self.lamb*logJ**2/2
339 def _potential_3d(self, F:FloatArray)-> FloatVars:
340 J = det(F)
341 I_C = trace2(F)
342 return self._potential(J, I_C)
344 def _potential_planestrain(self, F2d:FloatArray)-> FloatVars:
345 if F2d.shape[-2:] == (3,2):
346 m = metric(F2d)
347 J = np.sqrt(det(m))
348 I_C = 1.+trace(m)
349 elif F2d.shape[-2:] == (2,2):
350 J = det(F2d)
351 I_C = trace2(F2d) +1
352 else :
353 raise ValueError('_potential_planestrain, F2d must be of shape (...,3,2) or (...,2,2)')
354 return self._potential(J, I_C)
356 def _potential_planestress(self, F2d:FloatArray)-> FloatVars:
357 if F2d.shape[-2:] == (3,2):
358 m = metric(F2d)
359 J2d = np.sqrt(det(m))
360 I_C2d = trace(m)
361 elif F2d.shape[-2:] == (2,2):
362 J2d = det(F2d)
363 I_C2d = trace2(F2d)
364 else :
365 raise ValueError('_potential_planestress, F2d must be of shape (...,3,2) or (...,2,2)')
366 l3 = self.l3_planestress(J2d)
367 J = l3*J2d
368 I_C = l3**2+I_C2d
369 return self._potential(J, I_C)
371 def _PKI_3d(self, F:FloatArray)->FloatArray:
372 '''Compute Piola Kirchhoff I Stress tensor as a function of F: 3D case (and Plane Strain)'''
373 if F.shape[-2:] not in [(3,3), (3,2), (2,2)]:
374 raise ValueError("F must be a nd.array of shape[-2:] in [(3,3), (3,2), (2,2)")
375 J, invFT = det_invt(F)
376 #print("### ", J, F, invFT)
377 return self._PKI_exp(J, F, invFT)
379 def _PKI_exp_v0(self, J, F, invFT)->FloatArray:
380 ''' Compute Piola Kirchhoff I Stress tensor as a function of J, F, invFT'''
381 return self.mu*(F-invFT) + self.lamb*np.log(J)[..., np.newaxis, np.newaxis]*invFT
383 def _PKI_exp(self, J, F, invFT) -> FloatArray:
384 ''' Compute Piola Kirchhoff I Stress tensor as a function of J, F, invFT'''
385 # 1. Factor out the scalar coefficients to eliminate array subtractions
386 coef = (self.lamb * np.log(J) - self.mu)[..., np.newaxis, np.newaxis]
387 # 2. Allocate exact output memory once
388 ret = np.empty_like(F)
389 # 3. Compute in-place to bypass temporary array copies
390 np.multiply(F, self.mu, out=ret)
391 ret += coef * invFT
392 return ret
394 def _PKI_planestress(self, F2d:FloatArray)->FloatArray:
395 ''' Compute Piola-Kirchhoff I Stress tensor as a function of F : plane stress case)'''
396 if F2d.shape[-2:] not in[(2,2), (3,2)]:
397 raise ValueError("F must be a nd.array of shape (...,3,3)")
398 J2d, invFT2d = det_invt(F2d)
399 if np.any(J2d<=0.):
400 raise ValueError("det(F2d) must be > 0.")
401 l3 = self.l3_planestress(J2d)
402 J = l3*J2d
403 return self._PKI_exp(J, F2d, invFT2d)
407 def _dPKIdF_exp_v0(self, c0, c1, c2, invFT)->FloatArray:
408 c0 = np.atleast_1d(c0)[(...,) + (np.newaxis,)*4]
409 c1 = np.atleast_1d(c1)[(...,) + (np.newaxis,)*4]
410 c2 = np.atleast_1d(c2)[(...,) + (np.newaxis,)*4]
411 tshape = invFT.shape[-2:]
412 I4_ = np.eye(tshape[0]*tshape[1]).reshape(tshape+tshape)
413 I4_ = I4_.reshape((1,)*(invFT.ndim-2)+I4_.shape, copy=False)
414 c0_I4 = c0*I4_
415 invFTxinvFT = cross_ijkl(invFT, invFT) # iFT x d(log(J)) = iFTxiFT
416 c2_invFTxinvFT = c2*invFTxinvFT
417 c1_dinvFTdF = -c1*(invFTxinvFT.transpose(list(range(invFTxinvFT.ndim-4)) + [-4,-1,-2,-3]))
418 #c1_dinvFTdF = c1*dinvt(invFT, invFT)
419 ret = c0_I4 + c1_dinvFTdF +c2_invFTxinvFT
420 if invFT.shape[:-2] == ():
421 return ret.squeeze()
422 return ret
424 def _dPKIdF_exp_v1(self, c0, c1, c2, invFT) -> FloatArray:
425 # 1. Expand scalars/arrays for 4D tensor broadcasting
426 c0 = np.atleast_1d(c0)[..., None, None, None, None]
427 c1 = np.atleast_1d(c1)[..., None, None, None, None]
428 c2 = np.atleast_1d(c2)[..., None, None, None, None]
430 # 2. Build the 4th-order identity tensor I4
431 d = invFT.shape[-1]
432 I4 = np.eye(d * d).reshape((d, d, d, d))
434 # 3. Allocate full memory block to match the number of quadrature points
435 full_shape = invFT.shape[:-2] + (d, d, d, d)
436 ret = np.empty(full_shape, dtype=invFT.dtype)
437 ret[...] = c0 * I4
439 # 4. Compute and subtract the c1 * dinvFTdF term in-place
440 ret -= c1 * (invFT[..., :, None, None, :] * invFT.swapaxes(-1, -2)[..., None, :, :, None])
442 # 5. Compute and add the c2 * invFTxinvFT term in-place
443 ret += c2 * (invFT[..., :, :, None, None] * invFT[..., None, None, :, :])
445 if invFT.ndim == 2:
446 return ret.squeeze()
447 return ret
449 def _dPKIdF_exp_v2(self, F)->FloatArray:
450 mu, lamb = self.mu, self.lamb
451 tshape = F.shape[-2:]
452 I4_ = np.eye(tshape[0]*tshape[1]).reshape(tshape+tshape)
453 I4_.shape = (1,)*(F.ndim-2)+I4_.shape
454 J, iFT = det_invt(F)
455 J_1111 = np.atleast_1d(J)[(...,) + (np.newaxis,)*4]
456 return mu*I4_ + lamb*cross_ijkl(iFT, iFT) + (lamb*np.log(J_1111)-mu)*dinvt(F)
458 def _dPKIdF_exp_v3(self, c0, c1, c2, invFT) -> FloatArray:
459 d = invFT.shape[-1]
460 I4 = np.eye(d * d).reshape((d, d, d, d))
462 # Expand constants for broadcasting
463 c0_b = np.atleast_1d(c0)[..., None, None, None, None]
464 c1_b = np.atleast_1d(c1)[..., None, None, None, None]
465 c2_b = np.atleast_1d(c2)[..., None, None, None, None]
467 # 1. Base allocation (c0 * I4)
468 full_shape = invFT.shape[:-2] + (d, d, d, d)
469 ret = np.empty(full_shape, dtype=invFT.dtype)
470 ret[...] = c0_b * I4
472 # 2. Term: -c1 * invFT_il * invFT_jk
473 # Directly computes the transposed term in C without copying memory
474 ret -= c1_b * np.einsum('...il,...kj->...ijkl', invFT, invFT)
476 # 3. Term: c2 * invFT_ij * invFT_kl
477 # Replaces cross_ijkl with direct C-loop
478 ret += c2_b * np.einsum('...ij,...kl->...ijkl', invFT, invFT)
480 if invFT.ndim == 2:
481 return ret.squeeze()
482 return ret
484 @staticmethod
485 @njit(parallel=True, fastmath=True, cache = True)
486 def _dPKIdF_numba_core(c0, c1, c2, invFT, ret):
487 N = invFT.shape[0]
488 d = invFT.shape[1]
490 for g in prange(N):
491 c0_g = c0[g]
492 c1_g = c1[g]
493 c2_g = c2[g]
495 for i in range(d):
496 for j in range(d):
497 for k in range(d):
498 for l in range(d):
499 val = 0.0
500 if i == k and j == l:
501 val += c0_g
502 val -= c1_g * invFT[g, i, l] * invFT[g, k, j]
503 val += c2_g * invFT[g, i, j] * invFT[g, k, l]
505 ret[g, i, j, k, l] = val
507 def _dPKIdF_exp(self, c0, c1, c2, invFT):
508 orig_shape = invFT.shape
509 d = orig_shape[-1]
511 if invFT.ndim > 2:
512 N = np.prod(orig_shape[:-2], dtype=int)
513 else:
514 N = 1
516 invFT_flat = invFT.reshape(N, d, d)
518 c0_flat = np.broadcast_to(c0, (N,)).astype(invFT.dtype)
519 c1_flat = np.broadcast_to(c1, (N,)).astype(invFT.dtype)
520 c2_flat = np.broadcast_to(c2, (N,)).astype(invFT.dtype)
522 ret_flat = np.empty((N, d, d, d, d), dtype=invFT.dtype)
524 # Call the static method via self
525 self._dPKIdF_numba_core(c0_flat, c1_flat, c2_flat, invFT_flat, ret_flat)
527 if invFT.ndim == 2:
528 return ret_flat.squeeze()
530 return ret_flat.reshape(orig_shape[:-2] + (d, d, d, d))
532 def potentialF(self, F:FloatArray)-> FloatVars:
533 """ return the Hyperelastic potential as a function of F"""
534 self._checktensorsize(F)
535 if self.hyp == DIM3:
536 return self._potential_3d(F)
537 if self.hyp == PLANESTRAIN:
538 return self._potential_planestrain(F)
539 if self.hyp == PLANESTRESS:
540 return self._potential_planestress(F)
541 raise ValueError("MaterialHypothesis" + str(self.hyp) +"Not implemented")
543 def PKI(self, F:FloatArray)->FloatArray:
544 ''' Compute Piola Kirchhoff I Stress tensor as a function of F'''
545 if self.hyp in [PLANESTRAIN, DIM3]:
546 return self._PKI_3d(F)
547 return self._PKI_planestress(F)
549 def dPKIdF(self, F:FloatArray)->FloatArray:
550 ''' Compute the deriviatives of the Piola Kirchhoff I Stress tensor
551 with regard to F at
552 '''
553 self._checktensorsize(F)
554 J, invFT = det_invt(F)
555 #print('invFT shap', invFT.shape)
556 mu =self.mu
557 lamb = self.lamb
558 c0 = mu
559 if self.hyp in [DIM3, PLANESTRAIN]:
560 c1 = lamb*np.log(J)-mu
561 c2 = lamb
562 elif self.hyp == PLANESTRESS :
563 J2d = J
564 l3 = self.l3_planestress(J2d)
565 J = l3*J2d
566 c1 = lamb*np.log(J)-mu
567 c2 = lamb*(1.-lamb/(lamb+2*mu*l3*l3))
568 else:
569 raise ValueError(f"dPKIdF not implemented for {str(self.hyp):}")
570 if F.shape[-2:] == (3,2):
571 I4_ = np.eye(6).reshape((1,)*(F.ndim-2)+(3,2,3,2), copy = False)
572 c1 = np.atleast_1d(c1)[(...,)+ (np.newaxis,)*4]
573 c2 = np.atleast_1d(c2)[(...,)+ (np.newaxis,)*4]
574 return c0*I4_ + c2*cross_ijkl(invFT, invFT) + c1*dinvt(F)
575 return self._dPKIdF_exp(c0, c1, c2, invFT)
577 def C33_planestress(self, II_C2D:FloatVars) -> FloatVars:
578 """ Return the square of the thichkess dilatation in plane stress Case,
579 knowing the secons invariant of the planar componant of C
580 mu(C33 -1) + kappa ln(J) = sig33
581 """
582 mu, lamb= (self.mu, self.lamb)
583 a = 0.5*lamb/mu
584 ia = 1/a
585 #b = np.log(II_C2D)
586 # r0 = a*np.real(lambertw( np.exp(ia-b)*ia ))
587 # r1 = a*np.real(lambertw( ia/II_C2D*np.exp(ia))
588 # print("ga", r0-r1)
589 return a*np.real(lambertw( ia/II_C2D*np.exp(ia)))
591 def dC33dII_planestress(self, II_C2D:FloatVars, C33:FloatVars|None = None) -> FloatVars:
592 """ Return the square of the thichkess dilatation in plane stress Case,
593 knowing the secons invariant of the planar componant of C """
594 mu, lamb= (self.mu, self.lamb)
596 II = II_C2D
597 if C33 is None :
598 C33 = self.C33_planestress(II)
599 a = cast(npfloat, 0.5*lamb/mu)
600 dC33dII = -a*C33/(II*(a+C33))
601 return dC33dII
603 def potentialC(self, C:FloatArray) -> FloatArray:
604 """ Compute The elastic potential as a function of C"""
605 mu, lamb= (self.mu, self.lamb)
606 if C.shape[-2:] == (3,3):
607 I = trace(C)
608 logIII = np.log(np.linalg.det(C))
609 elif C.shape[-2:] == (2,2):
610 I_C2D = trace(C)
611 II_C2D = det(C)
612 if self.hyp == MaterialHypothesis.PLANESTRESS :
613 C33 = self.C33_planestress(II_C2D)
614 elif self.hyp == MaterialHypothesis.PLANESTRAIN:
615 C33 = 1.
616 else:
617 raise ValueError("self.hyp must be either PLANESTRESS or PLANESTRAIN")
618 I = I_C2D + C33
619 logIII = np.log(II_C2D*C33)
620 else :
621 raise ValueError("C must be of shape (..., 2,2) or (..., 3, 3)")
622 return mu*(I - 3 - logIII)/2. + lamb*logIII**2/8
624 def PKII(self, C:FloatArray) -> FloatArray:
625 """ Return second Piola Kirchoff Stres as a function of C """
626 lamb, mu = (self.lamb, self.mu)
627 if C.shape[-2:] == (3,3) or \
628 (C.shape[-2:] == (2,2) and self.hyp == MaterialHypothesis.PLANESTRAIN) :
629 dim = C.shape[-1]
630 III = np.linalg.det(C).reshape( C.shape[:-2] + (1,1))
631 logIII = np.log(III)
632 Id = np.eye(dim).reshape( (1,)*(C.ndim-2) + (dim, dim) )
633 iCT = invt(C, III)
635 kappa = 0.5*lamb*logIII - mu
636 # print("fff", mu*Id[0,2,2], kappa[0], iCT[0,2,2])
637 elif C.shape[-2:] == (2,2) and self.hyp == MaterialHypothesis.PLANESTRESS :
638 Id = np.eye(2).reshape ((1,)*(C.ndim-2) + (2,2))
639 II_2D = det(C).reshape(C.shape[:-2]+(1,1))
640 iCT = invt(C, II_2D)
641 C33 = self.C33_planestress(II_2D)
642 logIII = np.log(C33*II_2D)
643 kappa = mu*C33*(lamb*(logIII-1) - 2*mu )/(2*mu*C33+lamb)
644 else :
645 raise NotImplementedError()
648 return mu*Id + kappa*iCT
650 def dPKIIdC(self, C:FloatArray) -> FloatArray:
651 lamb, mu = (self.lamb, self.mu)
652 iC = inv(C)
653 #iCxiC = cross_ijkl(iC, iC)
654 #diCdC = dinv(C, iA= iC)
655 diCTdC = dinvt(C)
656 iCT = smallsize_tensor.trans(iC)
657 iCTxiCT = cross_ijkl(iCT,iCT)
659 #diCdC = dinv_sym(C, iA= iC)
661 #diCTdC = dinvt(C, iCT)
662 if C.shape[-2:] == (3,3) or \
663 (C.shape[-2:] == (2,2) and self.hyp == MaterialHypothesis.PLANESTRAIN) :
664 III = det(C).reshape( C.shape[:-2] + (1,1,1,1))
665 logIII = np.log(III)
666 # kappa = 0.5*lamb*logIII- mu
667 # IIIdkappadIII = 0.5*lamb
668 # dkappa = IIIdkappadIII
669 return (0.5*lamb)*iCTxiCT + (0.5*lamb*logIII- mu)*diCTdC
670 elif C.shape[-2:] == (2,2) and self.hyp == MaterialHypothesis.PLANESTRESS:
671 II_2D = det(C).reshape(C.shape[:-2]+(1,1,1,1))
672 C33 = self.C33_planestress(II_2D)
673 III = C33*II_2D
674 logIII = np.log(III)
675 a = 2*mu*C33+lamb
676 kappa = mu*C33*(lamb*(logIII-1) - 2*mu)/a
677 II_2DxdkappadII_2D = mu*lamb*C33*( a**2 - 2*mu*lamb*(C33-1) -lamb**2*(logIII) ) /a**3
678 dkappa = II_2DxdkappadII_2D
679 #II_2DxdkappadII_2D = mu*lamb*C33*( a**2 - 2*mu*lamb*(C33-1) -lamb**2*(logIII) ) /a**3
680 return dkappa*iCTxiCT + kappa*diCTdC
681 #return dkappa*iCxiC + kappa*diCdC
683 else :
684 raise ValueError("incorrect size for C (..., 2,2) or (..., (3,3))")
686 def Cauchy(self, b:FloatArray, out_sigma:FloatArray|None = None)->FloatArray:
687 lamb, mu = (self.lamb, self.mu)
688 if b.shape[-2:] == (2,2) and self.hyp == MaterialHypothesis.PLANESTRAIN:
689 if out_sigma is None :
690 out_sigma = np.empty(b.shape[:-2] + (2,2))
691 _NH_cauchy_plane_strain_guvec(b, lamb, mu, out_sigma)
692 return out_sigma
693 if b.shape[-2:] == (3,3):
694 dim = b.shape[-1]
695 J = np.sqrt(det(b).reshape( b.shape[:-2] + (1,1)))
696 Id = np.eye(dim).reshape( (1,)*(b.ndim-2) + (dim, dim) )
697 elif b.shape[-2:] == (2,2) and self.hyp == MaterialHypothesis.PLANESTRESS :
698 II_2D = det(b).reshape(b.shape[:-2]+(1,1))
699 Id = np.eye(2).reshape ((1,)*(b.ndim-2) + (2,2))
700 b33 = self.C33_planestress(II_2D)
701 J = np.sqrt(b33*II_2D)
702 else :
703 raise NotImplementedError(f"Cauchy not implemented for \
704 {str(b.shape[-2:])} and {self.hyp.name}" )
705 return mu/J*(b-Id) + lamb*np.log(J)/J*Id
707 def spacial_elasticty_tensor(self, b:FloatArray, out_c: FloatArray|None = None)->FloatArray:
708 """
709 Compute the spatial (Eulerian) elasticity tensor for a Compressible Neo-Hookean material.
711 This tensor (c) relates the Lie derivative of the Kirchhoff stress to the
712 rate of deformation tensor. It is derived from the push-forward of the
713 material elasticity tensor to the current configuration.
715 The formula used is:
716 c = (lamb / J) * (I ⊗ I) + 2 * (mu - lamb * ln(J)) / J * Is
717 where J = det(F) = sqrt(det(b)), I ⊗ I is the identity for the volumetric
718 part, and Is is the symmetric fourth-order identity tensor.
720 Args:
721 b (ndarray): Left Cauchy-Green deformation tensor (Finger tensor).
722 Shape: (dim, dim) or (batch_dims, dim, dim).
723 Supported dims: 3D or 2D (under Plane Strain hypothesis).
725 Returns:
726 ndarray: Fourth-order spatial elasticity tensor.
727 Shape: (dim, dim, dim, dim) or (batch_dims, dim, dim, dim, dim).
729 Raises:
730 NotImplementedError: If the dimensions are not (3,3) or if (2,2) is
731 provided without Plane Strain hypothesis.
732 """
733 lamb, mu = (self.lamb, self.mu)
734 if (b.shape[-2:] == (2,2) and self.hyp == MaterialHypothesis.PLANESTRAIN) :
735 if out_c is None:
736 batch_shape = b.shape[:-2]
737 out_c = np.empty(batch_shape + (2,2,2,2), dtype=b.dtype)
738 _NH_spacial_elasticity_plane_strain_guvec(b, self.lamb, self.mu, \
739 out=out_c) # type: ignore
740 return out_c
741 if b.shape[-2:] == (3,3):
742 nbtens = b.shape[:-2]
743 bcshape = (1,)*len(nbtens)
744 J = np.sqrt(det(b)).reshape((-1,1,1,1,1))
745 IxI = Sp4[3].reshape(bcshape+(3,3,3,3))
746 Is = Sy4[3].reshape(bcshape+(3,3,3,3))
747 c = lamb/J*IxI + 2*(mu-lamb*np.log(J))/J*Is
748 if b.ndim == 2 :
749 return c.squeeze()
750 return c
751 else :
752 raise NotImplementedError()
754 def get_properties(self):
755 """ return material properies mu and lambda """
756 return np.array([self.mu, self.lamb], dtype = npfloat)
758 @staticmethod
759 @njit(inline='always', fastmath=True, cache=True)
760 def sigma_c_kernel(F, props, out_sigma, out_c):
761 """
762 Calcule la contrainte de Cauchy et le tenseur d'élasticité spatial 2D.
763 Zéro allocation, remplissage en place.
764 """
765 mu = props[0]
766 lamb = props[1]
767 # --- Cinématique ---
768 # J = det(F)
769 J = F[0, 0] * F[1, 1] - F[0, 1] * F[1, 0]
770 # b = F * F^T (Tenseur de Cauchy-Green gauche)
771 b00 = F[0, 0]*F[0, 0] + F[0, 1]*F[0, 1]
772 b01 = F[0, 0]*F[1, 0] + F[0, 1]*F[1, 1]
773 b11 = F[1, 0]*F[1, 0] + F[1, 1]*F[1, 1]
774 fill_nh_sigma_2d(J, b00, b11, b01, mu, lamb, out_sigma)
775 fill_nh_c_2d(J, mu, lamb, out_c)
777def check_material_law_F(law:HyperElasticMaterial, F:FloatArray, verb = True)->Tuple[bool, bool]:
778 ''' check PKI(F) and dPKIdF(F) return PKI_ok and dPKIDF_ok
779 PKI_ok and dPKIdF_ok, respectively True if PKI (dPKIdF) is close to
780 the finite difference approximation, False otherwise.
781 '''
782 # pylint: disable=C0103
783 #law = law
784 PKI = law.PKI(F)
785 PKI_num = df_num(law.potentialF, F, F.shape[-2:], ())
786 PKI_ok = check_relerr(PKI, PKI_num, PKI.shape[-2:] )
787 dPKIdF = law.dPKIdF(F)
788 dPKIdF_num = df_num(law.PKI, F, F.shape[-2:], F.shape[-2:])
789 dPKIdF_ok = check_relerr(dPKIdF, dPKIdF_num, PKI.shape[-2:] + F.shape[-2:], verbose= verb)
790 return PKI_ok, dPKIdF_ok
792def check_material_law_C(law:HyperElasticMaterial, C:FloatArray, verb = True)->Tuple[bool, bool]:
793 ''' check PKII(C) and dPKIIdC(C) return PKII_ok and dPKIIDC_ok
794 PKII_ok and dPKIIdC_ok, respectively True if PKII (dPKIIdC) is close to
795 the finite difference approximation, False otherwise.
796 '''
797 # pylint: disable=C0103
798 #verb = True
799 PKII = law.PKII(C)
800 shape = C.shape[-2:]
801 PKII_num = 2.*df_num(law.potentialC, C, shape, ())
802 PKII_ok = check_relerr(PKII, PKII_num, PKII.shape[-2:] , verbose= verb)
803 dPKIIdC = law.dPKIIdC(C)
804 dPKIIdC_num = df_num(law.PKII, C, shape, shape)
805 #dPKIIdC_num2 = df_num(lambda C : law.PKII( C.swapaxes(-1,-2)), C, shape, shape)
806 #dPKIIdC_num = 0.5*(dPKIIdC_num + dPKIIdC_num2)
807 dPKIIdC_ok = check_relerr(dPKIIdC, dPKIIdC_num, PKII.shape[-2:] + C.shape[-2:], verbose = verb)
808 return PKII_ok, dPKIIdC_ok
810def benchmark(law:HyperElasticMaterial, n:int, mem = False, prof = True):
811 """ benchmark : compute and profile law.PKI and law.dPKIDF """
812 # pylint: disable=C0103
813 dim = HYP2DIM[law.hyp]
814 if law.hyp in [PLANESTRESS, PLANESTRAIN] and mem:
815 F = random_lintransform2dto3d(n)
816 print(f'computing PKI and dPKIDF for {n:.3e} F32, in in {law.hyp}')
817 else :
818 gen = random_lintransform
819 F = gen(n, dim)
820 print(f'computing PKI and dPKIDF for {n:.3e} F, in {law.hyp}')
821 Profile.doprofile = prof
822 Profile.enable()
823 start = time.time()
824 __ = law.PKI(F)
825 __ = law.dPKIdF(F)
826 end = time.time()
827 Profile.disable()
828 Profile.print_stats()
829 print('time : ', end - start)
831if __name__ == "__main__":
832 lawb = CompressibleNeoHookeen(E = 1., nu = 0.3, hyp = PLANESTRESS)
833 benchmark(lawb, 1000000, prof=True)
835 # computing PKI and dPKIDF for 1.000e+06 F, in MaterialHypothesis.DIM3
836 # 550 function calls (532 primitive calls) in 1.540 seconds
837 # lawPSTRAIN = CompressibleNeoHookeen(E=1.2, nu = 0.3, hyp=PLANESTRAIN)
838 # lawPSTRESS = CompressibleNeoHookeen(E=1.2, nu = 0.3, hyp=PLANESTRESS)
839 # # benchmark(lawPSTRESS, 10000, mem =True, prof =True)
840 # law1 = SimpleLaw0()
841 # law2 = SimpleLaw1()
842 # law3 = SimpleLaw2()
843 # F2d0 = np.array([[2.,0.1],[1.,1.3]])
844 # F2d1 = np.array([[[2.,0.1],[1.,1.3]]])
845 # F2d2 = np.array([[[2.,0.1],[1.,1.3]],
846 # [[2.,0.5],[3.,1.1]]
847 # ])
848 # Fm = np.array( [
849 # [[1.,0.],[0.,1.],[0,0]],
850 # [[1.,0.],[0.,0.],[0,1]]
851 # ])
852 # #Fm = np.array([[1.,0.],[0.,1.],[0,0]])
853 # Fm2d = np.array([[2.,0.1],[1.,1.3],[0.,0.]])
854 # Fm = np.array([[[2.,0.], [0.,1.], [0,0]],
855 # [[2.,0.1],[0.1,1.],[0,0]]])
856 # Fm0= np.array([[[3.,0.1], [1., 0.1], [1.,1.]],
857 # [[2.,0.1],[0.1,1.],[0.1,0.2]]])
858 # Fm1= np.array([[[3.,0.1], [1., 0.1], [1.,1.]],
859 # [[2.,0.1],[0.1,1.],[0.1,0.2]]])
860 # lawu = lawPSTRESS
861 # Fu = Fm1
862 # print(test_material_law(lawu, Fu))
864# computing PKI and dPKIDF for 1.000e+06 F, in MaterialHypothesis.DIM3
865# 62 function calls in 1.243 seconds
866# Ordered by: cumulative time
867# ncalls tottime percall cumtime percall filename:lineno(function)
868# 1 0.013 0.013 1.031 1.031 hyperelasticity.py:221(dPKIdF)
869# 1 0.682 0.682 0.865 0.865 hyperelasticity.py:186(_dPKIdF_exp)
870# 2 0.035 0.017 0.310 0.155 smallsize_tensor.py:69(det_invt)
871# 2 0.262 0.131 0.262 0.131 smallsize_tensor.py:25(_cofac)
872# 1 0.000 0.000 0.211 0.211 hyperelasticity.py:214(PKI)
873# 1 0.003 0.003 0.211 0.211 hyperelasticity.py:164(_PKI_3d)
874# 1 0.182 0.182 0.182 0.182 smallsize_tensor.py:79(cross_ijkl)
875# 1 0.052 0.052 0.052 0.052 hyperelasticity.py:171(_PKI_exp)
876# 2 0.014 0.007 0.014 0.007 smallsize_tensor.py:41(_det)