"""
Plot TriSeg.
By Nick van Osta
The function triseg2022(model, ax, idx) plots triseg on ax at time-iteration
idx.
"""
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Wedge
# TODO: Generalize plot for different wall types
# now only model VanOsta2023 is compatible
[docs]
def update_triseg(model, ax, idx):
"""
Updates the already plotted TriSeg.
init_plot has to be called first.
Parameters
----------
model: CircAdapt
CircAdapt model with TriSeg2022 and Wall2022.
ax: ax
Pyplot axis.
idx: int
Time index.
"""
Y = model['TriSeg']['Y'][idx, 0]
V_sv = model['TriSeg']['V'][idx, 0]
V_wall = model['Wall']['V_wall'][2:5]
Am = model['Wall']['Am'][:, 2:5][idx, :]
Cm = model['Wall']['Cm'][:, 2:5][idx, :]
Vm = np.array([
-model['Cavity']['V'][:, 'cLv'][idx] - 0.5*V_wall[0] - 0.5*V_wall[1],
0,
model['Cavity']['V'][:, 'cRv'][idx] + 0.5*V_wall[2] + 0.5*V_wall[1],
]) + V_sv
SignVm = np.sign(Vm)
Vm = np.abs(Vm)
Ym3 = np.array([Y, Y, Y])
V = (3 / np.pi) * Vm
Q = (V + np.sqrt(V**2 + Ym3**6))**(1/3)
Xm = SignVm * (Q - Ym3**2 / Q)
_update_wall(ax.patches[0], Y, Xm[0], V_wall[0], Am[0], Cm[0])
_update_wall(ax.patches[1], Y, Xm[1], V_wall[1], Am[1], Cm[1])
_update_wall(ax.patches[2], Y, Xm[2], V_wall[2], Am[2], Cm[2])
node_width = 0.5 * np.max(V_wall/Am)
_update_nodes(ax.patches[3], Y, node_width)
_update_nodes(ax.patches[4], -Y, node_width)
def _update_nodes(wedge, Y, d):
wedge.set_center((0, Y))
wedge.set_radius(d)
wedge.set_width(d)
def _update_wall(wedge, Y, Xm, V_wall, Am, Cm):
center = _update_wall_center(wedge, Cm, Xm, V_wall, Am)
_update_wall_thickness(wedge, V_wall, Am)
_update_wall_radius(wedge, Cm, V_wall, Am)
_update_wall_theta(wedge, center, Y, Cm, Xm, V_wall, Am)
def _update_wall_theta(wedge, center, Y, Cm, Xm, V_wall, Am):
# calculate and set angle
if np.sign(center) != np.sign(Xm) and Xm < 0:
theta1 = -90 - np.arccos(Y * abs(Cm)) * 180 / np.pi
theta2 = 90 + np.arccos(Y * abs(Cm)) * 180 / np.pi
elif np.sign(center) != np.sign(Xm):
theta1 = np.arcsin(Y * abs(Cm)) * 180 / np.pi
theta2 = -np.arcsin(Y * abs(Cm)) * 180 / np.pi
elif Xm < 0:
theta1 = -90 + np.arccos(Y * abs(Cm)) * 180 / np.pi
theta2 = 90 - np.arccos(Y * abs(Cm)) * 180 / np.pi
else:
theta1 = 180 - np.arcsin(Y * abs(Cm)) * 180 / np.pi
theta2 = 180 + np.arcsin(Y * abs(Cm)) * 180 / np.pi
wedge.set_theta1(theta1)
wedge.set_theta2(theta2)
def _update_wall_center(wedge, Cm, Xm, V_wall, Am):
if Xm < 0:
center = 1/Cm - Xm + 0*V_wall/Am
else:
center = -Xm + 1/Cm - 0.*V_wall/Am
wedge.set_center((center, 0))
return center
def _update_wall_thickness(wedge, V_wall, Am):
wedge.set_width(V_wall/Am)
def _update_wall_radius(wedge, Cm, V_wall, Am):
R = np.abs(1/Cm)
wedge.set_radius(R+0.5*V_wall/Am)
[docs]
def init_plot(model, ax, colors=None):
if colors is None:
colors = [None, None, None]
# Define the initial radius and the number of frames
initial_radius = 1.0
wall_wedges = [
# LV, SV, RV
Wedge((0, 0), initial_radius, 0, 360, alpha=1.0, width=0.,
fc=colors[0],
),
Wedge((1, 0), initial_radius, 0, 360, alpha=1.0, width=0.,
fc=colors[1],
),
Wedge((0, 0), initial_radius, 0, 360, alpha=1.0, width=0.,
fc=colors[2],
),
# nodes
Wedge((0, 0), initial_radius, 0, 360, alpha=1.0, width=0.,
fc=[1, 1, 1],
ec=[0, 0, 0],
lw=2,
),
Wedge((0, 0), initial_radius, 0, 360, alpha=1.0, width=0.,
fc=[1, 1, 1],
ec=[0, 0, 0],
lw=2,
)
]
for w in wall_wedges:
ax.add_patch(w)
Y = model['TriSeg']['Y'][:, 0]
V_sv = model['TriSeg']['V'][:, 0]
V_wall = model['Wall']['V_wall'][2:5]
Am = model['Wall']['Am'][:, 2:5][:, :]
Vm = [
-model['Cavity']['V'][:, 'cLv'] - 0.5*V_wall[0] - 0.5*V_wall[1] + V_sv,
V_sv,
model['Cavity']['V'][:, 'cRv'][:] + 0.5*V_wall[2] + 0.5*V_wall[1] + V_sv,
]
SignVm = np.sign(Vm)
Vm = np.abs(Vm)
Ym3 = np.array([Y, Y, Y])
V = (3/np.pi)*Vm
Q = (V + np.sqrt(V**2 + Ym3**6))**(1/3)
Xm = SignVm * ( Q - Ym3**2 / Q )
lim = np.max([np.max(Y), np.max(np.abs(Xm))]) + np.max(V_wall/Am)
plt.axis('equal')
ax.set_xlim(-lim, lim)
ax.set_ylim(-lim, lim)
[docs]
def triseg2022(model, ax, idx, colors=None):
init_plot(model, ax, colors)
update_triseg(model, ax, idx)
[docs]
def mmode(model, ax):
"""
Plot MMode for TriSeg Module
Parameters
----------
model : CircAdapt model
Model needs to have an TriSeg2022 object.
ax : Matplotlib ax object
Used to plot mmode to
"""
# Calculate
Y = model['TriSeg']['Y'][:, 0]
V_sv = model['TriSeg']['V'][:, 0]
V_wall = model['Wall']['V_wall'][2:5]
Am = model['Wall']['Am'][:, 2:5][:, :]
d = V_wall / Am
Vm = [
-model['Cavity']['V'][:, 'cLv'] - 0.5*V_wall[0] - 0.5*V_wall[1] + V_sv,
V_sv,
model['Cavity']['V'][:, 'cRv'][:] + 0.5*V_wall[2] + 0.5*V_wall[1] + V_sv,
]
SignVm = np.sign(Vm)
Vm = np.abs(Vm)
Ym3 = np.array([Y, Y, Y])
V = (3/np.pi)*Vm
Q = (V + np.sqrt(V**2 + Ym3**6))**(1/3)
Xm = SignVm * ( Q - Ym3**2 / Q ) * 1e3
d *= 1e3
t = (model['Solver']['t']-model['Solver']['t'][0]) * 1e3
def split_blocks(arr):
true_blocks = []
current_block = []
original_indices = []
array_to_look_true = np.array(arr[:])
array_to_look_false = np.array(arr[:])
array_to_look_true[:np.argmax(array_to_look_true)] = True
while np.any(array_to_look_false):
array_to_add = np.zeros(len(arr), dtype=bool)
idx0 = np.argmax(array_to_look_false)
array_to_look_true[:idx0] = True
idx1 = np.argmin(array_to_look_true)
if (idx1 == 0):
idx1 = None
array_to_add[idx0:idx1] = True
true_blocks.append(array_to_add)
# prepare for next
array_to_look_false[:idx1] = False
return true_blocks
# LV phases
t_cycle = model['General']['t_cycle']
n_beats = round((t[-1] - t[0]) / t_cycle)
for i_beat in range(n_beats):
idx_beat = (t > (t[0] + i_beat * t_cycle*1e3)) & (t < (t[0] + (i_beat + 1) * t_cycle*1e3))
max_flow = np.max(model['Valve']['q'][:, ['LvSyArt', 'RvPuArt']], axis=0)
idx_ejection = (model['Valve']['q'][:, ['LvSyArt', 'RvPuArt']] > 0.01*max_flow) & np.array([idx_beat, idx_beat]).T
for i in range(2):
idx_ejection_plot = idx_ejection[:, i]
filled_poly = plt.fill_between(t[idx_ejection_plot],
Xm[i, idx_ejection_plot]+0.5*d[idx_ejection_plot, i],
Xm[i+1, idx_ejection_plot]-0.5*d[idx_ejection_plot, i+1],
color=[0.5,0.5,0.5])
# (x0, y0), (x1, y1) = filled_poly.get_paths()[0].get_extents().get_points()
# plt.text((x0 + x1) / 2, (y0 + y1) / 2, "Ejection", ha='center', va='center', fontsize=16)
# filling
max_flow = np.max(model['Valve']['q'][:, ['LaLv', 'RaRv']], axis=0)
idx_filling = (model['Valve']['q'][:, ['LaLv', 'RaRv']] > 0.01*max_flow) & np.array([idx_beat, idx_beat]).T
for i in range(2):
blocks_filling = split_blocks(idx_filling[:, i])
for block in blocks_filling:
filled_poly = plt.fill_between(t[block],
Xm[i, block]+0.5*d[block, i],
Xm[i+1, block]-0.5*d[block, i+1],
color=[0.9, 0.9, .9])
for i_wall in range(3):
plt.fill_between(t, Xm[i_wall, :]-0.5*d[:, i_wall], Xm[i_wall, :]+0.5*d[:, i_wall])
plt.plot(t, Xm[i_wall, :])