.. 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-332 .. 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.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 | .. 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.ScalarParameter() 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( groupElem, displacement, 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): step0A = "step0A" # active_stress step0B = "step0B" # pressure step1 = "step1" # active_stress + pressure if __name__ == "__main__": Display.Clear() # ---------------------------------------------- # Config # ---------------------------------------------- useCoarseConfig = True ellipsoid = "ellipsoid0.03" if useCoarseConfig else "ellipsoid0.005" config = Config.step1 fiberSource = "analytic" # fiberSource = "vtu" matrixType = MatrixType.rigi # matrixType = MatrixType.mass # matrixType = 15 results_dir = Folder.Join(RESULTS_DIR, config.name, ellipsoid) 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.step0B: activeStress_values *= 0 if config is Config.step0A: 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 13.501 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 `_