.. 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 ``following 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-342 .. 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.195009922818e-09 ===== hyperelastic problem at iteration 1 ===== At Newton iteration 1 norm is 8.636333907881e-02 At Newton iteration 2 norm is 8.112107298115e-04 At Newton iteration 3 norm is 3.260029883637e-08 ===== hyperelastic problem at iteration 2 ===== At Newton iteration 1 norm is 6.667722093494e-02 At Newton iteration 2 norm is 6.598459929744e-04 At Newton iteration 3 norm is 1.979972955817e-08 ===== hyperelastic problem at iteration 3 ===== At Newton iteration 1 norm is 5.443853101852e-02 At Newton iteration 2 norm is 3.574030065331e-04 At Newton iteration 3 norm is 5.100387803189e-09 ===== hyperelastic problem at iteration 4 ===== At Newton iteration 1 norm is 3.568161669214e-02 At Newton iteration 2 norm is 1.879472302824e-04 At Newton iteration 3 norm is 7.634436538104e-10 ===== hyperelastic problem at iteration 5 ===== At Newton iteration 1 norm is 3.345855815030e-02 At Newton iteration 2 norm is 1.104480380067e-04 At Newton iteration 3 norm is 1.727237827593e-10 ===== hyperelastic problem at iteration 6 ===== At Newton iteration 1 norm is 2.145570404913e-02 At Newton iteration 2 norm is 5.743376125986e-05 At Newton iteration 3 norm is 4.241006735523e-11 ===== hyperelastic problem at iteration 7 ===== At Newton iteration 1 norm is 2.002055496982e-02 At Newton iteration 2 norm is 2.858175613789e-05 At Newton iteration 3 norm is 9.329252410960e-12 ===== hyperelastic problem at iteration 8 ===== At Newton iteration 1 norm is 1.215316851397e-02 At Newton iteration 2 norm is 1.697272618343e-05 At Newton iteration 3 norm is 3.483195564703e-12 ===== hyperelastic problem at iteration 9 ===== At Newton iteration 1 norm is 1.284612269409e-02 At Newton iteration 2 norm is 1.004326521727e-05 At Newton iteration 3 norm is 1.236643848330e-12 ==== hyperelastic problem at iteration 10 ==== At Newton iteration 1 norm is 7.358674634723e-03 At Newton iteration 2 norm is 4.663340717413e-06 ==== hyperelastic problem at iteration 11 ==== At Newton iteration 1 norm is 7.987532659065e-03 At Newton iteration 2 norm is 3.232423353815e-06 ==== hyperelastic problem at iteration 12 ==== At Newton iteration 1 norm is 4.809825935357e-03 At Newton iteration 2 norm is 1.536260502767e-06 ==== hyperelastic problem at iteration 13 ==== At Newton iteration 1 norm is 1.039439603744e-02 At Newton iteration 2 norm is 3.486662051805e-06 ==== hyperelastic problem at iteration 14 ==== At Newton iteration 1 norm is 1.766079917895e+00 At Newton iteration 2 norm is 9.125789691990e-02 At Newton iteration 3 norm is 1.087619504333e-03 At Newton iteration 4 norm is 2.420674925161e-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.637389813391e-04 At Newton iteration 4 norm is 1.322538469399e-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.109056791491e-04 At Newton iteration 4 norm is 1.090220825340e-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.687518365728e-05 At Newton iteration 4 norm is 7.305379903751e-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.438425390513e-05 At Newton iteration 4 norm is 1.614494114363e-10 ==== hyperelastic problem at iteration 19 ==== At Newton iteration 1 norm is 1.528855945657e+00 At Newton iteration 2 norm is 2.503221856980e-01 At Newton iteration 3 norm is 6.025836142854e-05 At Newton iteration 4 norm is 1.927383623544e-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.371244262386e-05 At Newton iteration 4 norm is 8.945000499470e-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.447758718910e-05 At Newton iteration 4 norm is 1.275834939250e-10 ==== hyperelastic problem at iteration 22 ==== At Newton iteration 1 norm is 1.263965377752e+00 At Newton iteration 2 norm is 1.478516823518e-01 At Newton iteration 3 norm is 3.463498059109e-05 At Newton iteration 4 norm is 9.066723230797e-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.773930957736e-05 At Newton iteration 4 norm is 1.201000534770e-10 ==== hyperelastic problem at iteration 24 ==== At Newton iteration 1 norm is 1.138920590138e+00 At Newton iteration 2 norm is 9.680485073072e-02 At Newton iteration 3 norm is 2.354890713975e-05 At Newton iteration 4 norm is 5.134185464888e-11 ==== hyperelastic problem at iteration 25 ==== At Newton iteration 1 norm is 1.150907660345e+00 At Newton iteration 2 norm is 8.160509659659e-02 At Newton iteration 3 norm is 2.098969821691e-05 At Newton iteration 4 norm is 3.973533050176e-11 ==== hyperelastic problem at iteration 26 ==== At Newton iteration 1 norm is 1.026993249552e+00 At Newton iteration 2 norm is 6.433957030618e-02 At Newton iteration 3 norm is 1.221062639626e-05 At Newton iteration 4 norm is 1.357096359652e-11 ==== hyperelastic problem at iteration 27 ==== At Newton iteration 1 norm is 1.016718457051e+00 At Newton iteration 2 norm is 5.201368367357e-02 At Newton iteration 3 norm is 1.036869328330e-05 At Newton iteration 4 norm is 9.039831653444e-12 ==== hyperelastic problem at iteration 28 ==== At Newton iteration 1 norm is 9.079659358230e-01 At Newton iteration 2 norm is 4.111574615343e-02 At Newton iteration 3 norm is 6.121538198192e-06 ==== hyperelastic problem at iteration 29 ==== At Newton iteration 1 norm is 8.911742704472e-01 At Newton iteration 2 norm is 3.252948938893e-02 At Newton iteration 3 norm is 5.246258481852e-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.138421966246e-06 ==== hyperelastic problem at iteration 31 ==== At Newton iteration 1 norm is 7.774807157865e-01 At Newton iteration 2 norm is 2.033244512176e-02 At Newton iteration 3 norm is 2.720350001102e-06 ==== hyperelastic problem at iteration 32 ==== At Newton iteration 1 norm is 6.966625324298e-01 At Newton iteration 2 norm is 1.638331490730e-02 At Newton iteration 3 norm is 1.672881886969e-06 ==== hyperelastic problem at iteration 33 ==== At Newton iteration 1 norm is 6.777500156535e-01 At Newton iteration 2 norm is 1.304078137851e-02 At Newton iteration 3 norm is 1.444717982399e-06 ==== hyperelastic problem at iteration 34 ==== At Newton iteration 1 norm is 6.130620615770e-01 At Newton iteration 2 norm is 1.077692758014e-02 At Newton iteration 3 norm is 9.146600434499e-07 ==== hyperelastic problem at iteration 35 ==== At Newton iteration 1 norm is 5.902687207387e-01 At Newton iteration 2 norm is 8.705891529260e-03 At Newton iteration 3 norm is 7.830889352535e-07 ==== hyperelastic problem at iteration 36 ==== At Newton iteration 1 norm is 5.422373569372e-01 At Newton iteration 2 norm is 7.351341456339e-03 At Newton iteration 3 norm is 5.100601649949e-07 ==== hyperelastic problem at iteration 37 ==== At Newton iteration 1 norm is 5.124490985705e-01 At Newton iteration 2 norm is 6.024003824768e-03 At Newton iteration 3 norm is 4.294865320165e-07 ==== hyperelastic problem at iteration 38 ==== At Newton iteration 1 norm is 4.155515550564e-01 At Newton iteration 2 norm is 4.091144051917e-03 At Newton iteration 3 norm is 1.715786262025e-07 ==== hyperelastic problem at iteration 39 ==== At Newton iteration 1 norm is 7.223803002004e+00 At Newton iteration 2 norm is 3.374928378035e-01 At Newton iteration 3 norm is 1.512215943614e-03 At Newton iteration 4 norm is 1.209413613077e-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.188459381137e-02 At Newton iteration 4 norm is 9.884492877704e-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.626756732999e-03 At Newton iteration 4 norm is 1.977426702712e-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.294449634168e-03 At Newton iteration 4 norm is 1.571831071135e-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.286824519142e-04 At Newton iteration 4 norm is 8.291984702509e-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.113008822774e-04 At Newton iteration 4 norm is 3.912718376848e-09 ==== hyperelastic problem at iteration 45 ==== At Newton iteration 1 norm is 1.334126825250e+00 At Newton iteration 2 norm is 1.871344231620e-01 At Newton iteration 3 norm is 7.617148575375e-05 At Newton iteration 4 norm is 4.128086973086e-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.908254974208e-05 At Newton iteration 4 norm is 2.065097539425e-10 ==== hyperelastic problem at iteration 47 ==== At Newton iteration 1 norm is 1.105992962329e+00 At Newton iteration 2 norm is 7.189486787907e-02 At Newton iteration 3 norm is 1.594932226076e-05 At Newton iteration 4 norm is 8.947918140863e-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.042709878866e-05 At Newton iteration 4 norm is 5.660691005643e-12 ==== hyperelastic problem at iteration 49 ==== At Newton iteration 1 norm is 9.805659860480e-01 At Newton iteration 2 norm is 2.521242433352e-02 At Newton iteration 3 norm is 3.338387369209e-06 ==== hyperelastic problem at iteration 50 ==== At Newton iteration 1 norm is 9.015430803165e-01 At Newton iteration 2 norm is 1.236194514671e-02 At Newton iteration 3 norm is 2.147490516882e-06 ==== hyperelastic problem at iteration 51 ==== At Newton iteration 1 norm is 8.888258892399e-01 At Newton iteration 2 norm is 9.729822282724e-03 At Newton iteration 3 norm is 7.514945099176e-07 ==== hyperelastic problem at iteration 52 ==== At Newton iteration 1 norm is 8.174872022538e-01 At Newton iteration 2 norm is 4.528749982911e-03 At Newton iteration 3 norm is 5.045362567222e-07 ==== hyperelastic problem at iteration 53 ==== At Newton iteration 1 norm is 8.132221796678e-01 At Newton iteration 2 norm is 3.991324564921e-03 At Newton iteration 3 norm is 2.558138578404e-07 ==== hyperelastic problem at iteration 54 ==== At Newton iteration 1 norm is 7.548929810226e-01 At Newton iteration 2 norm is 1.860699119017e-03 At Newton iteration 3 norm is 8.104728099265e-08 ==== hyperelastic problem at iteration 55 ==== At Newton iteration 1 norm is 7.481044324089e-01 At Newton iteration 2 norm is 1.727590418353e-03 At Newton iteration 3 norm is 5.737153373598e-08 ==== hyperelastic problem at iteration 56 ==== At Newton iteration 1 norm is 7.008317160020e-01 At Newton iteration 2 norm is 9.127256508083e-04 At Newton iteration 3 norm is 6.812612494630e-08 ==== hyperelastic problem at iteration 57 ==== At Newton iteration 1 norm is 6.880730859717e-01 At Newton iteration 2 norm is 7.777542154026e-04 At Newton iteration 3 norm is 2.432568218601e-08 ==== hyperelastic problem at iteration 58 ==== At Newton iteration 1 norm is 6.497186901320e-01 At Newton iteration 2 norm is 4.614777707682e-04 At Newton iteration 3 norm is 2.275536183043e-08 ==== hyperelastic problem at iteration 59 ==== At Newton iteration 1 norm is 6.317440138225e-01 At Newton iteration 2 norm is 3.629149331610e-04 At Newton iteration 3 norm is 4.184008359101e-09 ==== hyperelastic problem at iteration 60 ==== At Newton iteration 1 norm is 5.995733407555e-01 At Newton iteration 2 norm is 2.346991735229e-04 At Newton iteration 3 norm is 3.682131550407e-09 ==== hyperelastic problem at iteration 61 ==== At Newton iteration 1 norm is 5.785810082564e-01 At Newton iteration 2 norm is 1.973183190452e-04 At Newton iteration 3 norm is 1.075611960852e-09 ==== hyperelastic problem at iteration 62 ==== At Newton iteration 1 norm is 5.500050438403e-01 At Newton iteration 2 norm is 1.313525537202e-04 At Newton iteration 3 norm is 8.066876155885e-10 ==== hyperelastic problem at iteration 63 ==== At Newton iteration 1 norm is 5.275749142850e-01 At Newton iteration 2 norm is 1.234968809272e-04 At Newton iteration 3 norm is 5.010300824091e-10 ==== hyperelastic problem at iteration 64 ==== At Newton iteration 1 norm is 5.011748613697e-01 At Newton iteration 2 norm is 8.399457910998e-05 At Newton iteration 3 norm is 2.632478249504e-10 ==== hyperelastic problem at iteration 65 ==== At Newton iteration 1 norm is 4.782874190070e-01 At Newton iteration 2 norm is 8.538851801639e-05 At Newton iteration 3 norm is 3.156901250987e-10 ==== hyperelastic problem at iteration 66 ==== At Newton iteration 1 norm is 4.534428032311e-01 At Newton iteration 2 norm is 6.310035831464e-05 At Newton iteration 3 norm is 1.597719310355e-10 ==== hyperelastic problem at iteration 67 ==== At Newton iteration 1 norm is 4.307005851754e-01 At Newton iteration 2 norm is 6.570998136604e-05 At Newton iteration 3 norm is 2.207327611751e-10 ==== hyperelastic problem at iteration 68 ==== At Newton iteration 1 norm is 4.071546078103e-01 At Newton iteration 2 norm is 5.213994258629e-05 At Newton iteration 3 norm is 1.363506384678e-10 ==== hyperelastic problem at iteration 69 ==== At Newton iteration 1 norm is 3.849569253036e-01 At Newton iteration 2 norm is 5.431111795866e-05 At Newton iteration 3 norm is 1.550013076227e-10 ==== hyperelastic problem at iteration 70 ==== At Newton iteration 1 norm is 3.626294460454e-01 At Newton iteration 2 norm is 4.555747806763e-05 At Newton iteration 3 norm is 1.070559539459e-10 ==== hyperelastic problem at iteration 71 ==== At Newton iteration 1 norm is 3.412633887488e-01 At Newton iteration 2 norm is 4.639574951534e-05 At Newton iteration 3 norm is 1.029522747784e-10 ==== hyperelastic problem at iteration 72 ==== At Newton iteration 1 norm is 3.201721472453e-01 At Newton iteration 2 norm is 4.053811622740e-05 At Newton iteration 3 norm is 7.348528934907e-11 ==== hyperelastic problem at iteration 73 ==== At Newton iteration 1 norm is 2.998699388741e-01 At Newton iteration 2 norm is 3.989850739047e-05 At Newton iteration 3 norm is 6.344579930693e-11 ==== hyperelastic problem at iteration 74 ==== At Newton iteration 1 norm is 2.801031257124e-01 At Newton iteration 2 norm is 3.548282531170e-05 At Newton iteration 3 norm is 4.525626213656e-11 ==== hyperelastic problem at iteration 75 ==== At Newton iteration 1 norm is 2.610859087223e-01 At Newton iteration 2 norm is 3.388145461169e-05 At Newton iteration 3 norm is 3.593471596824e-11 ==== hyperelastic problem at iteration 76 ==== At Newton iteration 1 norm is 2.427962193749e-01 At Newton iteration 2 norm is 3.024139519656e-05 At Newton iteration 3 norm is 2.509876102040e-11 ==== hyperelastic problem at iteration 77 ==== At Newton iteration 1 norm is 2.253105891341e-01 At Newton iteration 2 norm is 2.819382664981e-05 At Newton iteration 3 norm is 1.866126941999e-11 ==== hyperelastic problem at iteration 78 ==== At Newton iteration 1 norm is 2.087297394796e-01 At Newton iteration 2 norm is 2.519949220411e-05 At Newton iteration 3 norm is 1.285995619943e-11 ==== hyperelastic problem at iteration 79 ==== At Newton iteration 1 norm is 1.930908341347e-01 At Newton iteration 2 norm is 2.318537656879e-05 At Newton iteration 3 norm is 9.632320228728e-12 ==== hyperelastic problem at iteration 80 ==== At Newton iteration 1 norm is 1.785677476811e-01 At Newton iteration 2 norm is 2.083534787055e-05 At Newton iteration 3 norm is 7.422819137161e-12 Saved mesh data in: /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v3.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/v3.0.0/examples/CardiacElastoDynamics/results/step1/ellipsoid0.03_dt0.0125_analytic_rigi Generate movie 01/20 (5.00 %) 2.83 s Generate movie 02/20 (10.00 %) 2.37 s Generate movie 03/20 (15.00 %) 2.21 s Generate movie 04/20 (20.00 %) 2.10 s Generate movie 05/20 (25.00 %) 1.94 s Generate movie 06/20 (30.00 %) 1.77 s Generate movie 07/20 (35.00 %) 1.65 s Generate movie 08/20 (40.00 %) 1.51 s Generate movie 09/20 (45.00 %) 1.37 s Generate movie 10/20 (50.00 %) 1.25 s Generate movie 11/20 (55.00 %) 1.16 s Generate movie 12/20 (60.00 %) 1.04 s Generate movie 13/20 (65.00 %) 912.10 ms Generate movie 14/20 (70.00 %) 787.56 ms Generate movie 15/20 (75.00 %) 652.88 ms Generate movie 16/20 (80.00 %) 527.90 ms Generate movie 17/20 (85.00 %) 391.01 ms Generate movie 18/20 (90.00 %) 260.55 ms Generate movie 19/20 (95.00 %) 129.62 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 ( Terminal, Matplotlib, 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 — tracks the deformed normal, so it depends on # the current iterate / pressure and is rebuilt every Newton step. 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__": Terminal.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 = Matplotlib.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" Matplotlib.Save_fig(results_dir, name) ax = Matplotlib.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" Matplotlib.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 = Matplotlib.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: Matplotlib.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) Matplotlib.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, ) Matplotlib.plt.show() .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 12.879 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 `_