.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/CardiacElastoDynamics/MonoVentricular.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_CardiacElastoDynamics_MonoVentricular.py: 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 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. .. GENERATED FROM PYTHON SOURCE LINES 20-331 .. image-sg:: /examples/CardiacElastoDynamics/images/sphx_glr_MonoVentricular_001.gif :alt: MonoVentricular :srcset: /examples/CardiacElastoDynamics/images/sphx_glr_MonoVentricular_001.gif :class: sphx-glr-single-img .. image-sg:: /examples/CardiacElastoDynamics/images/sphx_glr_MonoVentricular_002.png :alt: Particle p0, Particle p1 :srcset: /examples/CardiacElastoDynamics/images/sphx_glr_MonoVentricular_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none ===== 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 | .. code-block:: Python :lineno-start: 21 from enum import Enum import numpy as np from EasyFEA import Display, Folder, PyVista, MatrixType, Models, Simulations, AlgoType from EasyFEA.FEM import Operators from EasyFEA.Utilities import _params from utils import RESULTS_DIR, DATA_DIR, Get_config, Get_values class CardiacElastoDynamics(Simulations.HyperElastic): pressure = _params.PositiveScalarParameter() def __init__( self, mesh, model, folder="", tolConv=0.00001, maxIter=20, verbosity=False ): super().__init__(mesh, model, folder, tolConv, maxIter, verbosity) self.pressure = 0.0 def Construct_local_matrix_system(self, problemType): assert isinstance(self.material, Models.HyperElastic.HolzapfelOgden) nPg = self.material.T1.shape[1] results = super().Construct_local_matrix_system(problemType, nPg) # current Newton-Raphson iterate (updated via u += delta_u) displacement = self._Solver_Get_Newton_Raphson_current_solution() if self.algo in AlgoType.Get_Hyperbolic_Types(): displacement, _, _ = self._Solver_Evaluate_u_v_a_for_time_scheme( problemType, displacement ) for groupElem in self.mesh.Get_list_groupElem(self.dim - 1): # Following pressure (operator returns slot-ready values) tangent_e, residual_e = Operators.NonLinear.FollowingPressure( displacement, groupElem, self.pressure, groupElem.Get_Elements_Tag("endo"), MatrixType.mass, ) # top — isotropic surface mass penalty (Robin α·u + β·u̇ = 0) M_e = Operators.Bilinear.UV(groupElem, dof_n=3) Ktop_e = np.zeros_like(M_e) Ctop_e = np.zeros_like(M_e) if "top" in groupElem.elementTags: top_e = groupElem.Get_Elements_Tag("top") alpha_top = 1e5 Ktop_e[top_e] = alpha_top * M_e[top_e] beta_top = 5e3 Ctop_e[top_e] = beta_top * M_e[top_e] # epi — normal-direction mass penalty (Robin α·(u·n̂) + β·(u̇·n̂) = 0) Ms_e = Operators.Bilinear.MassAlongNormal(groupElem) Kepi_e = np.zeros_like(Ms_e) Cepi_e = np.zeros_like(Ms_e) if "epi" in groupElem.elementTags: epi_e = groupElem.Get_Elements_Tag("epi") alpha_epi = 1e8 Kepi_e[epi_e] = alpha_epi * Ms_e[epi_e] beta_epi = 5e3 Cepi_e[epi_e] = beta_epi * Ms_e[epi_e] # Penalty residual contribution: −K_penalty · u_t at current iterate K_penalty_e = Ktop_e + Kepi_e assembly_e = groupElem.Get_assembly_e(self.dim) u_e = displacement[assembly_e] # (Ne_surf, nPe·3) f_penalty_e = np.einsum("eij,ej->ei", K_penalty_e, u_e) results[groupElem] = ( tangent_e + K_penalty_e, Ctop_e + Cepi_e, None, residual_e - f_penalty_e, ) return results class Config(str, Enum): D = "D" # active_stress + pressure A = "A" # active_stress B = "B" # pressure if __name__ == "__main__": Display.Clear() # ---------------------------------------------- # Config # ---------------------------------------------- useCoarseConfig = True ellipsoid = "ellipsoid_0.03" if useCoarseConfig else "ellipsoid_0.005" config = Config.D fiberSource = "analytic" # fiberSource = "vtu" matrixType = MatrixType.mass # matrixType = 15 results_dir = Folder.Join(RESULTS_DIR, ellipsoid + f"_{config.name}") doSimu = True plotGraph = False plotParticles = True saveParticles = True makeMovie = True # ---------------------------------------------- # time-history needed for plotting in both doSimu / Load_Simu flows Nt = 80 if useCoarseConfig else 1000 t_values, activeStress_values, pressure_values = Get_values(Tmax=1.0, Nt=Nt) dt = t_values[1] - t_values[0] results_dir += f"_dt{dt}_{fiberSource}_{matrixType}" if plotGraph: ax = Display.Init_Axes() ax.grid() ax.set_xlabel(r"$t$ [s]") ax.set_ylabel(r"$\tau(t)$ [Pa]") ax.plot(t_values, activeStress_values) name = "active_pressure" Display.Save_fig(results_dir, name) ax = Display.Init_Axes() ax.grid() ax.set_xlabel(r"$t$ [s]") ax.set_ylabel(r"$p(t)$ [Pa]") ax.plot(t_values, pressure_values) name = "pressure" Display.Save_fig(results_dir, name) if config is Config.B: activeStress_values *= 0 if config is Config.A: pressure_values *= 0 if doSimu: # ---------------------------------------------- # Mesh, fibers and sheets # ---------------------------------------------- mesh, fibers_e_pg, sheets_e_pg = Get_config( Folder.Join(DATA_DIR, ellipsoid), matrixType=matrixType, fiberSource=fiberSource, plotMesh=False, plotTags=False, plotFibers=False, ) # ---------------------------------------------- # Material # ---------------------------------------------- # solid a = 59.0 a_f = 18472.0 a_fs = 216.0 a_s = 2481.0 b = 8.023 b_f = 16.026 b_fs = 11.436 b_s = 11.12 material = Models.HyperElastic.HolzapfelOgden( dim=3, C0=a / 2 / b, C1=b, C2=a_f / 2 / b_f, C3=b_f, C4=a_s / 2 / b_s, C5=b_s, C6=a_fs / 2 / b_fs, C7=b_fs, K=1e6, Mu1=0.0, Mu2=0.0, T1=fibers_e_pg, T2=sheets_e_pg, ks=100, ) material.eta = 100.0 material.Set_active_stress_vec(material.T1) # ---------------------------------------------- # Simulation # ---------------------------------------------- simu = CardiacElastoDynamics(mesh, material, folder=results_dir) simu.Solver_Set_Hyperbolic_Algorithm(dt, algo=AlgoType.midpoint) simu.rho = 1000 for t in t_values: simu.Bc_Init() simu.pressure = np.interp(t + dt / 2, t_values, pressure_values) material.active_stress = np.interp( t + dt / 2, t_values, activeStress_values ) simu.Solve() simu.Save_Iter() simu.Save(results_dir) else: simu = Simulations.Load_Simu(results_dir) simu._Gather() if plotParticles and simu.isGathered: coords = [(0.025, 0.03, 0), (0, 0.03, 0)] evalCoords = np.array(coords) evalElements = simu.mesh.groupElem._Get_nearby_elements(evalCoords) Niter = simu.Niter values = np.empty((Niter, len(coords), 3)) for i in range(Niter): simu.Set_Iter(i) values[i] = simu.mesh.Evaluate_dofsValues_at_coordinates( evalCoords, simu.displacement, evalElements ) times = t_values[:Niter] axs = Display.plt.subplots(3, 2, sharex=True)[1] for p, (particle, coord) in enumerate(zip(["p0", "p1"], coords)): for c, component in enumerate(["x", "y", "z"]): ax: Display.plt.Axes = axs[c, p] ax.grid() if c == 2: ax.set_xlabel("Time [s]") if p == 0: ax.set_ylabel(f"Displacement {component}-component [m]") if c == 0: ax.set_title(f"Particle {particle}") ax.plot(times, values[:, p, c]) width, height = ax.figure.get_size_inches() ax.figure.set_size_inches(width * 1.5, height * 2.5) Display.Save_fig(results_dir, "particles") if saveParticles and simu.isGathered: # per-iteration deformed volume volumes = np.empty(Niter) for i in range(Niter): simu.Set_Iter(i) deformed = simu.mesh.copy() deformed.coord += simu.displacement.reshape(-1, 3) volumes[i] = deformed.volume dict_particles = { "time": times, "displacement": { f"p{p}": { "ux": values[:, p, 0], "uy": values[:, p, 1], "uz": values[:, p, 2], "magnitude": np.linalg.norm(values[:, p, :], axis=1), } for p in range(2) }, "stress": { "time": None, "p0": {"magnitude": None}, "p1": {"magnitude": None}, }, "volume": volumes, } Simulations.Save_pickle(dict_particles, results_dir, "particles") if makeMovie: values = [simu.Result("ux", iter=i) for i in range(simu.Niter)] clim = (np.min(values), np.max(values)) PyVista.Movie_simu( simu, "ux", results_dir, "ux.gif", N=20, deformFactor=1.0, clim=clim, plotMesh=True, ) Display.plt.show() .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 19.953 seconds) .. _sphx_glr_download_examples_CardiacElastoDynamics_MonoVentricular.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: MonoVentricular.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: MonoVentricular.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: MonoVentricular.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_