diff --git a/functions.py b/functions.py index 74d5e5a..5577723 100644 --- a/functions.py +++ b/functions.py @@ -5,14 +5,49 @@ Created on Thu Jul 16 15:08:27 2020 @author: Aloma Blanch """ - import numpy as np -import matplotlib.pyplot as plt import pandas as pd -from statistics import mean -# from scipy.signal import find_peaks +import matplotlib.pyplot as plt import matplotlib.transforms as mtransforms +from math import floor +from statistics import mean +from scipy.signal import find_peaks + +def cycle(folder,dt,save_path): + pressure = np.loadtxt(folder + '/PHistRCR.dat',skiprows=2) + N_ts = pressure.shape[0] + T = round(dt*N_ts,3) + time = np.linspace(0,T,N_ts) + # Find peaks, keep only the maximums + peaks, _ = find_peaks(pressure[:,-1]) + idx = np.where(pressure[peaks,-1] >= np.mean(pressure)) + peaks_loc = peaks[idx] + P_peaks = pressure[peaks,-1][pressure[peaks,-1] >= np.mean(pressure)] + # Rmove the multiple picks found in the same location + idx_del = [] + t_pk_loc = time[peaks_loc].tolist() + peak_tdiff = [t_pk_loc[n]-t_pk_loc[n-1] for n in range(1,len(t_pk_loc))] + for i in range(0,len(peak_tdiff)): + if peak_tdiff[i]<= dt*10: + idx_del.append(i) + t_pk_loc[i] = [] + peak_tdiff[i] = [] + t_pk_loc[:] = [x for x in t_pk_loc if x] + peak_tdiff[:] = [x for x in peak_tdiff if x] + P_peaks = np.delete(P_peaks,idx_del) + + fig, ax = plt.subplots() + ax.plot(time,pressure[:,-1]/1333.22,'b') + ax.plot(t_pk_loc, P_peaks/1333.22, "or",label='peak location') + ax.set(xlabel='time [s]', ylabel='Pressure [mmHg]', + title='Pressure at last outlet') + ax.spines['right'].set_visible(False) + ax.spines['top'].set_visible(False) + ax.legend(loc=0) + # plt.show() + plt.savefig(save_path + '/Pressure_max_n_cycles.pdf') + return (floor(mean(peak_tdiff)*1000)/1000,len(t_pk_loc)) def error_plot(folder,t_step,r_criteria,save_path): # Load iterations and residual error diff --git a/main.py b/main.py index 786845d..0105028 100644 --- a/main.py +++ b/main.py @@ -20,7 +20,7 @@ import os.path # from fpdf import FPDF -from functions import error_plot, periodicity, pressure, flow, inlet_flow_waveform, rcr +from functions import error_plot, periodicity, pressure, flow, inlet_flow_waveform, rcr, cycle from generatePDF import generatePDF @@ -33,11 +33,6 @@ save_path = folder+'/'+project+'-report' os.mkdir(save_path) save_pdf = save_path + '/' + project + '-report.pdf' -# Input parameters -T_cyc = 0.477 -n_cyc = 5 - - # Read important parameters from - solver.inp file mylines = [] # Declare an empty list named mylines. with open (project_folder + '/solver.inp', 'rt') as myfile: # Open lorem.txt for reading text data. @@ -55,6 +50,9 @@ t_btw_rst = int(mylines[6][37:-1]) # Extracting Outlet Boundary Conditions from - rcrt.dat fle Rc_C_Rd = rcr(project_folder) + +# Extracting number of cycles and period of cardiac cycle +(T_cyc,n_cyc) = cycle(folder,dt,save_path) # Cehcking convergency and periodicity error_plot(folder,dt,rc,save_path)