Define materials

This example shows how to create a mechanical material and assign it to a heart part.

Import material module

from pathlib import Path

import matplotlib.pyplot as plt

from ansys.health.heart.settings.material.curve import (
    ActiveCurve,
    constant_ca2,
    kumaraswamy_active,
)
from ansys.health.heart.settings.material.ep_material import CellModel, EPMaterial
from ansys.health.heart.settings.material.material import (
    ACTIVE,
    ANISO,
    ISO,
    ActiveModel,
    Mat295,
    NeoHookean,
)

Note

The unit system for heart modeling in LS-DYNA is ["MPa", "mm", "N", "ms", "g"].

Create a material

Create a Neo-Hookean material as follows.

neo = NeoHookean(rho=0.001, c10=1, nu=0.499)

C:\Users\ansys\actions-runner\_work\pyansys-heart\pyansys-heart\examples\preprocessor\demo-material_pr.py:61: DeprecationWarning: Call to deprecated class NeoHookean. (Use *MAT_295 with the ISO module instead.)
  neo = NeoHookean(rho=0.001, c10=1, nu=0.499)

The recommended way to create a Neo-Hookean material is by activating only the isotropic module in MAT_295.

neo2 = Mat295(rho=0.001, iso=ISO(itype=1, beta=2, kappa=1, mu1=0.05, alpha1=2))

Note

For more information on MAT_295, see the LS-DYNA manuals.

# Additional steps follow for creating MAT_295, which is used for myocardium.

# step 1: create an isotropic module
iso = ISO(k1=1, k2=1, kappa=100)

# step 2: create an anisotropic module
fiber = ANISO.HGOFiber(k1=1, k2=1)
aniso1 = ANISO(fibers=[fiber])

# Create fiber with sheet, and their interactions
sheet = ANISO.HGOFiber(k1=1, k2=1)
aniso2 = ANISO(fibers=[fiber, sheet], k1fs=1, k2fs=1)

# step3: create the active module

# example 1:
# create active model 1
ac_model1 = ActiveModel.Model1()
# create Ca2+ curve
ac_curve1 = ActiveCurve(constant_ca2(tb=800, ca2ionm=ac_model1.ca2ionm), type="ca2", threshold=0.5)
# build active module
active = ACTIVE(model=ac_model1, ca2_curve=ac_curve1)

## Active model 1 must have a constant ca2ion,
# but the curve must cross the threshold at every start of the heart beat.

Plot Ca2+ with threshold

Plot Ca2+ with the threshold.

fig = active.ca2_curve.plot_time_vs_ca2()
plt.show()

# example 2
# create active model 3
ac_model3 = ActiveModel.Model3()
# create a stress curve and show
ac_curve3 = ActiveCurve(kumaraswamy_active(t_end=800), type="stress")
fig = ac_curve3.plot_time_vs_stress()
plt.show()

• 
• 

Note

When eta=0 in model 3, the stress curve is the active stress for all elements. If eta!=0, this is the idealized active stress when fiber stretch stays at 1.

# PyAnsys Heart converts the stress curve to Ca2+ curve (input of MAT_295)
fig = ac_curve3.plot_time_vs_ca2()
plt.show()

# build active module
active3 = ACTIVE(model=ac_model3, ca2_curve=ac_curve3)

Create MAT_295 with modules

Create MAT_295 with the preceding modules.

iso_mat = Mat295(rho=1, iso=iso, aniso=None, active=None)
passive_mat = Mat295(rho=1, iso=iso, aniso=aniso1, active=None)
active_mat = Mat295(rho=1, iso=iso, aniso=aniso1, active=active)

Create EP materials

ep_mat_active = EPMaterial.Active(
    sigma_fiber=1, sigma_sheet=0.5, beta=140, cm=0.01, cell_model=CellModel.Tentusscher()
)
epinsulator = EPMaterial.Insulator()

# import pyvista as pv


# ##############################################################################
# .. note::
#    The Ca2+ curve is ignored if the simulation is coupled with electrophysiology.

Assign materials to a part

Assign the materials to the heart model.

Load a heart model

Note

You must complete the full heart preprocessing example first.

import ansys.health.heart.examples as examples
import ansys.health.heart.models as models

heart_model_vtu, heart_model_partinfo, _ = examples.get_preprocessed_fullheart()
workdir = str(Path.home() / "pyansys-heart" / "Rodero2021")

# Load a full-heart model.
heartmodel: models.FullHeart = models.FullHeart(working_directory=workdir)
heartmodel.load_model_from_mesh(heart_model_vtu, heart_model_partinfo)

heartmodel.plot_mesh(show_edges=False)

# Print the default material. You should see that the material is empty.
print(heartmodel.left_ventricle.meca_material)
print(heartmodel.left_ventricle.ep_material)

Static Scene

Interactive Scene

Material is empty.
Material is empty.

Note

If no material is set before writing k files, the default material from the settings object is used.

# Assign the material that you just created.
heartmodel.left_ventricle.meca_material = active_mat
heartmodel.left_ventricle.ep_material = ep_mat_active

# Print it. You should see the following:
# MAT295(rho=1, iso=ISO(itype=-3, beta=0.0, nu=0.499, k1=1, k2=1), aopt=2.0, aniso=ANISO(atype=-1, fibers=[ANISO.HGOFiber(k1=1, k2=1, a=0.0, b=1.0, _theta=0.0, _ftype=1, _fcid=0)], k1fs=None, k2fs=None, vec_a=(1.0, 0.0, 0.0), vec_d=(0.0, 1.0, 0.0), nf=1, intype=0), active=ActiveModel.Model1(t0=None, ca2ion=None, ca2ionm=4.35, n=2, taumax=0.125, stf=0.0, b=4.75, l0=1.58, l=1.85, dtmax=150, mr=1048.9, tr=-1629.0))  # noqa
print(heartmodel.left_ventricle.meca_material)
print(heartmodel.left_ventricle.ep_material)

Mat295(rho=1, iso=ISO(itype=-3, beta=0.0, nu=0.4950166112956811, k1=1, k2=1, mu1=None, alpha1=None, kappa=100), aopt=2.0, aniso=ANISO(atype=-1, fibers=[ANISO.HGOFiber(k1=1, k2=1, a=0.0, b=1.0, _theta=0.0, _ftype=1, _fcid=0)], k1fs=None, k2fs=None, vec_a=(1.0, 0.0, 0.0), vec_d=(0.0, 1.0, 0.0), nf=1, intype=0), active=ACTIVE(acid=None, actype=1, acthr=0.5, acdir=1, sf=1.0, ss=0.0, sn=0.0, model=ActiveModel.Model1(t0=None, ca2ion=None, ca2ionm=4.35, n=2, taumax=0.125, stf=0.0, b=4.75, l0=1.58, l=1.85, dtmax=150, mr=1048.9, tr=-1629.0), ca2_curve=<ansys.health.heart.settings.material.curve.ActiveCurve object at 0x0000023D90A92890>))
EPMaterial.Active(sigma_fiber=1, sigma_sheet=0.5, sigma_sheet_normal=0.1, beta=140, cm=0.01, lambda_=None, cell_model=CellModel.Tentusscher(gas_constant=8314.472, t=310, faraday_constant=96485.3415, cm=0.185, vc=0.016404, vsr=0.001094, vss=5.468e-05, pkna=0.03, ko=5.4, nao=140.0, cao=2.0, gk1=5.405, gkr=0.153, gna=14.838, gbna=0.0002, gcal=3.98e-05, gbca=0.000592, gpca=0.1238, gpk=0.0146, pnak=2.724, km=1.0, kmna=40.0, knaca=1000.0, ksat=0.1, alpha=2.5, gamma=0.35, kmca=1.38, kmnai=87.5, kpca=0.0005, k1=0.15, k2=0.045, k3=0.06, k4=0.005, ec=1.5, maxsr=2.5, minsr=1.0, vrel=0.102, vleak=0.00036, vxfer=0.0038, vmaxup=0.006375, kup=0.00025, bufc=0.2, kbufc=0.001, bufsr=10.0, kbufsf=0.3, bufss=0.4, kbufss=0.00025, gks=0.392, gto=0.294, v=-85.23, ki=136.89, nai=8.604, cai=0.000126, cass=0.00036, casr=3.64, rpri=0.9073, xr1=0.00621, xr2=0.4712, xs=0.0095, m=0.00172, h=0.7444, j=0.7045, d=3.373e-05, f=0.7888, f2=0.9755, fcass=0.9953, s=0.999998, r=2.42e-08))

Create a new part and set material

# # A new part can be created by elements IDs
# ids = np.where(heartmodel.mesh.point_data_to_cell_data()["uvc_longitudinal"] > 0.9)[0]
# new_part: Part = heartmodel.create_part_by_ids(ids, "new_part")

# # Show the part
# plotter = heartmodel.plot_part(new_part)