Note
Go to the end to download the full example code.
MonoVentricular#
Passive + active hyperelastic simulation of an ellipsoidal left-ventricle model.
Combines a Holzapfel-Ogden orthotropic law (fiber + sheet directions), an active stress along the fiber direction, and a follower pressure on the endocardial surface. Time integration uses the midpoint hyperbolic scheme.
Reproduces Benchmark 1: monoventricular mechanics (§3) of the cardiac elastodynamics benchmark published in Comput. Methods Appl. Mech. Engrg.: https://www.sciencedirect.com/science/article/pii/S0045782524007394
The mesh.msh / fiber.vtu / sheet.vtu files read for fiberSource="vtu" are generated beforehand with the cardiac_benchmark_toolkit — see the module docstring of utils.py for the exact procedure. (fiberSource="analytic" builds the fibers/sheets directly in EasyFEA and needs no external data.)
With useCoarseConfig=False this becomes a large 3D, non-linear, transient problem for which the default direct solver is slow. In that case it is recommended to run it either in parallel with MPI and PETSc (e.g. mpiexec -n <N> python MonoVentricular.py with a PETSc-backed solver), or, on a single process, with the pypardiso solver — both markedly cut the solve time. The default useCoarseConfig=True is light enough to run as-is.
===== hyperelastic problem at iteration 0 =====
At Newton iteration 1 norm is 8.686270292626e-02
At Newton iteration 2 norm is 2.083382951295e-04
At Newton iteration 3 norm is 2.195014942968e-09
===== hyperelastic problem at iteration 1 =====
At Newton iteration 1 norm is 8.636333907881e-02
At Newton iteration 2 norm is 8.112107297989e-04
At Newton iteration 3 norm is 3.260032262147e-08
===== hyperelastic problem at iteration 2 =====
At Newton iteration 1 norm is 6.667722093496e-02
At Newton iteration 2 norm is 6.598459930140e-04
At Newton iteration 3 norm is 1.979971164883e-08
===== hyperelastic problem at iteration 3 =====
At Newton iteration 1 norm is 5.443853101853e-02
At Newton iteration 2 norm is 3.574030065226e-04
At Newton iteration 3 norm is 5.100376790617e-09
===== hyperelastic problem at iteration 4 =====
At Newton iteration 1 norm is 3.568161669214e-02
At Newton iteration 2 norm is 1.879472302928e-04
At Newton iteration 3 norm is 7.634289365495e-10
===== hyperelastic problem at iteration 5 =====
At Newton iteration 1 norm is 3.345855815028e-02
At Newton iteration 2 norm is 1.104480379896e-04
At Newton iteration 3 norm is 1.727627710512e-10
===== hyperelastic problem at iteration 6 =====
At Newton iteration 1 norm is 2.145570404913e-02
At Newton iteration 2 norm is 5.743376127188e-05
At Newton iteration 3 norm is 4.241768320603e-11
===== hyperelastic problem at iteration 7 =====
At Newton iteration 1 norm is 2.002055496981e-02
At Newton iteration 2 norm is 2.858175612514e-05
At Newton iteration 3 norm is 9.345753106690e-12
===== hyperelastic problem at iteration 8 =====
At Newton iteration 1 norm is 1.215316851402e-02
At Newton iteration 2 norm is 1.697272617536e-05
At Newton iteration 3 norm is 3.527291499078e-12
===== hyperelastic problem at iteration 9 =====
At Newton iteration 1 norm is 1.284612269410e-02
At Newton iteration 2 norm is 1.004326528933e-05
At Newton iteration 3 norm is 1.235996777177e-12
==== hyperelastic problem at iteration 10 ====
At Newton iteration 1 norm is 7.358674634754e-03
At Newton iteration 2 norm is 4.663340737393e-06
==== hyperelastic problem at iteration 11 ====
At Newton iteration 1 norm is 7.987532659073e-03
At Newton iteration 2 norm is 3.232423375926e-06
==== hyperelastic problem at iteration 12 ====
At Newton iteration 1 norm is 4.809825935390e-03
At Newton iteration 2 norm is 1.536260526938e-06
==== hyperelastic problem at iteration 13 ====
At Newton iteration 1 norm is 1.039439603742e-02
At Newton iteration 2 norm is 3.486662015805e-06
==== hyperelastic problem at iteration 14 ====
At Newton iteration 1 norm is 1.766079917895e+00
At Newton iteration 2 norm is 9.125789691989e-02
At Newton iteration 3 norm is 1.087619504319e-03
At Newton iteration 4 norm is 2.420676246989e-08
==== hyperelastic problem at iteration 15 ====
At Newton iteration 1 norm is 2.001628154789e+00
At Newton iteration 2 norm is 3.964226109087e-01
At Newton iteration 3 norm is 2.637389813702e-04
At Newton iteration 4 norm is 1.322534260539e-09
==== hyperelastic problem at iteration 16 ====
At Newton iteration 1 norm is 1.719143781265e+00
At Newton iteration 2 norm is 3.830356713260e-01
At Newton iteration 3 norm is 1.109056791564e-04
At Newton iteration 4 norm is 1.090225288152e-09
==== hyperelastic problem at iteration 17 ====
At Newton iteration 1 norm is 1.641029957523e+00
At Newton iteration 2 norm is 2.661771005199e-01
At Newton iteration 3 norm is 8.687518365251e-05
At Newton iteration 4 norm is 7.305458876243e-10
==== hyperelastic problem at iteration 18 ====
At Newton iteration 1 norm is 1.494367369726e+00
At Newton iteration 2 norm is 2.463078554385e-01
At Newton iteration 3 norm is 5.438425389559e-05
At Newton iteration 4 norm is 1.614534352445e-10
==== hyperelastic problem at iteration 19 ====
At Newton iteration 1 norm is 1.528855945657e+00
At Newton iteration 2 norm is 2.503221856979e-01
At Newton iteration 3 norm is 6.025836143964e-05
At Newton iteration 4 norm is 1.927372767715e-10
==== hyperelastic problem at iteration 20 ====
At Newton iteration 1 norm is 1.375005604581e+00
At Newton iteration 2 norm is 2.018573037972e-01
At Newton iteration 3 norm is 4.371244263022e-05
At Newton iteration 4 norm is 8.944786097194e-11
==== hyperelastic problem at iteration 21 ====
At Newton iteration 1 norm is 1.394964362988e+00
At Newton iteration 2 norm is 1.714802611816e-01
At Newton iteration 3 norm is 4.447758721415e-05
At Newton iteration 4 norm is 1.276317202623e-10
==== hyperelastic problem at iteration 22 ====
At Newton iteration 1 norm is 1.263965377752e+00
At Newton iteration 2 norm is 1.478516823517e-01
At Newton iteration 3 norm is 3.463498057017e-05
At Newton iteration 4 norm is 9.065879147737e-11
==== hyperelastic problem at iteration 23 ====
At Newton iteration 1 norm is 1.285731015645e+00
At Newton iteration 2 norm is 1.217372426517e-01
At Newton iteration 3 norm is 3.773930959185e-05
At Newton iteration 4 norm is 1.201147139048e-10
==== hyperelastic problem at iteration 24 ====
At Newton iteration 1 norm is 1.138920590138e+00
At Newton iteration 2 norm is 9.680485073069e-02
At Newton iteration 3 norm is 2.354890714584e-05
At Newton iteration 4 norm is 5.133215014504e-11
==== hyperelastic problem at iteration 25 ====
At Newton iteration 1 norm is 1.150907660345e+00
At Newton iteration 2 norm is 8.160509659664e-02
At Newton iteration 3 norm is 2.098969826639e-05
At Newton iteration 4 norm is 3.969201599009e-11
==== hyperelastic problem at iteration 26 ====
At Newton iteration 1 norm is 1.026993249552e+00
At Newton iteration 2 norm is 6.433957030611e-02
At Newton iteration 3 norm is 1.221062638004e-05
At Newton iteration 4 norm is 1.355694167085e-11
==== hyperelastic problem at iteration 27 ====
At Newton iteration 1 norm is 1.016718457051e+00
At Newton iteration 2 norm is 5.201368367350e-02
At Newton iteration 3 norm is 1.036869328924e-05
At Newton iteration 4 norm is 9.053735305167e-12
==== hyperelastic problem at iteration 28 ====
At Newton iteration 1 norm is 9.079659358231e-01
At Newton iteration 2 norm is 4.111574615345e-02
At Newton iteration 3 norm is 6.121538205175e-06
==== hyperelastic problem at iteration 29 ====
At Newton iteration 1 norm is 8.911742704472e-01
At Newton iteration 2 norm is 3.252948938897e-02
At Newton iteration 3 norm is 5.246258504945e-06
==== hyperelastic problem at iteration 30 ====
At Newton iteration 1 norm is 7.949198299427e-01
At Newton iteration 2 norm is 2.573282611816e-02
At Newton iteration 3 norm is 3.138421962920e-06
==== hyperelastic problem at iteration 31 ====
At Newton iteration 1 norm is 7.774807157864e-01
At Newton iteration 2 norm is 2.033244512179e-02
At Newton iteration 3 norm is 2.720349999965e-06
==== hyperelastic problem at iteration 32 ====
At Newton iteration 1 norm is 6.966625324299e-01
At Newton iteration 2 norm is 1.638331490727e-02
At Newton iteration 3 norm is 1.672881898897e-06
==== hyperelastic problem at iteration 33 ====
At Newton iteration 1 norm is 6.777500156534e-01
At Newton iteration 2 norm is 1.304078137846e-02
At Newton iteration 3 norm is 1.444718014665e-06
==== hyperelastic problem at iteration 34 ====
At Newton iteration 1 norm is 6.130620615771e-01
At Newton iteration 2 norm is 1.077692758011e-02
At Newton iteration 3 norm is 9.146600842569e-07
==== hyperelastic problem at iteration 35 ====
At Newton iteration 1 norm is 5.902687207387e-01
At Newton iteration 2 norm is 8.705891529224e-03
At Newton iteration 3 norm is 7.830888659807e-07
==== hyperelastic problem at iteration 36 ====
At Newton iteration 1 norm is 5.422373569373e-01
At Newton iteration 2 norm is 7.351341456356e-03
At Newton iteration 3 norm is 5.100601098657e-07
==== hyperelastic problem at iteration 37 ====
At Newton iteration 1 norm is 5.124490985705e-01
At Newton iteration 2 norm is 6.024003824795e-03
At Newton iteration 3 norm is 4.294865641252e-07
==== hyperelastic problem at iteration 38 ====
At Newton iteration 1 norm is 4.155515550565e-01
At Newton iteration 2 norm is 4.091144051898e-03
At Newton iteration 3 norm is 1.715786348287e-07
==== hyperelastic problem at iteration 39 ====
At Newton iteration 1 norm is 7.223803002004e+00
At Newton iteration 2 norm is 3.374928378036e-01
At Newton iteration 3 norm is 1.512215943606e-03
At Newton iteration 4 norm is 1.209413876421e-07
==== hyperelastic problem at iteration 40 ====
At Newton iteration 1 norm is 6.790874276672e+00
At Newton iteration 2 norm is 1.262651254492e+00
At Newton iteration 3 norm is 1.188459381142e-02
At Newton iteration 4 norm is 9.884492848622e-06
==== hyperelastic problem at iteration 41 ====
At Newton iteration 1 norm is 4.406861200099e+00
At Newton iteration 2 norm is 1.225542894193e+00
At Newton iteration 3 norm is 4.626756732993e-03
At Newton iteration 4 norm is 1.977426693459e-06
==== hyperelastic problem at iteration 42 ====
At Newton iteration 1 norm is 3.378202669326e+00
At Newton iteration 2 norm is 6.716046814598e-01
At Newton iteration 3 norm is 1.294449634204e-03
At Newton iteration 4 norm is 1.571830950457e-07
==== hyperelastic problem at iteration 43 ====
At Newton iteration 1 norm is 1.993432017108e+00
At Newton iteration 2 norm is 4.776884593473e-01
At Newton iteration 3 norm is 3.286824519224e-04
At Newton iteration 4 norm is 8.291959281829e-09
==== hyperelastic problem at iteration 44 ====
At Newton iteration 1 norm is 1.944264247387e+00
At Newton iteration 2 norm is 2.955186274079e-01
At Newton iteration 3 norm is 2.113008822716e-04
At Newton iteration 4 norm is 3.912655961746e-09
==== hyperelastic problem at iteration 45 ====
At Newton iteration 1 norm is 1.334126825250e+00
At Newton iteration 2 norm is 1.871344231619e-01
At Newton iteration 3 norm is 7.617148573999e-05
At Newton iteration 4 norm is 4.128247384862e-10
==== hyperelastic problem at iteration 46 ====
At Newton iteration 1 norm is 1.327886215182e+00
At Newton iteration 2 norm is 1.041889384726e-01
At Newton iteration 3 norm is 4.908254976213e-05
At Newton iteration 4 norm is 2.065514114021e-10
==== hyperelastic problem at iteration 47 ====
At Newton iteration 1 norm is 1.105992962329e+00
At Newton iteration 2 norm is 7.189486787911e-02
At Newton iteration 3 norm is 1.594932229102e-05
At Newton iteration 4 norm is 8.968784050886e-12
==== hyperelastic problem at iteration 48 ====
At Newton iteration 1 norm is 1.045556999382e+00
At Newton iteration 2 norm is 3.645902679632e-02
At Newton iteration 3 norm is 1.042709877240e-05
At Newton iteration 4 norm is 5.680974521649e-12
==== hyperelastic problem at iteration 49 ====
At Newton iteration 1 norm is 9.805659860478e-01
At Newton iteration 2 norm is 2.521242433350e-02
At Newton iteration 3 norm is 3.338387370888e-06
==== hyperelastic problem at iteration 50 ====
At Newton iteration 1 norm is 9.015430803163e-01
At Newton iteration 2 norm is 1.236194514669e-02
At Newton iteration 3 norm is 2.147490501810e-06
==== hyperelastic problem at iteration 51 ====
At Newton iteration 1 norm is 8.888258892398e-01
At Newton iteration 2 norm is 9.729822282684e-03
At Newton iteration 3 norm is 7.514945046735e-07
==== hyperelastic problem at iteration 52 ====
At Newton iteration 1 norm is 8.174872022537e-01
At Newton iteration 2 norm is 4.528749982887e-03
At Newton iteration 3 norm is 5.045362407583e-07
==== hyperelastic problem at iteration 53 ====
At Newton iteration 1 norm is 8.132221796677e-01
At Newton iteration 2 norm is 3.991324564933e-03
At Newton iteration 3 norm is 2.558138833736e-07
==== hyperelastic problem at iteration 54 ====
At Newton iteration 1 norm is 7.548929810225e-01
At Newton iteration 2 norm is 1.860699119020e-03
At Newton iteration 3 norm is 8.104726244957e-08
==== hyperelastic problem at iteration 55 ====
At Newton iteration 1 norm is 7.481044324088e-01
At Newton iteration 2 norm is 1.727590418354e-03
At Newton iteration 3 norm is 5.737151473020e-08
==== hyperelastic problem at iteration 56 ====
At Newton iteration 1 norm is 7.008317160020e-01
At Newton iteration 2 norm is 9.127256507902e-04
At Newton iteration 3 norm is 6.812613363506e-08
==== hyperelastic problem at iteration 57 ====
At Newton iteration 1 norm is 6.880730859717e-01
At Newton iteration 2 norm is 7.777542153835e-04
At Newton iteration 3 norm is 2.432568018252e-08
==== hyperelastic problem at iteration 58 ====
At Newton iteration 1 norm is 6.497186901320e-01
At Newton iteration 2 norm is 4.614777707713e-04
At Newton iteration 3 norm is 2.275538484017e-08
==== hyperelastic problem at iteration 59 ====
At Newton iteration 1 norm is 6.317440138225e-01
At Newton iteration 2 norm is 3.629149331827e-04
At Newton iteration 3 norm is 4.183972630005e-09
==== hyperelastic problem at iteration 60 ====
At Newton iteration 1 norm is 5.995733407556e-01
At Newton iteration 2 norm is 2.346991735432e-04
At Newton iteration 3 norm is 3.682116152617e-09
==== hyperelastic problem at iteration 61 ====
At Newton iteration 1 norm is 5.785810082564e-01
At Newton iteration 2 norm is 1.973183190447e-04
At Newton iteration 3 norm is 1.075630130725e-09
==== hyperelastic problem at iteration 62 ====
At Newton iteration 1 norm is 5.500050438403e-01
At Newton iteration 2 norm is 1.313525537067e-04
At Newton iteration 3 norm is 8.066794410783e-10
==== hyperelastic problem at iteration 63 ====
At Newton iteration 1 norm is 5.275749142850e-01
At Newton iteration 2 norm is 1.234968808907e-04
At Newton iteration 3 norm is 5.010408362634e-10
==== hyperelastic problem at iteration 64 ====
At Newton iteration 1 norm is 5.011748613698e-01
At Newton iteration 2 norm is 8.399457911006e-05
At Newton iteration 3 norm is 2.632432200953e-10
==== hyperelastic problem at iteration 65 ====
At Newton iteration 1 norm is 4.782874190071e-01
At Newton iteration 2 norm is 8.538851797859e-05
At Newton iteration 3 norm is 3.156883776320e-10
==== hyperelastic problem at iteration 66 ====
At Newton iteration 1 norm is 4.534428032312e-01
At Newton iteration 2 norm is 6.310035833323e-05
At Newton iteration 3 norm is 1.597833567116e-10
==== hyperelastic problem at iteration 67 ====
At Newton iteration 1 norm is 4.307005851755e-01
At Newton iteration 2 norm is 6.570998134897e-05
At Newton iteration 3 norm is 2.207460118586e-10
==== hyperelastic problem at iteration 68 ====
At Newton iteration 1 norm is 4.071546078104e-01
At Newton iteration 2 norm is 5.213994257418e-05
At Newton iteration 3 norm is 1.363626068607e-10
==== hyperelastic problem at iteration 69 ====
At Newton iteration 1 norm is 3.849569253037e-01
At Newton iteration 2 norm is 5.431111795198e-05
At Newton iteration 3 norm is 1.550038396825e-10
==== hyperelastic problem at iteration 70 ====
At Newton iteration 1 norm is 3.626294460455e-01
At Newton iteration 2 norm is 4.555747807564e-05
At Newton iteration 3 norm is 1.070495159118e-10
==== hyperelastic problem at iteration 71 ====
At Newton iteration 1 norm is 3.412633887489e-01
At Newton iteration 2 norm is 4.639574949926e-05
At Newton iteration 3 norm is 1.029586294596e-10
==== hyperelastic problem at iteration 72 ====
At Newton iteration 1 norm is 3.201721472454e-01
At Newton iteration 2 norm is 4.053811618939e-05
At Newton iteration 3 norm is 7.351764730501e-11
==== hyperelastic problem at iteration 73 ====
At Newton iteration 1 norm is 2.998699388742e-01
At Newton iteration 2 norm is 3.989850739037e-05
At Newton iteration 3 norm is 6.345489510249e-11
==== hyperelastic problem at iteration 74 ====
At Newton iteration 1 norm is 2.801031257125e-01
At Newton iteration 2 norm is 3.548282532361e-05
At Newton iteration 3 norm is 4.524534475363e-11
==== hyperelastic problem at iteration 75 ====
At Newton iteration 1 norm is 2.610859087224e-01
At Newton iteration 2 norm is 3.388145462972e-05
At Newton iteration 3 norm is 3.594736322019e-11
==== hyperelastic problem at iteration 76 ====
At Newton iteration 1 norm is 2.427962193749e-01
At Newton iteration 2 norm is 3.024139515013e-05
At Newton iteration 3 norm is 2.509612755194e-11
==== hyperelastic problem at iteration 77 ====
At Newton iteration 1 norm is 2.253105891342e-01
At Newton iteration 2 norm is 2.819382665239e-05
At Newton iteration 3 norm is 1.863969815643e-11
==== hyperelastic problem at iteration 78 ====
At Newton iteration 1 norm is 2.087297394796e-01
At Newton iteration 2 norm is 2.519949220164e-05
At Newton iteration 3 norm is 1.285121661837e-11
==== hyperelastic problem at iteration 79 ====
At Newton iteration 1 norm is 1.930908341348e-01
At Newton iteration 2 norm is 2.318537658185e-05
At Newton iteration 3 norm is 9.618019872682e-12
==== hyperelastic problem at iteration 80 ====
At Newton iteration 1 norm is 1.785677476812e-01
At Newton iteration 2 norm is 2.083534784179e-05
At Newton iteration 3 norm is 7.439353259071e-12
Saved mesh data in:
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v2.0.0/examples/CardiacElastoDynamics/results/step1/ellipsoid0.03_dt0.0125_analytic_rigi/Meshes
Saved simulation and summary in:
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v2.0.0/examples/CardiacElastoDynamics/results/step1/ellipsoid0.03_dt0.0125_analytic_rigi
Generate movie 01/20 (5.00 %) 2.84 s
Generate movie 02/20 (10.00 %) 2.35 s
Generate movie 03/20 (15.00 %) 2.20 s
Generate movie 04/20 (20.00 %) 2.05 s
Generate movie 05/20 (25.00 %) 1.92 s
Generate movie 06/20 (30.00 %) 1.75 s
Generate movie 07/20 (35.00 %) 1.60 s
Generate movie 08/20 (40.00 %) 1.47 s
Generate movie 09/20 (45.00 %) 1.36 s
Generate movie 10/20 (50.00 %) 1.23 s
Generate movie 11/20 (55.00 %) 1.14 s
Generate movie 12/20 (60.00 %) 1.02 s
Generate movie 13/20 (65.00 %) 909.85 ms
Generate movie 14/20 (70.00 %) 777.42 ms
Generate movie 15/20 (75.00 %) 639.16 ms
Generate movie 16/20 (80.00 %) 516.13 ms
Generate movie 17/20 (85.00 %) 388.85 ms
Generate movie 18/20 (90.00 %) 260.32 ms
Generate movie 19/20 (95.00 %) 128.72 ms
Generate movie 20/20 (100.00 %) 0.00 µs
21 from enum import Enum
22
23 import numpy as np
24
25 from EasyFEA import Display, Folder, PyVista, MatrixType, Models, Simulations, AlgoType
26 from EasyFEA.FEM import Operators
27 from EasyFEA.Utilities import _params
28
29 from utils import RESULTS_DIR, DATA_DIR, Get_config, Get_values
30
31
32 class CardiacElastoDynamics(Simulations.HyperElastic):
33
34 pressure = _params.ScalarParameter()
35
36 def __init__(
37 self, mesh, model, folder="", tolConv=0.00001, maxIter=20, verbosity=False
38 ):
39 super().__init__(mesh, model, folder, tolConv, maxIter, verbosity)
40 self.pressure = 0.0
41
42 def Construct_local_matrix_system(self, problemType):
43
44 assert isinstance(self.material, Models.HyperElastic.HolzapfelOgden)
45 nPg = self.material.T1.shape[1]
46
47 results = super().Construct_local_matrix_system(problemType, nPg)
48
49 # current Newton-Raphson iterate (updated via u += delta_u)
50 displacement = self._Solver_Get_Newton_Raphson_current_solution()
51 if self.algo in AlgoType.Get_Hyperbolic_Types():
52 displacement, _, _ = self._Solver_Evaluate_u_v_a_for_time_scheme(
53 problemType, displacement
54 )
55
56 for groupElem in self.mesh.Get_list_groupElem(self.dim - 1):
57
58 # Following pressure (operator returns slot-ready values)
59 tangent_e, residual_e = Operators.NonLinear.FollowingPressure(
60 groupElem,
61 displacement,
62 self.pressure,
63 groupElem.Get_Elements_Tag("endo"),
64 MatrixType.mass,
65 )
66
67 # top — isotropic surface mass penalty (Robin α·u + β·u̇ = 0)
68 M_e = Operators.Bilinear.UV(groupElem, dof_n=3)
69
70 Ktop_e = np.zeros_like(M_e)
71 Ctop_e = np.zeros_like(M_e)
72 if "top" in groupElem.elementTags:
73 top_e = groupElem.Get_Elements_Tag("top")
74 alpha_top = 1e5
75 Ktop_e[top_e] = alpha_top * M_e[top_e]
76
77 beta_top = 5e3
78 Ctop_e[top_e] = beta_top * M_e[top_e]
79
80 # epi — normal-direction mass penalty (Robin α·(u·n̂) + β·(u̇·n̂) = 0)
81 Ms_e = Operators.Bilinear.MassAlongNormal(groupElem)
82
83 Kepi_e = np.zeros_like(Ms_e)
84 Cepi_e = np.zeros_like(Ms_e)
85 if "epi" in groupElem.elementTags:
86 epi_e = groupElem.Get_Elements_Tag("epi")
87 alpha_epi = 1e8
88 Kepi_e[epi_e] = alpha_epi * Ms_e[epi_e]
89
90 beta_epi = 5e3
91 Cepi_e[epi_e] = beta_epi * Ms_e[epi_e]
92
93 # Penalty residual contribution: −K_penalty · u_t at current iterate
94 K_penalty_e = Ktop_e + Kepi_e
95 assembly_e = groupElem.Get_assembly_e(self.dim)
96 u_e = displacement[assembly_e] # (Ne_surf, nPe·3)
97 f_penalty_e = np.einsum("eij,ej->ei", K_penalty_e, u_e)
98
99 results[groupElem] = (
100 tangent_e + K_penalty_e,
101 Ctop_e + Cepi_e,
102 None,
103 residual_e - f_penalty_e,
104 )
105
106 return results
107
108
109 class Config(str, Enum):
110 step0A = "step0A" # active_stress
111 step0B = "step0B" # pressure
112 step1 = "step1" # active_stress + pressure
113
114
115 if __name__ == "__main__":
116
117 Display.Clear()
118
119 # ----------------------------------------------
120 # Config
121 # ----------------------------------------------
122
123 useCoarseConfig = True
124
125 ellipsoid = "ellipsoid0.03" if useCoarseConfig else "ellipsoid0.005"
126
127 config = Config.step1
128
129 fiberSource = "analytic"
130 # fiberSource = "vtu"
131
132 matrixType = MatrixType.rigi
133 # matrixType = MatrixType.mass
134 # matrixType = 15
135
136 results_dir = Folder.Join(RESULTS_DIR, config.name, ellipsoid)
137
138 doSimu = True
139 plotGraph = False
140 plotParticles = True
141 saveParticles = True
142 makeMovie = True
143
144 # ----------------------------------------------
145 # time-history needed for plotting in both doSimu / Load_Simu flows
146
147 Nt = 80 if useCoarseConfig else 1000
148
149 t_values, activeStress_values, pressure_values = Get_values(Tmax=1.0, Nt=Nt)
150 dt = t_values[1] - t_values[0]
151 results_dir += f"_dt{dt}_{fiberSource}_{matrixType}"
152
153 if plotGraph:
154 ax = Display.Init_Axes()
155 ax.grid()
156 ax.set_xlabel(r"$t$ [s]")
157 ax.set_ylabel(r"$\tau(t)$ [Pa]")
158 ax.plot(t_values, activeStress_values)
159 name = "active_pressure"
160 Display.Save_fig(results_dir, name)
161
162 ax = Display.Init_Axes()
163 ax.grid()
164 ax.set_xlabel(r"$t$ [s]")
165 ax.set_ylabel(r"$p(t)$ [Pa]")
166 ax.plot(t_values, pressure_values)
167 name = "pressure"
168 Display.Save_fig(results_dir, name)
169
170 if config is Config.step0B:
171 activeStress_values *= 0
172 if config is Config.step0A:
173 pressure_values *= 0
174
175 if doSimu:
176
177 # ----------------------------------------------
178 # Mesh, fibers and sheets
179 # ----------------------------------------------
180
181 mesh, fibers_e_pg, sheets_e_pg = Get_config(
182 Folder.Join(DATA_DIR, ellipsoid),
183 matrixType=matrixType,
184 fiberSource=fiberSource,
185 plotMesh=False,
186 plotTags=False,
187 plotFibers=False,
188 )
189
190 # ----------------------------------------------
191 # Material
192 # ----------------------------------------------
193
194 # solid
195 a = 59.0
196 a_f = 18472.0
197 a_fs = 216.0
198 a_s = 2481.0
199 b = 8.023
200 b_f = 16.026
201 b_fs = 11.436
202 b_s = 11.12
203
204 material = Models.HyperElastic.HolzapfelOgden(
205 dim=3,
206 C0=a / 2 / b,
207 C1=b,
208 C2=a_f / 2 / b_f,
209 C3=b_f,
210 C4=a_s / 2 / b_s,
211 C5=b_s,
212 C6=a_fs / 2 / b_fs,
213 C7=b_fs,
214 K=1e6,
215 Mu1=0.0,
216 Mu2=0.0,
217 T1=fibers_e_pg,
218 T2=sheets_e_pg,
219 ks=100,
220 )
221 material.eta = 100.0
222 material.Set_active_stress_vec(material.T1)
223
224 # ----------------------------------------------
225 # Simulation
226 # ----------------------------------------------
227
228 simu = CardiacElastoDynamics(mesh, material, folder=results_dir)
229
230 simu.Solver_Set_Hyperbolic_Algorithm(dt, algo=AlgoType.midpoint)
231 simu.rho = 1000
232
233 for t in t_values:
234 simu.Bc_Init()
235 simu.pressure = np.interp(t + dt / 2, t_values, pressure_values)
236 material.active_stress = np.interp(
237 t + dt / 2, t_values, activeStress_values
238 )
239 simu.Solve()
240 simu.Save_Iter()
241
242 simu.Save(results_dir)
243
244 else:
245 simu = Simulations.Load_Simu(results_dir)
246
247 simu._Gather()
248
249 if plotParticles and simu.isGathered:
250
251 coords = [(0.025, 0.03, 0), (0, 0.03, 0)]
252 evalCoords = np.array(coords)
253 evalElements = simu.mesh.groupElem._Get_nearby_elements(evalCoords)
254
255 Niter = simu.Niter
256 values = np.empty((Niter, len(coords), 3))
257 for i in range(Niter):
258 simu.Set_Iter(i)
259 values[i] = simu.mesh.Evaluate_dofsValues_at_coordinates(
260 evalCoords, simu.displacement, evalElements
261 )
262
263 times = t_values[:Niter]
264 axs = Display.plt.subplots(3, 2, sharex=True)[1]
265
266 for p, (particle, coord) in enumerate(zip(["p0", "p1"], coords)):
267
268 for c, component in enumerate(["x", "y", "z"]):
269
270 ax: Display.plt.Axes = axs[c, p]
271
272 ax.grid()
273
274 if c == 2:
275 ax.set_xlabel("Time [s]")
276 if p == 0:
277 ax.set_ylabel(f"Displacement {component}-component [m]")
278 if c == 0:
279 ax.set_title(f"Particle {particle}")
280
281 ax.plot(times, values[:, p, c])
282
283 width, height = ax.figure.get_size_inches()
284 ax.figure.set_size_inches(width * 1.5, height * 2.5)
285 Display.Save_fig(results_dir, "particles")
286
287 if saveParticles and simu.isGathered:
288
289 # per-iteration deformed volume
290 volumes = np.empty(Niter)
291 for i in range(Niter):
292 simu.Set_Iter(i)
293 deformed = simu.mesh.copy()
294 deformed.coord += simu.displacement.reshape(-1, 3)
295 volumes[i] = deformed.volume
296
297 dict_particles = {
298 "time": times,
299 "displacement": {
300 f"p{p}": {
301 "ux": values[:, p, 0],
302 "uy": values[:, p, 1],
303 "uz": values[:, p, 2],
304 "magnitude": np.linalg.norm(values[:, p, :], axis=1),
305 }
306 for p in range(2)
307 },
308 "stress": {
309 "time": None,
310 "p0": {"magnitude": None},
311 "p1": {"magnitude": None},
312 },
313 "volume": volumes,
314 }
315 Simulations.Save_pickle(dict_particles, results_dir, "particles")
316
317 if makeMovie:
318 values = [simu.Result("ux", iter=i) for i in range(simu.Niter)]
319 clim = (np.min(values), np.max(values))
320 PyVista.Movie_simu(
321 simu,
322 "ux",
323 results_dir,
324 "ux.gif",
325 N=20,
326 deformFactor=1.0,
327 clim=clim,
328 plotMesh=True,
329 )
330
331 Display.plt.show()
Total running time of the script: (0 minutes 13.501 seconds)