Coverage for physioblocks / library / model_components / dynamics.py: 100%

1# SPDX-FileCopyrightText: Copyright INRIA 

2# 

3# SPDX-License-Identifier: LGPL-3.0-only 

4# 

5# Copyright INRIA 

6# 

7# This file is part of PhysioBlocks, a library mostly developed by the 

8# [Ananke project-team](https://team.inria.fr/ananke) at INRIA. 

9# 

10# Authors: 

11# - Colin Drieu 

12# - Dominique Chapelle 

13# - François Kimmig 

14# - Philippe Moireau 

15# 

16# PhysioBlocks is free software: you can redistribute it and/or modify it under the 

17# terms of the GNU Lesser General Public License as published by the Free Software 

18# Foundation, version 3 of the License. 

19# 

20# PhysioBlocks is distributed in the hope that it will be useful, but WITHOUT ANY 

21# WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A 

22# PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details. 

23# 

24# You should have received a copy of the GNU Lesser General Public License along with 

25# PhysioBlocks. If not, see <https://www.gnu.org/licenses/>. 

26 

27""" 

28Describes dynamics models. 

29""" 

30 

31from dataclasses import dataclass 

32from typing import Any 

33 

34import numpy as np 

35from numpy.typing import NDArray 

36 

37import physioblocks.utils.math_utils as math_utils 

38from physioblocks.computing import Expression, ModelComponent, Quantity, diff, mid_point 

39from physioblocks.computing.assembling import EqSystem 

40from physioblocks.registers import register_type 

41from physioblocks.simulation import State, Time 

42from physioblocks.simulation.solvers import NewtonSolver 

43 

44# Constant for the spherical dynamics type id 

45SPHERICAL_DYNAMICS_TYPE_ID = "spherical_dynamics" 

46 

47SPHERICAL_DYNAMICS_STATIC_DISP_LOCAL_ID = "disp" 

48 

49_DISP_EPSILON = 1.0e-6 

50_STATIC_PROBLEM_TOL = 1e-6 

51_STATIC_PROBLEM_MAX_IT = 10 

52_STATIC_PROBLEM_DISP_MAG = 0.010 

53_STATIC_PROBLEM_MIN_PRESSURE_STEP = 1.0 

54 

55 

56@dataclass 

57class _SphericalDynamicsStaticModelComponent(ModelComponent): 

58 """ 

59 ModelComponent containing the quantities and the expressions to solve the 

60 static problem in the Spherical Dynamics Model in order to initialize the 

61 displacement and fiber deformation. 

62 """ 

63 

64 disp: Quantity[np.float64] 

65 """Displacement""" 

66 

67 pressure: Quantity[np.float64] 

68 """pressure""" 

69 

70 thickness: Quantity[np.float64] 

71 """thickness""" 

72 

73 hyperelastic_cst: Quantity[NDArray[np.float64]] 

74 """hyperelastic constant""" 

75 

76 inv_radius: Quantity[np.float64] 

77 """sphere radius""" 

78 

79 pressure_external: Quantity[np.float64] 

80 """external pressure""" 

81 

82 def initialize(self) -> None: 

83 """ 

84 Initialize block attributes from current quantity values 

85 """ 

86 self.hyperelastic_cst_01 = ( 

87 self.hyperelastic_cst.current[0] * self.hyperelastic_cst.current[1] 

88 ) 

89 self.hyperelastic_cst_23 = ( 

90 self.hyperelastic_cst.current[2] * self.hyperelastic_cst.current[3] 

91 ) 

92 self.thickness_radius_ratio = self.thickness.current * self.inv_radius.current 

93 self.half_thickness_radius_ratio = 0.5 * self.thickness_radius_ratio 

94 

95 def dynamics_static_residual(self) -> Any: 

96 """ 

97 Compute the residual for the dynamics static problem. 

98 

99 :return: the static problem residual 

100 :rtype: np.float64 

101 """ 

102 pressure_external = mid_point(self.pressure_external) 

103 

104 disp_new_ratio = 1.0 + self.disp.new * self.inv_radius.current 

105 fluid_volume_diff = np.pow( 

106 disp_new_ratio 

107 - np.pow(disp_new_ratio, -2) * self.half_thickness_radius_ratio, 

108 2, 

109 ) * (1.0 + np.pow(disp_new_ratio, -3) * self.thickness_radius_ratio) 

110 

111 j4_new = np.pow(disp_new_ratio, 2) 

112 j1_new = 2.0 * j4_new + np.pow(disp_new_ratio, -4) 

113 

114 j4_diff_new = 2.0 * self.inv_radius.current * (disp_new_ratio) 

115 j1_diff_new = 2.0 * ( 

116 j4_diff_new - 2.0 * self.inv_radius.current * np.pow(disp_new_ratio, -5) 

117 ) 

118 

119 hyperelastic_potential_diff = 2.0 * ( 

120 ( 

121 self.hyperelastic_cst_01 

122 * j1_diff_new 

123 * (j1_new - 3.0) 

124 * np.exp(self.hyperelastic_cst.current[1] * np.pow(j1_new - 3.0, 2)) 

125 ) 

126 + ( 

127 self.hyperelastic_cst_23 

128 * j4_diff_new 

129 * (j4_new - 1.0) 

130 * np.exp(self.hyperelastic_cst.current[3] * np.pow(j4_new - 1.0, 2)) 

131 ) 

132 ) 

133 

134 passive_stress = self.thickness.current * hyperelastic_potential_diff 

135 external_stress = (self.pressure.new - pressure_external) * fluid_volume_diff 

136 

137 return passive_stress - external_stress 

138 

139 def dynamics_static_residual_ddisp(self) -> Any: 

140 """ 

141 Compute the partial derivative for disp for the the static problem residual. 

142 

143 :return: the residual partial derivative for disp 

144 :rtype: np.float64 

145 """ 

146 disp_new_ratio = 1.0 + self.disp.new * self.inv_radius.current 

147 pressure_external = mid_point(self.pressure_external) 

148 

149 j4_new = np.pow(disp_new_ratio, 2) 

150 j1_new = 2.0 * j4_new + np.pow(disp_new_ratio, -4) 

151 

152 j4_diff_new = 2.0 * self.inv_radius.current * (disp_new_ratio) 

153 j1_diff_new = 2.0 * ( 

154 j4_diff_new - 2.0 * self.inv_radius.current * np.pow(disp_new_ratio, -5) 

155 ) 

156 

157 j1_diff_diff_new = 4.0 * np.pow(self.inv_radius.current, 2) + 20.0 * np.pow( 

158 self.inv_radius.current, 2 

159 ) * np.pow(disp_new_ratio, -6) 

160 j4_diff_diff_new = 2.0 * np.pow(self.inv_radius.current, 2) 

161 

162 hyperelastic_potential_diff_diff = 2.0 * ( 

163 self.hyperelastic_cst.current[0] 

164 * self.hyperelastic_cst.current[1] 

165 * ( 

166 j1_diff_diff_new * (j1_new - 3.0) 

167 + np.pow(j1_diff_new, 2) 

168 + ( 

169 2.0 

170 * self.hyperelastic_cst.current[1] 

171 * np.pow(j1_diff_new * (j1_new - 3.0), 2) 

172 ) 

173 ) 

174 * np.exp(self.hyperelastic_cst.current[1] * np.pow(j1_new - 3.0, 2)) 

175 + self.hyperelastic_cst.current[2] 

176 * self.hyperelastic_cst.current[3] 

177 * ( 

178 j4_diff_diff_new * (j4_new - 1.0) 

179 + np.pow(j4_diff_new, 2) 

180 + ( 

181 2.0 

182 * self.hyperelastic_cst.current[3] 

183 * np.pow(j4_diff_new * (j4_new - 1.0), 2) 

184 ) 

185 ) 

186 * np.exp(self.hyperelastic_cst.current[3] * np.pow(j4_new - 1.0, 2)) 

187 ) 

188 

189 passive_stress_diff = self.thickness.current * hyperelastic_potential_diff_diff 

190 

191 fluid_volume_diff_diff = self.inv_radius.current * ( 

192 2.0 

193 * ( 

194 1.0 

195 + self.disp.new * self.inv_radius.current 

196 - np.pow(disp_new_ratio, -2) * self.half_thickness_radius_ratio 

197 ) 

198 * np.pow(1.0 + np.pow(disp_new_ratio, -3) * self.thickness_radius_ratio, 2) 

199 - 3.0 

200 * self.thickness_radius_ratio 

201 * np.pow(disp_new_ratio, -4) 

202 * np.pow( 

203 1.0 

204 + self.disp.new * self.inv_radius.current 

205 - np.pow(disp_new_ratio, -2) * self.half_thickness_radius_ratio, 

206 2, 

207 ) 

208 ) 

209 

210 ext_stress_diff = ( 

211 self.pressure.new - pressure_external 

212 ) * fluid_volume_diff_diff 

213 

214 return passive_stress_diff - ext_stress_diff 

215 

216 

217_spherical_dynamics_static_residual_expr = Expression( 

218 1, 

219 _SphericalDynamicsStaticModelComponent.dynamics_static_residual, 

220 { 

221 SPHERICAL_DYNAMICS_STATIC_DISP_LOCAL_ID: _SphericalDynamicsStaticModelComponent.dynamics_static_residual_ddisp # noqa: E501 

222 }, 

223) 

224 

225_SphericalDynamicsStaticModelComponent.declares_internal_expression( 

226 SPHERICAL_DYNAMICS_STATIC_DISP_LOCAL_ID, _spherical_dynamics_static_residual_expr 

227) 

228 

229 

230@dataclass 

231@register_type(SPHERICAL_DYNAMICS_TYPE_ID) 

232class SphericalDynamicsModelComponent(ModelComponent): 

233 r""" 

234 Implementation of the spherical dynamics model. 

235 

236 It provides a dynamics on :math:`y` the displacement of a sphere radius. 

237 

238 **Internal Equation:** 

239 

240 .. math:: 

241 

242 \frac{P - P_{ext}}{|\Omega_0|} \frac{\partial V(y)}{\partial y} 

243 - \rho_0 \ddot{y} 

244 - \frac{\partial W_e(y)}{\partial y} 

245 - \frac{1}{R_0} \Sigma_v(y,\dot{y}) 

246 - \frac{k_s}{R_0} \Bigl( \frac{y}{R_0}-e_c \Bigr) = 0 

247 

248 With: 

249 

250 * :math:`W_e(y)` represents an elastic energy per unit volume given by 

251 

252 .. math:: 

253 

254 W_e(y) = 

255 C_0 \exp \Bigl[ C_1 \bigl(2C(y)+C(y)^{-2}-3\bigr)^2 \Bigr] 

256 + C_2 \exp \Bigl[ C_3 \bigl(C(y)-1\bigr)^2 \Bigr] 

257 

258 where :math:`(C_0,C_1,C_2,C_3)` denote material parameters of the passive 

259 constitutive law, and :math:`C(y)= (1+y/R_0)^2` represents the component 

260 of the right Cauchy-Green deformation tensor along the fiber direction 

261 

262 * :math:`\Sigma_v(y,\dot{y})` denotes a viscous stress defined by 

263 

264 .. math:: 

265 

266 \Sigma_v(y,\dot{y}) = 

267 2 \eta \, \bigl( C(y) + 2 C(y)^{-5} \bigr) \frac{\dot{y}}{R_0} 

268 

269 with :math:`\eta` a solid viscosity parameter 

270 

271 * :math:`k_s` denotes a passive elastic modulus 

272 

273 * :math:`|\Omega_0|` is the total volume of sphere tissue in the reference 

274 configuration 

275 

276 .. math:: 

277 

278 |\Omega_0| = 

279 \frac{4}{3}\pi \Biggl[ \Bigl( R_0+\frac{d_0}{2} \Bigr)^3 

280 - \Bigl( R_0-\frac{d_0}{2} \Bigr)^3 \Biggr] 

281 

282 with :math:`d_0` the wall thickness in the reference configuration, and 

283 factorizing :math:`R_0` from the bracket a Taylor expansion readily gives 

284 the following approximation that we use in the implementation up to the 

285 third order in :math:`d_0/2R_0`: 

286 

287 .. math:: 

288 

289 |\Omega_0| = 4\pi R_0^2 d_0 

290 

291 * :math:`V(y)` denotes the volume of the sphere, i.e. 

292 

293 .. math:: 

294 

295 V(y) = \frac{4}{3}\pi \Bigl( R(y)-\frac{d(y)}{2} \Bigr)^3 

296 = \frac{4}{3}\pi \Bigl( R_0+y-\frac{d_0}{2}C(y)^{-1} \Bigr)^3, 

297 

298 with :math:`d(y)` the wall thickness in the deformed configuration. 

299 

300 **Discretised:** 

301 

302 .. math:: 

303 

304 \frac{P^{n + \frac{1}{2}} 

305 - P_{ext}^{n + \frac{1}{2}}}{|\Omega_0|} DV^{{n + \frac{1}{2}}\sharp} 

306 - \rho_0 \Bigl[ \frac{\dot{y}^{n+1}-\dot{y}^n}{\Delta t^n} \Bigr] 

307 - DW^{{n + \frac{1}{2}}\sharp} 

308 - \frac{1}{R_0} \Sigma_v(y^{n + \frac{1}{2}},\dot{y}^{n + \frac{1}{2}}) 

309 - \frac{k_s}{R_0} \Bigl( 

310 \frac{y^{n + \frac{1}{2}}}{R_0}-e_c^{n + \frac{1}{2}} \Bigr) = 0 

311 """ 

312 

313 disp: Quantity[np.float64] 

314 """:math:`y` the displacement along the radius.""" 

315 

316 fib_deform: Quantity[np.float64] 

317 """:math:`e_c` the sphere fiber deformation""" 

318 

319 pressure: Quantity[np.float64] 

320 """:math:`P` the pressure in the sphere""" 

321 

322 pressure_external: Quantity[np.float64] 

323 """:math:`P_{ext}` the pressure outside the sphere""" 

324 

325 vel: Quantity[np.float64] 

326 """:math:`\\dot{y}` the displacement velocity""" 

327 

328 radius: Quantity[np.float64] 

329 """:math:`R_0` the initial sphere radius""" 

330 

331 vol_mass: Quantity[np.float64] 

332 """:math:`\\rho_0` the volumic mass of the sphere tissu""" 

333 

334 thickness: Quantity[np.float64] 

335 """:math:`d_0` the initial thickness of the sphere""" 

336 

337 damping_coef: Quantity[np.float64] 

338 """Damping coeff for the viscous coeff :math:`\\eta`""" 

339 

340 series_stiffness: Quantity[np.float64] 

341 """:math:`k_s` the series stiffness""" 

342 

343 hyperelastic_cst: Quantity[NDArray[np.float64]] 

344 """:math:`(C_0,C_1,C_2,C_3)` the hyperelastic constants""" 

345 

346 time: Time 

347 """Simulation time""" 

348 

349 def initialize(self) -> None: 

350 """ 

351 Compute initial block quantities (sphere radius, surface, volume, ...) and solve 

352 the associated static problem to get initial displacement and fiber deformation. 

353 """ 

354 self.inv_radius = 1.0 / self.radius.current 

355 self.sphere_surface = 4.0 * np.pi * np.power(self.radius.current, 2) 

356 self.sphere_volume = 4.0 / 3.0 * np.pi * np.power(self.radius.current, 3) 

357 self.thickness_radius_ratio = self.thickness.current * self.inv_radius 

358 self.half_thickness_radius_ratio = 0.5 * self.thickness_radius_ratio 

359 self.viscous_coef = ( 

360 2.0 

361 * self.damping_coef.current 

362 * self.thickness.current 

363 * np.pow(self.inv_radius, 2) 

364 ) 

365 self.hyperelastic_cst_01 = ( 

366 self.hyperelastic_cst.current[0] * self.hyperelastic_cst.current[1] 

367 ) 

368 self.hyperelastic_cst_23 = ( 

369 self.hyperelastic_cst.current[2] * self.hyperelastic_cst.current[3] 

370 ) 

371 

372 disp_init, fib_deform_init = self._solve_static_problem() 

373 

374 self.disp.initialize(disp_init) 

375 self.fib_deform.initialize(fib_deform_init) 

376 

377 def dynamics_residual(self) -> Any: 

378 """ 

379 Compute the residual giving a dynamics on `disp` 

380 

381 :return: the residual 

382 :rtype: np.float64 

383 """ 

384 

385 disp_mid = mid_point(self.disp) 

386 fib_deform_mid = mid_point(self.fib_deform) 

387 p_mid_point = mid_point(self.pressure) 

388 pressure_external = mid_point(self.pressure_external) 

389 

390 disp_diff = diff(self.disp) 

391 

392 disp_mid_ratio = 1.0 + disp_mid * self.inv_radius 

393 disp_new_ratio = 1.0 + self.disp.new * self.inv_radius 

394 disp_cur_ratio = 1.0 + self.disp.current * self.inv_radius 

395 

396 j4_mid = np.pow(disp_mid_ratio, 2) 

397 j4_diff_mid = 2.0 * self.inv_radius * disp_mid_ratio 

398 

399 j1_mid = 2.0 * j4_mid + np.pow(disp_mid_ratio, -4) 

400 

401 # inertia stress 

402 inertia_stress = ( 

403 self.vol_mass.current 

404 * self.thickness.current 

405 * diff(self.vel) 

406 * self.time.inv_dt 

407 ) 

408 

409 # viscous stress 

410 viscous_diff_vel = ( 

411 self.viscous_coef * j4_mid * (1.0 + 2.0 * np.pow(disp_mid_ratio, -12)) 

412 ) 

413 viscous_stress = viscous_diff_vel * disp_diff * self.time.inv_dt 

414 

415 # active stress 

416 active_stress = ( 

417 self.thickness_radius_ratio 

418 * self.series_stiffness.current 

419 * (disp_mid * self.inv_radius - fib_deform_mid) 

420 ) 

421 if np.abs(disp_diff / self.disp.current) > _DISP_EPSILON: 

422 disp_diff_ratio = disp_diff * self.inv_radius 

423 

424 # external stress 

425 fluid_volume_cur = ( 

426 disp_cur_ratio 

427 - np.pow(disp_cur_ratio, -2) * self.half_thickness_radius_ratio 

428 ) 

429 fluid_volume_new = ( 

430 disp_new_ratio 

431 - np.pow(disp_new_ratio, -2) * self.half_thickness_radius_ratio 

432 ) 

433 fluid_volume_diff3 = ( 

434 disp_diff_ratio 

435 - math_utils.power_diff( 

436 disp_new_ratio, disp_cur_ratio, disp_diff_ratio, -2 

437 ) 

438 * self.half_thickness_radius_ratio 

439 ) 

440 fluid_volume_diff = self.sphere_volume * math_utils.power_diff( 

441 fluid_volume_new, fluid_volume_cur, fluid_volume_diff3, 3 

442 ) 

443 

444 external_stress = ( 

445 (p_mid_point - pressure_external) 

446 * fluid_volume_diff 

447 / (self.sphere_surface * disp_diff) 

448 ) 

449 

450 # passive stress 

451 j4_new = np.pow(disp_new_ratio, 2) 

452 j4_cur = np.pow(disp_cur_ratio, 2) 

453 

454 j4_diff = math_utils.power_diff( 

455 disp_new_ratio, disp_cur_ratio, disp_diff_ratio, 2 

456 ) 

457 

458 j1_new = 2.0 * j4_new + np.pow(disp_new_ratio, -4) 

459 j1_cur = 2.0 * j4_cur + np.pow(disp_cur_ratio, -4) 

460 j1_diff = 2.0 * j4_diff + math_utils.power_diff( 

461 disp_new_ratio, disp_cur_ratio, disp_diff_ratio, -4 

462 ) 

463 

464 arg_exp_j1_new = self.hyperelastic_cst.current[1] * np.pow(j1_new - 3.0, 2) 

465 arg_exp_j1_cur = self.hyperelastic_cst.current[1] * np.pow(j1_cur - 3.0, 2) 

466 arg_exp_j1_diff = self.hyperelastic_cst.current[1] * math_utils.power_diff( 

467 j1_new - 3.0, j1_cur - 3.0, j1_diff, 2 

468 ) 

469 

470 arg_exp_j4_new = self.hyperelastic_cst.current[3] * np.pow(j4_new - 1.0, 2) 

471 arg_exp_j4_cur = self.hyperelastic_cst.current[3] * np.pow(j4_cur - 1.0, 2) 

472 arg_exp_j4_diff = self.hyperelastic_cst.current[3] * math_utils.power_diff( 

473 j4_new - 1.0, j4_cur - 1.0, j4_diff, 2 

474 ) 

475 

476 exp_j1_diff = math_utils.exp_diff( 

477 arg_exp_j1_new, arg_exp_j1_cur, arg_exp_j1_diff 

478 ) 

479 exp_j4_diff = math_utils.exp_diff( 

480 arg_exp_j4_new, arg_exp_j4_cur, arg_exp_j4_diff 

481 ) 

482 

483 hyperelastic_potential_diff = ( 

484 self.hyperelastic_cst.current[0] * exp_j1_diff 

485 + self.hyperelastic_cst.current[2] * exp_j4_diff 

486 ) 

487 

488 passive_stress = ( 

489 self.thickness.current * hyperelastic_potential_diff / disp_diff 

490 ) 

491 

492 else: 

493 # external stress 

494 disp_mid_ratio_adj = np.pow( 

495 disp_mid_ratio 

496 - np.pow(disp_mid_ratio, -2) * self.half_thickness_radius_ratio, 

497 2, 

498 ) 

499 external_stress = ( 

500 (p_mid_point - pressure_external) 

501 * disp_mid_ratio_adj 

502 * (1.0 + np.pow(disp_mid_ratio, -3) * self.thickness_radius_ratio) 

503 ) 

504 

505 # passive stress 

506 exp_j1_mid = np.exp( 

507 self.hyperelastic_cst.current[1] * np.pow(j1_mid - 3.0, 2) 

508 ) 

509 exp_j4_mid = np.exp( 

510 self.hyperelastic_cst.current[3] * np.pow(j4_mid - 1.0, 2) 

511 ) 

512 

513 j1_diff_mid = 2.0 * j4_diff_mid - 4.0 * self.inv_radius * np.pow( 

514 disp_mid_ratio, -5 

515 ) 

516 

517 hyperelastic_potential_diff_mid = ( 

518 2.0 

519 * self.hyperelastic_cst_01 

520 * j1_diff_mid 

521 * (j1_mid - 3.0) 

522 * exp_j1_mid 

523 ) + ( 

524 2.0 

525 * self.hyperelastic_cst_23 

526 * j4_diff_mid 

527 * (j4_mid - 1.0) 

528 * exp_j4_mid 

529 ) 

530 

531 passive_stress = self.thickness.current * hyperelastic_potential_diff_mid 

532 

533 return ( 

534 inertia_stress 

535 + passive_stress 

536 + viscous_stress 

537 + active_stress 

538 - external_stress 

539 ) 

540 

541 def dynamics_residual_ddisp(self) -> Any: 

542 """ 

543 Compute the partial derivative of the residual for `disp`. 

544 

545 :return: the residual partial derivative for `disp` 

546 :rtype: np.float64 

547 """ 

548 

549 disp_mid = mid_point(self.disp) 

550 disp_diff = diff(self.disp) 

551 disp_mid_ratio = 1.0 + disp_mid * self.inv_radius 

552 disp_cur_ratio = 1.0 + self.disp.current * self.inv_radius 

553 disp_new_ratio = 1.0 + self.disp.new * self.inv_radius 

554 pressure_external = mid_point(self.pressure_external) 

555 

556 p_mid_point = mid_point(self.pressure) 

557 

558 j4_mid = np.pow(disp_mid_ratio, 2) 

559 j4_diff_mid = 2.0 * self.inv_radius * (disp_mid_ratio) 

560 j1_mid = 2.0 * j4_mid + np.pow(disp_mid_ratio, -4) 

561 

562 # viscous 

563 viscous_diff_vel = ( 

564 self.viscous_coef * j4_mid * (1.0 + 2.0 * np.pow(disp_mid_ratio, -12)) 

565 ) 

566 viscous_potential_diff_disp_mid = ( 

567 self.viscous_coef 

568 * disp_diff 

569 * self.time.inv_dt 

570 * ( 

571 j4_diff_mid * (1.0 + 2.0 * np.pow(disp_mid_ratio, -12)) 

572 - 24.0 * self.inv_radius * j4_mid * np.pow(disp_mid_ratio, -13) 

573 ) 

574 ) 

575 viscous_ddisp = ( 

576 0.5 * viscous_potential_diff_disp_mid + self.time.inv_dt * viscous_diff_vel 

577 ) 

578 

579 # active 

580 active_ddisp = ( 

581 self.half_thickness_radius_ratio 

582 * self.series_stiffness.current 

583 * self.inv_radius 

584 ) 

585 

586 if np.fabs(disp_diff / self.disp.current) > _DISP_EPSILON: 

587 # external 

588 disp_diff_ratio = disp_diff * self.inv_radius 

589 fluid_volume_term_cur = ( 

590 disp_cur_ratio 

591 - np.pow(disp_cur_ratio, -2) * self.half_thickness_radius_ratio 

592 ) 

593 fluid_volume_term_new = ( 

594 disp_new_ratio 

595 - np.pow(disp_new_ratio, -2) * self.half_thickness_radius_ratio 

596 ) 

597 fluid_volume_term_diff3 = ( 

598 disp_diff_ratio 

599 - math_utils.power_diff( 

600 disp_new_ratio, disp_cur_ratio, disp_diff_ratio, -2 

601 ) 

602 * self.half_thickness_radius_ratio 

603 ) 

604 fluid_volume_diff = self.sphere_volume * math_utils.power_diff( 

605 fluid_volume_term_new, fluid_volume_term_cur, fluid_volume_term_diff3, 3 

606 ) 

607 

608 fluid_volume_diff_new = ( 

609 self.sphere_surface 

610 * np.pow( 

611 disp_new_ratio 

612 - np.pow(disp_new_ratio, -2) * self.half_thickness_radius_ratio, 

613 2, 

614 ) 

615 * (1.0 + np.pow(disp_new_ratio, -3) * self.thickness_radius_ratio) 

616 ) 

617 

618 fluid_volume_diff_scheme_diff = ( 

619 (1.0 / self.sphere_surface) 

620 * (fluid_volume_diff_new * disp_diff - fluid_volume_diff) 

621 * np.pow(disp_diff, -2) 

622 ) 

623 

624 external_ddisp = ( 

625 p_mid_point - pressure_external 

626 ) * fluid_volume_diff_scheme_diff 

627 

628 # passive 

629 j4_cur = np.pow(disp_cur_ratio, 2) 

630 j4_new = np.pow(disp_new_ratio, 2) 

631 j4_diff_new = 2.0 * self.inv_radius * disp_new_ratio 

632 j4_diff = math_utils.power_diff( 

633 disp_new_ratio, disp_cur_ratio, disp_diff_ratio, 2 

634 ) 

635 

636 j1_new = 2.0 * j4_new + np.pow(disp_new_ratio, -4) 

637 j1_cur = 2.0 * j4_cur + np.pow(disp_cur_ratio, -4) 

638 j1_diff_new = 2.0 * j4_diff_new - 4.0 * self.inv_radius * np.pow( 

639 disp_new_ratio, -5 

640 ) 

641 j1_diff = 2.0 * j4_diff + math_utils.power_diff( 

642 disp_new_ratio, disp_cur_ratio, disp_diff_ratio, -4 

643 ) 

644 

645 arg_exp_j1_cur = self.hyperelastic_cst.current[1] * np.pow(j1_cur - 3.0, 2) 

646 arg_exp_j1_new = self.hyperelastic_cst.current[1] * np.pow(j1_new - 3.0, 2) 

647 arg_exp_j1_diff = self.hyperelastic_cst.current[1] * math_utils.power_diff( 

648 j1_new - 3.0, j1_cur - 3.0, j1_diff, 2 

649 ) 

650 exp_j1_new = np.exp(arg_exp_j1_new) 

651 

652 arg_exp_j4_cur = self.hyperelastic_cst.current[3] * np.pow(j4_cur - 1.0, 2) 

653 arg_exp_j4_new = self.hyperelastic_cst.current[3] * np.pow(j4_new - 1.0, 2) 

654 arg_exp_j4_diff = self.hyperelastic_cst.current[3] * math_utils.power_diff( 

655 j4_new - 1.0, j4_cur - 1.0, j4_diff, 2 

656 ) 

657 exp_j4_new = np.exp(arg_exp_j4_new) 

658 

659 exp_j4_diff = math_utils.exp_diff( 

660 arg_exp_j4_new, arg_exp_j4_cur, arg_exp_j4_diff 

661 ) 

662 exp_j1_diff = math_utils.exp_diff( 

663 arg_exp_j1_new, arg_exp_j1_cur, arg_exp_j1_diff 

664 ) 

665 

666 hyperelastic_potential_diff = ( 

667 self.hyperelastic_cst.current[0] * exp_j1_diff 

668 + self.hyperelastic_cst.current[2] * exp_j4_diff 

669 ) 

670 hyperelastic_potential_diff_new = ( 

671 2.0 

672 * self.hyperelastic_cst_01 

673 * j1_diff_new 

674 * (j1_new - 3.0) 

675 * exp_j1_new 

676 ) + ( 

677 2.0 

678 * self.hyperelastic_cst_23 

679 * j4_diff_new 

680 * (j4_new - 1.0) 

681 * exp_j4_new 

682 ) 

683 passive_ddisp = ( 

684 self.thickness.current 

685 * ( 

686 hyperelastic_potential_diff_new * disp_diff 

687 - hyperelastic_potential_diff 

688 ) 

689 / np.pow(disp_diff, 2) 

690 ) 

691 

692 else: 

693 # external 

694 disp_mid_ratio_adj = np.pow( 

695 disp_mid_ratio 

696 - np.pow(disp_mid_ratio, -2) * self.half_thickness_radius_ratio, 

697 2, 

698 ) 

699 fluid_volume_diff_diff_mid = ( 

700 4.0 

701 * np.pi 

702 * self.radius.current 

703 * ( 

704 2.0 

705 * ( 

706 disp_mid_ratio 

707 - np.pow(disp_mid_ratio, -2) * self.half_thickness_radius_ratio 

708 ) 

709 * np.pow( 

710 1.0 + np.pow(disp_mid_ratio, -3) * self.thickness_radius_ratio, 

711 2, 

712 ) 

713 - ( 

714 3.0 

715 * self.thickness_radius_ratio 

716 * np.pow(disp_mid_ratio, -4) 

717 * disp_mid_ratio_adj 

718 ) 

719 ) 

720 ) 

721 fluid_volume_diff_scheme_diff = ( 

722 1.0 / self.sphere_surface * 0.5 * fluid_volume_diff_diff_mid 

723 ) 

724 external_ddisp = ( 

725 p_mid_point - pressure_external 

726 ) * fluid_volume_diff_scheme_diff 

727 

728 # passive 

729 exp_j1_mid = np.exp( 

730 self.hyperelastic_cst.current[1] * np.pow(j1_mid - 3.0, 2) 

731 ) 

732 exp_j4_mid = np.exp( 

733 self.hyperelastic_cst.current[3] * np.pow(j4_mid - 1.0, 2) 

734 ) 

735 j1_diff_mid = 2.0 * j4_diff_mid - 4.0 * self.inv_radius * np.pow( 

736 disp_mid_ratio, -5 

737 ) 

738 j4_diff_diff = 2.0 * np.pow(self.inv_radius, 2) 

739 j1_diff_diff_mid = 4.0 * np.pow(self.inv_radius, 2) + 20.0 * np.pow( 

740 self.inv_radius, 2 

741 ) * np.pow(disp_mid_ratio, -6) 

742 hyperelastic_potential_diff_diff_mid = ( 

743 2.0 

744 * self.hyperelastic_cst_01 

745 * ( 

746 j1_diff_diff_mid * (j1_mid - 3.0) 

747 + np.pow(j1_diff_mid, 2) 

748 + 2.0 

749 * self.hyperelastic_cst.current[1] 

750 * np.pow(j1_diff_mid * (j1_mid - 3.0), 2) 

751 ) 

752 * exp_j1_mid 

753 + 2.0 

754 * self.hyperelastic_cst_23 

755 * ( 

756 j4_diff_diff * (j4_mid - 1.0) 

757 + np.pow(j4_diff_mid, 2) 

758 + 2.0 

759 * self.hyperelastic_cst.current[3] 

760 * np.pow(j4_diff_mid * (j4_mid - 1.0), 2) 

761 ) 

762 * exp_j4_mid 

763 ) 

764 passive_ddisp = ( 

765 self.thickness.current * 0.5 * hyperelastic_potential_diff_diff_mid 

766 ) 

767 

768 return viscous_ddisp + active_ddisp - external_ddisp + passive_ddisp 

769 

770 def dynamics_residual_dfib_deform(self) -> Any: 

771 """ 

772 

773 Compute the partial derivative of the residual for `fib_deform`. 

774 

775 :return: the residual partial derivative for `fib_deform` 

776 :rtype: np.float64 

777 """ 

778 # activ stress 

779 return -self.half_thickness_radius_ratio * self.series_stiffness.current 

780 

781 def dynamics_residual_dvel(self) -> Any: 

782 """ 

783 Compute the partial derivative of the residual for `vel`. 

784 

785 :return: the residual partial derivative for `vel` 

786 :rtype: np.float64 

787 """ 

788 # inertia 

789 return self.vol_mass.current * self.thickness.current * self.time.inv_dt 

790 

791 def dynamics_residual_dpressure(self) -> Any: 

792 """ 

793 Compute the partial derivative of the residual for the `pressure`. 

794 

795 :return: the residual partial derivative for `pressure` 

796 :rtype: np.float64 

797 """ 

798 

799 # external 

800 disp_diff = diff(self.disp) 

801 disp_mid = mid_point(self.disp) 

802 disp_mid_ratio = 1.0 + disp_mid * self.inv_radius 

803 disp_new_ratio = 1.0 + self.disp.new * self.inv_radius 

804 disp_cur_ratio = 1.0 + self.disp.current * self.inv_radius 

805 

806 if np.fabs(disp_diff / self.disp.current) > _DISP_EPSILON: 

807 disp_diff_ratio = disp_diff * self.inv_radius 

808 fluid_volume_term_cur = ( 

809 disp_cur_ratio 

810 - np.pow(disp_cur_ratio, -2) * self.half_thickness_radius_ratio 

811 ) 

812 fluid_volume_term_new = ( 

813 disp_new_ratio 

814 - np.pow(disp_new_ratio, -2) * self.half_thickness_radius_ratio 

815 ) 

816 

817 fluid_volume_term_diff3 = ( 

818 disp_diff_ratio 

819 - math_utils.power_diff( 

820 disp_new_ratio, disp_cur_ratio, disp_diff_ratio, -2 

821 ) 

822 * self.half_thickness_radius_ratio 

823 ) 

824 

825 fluid_volume_diff = self.sphere_volume * math_utils.power_diff( 

826 fluid_volume_term_new, fluid_volume_term_cur, fluid_volume_term_diff3, 3 

827 ) 

828 external_dpressure = -( 

829 (0.5 * fluid_volume_diff) / (self.sphere_surface * disp_diff) 

830 ) 

831 else: 

832 disp_mid_ratio_adj = np.pow( 

833 disp_mid_ratio 

834 - np.pow(disp_mid_ratio, -2) * self.half_thickness_radius_ratio, 

835 2, 

836 ) 

837 external_dpressure = ( 

838 -0.5 

839 * disp_mid_ratio_adj 

840 * (1.0 + np.pow(disp_mid_ratio, -3) * self.thickness_radius_ratio) 

841 ) 

842 

843 return external_dpressure 

844 

845 def dynamics_residual_dpressure_external(self) -> Any: 

846 """ 

847 Compute the partial derivative of the residual `pressure_external`. 

848 

849 :return: the residual partial derivative `pressure_external`. 

850 :rtype: np.float64 

851 """ 

852 

853 disp_diff = diff(self.disp) 

854 

855 if np.abs(disp_diff / self.disp.current) > _DISP_EPSILON: 

856 disp_diff_ratio = disp_diff * self.inv_radius