Coverage for pyxfem/solidmechanics/largestrain.py: 15%

1""" 

2Continuum Mechanics : Tensor Transformations and Acoustic Analysis Module. 

3 

4This module provides a suite of functions to perform pull-back and push-forward 

5operations between different stress measures (Cauchy, PK1, PK2) and their 

6associated elasticity tangents. It also includes tools for computing acoustic 

7tensors and wave speeds in finite strain frameworks. 

8 

9Supported measures: 

10 - Deformation gradients (F) and Cauchy-Green tensors (C, b). 

11 - Stress tensors: Cauchy (s), First Piola-Kirchhoff (P), Second Piola-Kirchhoff (S). 

12 - Elasticity tensors: Material (C), Spatial (c), and Mixed (dPdF). 

13 

14Compatibility: 

15 All functions (but the acoustic tensors) are vectorized. 

16 They support single tensors of shape (dim, dim) 

17 or batches of tensors of shape (..., dim, dim) using NumPy broadcasting 

18 via the `smallsize_tensor` backend. 

19""" 

20# pylint: disable = C0103 

21import numpy as np 

22from typing import Optional 

23from numba import guvectorize, float64, float32, objmode 

24import numba as nb 

25 

26from pyxfem.core.config import FASTMATH, FloatArray 

27from pyxfem.core.smallsize_tensor import metric, sqrt22, trans, inv, det, \ 

28 det_invt, cofac, cross_ikjl, Id2, eigval22sym, AAT22, sqrt22_kernel 

29 

30 

31def rightCauchyGreen(F): 

32 """ 

33 Compute the Right Cauchy-Green deformation tensor C = F^T @ F. 

34 

35 Args: 

36 F (ndarray): Deformation gradient tensor(s) of shape (..., dim, dim). 

37 Returns: 

38 ndarray: The material deformation tensor(s) C. 

39 """ 

40 return metric(F) 

41 

42def leftCauchyGreen(F:FloatArray, out:FloatArray|None = None) -> FloatArray: 

43 """ 

44 Compute the Left Cauchy-Green deformation tensor b = F @ F^T. 

45 

46 Args: 

47 F (ndarray): Deformation gradient tensor(s) of shape (..., dim, dim). 

48 Returns: 

49 ndarray: The spatial deformation tensor(s) b. 

50 """ 

51 if out is None : 

52 out = np.empty_like(F) 

53 if F.shape[-2:] == (2,2): 

54 return AAT22(F, out) 

55 out[:] = F@trans(F) 

56 return out 

57 

58def PKII_to_PKI(PKII, F): 

59 """ 

60 Transform Second Piola-Kirchhoff stress (S) to First Piola-Kirchhoff stress (P). 

61 Formula: P = F @ S 

62 

63 Args: 

64 PKII (ndarray): PK2 stress tensor(s) of shape (..., dim, dim). 

65 F (ndarray): Deformation gradient tensor(s) of shape (..., dim, dim). 

66 Returns: 

67 ndarray: First Piola-Kirchhoff stress tensor(s). 

68 """ 

69 return F@PKII 

70 

71def PKI_to_PKII(PKI, F): 

72 """ 

73 Transform First Piola-Kirchhoff stress (P) to Second Piola-Kirchhoff stress (S). 

74 Formula: S = inv(F) @ P 

75 

76 Args: 

77 PKI (ndarray): PK1 stress tensor(s) of shape (..., dim, dim). 

78 F (ndarray): Deformation gradient tensor(s) of shape (..., dim, dim). 

79 Returns: 

80 ndarray: Second Piola-Kirchhoff stress tensor(s). 

81 """ 

82 return inv(F)@PKI 

83 

84def PKI_to_Cauchy(PKI, F): 

85 """ 

86 Push-forward First Piola-Kirchhoff stress (P) to Cauchy stress (sigma). 

87 Formula: sigma = (1/J) * P @ F^T 

88 

89 Args: 

90 PKI (ndarray): PK1 stress tensor(s) of shape (..., dim, dim). 

91 F (ndarray): Deformation gradient tensor(s) of shape (..., dim, dim). 

92 Returns: 

93 ndarray: Cauchy stress tensor(s). 

94 """ 

95 J = det(F)[..., np.newaxis, np.newaxis] 

96 return PKI@trans(F)/J 

97 

98def Cauchy_to_PKI(Cauchy, F): 

99 """ 

100 Pull-back Cauchy stress (sigma) to First Piola-Kirchhoff stress (P). 

101 Formula: P = J * sigma @ F^-T = sigma @ cofac(F) 

102 

103 Args: 

104 Cauchy (ndarray): Cauchy stress tensor(s) of shape (..., dim, dim). 

105 F (ndarray): Deformation gradient tensor(s) of shape (..., dim, dim). 

106 Returns: 

107 ndarray: First Piola-Kirchhoff stress tensor(s). 

108 """ 

109 #return det(F)*Cauchy@invt(F) 

110 print("## ", Cauchy.shape, F.shape, cofac(F).shape) 

111 return Cauchy@cofac(F) 

112 

113def PKII_to_Cauchy(PKII, F): 

114 """ 

115 Push-forward Second Piola-Kirchhoff stress (S) to Cauchy stress (sigma). 

116 Formula: sigma = (1/J) * F @ S @ F^T 

117 

118 Args: 

119 PKII (ndarray): PK2 stress tensor(s) of shape (..., dim, dim). 

120 F (ndarray): Deformation gradient tensor(s) of shape (..., dim, dim). 

121 Returns: 

122 ndarray: Cauchy stress tensor(s). 

123 """ 

124 J = det(F)[..., np.newaxis, np.newaxis] 

125 return F@PKII@trans(F)/J 

126 

127def Cauchy_to_PKII(Cauchy, F): 

128 """ 

129 Pull-back Cauchy stress (sigma) to Second Piola-Kirchhoff stress (S). 

130 Formula: S = J * F^-1 @ sigma @ F^-T 

131 

132 Args: 

133 Cauchy (ndarray): Cauchy stress tensor(s) of shape (..., dim, dim). 

134 F (ndarray): Deformation gradient tensor(s) of shape (..., dim, dim). 

135 Returns: 

136 ndarray: Second Piola-Kirchhoff stress tensor(s). 

137 """ 

138 J, FmT = det_invt(F) 

139 J_view = np.reshape(J, (-1, 1, 1)) if np.ndim(J) > 0 else J 

140 return J_view*trans(FmT)@Cauchy@FmT 

141 

142def dPKIIdE_to_dPKIdF(dPKIIdE, PKII, F): 

143 """ 

144 Transform the material tangent (dS/dE) to the mixed tangent (dP/dF). 

145 

146 Args: 

147 dPKIIdE (ndarray): Fourth-order material elasticity tensor(s). 

148 PKII (ndarray): PK2 stress tensor(s). 

149 F (ndarray): Deformation gradient tensor(s). 

150 Returns: 

151 ndarray: Mixed elasticity tensor(s) dP/dF. 

152 """ 

153 dPKIdF = np.einsum("...iI, ...IJKL, ...kK->...iJkL", F, dPKIIdE, F) 

154 dim = F.shape[-1] 

155 Id = Id2[dim].reshape( (1,)*(F.ndim-2) + (dim, dim)) 

156 dPKIdF += cross_ikjl(Id, PKII) 

157 return dPKIdF 

158 

159def dPKIdF_to_dPKIIdE(dPKIdF, PKI, F): 

160 """ 

161 Transform the mixed tangent (dP/dF) to the material tangent (dS/dE). 

162 

163 Args: 

164 dPKIdF (ndarray): Mixed elasticity tensor(s) dP/dF. 

165 PKI (ndarray): PK1 stress tensor(s). 

166 F (ndarray): Deformation gradient tensor(s). 

167 Returns: 

168 ndarray: Material elasticity tensor(s) dS/dE. 

169 """ 

170 iF = inv(F) 

171 PkII = iF@PKI 

172 dim = F.shape[-1] 

173 Id = Id2[dim].reshape( (1,)*(F.ndim-2) + (dim, dim)) 

174 tmp = dPKIdF - cross_ikjl(Id, PkII) 

175 dPKIIdE = np.einsum("...Ii, ...iJkL, ...Kk->...IJKL", iF, tmp, iF) 

176 return dPKIIdE 

177 

178def dPKIIdE_to_c(dPKIIdE, F): 

179 """ 

180 Push-forward the material tangent (dS/dE) to the spatial elasticity tensor (c). 

181 

182 Args: 

183 dPKIIdE (ndarray): Material elasticity tensor(s) (dPKII/dE). 

184 F (ndarray): Deformation gradient tensor(s). 

185 Returns: 

186 ndarray: Spatial elasticity tensor(s) c. 

187 """ 

188 J = det(F).reshape((-1,1,1,1,1)) 

189 c = np.einsum("...iI, ...jJ, ...kK, ...lL, ...IJKL -> ...ijkl", 

190 F,F,F,F,dPKIIdE)/J 

191 return c 

192 

193def c_to_dPKIIdE(c, F): 

194 """ 

195 Pull-back the spatial elasticity tensor (c) to the material tangent (dS/dE). 

196 

197 Args: 

198 c (ndarray): Spatial elasticity tensor(s). 

199 F (ndarray): Deformation gradient tensor(s). 

200 Returns: 

201 ndarray: Material elasticity tensor(s) dS/dE. 

202 """ 

203 J = det(F).reshape((-1,1,1,1,1)) 

204 iF = inv(F) 

205 C = J*np.einsum("...Ii, ...Jj, ...Kk, ...Ll, ...ijkl -> ...IJKL", 

206 iF,iF,iF,iF,c) 

207 return C 

208 

209def dPKIdF_to_c(dPKIdF, PKI, F): 

210 """ 

211 Convert the mixed tangent (dP/dF) to the spatial elasticity tensor (c). 

212 

213 Args: 

214 dPKIdF (ndarray): Mixed elasticity tensor(s). 

215 PKI (ndarray): PK1 stress tensor(s). 

216 F (ndarray): Deformation gradient tensor(s). 

217 Returns: 

218 ndarray: Spatial elasticity tensor(s) c. 

219 """ 

220 J = det(F) 

221 dim = F.shape[-1] 

222 s = PKI@trans(F)/J.reshape((-1,1,1)) 

223 Id = Id2[dim].reshape( (1,)*(F.ndim-2) + (dim, dim)) 

224 c = np.einsum("...jJ, ...iJkL, ...lL->...ijkl", F, dPKIdF , F)/J.reshape((-1,1,1,1,1)) 

225 c -= cross_ikjl(Id, s) 

226 return c 

227 

228def c_dPKIdF(c, s, F): 

229 """ 

230 Convert the spatial elasticity tensor (c) to the mixed tangent (dP/dF). 

231 

232 Args: 

233 c (ndarray): Spatial elasticity tensor(s). 

234 s (ndarray): Cauchy stress tensor(s). 

235 F (ndarray): Deformation gradient tensor(s). 

236 Returns: 

237 ndarray: Mixed elasticity tensor(s) dP/dF. 

238 """ 

239 J = det(F) 

240 dim = F.shape[-1] 

241 Id = Id2[dim].reshape( (1,)*(F.ndim-2) + (dim, dim)) 

242 tmp = J.reshape((-1,1,1,1,1))*(c + cross_ikjl(Id, s)) 

243 iF = inv(F) 

244 dPdF = np.einsum("...Jj, ...Ll, ...ijkl -> ...iJkL", iF,iF, tmp) 

245 return dPdF 

246 

247def acoustic_tensor_mixte(dPdF, N, J): 

248 """ 

249 Compute the acoustic tensor using the mixed tangent dP/dF. 

250 

251 Args: 

252 dPdF (ndarray): Mixed elasticity tensor(s). 

253 N (ndarray): normal vector(s) in reference config (..., dim). 

254 N is the piola pull back of a unit normal vector in current configuration 

255 N = (1/J) F^T n with |n| = 1 

256 J (float/ndarray): Determinant of F. 

257 Returns: 

258 ndarray: The acoustic tensor(s) of shape (..., dim, dim). 

259 Qik = J AiJkL NJ N L 

260 """ 

261 return np.einsum( "I,iIjJ,J->ij", N, dPdF, N)*J 

262 

263def acoustic_tensor_space_v0(c, s, n): 

264 """ 

265 Compute the acoustic tensor using spatial variables (c, sigma). 

266 

267 Args: 

268 c (ndarray): Spatial elasticity tensor(s). 

269 s (ndarray): Cauchy stress tensor(s). 

270 n (ndarray): Direction vector(s) in current config (..., dim). 

271 Returns: 

272 ndarray: The spatial acoustic tensor(s). 

273 """ 

274 l = len(n.shape)-1 

275 dim = n.shape[-1] 

276 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] 

277 

278 

279 

280def acoustic_tensor_space_v1(c, s, n, out=None): 

281 """ Version déroulée spécifiquement pour la 2D """ 

282 if n.shape[-1] != 2: 

283 raise NotImplementedError() 

284 if out is None: 

285 out = np.empty(n.shape+(2,)) 

286 #out = np.empty_like(s) 

287 

288 n0, n1 = n[..., 0], n[..., 1] 

289 

290 # --- Calcul de la contrainte normale (scalaire) --- 

291 # sigma_n = n0*(s00*n0 + s01*n1) + n1*(s10*n0 + s11*n1) 

292 # Note : On utilise des vues temporaires, Python gère ça très vite 

293 sigma_n = n0 * (s[..., 0, 0]*n0 + s[..., 0, 1]*n1) + \ 

294 n1 * (s[..., 1, 0]*n0 + s[..., 1, 1]*n1) 

295 

296 # --- Calcul de A00 --- 

297 out[..., 0, 0] = (c[..., 0, 0, 0, 0] * n0 * n0 + 

298 c[..., 0, 1, 0, 0] * n1 * n0 + 

299 c[..., 0, 0, 0, 1] * n0 * n1 + 

300 c[..., 0, 1, 0, 1] * n1 * n1) + sigma_n 

301 

302 # --- Calcul de A01 --- 

303 out[..., 0, 1] = (c[..., 0, 0, 1, 0] * n0 * n0 + 

304 c[..., 0, 1, 1, 0] * n1 * n0 + 

305 c[..., 0, 0, 1, 1] * n0 * n1 + 

306 c[..., 0, 1, 1, 1] * n1 * n1) 

307 

308 # --- Calcul de A10 --- 

309 out[..., 1, 0] = (c[..., 1, 0, 0, 0] * n0 * n0 + 

310 c[..., 1, 1, 0, 0] * n1 * n0 + 

311 c[..., 1, 0, 0, 1] * n0 * n1 + 

312 c[..., 1, 1, 0, 1] * n1 * n1) 

313 

314 # --- Calcul de A11 --- 

315 out[..., 1, 1] = (c[..., 1, 0, 1, 0] * n0 * n0 + 

316 c[..., 1, 1, 1, 0] * n1 * n0 + 

317 c[..., 1, 0, 1, 1] * n0 * n1 + 

318 c[..., 1, 1, 1, 1] * n1 * n1) + sigma_n 

319 

320 return out 

321 

322def acoustic_tensor_space(c: np.ndarray, s: np.ndarray, n: np.ndarray, out=None) -> np.ndarray: 

323 # ... 

324 if out is None: 

325 batch_out = np.broadcast_shapes(c.shape[:-4], s.shape[:-2], n.shape[:-1]) 

326 out_dtype = n.dtype 

327 out = np.empty(batch_out + (2, 2), dtype=out_dtype) 

328 _acoustic_tensor_2d_guvec(c, s, n, out=out) 

329 

330 return out 

331 

332#acoustic_tensor_space = acoustic_tensor_space 

333 

334@nb.njit(cache=True, fastmath = FASTMATH) 

335def acoustic_tensor_2d_kernel_njit(c, s, n, out): 

336 n0 = n[0] 

337 n1 = n[1] 

338 sigma_n = n0 * (s[0, 0]*n0 + s[0, 1]*n1) + \ 

339 n1 * (s[1, 0]*n0 + s[1, 1]*n1) 

340 

341 out[0, 0] = (c[0, 0, 0, 0] * n0 * n0 + c[0, 1, 0, 0] * n1 * n0 + 

342 c[0, 0, 0, 1] * n0 * n1 + c[0, 1, 0, 1] * n1 * n1) + sigma_n 

343 out[0, 1] = (c[0, 0, 1, 0] * n0 * n0 + c[0, 1, 1, 0] * n1 * n0 + 

344 c[0, 0, 1, 1] * n0 * n1 + c[0, 1, 1, 1] * n1 * n1) 

345 out[1, 0] = (c[1, 0, 0, 0] * n0 * n0 + c[1, 1, 0, 0] * n1 * n0 + 

346 c[1, 0, 0, 1] * n0 * n1 + c[1, 1, 0, 1] * n1 * n1) 

347 out[1, 1] = (c[1, 0, 1, 0] * n0 * n0 + c[1, 1, 1, 0] * n1 * n0 + 

348 c[1, 0, 1, 1] * n0 * n1 + c[1, 1, 1, 1] * n1 * n1) + sigma_n 

349 

350@guvectorize([ 

351 (float64[:,:,:,:], float64[:,:], float64[:], float64[:,:]), 

352 (float32[:,:,:,:], float32[:,:], float32[:], float32[:,:]) 

353], '(m,m,m,m),(m,m),(m)->(m,m)', nopython=True, cache=True) 

354def _acoustic_tensor_2d_guvec(c, s, n, out): 

355 acoustic_tensor_2d_kernel_njit(c, s, n, out) 

356 

357 

358 

359 

360 

361def acoustic_impedance_tensor(rho, c, s, n, outZ = None): 

362 """ compute the acoustic tensor impedance in direction n 

363 c is the 4 order tangent spacial material tensor 

364 s is the cauchy stress 

365 n is a direction (norm 1 ..) 

366 rho is the density 

367 return OutZ : acoustic tensor in direction n 

368 """ 

369 #print(rho.shape, c.shape, s.shape ,n.shape, outZ.shape) 

370 

371 #print(c.shape, s.shape, n.shape) 

372 if outZ is None : 

373 batch_out = np.broadcast_shapes(c.shape[:-4], s.shape[:-2], n.shape[:-1]) 

374 out_dtype = n.dtype 

375 outZ = np.empty(batch_out + (2, 2), dtype=out_dtype) 

376 

377 acoustic_tensor_space(c,s,n, outZ) 

378 #print("G", G.shape) 

379 outZ *= rho[..., None, None] 

380 return sqrt22(outZ, outZ) 

381 

382 

383 #Z = sqrt22(rho.reshape(rho.shape+(1,1))*A) 

384 #return outZ 

385 

386@nb.njit(parallel = True, cache = True, fastmath = FASTMATH) 

387def acoustic_impedance_tensor_pnjit(rho, c, s, n, out): 

388 #print(rho.shape, c.shape, s.shape ,n.shape, out.shape) 

389 N = rho.shape[0] 

390 assert(N == c.shape[0] and s.shape[0] == N and n.shape[0] == N and out.shape[0] == N) 

391 assert( len(n.shape) == 3) # typically we have more than one normal for each 

392 for j in nb.prange(N): 

393 jrho = rho[j, 0] 

394 jc = c[j, 0] 

395 js = s[j, 0] 

396 jvn = n[j] 

397 jvout = out[j] 

398 for r in range(n.shape[1]) : 

399 jrout = jvout[r] 

400 acoustic_tensor_2d_kernel_njit(jc, js, jvn[r], jrout ) 

401 jrout[0,0]*=jrho 

402 jrout[0,1]*=jrho 

403 jrout[1,0]*=jrho 

404 jrout[1,1]*=jrho 

405 sqrt22_kernel(jrout, jrout) 

406 return out 

407 

408def acoustic_tensor_material(C, S, N, F): 

409 """ 

410 Compute the acoustic tensor using material variables (C, S). 

411 Args: 

412 C (ndarray): Material elasticity tensor(s). 

413 S (ndarray): PK2 stress tensor(s). 

414 N (ndarray): normal vector(s) in reference config (..., dim). 

415 N is the piola pull back of a unit normal vector in current configuration 

416 N = (1/J) F^T n with |n| = 1 

417 F (ndarray): Deformation gradient tensor(s). 

418 Returns: 

419 ndarray: The material acoustic tensor(s). 

420 """ 

421 

422 dim = len(N) 

423 J =det(F) 

424 

425 A = J*(np.einsum( "iI,J,IJKL,kK,L -> ik", F, N, C, F, N ) + N.dot(S@N)*Id2[dim]) 

426 return A 

427 

428def wavespeed(A, rho): 

429 """ 

430 Calculate wave speeds from the acoustic tensor and density. 

431 For now only work for 2x2 acoustic tensor .... 

432 Args: 

433 A (ndarray): Acoustic tensor(s). 

434 rho (float/ndarray): Material density. 

435 Returns: 

436 ndarray: Wave speeds (sqrt of eigenvalues of A/rho). 

437 """ 

438 if A.shape[-1] != 2 : 

439 raise NotImplementedError("Wavespeed only works for acoustic tensor in 2d space") 

440 eigs = eigval22sym(A) 

441 np.sqrt(eigs/rho) 

442 return eigs 

443 

444 

445def test_transport_formula(): 

446 """ checking the transformation formulas """ 

447 # pylint: disable=import-outside-toplevel 

448 from pyxfem.solidmechanics.hyperelasticity import CompressibleNeoHookeen 

449 from pyxfem.solidmechanics.elasticity import MaterialHypothesis 

450 law_ps = CompressibleNeoHookeen(lamb = 1., mu=0.5, hyp=MaterialHypothesis.PLANESTRAIN) 

451 law_3D = CompressibleNeoHookeen(lamb = 1., mu=0.5, hyp=MaterialHypothesis.DIM3) 

452 

453 F_2D = np.array([[2.,0.1], [0.2, 1.]]) 

454 array_of_F_2D = np.array([[[2.,0.3], [0.22, 1.]], 

455 [[1.2,0.3], [0.21, 0.8]], 

456 [[2.,0.3], [0.22, 1.]]] ) 

457 F_3D = np.array([[2.,0.1, -0.1], [0.05, 1.1, 0.2], [0.15, -0.13, 0.8] ]) 

458 array_of_F_3D = np.array([[[2.,0.1, -0.1], [0.05, 1.1, 0.2], [0.15, -0.13, 0.8] ], 

459 [[3.,0.5, -0.1], [0.05, 1.1, 0.2], [-0.15, -0.13, 0.8] ] ]) 

460 

461 

462 for F, name, law in zip([F_2D, array_of_F_2D, F_3D, array_of_F_3D], 

463 ["F 2D", "array of F 2D", "F3D", "array of F 3D"], 

464 [law_ps, law_ps, law_3D, law_3D]): 

465 print("### Start test for ", name ) 

466 C = rightCauchyGreen(F) 

467 b = leftCauchyGreen(F) 

468 P = law.PKI(F) 

469 S = law.PKII(C) 

470 s = law.Cauchy(b) 

471 print(P, S, s) 

472 print(s.shape) 

473 P2 = PKII_to_PKI(S, F) 

474 print('S-> P error ', np.linalg.norm(P-P2)) 

475 P3 = Cauchy_to_PKI(s, F) 

476 print('s-> P error ', np.linalg.norm(P-P3)) 

477 S2 = PKI_to_PKII(P, F) 

478 print('P-> S error ', np.linalg.norm(S-S2)) 

479 S3 = Cauchy_to_PKII(s, F) 

480 print('s-> S error ', np.linalg.norm(S-S3)) 

481 s2 = PKI_to_Cauchy(P,F) 

482 s3 = PKII_to_Cauchy(S,F) 

483 print('P-> s error ', np.linalg.norm(s-s2)) 

484 print('S-> s error ', np.linalg.norm(s-s3)) 

485 

486 dPdF = law.dPKIdF(F) 

487 dSdE = 2.*law.dPKIIdC(C) 

488 c = law.spacial_elasticty_tensor(b) 

489 

490 dPdF2 = dPKIIdE_to_dPKIdF(dSdE, S, F) 

491 print("dSdE-> dPdF error", np.linalg.norm(dPdF- dPdF2) ) 

492 dPdF3 = c_dPKIdF(c, s, F) 

493 print("c-> dPdF error", np.linalg.norm(dPdF- dPdF3) ) 

494 dSdE2 = dPKIdF_to_dPKIIdE(dPdF, P, F) 

495 print("dPdF-> dSdE error", np.linalg.norm(dSdE- dSdE2) ) 

496 dSdE3 = c_to_dPKIIdE(c, F) 

497 print("c-> dSdE error", np.linalg.norm(dSdE- dSdE3) ) 

498 c2 = dPKIIdE_to_c(dSdE, F) 

499 print("dSdE-> c error", np.linalg.norm(c- c2) ) 

500 c3 = dPKIdF_to_c(dPdF, P, F) 

501 print("dPdF-> c error", np.linalg.norm(c- c3) ) 

502 

503def test_acoustic_tensor(): 

504 """ Testing Acoustic tensor and wave seed computation """ 

505 # pylint: disable=import-outside-toplevel 

506 from pyxfem.solidmechanics.hyperelasticity import CompressibleNeoHookeen 

507 from pyxfem.solidmechanics.elasticity import MaterialHypothesis 

508 law = CompressibleNeoHookeen(lamb = 3., mu= 7., hyp=MaterialHypothesis.PLANESTRAIN) 

509 rho0 = 1. 

510 #F = np.array([[1.,0.], [0., 1.]]) 

511 F = np.array([[3.,1.], [0.8, 2]]) 

512 n = np.array([1.,0.8]) 

513 n /= np.linalg.norm(n) 

514 dPdF = law.dPKIdF(F) 

515 J = det(F) 

516 N = trans(F)@n/J 

517 A1 = acoustic_tensor_mixte(dPdF, N, J) 

518 #print("A1", A1) 

519 

520 b = leftCauchyGreen(F) 

521 s = law.Cauchy(b) 

522 

523 

524 c = law.spacial_elasticty_tensor(b) 

525 A2 = acoustic_tensor_space(c, s, n) 

526 

527 S = law.PKII(rightCauchyGreen(F)) 

528 C = law.dPKIIdE(rightCauchyGreen(F)) 

529 A3 = acoustic_tensor_material(C, S, N, F) 

530 print("A1", A1) 

531 print("A2", A2) 

532 print("A3", A3) 

533 

534 J = det(F) 

535 rho = rho0/J 

536 

537 print('wavespeed', wavespeed(A1, rho)) 

538 print('wavespeed', wavespeed(A2, rho)) 

539 print('wavespeed', wavespeed(A3, rho)) 

540 

541 print("test Acoustic tensor for multiple direction") 

542 ns = np.array([[1.,0.],[0,1], [np.sqrt(2)/2, np.sqrt(2)/2], [n[0], n[1]]]) 

543 As = acoustic_tensor_space(c[np.newaxis], s[np.newaxis], ns) 

544 

545 

546 

547 print(As.shape) 

548 print(As) 

549 c = np.stack((c, 2*c), axis=0) 

550 s = np.stack((s, 2*s), axis=0) 

551 

552 

553 As = acoustic_tensor_space(c[:,np.newaxis], s[:,np.newaxis], ns[np.newaxis]) 

554 print("2 states, 4 vectors ") 

555 print(As.shape) 

556 print(As) 

557 rho = np.array([rho, rho]) 

558 print('wavespeed', wavespeed(As, rho[np.newaxis,:])) 

559 

560 

561if __name__ == "__main__": 

562 test_transport_formula() 

563 test_acoustic_tensor()