from pathlib import Path
import dolfinx
import numpy as np
from ..utils import space_from_string
from . import utils
[docs]
def mu_theta(
x: np.ndarray,
y: np.ndarray,
z: np.ndarray,
rs: np.ndarray,
rl: np.ndarray,
long_axis: int = 0,
) -> tuple[np.ndarray, np.ndarray, list[int]]:
"""Get the angles mu and theta from the coordinates x, y, z
given the long axis.
Parameters
----------
x : np.ndarray
The x-coordinates
y : np.ndarray
The y-coordinates
z : np.ndarray
The z-coordinates
rs : np.ndarray
The short radius
rl : np.ndarray
The long radius
long_axis : int, optional
The long axis, by default 0 (x-axis)
Returns
-------
tuple[np.ndarray, np.ndarray, list[int]]
The angles mu and theta and the permutation of the axes
Raises
------
ValueError
If the long axis is not 0, 1 or 2
"""
if long_axis == 0:
a = np.sqrt(y**2 + z**2) / rs
b = x / rl
theta = np.pi - np.arctan2(z, -y)
perm = [0, 1, 2]
elif long_axis == 1:
a = np.sqrt(x**2 + z**2) / rs
b = y / rl
theta = np.pi - np.arctan2(z, -x)
perm = [1, 0, 2]
elif long_axis == 2:
a = np.sqrt(x**2 + y**2) / rs
b = z / rl
theta = np.pi - np.arctan2(x, -y)
perm = [2, 1, 0]
else:
raise ValueError("Invalid long_axis")
mu = np.arctan2(a, b)
theta[mu < 1e-7] = 0.0
return mu, theta, perm
[docs]
def transmural_coordinate(
V: dolfinx.fem.FunctionSpace,
r_short_endo=0.025,
r_short_epi=0.035,
r_long_endo=0.09,
r_long_epi=0.097,
long_axis: int = 0,
num_iterations: int = 60,
) -> dolfinx.fem.Function:
"""The transmural coordinate t in [0, 1] of the ellipsoid, as defined in
the benchmark of Land et al. (2015): the point lies on the ellipsoid with
radii r_short(t) and r_long(t), interpolated linearly between the
endocardium (t = 0) and the epicardium (t = 1).
The root is found by bisection, as
(r^2 - x_long^2) / r_short(t)^2 + x_long^2 / r_long(t)^2 - 1
decreases with t. Points just outside the wall get t = 0 or t = 1.
Parameters
----------
V : dolfinx.fem.FunctionSpace
Scalar function space of the transmural coordinate
r_short_endo : float, optional
Short radius at the endocardium, by default 0.025
r_short_epi : float, optional
Short radius at the epicardium, by default 0.035
r_long_endo : float, optional
Long radius at the endocardium, by default 0.09
r_long_epi : float, optional
Long radius at the epicardium, by default 0.097
long_axis : int, optional
Long axis, by default 0 (x-axis)
num_iterations : int, optional
Number of bisection steps, by default 60
Returns
-------
dolfinx.fem.Function
The transmural coordinate
"""
x = V.tabulate_dof_coordinates()
b2 = x[:, long_axis] ** 2
a2 = (x**2).sum(axis=1) - b2
def g(t):
rs = r_short_endo + (r_short_epi - r_short_endo) * t
rl = r_long_endo + (r_long_epi - r_long_endo) * t
return a2 / rs**2 + b2 / rl**2 - 1.0
lo = np.zeros_like(b2)
hi = np.ones_like(b2)
for _ in range(num_iterations):
mid = 0.5 * (lo + hi)
outside = g(mid) > 0.0
lo = np.where(outside, mid, lo)
hi = np.where(outside, hi, mid)
t = dolfinx.fem.Function(V, name="t")
t.x.array[:] = 0.5 * (lo + hi)
return t
[docs]
def compute_system(
t_func: dolfinx.fem.Function,
r_short_endo=0.025,
r_short_epi=0.035,
r_long_endo=0.09,
r_long_epi=0.097,
alpha_endo: float = -60,
alpha_epi: float = 60,
long_axis: int = 0,
**kwargs,
) -> utils.Microstructure:
"""Compute the microstructure for the given time function.
Parameters
----------
t_func : dolfinx.fem.Function
Solution of the Laplace equation
r_short_endo : float, optional
Short radius at the endocardium, by default 0.025
r_short_epi : float, optional
Short radius at the epicardium, by default 0.035
r_long_endo : float, optional
Long radius at the endocardium, by default 0.09
r_long_epi : float, optional
Long radius at the epicardium, by default 0.097
alpha_endo : float, optional
Angle at the endocardium, by default -60
alpha_epi : float, optional
Angle at the epicardium, by default 60
long_axis : int, optional
Long axis, by default 0 (x-axis)
Returns
-------
utils.Microstructure
The microstructure
"""
V = t_func.function_space
element = V.ufl_element()
mesh = V.mesh
dof_coordinates = V.tabulate_dof_coordinates()
alpha = lambda x: (alpha_endo + (alpha_epi - alpha_endo) * x) * (np.pi / 180)
r_long = lambda x: r_long_endo + (r_long_epi - r_long_endo) * x
r_short = lambda x: r_short_endo + (r_short_epi - r_short_endo) * x
drl_dt = r_long_epi - r_long_endo
drs_dt = r_short_epi - r_short_endo
t = t_func.x.array
rl = r_long(t)
rs = r_short(t)
al = alpha(t)
x = dof_coordinates[:, 0]
y = dof_coordinates[:, 1]
z = dof_coordinates[:, 2]
mu, theta, perm = mu_theta(x, y, z, rs, rl, long_axis=long_axis)
e_t = np.array(
[
drl_dt * np.cos(mu),
drs_dt * np.sin(mu) * np.cos(theta),
drs_dt * np.sin(mu) * np.sin(theta),
],
)[perm]
e_t = utils.normalize(e_t)
e_mu = np.array(
[
-rl * np.sin(mu),
rs * np.cos(mu) * np.cos(theta),
rs * np.cos(mu) * np.sin(theta),
],
)[perm]
e_mu = utils.normalize(e_mu)
e_theta = np.array(
[
np.zeros_like(t),
-rs * np.sin(mu) * np.sin(theta),
rs * np.sin(mu) * np.cos(theta),
],
)[perm]
e_theta = utils.normalize(e_theta)
f0 = np.sin(al) * e_mu + np.cos(al) * e_theta
f0 = utils.normalize(f0)
n0 = np.cross(e_mu, e_theta, axis=0)
n0 = utils.normalize(n0)
s0 = np.cross(f0, n0, axis=0)
s0 = utils.normalize(s0)
# The fibers have three components also on the meridian of an axisymmetric
# geometry, where the third one is the circumferential direction
Vv = space_from_string(
space_string=f"{element.family_name}_{element.degree}",
mesh=mesh,
dim=3,
discontinuous=element.discontinuous,
)
fiber = dolfinx.fem.Function(Vv)
norm_f = np.linalg.norm(f0, axis=0)
fiber.x.array[:] = (f0 / norm_f).T.reshape(-1)
fiber.name = "f0"
sheet = dolfinx.fem.Function(Vv)
sheet.x.array[:] = s0.T.reshape(-1)
sheet.name = "s0"
sheet_normal = dolfinx.fem.Function(Vv)
sheet_normal.x.array[:] = n0.T.reshape(-1)
sheet_normal.name = "n0"
return utils.Microstructure(f0=fiber, s0=sheet, n0=sheet_normal)
[docs]
def create_microstructure(
mesh,
ffun,
markers,
function_space="P_1",
r_short_endo=0.025,
r_short_epi=0.035,
r_long_endo=0.09,
r_long_epi=0.097,
alpha_endo: float = -60,
alpha_epi: float = 60,
long_axis: int = 0,
outdir: str | Path | None = None,
transmural: str = "laplace",
**kwargs,
):
"""Create the microstructure of an LV ellipsoid.
`transmural` selects the transmural coordinate the fiber angle varies
linearly in: ``"laplace"``, the solution of a Laplace problem between ENDO
and EPI, or ``"analytic"``, the coordinate of the ellipsoid itself (see
`transmural_coordinate`), which is the one of the benchmark of Land et al.
(2015). Use ``"analytic"`` on the meridian of an axisymmetric geometry,
where a planar Laplace problem is not the axisymmetric one.
"""
if transmural == "laplace":
t = utils.laplace(
mesh,
ffun,
endo_marker=markers["ENDO"][0],
epi_marker=markers["EPI"][0],
function_space=function_space,
)
elif transmural == "analytic":
t = transmural_coordinate(
space_from_string(function_space, mesh, dim=1),
r_short_endo=r_short_endo,
r_short_epi=r_short_epi,
r_long_endo=r_long_endo,
r_long_epi=r_long_epi,
long_axis=long_axis,
)
else:
raise ValueError(
f"Unknown transmural coordinate {transmural!r}, expected 'laplace' or 'analytic'"
)
if outdir is not None:
try:
with dolfinx.io.VTXWriter(mesh.comm, Path(outdir) / f"{transmural}.bp", [t]) as file:
file.write(0.0)
except (RuntimeError, ValueError):
pass
system = compute_system(
t,
function_space=function_space,
r_short_endo=r_short_endo,
r_short_epi=r_short_epi,
r_long_endo=r_long_endo,
r_long_epi=r_long_epi,
alpha_endo=alpha_endo,
alpha_epi=alpha_epi,
long_axis=long_axis,
)
if outdir is not None:
utils.save_microstructure(mesh, system, path=Path(outdir) / "geometry.bp")
return system
