###
# Run_model.py (C) of thesis research "Modeling and Simulating Three Phases of Steel: Austenite, Feritte and Cementite" 2018 by Thijs Verbeek
# is licensed under a Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License.
#
# You should have received a copy of the license along with this
# work. If not, see <http://creativecommons.org/licenses/by-nc-sa/4.0/>.
###
 
if __name__ == "__main__":
    from copy import deepcopy
    import sys
    import scipy.io as sio
    import pickle
    import traceback
    import numpy as np
    import scipy.linalg as la
    import matplotlib.pyplot as plt
#    from matplotlib.patches import Ellipse
    from set_layout_plot import plot_settings
    plot_settings()

    from Grid import cGrid2D
    import do_time_step2D

    import parameters as param
    import equilibrium2D

    import time
    from datetime import timedelta


    def H(grid,phi_x):
        BP_points = np.where(grid.is_point_BP[grid.diffusion_mesh_gamma.points_index,0] == 1)[0]
        sigma = np.abs(np.max(phi_x[BP_points]))/2
        H = np.zeros(len(phi_x))
        where_smaller_sigma = np.where(np.abs(- phi_x - sigma) <= sigma)[0]
        H[where_smaller_sigma] = (1/2)*(1 - np.sin((phi_x[where_smaller_sigma] - sigma)*np.pi/(2*sigma)))
        H[np.where(phi_x < -2*sigma)[0]] = 0
        return H
    plt.close('all')


    prefix_dim = 0
    x0 = 0*10**(-prefix_dim)
    xN = 25*10**(-prefix_dim)
    y0 = 0*10**(-prefix_dim)
    yN = 25*10**(-prefix_dim)
    d=2
    t=0
    t_end = param.t_end() - t

    T = param.T0()
    interface_gammaalpha0 = -np.ones(2)*10**(-prefix_dim)
    interface_alphatheta0 =  10+np.array([2.5, 0])*10**(-prefix_dim)

    x03 = 10+np.array([4.0, 4.6]) * 10**(-prefix_dim)
    x13 = 10+np.array([1, 4.6]) * 10**(-prefix_dim) #cause of symmetry not needed
    t0_t = 5*np.pi/12
    t1_t = np.pi - t0_t
    a_t = (x03[0] - interface_alphatheta0[0])/np.cos(t0_t)
    b_t = (x03[1] - interface_alphatheta0[1])/np.sin(t0_t)
    omega_t0 = np.arctan(b_t*np.cos(t0_t)/(-a_t*np.sin(t0_t)))
    omega = np.arccos((1.05 - .94)/.39)

    interface_gammaalpha0[0] = interface_alphatheta0[0]

#    t0_ga = np.arctan(-np.cos(omega_t0 + omega)/(x03[0] - interface_gammaalpha0[0]))
#    t1_ga = np.arctan(-np.cos(omega_t0 + omega)/(x03[0] - interface_gammaalpha0[0]))
#    a_ga = ((x03[0] - interface_gammaalpha0[0])/np.cos(t0_ga))
#    interface_gammaalpha0[1] = x03[1] - np.sin(omega + omega_t0) * np.tan(t0_ga)
#    b_ga = ((x03[1] - interface_gammaalpha0[1])/np.sin(t0_ga))
#    t1_ga = np.arccos((x13[0]-interface_gammaalpha0[0] )/a_ga)

    a_ga = 1.39 * (x03[0] - interface_alphatheta0[0])
    t0_ga = np.arccos((x03[0] - interface_gammaalpha0[0])/a_ga)
    b_ga = - a_ga * np.tan(np.pi - omega - omega_t0) * np.tan(t0_ga)
#    b_ga = a_ga
    interface_gammaalpha0[1] = x03[1] - b_ga * np.sin(t0_ga)



    N = 7**2
    ext = 'Interpolate' #'L2' 'Interpolate' 'Shift'

    area = (xN-x0) * (yN - y0) *(1e-6)**2

### Loading saved time steps
#    name = 'gridhmax025_ind_shift_99dtk1_ind_shift_no_reinit'
#    pkl_file = open('Results/Triple_first/'+name+'.pkl', 'rb')
#    data = pickle.load(pkl_file)
#    grid = data['grid']
#    c_gamma = data['c_gamma']
#    c_alpha = data['c_alpha']
#    dcdn_alpha_ga = data['dcdn_alpha_ga']
#    t = data['t']
#    t_old = data['t_old']
#    m_C_gamma = data['m_C_gamma']
#    m_C_alpha = data['m_C_alpha']
#    m_C_theta = data['m_C_theta']
#    area_gamma = data['area_gamma']
#    area_alpha = data['area_alpha']
#    area_theta = data['area_theta']
#    t_vec = data['t_vec']
#    loaded = 1
#    m_C0 = m_C_gamma[0] + m_C_alpha[0] + m_C_theta[0]
#    pkl_file.close()
#
#    grid.two_particles = None


## Starting from t0, mesh loaded from already made mesh (.pkl file) 
    name = 'gridhmax025'
    pkl_file = open('Results/'+name+'.pkl', 'rb')
    data = pickle.load(pkl_file)
    grid = data['grid']
    loaded = 0
    pkl_file.close()
    t_old = 0

### Starting from t0, mesh still has to be made from matlab mesh (.mat file)
#    name = 'gridhmax025'
#    loaded = 0
#    mesh = sio.loadmat('Grid/gridhmax025.mat')  #gridhmax1.mat 05 025 0125 00625 003125                        #gridN102.mat 200 401 801
#    start_time = time.monotonic()
#    ## uniform (Eqidistant with small offset (1e-5)) random (totally random points in given domain) regular(matlab using Delaunay) 
#    grid = cGrid2D.cGrid(N,x0,x0,xN, xN,'uniform',d, interface_gammaalpha0, interface_alphatheta0, mesh) #uniform random regular(matlab)
#    end_time = time.monotonic()
#    print('Time taken:',timedelta(seconds=end_time - start_time),'\n')
#    t_old = 0
#    
#    
### Saving it to .pkl file
#    data = {'grid': grid}
#    output = open('./Results/'+name+'.pkl', 'wb')
#    pickle.dump(data, output)
#    output.close()

### Plot T_E, T and interface edges in T_E (at = red, gt = green, ga = blue, ga_ghost = lightblue)
    f, ax2 = plt.subplots(1,1)
    ax2.triplot(grid.enriched_mesh.points[:,0],grid.enriched_mesh.points[:,1], grid.enriched_mesh.triangles_points_index, linewidth = 1, color = 'black', alpha = 1)
    ax2.triplot(grid.background_mesh.points[:,0],grid.background_mesh.points[:,1], grid.background_mesh.triangles_points_index, linewidth = 1, color = 'black', alpha = .3)
    [ax2.plot(grid.enriched_mesh.points[grid.interface_alphatheta.edges_index[i,:],0],grid.enriched_mesh.points[grid.interface_alphatheta.edges_index[i,:],1], linewidth = 2, color = 'red') for i in range(0,len(grid.interface_alphatheta.edges))]
    [ax2.plot(grid.enriched_mesh.points[grid.interface_gammatheta.edges_index[i,:],0],grid.enriched_mesh.points[grid.interface_gammatheta.edges_index[i,:],1], linewidth = 2, color = 'green') for i in range(0,len(grid.interface_gammatheta.edges))]
    [ax2.plot(grid.enriched_mesh.points[grid.interface_gammaalpha.edges_index[i,:],0],grid.enriched_mesh.points[grid.interface_gammaalpha.edges_index[i,:],1], linewidth = 2, color = 'blue') for i in range(0,len(grid.interface_gammaalpha.edges))]
    [ax2.plot(grid.enriched_mesh.points[grid.interface_gammaalpha.ghost_edges_index[i,:],0],grid.enriched_mesh.points[grid.interface_gammaalpha.ghost_edges_index[i,:],1], linewidth = 2, color = 'lightblue') for i in range(0,len(grid.interface_gammaalpha.ghost_edges_index))]
    #Real ellipses as comparison
#    ellipse_at = Ellipse(xy=(grid.interface_alphatheta0[0], grid.interface_alphatheta0[1]), width=2*a_t, height=2*b_t,
#                            edgecolor='red', fc='None', lw=1.5)
#    ax2.add_patch(ellipse_at)
#    ellipse_ga = Ellipse(xy=(grid.interface_gammaalpha0[0], grid.interface_gammaalpha0[1]), width=2*a_ga, height=2*b_ga,
#                            edgecolor='blue', fc='None', lw=1.5)
#    ax2.add_patch(ellipse_ga)


    if loaded == 0:
        ### Calculating initial solution        
        H_ga = H(grid, grid.level_set_gammaalpha_enriched.phi[grid.diffusion_mesh_gamma.points_index])
        H_t  = H(grid, grid.level_set_alphatheta_enriched.phi[grid.diffusion_mesh_gamma.points_index])
        H_gat = (np.max(np.array([H_ga, H_t]),axis = 0))

        c_alpha = equilibrium2D.equilibrium_alpha(grid, T)
        c0_gamma = param.c_gammatheta(T) #(param.c_gammatheta(T) + param.c_gammaalpha(T))/2 #
        c_gamma_eq, dcdn_alpha_ga = equilibrium2D.equilibrium_gamma(grid, c_alpha, T, c0_gamma, Gx_vec = None, Curvature = None)

        c_gamma_ga = (c_gamma_eq - c0_gamma)*H_ga + c0_gamma
        c_gamma_t  = (c_gamma_eq - c0_gamma)*H_t + c0_gamma
        c_gamma  = (c_gamma_eq - c0_gamma)*H_gat + c0_gamma
        c_theta = np.ones(grid.precipitate_mesh_theta.n) * param.c0_theta()

    if t == 0:
        ### Mass and area calculation
        m_C_gamma = 0
        for i in grid.diffusion_mesh_gamma.points_index:
            m_C_gamma = m_C_gamma + np.sum(np.abs(grid.enriched_mesh.area_triangle[np.nonzero(grid.enriched_mesh.point_in_tri[i])[1]])) / 6 * c_gamma[grid.diffusion_mesh_gamma.enriched_to_gamma[i]]
    
        m_C_alpha = 0
        for i in grid.diffusion_mesh_alpha.points_index:
            m_C_alpha = m_C_alpha + np.sum(np.abs(grid.enriched_mesh.area_triangle[np.nonzero(grid.enriched_mesh.point_in_tri[i])[1]])) / 6 * c_alpha[grid.diffusion_mesh_alpha.enriched_to_alpha[i]]
    
        area_gamma = np.sum(np.abs(grid.enriched_mesh.area_triangle[grid.diffusion_mesh_gamma.triangles_index]))*(1e-6)**2
        area_alpha = np.sum(np.abs(grid.enriched_mesh.area_triangle[grid.diffusion_mesh_alpha.triangles_index]))*(1e-6)**2
        area_theta = (xN-x0) * (yN - y0) *(1e-6)**2 - area_gamma - area_alpha
        area = (xN-x0) * (yN - y0) *(1e-6)**2
        m_C_theta = area_theta * param.c0_theta()
        m_C = (m_C_gamma + m_C_alpha)*(1e-6)**2 + m_C_theta
        N_blocks_gamma = area_gamma/(.36e-9)**3
        N_blocks_alpha = area_alpha/(.29e-9)**3
        N_Fe = N_blocks_gamma*4 + N_blocks_alpha*2
        NA = 6.022140857e23 #1/mol
        mol_Fe = N_Fe/NA
        m_Fe = mol_Fe*55.845 #g
        m_Fe_theta = (m_C_theta / 12.01) * 3 * 55.845
        m = m_Fe + m_C + m_Fe_theta
        m_C0 = m_C
        print('total carbon mass at t = %.4e' %t, '%.4e' %m_C)
        print('total carbon mass percentage at t = %.4e' %t, '%.6e' %(m_C/m_C0))
        print('area gamma/alpha/theta % %.4e' %(area_gamma/area*100), '%.4e' %(area_alpha/area*100), '%.4e' %(area_theta/area*100))
        print('wt % C', m_C/m * 100)

#   ### Plot interfaces before and after
#    plt.close('all')
#
#    f, (ax22, ax32) = plt.subplots(1,2, sharex = True, sharey = True)
#    [ax22.plot(grid.enriched_mesh.points[grid.interface_alphatheta.edges_index[i,:],0],grid.enriched_mesh.points[grid.interface_alphatheta.edges_index[i,:],1], linewidth = 2, color = 'red') for i in range(0,len(grid.interface_alphatheta.edges))]
#    [ax22.plot(grid.enriched_mesh.points[grid.interface_gammatheta.edges_index[i,:],0],grid.enriched_mesh.points[grid.interface_gammatheta.edges_index[i,:],1], linewidth = 2, color = 'green') for i in range(0,len(grid.interface_gammatheta.edges))]
#    [ax22.plot(grid.enriched_mesh.points[grid.interface_gammaalpha.edges_index[i,:],0],grid.enriched_mesh.points[grid.interface_gammaalpha.edges_index[i,:],1], linewidth = 2, color = 'blue') for i in range(0,len(grid.interface_gammaalpha.edges))]
#    [ax22.plot(grid.enriched_mesh.points[grid.interface_gammaalpha.ghost_edges_index[i,:],0],grid.enriched_mesh.points[grid.interface_gammaalpha.ghost_edges_index[i,:],1], linewidth = 1, color = 'blue') for i in range(0,len(grid.interface_gammaalpha.ghost_edges_index))]
#
#    ax22.set_title('Interfaces at t =' + str(t)+'s')
#    ax22.set_xlim([grid.x0-.1, grid.xN + .1])
#    ax22.set_ylim([grid.y0-.1, grid.yN + .1])
#
#
#    f, (ax2, ax3) = plt.subplots(1,2, sharex = True, sharey = True)
#    ax2.triplot(grid.enriched_mesh.points[:,0],grid.enriched_mesh.points[:,1], grid.enriched_mesh.triangles_points_index, linewidth = 1, color = 'black', alpha = 1)
#    ax2.triplot(grid.background_mesh.points[:,0],grid.background_mesh.points[:,1], grid.background_mesh.triangles_points_index, linewidth = 1, color = 'black', alpha = .3)
#    ax3.triplot(grid.background_mesh.points[:,0],grid.background_mesh.points[:,1], grid.background_mesh.triangles_points_index, linewidth = 1, color = 'black', alpha = .3)
#    [ax2.plot(grid.enriched_mesh.points[grid.interface_alphatheta.edges_index[i,:],0],grid.enriched_mesh.points[grid.interface_alphatheta.edges_index[i,:],1], linewidth = 2, color = 'red') for i in range(0,len(grid.interface_alphatheta.edges))]
#    [ax2.plot(grid.enriched_mesh.points[grid.interface_gammatheta.edges_index[i,:],0],grid.enriched_mesh.points[grid.interface_gammatheta.edges_index[i,:],1], linewidth = 2, color = 'green') for i in range(0,len(grid.interface_gammatheta.edges))]
#    [ax2.plot(grid.enriched_mesh.points[grid.interface_gammaalpha.edges_index[i,:],0],grid.enriched_mesh.points[grid.interface_gammaalpha.edges_index[i,:],1], linewidth = 2, color = 'blue') for i in range(0,len(grid.interface_gammaalpha.edges))]
#    [ax2.plot(grid.enriched_mesh.points[grid.interface_gammaalpha.ghost_edges_index[i,:],0],grid.enriched_mesh.points[grid.interface_gammaalpha.ghost_edges_index[i,:],1], linewidth = 2, color = 'lightblue') for i in range(0,len(grid.interface_gammaalpha.ghost_edges_index))]
#
#    if grid.background_mesh.i_ghost_fixed != -1:
#        ax2.plot(grid.background_mesh.points[grid.background_mesh.i_ghost_fixed,0], grid.background_mesh.points[grid.background_mesh.i_ghost_fixed,1], 'og')

    ### Calculating dc/dn of alpha on ga-interface
    dcdxy_alpha_ga = np.zeros((grid.interface_gammaalpha.n,2))
    k=0
    ga_points = grid.interface_gammaalpha.points_index
    for j in ga_points:
        xp = grid.enriched_mesh.points[j]
        other_ii_all = np.unique(grid.enriched_mesh.triangles_points_index[grid.enriched_mesh.point_in_tri[j].nonzero()[1]].flatten())
        #alpha
        other_ii = np.intersect1d(other_ii_all[np.where(grid.level_set_gammaalpha_enriched.phi[other_ii_all] >= 0) ], other_ii_all[np.where(grid.level_set_alphatheta_enriched.phi[other_ii_all] <= 0) ])
        other_i = grid.enriched_mesh.points[other_ii]
        hp = np.max(la.norm(other_i - xp, axis = 1))
        n_row = len(other_ii)
        interface_points = grid.has_point_been_shifted[other_ii,1] == -1
        BP_points_bool = grid.is_point_BP[other_ii,0] == 0
        has_internal = any(np.logical_and(interface_points, BP_points_bool))
        new_other_ii = other_ii.copy()
        r = 2
        while not(has_internal):
            for l in range(0,len(new_other_ii)):
                for tri_i in grid.enriched_mesh.point_in_tri[new_other_ii[l]].nonzero()[1]:
                    trii = grid.enriched_mesh.triangles_points_index[tri_i]
                    trii_alpha_coord = grid.diffusion_mesh_alpha.enriched_to_alpha[trii]
                    if all(trii_alpha_coord != -1):
                        for i in np.setdiff1d(trii, new_other_ii[l]):
                            if la.norm(grid.enriched_mesh.points[i] - xp, 2) <= r * hp:
                                other_ii = np.append(other_ii, i)
            new_other_ii = np.unique(other_ii)
            interface_points = grid.has_point_been_shifted[new_other_ii,1] == -1
            BP_points_bool = grid.is_point_BP[new_other_ii,0] == 0
            has_internal = any(np.logical_and(interface_points, BP_points_bool))
            r = r+1
        where_internal = np.where(np.logical_and(interface_points, BP_points_bool) == True)
        internal_patches = grid.enriched_mesh.patch_matrix[new_other_ii[where_internal]].nonzero()[1]
        internal_patches_alpha = np.intersect1d(internal_patches[np.where(grid.level_set_gammaalpha_enriched.phi[internal_patches] >= 0) ], internal_patches[np.where(grid.level_set_alphatheta_enriched.phi[internal_patches] <= 0) ])
        other_ii = np.union1d(other_ii, internal_patches_alpha)
        n_row = len(other_ii)
        other_i = grid.enriched_mesh.points[other_ii]
        eta = (other_i - xp)/hp
        M_eta = np.ones((n_row, 6))
        c_a  = np.zeros(n_row)
        for i in range(0,n_row):
            etaj = eta[i]
            M_eta[i,1:6] = np.array([etaj[0], etaj[1], etaj[0]**2, etaj[0] * etaj[1], etaj[1]**2])
            c_a[i ] = c_alpha[grid.diffusion_mesh_alpha.enriched_to_alpha[other_ii[i]]]
        B = la.lstsq(M_eta, c_a)[0][[1,2]]
        dcdxy_alpha_ga[k] = B/hp
        k=k+1
    if loaded == 0:
        eq_dcdn_alpha_ga = dcdn_alpha_ga
        dcdn_alpha_ga = np.sum(dcdxy_alpha_ga * grid.n_ga[ga_points], axis = 1)
        if any(eq_dcdn_alpha_ga - dcdn_alpha_ga != 0):
            print('dcdn_alpha different than loaded dcdn?')

    t_vec = t
    start_time = time.monotonic()
    for k in range(0, 5):  # Amount of time steps can also use while t < t_end
        t_old = t
        gelukt = 0
        grid_old = deepcopy(grid)
        c_gamma_old = c_gamma.copy()
        c_alpha_old = c_alpha.copy()
        dcdn_alpha_ga_old = dcdn_alpha_ga.copy()
#       ### Do a timestep without testing for errors, extra option (do reinitialization, cooling rate, if there is a end time, if Gibbs-Thomson effect is used) set to default
#        c_gamma, c_alpha, dcdn_alpha_ga, t = do_time_step2D.do_time_step2D(c_gamma, c_alpha, dcdn_alpha_ga, param.c0_theta(), grid, t, t_old, t_end, k, ext, Reinit = False, Cooling_rate = 0, Curvature = False)   
        try:
            c_gamma, c_alpha, dcdn_alpha_ga, t = do_time_step2D.do_time_step2D(c_gamma, c_alpha, dcdn_alpha_ga, param.c0_theta(), grid, t, t_old, t_end, k, ext, Reinit = False, Cooling_rate = 0, Curvature = False)                                                                              
            gelukt = 1
        except Exception as e:
            type_, value_, traceback_ = sys.exc_info()
            print('k ',k,'Error: ',e)
            traceback.format_tb(traceback_)
            data = {'grid': grid,
                    'c_gamma': c_gamma,
                    'c_alpha': c_alpha,
                    'dcdn_alpha_ga': dcdn_alpha_ga,
                    't': t,
                    't_old': t_old,
                    'm_C_gamma' : m_C_gamma,
                    'm_C_alpha' : m_C_alpha,
                    'm_C_theta' : m_C_theta,
                    'area_gamma' : area_gamma,
                    'area_alpha' : area_alpha,
                    'area_theta' : area_theta,
                    't_vec' : t_vec}
            output = open('./Results/'+name+'_'+str(k)+'dtk1_error.pkl', 'wb')
            pickle.dump(data, output)
            output.close()
            break
        ### Check if area is conserved, else probably error by overlapping triangles
        if np.abs(np.sum(np.abs(grid.enriched_mesh.area_triangle)) - (grid.xN - grid.x0) * (grid.yN-grid.y0)) > 1e-10:
            data = {'grid': grid,
                    'c_gamma': c_gamma,
                    'c_alpha': c_alpha,
                    'dcdn_alpha_ga': dcdn_alpha_ga,
                    't': t,
                    't_old': t_old,
                    'm_C_gamma' : m_C_gamma,
                    'm_C_alpha' : m_C_alpha,
                    'm_C_theta' : m_C_theta,
                    'area_gamma' : area_gamma,
                    'area_alpha' : area_alpha,
                    'area_theta' : area_theta,
                    't_vec' : t_vec}
            output = open('./Results/'+name+'_'+str(k)+'dtk1_area_triangles.pkl', 'wb')
            pickle.dump(data, output)
            output.close()
            print('area triangles not sum to 625')
            break
        ### Check if time step is not too small (velocity too high)
        if t - t_old < 1e-12:
            print('dt is really small, vn_t ext bad')
            data = {'grid': grid,
                    'c_gamma': c_gamma,
                    'c_alpha': c_alpha,
                    'dcdn_alpha_ga': dcdn_alpha_ga,
                    't': t,
                    't_old': t_old,
                    'm_C_gamma' : m_C_gamma,
                    'm_C_alpha' : m_C_alpha,
                    'm_C_theta' : m_C_theta,
                    'area_gamma' : area_gamma,
                    'area_alpha' : area_alpha,
                    'area_theta' : area_theta,
                    't_vec' : t_vec}
            output = open('./Results/'+name+'_'+str(k)+'dt_dt_small.pkl', 'wb')
            pickle.dump(data, output)
            output.close()
            gelukt = 0
            break


        print('k = ',k,'  n_gamma ', grid.diffusion_mesh_gamma.n, 'n_alpha ', grid.diffusion_mesh_alpha.n)

        ### Check maximum principle on ga-interface
        ga_points_no_x = grid.diffusion_mesh_gamma.enriched_to_gamma[np.setdiff1d(grid.interface_gammaalpha.points_index, grid.interface_gammaalpha.triple_points_index)]
        if np.min(c_gamma[ga_points_no_x]) < 0.9999*param.c_gammaalpha(T):
            print('a c_g point on ga < c_ga. min ', np.min(c_gamma[ga_points_no_x]), 'c_ga ', param.c_gammaalpha(T))
            
        ### Calculate masses, areas
        m_C_gamma_new = 0
        for i in grid.diffusion_mesh_gamma.points_index:
            m_C_gamma_new = m_C_gamma_new+ np.sum(np.abs(grid.enriched_mesh.area_triangle[np.nonzero(grid.enriched_mesh.point_in_tri[i])[1]])) / 6 * c_gamma[grid.diffusion_mesh_gamma.enriched_to_gamma[i]]
        m_C_gamma  = np.append(m_C_gamma , m_C_gamma_new)
    
        m_C_alpha_new = 0
        for i in grid.diffusion_mesh_alpha.points_index:
            m_C_alpha_new = m_C_alpha_new + np.sum(np.abs(grid.enriched_mesh.area_triangle[np.nonzero(grid.enriched_mesh.point_in_tri[i])[1]])) / 6 * c_alpha[grid.diffusion_mesh_alpha.enriched_to_alpha[i]]
        m_C_alpha = np.append(m_C_alpha, m_C_alpha_new)
            
        area_gamma_new = np.sum(np.abs(grid.enriched_mesh.area_triangle[grid.diffusion_mesh_gamma.triangles_index]))*(1e-6)**2
        area_alpha_new = np.sum(np.abs(grid.enriched_mesh.area_triangle[grid.diffusion_mesh_alpha.triangles_index]))*(1e-6)**2
        area_theta_new = (xN-x0) * (yN - y0) *(1e-6)**2 - area_gamma_new - area_alpha_new
        area_gamma = np.append(area_gamma, area_gamma_new)
        area_alpha = np.append(area_alpha, area_alpha_new)
        area_theta = np.append(area_theta, area_theta_new)
#        area = (xN-x0) * (yN - y0) *(1e-6)**2
        m_C_theta_new = area_theta_new * param.c0_theta()
        m_C_theta = np.append(m_C_theta, m_C_theta_new)
        m_C = (m_C_gamma_new + m_C_alpha_new)*(1e-6)**2 + m_C_theta_new
        N_blocks_gamma = area_gamma_new/(.36e-9)**3
        N_blocks_alpha = area_alpha_new/(.29e-9)**3
        N_Fe = N_blocks_gamma*4 + N_blocks_alpha*2
        NA = 6.022140857e23 #1/mol
        mol_Fe = N_Fe/NA
        m_Fe = mol_Fe*55.845 #g
        m_Fe_theta = (m_C_theta_new / 12.01) * 3 * 55.845
        m = m_Fe + m_C + m_C_theta_new
        print('total carbon mass at t = %.4e' %t, '%.4e' %m_C)
        print('total carbon mass percentage at t = %.4e' %t, '%.6e' %(m_C/m_C0))
        print('area gamma/alpha/theta % %.4e' %(area_gamma_new/area*100), '%.4e' %(area_alpha_new/area*100), '%.4e' %(area_theta_new/area*100))
        t_vec =  np.append(t_vec, t)


    end_time = time.monotonic()
    print('Time taken:',timedelta(seconds=end_time - start_time),'\n')

#    ### After plot
#    ax3.triplot(grid.enriched_mesh.points[:,0],grid.enriched_mesh.points[:,1], grid.enriched_mesh.triangles_points_index, linewidth = 1, color = 'black', alpha = 1)
#    [ax3.plot(grid.enriched_mesh.points[grid.interface_alphatheta.edges_index[i,:],0],grid.enriched_mesh.points[grid.interface_alphatheta.edges_index[i,:],1], linewidth = 2, color = 'red') for i in range(0,len(grid.interface_alphatheta.edges))]
#    [ax3.plot(grid.enriched_mesh.points[grid.interface_gammatheta.edges_index[i,:],0],grid.enriched_mesh.points[grid.interface_gammatheta.edges_index[i,:],1], linewidth = 2, color = 'green') for i in range(0,len(grid.interface_gammatheta.edges))]
#    [ax3.plot(grid.enriched_mesh.points[grid.interface_gammaalpha.edges_index[i,:],0],grid.enriched_mesh.points[grid.interface_gammaalpha.edges_index[i,:],1], linewidth = 2, color = 'blue', alpha = 1) for i in range(0,len(grid.interface_gammaalpha.edges))]
#    [ax3.plot(grid.enriched_mesh.points[grid.interface_gammaalpha.ghost_edges_index[i,:],0],grid.enriched_mesh.points[grid.interface_gammaalpha.ghost_edges_index[i,:],1], linewidth = 2, color = 'lightblue') for i in range(0,len(grid.interface_gammaalpha.ghost_edges_index))]
#
#    [ax32.plot(grid.enriched_mesh.points[grid.interface_alphatheta.edges_index[i,:],0],grid.enriched_mesh.points[grid.interface_alphatheta.edges_index[i,:],1], linewidth = 2, color = 'red') for i in range(0,len(grid.interface_alphatheta.edges))]
#    [ax32.plot(grid.enriched_mesh.points[grid.interface_gammatheta.edges_index[i,:],0],grid.enriched_mesh.points[grid.interface_gammatheta.edges_index[i,:],1], linewidth = 2, color = 'green') for i in range(0,len(grid.interface_gammatheta.edges))]
#    [ax32.plot(grid.enriched_mesh.points[grid.interface_gammaalpha.edges_index[i,:],0],grid.enriched_mesh.points[grid.interface_gammaalpha.edges_index[i,:],1], linewidth = 2, color = 'blue') for i in range(0,len(grid.interface_gammaalpha.edges))]
#    [ax32.plot(grid.enriched_mesh.points[grid.interface_gammaalpha.ghost_edges_index[i,:],0],grid.enriched_mesh.points[grid.interface_gammaalpha.ghost_edges_index[i,:],1], linewidth = 1, color = 'blue') for i in range(0,len(grid.interface_gammaalpha.ghost_edges_index))]
#    ax32.set_title('Interfaces at t =' + str(t) + 's')
#    ax32.set_xlim([grid.x0-.1, grid.xN + .1])
#    ax32.set_ylim([grid.y0-.1, grid.yN + .1])


