aloma/SimVascular_report_PostProcess
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Creating the PDF output

Browse Source
This commit is contained in:
Aloma Blanch 2020-07-16 23:33:48 -05:00
parent e079927b84
commit 1e4ef53bc6
2 changed files with 82 additions and 42 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

88
functions.py
View File

@ -12,7 +12,7 @@ from scipy.signal import find_peaks
import matplotlib.transforms as mtransforms import matplotlib.transforms as mtransforms
def error_plot(folder,t_step,r_criteria,save): def error_plot(folder,t_step,r_criteria,save_path):
# Load iterations and residual error # Load iterations and residual error
histor = folder + '/histor.dat' histor = folder + '/histor.dat'
input = open(histor, 'r') input = open(histor, 'r')
@ -29,30 +29,31 @@ def error_plot(folder,t_step,r_criteria,save):
error.append(error_info[2][n-1]) error.append(error_info[2][n-1])
time = np.linspace(start = t_step, stop = len(error)*t_step, num = len(error)) time = np.linspace(start = t_step, stop = len(error)*t_step, num = len(error))
# Plots of interest
# Liniar Scale # Liniar Scale
plt.figure() fig, ax = plt.subplots()
plt.plot(time,error) ax.plot(time,error)
plt.plot(time,r_criteria*np.ones(len(error)),'r') ax.plot(time,r_criteria*np.ones(len(error)),'r')
plt.ylabel('Residual error') ax.set(xlabel='Time steps', ylabel='Residual error',
plt.xlabel('Time steps') title='Last nonlinear residual error for each time step')
plt.title('Last nonlinear residual error for each time step') ax.spines['right'].set_visible(False)
plt.grid(True) ax.spines['top'].set_visible(False)
if save: plt.savefig(plt_folder + case + '_Last_nonlin_res_error.pdf') plt.show()
plt.savefig(save_path + '/Last_nonlin_res_error.pdf')
# Semilog scale # Semilog scale
plt.figure() fig, ax = plt.subplots()
plt.semilogy(time,error) ax.semilogy(time,error)
plt.semilogy(time,r_criteria*np.ones(len(error)),'r') ax.semilogy(time,r_criteria*np.ones(len(error)),'r')
plt.ylabel('Residual error') ax.set(xlabel='Time steps', ylabel='Residual error',
plt.xlabel('Time steps') title='Log - Last nonlinear residual error for each time step')
plt.title('Log - Last nonlinear residual error for each time step') ax.spines['right'].set_visible(False)
plt.grid(True) ax.spines['top'].set_visible(False)
if save: plt.savefig(plt_folder + case + '_Log_Last_nonlin_res_error.pdf') plt.show()
plt.savefig(save_path + '/Log_Last_nonlin_res_error.pdf')
fig.savefig(save_path + '/Log_Last_nonlin_res_error.jpg',dpi=150)
def periodicity(project,folder,dt,T_cyc,n_cyc):
def periodicity(project,folder,dt,T_cyc,n_cyc,save_path):
pressure = np.loadtxt(folder+'/PHistRCR.dat',skiprows=2,) pressure = np.loadtxt(folder+'/PHistRCR.dat',skiprows=2,)
time = np.linspace(0,T_cyc*n_cyc,round(T_cyc/dt*n_cyc)) time = np.linspace(0,T_cyc*n_cyc,round(T_cyc/dt*n_cyc))
peak_P = [] peak_P = []
@ -66,19 +67,23 @@ def periodicity(project,folder,dt,T_cyc,n_cyc):
peak_Pdiff = list(map(abs, peak_Pdiff)) peak_Pdiff = list(map(abs, peak_Pdiff))
fig, ax = plt.subplots() fig, ax = plt.subplots()
ax.plot(time,pressure[:,-1]/1333.22) ax.plot(time,pressure[:,-1]/1333.22,'b')
ax.plot(time[peak_P_pos], peak_P,'ro',label='Cylce pike') ax.plot(time[peak_P_pos], peak_P,'ro',label='Cylce pike')
ax.set(xlabel='time [s]', ylabel='Pressure [mmHg]', ax.set(xlabel='time [s]', ylabel='Pressure [mmHg]',
title='Pressure @ last outlet') title='Pressure at last outlet')
ax.spines['right'].set_visible(False) ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False) ax.spines['top'].set_visible(False)
ax.legend(loc=0) ax.legend(loc=0)
plt.show() plt.show()
plt.savefig(save_path + '/periodicity.pdf')
fig.savefig(save_path + '/periodicity.jpg',dpi=150)
if (peak_Pdiff[-1]<=1): print('The numerical simulation \'{0}\' has achieve periodicity!\nSystolic Blood Pressure (SBP):\nsecond-last cycle = {1:.2f} mmHg,\nlast cycle = {2:.2f} mmHg,\n\u0394mmHg = {3:.2f} mmHg'.format(project,peak_P[-2],peak_P[-1],peak_Pdiff[-1])) if (peak_Pdiff[-1]<=1):
print('The numerical simulation \'{0}\' has achieve periodicity!\nSystolic Blood Pressure (SBP):\nsecond-last cycle = {1:.2f} mmHg,\nlast cycle = {2:.2f} mmHg,\n\u0394mmHg = {3:.2f} mmHg'.format(project,peak_P[-2],peak_P[-1],peak_Pdiff[-1]))
txt = 'The numerical simulation \'{0}\' has achieve periodicity!\nSystolic Blood Pressure (SBP):\nsecond-last cycle = {1:.2f} mmHg,\nlast cycle = {2:.2f} mmHg,\nDelta_mmHg = {3:.2f} mmHg'.format(project,peak_P[-2],peak_P[-1],peak_Pdiff[-1])
def pressure(folder,N_ts,T_cyc,dt,n_cyc): return txt
def pressure(folder,N_ts,T_cyc,dt,n_cyc,save_path):
pressure = np.loadtxt(folder+'/PHistRCR.dat',skiprows=2,) pressure = np.loadtxt(folder+'/PHistRCR.dat',skiprows=2,)
Nc = round(T_cyc/dt) Nc = round(T_cyc/dt)
time = np.linspace(0,T_cyc,Nc) time = np.linspace(0,T_cyc,Nc)
@ -93,15 +98,17 @@ def pressure(folder,N_ts,T_cyc,dt,n_cyc):
MBP[i] = (mean(pressure[N_ts-Nc:N_ts,i]/1333.22)) MBP[i] = (mean(pressure[N_ts-Nc:N_ts,i]/1333.22))
PP = SBP-DBP PP = SBP-DBP
ax.set(xlabel='time [s]', ylabel='Pressure [mmHg]', ax.set(xlabel='time [s]', ylabel='Pressure [mmHg]',
title='Pressure @ each outlet') title='Pressure at each outlet')
ax.spines['right'].set_visible(False) ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False) ax.spines['top'].set_visible(False)
ax.legend(loc=0) ax.legend(loc=0)
plt.show() plt.show()
plt.savefig(save_path + '/pressure.pdf')
fig.savefig(save_path + '/pressure.jpg',dpi=150)
return (DBP,MBP,SBP,PP) return (DBP,MBP,SBP,PP)
def flow(folder,N_ts,T_cyc,dt,n_cyc): def flow(folder,N_ts,T_cyc,dt,n_cyc,save_path):
flow = np.loadtxt(folder+'/QHistRCR.dat',skiprows=2,) flow = np.loadtxt(folder+'/QHistRCR.dat',skiprows=2,)
Nc = round(T_cyc/dt) Nc = round(T_cyc/dt)
time = np.linspace(0,T_cyc,Nc) time = np.linspace(0,T_cyc,Nc)
@ -112,14 +119,16 @@ def flow(folder,N_ts,T_cyc,dt,n_cyc):
Q[i] = (mean(flow[N_ts-Nc:N_ts,i])) Q[i] = (mean(flow[N_ts-Nc:N_ts,i]))
ax.set(xlabel='time [s]', ylabel='Flow [mL/s]', ax.set(xlabel='time [s]', ylabel='Flow [mL/s]',
title='Flow @ each outlet') title='Flow at each outlet')
ax.spines['right'].set_visible(False) ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False) ax.spines['top'].set_visible(False)
ax.legend(loc=0) ax.legend(loc=0)
plt.show() plt.show()
plt.savefig(save_path + '/flow.pdf')
fig.savefig(save_path + '/flow.jpg',dpi=150)
return Q return Q
def inlet_flow_waveform(project_folder,t_btw_rst,N_ts,dt,T_cyc,n_cyc): def inlet_flow_waveform(project_folder,t_btw_rst,N_ts,dt,T_cyc,n_cyc,save_path):
x = np.loadtxt(project_folder+'/ROI-1.flow') x = np.loadtxt(project_folder+'/ROI-1.flow')
t = x[:,0] t = x[:,0]
Q = -x[:,1] Q = -x[:,1]
@ -136,12 +145,11 @@ def inlet_flow_waveform(project_folder,t_btw_rst,N_ts,dt,T_cyc,n_cyc):
fig, ax = plt.subplots() fig, ax = plt.subplots()
ax.plot(t, Q, 'r') ax.plot(t, Q, 'r')
ax.plot(t_pts, Q_pts, 'ob') ax.plot(t_pts, Q_pts, 'ob',label='Time steps saved')
trans_offset = mtransforms.offset_copy(ax.transData, fig=fig, trans_offset = mtransforms.offset_copy(ax.transData, fig=fig,
x=-0, y=0.15, units='inches') x=-0, y=0.15, units='inches')
ax.set(xlabel='Time [s]', ylabel='Flow Rate - Q [mL/s]', ax.set(xlabel='Time [s]', ylabel='Flow Rate - Q [mL/s]',
title='Inlet Flow rate Waveform - 1 cycle') title='Inlet Flow Rate Waveform')
ax.set_ylim([-10, 90]) ax.set_ylim([-10, 90])
ax.spines['right'].set_visible(False) ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False) ax.spines['top'].set_visible(False)
@ -151,8 +159,14 @@ def inlet_flow_waveform(project_folder,t_btw_rst,N_ts,dt,T_cyc,n_cyc):
for i in range(0,np.unique(np.round(t_pts,3)).shape[0]): for i in range(0,np.unique(np.round(t_pts,3)).shape[0]):
time.append('$t_'+str(i+1)+'$') time.append('$t_'+str(i+1)+'$')
for x, y, t in zip(t_pts[(-n_cyc-1):], Q_pts[(-n_cyc-1):], time): for x, y, t in zip(t_pts[(-n_cyc-1):], Q_pts[(-n_cyc-1):], time):
plt.text(x, y, t, transform=trans_offset, fontsize=12) plt.text(x, y, t, transform=trans_offset, fontsize=12)
ax.legend(loc=0)
plt.show() plt.show()
plt.savefig(save_path + '/inlet_waveform.pdf')
fig.savefig(save_path + '/inlet_waveform.jpg',dpi=150)
print('{0} time steps saved, available to visualize in ParaView.'.format(np.unique(np.round(t_pts,3)).shape[0]))
txt = '{0} time steps saved, available to visualize in ParaView.'.format(np.unique(np.round(t_pts,3)).shape[0])
return txt

36
main.py
View File

@ -11,6 +11,7 @@ import os.path
from scipy.signal import find_peaks_cwt from scipy.signal import find_peaks_cwt
from scipy import signal from scipy import signal
import statistics import statistics
from fpdf import FPDF
from functions import error_plot, periodicity, pressure, flow, inlet_flow_waveform from functions import error_plot, periodicity, pressure, flow, inlet_flow_waveform
@ -20,6 +21,9 @@ Tk().withdraw()
folder = askdirectory() folder = askdirectory()
project_folder = os.path.dirname(folder) project_folder = os.path.dirname(folder)
project = os.path.basename(project_folder) project = os.path.basename(project_folder)
save_path = folder+'/'+project+'-report'
# os.mkdir(save_path)
save_pdf = save_path + '/' + project + '-report.pdf'
# Input parameters # Input parameters
T_cyc = 0.477 T_cyc = 0.477
@ -42,14 +46,36 @@ rc = float(mylines[26][18:-1])
t_btw_rst = int(mylines[6][37:-1]) t_btw_rst = int(mylines[6][37:-1])
# Cehcking convergency and periodicity # Cehcking convergency and periodicity
error_plot(folder,dt,rc,False) error_plot(folder,dt,rc,save_path)
periodicity(project,folder,dt,T_cyc,n_cyc) txt1 = periodicity(project,folder,dt,T_cyc,n_cyc,save_path)
# Pressure - Outlets # Pressure - Outlets
(DBP,MBP,SBP,PP) = pressure(folder,N_ts,T_cyc,dt,n_cyc) (DBP,MBP,SBP,PP) = pressure(folder,N_ts,T_cyc,dt,n_cyc,save_path)
# Flow Rate - Outlets # Flow Rate - Outlets
(Q_avg) = flow(folder,N_ts,T_cyc,dt,n_cyc) (Q_avg) = flow(folder,N_ts,T_cyc,dt,n_cyc,save_path)
# Inlet Flow Rate + and t saved # Inlet Flow Rate + and t saved
inlet_flow_waveform(project_folder,t_btw_rst,N_ts,dt,T_cyc,n_cyc) txt2 = inlet_flow_waveform(project_folder,t_btw_rst,N_ts,dt,T_cyc,n_cyc,save_path)
# Create PDF report
pdf = FPDF('P','mm','Letter')
pdf.add_page()
pdf.set_font('Times', 'B', 18)
pdf.cell(200, 20, 'Project name: '+ project, 0, 1, 'C')
# Keep this ratio - 5.42/3.86
pdf.set_font('Times', '', 14)
pdf.cell(200, 30, 'Cehcking convergency and periodicity:', 0, 1, 'L')
pdf.image(save_path +'/Log_Last_nonlin_res_error.jpg', x = None, y = None, w = 140, h = 100, type = '', link = '')
pdf.image(save_path +'/periodicity.jpg', x = None, y = None, w = 140, h = 100, type = '', link = '')
pdf.cell(200, 20, 'Cehcking Pressures at each outlet:', 0, 1, 'L')
pdf.image(save_path +'/pressure.jpg', x = None, y = None, w = 140, h = 100, type = '', link = '')
pdf.cell(200, 20, 'Cehcking Flow Rate at each outlet:', 0, 1, 'L')
pdf.image(save_path +'/flow.jpg', x = None, y = None, w = 140, h = 100, type = '', link = '')
pdf.image(save_path +'/inlet_waveform.jpg', x = None, y = None, w = 140, h = 100, type = '', link = '')
pdf.output(save_pdf, 'F')