diff --git a/tests/cases/struct/LV_HolzapfelOgden_active/LV_stress.dat b/tests/cases/struct/LV_HolzapfelOgden_active/LV_stress.dat new file mode 100644 index 000000000..2bd2d3038 --- /dev/null +++ b/tests/cases/struct/LV_HolzapfelOgden_active/LV_stress.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:18f506b3429a491b57ee8073631ba698cf4b79b5fc68a19ededefe12b5b52f6f +size 20239 diff --git a/tests/cases/struct/LV_HolzapfelOgden_active/README.md b/tests/cases/struct/LV_HolzapfelOgden_active/README.md new file mode 100644 index 000000000..d6438600b --- /dev/null +++ b/tests/cases/struct/LV_HolzapfelOgden_active/README.md @@ -0,0 +1,21 @@ +This test case simulates an idealized left ventricle with a Holzapfel-Ogden material model +contracting due to time-dependent active stress, and subject to a time-dependent +pressure load on the endocardial surface. The full problem is described in +case 1A of the cardiac elastodynamcis benchmark paper by Aróstica et al. (2025)[1]. A comparison of the displacement of two points throughout the cardiac cycle as computed by multiple solvers including svMultiphysics is shown below: + +![Displacement Benchmark](comparison_plots_p0_p1_step_1_nonblinded.png) + +Aditionally, we present a pressure-volume loop for the idealized left ventricle. Note that the time-dependent pressure load in this problem is not intended to reflect physiological conditions. + +![P-V loop](p-v_loop.png) + +These plots can be generated in figures.py by dowloading the benchmark dataset from https://zenodo.org/records/14260459 and running solver.xml with 1000 time steps to obtain the svMultiphysics results for the entire cycle. + +[1]: Reidmen Aróstica, David Nolte, Aaron Brown, Amadeus Gebauer, Elias Karabelas, Javiera Jilberto, Matteo Salvador, Michele Bucelli, Roberto Piersanti, Kasra Osouli, Christoph Augustin, Henrik Finsberg, Lei Shi, Marc Hirschvogel, Martin Pfaller, Pasquale Claudio Africa, Matthias Gsell, Alison Marsden, David Nordsletten, Francesco Regazzoni, Gernot Plank, Joakim Sundnes, Luca Dede’, Mathias Peirlinck, Vijay Vedula, Wolfgang Wall, Cristóbal Bertoglio, +A software benchmark for cardiac elastodynamics, +Computer Methods in Applied Mechanics and Engineering, +Volume 435, +2025, +117485, +ISSN 0045-7825, +https://doi.org/10.1016/j.cma.2024.117485. \ No newline at end of file diff --git a/tests/cases/struct/LV_HolzapfelOgden_active/comparison_plots_p0_p1_step_1_nonblinded.png b/tests/cases/struct/LV_HolzapfelOgden_active/comparison_plots_p0_p1_step_1_nonblinded.png new file mode 100644 index 000000000..4dd42c905 --- /dev/null +++ b/tests/cases/struct/LV_HolzapfelOgden_active/comparison_plots_p0_p1_step_1_nonblinded.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:87b5271005953d00271a71467688b5bae80ab366ed245df6265098b0b46e34f2 +size 383754 diff --git a/tests/cases/struct/LV_HolzapfelOgden_active/figures.py b/tests/cases/struct/LV_HolzapfelOgden_active/figures.py new file mode 100644 index 000000000..f6a9ad757 --- /dev/null +++ b/tests/cases/struct/LV_HolzapfelOgden_active/figures.py @@ -0,0 +1,502 @@ +""" +This script plots the displacement and pressure-volume loop for the LV_HolzapfelOgden_active case +and compares results to those found in the benchmark paper by Aróstica et al. (2025). + +The resulting plots are included in this directory. If the user wishes to generate the plots, +they can do so by downloading the benchmark dataset from https://zenodo.org/records/14260459 +and running solver.xml with 1000 time steps to obtain the svMultiPhysics results. +""" + +import pickle +from pathlib import Path +import numpy as np +import matplotlib.pyplot as plt +import seaborn as sns +import os +import pyvista as pv +from scipy.interpolate import interp1d + +sns.set_theme() + +DATA_PATH = Path(__file__).parent / "data" + +# if DATA_PATH doesn't exist, print error messsage and exit +if not DATA_PATH.exists(): + print(f"Error: Benchmark path {DATA_PATH} does not exist. Benchmark results can be found at https://zenodo.org/records/14260459") + exit(1) + +# Path to svMultiPhysics results +results_path = 'results' +if not Path(results_path).exists(): + print(f"Error: svMultiPhysics results path {results_path} does not exist. Run solver.xml with 1000 time steps to generate results.") + exit(1) + + +def load_dataset_from_pickle(filename: str | Path) -> dict[str, np.ndarray]: + if not isinstance(filename, Path): + filename = Path(filename) + + if not filename.exists(): + raise FileNotFoundError(f"File {filename} not found") + + with open(filename.as_posix(), "rb") as fl: + return pickle.load(fl) + + +def get_array_at_point(mesh, point, array_name): + """Interpolate array value at an arbitrary point using PyVista's sample method.""" + # Create a PolyData object for the query point + point_poly = pv.PolyData(np.array([point])) + # Interpolate the array at the point + sampled = point_poly.sample(mesh) + if array_name in sampled.point_data: + return sampled.point_data[array_name][0] + else: + raise ValueError(f"Array '{array_name}' not found in mesh") + + +def get_displacements_at_points(start_timestep, end_timestep, step, timestep_size, results_folder, sample_points): + """Get displacements at specified points over time.""" + t = [] + displacements = [] + + for k in range(start_timestep, end_timestep + 1, step): + # Load results VTU mesh + result = pv.read(os.path.join(results_folder, f"result_{k:03d}.vtu")) + + # Get displacement at each sample point + point_displacements = [] + for point in sample_points: + if k == 0: + # At time 0, displacement is zero + disp = np.array([0.0, 0.0, 0.0]) + else: + # Get displacement at the point + disp = get_array_at_point(result, point, 'Displacement') + point_displacements.append(disp) + + t.append(k * timestep_size) + displacements.append(point_displacements) + + return np.array(t), np.array(displacements) + + +def calc_volume_3D(start_timestep, end_timestep, step, timestep_size, results_folder, reference_surface, save_intermediate_data=False, intermediate_output_folder=None): + """Calculate volume over time using 3D mesh.""" + t = [] + vol = [] + + # Load reference surface + ref_lumen = pv.read(reference_surface) + + for k in range(start_timestep, end_timestep + 1, step): + # Load results VTU mesh + result = pv.read(os.path.join(results_folder, f"result_{k:03d}.vtu")) + + # Sample result onto ref_lumen + resampled_lumen = ref_lumen.sample(result) + + # Warp resampled surface by displacement + warped_lumen = resampled_lumen.warp_by_vector('Displacement') + + # Save warped and filled lumen if requested + if save_intermediate_data: + warped_lumen.save(os.path.join(intermediate_output_folder, f'resampled_warped_and_filled_{k:03d}.vtp')) + + # Add time and volume to arrays + t.append(k * timestep_size) + vol.append(warped_lumen.volume) + + print(f"Iteration: {k}, Volume: {warped_lumen.volume}") + + return (t, vol) + + +def get_svmultiphysics_displacement_data(): + """Get svMultiPhysics displacement data in the same format as benchmark datasets.""" + # Set points needed to calculate the displacements + p0 = np.array([0.025, 0.03, 0.0]) + p1 = np.array([0.0, 0.03, 0.0]) + sample_points = [p0, p1] + results_path = 'results' + + # Get the displacements at the sample points + t_displacements, displacements = get_displacements_at_points(10, 1000, 10, 1e-3, results_path, sample_points) + + # Format data to match benchmark dataset structure + svmultiphysics_data = { + "time": t_displacements, + "displacement": { + "p0": { + "ux": displacements[:, 0, 0], + "uy": displacements[:, 0, 1], + "uz": displacements[:, 0, 2] + }, + "p1": { + "ux": displacements[:, 1, 0], + "uy": displacements[:, 1, 1], + "uz": displacements[:, 1, 2] + } + } + } + + return svmultiphysics_data + + +def plot_pv_loop_only(): + """Plot only the pressure-volume loop.""" + # Plot P-V loop + results_path = 'results' + t, vol = calc_volume_3D(10, 1000, 10, 1e-3, results_path, 'mesh/mesh-surfaces/endo.vtp') + vol = np.array(vol) * 1e6 # Convert to mL + + # Load pressure data from the second column of pressure.dat + pressure_data = np.loadtxt('pressure.dat') + pressure_time = pressure_data[1:, 0] # First column: time + pressure_values = pressure_data[1:, 1] / 133.322387415 # Second column: pressure + + # Interpolate pressure to match volume time points + pressure_interp = interp1d(pressure_time, pressure_values, kind='linear', bounds_error=False, fill_value='extrapolate') + pressure_at_volume_times = pressure_interp(t) + + plt.figure(figsize=(10, 8)) + plt.plot(vol, pressure_at_volume_times, 'b-', linewidth=2, label='Pressure-Volume Loop') + plt.xlabel('Volume [mL]', fontsize=12) + plt.ylabel('Pressure [mmHg]', fontsize=12) + plt.title('Pressure-Volume Loop', fontsize=14) + plt.grid(True, alpha=0.3) + + plt.tight_layout() + plt.savefig('p-v_loop.png', dpi=300, bbox_inches='tight') + plt.show() + + # Print some statistics + print(f"Maximum volume: {np.max(vol):.2f} mL") + print(f"Minimum volume: {np.min(vol):.2f} mL") + print(f"Stroke volume: {np.max(vol) - np.min(vol):.2f} mL") + print(f"Maximum pressure: {np.max(pressure_at_volume_times):.2f} mmHg") + print(f"Minimum pressure: {np.min(pressure_at_volume_times):.2f} mmHg") + + +def plot_displacements_and_pv_loop(): + """Plot displacements and pressure-volume loop from process_results.py.""" + # Get svMultiPhysics data + svmultiphysics_data = get_svmultiphysics_displacement_data() + + # Load benchmark data + data = np.load('Step_1_US_P1_h5.npz') + coords = ['x', 'y', 'z'] + + # Plot the displacements for both points and for each of three coordinates on 2 x 3 grid + plt.figure(figsize=(12, 12)) + for i in range(3): + plt.subplot(3, 2, 2*i+1) + plt.plot(svmultiphysics_data['time'], svmultiphysics_data['displacement']['p0'][f'u{coords[i]}'], label='svMultiPhysics') + plt.plot(data['time'], data['u_0'][:, i], label='benchmark') + plt.xlabel('Time [s]') + plt.ylabel(coords[i] + ' Displacement [m]') + plt.legend() + plt.tight_layout() + + plt.subplot(3, 2, 2*i+2) + plt.plot(svmultiphysics_data['time'], svmultiphysics_data['displacement']['p1'][f'u{coords[i]}'], label='svMultiPhysics') + plt.plot(data['time'], data['u_1'][:, i], label='benchmark') + plt.xlabel('Time [s]') + plt.ylabel(coords[i] + ' Displacement [m]') + plt.legend() + plt.tight_layout() + + # Label column 1 'p_0' and column 2 'p_1' + plt.subplot(3, 2, 1) + plt.title(r'$p_0$') + plt.subplot(3, 2, 2) + plt.title(r'$p_1$') + + plt.tight_layout() + plt.savefig('displacements.png', dpi=300, bbox_inches='tight') + plt.show() + + # Plot P-V loop + results_path = 'results' + t, vol = calc_volume_3D(0, 1000, 10, 1e-3, results_path, 'mesh/mesh-surfaces/endo.vtp') + vol = np.array(vol) * 1e6 # Convert to mL + + # Load pressure data from the second column of pressure.dat + pressure_data = np.loadtxt('pressure.dat') + pressure_time = pressure_data[1:, 0] # First column: time + pressure_values = pressure_data[1:, 1] / 133.322387415 # Second column: pressure + + # Interpolate pressure to match volume time points + pressure_interp = interp1d(pressure_time, pressure_values, kind='linear', bounds_error=False, fill_value='extrapolate') + pressure_at_volume_times = pressure_interp(t) + + plt.figure(figsize=(10, 8)) + plt.plot(vol, pressure_at_volume_times, 'b-', linewidth=2, label='Pressure-Volume Loop') + plt.xlabel('Volume [mL]', fontsize=12) + plt.ylabel('Pressure [mmHg]', fontsize=12) + plt.title('Pressure-Volume Loop', fontsize=14) + plt.grid(True, alpha=0.3) + + plt.tight_layout() + plt.savefig('p-v_loop.png', dpi=300, bbox_inches='tight') + plt.show() + + # Print some statistics + print(f"Maximum volume: {np.max(vol):.2f} mL") + print(f"Minimum volume: {np.min(vol):.2f} mL") + print(f"Stroke volume: {np.max(vol) - np.min(vol):.2f} mL") + print(f"Maximum pressure: {np.max(pressure_at_volume_times):.2f} mmHg") + print(f"Minimum pressure: {np.min(pressure_at_volume_times):.2f} mmHg") + + +LABEL_NAMES = [ + "CARPentry", + "Ambit", + "4C", + "Simula", + "CHimeRA", + r"$\mathcal{C}$Heart", + r"life$^{\mathbf{X}}$", + r"SimVascular $\mathbb{P}_1$", + r"SimVascular $\mathbb{P}_2$", + "COMSOL", +] + +LABEL_NAMES_BIV = [ + "CARPentry", + "Ambit", + "4C", + "Simula", + "CHimeRA", + r"$\mathcal{C}$Heart", + r"life$^{\mathbf{X}}$", + r"SimVascular $\mathbb{P}_1$", + "COMSOL", +] + +COLORS = [ + "tab:blue", + "tab:orange", + "tab:green", + "tab:red", + "tab:purple", + "tab:brown", + "tab:pink", + "tab:gray", + "tab:cyan", + "tab:olive", +] + +COLORS_BIV = [ + "tab:blue", + "tab:orange", + "tab:green", + "tab:red", + "tab:purple", + "tab:brown", + "tab:pink", + "tab:gray", + "tab:olive", +] + +LABELS_P1 = [ + r"Simula-$\mathbb{P}_1$", + r"CHimeRA-$\mathbb{P}_1$", + r"$\mathcal{C}$Heart-$\mathbb{P}_1$", + r"life$^{\mathbf{X}}$-$\mathbb{P}_1$", + r"COMSOL-$\mathbb{P}_1$", +] + +LABELS_P2 = [ + r"Simula-$\mathbb{P}_2$", + r"CHimeRA-$\mathbb{P}_2$", + r"$\mathcal{C}$Heart-$\mathbb{P}_2$", + r"life$^{\mathbf{X}}$-$\mathbb{P}_2$", + r"COMSOL-$\mathbb{P}_2$", +] + +COLORS_P1_P2 = ["tab:red", "tab:purple", "tab:brown", "tab:pink", "tab:olive"] + + +TEAMS_DATASETS_1 = [ + load_dataset_from_pickle( + DATA_PATH / "monoventricular_nonblinded_step_1_group_carpentry.pickle" + ), + load_dataset_from_pickle( + DATA_PATH / "monoventricular_nonblinded_step_1_group_ambit.pickle" + ), + load_dataset_from_pickle( + DATA_PATH / "monoventricular_nonblinded_step_1_group_4c.pickle" + ), + load_dataset_from_pickle( + DATA_PATH / "monoventricular_nonblinded_step_1_group_simula.pickle" + ), + load_dataset_from_pickle( + DATA_PATH / "monoventricular_nonblinded_step_1_group_chimera.pickle" + ), + load_dataset_from_pickle( + DATA_PATH / "monoventricular_nonblinded_step_1_group_cheart.pickle" + ), + load_dataset_from_pickle( + DATA_PATH / "monoventricular_nonblinded_step_1_group_lifex.pickle" + ), + load_dataset_from_pickle( + DATA_PATH / "monoventricular_nonblinded_step_1_group_simvascular_p1p1.pickle" + ), + load_dataset_from_pickle( + DATA_PATH / "monoventricular_nonblinded_step_1_group_simvascular_p2.pickle" + ), + load_dataset_from_pickle( + DATA_PATH / "monoventricular_nonblinded_step_1_group_comsol.pickle" + ), +] + + +def compute_plot_displacement_monoventricular( + teams_datasets: list[dict[str, np.ndarray]], + filename: str | Path, + labels_names: list[str], + colors: list[str], +) -> None: + """Computes the displacement figures based on the groups datasets.""" + if isinstance(filename, str): + filename = Path(filename) + + if not isinstance(teams_datasets, list): + raise Exception("teams_datasets must be a list") + + if not len(teams_datasets) == len(labels_names): + raise Exception("teams_datasets and labels_names must have the same length") + + if not len(teams_datasets) == len(colors): + raise Exception("teams_datasets and colors must have the same length") + + filename.parent.mkdir(parents=True, exist_ok=True) + + fig, axs = plt.subplots(3, 2, sharex=True, figsize=(12, 12)) + + axs[0, 0].set_title(r"Particle $\mathbf{p}_0$") + axs[0, 0].set_ylabel("Displacement x-component [m]") + axs[1, 0].set_ylabel("Displacement y-component [m]") + axs[2, 0].set_ylabel("Displacement z-component [m]") + axs[2, 0].set_xlabel("Time [s]") + + axs[0, 1].set_title(r"Particle $\mathbf{p}_1$") + axs[2, 1].set_xlabel("Time [s]") + + for data, lbl, color in zip(teams_datasets, labels_names, colors): + for i, u_type in zip([0, 1, 2], ["ux", "uy", "uz"]): + if lbl == "svMultiPhysics": + axs[i, 0].plot( + data["time"], data["displacement"]["p0"][u_type], label=lbl, color="k", linestyle="--" + ) + axs[i, 1].plot( + data["time"], data["displacement"]["p1"][u_type], label=lbl, color="k", linestyle="--" + ) + else: + axs[i, 0].plot( + data["time"], data["displacement"]["p0"][u_type], label=lbl, color=color + ) + axs[i, 1].plot( + data["time"], data["displacement"]["p1"][u_type], label=lbl, color=color + ) + + plt.legend(loc="best", fancybox=True, shadow=True) + fig.tight_layout() + fig.savefig(filename.as_posix(), bbox_inches="tight", dpi=120) + # plt.show() + + +def compute_statistics( + datasets: list[dict[str, np.ndarray]], point: str = "p0" +) -> dict[str, np.ndarray]: + time = datasets[0]["time"] + number_of_element_per_dataset = datasets[0]["time"].shape[0] + number_of_datasets = len(datasets) + condensated_dataset_ux = np.zeros( + (number_of_element_per_dataset, number_of_datasets) + ) + condensated_dataset_uy = np.zeros_like(condensated_dataset_ux) + condensated_dataset_uz = np.zeros_like(condensated_dataset_ux) + red_dataset_u = np.zeros((number_of_datasets,)) + + for i, dataset in enumerate(datasets): + condensated_dataset_ux[:, i] = dataset["displacement"][point]["ux"] + condensated_dataset_uy[:, i] = dataset["displacement"][point]["uy"] + condensated_dataset_uz[:, i] = dataset["displacement"][point]["uz"] + + mean_ux = np.mean(condensated_dataset_ux, axis=1) + mean_uy = np.mean(condensated_dataset_uy, axis=1) + mean_uz = np.mean(condensated_dataset_uz, axis=1) + + for i in range(number_of_datasets): + diff_to_mean_norm = np.sqrt( + np.abs(condensated_dataset_ux[:, i] - mean_ux) ** 2 + + np.abs(condensated_dataset_uy[:, i] - mean_uy) ** 2 + + np.abs(condensated_dataset_uz[:, i] - mean_uz) ** 2 + ) + + mean_norm = np.sqrt( + np.abs(mean_ux) ** 2 + np.abs(mean_uy) ** 2 + np.abs(mean_uz) ** 2 + ) + + red_dataset_u[i] = np.mean(diff_to_mean_norm / mean_norm, axis=0) + + statistics_dataset = { + "time": time, + "mean_ux": mean_ux, + "mean_uy": mean_uy, + "mean_uz": mean_uz, + "std_ux": np.std(condensated_dataset_ux, axis=1), + "std_uy": np.std(condensated_dataset_uy, axis=1), + "std_uz": np.std(condensated_dataset_uz, axis=1), + "red_u": red_dataset_u, + } + + return statistics_dataset + + +def compute_plots_monoventricular_nonblinded_step_1() -> None: + """Computes displacement curves for the monoventricular non-blinded step 1""" + # Get svMultiPhysics data + svmultiphysics_data = get_svmultiphysics_displacement_data() + + # Add svMultiPhysics to the datasets, labels, and colors + all_datasets = TEAMS_DATASETS_1 + [svmultiphysics_data] + all_labels = LABEL_NAMES + ["svMultiPhysics"] + all_colors = COLORS + ["tab:red"] # Add red color for svMultiPhysics + + compute_plot_displacement_monoventricular( + all_datasets, + "./comparison_plots_p0_p1_step_1_nonblinded.png", + all_labels, + all_colors, + ) + + +def compute_statistics_per_dataset( + labels: list[str], datasets: list[dict[str, np.ndarray]], header: str = "" +) -> None: + if len(labels) != len(datasets): + raise ValueError("labels and datasets do not match in length") + + print(header) + statistics_p0 = compute_statistics(datasets, point="p0") + statistics_p1 = compute_statistics(datasets, point="p1") + + red_p0 = np.round(statistics_p0["red_u"], decimals=3) + red_p1 = np.round(statistics_p1["red_u"], decimals=3) + + print(f"Team {'TEAM NAME':25s} - {'p0':4s}, {'p1':4s}") + for label, i in zip(labels, range(len(red_p0))): + print(f"Team {label:25s} - {red_p0[i]:4.3f}, {red_p1[i]:4.3f}") + + print("-----------------------------------\n") + + +if __name__ == "__main__": + compute_plots_monoventricular_nonblinded_step_1() + + # Also plot PV loop separately + print("Plotting pressure-volume loop...") + plot_pv_loop_only() diff --git a/tests/cases/struct/LV_HolzapfelOgden_active/generate_pressure.py b/tests/cases/struct/LV_HolzapfelOgden_active/generate_pressure.py new file mode 100644 index 000000000..7d3c3a984 --- /dev/null +++ b/tests/cases/struct/LV_HolzapfelOgden_active/generate_pressure.py @@ -0,0 +1,131 @@ +import numpy as np +import matplotlib.pyplot as plt +from scipy.integrate import solve_ivp + +params = {} + +params['alpha_min'] = -30 +params['alpha_max'] = 5 +params['alpha_pre'] = 5 +params['alpha_mid'] = 1 +params['sigma_pre'] = 7000 +params['sigma_mid'] = 16000 +params['sigma_0'] = 1.5e5 +params['t_sys_pre'] = 0.17 +params['t_dias_pre'] = 0.484 +params['t_sys'] = 0.16 +params['t_dias'] = 0.484 +params['gamma'] = 0.005 + +def S_plus(delta_t, params): + return 0.5* (1 + np.tanh(delta_t/params['gamma'])) + +def S_minus(delta_t, params): + return 0.5* (1 - np.tanh(delta_t/params['gamma'])) + +def g_pre(t, params): + return S_minus(t - params['t_dias_pre'], params) + +def f(t, params): + return S_plus(t - params['t_sys'], params) * S_minus(t - params['t_dias'], params) + +def f_pre(t, params): + return S_plus(t - params['t_sys_pre'], params) * S_minus(t - params['t_dias_pre'], params) + +def a(t, params): + return params['alpha_max'] * f(t, params) + params['alpha_min'] * (1-f(t, params)) + +def a_pre(t, params): + return params['alpha_max'] * f_pre(t, params) + params['alpha_min'] * (1-f_pre(t, params)) + +def b(t, params): + return a_pre(t, params) + params['alpha_pre'] * g_pre(t, params) + params['alpha_mid'] + +# Differential equation for pressure +def dpdt(t, p): + b_abs = abs(b(t, params)) + b_max = max(b(t, params), 0) + g_pre_max = max(g_pre(t, params), 0) + return -b_abs*p + params['sigma_mid']*b_max + params['sigma_pre']*g_pre_max + +# Differential equation for stress +def dtaudt(t, tau): + a_abs = abs(a(t, params)) + a_max = max(a(t, params), 0) + return -a_abs*tau + params['sigma_0']*a_max + + +def save_sol(sol, filename, f_modes = 512): + # save data to .dat file + # first row is [len(sol.y[0]) f_modes] + # first column is sol.t rounded to six decimal places + # second column is sol.y[0] rounded to six decimal places + + data = np.column_stack((np.round(sol.t, 6), np.round(sol.y[0], 6))) + + data = np.vstack((np.array([len(sol.y[0]), f_modes]), data)) + + # Save to file with the first row as integers and the rest as floats + with open(filename, 'w') as f: + # Write the first row as integers + np.savetxt(f, [data[0]], fmt='%d', delimiter=' ') + # Write the remaining rows as floats with 6 decimal places + np.savetxt(f, data[1:], fmt='%1.6f', delimiter=' ') + +# Solve the ODE for pressure + +p0 = 0 +t_eval = np.linspace(0, 1, 1001) +p_sol = solve_ivp(dpdt, [0, 1], [p0], t_eval=t_eval, method='DOP853') +print('Max pressure:', max(p_sol.y[0])) + +plt.plot(p_sol.t, p_sol.y[0]) +plt.xlabel('Time') +plt.ylabel('Pressure') +''' +# load pressure_original.dat and convert string to float +p_original = np.loadtxt('pressure_original.dat', delimiter=',') +p_original[1:,1] = 0.1*p_original[1:,1] +plt.plot(p_original[1:,0], p_original[1:,1]) + +with open('pressure_original_Pa.dat', 'w') as f: + # Write the first row as integers + np.savetxt(f, [p_original[0]], fmt='%d', delimiter=' ') + # Write the remaining rows as floats with 6 decimal places + np.savetxt(f, p_original[1:], fmt='%1.6f', delimiter=' ') +''' +plt.show() + + +# Save the solution to a file +save_sol(p_sol, 'pressure.dat', 512) + +# Solve the ODE for stress + +tau0 = 0 +t_eval = np.linspace(0, 1, 1001) +tau_sol = solve_ivp(dtaudt, [0, 1], [tau0], t_eval=t_eval) + +print('Max stress:', max(tau_sol.y[0])) + +plt.plot(tau_sol.t, tau_sol.y[0]) +plt.xlabel('Time') +plt.ylabel('Stress') +''' +# load stress_original.dat and convert string to float +tau_original = np.loadtxt('stress_original.dat', delimiter=' ') +tau_original[1:,1] = 0.1*tau_original[1:,1] +plt.plot(tau_original[1:,0], tau_original[1:,1]) + +with open('stress_original_Pa.dat', 'w') as f: + # Write the first row as integers + np.savetxt(f, [tau_original[0]], fmt='%d', delimiter=' ') + # Write the remaining rows as floats with 6 decimal places + np.savetxt(f, tau_original[1:], fmt='%1.6f', delimiter=' ') +''' +plt.show() + +# Save the solution to a file +save_sol(tau_sol, 'LV_stress.dat', 512) + + diff --git a/tests/cases/struct/LV_HolzapfelOgden_active/mesh/fibersLong.vtu b/tests/cases/struct/LV_HolzapfelOgden_active/mesh/fibersLong.vtu new file mode 100644 index 000000000..28e2d5cf4 --- /dev/null +++ b/tests/cases/struct/LV_HolzapfelOgden_active/mesh/fibersLong.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:9c6570939679034cdbba8a2577d5c84baa68339aa8b38da3e2960356bbfcc5d9 +size 868063 diff --git a/tests/cases/struct/LV_HolzapfelOgden_active/mesh/fibersSheet.vtu b/tests/cases/struct/LV_HolzapfelOgden_active/mesh/fibersSheet.vtu new file mode 100644 index 000000000..425397698 --- /dev/null +++ b/tests/cases/struct/LV_HolzapfelOgden_active/mesh/fibersSheet.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1939854b0527c6f78353b12b25aee37c0b82a1202c40e10f4f48d0d2c46ca7a2 +size 868763 diff --git a/tests/cases/struct/LV_HolzapfelOgden_active/mesh/mesh-complete.mesh.vtu b/tests/cases/struct/LV_HolzapfelOgden_active/mesh/mesh-complete.mesh.vtu new file mode 100644 index 000000000..451bf197a --- /dev/null +++ b/tests/cases/struct/LV_HolzapfelOgden_active/mesh/mesh-complete.mesh.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:17f70504f6760f7e836373181fc3487ecc35312bc286fe9b2962e877fe55e517 +size 834394 diff --git a/tests/cases/struct/LV_HolzapfelOgden_active/mesh/mesh-surfaces/endo.vtp b/tests/cases/struct/LV_HolzapfelOgden_active/mesh/mesh-surfaces/endo.vtp new file mode 100644 index 000000000..8bf3ef696 --- /dev/null +++ b/tests/cases/struct/LV_HolzapfelOgden_active/mesh/mesh-surfaces/endo.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:eefb1443d58a90414d3c94b0e6128fc8b3698dc8c25ea33ffeeb573b263de023 +size 68199 diff --git a/tests/cases/struct/LV_HolzapfelOgden_active/mesh/mesh-surfaces/epi.vtp b/tests/cases/struct/LV_HolzapfelOgden_active/mesh/mesh-surfaces/epi.vtp new file mode 100644 index 000000000..ae626668d --- /dev/null +++ b/tests/cases/struct/LV_HolzapfelOgden_active/mesh/mesh-surfaces/epi.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:f08544d5b0cf276d05e999e6116986bea26442506ae51904ce386530a38cf5d1 +size 91341 diff --git a/tests/cases/struct/LV_HolzapfelOgden_active/mesh/mesh-surfaces/top.vtp b/tests/cases/struct/LV_HolzapfelOgden_active/mesh/mesh-surfaces/top.vtp new file mode 100644 index 000000000..40eb6411c --- /dev/null +++ b/tests/cases/struct/LV_HolzapfelOgden_active/mesh/mesh-surfaces/top.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:685fe50130b4f5b2a8a3f4b66cbd381a1741b4bc0f3b5e46168a21553b32215f +size 8884 diff --git a/tests/cases/struct/LV_HolzapfelOgden_active/p-v_loop.png b/tests/cases/struct/LV_HolzapfelOgden_active/p-v_loop.png new file mode 100644 index 000000000..f6d6a597b --- /dev/null +++ b/tests/cases/struct/LV_HolzapfelOgden_active/p-v_loop.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:2b5b11136a27fbc74cdf65b4ffeac17ecc2c491d6480e37683e1c295419df262 +size 190281 diff --git a/tests/cases/struct/LV_HolzapfelOgden_active/pressure.dat b/tests/cases/struct/LV_HolzapfelOgden_active/pressure.dat new file mode 100644 index 000000000..680a13c97 --- /dev/null +++ b/tests/cases/struct/LV_HolzapfelOgden_active/pressure.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:564f6e7c91faa9a4f9234f1d7f4d0ae96997c5b36a072879b6523668dcfd555d +size 20059 diff --git a/tests/cases/struct/LV_HolzapfelOgden_active/result_001.vtu b/tests/cases/struct/LV_HolzapfelOgden_active/result_001.vtu new file mode 100644 index 000000000..8c9c466d6 --- /dev/null +++ b/tests/cases/struct/LV_HolzapfelOgden_active/result_001.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:28755d0b759cc113191b1f54abee26fa67e98e15d1456ee55f79f34e4833b033 +size 1869197 diff --git a/tests/cases/struct/LV_HolzapfelOgden_active/solver.xml b/tests/cases/struct/LV_HolzapfelOgden_active/solver.xml new file mode 100644 index 000000000..fd9632f33 --- /dev/null +++ b/tests/cases/struct/LV_HolzapfelOgden_active/solver.xml @@ -0,0 +1,128 @@ + + + + + 0 + 3 + 1 + 1e-3 + 0.50 + STOP_SIM + + 1 + result + 1 + 1 + + 50 + 0 + + 1 + 1 + 1 + + + + + + + mesh/mesh-complete.mesh.vtu + + + mesh/mesh-surfaces/endo.vtp + + + + mesh/mesh-surfaces/epi.vtp + + + + mesh/mesh-surfaces/top.vtp + + + mesh/fibersLong.vtu + mesh/fibersSheet.vtu + + + + + + + true + 3 + 10 + 1e-9 + + 1000.0 + 1.0e6 + 0.483333 + + + 100.0 + + + + 59.0 + 8.023 + 18472.0 + 16.026 + 2481.0 + 11.12 + 216.0 + 11.436 + 100.0 + + + ST91 + 1e6 + + + LV_stress.dat + false + + + + true + true + true + true + true + true + true + true + + + + + fsils + + 1e-12 + 1e-12 + 1000 + 100 + + + + Robin + 1.0e8 + 5.0e3 + 1 + + + + Robin + 1.0e5 + 5.0e3 + + + + Neu + Unsteady + pressure.dat + false + true + + + + + diff --git a/tests/cases/ustruct/LV_HolzapfelOgden_active/LV_stress.dat b/tests/cases/ustruct/LV_HolzapfelOgden_active/LV_stress.dat new file mode 100644 index 000000000..2bd2d3038 --- /dev/null +++ b/tests/cases/ustruct/LV_HolzapfelOgden_active/LV_stress.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:18f506b3429a491b57ee8073631ba698cf4b79b5fc68a19ededefe12b5b52f6f +size 20239 diff --git a/tests/cases/ustruct/LV_HolzapfelOgden_active/README.md b/tests/cases/ustruct/LV_HolzapfelOgden_active/README.md new file mode 100644 index 000000000..d6438600b --- /dev/null +++ b/tests/cases/ustruct/LV_HolzapfelOgden_active/README.md @@ -0,0 +1,21 @@ +This test case simulates an idealized left ventricle with a Holzapfel-Ogden material model +contracting due to time-dependent active stress, and subject to a time-dependent +pressure load on the endocardial surface. The full problem is described in +case 1A of the cardiac elastodynamcis benchmark paper by Aróstica et al. (2025)[1]. A comparison of the displacement of two points throughout the cardiac cycle as computed by multiple solvers including svMultiphysics is shown below: + +![Displacement Benchmark](comparison_plots_p0_p1_step_1_nonblinded.png) + +Aditionally, we present a pressure-volume loop for the idealized left ventricle. Note that the time-dependent pressure load in this problem is not intended to reflect physiological conditions. + +![P-V loop](p-v_loop.png) + +These plots can be generated in figures.py by dowloading the benchmark dataset from https://zenodo.org/records/14260459 and running solver.xml with 1000 time steps to obtain the svMultiphysics results for the entire cycle. + +[1]: Reidmen Aróstica, David Nolte, Aaron Brown, Amadeus Gebauer, Elias Karabelas, Javiera Jilberto, Matteo Salvador, Michele Bucelli, Roberto Piersanti, Kasra Osouli, Christoph Augustin, Henrik Finsberg, Lei Shi, Marc Hirschvogel, Martin Pfaller, Pasquale Claudio Africa, Matthias Gsell, Alison Marsden, David Nordsletten, Francesco Regazzoni, Gernot Plank, Joakim Sundnes, Luca Dede’, Mathias Peirlinck, Vijay Vedula, Wolfgang Wall, Cristóbal Bertoglio, +A software benchmark for cardiac elastodynamics, +Computer Methods in Applied Mechanics and Engineering, +Volume 435, +2025, +117485, +ISSN 0045-7825, +https://doi.org/10.1016/j.cma.2024.117485. \ No newline at end of file diff --git a/tests/cases/ustruct/LV_HolzapfelOgden_active/comparison_plots_p0_p1_step_1_nonblinded.png b/tests/cases/ustruct/LV_HolzapfelOgden_active/comparison_plots_p0_p1_step_1_nonblinded.png new file mode 100644 index 000000000..f4afc4259 --- /dev/null +++ b/tests/cases/ustruct/LV_HolzapfelOgden_active/comparison_plots_p0_p1_step_1_nonblinded.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1cd01088f5a1eb1771217c8703b7feda9fefefcfa1c0be4243f0f247abd0fa3d +size 365690 diff --git a/tests/cases/ustruct/LV_HolzapfelOgden_active/figures.py b/tests/cases/ustruct/LV_HolzapfelOgden_active/figures.py new file mode 100644 index 000000000..f6a9ad757 --- /dev/null +++ b/tests/cases/ustruct/LV_HolzapfelOgden_active/figures.py @@ -0,0 +1,502 @@ +""" +This script plots the displacement and pressure-volume loop for the LV_HolzapfelOgden_active case +and compares results to those found in the benchmark paper by Aróstica et al. (2025). + +The resulting plots are included in this directory. If the user wishes to generate the plots, +they can do so by downloading the benchmark dataset from https://zenodo.org/records/14260459 +and running solver.xml with 1000 time steps to obtain the svMultiPhysics results. +""" + +import pickle +from pathlib import Path +import numpy as np +import matplotlib.pyplot as plt +import seaborn as sns +import os +import pyvista as pv +from scipy.interpolate import interp1d + +sns.set_theme() + +DATA_PATH = Path(__file__).parent / "data" + +# if DATA_PATH doesn't exist, print error messsage and exit +if not DATA_PATH.exists(): + print(f"Error: Benchmark path {DATA_PATH} does not exist. Benchmark results can be found at https://zenodo.org/records/14260459") + exit(1) + +# Path to svMultiPhysics results +results_path = 'results' +if not Path(results_path).exists(): + print(f"Error: svMultiPhysics results path {results_path} does not exist. Run solver.xml with 1000 time steps to generate results.") + exit(1) + + +def load_dataset_from_pickle(filename: str | Path) -> dict[str, np.ndarray]: + if not isinstance(filename, Path): + filename = Path(filename) + + if not filename.exists(): + raise FileNotFoundError(f"File {filename} not found") + + with open(filename.as_posix(), "rb") as fl: + return pickle.load(fl) + + +def get_array_at_point(mesh, point, array_name): + """Interpolate array value at an arbitrary point using PyVista's sample method.""" + # Create a PolyData object for the query point + point_poly = pv.PolyData(np.array([point])) + # Interpolate the array at the point + sampled = point_poly.sample(mesh) + if array_name in sampled.point_data: + return sampled.point_data[array_name][0] + else: + raise ValueError(f"Array '{array_name}' not found in mesh") + + +def get_displacements_at_points(start_timestep, end_timestep, step, timestep_size, results_folder, sample_points): + """Get displacements at specified points over time.""" + t = [] + displacements = [] + + for k in range(start_timestep, end_timestep + 1, step): + # Load results VTU mesh + result = pv.read(os.path.join(results_folder, f"result_{k:03d}.vtu")) + + # Get displacement at each sample point + point_displacements = [] + for point in sample_points: + if k == 0: + # At time 0, displacement is zero + disp = np.array([0.0, 0.0, 0.0]) + else: + # Get displacement at the point + disp = get_array_at_point(result, point, 'Displacement') + point_displacements.append(disp) + + t.append(k * timestep_size) + displacements.append(point_displacements) + + return np.array(t), np.array(displacements) + + +def calc_volume_3D(start_timestep, end_timestep, step, timestep_size, results_folder, reference_surface, save_intermediate_data=False, intermediate_output_folder=None): + """Calculate volume over time using 3D mesh.""" + t = [] + vol = [] + + # Load reference surface + ref_lumen = pv.read(reference_surface) + + for k in range(start_timestep, end_timestep + 1, step): + # Load results VTU mesh + result = pv.read(os.path.join(results_folder, f"result_{k:03d}.vtu")) + + # Sample result onto ref_lumen + resampled_lumen = ref_lumen.sample(result) + + # Warp resampled surface by displacement + warped_lumen = resampled_lumen.warp_by_vector('Displacement') + + # Save warped and filled lumen if requested + if save_intermediate_data: + warped_lumen.save(os.path.join(intermediate_output_folder, f'resampled_warped_and_filled_{k:03d}.vtp')) + + # Add time and volume to arrays + t.append(k * timestep_size) + vol.append(warped_lumen.volume) + + print(f"Iteration: {k}, Volume: {warped_lumen.volume}") + + return (t, vol) + + +def get_svmultiphysics_displacement_data(): + """Get svMultiPhysics displacement data in the same format as benchmark datasets.""" + # Set points needed to calculate the displacements + p0 = np.array([0.025, 0.03, 0.0]) + p1 = np.array([0.0, 0.03, 0.0]) + sample_points = [p0, p1] + results_path = 'results' + + # Get the displacements at the sample points + t_displacements, displacements = get_displacements_at_points(10, 1000, 10, 1e-3, results_path, sample_points) + + # Format data to match benchmark dataset structure + svmultiphysics_data = { + "time": t_displacements, + "displacement": { + "p0": { + "ux": displacements[:, 0, 0], + "uy": displacements[:, 0, 1], + "uz": displacements[:, 0, 2] + }, + "p1": { + "ux": displacements[:, 1, 0], + "uy": displacements[:, 1, 1], + "uz": displacements[:, 1, 2] + } + } + } + + return svmultiphysics_data + + +def plot_pv_loop_only(): + """Plot only the pressure-volume loop.""" + # Plot P-V loop + results_path = 'results' + t, vol = calc_volume_3D(10, 1000, 10, 1e-3, results_path, 'mesh/mesh-surfaces/endo.vtp') + vol = np.array(vol) * 1e6 # Convert to mL + + # Load pressure data from the second column of pressure.dat + pressure_data = np.loadtxt('pressure.dat') + pressure_time = pressure_data[1:, 0] # First column: time + pressure_values = pressure_data[1:, 1] / 133.322387415 # Second column: pressure + + # Interpolate pressure to match volume time points + pressure_interp = interp1d(pressure_time, pressure_values, kind='linear', bounds_error=False, fill_value='extrapolate') + pressure_at_volume_times = pressure_interp(t) + + plt.figure(figsize=(10, 8)) + plt.plot(vol, pressure_at_volume_times, 'b-', linewidth=2, label='Pressure-Volume Loop') + plt.xlabel('Volume [mL]', fontsize=12) + plt.ylabel('Pressure [mmHg]', fontsize=12) + plt.title('Pressure-Volume Loop', fontsize=14) + plt.grid(True, alpha=0.3) + + plt.tight_layout() + plt.savefig('p-v_loop.png', dpi=300, bbox_inches='tight') + plt.show() + + # Print some statistics + print(f"Maximum volume: {np.max(vol):.2f} mL") + print(f"Minimum volume: {np.min(vol):.2f} mL") + print(f"Stroke volume: {np.max(vol) - np.min(vol):.2f} mL") + print(f"Maximum pressure: {np.max(pressure_at_volume_times):.2f} mmHg") + print(f"Minimum pressure: {np.min(pressure_at_volume_times):.2f} mmHg") + + +def plot_displacements_and_pv_loop(): + """Plot displacements and pressure-volume loop from process_results.py.""" + # Get svMultiPhysics data + svmultiphysics_data = get_svmultiphysics_displacement_data() + + # Load benchmark data + data = np.load('Step_1_US_P1_h5.npz') + coords = ['x', 'y', 'z'] + + # Plot the displacements for both points and for each of three coordinates on 2 x 3 grid + plt.figure(figsize=(12, 12)) + for i in range(3): + plt.subplot(3, 2, 2*i+1) + plt.plot(svmultiphysics_data['time'], svmultiphysics_data['displacement']['p0'][f'u{coords[i]}'], label='svMultiPhysics') + plt.plot(data['time'], data['u_0'][:, i], label='benchmark') + plt.xlabel('Time [s]') + plt.ylabel(coords[i] + ' Displacement [m]') + plt.legend() + plt.tight_layout() + + plt.subplot(3, 2, 2*i+2) + plt.plot(svmultiphysics_data['time'], svmultiphysics_data['displacement']['p1'][f'u{coords[i]}'], label='svMultiPhysics') + plt.plot(data['time'], data['u_1'][:, i], label='benchmark') + plt.xlabel('Time [s]') + plt.ylabel(coords[i] + ' Displacement [m]') + plt.legend() + plt.tight_layout() + + # Label column 1 'p_0' and column 2 'p_1' + plt.subplot(3, 2, 1) + plt.title(r'$p_0$') + plt.subplot(3, 2, 2) + plt.title(r'$p_1$') + + plt.tight_layout() + plt.savefig('displacements.png', dpi=300, bbox_inches='tight') + plt.show() + + # Plot P-V loop + results_path = 'results' + t, vol = calc_volume_3D(0, 1000, 10, 1e-3, results_path, 'mesh/mesh-surfaces/endo.vtp') + vol = np.array(vol) * 1e6 # Convert to mL + + # Load pressure data from the second column of pressure.dat + pressure_data = np.loadtxt('pressure.dat') + pressure_time = pressure_data[1:, 0] # First column: time + pressure_values = pressure_data[1:, 1] / 133.322387415 # Second column: pressure + + # Interpolate pressure to match volume time points + pressure_interp = interp1d(pressure_time, pressure_values, kind='linear', bounds_error=False, fill_value='extrapolate') + pressure_at_volume_times = pressure_interp(t) + + plt.figure(figsize=(10, 8)) + plt.plot(vol, pressure_at_volume_times, 'b-', linewidth=2, label='Pressure-Volume Loop') + plt.xlabel('Volume [mL]', fontsize=12) + plt.ylabel('Pressure [mmHg]', fontsize=12) + plt.title('Pressure-Volume Loop', fontsize=14) + plt.grid(True, alpha=0.3) + + plt.tight_layout() + plt.savefig('p-v_loop.png', dpi=300, bbox_inches='tight') + plt.show() + + # Print some statistics + print(f"Maximum volume: {np.max(vol):.2f} mL") + print(f"Minimum volume: {np.min(vol):.2f} mL") + print(f"Stroke volume: {np.max(vol) - np.min(vol):.2f} mL") + print(f"Maximum pressure: {np.max(pressure_at_volume_times):.2f} mmHg") + print(f"Minimum pressure: {np.min(pressure_at_volume_times):.2f} mmHg") + + +LABEL_NAMES = [ + "CARPentry", + "Ambit", + "4C", + "Simula", + "CHimeRA", + r"$\mathcal{C}$Heart", + r"life$^{\mathbf{X}}$", + r"SimVascular $\mathbb{P}_1$", + r"SimVascular $\mathbb{P}_2$", + "COMSOL", +] + +LABEL_NAMES_BIV = [ + "CARPentry", + "Ambit", + "4C", + "Simula", + "CHimeRA", + r"$\mathcal{C}$Heart", + r"life$^{\mathbf{X}}$", + r"SimVascular $\mathbb{P}_1$", + "COMSOL", +] + +COLORS = [ + "tab:blue", + "tab:orange", + "tab:green", + "tab:red", + "tab:purple", + "tab:brown", + "tab:pink", + "tab:gray", + "tab:cyan", + "tab:olive", +] + +COLORS_BIV = [ + "tab:blue", + "tab:orange", + "tab:green", + "tab:red", + "tab:purple", + "tab:brown", + "tab:pink", + "tab:gray", + "tab:olive", +] + +LABELS_P1 = [ + r"Simula-$\mathbb{P}_1$", + r"CHimeRA-$\mathbb{P}_1$", + r"$\mathcal{C}$Heart-$\mathbb{P}_1$", + r"life$^{\mathbf{X}}$-$\mathbb{P}_1$", + r"COMSOL-$\mathbb{P}_1$", +] + +LABELS_P2 = [ + r"Simula-$\mathbb{P}_2$", + r"CHimeRA-$\mathbb{P}_2$", + r"$\mathcal{C}$Heart-$\mathbb{P}_2$", + r"life$^{\mathbf{X}}$-$\mathbb{P}_2$", + r"COMSOL-$\mathbb{P}_2$", +] + +COLORS_P1_P2 = ["tab:red", "tab:purple", "tab:brown", "tab:pink", "tab:olive"] + + +TEAMS_DATASETS_1 = [ + load_dataset_from_pickle( + DATA_PATH / "monoventricular_nonblinded_step_1_group_carpentry.pickle" + ), + load_dataset_from_pickle( + DATA_PATH / "monoventricular_nonblinded_step_1_group_ambit.pickle" + ), + load_dataset_from_pickle( + DATA_PATH / "monoventricular_nonblinded_step_1_group_4c.pickle" + ), + load_dataset_from_pickle( + DATA_PATH / "monoventricular_nonblinded_step_1_group_simula.pickle" + ), + load_dataset_from_pickle( + DATA_PATH / "monoventricular_nonblinded_step_1_group_chimera.pickle" + ), + load_dataset_from_pickle( + DATA_PATH / "monoventricular_nonblinded_step_1_group_cheart.pickle" + ), + load_dataset_from_pickle( + DATA_PATH / "monoventricular_nonblinded_step_1_group_lifex.pickle" + ), + load_dataset_from_pickle( + DATA_PATH / "monoventricular_nonblinded_step_1_group_simvascular_p1p1.pickle" + ), + load_dataset_from_pickle( + DATA_PATH / "monoventricular_nonblinded_step_1_group_simvascular_p2.pickle" + ), + load_dataset_from_pickle( + DATA_PATH / "monoventricular_nonblinded_step_1_group_comsol.pickle" + ), +] + + +def compute_plot_displacement_monoventricular( + teams_datasets: list[dict[str, np.ndarray]], + filename: str | Path, + labels_names: list[str], + colors: list[str], +) -> None: + """Computes the displacement figures based on the groups datasets.""" + if isinstance(filename, str): + filename = Path(filename) + + if not isinstance(teams_datasets, list): + raise Exception("teams_datasets must be a list") + + if not len(teams_datasets) == len(labels_names): + raise Exception("teams_datasets and labels_names must have the same length") + + if not len(teams_datasets) == len(colors): + raise Exception("teams_datasets and colors must have the same length") + + filename.parent.mkdir(parents=True, exist_ok=True) + + fig, axs = plt.subplots(3, 2, sharex=True, figsize=(12, 12)) + + axs[0, 0].set_title(r"Particle $\mathbf{p}_0$") + axs[0, 0].set_ylabel("Displacement x-component [m]") + axs[1, 0].set_ylabel("Displacement y-component [m]") + axs[2, 0].set_ylabel("Displacement z-component [m]") + axs[2, 0].set_xlabel("Time [s]") + + axs[0, 1].set_title(r"Particle $\mathbf{p}_1$") + axs[2, 1].set_xlabel("Time [s]") + + for data, lbl, color in zip(teams_datasets, labels_names, colors): + for i, u_type in zip([0, 1, 2], ["ux", "uy", "uz"]): + if lbl == "svMultiPhysics": + axs[i, 0].plot( + data["time"], data["displacement"]["p0"][u_type], label=lbl, color="k", linestyle="--" + ) + axs[i, 1].plot( + data["time"], data["displacement"]["p1"][u_type], label=lbl, color="k", linestyle="--" + ) + else: + axs[i, 0].plot( + data["time"], data["displacement"]["p0"][u_type], label=lbl, color=color + ) + axs[i, 1].plot( + data["time"], data["displacement"]["p1"][u_type], label=lbl, color=color + ) + + plt.legend(loc="best", fancybox=True, shadow=True) + fig.tight_layout() + fig.savefig(filename.as_posix(), bbox_inches="tight", dpi=120) + # plt.show() + + +def compute_statistics( + datasets: list[dict[str, np.ndarray]], point: str = "p0" +) -> dict[str, np.ndarray]: + time = datasets[0]["time"] + number_of_element_per_dataset = datasets[0]["time"].shape[0] + number_of_datasets = len(datasets) + condensated_dataset_ux = np.zeros( + (number_of_element_per_dataset, number_of_datasets) + ) + condensated_dataset_uy = np.zeros_like(condensated_dataset_ux) + condensated_dataset_uz = np.zeros_like(condensated_dataset_ux) + red_dataset_u = np.zeros((number_of_datasets,)) + + for i, dataset in enumerate(datasets): + condensated_dataset_ux[:, i] = dataset["displacement"][point]["ux"] + condensated_dataset_uy[:, i] = dataset["displacement"][point]["uy"] + condensated_dataset_uz[:, i] = dataset["displacement"][point]["uz"] + + mean_ux = np.mean(condensated_dataset_ux, axis=1) + mean_uy = np.mean(condensated_dataset_uy, axis=1) + mean_uz = np.mean(condensated_dataset_uz, axis=1) + + for i in range(number_of_datasets): + diff_to_mean_norm = np.sqrt( + np.abs(condensated_dataset_ux[:, i] - mean_ux) ** 2 + + np.abs(condensated_dataset_uy[:, i] - mean_uy) ** 2 + + np.abs(condensated_dataset_uz[:, i] - mean_uz) ** 2 + ) + + mean_norm = np.sqrt( + np.abs(mean_ux) ** 2 + np.abs(mean_uy) ** 2 + np.abs(mean_uz) ** 2 + ) + + red_dataset_u[i] = np.mean(diff_to_mean_norm / mean_norm, axis=0) + + statistics_dataset = { + "time": time, + "mean_ux": mean_ux, + "mean_uy": mean_uy, + "mean_uz": mean_uz, + "std_ux": np.std(condensated_dataset_ux, axis=1), + "std_uy": np.std(condensated_dataset_uy, axis=1), + "std_uz": np.std(condensated_dataset_uz, axis=1), + "red_u": red_dataset_u, + } + + return statistics_dataset + + +def compute_plots_monoventricular_nonblinded_step_1() -> None: + """Computes displacement curves for the monoventricular non-blinded step 1""" + # Get svMultiPhysics data + svmultiphysics_data = get_svmultiphysics_displacement_data() + + # Add svMultiPhysics to the datasets, labels, and colors + all_datasets = TEAMS_DATASETS_1 + [svmultiphysics_data] + all_labels = LABEL_NAMES + ["svMultiPhysics"] + all_colors = COLORS + ["tab:red"] # Add red color for svMultiPhysics + + compute_plot_displacement_monoventricular( + all_datasets, + "./comparison_plots_p0_p1_step_1_nonblinded.png", + all_labels, + all_colors, + ) + + +def compute_statistics_per_dataset( + labels: list[str], datasets: list[dict[str, np.ndarray]], header: str = "" +) -> None: + if len(labels) != len(datasets): + raise ValueError("labels and datasets do not match in length") + + print(header) + statistics_p0 = compute_statistics(datasets, point="p0") + statistics_p1 = compute_statistics(datasets, point="p1") + + red_p0 = np.round(statistics_p0["red_u"], decimals=3) + red_p1 = np.round(statistics_p1["red_u"], decimals=3) + + print(f"Team {'TEAM NAME':25s} - {'p0':4s}, {'p1':4s}") + for label, i in zip(labels, range(len(red_p0))): + print(f"Team {label:25s} - {red_p0[i]:4.3f}, {red_p1[i]:4.3f}") + + print("-----------------------------------\n") + + +if __name__ == "__main__": + compute_plots_monoventricular_nonblinded_step_1() + + # Also plot PV loop separately + print("Plotting pressure-volume loop...") + plot_pv_loop_only() diff --git a/tests/cases/ustruct/LV_HolzapfelOgden_active/generate_pressure.py b/tests/cases/ustruct/LV_HolzapfelOgden_active/generate_pressure.py new file mode 100644 index 000000000..7d3c3a984 --- /dev/null +++ b/tests/cases/ustruct/LV_HolzapfelOgden_active/generate_pressure.py @@ -0,0 +1,131 @@ +import numpy as np +import matplotlib.pyplot as plt +from scipy.integrate import solve_ivp + +params = {} + +params['alpha_min'] = -30 +params['alpha_max'] = 5 +params['alpha_pre'] = 5 +params['alpha_mid'] = 1 +params['sigma_pre'] = 7000 +params['sigma_mid'] = 16000 +params['sigma_0'] = 1.5e5 +params['t_sys_pre'] = 0.17 +params['t_dias_pre'] = 0.484 +params['t_sys'] = 0.16 +params['t_dias'] = 0.484 +params['gamma'] = 0.005 + +def S_plus(delta_t, params): + return 0.5* (1 + np.tanh(delta_t/params['gamma'])) + +def S_minus(delta_t, params): + return 0.5* (1 - np.tanh(delta_t/params['gamma'])) + +def g_pre(t, params): + return S_minus(t - params['t_dias_pre'], params) + +def f(t, params): + return S_plus(t - params['t_sys'], params) * S_minus(t - params['t_dias'], params) + +def f_pre(t, params): + return S_plus(t - params['t_sys_pre'], params) * S_minus(t - params['t_dias_pre'], params) + +def a(t, params): + return params['alpha_max'] * f(t, params) + params['alpha_min'] * (1-f(t, params)) + +def a_pre(t, params): + return params['alpha_max'] * f_pre(t, params) + params['alpha_min'] * (1-f_pre(t, params)) + +def b(t, params): + return a_pre(t, params) + params['alpha_pre'] * g_pre(t, params) + params['alpha_mid'] + +# Differential equation for pressure +def dpdt(t, p): + b_abs = abs(b(t, params)) + b_max = max(b(t, params), 0) + g_pre_max = max(g_pre(t, params), 0) + return -b_abs*p + params['sigma_mid']*b_max + params['sigma_pre']*g_pre_max + +# Differential equation for stress +def dtaudt(t, tau): + a_abs = abs(a(t, params)) + a_max = max(a(t, params), 0) + return -a_abs*tau + params['sigma_0']*a_max + + +def save_sol(sol, filename, f_modes = 512): + # save data to .dat file + # first row is [len(sol.y[0]) f_modes] + # first column is sol.t rounded to six decimal places + # second column is sol.y[0] rounded to six decimal places + + data = np.column_stack((np.round(sol.t, 6), np.round(sol.y[0], 6))) + + data = np.vstack((np.array([len(sol.y[0]), f_modes]), data)) + + # Save to file with the first row as integers and the rest as floats + with open(filename, 'w') as f: + # Write the first row as integers + np.savetxt(f, [data[0]], fmt='%d', delimiter=' ') + # Write the remaining rows as floats with 6 decimal places + np.savetxt(f, data[1:], fmt='%1.6f', delimiter=' ') + +# Solve the ODE for pressure + +p0 = 0 +t_eval = np.linspace(0, 1, 1001) +p_sol = solve_ivp(dpdt, [0, 1], [p0], t_eval=t_eval, method='DOP853') +print('Max pressure:', max(p_sol.y[0])) + +plt.plot(p_sol.t, p_sol.y[0]) +plt.xlabel('Time') +plt.ylabel('Pressure') +''' +# load pressure_original.dat and convert string to float +p_original = np.loadtxt('pressure_original.dat', delimiter=',') +p_original[1:,1] = 0.1*p_original[1:,1] +plt.plot(p_original[1:,0], p_original[1:,1]) + +with open('pressure_original_Pa.dat', 'w') as f: + # Write the first row as integers + np.savetxt(f, [p_original[0]], fmt='%d', delimiter=' ') + # Write the remaining rows as floats with 6 decimal places + np.savetxt(f, p_original[1:], fmt='%1.6f', delimiter=' ') +''' +plt.show() + + +# Save the solution to a file +save_sol(p_sol, 'pressure.dat', 512) + +# Solve the ODE for stress + +tau0 = 0 +t_eval = np.linspace(0, 1, 1001) +tau_sol = solve_ivp(dtaudt, [0, 1], [tau0], t_eval=t_eval) + +print('Max stress:', max(tau_sol.y[0])) + +plt.plot(tau_sol.t, tau_sol.y[0]) +plt.xlabel('Time') +plt.ylabel('Stress') +''' +# load stress_original.dat and convert string to float +tau_original = np.loadtxt('stress_original.dat', delimiter=' ') +tau_original[1:,1] = 0.1*tau_original[1:,1] +plt.plot(tau_original[1:,0], tau_original[1:,1]) + +with open('stress_original_Pa.dat', 'w') as f: + # Write the first row as integers + np.savetxt(f, [tau_original[0]], fmt='%d', delimiter=' ') + # Write the remaining rows as floats with 6 decimal places + np.savetxt(f, tau_original[1:], fmt='%1.6f', delimiter=' ') +''' +plt.show() + +# Save the solution to a file +save_sol(tau_sol, 'LV_stress.dat', 512) + + diff --git a/tests/cases/ustruct/LV_HolzapfelOgden_active/mesh/fibersLong.vtu b/tests/cases/ustruct/LV_HolzapfelOgden_active/mesh/fibersLong.vtu new file mode 100644 index 000000000..28e2d5cf4 --- /dev/null +++ b/tests/cases/ustruct/LV_HolzapfelOgden_active/mesh/fibersLong.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:9c6570939679034cdbba8a2577d5c84baa68339aa8b38da3e2960356bbfcc5d9 +size 868063 diff --git a/tests/cases/ustruct/LV_HolzapfelOgden_active/mesh/fibersSheet.vtu b/tests/cases/ustruct/LV_HolzapfelOgden_active/mesh/fibersSheet.vtu new file mode 100644 index 000000000..425397698 --- /dev/null +++ b/tests/cases/ustruct/LV_HolzapfelOgden_active/mesh/fibersSheet.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1939854b0527c6f78353b12b25aee37c0b82a1202c40e10f4f48d0d2c46ca7a2 +size 868763 diff --git a/tests/cases/ustruct/LV_HolzapfelOgden_active/mesh/mesh-complete.mesh.vtu b/tests/cases/ustruct/LV_HolzapfelOgden_active/mesh/mesh-complete.mesh.vtu new file mode 100644 index 000000000..451bf197a --- /dev/null +++ b/tests/cases/ustruct/LV_HolzapfelOgden_active/mesh/mesh-complete.mesh.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:17f70504f6760f7e836373181fc3487ecc35312bc286fe9b2962e877fe55e517 +size 834394 diff --git a/tests/cases/ustruct/LV_HolzapfelOgden_active/mesh/mesh-surfaces/endo.vtp b/tests/cases/ustruct/LV_HolzapfelOgden_active/mesh/mesh-surfaces/endo.vtp new file mode 100644 index 000000000..8bf3ef696 --- /dev/null +++ b/tests/cases/ustruct/LV_HolzapfelOgden_active/mesh/mesh-surfaces/endo.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:eefb1443d58a90414d3c94b0e6128fc8b3698dc8c25ea33ffeeb573b263de023 +size 68199 diff --git a/tests/cases/ustruct/LV_HolzapfelOgden_active/mesh/mesh-surfaces/epi.vtp b/tests/cases/ustruct/LV_HolzapfelOgden_active/mesh/mesh-surfaces/epi.vtp new file mode 100644 index 000000000..ae626668d --- /dev/null +++ b/tests/cases/ustruct/LV_HolzapfelOgden_active/mesh/mesh-surfaces/epi.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:f08544d5b0cf276d05e999e6116986bea26442506ae51904ce386530a38cf5d1 +size 91341 diff --git a/tests/cases/ustruct/LV_HolzapfelOgden_active/mesh/mesh-surfaces/top.vtp b/tests/cases/ustruct/LV_HolzapfelOgden_active/mesh/mesh-surfaces/top.vtp new file mode 100644 index 000000000..40eb6411c --- /dev/null +++ b/tests/cases/ustruct/LV_HolzapfelOgden_active/mesh/mesh-surfaces/top.vtp @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:685fe50130b4f5b2a8a3f4b66cbd381a1741b4bc0f3b5e46168a21553b32215f +size 8884 diff --git a/tests/cases/ustruct/LV_HolzapfelOgden_active/p-v_loop.png b/tests/cases/ustruct/LV_HolzapfelOgden_active/p-v_loop.png new file mode 100644 index 000000000..4975c4d28 --- /dev/null +++ b/tests/cases/ustruct/LV_HolzapfelOgden_active/p-v_loop.png @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:63bf0c47ba2848c9bad99a5d7e74cd35152e43ea95bb4a44500f35fdb5e0548f +size 187082 diff --git a/tests/cases/ustruct/LV_HolzapfelOgden_active/pressure.dat b/tests/cases/ustruct/LV_HolzapfelOgden_active/pressure.dat new file mode 100644 index 000000000..680a13c97 --- /dev/null +++ b/tests/cases/ustruct/LV_HolzapfelOgden_active/pressure.dat @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:564f6e7c91faa9a4f9234f1d7f4d0ae96997c5b36a072879b6523668dcfd555d +size 20059 diff --git a/tests/cases/ustruct/LV_HolzapfelOgden_active/result_001.vtu b/tests/cases/ustruct/LV_HolzapfelOgden_active/result_001.vtu new file mode 100644 index 000000000..f44f798ec --- /dev/null +++ b/tests/cases/ustruct/LV_HolzapfelOgden_active/result_001.vtu @@ -0,0 +1,3 @@ +version https://git-lfs.github.com/spec/v1 +oid sha256:1adf6fc288b6e02bcc4f629fbb5620c8e96c40775b85fcc7b2bbe2ea919cac4e +size 1915077 diff --git a/tests/cases/ustruct/LV_HolzapfelOgden_active/solver.xml b/tests/cases/ustruct/LV_HolzapfelOgden_active/solver.xml new file mode 100644 index 000000000..24c553250 --- /dev/null +++ b/tests/cases/ustruct/LV_HolzapfelOgden_active/solver.xml @@ -0,0 +1,132 @@ + + + + + 0 + 3 + 1 + 1e-3 + 0.50 + STOP_SIM + + 1 + result + 1 + 1 + + 50 + 0 + + 1 + 1 + 1 + + + + + + + mesh/mesh-complete.mesh.vtu + + + mesh/mesh-surfaces/endo.vtp + + + + mesh/mesh-surfaces/epi.vtp + + + + mesh/mesh-surfaces/top.vtp + + + mesh/fibersLong.vtu + mesh/fibersSheet.vtu + + + + + + + true + 3 + 10 + 1e-9 + + 1000.0 + 1.0e6 + 0.483333 + + + 100.0 + + + + 59.0 + 8.023 + 18472.0 + 16.026 + 2481.0 + 11.12 + 216.0 + 11.436 + 100.0 + + + ST91 + 1e6 + + 1e-5 + 1e-5 + + + LV_stress.dat + false + + + + true + true + true + true + true + true + true + true + true + + + + + fsils + + 1e-12 + 1e-12 + 1000 + 100 + + + + Robin + 1.0e8 + 5.0e3 + 1 + + + + Robin + 1.0e5 + 5.0e3 + + + + Neu + Unsteady + pressure.dat + false + true + + + + + diff --git a/tests/test_struct.py b/tests/test_struct.py index e0785e6be..6ba9db89f 100644 --- a/tests/test_struct.py +++ b/tests/test_struct.py @@ -23,6 +23,9 @@ def test_LV_Guccione_passive(n_proc): test_folder = "LV_Guccione_passive" run_with_reference(base_folder, test_folder, fields, n_proc) +def test_LV_HolzapfelOgden_active(n_proc): + test_folder = "LV_HolzapfelOgden_active" + run_with_reference(base_folder, test_folder, fields, n_proc) def test_LV_HolzapfelOgden_passive(n_proc): test_folder = "LV_HolzapfelOgden_passive" diff --git a/tests/test_ustruct.py b/tests/test_ustruct.py index 905dada21..e5d96993c 100644 --- a/tests/test_ustruct.py +++ b/tests/test_ustruct.py @@ -34,6 +34,10 @@ def test_LV_Guccione_active(n_proc): test_folder = "LV_Guccione_active" run_with_reference(base_folder, test_folder, fields, n_proc) +def test_LV_HolzapfelOgden_active(n_proc): + test_folder = "LV_HolzapfelOgden_active" + run_with_reference(base_folder, test_folder, fields, n_proc) + def test_LV_HolzapfelOgden_passive(n_proc): test_folder = "LV_HolzapfelOgden_passive" run_with_reference(base_folder, test_folder, fields, n_proc)