Note
Go to the end to download the full example code.
MonoVentricular#
Passive + active hyperelastic simulation of an ellipsoidal left-ventricle model.
Combines a Holzapfel-Ogden orthotropic law (fiber + sheet directions), an active stress along the fiber direction, and a 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 <N> python MonoVentricular.py with a PETSc-backed solver), or, on a single process, with the pypardiso solver — both markedly cut the solve time. The default useCoarseConfig=True is light enough to run as-is.
===== hyperelastic problem at iteration 0 =====
At Newton iteration 1 norm is 8.686270292626e-02
At Newton iteration 2 norm is 2.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
21 from enum import Enum
22
23 import numpy as np
24
25 from EasyFEA import (
26 Terminal,
27 Matplotlib,
28 Folder,
29 PyVista,
30 MatrixType,
31 Models,
32 Simulations,
33 AlgoType,
34 )
35 from EasyFEA.FEM import Operators
36 from EasyFEA.Utilities import _params
37
38 from utils import RESULTS_DIR, DATA_DIR, Get_config, Get_values
39
40
41 class CardiacElastoDynamics(Simulations.HyperElastic):
42
43 pressure = _params.ScalarParameter()
44
45 def __init__(
46 self, mesh, model, folder="", tolConv=0.00001, maxIter=20, verbosity=False
47 ):
48 super().__init__(mesh, model, folder, tolConv, maxIter, verbosity)
49 self.pressure = 0.0
50
51 def Construct_local_matrix_system(self, problemType):
52
53 assert isinstance(self.material, Models.HyperElastic.HolzapfelOgden)
54 nPg = self.material.T1.shape[1]
55
56 results = super().Construct_local_matrix_system(problemType, nPg)
57
58 # current Newton-Raphson iterate (updated via u += delta_u)
59 displacement = self._Solver_Get_Newton_Raphson_current_solution()
60 if self.algo in AlgoType.Get_Hyperbolic_Types():
61 displacement, _, _ = self._Solver_Evaluate_u_v_a_for_time_scheme(
62 problemType, displacement
63 )
64
65 for groupElem in self.mesh.Get_list_groupElem(self.dim - 1):
66
67 # Following pressure — tracks the deformed normal, so it depends on
68 # the current iterate / pressure and is rebuilt every Newton step.
69 tangent_e, residual_e = Operators.NonLinear.FollowingPressure(
70 groupElem,
71 displacement,
72 self.pressure,
73 groupElem.Get_Elements_Tag("endo"),
74 MatrixType.mass,
75 )
76
77 # top — isotropic surface mass penalty (Robin α·u + β·u̇ = 0)
78 M_e = Operators.Bilinear.UV(groupElem, dof_n=3)
79
80 Ktop_e = np.zeros_like(M_e)
81 Ctop_e = np.zeros_like(M_e)
82 if "top" in groupElem.elementTags:
83 top_e = groupElem.Get_Elements_Tag("top")
84 alpha_top = 1e5
85 Ktop_e[top_e] = alpha_top * M_e[top_e]
86
87 beta_top = 5e3
88 Ctop_e[top_e] = beta_top * M_e[top_e]
89
90 # epi — normal-direction mass penalty (Robin α·(u·n̂) + β·(u̇·n̂) = 0)
91 Ms_e = Operators.Bilinear.MassAlongNormal(groupElem)
92
93 Kepi_e = np.zeros_like(Ms_e)
94 Cepi_e = np.zeros_like(Ms_e)
95 if "epi" in groupElem.elementTags:
96 epi_e = groupElem.Get_Elements_Tag("epi")
97 alpha_epi = 1e8
98 Kepi_e[epi_e] = alpha_epi * Ms_e[epi_e]
99
100 beta_epi = 5e3
101 Cepi_e[epi_e] = beta_epi * Ms_e[epi_e]
102
103 # Penalty residual contribution: −K_penalty · u_t at current iterate
104 K_penalty_e = Ktop_e + Kepi_e
105 assembly_e = groupElem.Get_assembly_e(self.dim)
106 u_e = displacement[assembly_e] # (Ne_surf, nPe·3)
107 f_penalty_e = np.einsum("eij,ej->ei", K_penalty_e, u_e)
108
109 results[groupElem] = (
110 tangent_e + K_penalty_e,
111 Ctop_e + Cepi_e,
112 None,
113 residual_e - f_penalty_e,
114 )
115
116 return results
117
118
119 class Config(str, Enum):
120 step0A = "step0A" # active_stress
121 step0B = "step0B" # pressure
122 step1 = "step1" # active_stress + pressure
123
124
125 if __name__ == "__main__":
126
127 Terminal.Clear()
128
129 # ----------------------------------------------
130 # Config
131 # ----------------------------------------------
132
133 useCoarseConfig = True
134
135 ellipsoid = "ellipsoid0.03" if useCoarseConfig else "ellipsoid0.005"
136
137 config = Config.step1
138
139 fiberSource = "analytic"
140 # fiberSource = "vtu"
141
142 matrixType = MatrixType.rigi
143 # matrixType = MatrixType.mass
144 # matrixType = 15
145
146 results_dir = Folder.Join(RESULTS_DIR, config.name, ellipsoid)
147
148 doSimu = True
149 plotGraph = False
150 plotParticles = True
151 saveParticles = True
152 makeMovie = True
153
154 # ----------------------------------------------
155 # time-history needed for plotting in both doSimu / Load_Simu flows
156
157 Nt = 80 if useCoarseConfig else 1000
158
159 t_values, activeStress_values, pressure_values = Get_values(Tmax=1.0, Nt=Nt)
160 dt = t_values[1] - t_values[0]
161 results_dir += f"_dt{dt}_{fiberSource}_{matrixType}"
162
163 if plotGraph:
164 ax = Matplotlib.Init_Axes()
165 ax.grid()
166 ax.set_xlabel(r"$t$ [s]")
167 ax.set_ylabel(r"$\tau(t)$ [Pa]")
168 ax.plot(t_values, activeStress_values)
169 name = "active_pressure"
170 Matplotlib.Save_fig(results_dir, name)
171
172 ax = Matplotlib.Init_Axes()
173 ax.grid()
174 ax.set_xlabel(r"$t$ [s]")
175 ax.set_ylabel(r"$p(t)$ [Pa]")
176 ax.plot(t_values, pressure_values)
177 name = "pressure"
178 Matplotlib.Save_fig(results_dir, name)
179
180 if config is Config.step0B:
181 activeStress_values *= 0
182 if config is Config.step0A:
183 pressure_values *= 0
184
185 if doSimu:
186
187 # ----------------------------------------------
188 # Mesh, fibers and sheets
189 # ----------------------------------------------
190
191 mesh, fibers_e_pg, sheets_e_pg = Get_config(
192 Folder.Join(DATA_DIR, ellipsoid),
193 matrixType=matrixType,
194 fiberSource=fiberSource,
195 plotMesh=False,
196 plotTags=False,
197 plotFibers=False,
198 )
199
200 # ----------------------------------------------
201 # Material
202 # ----------------------------------------------
203
204 # solid
205 a = 59.0
206 a_f = 18472.0
207 a_fs = 216.0
208 a_s = 2481.0
209 b = 8.023
210 b_f = 16.026
211 b_fs = 11.436
212 b_s = 11.12
213
214 material = Models.HyperElastic.HolzapfelOgden(
215 dim=3,
216 C0=a / 2 / b,
217 C1=b,
218 C2=a_f / 2 / b_f,
219 C3=b_f,
220 C4=a_s / 2 / b_s,
221 C5=b_s,
222 C6=a_fs / 2 / b_fs,
223 C7=b_fs,
224 K=1e6,
225 Mu1=0.0,
226 Mu2=0.0,
227 T1=fibers_e_pg,
228 T2=sheets_e_pg,
229 ks=100,
230 )
231 material.eta = 100.0
232 material.Set_active_stress_vec(material.T1)
233
234 # ----------------------------------------------
235 # Simulation
236 # ----------------------------------------------
237
238 simu = CardiacElastoDynamics(mesh, material, folder=results_dir)
239
240 simu.Solver_Set_Hyperbolic_Algorithm(dt, algo=AlgoType.midpoint)
241 simu.rho = 1000
242
243 for t in t_values:
244 simu.Bc_Init()
245 simu.pressure = np.interp(t + dt / 2, t_values, pressure_values)
246 material.active_stress = np.interp(
247 t + dt / 2, t_values, activeStress_values
248 )
249 simu.Solve()
250 simu.Save_Iter()
251
252 simu.Save(results_dir)
253
254 else:
255 simu = Simulations.Load_Simu(results_dir)
256
257 simu._Gather()
258
259 if plotParticles and simu.isGathered:
260
261 coords = [(0.025, 0.03, 0), (0, 0.03, 0)]
262 evalCoords = np.array(coords)
263 evalElements = simu.mesh.groupElem._Get_nearby_elements(evalCoords)
264
265 Niter = simu.Niter
266 values = np.empty((Niter, len(coords), 3))
267 for i in range(Niter):
268 simu.Set_Iter(i)
269 values[i] = simu.mesh.Evaluate_dofsValues_at_coordinates(
270 evalCoords, simu.displacement, evalElements
271 )
272
273 times = t_values[:Niter]
274 axs = Matplotlib.plt.subplots(3, 2, sharex=True)[1]
275
276 for p, (particle, coord) in enumerate(zip(["p0", "p1"], coords)):
277
278 for c, component in enumerate(["x", "y", "z"]):
279
280 ax: Matplotlib.plt.Axes = axs[c, p]
281
282 ax.grid()
283
284 if c == 2:
285 ax.set_xlabel("Time [s]")
286 if p == 0:
287 ax.set_ylabel(f"Displacement {component}-component [m]")
288 if c == 0:
289 ax.set_title(f"Particle {particle}")
290
291 ax.plot(times, values[:, p, c])
292
293 width, height = ax.figure.get_size_inches()
294 ax.figure.set_size_inches(width * 1.5, height * 2.5)
295 Matplotlib.Save_fig(results_dir, "particles")
296
297 if saveParticles and simu.isGathered:
298
299 # per-iteration deformed volume
300 volumes = np.empty(Niter)
301 for i in range(Niter):
302 simu.Set_Iter(i)
303 deformed = simu.mesh.copy()
304 deformed.coord += simu.displacement.reshape(-1, 3)
305 volumes[i] = deformed.volume
306
307 dict_particles = {
308 "time": times,
309 "displacement": {
310 f"p{p}": {
311 "ux": values[:, p, 0],
312 "uy": values[:, p, 1],
313 "uz": values[:, p, 2],
314 "magnitude": np.linalg.norm(values[:, p, :], axis=1),
315 }
316 for p in range(2)
317 },
318 "stress": {
319 "time": None,
320 "p0": {"magnitude": None},
321 "p1": {"magnitude": None},
322 },
323 "volume": volumes,
324 }
325 Simulations.Save_pickle(dict_particles, results_dir, "particles")
326
327 if makeMovie:
328 values = [simu.Result("ux", iter=i) for i in range(simu.Niter)]
329 clim = (np.min(values), np.max(values))
330 PyVista.Movie_simu(
331 simu,
332 "ux",
333 results_dir,
334 "ux.gif",
335 N=20,
336 deformFactor=1.0,
337 clim=clim,
338 plotMesh=True,
339 )
340
341 Matplotlib.plt.show()
Total running time of the script: (0 minutes 12.879 seconds)