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 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

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()