Module PostProcessing

#-------------------------------------------------------------------------------
# Librairies
#-------------------------------------------------------------------------------

import numpy as np
import matplotlib.pyplot as plt
import pyvista as pv
import pickle, porespy, skimage, time

# Own
from main import index_to_str

#-------------------------------------------------------------------------------
# Functions
#-------------------------------------------------------------------------------

def pp(dict_user):
    '''
    Main function for the post-processing.    
    '''
    # compute performances
    tic = time.perf_counter() 

    # initialize the dict
    dict_tempo = {
        'L_M_Neck': [],
        'L_freeSurface': [],
        'L_M0': [],
        'L_M1': [],
        'L_M2': [],
        'L_M3': [],
        'map': None
    }

    # read and plot the .csv data.
    ReadCSV(dict_user, dict_tempo)

    print('\nread pvtu')
    # determine the zone of interest for the porosity
    # a box is generated from the lastest configuration (more dense)
    # read the .pvtu
    ReadPVTU(dict_user, dict_tempo, 'output/vtk/PF_Sintering_other_'+dict_user['last_j_str']+'.pvtu')     
    # compute the zone of interest
    RebuildMap_ZoneInterest(dict_user, dict_tempo)

    # iterate on the .pvtu files
    for i in range(dict_user['last_j']+1):
        # flag the initial map (for plot)
        if i == 0 :
            dict_tempo['ic_flag'] = True
        else:
            dict_tempo['ic_flag'] = False

        # read the .pvtu
        ReadPVTU(dict_user, dict_tempo, 'output/vtk/PF_Sintering_other_'+index_to_str(i)+'.pvtu')

        # Biswas, 2018
        # compute the neck area
        ComputeNeck(dict_tempo)
        # compute the free surface area
        ComputeFreeSurface(dict_tempo)
    
        # Guevel, 2022
        RebuildMap(dict_user, dict_tempo)
        ComputeMorphometers(dict_tempo)

    # plot final map
    PlotMap(dict_user, dict_tempo)

    # Biswas, 2018
    # plot the neck area
    PlotNeck(dict_tempo)
    # plot the free surface area
    PlotFreeSurface(dict_tempo)

    # Guevel, 2022
    PlotMorphometers(dict_tempo)

    # compute performances
    tac = time.perf_counter()
    hours = (tac-tic)//(60*60)
    minutes = (tac-tic - hours*60*60)//(60)
    seconds = int(tac-tic - hours*60*60 - minutes*60)
    print("\nPost-processing time : "+str(hours)+" hours "+str(minutes)+" minutes "+str(seconds)+" seconds")
    print('Post-processing ends')

#-------------------------------------------------------------------------------

def ReadCSV(dict_user, dict_tempo):
    '''
    Read the csv file generated with Moose (made by postprocessors).
    '''
    print('\nread csv')

    # read file
    f = open('output/PF_Sintering_csv.csv', "r")
    lines = f.readlines()
    f.close()
    # init data
    time_pp = []
    c_pp = []
    L_eta_pp = []
    for i_grain in range(dict_user['n_eta']):
        L_eta_pp.append([])

    # iterate on lines
    for line in lines[1:]:
        line = line.replace("\n", "")
        data = line.split(',')
        # read data
        time_pp.append(float(data[0]))
        c_pp.append(float(data[1]))
        for i_grain in range(dict_user['n_eta']):
            L_eta_pp[i_grain].append(float(data[2+i_grain]))

    # plot time-mass
    fig, ax1 = plt.subplots(1,1,figsize=(16,9))
    ax1.plot(time_pp, c_pp, linewidth=6)
    ax1.set_xlabel('time (-)', fontsize=25)
    ax1.set_ylabel('mean mass concentration (-)', fontsize=25)
    ax1.tick_params(axis='both', labelsize=20, width=3, length=3)
    fig.tight_layout()
    fig.savefig('output/evol_time_mass.png')
    plt.close(fig)
    # output
    print('conservation of the mass:')
    print('mass min:', round(min(c_pp),3), '(delta', int(100*(np.mean(c_pp)-min(c_pp))/np.mean(c_pp)),'% of mean value)')
    print('mass max:', round(max(c_pp),3), '(delta', int(100*(max(c_pp)-np.mean(c_pp))/np.mean(c_pp)),'% of mean value)')

    # plot time-etas
    fig, ax1 = plt.subplots(1,1,figsize=(16,9))
    for i_grain in range(dict_user['n_eta']):
        ax1.plot(time_pp, L_eta_pp[i_grain], linewidth=6)
    ax1.set_xlabel('time (-)', fontsize=25)
    ax1.set_ylabel('mean grain concentration (-)', fontsize=25)
    ax1.tick_params(axis='both', labelsize=20, width=3, length=3)
    fig.tight_layout()
    fig.savefig('output/evol_time_etas.png')
    plt.close(fig)

    # extract time from the iterations saved
    L_time_extracted = []
    for ite in dict_user['L_ite_saved']:
        L_time_extracted.append(time_pp[ite])
    # save
    dict_tempo['L_time_extracted'] = L_time_extracted

#-------------------------------------------------------------------------------

def ReadPVTU(dict_user, dict_tempo, filename):
    '''
    Read the PVTU file generated with Moose.
    '''
    # read the .pvtu file 
    dataset = pv.read(filename)

    # extract the points
    L_points = dataset.points

    # extract the concentration 
    L_c = dataset['c']
    
    # extract the etas
    L_L_eta = []
    for i_grain in range(dict_user['n_eta']):
        L_L_eta.append(dataset['eta'+str(i_grain)])

    # save
    dict_tempo['L_points'] = L_points
    dict_tempo['L_c'] = L_c
    dict_tempo['L_L_eta'] = L_L_eta

#-------------------------------------------------------------------------------

def ComputeNeck(dict_tempo):
    '''
    Compute the neck area.

    Defined as the integral of the function eta_i*eta_j.
    '''
    # initialization
    M_Neck = np.zeros((len(dict_tempo['L_L_eta']), len(dict_tempo['L_L_eta'])))
    # iterate on the pairs
    for i_grain in range(len(dict_tempo['L_L_eta'])-1):
        for j_grain in range(i_grain+1,len(dict_tempo['L_L_eta'])):
            L_eta_i_j = np.array(dict_tempo['L_L_eta'][i_grain])*np.array(dict_tempo['L_L_eta'][j_grain])
            M_Neck[i_grain, j_grain] = np.mean(L_eta_i_j)
    # save
    dict_tempo['L_M_Neck'].append(M_Neck)

#-------------------------------------------------------------------------------

def PlotNeck(dict_tempo):
    '''
    Plot the neck area.

    Defined as the integral of the function eta_i*eta_j.
    '''
    # open figure
    fig, (ax1) = plt.subplots(1,1,figsize=(16,9))
    # iterate on the pairs
    for i_grain in range(len(dict_tempo['L_L_eta'])-1):
        for j_grain in range(i_grain+1,len(dict_tempo['L_L_eta'])):
            L_Neck_i_j = []
            # extract the pair with time
            for ite in range(len(dict_tempo['L_M_Neck'])):
                L_Neck_i_j.append(dict_tempo['L_M_Neck'][ite][i_grain, j_grain])
            # plot
            ax1.plot(dict_tempo['L_time_extracted'], L_Neck_i_j, linewidth=6)
    # close figure
    ax1.set_xlabel('time (-)', fontsize=25)
    ax1.set_ylabel('neck area (-)', fontsize=25)
    ax1.tick_params(axis='both', labelsize=20, width=3, length=3)
    fig.tight_layout()
    fig.savefig('output/evol_time_neck.png')
    plt.close(fig)

#-------------------------------------------------------------------------------

def ComputeFreeSurface(dict_tempo):
    '''
    Compute the free surface area.

    Defined as the integral of the function c_int (=1 for  c in [0.45, 0.55]).
    '''
    # pp data
    L_c_int = []
    for c in dict_tempo['L_c']:
        if 0.45 <= c and c <= 0.55:
            L_c_int.append(1)
        else :
            L_c_int.append(0)
    # save
    dict_tempo['L_freeSurface'].append(np.mean(L_c_int))
    
#-------------------------------------------------------------------------------

def PlotFreeSurface(dict_tempo):
    '''
    Plot the free surface area.

    Defined as the integral of the function c_int (=1 for  c in [0.45, 0.55]).
    '''
    # open figure
    fig, (ax1) = plt.subplots(1,1,figsize=(16,9))
    ax1.plot(dict_tempo['L_time_extracted'], dict_tempo['L_freeSurface'], linewidth=6)
    ax1.set_xlabel('time (-)', fontsize=25)
    ax1.set_ylabel('free surface (-)', fontsize=25)
    ax1.tick_params(axis='both', labelsize=20, width=3, length=3)
    fig.tight_layout()
    fig.savefig('output/evol_time_freesurface.png')
    plt.close(fig)

#-------------------------------------------------------------------------------

def RebuildMap(dict_user, dict_tempo):
    '''
    Rebuild the map from the output of the pyvista.
    '''
    # initialization
    M_void = np.zeros((len(dict_user['L_y']), len(dict_user['L_x'])))
    M_solid = np.zeros((len(dict_user['L_y']), len(dict_user['L_x'])))

    # the map is not know
    if dict_tempo['map'] == None:
        # init map
        map = []
        # iterate on the points
        for i_point in range(len(dict_tempo['L_points'])):
            # search node in the mesh
            L_search = list(abs(np.array(dict_user['L_x']-list(dict_tempo['L_points'][i_point])[0])))
            i_x = L_search.index(min(L_search))
            L_search = list(abs(np.array(dict_user['L_y']-list(dict_tempo['L_points'][i_point])[1])))
            i_y = L_search.index(min(L_search))
            # save map
            map.append([i_x, i_y])
            # rebuild maps
            if dict_tempo['L_c'][i_point] < 0.5:
                M_void[-1-i_y, i_x] = 1
            else:
                M_solid[-1-i_y, i_x] = 1
        # save map
        dict_tempo['map'] = map

    # the map is know
    else:
        # iterate on the points
        for i_point in range(len(dict_tempo['L_points'])):
            # read the map
            i_x = dict_tempo['map'][i_point][0]
            i_y = dict_tempo['map'][i_point][1]
            # rebuild maps
            if dict_tempo['L_c'][i_point] < 0.5:
                M_void[-1-i_y, i_x] = 1
            else:
                M_solid[-1-i_y, i_x] = 1
    
    # extract maps (exclude void on the sides)
    i_x_min = dict_tempo['zoneInterest'][0]
    i_x_max = dict_tempo['zoneInterest'][1]
    i_y_min = dict_tempo['zoneInterest'][2]
    i_y_max = dict_tempo['zoneInterest'][3]

    # save
    dict_tempo['M_void'] = M_void.copy()
    dict_tempo['M_void_extracted'] = M_void.copy()[i_y_min: i_y_max+1, i_x_min: i_x_max+1]
    dict_tempo['M_solid'] = M_solid.copy()
    dict_tempo['M_solid_extracted'] = M_solid.copy()[i_y_min: i_y_max+1, i_x_min: i_x_max+1]

    # save ic 
    if dict_tempo['ic_flag']:
        dict_tempo['M_solid_ic'] = M_solid.copy()   

#-------------------------------------------------------------------------------

def RebuildMap_ZoneInterest(dict_user, dict_tempo):
    '''
    Rebuild the map from the output of the pyvista. And determine the zone of interest for the porosity estimation.
    '''    
    # initialization
    M_solid = np.zeros((len(dict_user['L_y']), len(dict_user['L_x'])))

    # the map is not know
    if dict_tempo['map'] == None:
        # init map
        map = []
        # iterate on the points
        for i_point in range(len(dict_tempo['L_points'])):
            # search node in the mesh
            L_search = list(abs(np.array(dict_user['L_x']-list(dict_tempo['L_points'][i_point])[0])))
            i_x = L_search.index(min(L_search))
            L_search = list(abs(np.array(dict_user['L_y']-list(dict_tempo['L_points'][i_point])[1])))
            i_y = L_search.index(min(L_search))
            # save map
            map.append([i_x, i_y])
            # rebuild maps
            if dict_tempo['L_c'][i_point] > 0.5:
                M_solid[-1-i_y, i_x] = 1
        # save map
        dict_tempo['map'] = map

    # the map is know
    else:
        # iterate on the points
        for i_point in range(len(dict_tempo['L_points'])):
            # read the map
            i_x = dict_tempo['map'][i_point][0]
            i_y = dict_tempo['map'][i_point][1]
            # rebuild maps
            if dict_tempo['L_c'][i_point] > 0.5:
                M_solid[-1-i_y, i_x] = 1
    
    # extract maps (exclude void on the sides)
    # find i_x_min
    #i_x_min = 0
    #while np.max(M_solid[:, i_x_min]) == 0:
    #    i_x_min = i_x_min + 1
    i_x_min = int(0.25*M_solid.shape[0])
    # find i_x_max
    #i_x_max = M_solid.shape[1]-1
    #while np.max(M_solid[:, i_x_max]) == 0:
    #    i_x_max = i_x_max - 1
    i_x_max = int(0.75*M_solid.shape[0])
    # find i_y_min
    #i_y_min = 0
    #while np.max(M_solid[i_y_min, :]) == 0:
    #    i_y_min = i_y_min + 1
    i_y_min = int(0.25*M_solid.shape[1])
    # find i_y_max
    #i_y_max = M_solid.shape[0]-1
    #while np.max(M_solid[i_y_max, :]) == 0:
    #    i_y_max = i_y_max - 1
    i_y_max = int(0.75*M_solid.shape[1])

    # save
    dict_tempo['zoneInterest'] = [i_x_min, i_x_max, i_y_min, i_y_max]

#-------------------------------------------------------------------------------

def PlotMap(dict_user, dict_tempo):
    '''
    plot the current configuration.
    '''
    # Plot maps
    fig, (ax1, ax2) = plt.subplots(1,2,figsize=(16,9))
    # parameters
    ax1.imshow(dict_tempo['M_solid_ic'], interpolation = 'nearest', extent=(dict_user['L_x'][0],dict_user['L_x'][-1],dict_user['L_y'][0],dict_user['L_x'][-1]))
    ax1.set_title('ic', fontsize=30)
    ax2.imshow(dict_tempo['M_solid'], interpolation = 'nearest', extent=(dict_user['L_x'][0],dict_user['L_x'][-1],dict_user['L_y'][0],dict_user['L_x'][-1]))
    ax2.set_title('final', fontsize=30)
    fig.tight_layout()
    fig.savefig('output/map_c.png')
    plt.close(fig)

    # Plot delta
    fig, ax1 = plt.subplots(1,1,figsize=(16,9))
    # parameters
    cax = ax1.imshow(dict_tempo['M_solid']-dict_tempo['M_solid_ic'], interpolation = 'nearest', extent=(dict_user['L_x'][0],dict_user['L_x'][-1],dict_user['L_y'][0],dict_user['L_x'][-1]), cmap='bwr')
    cbar = fig.colorbar(cax, ticks=[1, -1])
    cbar.ax.set_yticklabels(['precipitation', 'dissolution'])
    cbar.ax.tick_params(labelsize=20)
    ax1.set_title('final-ic', fontsize=30)
    fig.tight_layout()
    fig.savefig('output/map_deltac.png')
    plt.close(fig)

#-------------------------------------------------------------------------------

def ComputeMorphometers(dict_tempo):
    '''
    Compute the morphometers.

    see Guével, 2022
    '''
    # M0 porosity
    M0 = porespy.metrics.porosity(dict_tempo['M_void_extracted'])
    
    # M1 perimeter
    M1 = skimage.measure.perimeter(dict_tempo['M_solid'])
    
    # M2 grain size
    # this morphometer requires more lines
    # an algorithm is called to determine the size of the grains
    sizes = porespy.filters.local_thickness(dict_tempo['M_solid'])
    # then a size distribution is called  
    data = porespy.metrics.pore_size_distribution(sizes, log=False, bins=20)
    # finally the mean (weighted by the probability) is computed
    M2 = 0
    counter = 0
    for i_bin in range(len(data.bin_centers)):
        M2 = M2 + data.pdf[i_bin]*data.bin_centers[i_bin]
        counter = counter + data.pdf[i_bin]
    M2 = M2 / counter
    
    # M3 Euler
    M3 = skimage.measure.euler_number(dict_tempo['M_solid'])

    # save
    dict_tempo['L_M0'].append(M0)
    dict_tempo['L_M1'].append(M1)
    dict_tempo['L_M2'].append(M2)
    dict_tempo['L_M3'].append(M3)

#-------------------------------------------------------------------------------

def PlotMorphometers(dict_tempo):
    '''
    Plot the morphometers.

    see Guevel, 2022.
    '''
    # open figure
    fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2,2,figsize=(2*16,2*9))
    # M0
    ax1.plot(dict_tempo['L_time_extracted'], dict_tempo['L_M0'], linewidth=6)
    ax1.set_xlabel('time (-)', fontsize=25)
    ax1.set_ylabel('M0 porosity (-)', fontsize=25)
    ax1.tick_params(axis='both', labelsize=20, width=3, length=3)
    # M1
    ax2.plot(dict_tempo['L_time_extracted'], dict_tempo['L_M1'], linewidth=6)
    ax2.set_xlabel('time (-)', fontsize=25)
    ax2.set_ylabel('M1 perimeter (-)', fontsize=25)
    ax2.tick_params(axis='both', labelsize=20, width=3, length=3)
    # M2
    ax3.plot(dict_tempo['L_time_extracted'], dict_tempo['L_M2'], linewidth=6)
    ax3.set_xlabel('time (-)', fontsize=25)
    ax3.set_ylabel('M2 grain size (-)', fontsize=25)
    ax3.tick_params(axis='both', labelsize=20, width=3, length=3)
    # M3
    ax4.plot(dict_tempo['L_time_extracted'], dict_tempo['L_M3'], linewidth=6)
    ax4.set_xlabel('time (-)', fontsize=25)
    ax4.set_ylabel('M3 euler (-)', fontsize=25)
    ax4.tick_params(axis='both', labelsize=20, width=3, length=3)
    # close
    fig.tight_layout()
    fig.savefig('output/evol_time_morphometers.png')
    plt.close(fig)

#-------------------------------------------------------------------------------
# MAIN code
#-------------------------------------------------------------------------------

if __name__ == '__main__':

    # load dict_user
    with open('output/dict_user', 'rb') as handle:
        dict_user = pickle.load(handle)

    # call pp function
    pp(dict_user)