.. 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-364 .. 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 1.701414174823e-15 ===== hyperelastic problem at iteration 1 ===== At Newton iteration 1 norm is 2.561145112256e-15 ===== hyperelastic problem at iteration 2 ===== At Newton iteration 1 norm is 1.475069784143e-14 ===== hyperelastic problem at iteration 3 ===== At Newton iteration 1 norm is 5.425070561538e-14 ===== hyperelastic problem at iteration 4 ===== At Newton iteration 1 norm is 8.368265512706e-14 ===== hyperelastic problem at iteration 5 ===== At Newton iteration 1 norm is 1.113620833159e-13 ===== hyperelastic problem at iteration 6 ===== At Newton iteration 1 norm is 1.627322478676e-13 ===== hyperelastic problem at iteration 7 ===== At Newton iteration 1 norm is 2.248398965079e-13 ===== hyperelastic problem at iteration 8 ===== At Newton iteration 1 norm is 2.996327547142e-13 ===== hyperelastic problem at iteration 9 ===== At Newton iteration 1 norm is 3.549401174943e-13 ==== hyperelastic problem at iteration 10 ==== At Newton iteration 1 norm is 4.281129218549e-13 ==== hyperelastic problem at iteration 11 ==== At Newton iteration 1 norm is 5.210038574712e-13 ==== hyperelastic problem at iteration 12 ==== At Newton iteration 1 norm is 6.100341534836e-13 ==== hyperelastic problem at iteration 13 ==== At Newton iteration 1 norm is 6.639940554400e-13 ==== hyperelastic problem at iteration 14 ==== At Newton iteration 1 norm is 2.121145006475e+00 At Newton iteration 2 norm is 2.244664365385e-01 At Newton iteration 3 norm is 1.822228292771e-03 At Newton iteration 4 norm is 1.009136944641e-08 ==== hyperelastic problem at iteration 15 ==== At Newton iteration 1 norm is 2.956483984656e+00 At Newton iteration 2 norm is 1.425542179283e+00 At Newton iteration 3 norm is 8.173309904226e-04 At Newton iteration 4 norm is 2.970461139668e-08 ==== hyperelastic problem at iteration 16 ==== At Newton iteration 1 norm is 4.534701443792e+00 At Newton iteration 2 norm is 1.739437285044e+00 At Newton iteration 3 norm is 2.017654905897e-03 At Newton iteration 4 norm is 3.487782468133e-07 ==== hyperelastic problem at iteration 17 ==== At Newton iteration 1 norm is 3.895960153877e+00 At Newton iteration 2 norm is 1.176152564012e+00 At Newton iteration 3 norm is 2.603624618812e-03 At Newton iteration 4 norm is 3.745578279839e-07 ==== hyperelastic problem at iteration 18 ==== At Newton iteration 1 norm is 3.384061489690e+00 At Newton iteration 2 norm is 7.223055977667e-01 At Newton iteration 3 norm is 1.445157758984e-03 At Newton iteration 4 norm is 2.164179163907e-07 ==== hyperelastic problem at iteration 19 ==== At Newton iteration 1 norm is 3.107703747451e+00 At Newton iteration 2 norm is 4.831377484870e-01 At Newton iteration 3 norm is 4.674511724038e-04 At Newton iteration 4 norm is 1.370129532958e-08 ==== hyperelastic problem at iteration 20 ==== At Newton iteration 1 norm is 2.400301309016e+00 At Newton iteration 2 norm is 2.376509080746e-01 At Newton iteration 3 norm is 1.550920289754e-04 At Newton iteration 4 norm is 1.141343804945e-09 ==== hyperelastic problem at iteration 21 ==== At Newton iteration 1 norm is 2.053735648196e+00 At Newton iteration 2 norm is 1.407127798207e-01 At Newton iteration 3 norm is 7.457761090080e-05 At Newton iteration 4 norm is 2.686492103262e-10 ==== hyperelastic problem at iteration 22 ==== At Newton iteration 1 norm is 1.648632934792e+00 At Newton iteration 2 norm is 8.248569820926e-02 At Newton iteration 3 norm is 3.370626762255e-05 At Newton iteration 4 norm is 4.296388011456e-11 ==== hyperelastic problem at iteration 23 ==== At Newton iteration 1 norm is 1.546937778131e+00 At Newton iteration 2 norm is 5.138908864712e-02 At Newton iteration 3 norm is 2.361859617092e-05 At Newton iteration 4 norm is 2.592412978914e-11 ==== hyperelastic problem at iteration 24 ==== At Newton iteration 1 norm is 1.246231374608e+00 At Newton iteration 2 norm is 3.678501987527e-02 At Newton iteration 3 norm is 1.255733187249e-05 At Newton iteration 4 norm is 5.570717372593e-12 ==== hyperelastic problem at iteration 25 ==== At Newton iteration 1 norm is 1.286138575822e+00 At Newton iteration 2 norm is 3.042189935187e-02 At Newton iteration 3 norm is 1.132113422438e-05 At Newton iteration 4 norm is 5.158499377960e-12 ==== hyperelastic problem at iteration 26 ==== At Newton iteration 1 norm is 1.052777305722e+00 At Newton iteration 2 norm is 2.242681949620e-02 At Newton iteration 3 norm is 5.717998007093e-06 ==== hyperelastic problem at iteration 27 ==== At Newton iteration 1 norm is 1.103903434659e+00 At Newton iteration 2 norm is 2.040534283204e-02 At Newton iteration 3 norm is 5.352548039303e-06 ==== hyperelastic problem at iteration 28 ==== At Newton iteration 1 norm is 9.204032940145e-01 At Newton iteration 2 norm is 1.582456168409e-02 At Newton iteration 3 norm is 2.934412240132e-06 ==== hyperelastic problem at iteration 29 ==== At Newton iteration 1 norm is 9.622663633986e-01 At Newton iteration 2 norm is 1.410278235113e-02 At Newton iteration 3 norm is 2.552972081916e-06 ==== hyperelastic problem at iteration 30 ==== At Newton iteration 1 norm is 8.118259463013e-01 At Newton iteration 2 norm is 1.158913505552e-02 At Newton iteration 3 norm is 1.548246047467e-06 ==== hyperelastic problem at iteration 31 ==== At Newton iteration 1 norm is 8.425934270652e-01 At Newton iteration 2 norm is 1.013723656949e-02 At Newton iteration 3 norm is 1.291763321732e-06 ==== hyperelastic problem at iteration 32 ==== At Newton iteration 1 norm is 7.204393057062e-01 At Newton iteration 2 norm is 8.477600764943e-03 At Newton iteration 3 norm is 8.184823617384e-07 ==== hyperelastic problem at iteration 33 ==== At Newton iteration 1 norm is 7.377972487665e-01 At Newton iteration 2 norm is 7.446515414223e-03 At Newton iteration 3 norm is 6.809515971478e-07 ==== hyperelastic problem at iteration 34 ==== At Newton iteration 1 norm is 6.434522139823e-01 At Newton iteration 2 norm is 6.292283094605e-03 At Newton iteration 3 norm is 4.474450917052e-07 ==== hyperelastic problem at iteration 35 ==== At Newton iteration 1 norm is 6.449716445435e-01 At Newton iteration 2 norm is 5.449121003665e-03 At Newton iteration 3 norm is 3.612338277685e-07 ==== hyperelastic problem at iteration 36 ==== At Newton iteration 1 norm is 5.770992218766e-01 At Newton iteration 2 norm is 4.779379557433e-03 At Newton iteration 3 norm is 2.549535272615e-07 ==== hyperelastic problem at iteration 37 ==== At Newton iteration 1 norm is 5.617416644725e-01 At Newton iteration 2 norm is 3.969814965775e-03 At Newton iteration 3 norm is 1.928132138684e-07 ==== hyperelastic problem at iteration 38 ==== At Newton iteration 1 norm is 4.527668882620e-01 At Newton iteration 2 norm is 2.784520189617e-03 At Newton iteration 3 norm is 8.695694588041e-08 ==== hyperelastic problem at iteration 39 ==== At Newton iteration 1 norm is 7.708134092579e+00 At Newton iteration 2 norm is 5.033524590616e-01 At Newton iteration 3 norm is 2.830309336251e-03 At Newton iteration 4 norm is 2.071644420258e-07 ==== hyperelastic problem at iteration 40 ==== At Newton iteration 1 norm is 8.173676068906e+00 At Newton iteration 2 norm is 1.514269342892e+00 At Newton iteration 3 norm is 2.365636812884e-02 At Newton iteration 4 norm is 1.874702797534e-05 At Newton iteration 5 norm is 1.766643067913e-11 ==== hyperelastic problem at iteration 41 ==== At Newton iteration 1 norm is 5.774661522726e+00 At Newton iteration 2 norm is 1.158168367760e+00 At Newton iteration 3 norm is 1.114663730375e-02 At Newton iteration 4 norm is 8.048880942239e-06 ==== hyperelastic problem at iteration 42 ==== At Newton iteration 1 norm is 3.973974310219e+00 At Newton iteration 2 norm is 9.279348715137e-01 At Newton iteration 3 norm is 3.643992370980e-03 At Newton iteration 4 norm is 1.241935718851e-06 ==== hyperelastic problem at iteration 43 ==== At Newton iteration 1 norm is 2.869960624545e+00 At Newton iteration 2 norm is 8.104736221393e-01 At Newton iteration 3 norm is 1.216077166730e-03 At Newton iteration 4 norm is 1.319117267436e-07 ==== hyperelastic problem at iteration 44 ==== At Newton iteration 1 norm is 2.345337641331e+00 At Newton iteration 2 norm is 4.978589802734e-01 At Newton iteration 3 norm is 6.040745200315e-04 At Newton iteration 4 norm is 3.227187856604e-08 ==== hyperelastic problem at iteration 45 ==== At Newton iteration 1 norm is 1.838126443274e+00 At Newton iteration 2 norm is 3.259453883896e-01 At Newton iteration 3 norm is 2.575388933912e-04 At Newton iteration 4 norm is 5.416774725427e-09 ==== hyperelastic problem at iteration 46 ==== At Newton iteration 1 norm is 1.629863319305e+00 At Newton iteration 2 norm is 2.062340319492e-01 At Newton iteration 3 norm is 1.496742930423e-04 At Newton iteration 4 norm is 1.962778689313e-09 ==== hyperelastic problem at iteration 47 ==== At Newton iteration 1 norm is 1.426411463519e+00 At Newton iteration 2 norm is 1.393893368187e-01 At Newton iteration 3 norm is 6.612363435009e-05 At Newton iteration 4 norm is 3.071596735232e-10 ==== hyperelastic problem at iteration 48 ==== At Newton iteration 1 norm is 1.271951144871e+00 At Newton iteration 2 norm is 8.031681655918e-02 At Newton iteration 3 norm is 3.270248688238e-05 At Newton iteration 4 norm is 1.033602779131e-10 ==== hyperelastic problem at iteration 49 ==== At Newton iteration 1 norm is 1.193299744566e+00 At Newton iteration 2 norm is 5.313935202362e-02 At Newton iteration 3 norm is 1.299311699457e-05 At Newton iteration 4 norm is 8.706744314871e-12 ==== hyperelastic problem at iteration 50 ==== At Newton iteration 1 norm is 1.061410714356e+00 At Newton iteration 2 norm is 2.989828088849e-02 At Newton iteration 3 norm is 5.729961988905e-06 ==== hyperelastic problem at iteration 51 ==== At Newton iteration 1 norm is 1.040774409029e+00 At Newton iteration 2 norm is 1.979703671510e-02 At Newton iteration 3 norm is 2.268300864415e-06 ==== hyperelastic problem at iteration 52 ==== At Newton iteration 1 norm is 9.369639467442e-01 At Newton iteration 2 norm is 1.078177232293e-02 At Newton iteration 3 norm is 1.144628723979e-06 ==== hyperelastic problem at iteration 53 ==== At Newton iteration 1 norm is 9.295003368522e-01 At Newton iteration 2 norm is 7.678300257699e-03 At Newton iteration 3 norm is 4.727286959404e-07 ==== hyperelastic problem at iteration 54 ==== At Newton iteration 1 norm is 8.512946167850e-01 At Newton iteration 2 norm is 4.149009982224e-03 At Newton iteration 3 norm is 2.360240551911e-07 ==== hyperelastic problem at iteration 55 ==== At Newton iteration 1 norm is 8.420677727751e-01 At Newton iteration 2 norm is 3.148091499936e-03 At Newton iteration 3 norm is 1.652857568067e-07 ==== hyperelastic problem at iteration 56 ==== At Newton iteration 1 norm is 7.821065018809e-01 At Newton iteration 2 norm is 1.674287465544e-03 At Newton iteration 3 norm is 4.582995657601e-08 ==== hyperelastic problem at iteration 57 ==== At Newton iteration 1 norm is 7.674758886732e-01 At Newton iteration 2 norm is 1.326396895043e-03 At Newton iteration 3 norm is 2.333914294164e-08 ==== hyperelastic problem at iteration 58 ==== At Newton iteration 1 norm is 7.197511314937e-01 At Newton iteration 2 norm is 7.936822866902e-04 At Newton iteration 3 norm is 4.599855931125e-08 ==== hyperelastic problem at iteration 59 ==== At Newton iteration 1 norm is 6.996307456258e-01 At Newton iteration 2 norm is 5.905782810761e-04 At Newton iteration 3 norm is 1.195061886943e-08 ==== hyperelastic problem at iteration 60 ==== At Newton iteration 1 norm is 6.606632796233e-01 At Newton iteration 2 norm is 3.850977980078e-04 At Newton iteration 3 norm is 1.040892561972e-08 ==== hyperelastic problem at iteration 61 ==== At Newton iteration 1 norm is 6.371348335681e-01 At Newton iteration 2 norm is 2.936729195034e-04 At Newton iteration 3 norm is 1.938337284431e-09 ==== hyperelastic problem at iteration 62 ==== At Newton iteration 1 norm is 6.038981115209e-01 At Newton iteration 2 norm is 1.953221329294e-04 At Newton iteration 3 norm is 1.797402056268e-09 ==== hyperelastic problem at iteration 63 ==== At Newton iteration 1 norm is 5.790003413761e-01 At Newton iteration 2 norm is 1.695571680352e-04 At Newton iteration 3 norm is 5.815421539501e-10 ==== hyperelastic problem at iteration 64 ==== At Newton iteration 1 norm is 5.494457364906e-01 At Newton iteration 2 norm is 1.111039888199e-04 At Newton iteration 3 norm is 4.359748674653e-10 ==== hyperelastic problem at iteration 65 ==== At Newton iteration 1 norm is 5.245877433492e-01 At Newton iteration 2 norm is 1.099736070421e-04 At Newton iteration 3 norm is 2.654433659709e-10 ==== hyperelastic problem at iteration 66 ==== At Newton iteration 1 norm is 4.978052937194e-01 At Newton iteration 2 norm is 7.731081074140e-05 At Newton iteration 3 norm is 1.466760945309e-10 ==== hyperelastic problem at iteration 67 ==== At Newton iteration 1 norm is 4.738342956316e-01 At Newton iteration 2 norm is 8.336401844557e-05 At Newton iteration 3 norm is 1.513581218539e-10 ==== hyperelastic problem at iteration 68 ==== At Newton iteration 1 norm is 4.494083154957e-01 At Newton iteration 2 norm is 6.507252811327e-05 At Newton iteration 3 norm is 8.325370982536e-11 ==== hyperelastic problem at iteration 69 ==== At Newton iteration 1 norm is 4.268000319305e-01 At Newton iteration 2 norm is 7.214130941918e-05 At Newton iteration 3 norm is 9.730749796639e-11 ==== hyperelastic problem at iteration 70 ==== At Newton iteration 1 norm is 4.045434294606e-01 At Newton iteration 2 norm is 6.177243363777e-05 At Newton iteration 3 norm is 6.269236853081e-11 ==== hyperelastic problem at iteration 71 ==== At Newton iteration 1 norm is 3.836002435996e-01 At Newton iteration 2 norm is 6.619761620414e-05 At Newton iteration 3 norm is 6.754362195686e-11 ==== hyperelastic problem at iteration 72 ==== At Newton iteration 1 norm is 3.634479925583e-01 At Newton iteration 2 norm is 5.962238272375e-05 At Newton iteration 3 norm is 4.908986991402e-11 ==== hyperelastic problem at iteration 73 ==== At Newton iteration 1 norm is 3.443744681492e-01 At Newton iteration 2 norm is 6.089069793544e-05 At Newton iteration 3 norm is 4.968725596045e-11 ==== hyperelastic problem at iteration 74 ==== At Newton iteration 1 norm is 3.263339194711e-01 At Newton iteration 2 norm is 5.569283896651e-05 At Newton iteration 3 norm is 3.897789461268e-11 ==== hyperelastic problem at iteration 75 ==== At Newton iteration 1 norm is 3.093031037774e-01 At Newton iteration 2 norm is 5.466622157020e-05 At Newton iteration 3 norm is 3.797778261895e-11 ==== hyperelastic problem at iteration 76 ==== At Newton iteration 1 norm is 2.934514448583e-01 At Newton iteration 2 norm is 4.993056198137e-05 At Newton iteration 3 norm is 3.148435361632e-11 ==== hyperelastic problem at iteration 77 ==== At Newton iteration 1 norm is 2.786423447758e-01 At Newton iteration 2 norm is 4.743938503547e-05 At Newton iteration 3 norm is 3.037432798384e-11 ==== hyperelastic problem at iteration 78 ==== At Newton iteration 1 norm is 2.651156642838e-01 At Newton iteration 2 norm is 4.300524173056e-05 At Newton iteration 3 norm is 2.626562634164e-11 ==== hyperelastic problem at iteration 79 ==== At Newton iteration 1 norm is 2.527239189366e-01 At Newton iteration 2 norm is 3.985779202077e-05 At Newton iteration 3 norm is 2.517362514792e-11 ==== hyperelastic problem at iteration 80 ==== At Newton iteration 1 norm is 2.417092421370e-01 At Newton iteration 2 norm is 3.581933462041e-05 At Newton iteration 3 norm is 2.220310729555e-11 Saved mesh data in: /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v3.1.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.1.0/examples/CardiacElastoDynamics/results/step1/ellipsoid0.03_dt0.0125_analytic_rigi Generate movie 01/20 (5.00 %) 2.81 s Generate movie 02/20 (10.00 %) 2.39 s Generate movie 03/20 (15.00 %) 2.23 s Generate movie 04/20 (20.00 %) 2.09 s Generate movie 05/20 (25.00 %) 1.96 s Generate movie 06/20 (30.00 %) 1.76 s Generate movie 07/20 (35.00 %) 1.65 s Generate movie 08/20 (40.00 %) 1.48 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.15 s Generate movie 12/20 (60.00 %) 1.04 s Generate movie 13/20 (65.00 %) 918.82 ms Generate movie 14/20 (70.00 %) 787.71 ms Generate movie 15/20 (75.00 %) 665.98 ms Generate movie 16/20 (80.00 %) 531.73 ms Generate movie 17/20 (85.00 %) 390.92 ms Generate movie 18/20 (90.00 %) 263.55 ms Generate movie 19/20 (95.00 %) 130.45 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 utils import ( RESULTS_DIR, DATA_DIR, Get_config_ellipsoid, Get_stresses, Get_pressures, ) class CardiacElastoDynamics(Simulations.HyperElastic): def __init__( self, mesh, model, folder="", tolConv=0.00001, maxIter=20, verbosity=False, alpha_top=1e5, alpha_epi=1e8, beta_top=5e3, beta_epi=5e3, ): super().__init__(mesh, model, folder, tolConv, maxIter, verbosity) self.__dict_pressure: dict[str, float] = {} self.__alpha_top = alpha_top self.__alpha_epi = alpha_epi self.__beta_top = beta_top self.__beta_epi = beta_epi def Set_pressure(self, dict_pressure: dict[str, float]): self.__dict_pressure = dict_pressure 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. Kpressure_e, Rpressure_e = 0.0, 0.0 for tag, pressure in self.__dict_pressure.items(): if tag in groupElem.elementTags: tangent_e, residual_e = Operators.NonLinear.FollowingPressure( groupElem, displacement, pressure, groupElem.Get_Elements_Tag(tag), MatrixType.mass, ) Kpressure_e += tangent_e Rpressure_e += residual_e # 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") Ktop_e[top_e] = self.__alpha_top * M_e[top_e] Ctop_e[top_e] = self.__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") Kepi_e[epi_e] = self.__alpha_epi * Ms_e[epi_e] Cepi_e[epi_e] = self.__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] = ( Kpressure_e + K_penalty_e, Ctop_e + Cepi_e, None, Rpressure_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 meshName = "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, meshName) 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 times = np.linspace(0, 1, Nt + 1) dt = times[1] - times[0] stresses = Get_stresses(times) pressures = Get_pressures(times) 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(times, stresses) 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(times, pressures) name = "pressure" Matplotlib.Save_fig(results_dir, name) if config is Config.step0B: stresses *= 0 if config is Config.step0A: pressures *= 0 if doSimu: # ---------------------------------------------- # Mesh, fibers and sheets # ---------------------------------------------- mesh, fibers_e_pg, sheets_e_pg = Get_config_ellipsoid( Folder.Join(DATA_DIR, meshName), 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 times: simu.Bc_Init() simu.pressure = np.interp(t + dt / 2, times, pressures) material.active_stress = np.interp(t + dt / 2, times, stresses) 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 = times[: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 11.881 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 `_