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.026777061849e-04
At Newton iteration 3 norm is 2.476805657965e-09
===== hyperelastic problem at iteration 1 =====
At Newton iteration 1 norm is 8.627104418308e-02
At Newton iteration 2 norm is 7.783879848245e-04
At Newton iteration 3 norm is 3.817813043733e-08
===== hyperelastic problem at iteration 2 =====
At Newton iteration 1 norm is 6.626617918429e-02
At Newton iteration 2 norm is 6.155599929113e-04
At Newton iteration 3 norm is 2.220433209860e-08
===== hyperelastic problem at iteration 3 =====
At Newton iteration 1 norm is 5.365199085190e-02
At Newton iteration 2 norm is 3.217898859458e-04
At Newton iteration 3 norm is 5.280662666287e-09
===== hyperelastic problem at iteration 4 =====
At Newton iteration 1 norm is 3.448445477588e-02
At Newton iteration 2 norm is 1.651105212249e-04
At Newton iteration 3 norm is 7.904903585581e-10
===== hyperelastic problem at iteration 5 =====
At Newton iteration 1 norm is 3.232023636819e-02
At Newton iteration 2 norm is 9.636747040827e-05
At Newton iteration 3 norm is 1.695388706846e-10
===== hyperelastic problem at iteration 6 =====
At Newton iteration 1 norm is 2.005162696371e-02
At Newton iteration 2 norm is 4.861897171624e-05
At Newton iteration 3 norm is 4.367432338704e-11
===== hyperelastic problem at iteration 7 =====
At Newton iteration 1 norm is 1.879130633303e-02
At Newton iteration 2 norm is 2.350685330318e-05
At Newton iteration 3 norm is 8.781039636974e-12
===== hyperelastic problem at iteration 8 =====
At Newton iteration 1 norm is 1.079876021448e-02
At Newton iteration 2 norm is 1.367831305412e-05
At Newton iteration 3 norm is 2.527809298010e-12
===== hyperelastic problem at iteration 9 =====
At Newton iteration 1 norm is 1.186998995129e-02
At Newton iteration 2 norm is 7.941399724594e-06
==== hyperelastic problem at iteration 10 ====
At Newton iteration 1 norm is 6.296101809138e-03
At Newton iteration 2 norm is 3.410862824393e-06
==== hyperelastic problem at iteration 11 ====
At Newton iteration 1 norm is 7.292086420178e-03
At Newton iteration 2 norm is 2.415955353113e-06
==== hyperelastic problem at iteration 12 ====
At Newton iteration 1 norm is 4.146024277092e-03
At Newton iteration 2 norm is 1.043102533053e-06
==== hyperelastic problem at iteration 13 ====
At Newton iteration 1 norm is 1.012064090368e-02
At Newton iteration 2 norm is 2.954788400160e-06
==== hyperelastic problem at iteration 14 ====
At Newton iteration 1 norm is 1.668500585128e+00
At Newton iteration 2 norm is 8.144096288499e-02
At Newton iteration 3 norm is 4.910819037402e-04
At Newton iteration 4 norm is 1.606760612740e-08
==== hyperelastic problem at iteration 15 ====
At Newton iteration 1 norm is 1.864700262373e+00
At Newton iteration 2 norm is 4.799670120590e-01
At Newton iteration 3 norm is 3.843351422930e-04
At Newton iteration 4 norm is 6.334592679363e-09
==== hyperelastic problem at iteration 16 ====
At Newton iteration 1 norm is 1.561591222740e+00
At Newton iteration 2 norm is 4.993259773779e-01
At Newton iteration 3 norm is 1.623810601839e-04
At Newton iteration 4 norm is 2.416276960278e-09
==== hyperelastic problem at iteration 17 ====
At Newton iteration 1 norm is 1.486987763500e+00
At Newton iteration 2 norm is 3.506808946576e-01
At Newton iteration 3 norm is 1.277215927848e-04
At Newton iteration 4 norm is 1.762255552074e-09
==== hyperelastic problem at iteration 18 ====
At Newton iteration 1 norm is 1.402419150178e+00
At Newton iteration 2 norm is 2.899970761077e-01
At Newton iteration 3 norm is 5.801775436488e-05
At Newton iteration 4 norm is 1.865332615734e-10
==== hyperelastic problem at iteration 19 ====
At Newton iteration 1 norm is 1.370015019157e+00
At Newton iteration 2 norm is 2.645047452374e-01
At Newton iteration 3 norm is 6.059281399009e-05
At Newton iteration 4 norm is 2.377590994952e-10
==== hyperelastic problem at iteration 20 ====
At Newton iteration 1 norm is 1.266025514063e+00
At Newton iteration 2 norm is 2.057743195323e-01
At Newton iteration 3 norm is 3.505876167608e-05
At Newton iteration 4 norm is 6.949942387966e-11
==== hyperelastic problem at iteration 21 ====
At Newton iteration 1 norm is 1.232544598814e+00
At Newton iteration 2 norm is 1.619180149489e-01
At Newton iteration 3 norm is 3.632762760828e-05
At Newton iteration 4 norm is 9.585773023623e-11
==== hyperelastic problem at iteration 22 ====
At Newton iteration 1 norm is 1.129384135821e+00
At Newton iteration 2 norm is 1.351449983385e-01
At Newton iteration 3 norm is 2.312182284610e-05
At Newton iteration 4 norm is 4.275772660276e-11
==== hyperelastic problem at iteration 23 ====
At Newton iteration 1 norm is 1.137071588823e+00
At Newton iteration 2 norm is 1.068045114488e-01
At Newton iteration 3 norm is 2.712726181860e-05
At Newton iteration 4 norm is 6.155209824573e-11
==== hyperelastic problem at iteration 24 ====
At Newton iteration 1 norm is 1.002514755093e+00
At Newton iteration 2 norm is 8.252330081836e-02
At Newton iteration 3 norm is 1.515170364713e-05
At Newton iteration 4 norm is 2.177543656872e-11
==== hyperelastic problem at iteration 25 ====
At Newton iteration 1 norm is 1.020288334802e+00
At Newton iteration 2 norm is 6.770902107195e-02
At Newton iteration 3 norm is 1.556709343493e-05
At Newton iteration 4 norm is 2.182322492913e-11
==== hyperelastic problem at iteration 26 ====
At Newton iteration 1 norm is 8.994396136678e-01
At Newton iteration 2 norm is 5.177220568847e-02
At Newton iteration 3 norm is 8.393353337694e-06
==== hyperelastic problem at iteration 27 ====
At Newton iteration 1 norm is 9.030962138815e-01
At Newton iteration 2 norm is 4.106379099813e-02
At Newton iteration 3 norm is 8.120305110433e-06
==== hyperelastic problem at iteration 28 ====
At Newton iteration 1 norm is 7.944549454180e-01
At Newton iteration 2 norm is 3.153698026084e-02
At Newton iteration 3 norm is 4.402181396957e-06
==== hyperelastic problem at iteration 29 ====
At Newton iteration 1 norm is 7.945841772426e-01
At Newton iteration 2 norm is 2.484569345572e-02
At Newton iteration 3 norm is 4.232181201768e-06
==== hyperelastic problem at iteration 30 ====
At Newton iteration 1 norm is 6.978326590220e-01
At Newton iteration 2 norm is 1.916305658325e-02
At Newton iteration 3 norm is 2.340905302667e-06
==== hyperelastic problem at iteration 31 ====
At Newton iteration 1 norm is 6.959190103828e-01
At Newton iteration 2 norm is 1.535135817434e-02
At Newton iteration 3 norm is 2.229548738171e-06
==== hyperelastic problem at iteration 32 ====
At Newton iteration 1 norm is 6.148350779583e-01
At Newton iteration 2 norm is 1.218723825883e-02
At Newton iteration 3 norm is 1.295377321038e-06
==== hyperelastic problem at iteration 33 ====
At Newton iteration 1 norm is 6.087707482561e-01
At Newton iteration 2 norm is 9.981387860338e-03
At Newton iteration 3 norm is 1.198289818118e-06
==== hyperelastic problem at iteration 34 ====
At Newton iteration 1 norm is 5.440994521029e-01
At Newton iteration 2 norm is 8.226683861519e-03
At Newton iteration 3 norm is 7.387502323566e-07
==== hyperelastic problem at iteration 35 ====
At Newton iteration 1 norm is 5.318692836384e-01
At Newton iteration 2 norm is 6.882257891791e-03
At Newton iteration 3 norm is 6.580372741413e-07
==== hyperelastic problem at iteration 36 ====
At Newton iteration 1 norm is 4.838829358096e-01
At Newton iteration 2 norm is 5.853381989491e-03
At Newton iteration 3 norm is 4.294597536441e-07
==== hyperelastic problem at iteration 37 ====
At Newton iteration 1 norm is 4.631346903222e-01
At Newton iteration 2 norm is 4.948881865972e-03
At Newton iteration 3 norm is 3.653153313352e-07
==== hyperelastic problem at iteration 38 ====
At Newton iteration 1 norm is 3.736270457552e-01
At Newton iteration 2 norm is 3.373137985145e-03
At Newton iteration 3 norm is 1.494676499701e-07
==== hyperelastic problem at iteration 39 ====
At Newton iteration 1 norm is 6.704807451418e+00
At Newton iteration 2 norm is 3.058130447115e-01
At Newton iteration 3 norm is 1.163816551455e-03
At Newton iteration 4 norm is 7.605399175078e-08
==== hyperelastic problem at iteration 40 ====
At Newton iteration 1 norm is 6.165291742807e+00
At Newton iteration 2 norm is 1.234224662142e+00
At Newton iteration 3 norm is 1.072820113297e-02
At Newton iteration 4 norm is 7.189251973393e-06
==== hyperelastic problem at iteration 41 ====
At Newton iteration 1 norm is 4.057648321975e+00
At Newton iteration 2 norm is 1.204701586988e+00
At Newton iteration 3 norm is 4.177049374979e-03
At Newton iteration 4 norm is 1.421998616507e-06
==== hyperelastic problem at iteration 42 ====
At Newton iteration 1 norm is 3.079002105428e+00
At Newton iteration 2 norm is 7.034866918116e-01
At Newton iteration 3 norm is 1.233401229724e-03
At Newton iteration 4 norm is 1.254613192790e-07
==== hyperelastic problem at iteration 43 ====
At Newton iteration 1 norm is 1.884407486835e+00
At Newton iteration 2 norm is 5.230490771215e-01
At Newton iteration 3 norm is 3.210529552672e-04
At Newton iteration 4 norm is 8.691995924511e-09
==== hyperelastic problem at iteration 44 ====
At Newton iteration 1 norm is 1.805652277006e+00
At Newton iteration 2 norm is 3.317918404334e-01
At Newton iteration 3 norm is 2.374305946549e-04
At Newton iteration 4 norm is 5.106988780644e-09
==== hyperelastic problem at iteration 45 ====
At Newton iteration 1 norm is 1.333568056084e+00
At Newton iteration 2 norm is 2.217189655145e-01
At Newton iteration 3 norm is 8.713697103929e-05
At Newton iteration 4 norm is 6.365856600309e-10
==== hyperelastic problem at iteration 46 ====
At Newton iteration 1 norm is 1.262242782626e+00
At Newton iteration 2 norm is 1.260243805086e-01
At Newton iteration 3 norm is 6.096424367479e-05
At Newton iteration 4 norm is 3.609231695559e-10
==== hyperelastic problem at iteration 47 ====
At Newton iteration 1 norm is 1.136765084211e+00
At Newton iteration 2 norm is 8.978774810143e-02
At Newton iteration 3 norm is 2.005211249672e-05
At Newton iteration 4 norm is 2.050186738598e-11
==== hyperelastic problem at iteration 48 ====
At Newton iteration 1 norm is 1.023859872386e+00
At Newton iteration 2 norm is 4.541939147262e-02
At Newton iteration 3 norm is 1.357368093771e-05
At Newton iteration 4 norm is 1.388272020651e-11
==== hyperelastic problem at iteration 49 ====
At Newton iteration 1 norm is 1.012948847924e+00
At Newton iteration 2 norm is 3.209181262697e-02
At Newton iteration 3 norm is 4.608167756439e-06
==== hyperelastic problem at iteration 50 ====
At Newton iteration 1 norm is 9.011848606483e-01
At Newton iteration 2 norm is 1.514114435145e-02
At Newton iteration 3 norm is 2.649317509595e-06
==== hyperelastic problem at iteration 51 ====
At Newton iteration 1 norm is 9.147499543620e-01
At Newton iteration 2 norm is 1.182898822835e-02
At Newton iteration 3 norm is 9.585085960859e-07
==== hyperelastic problem at iteration 52 ====
At Newton iteration 1 norm is 8.262571045568e-01
At Newton iteration 2 norm is 5.116555475740e-03
At Newton iteration 3 norm is 5.523169250187e-07
==== hyperelastic problem at iteration 53 ====
At Newton iteration 1 norm is 8.340493842286e-01
At Newton iteration 2 norm is 4.510645600581e-03
At Newton iteration 3 norm is 3.220947802889e-07
==== hyperelastic problem at iteration 54 ====
At Newton iteration 1 norm is 7.673833114224e-01
At Newton iteration 2 norm is 1.893075604962e-03
At Newton iteration 3 norm is 1.017537400198e-07
==== hyperelastic problem at iteration 55 ====
At Newton iteration 1 norm is 7.658377469920e-01
At Newton iteration 2 norm is 1.806908734077e-03
At Newton iteration 3 norm is 1.232743151761e-07
==== hyperelastic problem at iteration 56 ====
At Newton iteration 1 norm is 7.143977434337e-01
At Newton iteration 2 norm is 8.793858776794e-04
At Newton iteration 3 norm is 5.046029910889e-08
==== hyperelastic problem at iteration 57 ====
At Newton iteration 1 norm is 7.035265861309e-01
At Newton iteration 2 norm is 7.682477952640e-04
At Newton iteration 3 norm is 3.358157699911e-08
==== hyperelastic problem at iteration 58 ====
At Newton iteration 1 norm is 6.629224135079e-01
At Newton iteration 2 norm is 4.467451354865e-04
At Newton iteration 3 norm is 3.860050859672e-08
==== hyperelastic problem at iteration 59 ====
At Newton iteration 1 norm is 6.449768305255e-01
At Newton iteration 2 norm is 3.411592777240e-04
At Newton iteration 3 norm is 6.323641182039e-09
==== hyperelastic problem at iteration 60 ====
At Newton iteration 1 norm is 6.115152876412e-01
At Newton iteration 2 norm is 2.043579222177e-04
At Newton iteration 3 norm is 3.770304606054e-09
==== hyperelastic problem at iteration 61 ====
At Newton iteration 1 norm is 5.895107293471e-01
At Newton iteration 2 norm is 1.758866635145e-04
At Newton iteration 3 norm is 9.623302568997e-10
==== hyperelastic problem at iteration 62 ====
At Newton iteration 1 norm is 5.598658875646e-01
At Newton iteration 2 norm is 1.091781553641e-04
At Newton iteration 3 norm is 4.444991781652e-10
==== hyperelastic problem at iteration 63 ====
At Newton iteration 1 norm is 5.359774646064e-01
At Newton iteration 2 norm is 1.067069369856e-04
At Newton iteration 3 norm is 3.663179626993e-10
==== hyperelastic problem at iteration 64 ====
At Newton iteration 1 norm is 5.083926956332e-01
At Newton iteration 2 norm is 7.128406133681e-05
At Newton iteration 3 norm is 1.401431295383e-10
==== hyperelastic problem at iteration 65 ====
At Newton iteration 1 norm is 4.838931730274e-01
At Newton iteration 2 norm is 7.282729287596e-05
At Newton iteration 3 norm is 2.249769469348e-10
==== hyperelastic problem at iteration 66 ====
At Newton iteration 1 norm is 4.577219107045e-01
At Newton iteration 2 norm is 5.496200771711e-05
At Newton iteration 3 norm is 1.184096267569e-10
==== hyperelastic problem at iteration 67 ====
At Newton iteration 1 norm is 4.333488299598e-01
At Newton iteration 2 norm is 5.584181487155e-05
At Newton iteration 3 norm is 1.599480318462e-10
==== hyperelastic problem at iteration 68 ====
At Newton iteration 1 norm is 4.084094019859e-01
At Newton iteration 2 norm is 4.546875915601e-05
At Newton iteration 3 norm is 1.064058454259e-10
==== hyperelastic problem at iteration 69 ====
At Newton iteration 1 norm is 3.846636417678e-01
At Newton iteration 2 norm is 4.610317427581e-05
At Newton iteration 3 norm is 1.126593035097e-10
==== hyperelastic problem at iteration 70 ====
At Newton iteration 1 norm is 3.610055281702e-01
At Newton iteration 2 norm is 3.976958232826e-05
At Newton iteration 3 norm is 8.161407068905e-11
==== hyperelastic problem at iteration 71 ====
At Newton iteration 1 norm is 3.382929501894e-01
At Newton iteration 2 norm is 3.967652130545e-05
At Newton iteration 3 norm is 7.476727071146e-11
==== hyperelastic problem at iteration 72 ====
At Newton iteration 1 norm is 3.160822917140e-01
At Newton iteration 2 norm is 3.551852763531e-05
At Newton iteration 3 norm is 5.526866323652e-11
==== hyperelastic problem at iteration 73 ====
At Newton iteration 1 norm is 2.947598112758e-01
At Newton iteration 2 norm is 3.462575531913e-05
At Newton iteration 3 norm is 4.659129433842e-11
==== hyperelastic problem at iteration 74 ====
At Newton iteration 1 norm is 2.742313667930e-01
At Newton iteration 2 norm is 3.155568012090e-05
At Newton iteration 3 norm is 3.430819861929e-11
==== hyperelastic problem at iteration 75 ====
At Newton iteration 1 norm is 2.546482912632e-01
At Newton iteration 2 norm is 3.023962475218e-05
At Newton iteration 3 norm is 2.732291361128e-11
==== hyperelastic problem at iteration 76 ====
At Newton iteration 1 norm is 2.360972905034e-01
At Newton iteration 2 norm is 2.777899516777e-05
At Newton iteration 3 norm is 1.987429435907e-11
==== hyperelastic problem at iteration 77 ====
At Newton iteration 1 norm is 2.186241117195e-01
At Newton iteration 2 norm is 2.628079129927e-05
At Newton iteration 3 norm is 1.525485890433e-11
==== hyperelastic problem at iteration 78 ====
At Newton iteration 1 norm is 2.024069659626e-01
At Newton iteration 2 norm is 2.424610478122e-05
At Newton iteration 3 norm is 1.116391676533e-11
==== hyperelastic problem at iteration 79 ====
At Newton iteration 1 norm is 1.874644769478e-01
At Newton iteration 2 norm is 2.277558533684e-05
At Newton iteration 3 norm is 8.841991736615e-12
==== hyperelastic problem at iteration 80 ====
At Newton iteration 1 norm is 1.740276406620e-01
At Newton iteration 2 norm is 2.106130764842e-05
At Newton iteration 3 norm is 7.088083490341e-12
Saved mesh data in:
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.11.0/examples/CardiacElastoDynamics/results/ellipsoid_0.03_D_dt0.0125_analytic_mass/Meshes
Saved simulation and summary in:
/home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.11.0/examples/CardiacElastoDynamics/results/ellipsoid_0.03_D_dt0.0125_analytic_mass
Generate movie 01/20 (5.00 %) 2.93 s
Generate movie 02/20 (10.00 %) 2.40 s
Generate movie 03/20 (15.00 %) 2.26 s
Generate movie 04/20 (20.00 %) 2.12 s
Generate movie 05/20 (25.00 %) 1.94 s
Generate movie 06/20 (30.00 %) 1.88 s
Generate movie 07/20 (35.00 %) 1.64 s
Generate movie 08/20 (40.00 %) 1.50 s
Generate movie 09/20 (45.00 %) 1.37 s
Generate movie 10/20 (50.00 %) 1.24 s
Generate movie 11/20 (55.00 %) 1.15 s
Generate movie 12/20 (60.00 %) 1.05 s
Generate movie 13/20 (65.00 %) 925.11 ms
Generate movie 14/20 (70.00 %) 788.25 ms
Generate movie 15/20 (75.00 %) 656.24 ms
Generate movie 16/20 (80.00 %) 530.01 ms
Generate movie 17/20 (85.00 %) 397.80 ms
Generate movie 18/20 (90.00 %) 265.50 ms
Generate movie 19/20 (95.00 %) 132.69 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.PositiveScalarParameter()
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 displacement,
61 groupElem,
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 D = "D" # active_stress + pressure
111 A = "A" # active_stress
112 B = "B" # pressure
113
114
115 if __name__ == "__main__":
116
117 Display.Clear()
118
119 # ----------------------------------------------
120 # Config
121 # ----------------------------------------------
122
123 useCoarseConfig = True
124
125 ellipsoid = "ellipsoid_0.03" if useCoarseConfig else "ellipsoid_0.005"
126
127 config = Config.D
128
129 fiberSource = "analytic"
130 # fiberSource = "vtu"
131
132 matrixType = MatrixType.mass
133 # matrixType = 15
134
135 results_dir = Folder.Join(RESULTS_DIR, ellipsoid + f"_{config.name}")
136
137 doSimu = True
138 plotGraph = False
139 plotParticles = True
140 saveParticles = True
141 makeMovie = True
142
143 # ----------------------------------------------
144 # time-history needed for plotting in both doSimu / Load_Simu flows
145
146 Nt = 80 if useCoarseConfig else 1000
147
148 t_values, activeStress_values, pressure_values = Get_values(Tmax=1.0, Nt=Nt)
149 dt = t_values[1] - t_values[0]
150 results_dir += f"_dt{dt}_{fiberSource}_{matrixType}"
151
152 if plotGraph:
153 ax = Display.Init_Axes()
154 ax.grid()
155 ax.set_xlabel(r"$t$ [s]")
156 ax.set_ylabel(r"$\tau(t)$ [Pa]")
157 ax.plot(t_values, activeStress_values)
158 name = "active_pressure"
159 Display.Save_fig(results_dir, name)
160
161 ax = Display.Init_Axes()
162 ax.grid()
163 ax.set_xlabel(r"$t$ [s]")
164 ax.set_ylabel(r"$p(t)$ [Pa]")
165 ax.plot(t_values, pressure_values)
166 name = "pressure"
167 Display.Save_fig(results_dir, name)
168
169 if config is Config.B:
170 activeStress_values *= 0
171 if config is Config.A:
172 pressure_values *= 0
173
174 if doSimu:
175
176 # ----------------------------------------------
177 # Mesh, fibers and sheets
178 # ----------------------------------------------
179
180 mesh, fibers_e_pg, sheets_e_pg = Get_config(
181 Folder.Join(DATA_DIR, ellipsoid),
182 matrixType=matrixType,
183 fiberSource=fiberSource,
184 plotMesh=False,
185 plotTags=False,
186 plotFibers=False,
187 )
188
189 # ----------------------------------------------
190 # Material
191 # ----------------------------------------------
192
193 # solid
194 a = 59.0
195 a_f = 18472.0
196 a_fs = 216.0
197 a_s = 2481.0
198 b = 8.023
199 b_f = 16.026
200 b_fs = 11.436
201 b_s = 11.12
202
203 material = Models.HyperElastic.HolzapfelOgden(
204 dim=3,
205 C0=a / 2 / b,
206 C1=b,
207 C2=a_f / 2 / b_f,
208 C3=b_f,
209 C4=a_s / 2 / b_s,
210 C5=b_s,
211 C6=a_fs / 2 / b_fs,
212 C7=b_fs,
213 K=1e6,
214 Mu1=0.0,
215 Mu2=0.0,
216 T1=fibers_e_pg,
217 T2=sheets_e_pg,
218 ks=100,
219 )
220 material.eta = 100.0
221 material.Set_active_stress_vec(material.T1)
222
223 # ----------------------------------------------
224 # Simulation
225 # ----------------------------------------------
226
227 simu = CardiacElastoDynamics(mesh, material, folder=results_dir)
228
229 simu.Solver_Set_Hyperbolic_Algorithm(dt, algo=AlgoType.midpoint)
230 simu.rho = 1000
231
232 for t in t_values:
233 simu.Bc_Init()
234 simu.pressure = np.interp(t + dt / 2, t_values, pressure_values)
235 material.active_stress = np.interp(
236 t + dt / 2, t_values, activeStress_values
237 )
238 simu.Solve()
239 simu.Save_Iter()
240
241 simu.Save(results_dir)
242
243 else:
244 simu = Simulations.Load_Simu(results_dir)
245
246 simu._Gather()
247
248 if plotParticles and simu.isGathered:
249
250 coords = [(0.025, 0.03, 0), (0, 0.03, 0)]
251 evalCoords = np.array(coords)
252 evalElements = simu.mesh.groupElem._Get_nearby_elements(evalCoords)
253
254 Niter = simu.Niter
255 values = np.empty((Niter, len(coords), 3))
256 for i in range(Niter):
257 simu.Set_Iter(i)
258 values[i] = simu.mesh.Evaluate_dofsValues_at_coordinates(
259 evalCoords, simu.displacement, evalElements
260 )
261
262 times = t_values[:Niter]
263 axs = Display.plt.subplots(3, 2, sharex=True)[1]
264
265 for p, (particle, coord) in enumerate(zip(["p0", "p1"], coords)):
266
267 for c, component in enumerate(["x", "y", "z"]):
268
269 ax: Display.plt.Axes = axs[c, p]
270
271 ax.grid()
272
273 if c == 2:
274 ax.set_xlabel("Time [s]")
275 if p == 0:
276 ax.set_ylabel(f"Displacement {component}-component [m]")
277 if c == 0:
278 ax.set_title(f"Particle {particle}")
279
280 ax.plot(times, values[:, p, c])
281
282 width, height = ax.figure.get_size_inches()
283 ax.figure.set_size_inches(width * 1.5, height * 2.5)
284 Display.Save_fig(results_dir, "particles")
285
286 if saveParticles and simu.isGathered:
287
288 # per-iteration deformed volume
289 volumes = np.empty(Niter)
290 for i in range(Niter):
291 simu.Set_Iter(i)
292 deformed = simu.mesh.copy()
293 deformed.coord += simu.displacement.reshape(-1, 3)
294 volumes[i] = deformed.volume
295
296 dict_particles = {
297 "time": times,
298 "displacement": {
299 f"p{p}": {
300 "ux": values[:, p, 0],
301 "uy": values[:, p, 1],
302 "uz": values[:, p, 2],
303 "magnitude": np.linalg.norm(values[:, p, :], axis=1),
304 }
305 for p in range(2)
306 },
307 "stress": {
308 "time": None,
309 "p0": {"magnitude": None},
310 "p1": {"magnitude": None},
311 },
312 "volume": volumes,
313 }
314 Simulations.Save_pickle(dict_particles, results_dir, "particles")
315
316 if makeMovie:
317 values = [simu.Result("ux", iter=i) for i in range(simu.Niter)]
318 clim = (np.min(values), np.max(values))
319 PyVista.Movie_simu(
320 simu,
321 "ux",
322 results_dir,
323 "ux.gif",
324 N=20,
325 deformFactor=1.0,
326 clim=clim,
327 plotMesh=True,
328 )
329
330 Display.plt.show()
Total running time of the script: (0 minutes 19.953 seconds)