arnau_op/post_Process.py

 #!/usr/bin/env python
import numpy as np
import matplotlib.pyplot as plt

T_cyc = 0.477
N_ts = 14280
cyc = 5

pressure_filename = "Z:\\Aloma\\Class\\final_BIEN6931\\Simulations\\MCS_m04\\96-procs_case\\PHistRCR.dat"
#pressure_filename = "Z:\\Aloma\\Class\\final_BIEN6931\\Simulations\\sim03_RCR_2cyc\\PHistRCR_ok.dat"
#pressure = np.loadtxt(pressure_filename)
pressure = np.loadtxt(pressure_filename, skiprows=1,)

time = np.linspace(0,T_cyc*cyc,pressure.shape[0])

# Ploting Pressure over iterations
# use this to have the plot in different window %matplotlib qt
n=2
plt.figure()
plt.plot(time,pressure[:,n]/1333.22)
plt.ylabel('Pressure [mmHg]')
plt.xlabel('Time [s]')
plt.title('Presure evolution ROI-2')
axes = plt.gca()
xmax = time[-1]
xmin = time[0]
axes.set_xlim([xmin,xmax])
ymin = min(pressure[:,n]/1333.22)
ymax = max(pressure[:,n]/1333.22+1)
axes.set_ylim([ymin,ymax])
plt.show()

# Ploting Pressure in all the vessels same plot over time
n = pressure.shape
i=0
plt.figure()
while (i<n[1]):
    plt.plot(time[round(N_ts-N_ts/cyc):N_ts],pressure[round(N_ts-N_ts/cyc):N_ts,i]/1333.22, label='ROI-'+str(i+2))
    i=i+1
    
plt.ylabel('Pressure [mmHg]')
plt.xlabel('Time [s]')
#plt.title('5th Cardiac Cylce')
axes = plt.gca()
#xmax = time[-1]
#xmin = time[round(N_ts-N_ts/cyc)]
#axes.set_xlim([xmin,xmax])
#ymin = np.amin(pressure[round(N_ts-N_ts/cyc):N_ts,:]/1333.22)
#ymax = np.amax(pressure[round(N_ts-N_ts/cyc):N_ts,:]/1333.22)+1
#axes.set_ylim([ymin,ymax])
plt.legend(loc='upper right')
plt.show()


# Ploting Pressure over time
i=0
while (i<n[1]):
    plt.plot(time,pressure[:,i]/1333.22)
    #plt.plot(round(N_ts/cyc), flowRate[round(N_ts/cyc),i], 'ro')
    plt.ylabel('Pressure [mmHg]')
    plt.xlabel('Time [s]')
    plt.title('Pressure waveform ROI-'+str(i+2)+' [cgs]')
    #axes = plt.gca()
    xmax = time[-1]
    xmin = time[0]
    axes.set_xlim([xmin,xmax])
    ymin = min(pressure[:,0]/1333.22)
    ymax = max(pressure[:,0]/1333.22)
    axes.set_ylim([ymin,ymax])
    plt.show()
    i = i + 1


# Plot Flow Rate to analyze that inlet is as expected all outlets all cycles
flowRate_filename = "Z:\\Aloma\\Class\\final_BIEN6931\\Simulations\\MCS_m04\\96-procs_case\\QHistRCR.dat"
#flowRate_filename = "Z:\\Aloma\\Class\\final_BIEN6931\\Simulations\\sim03_RCR_cyc\\QHistRCR_ok.dat"
#flowRate = np.loadtxt(flowRate_filename)
flowRate = np.loadtxt(flowRate_filename, skiprows=1,)

n = flowRate.shape
i=0

while (i<n[1]):
    plt.plot(time,flowRate[:,i])
    #plt.plot(round(N_ts/cyc), flowRate[round(N_ts/cyc),i], 'ro')
    plt.ylabel('Flow Rate Q [mL/s]')
    plt.xlabel('Time [s]')
    plt.title('Outlet Flow Rate ROI-'+str(i+2)+ ' [cgs]')
    #axes = plt.gca()
    xmax = time[-1]
    xmin = time[0]
    axes.set_xlim([xmin,xmax])
    ymin = min(flowRate[:,0])
    ymax = max(flowRate[:,0])
    axes.set_ylim([ymin,ymax])
    plt.show()
    i = i + 1
    
    
# Modifications need to be done in order to Plot only the cycle that user's want
cycle = int(input("Which cycle you whant to plot (velocity plot)? "))
n = flowRate.shape
i=0
while (i<n[1]):
    plt.plot(time[round(N_ts-N_ts/cyc):N_ts],flowRate[round(N_ts-N_ts/cyc):N_ts,i])
    i = i + 1

plt.ylabel('Flow Rate Q [mL/s]')
plt.xlabel('Time [s]')
plt.title('Out Flow Rate cylce '+str(cycle)+ '- all ROI [cgs]')
#axes = plt.gca() 
xmax = time[-1]
xmin = time[0]
axes.set_xlim([xmin,xmax])
ymin = min(flowRate[:,0])
ymax = max(flowRate[:,0])
axes.set_ylim([ymin,ymax])
plt.show()