###
# equilibrium2D.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/>.
###

import numpy as np
import scipy.linalg as la
from scipy.sparse import csr_matrix
from scipy.sparse import coo_matrix
from scipy.sparse import lil_matrix
from scipy.sparse import eye
import scipy.sparse.linalg as spla
from scipy.sparse import vstack, hstack

import parameters as param


def setcsrrow2id(amat, rowind):
    indptr = amat.indptr
    values = amat.data
    indxs = amat.indices

    # get the range of the data that is changed
    rowpa = indptr[rowind]
    rowpb = indptr[rowind+1]

    # new value and its new rowindex
    values[rowpa] = 1.0
    indxs[rowpa] = rowind

    # number of new zero values
    diffvals = rowpb - rowpa - 1

    # filter the data and indices and adjust the range
    values = np.r_[values[:rowpa+1], values[rowpb:]]
    indxs = np.r_[indxs[:rowpa+1], indxs[rowpb:]]
    indptr = np.r_[indptr[:rowind+1], indptr[rowind+1:]-diffvals]

    # hard set the new sparse data
    amat.indptr = indptr
    amat.data = values
    amat.indices = indxs



def equilibrium_alpha(grid, T, Curvature = False):
    c0_alpha = (param.c_alphatheta(T) + param.c_alphagamma(T))/2

    BP_points = np.where(grid.is_point_BP[grid.diffusion_mesh_alpha.points_index,0] == 1)[0]
    ag_points = grid.diffusion_mesh_alpha.enriched_to_alpha[grid.interface_gammaalpha.points_index]
    at_points = grid.diffusion_mesh_alpha.enriched_to_alpha[grid.interface_alphatheta.points_index]
    x_points = grid.diffusion_mesh_alpha.enriched_to_alpha[grid.triple_points_index]

    I = np.zeros(9*grid.diffusion_mesh_alpha.n_tri, dtype = np.int)
    J = np.zeros(9*grid.diffusion_mesh_alpha.n_tri, dtype = np.int)
    A_V = np.zeros(9*grid.diffusion_mesh_alpha.n_tri)
    l = 1
    for i in grid.diffusion_mesh_alpha.triangles_index:
        trii = grid.enriched_mesh.triangles_points_index[i]
        tri = grid.enriched_mesh.triangles[i]
        trii_alpha_coord = grid.diffusion_mesh_alpha.enriched_to_alpha[trii]
        Delta = 2*grid.enriched_mesh.area_triangle[i]

        beta0 = (tri[1,1] - tri[2,1])/Delta
        beta1 = (tri[2,1] - tri[0,1])/Delta
        beta2 = (tri[0,1] - tri[1,1])/Delta
        gamma0 = (tri[2,0] - tri[1,0])/Delta
        gamma1 = (tri[0,0] - tri[2,0])/Delta
        gamma2 = (tri[1,0] - tri[0,0])/Delta

        I[9*(l-1):9*(l-1)+9] = np.hstack((np.ones(3, dtype = np.int) * trii_alpha_coord[0], np.ones(3, dtype = np.int) * trii_alpha_coord[1], np.ones(3, dtype = np.int) * trii_alpha_coord[2]))
        J[9*(l-1):9*(l-1)+9] = np.hstack((trii_alpha_coord, trii_alpha_coord, trii_alpha_coord))


        A_V[9*(l-1):9*(l-1)+9] = - param.D_alpha(T) * (np.abs(Delta) / 2 * np.hstack((np.array([beta0**2 + gamma0**2           , beta0 * beta1 + gamma0 * gamma1, beta0 * beta2 + gamma0 * gamma2]) ,
                                                                                      np.array([beta1 * beta0 + gamma1 * gamma0, beta1**2 + gamma1**2           , beta1 * beta2 + gamma1 * gamma2]) ,
                                                                                      np.array([beta2 * beta0 + gamma2 * gamma0, beta2 * beta1 + gamma2 * gamma1, beta2**2 + gamma2**2           ]) )))
        l = l+1

    A_alpha = coo_matrix((A_V,(I,J)), shape = (grid.diffusion_mesh_alpha.n,grid.diffusion_mesh_alpha.n))
    A_alpha = A_alpha.tocsr()
    f = np.zeros(grid.diffusion_mesh_alpha.n)

#    #Dirichlet everywhere
#    A_alpha[BP_points] = eye(grid.diffusion_mesh_alpha.n).tolil()[BP_points]
#    A_alpha[ag_points] = eye(grid.diffusion_mesh_alpha.n).tolil()[ag_points]
#    A_alpha[at_points] = eye(grid.diffusion_mesh_alpha.n).tolil()[at_points]
#    A_alpha = A_alpha.tocsr()
#
#    f[BP_points] = c0_alpha
#    f[ag_points] = param.c_alphagamma(T)
#    f[at_points] = param.c_alphatheta(T)
#    f[x_points ] = (param.c_alphagamma(T) + param.c_alphatheta(T))/2
#
#    qr = spla.spsolve(A_alpha,f)

    [setcsrrow2id(A_alpha, BP_points[i]) for i in range(0, len(BP_points))]
    [setcsrrow2id(A_alpha, ag_points[i]) for i in range(0, len(ag_points))]

    f[BP_points] = c0_alpha
    if not(Curvature):
        f[ag_points] = param.c_alphagamma(T)
    else:
        for i in range(0, grid.interface_gammaalpha.n):
            f[ag_points[i]] = param.c_alphagamma_GT(T, grid.enriched_mesh.curvature_ga[i])

    I = np.zeros(4*(grid.interface_alphatheta.n_edges), dtype = np.int)
    J = np.zeros(4*(grid.interface_alphatheta.n_edges), dtype = np.int)
    S_BI_V = np.zeros(4*(grid.interface_alphatheta.n_edges))
    f_BI = np.zeros(grid.diffusion_mesh_alpha.n)
    l = 1
    for el in grid.interface_alphatheta.edges_index:
        L = la.norm(grid.enriched_mesh.points[el[0]] - grid.enriched_mesh.points[el[1]])
        I[4*(l-1):4*(l-1)+4] = grid.diffusion_mesh_alpha.enriched_to_alpha[np.array([el[0], el[0], el[1], el[1]])]
        J[4*(l-1):4*(l-1)+4] = grid.diffusion_mesh_alpha.enriched_to_alpha[np.array([el[0], el[1], el[0], el[1]])]
        S_BI_V[4*(l-1):4*(l-1)+4] = L/6*np.array([2, 1, 1, 2])
        l = l + 1
        f_BI[grid.diffusion_mesh_alpha.enriched_to_alpha[el]] = f_BI[grid.diffusion_mesh_alpha.enriched_to_alpha[el]] + L/2 * np.ones(2)
    S_BI = coo_matrix((S_BI_V,(I,J)), shape = (grid.diffusion_mesh_alpha.n,grid.diffusion_mesh_alpha.n))
    S_BI = S_BI.tocsr()
#    S_BI[x_points] = 0
#    f_BI[x_points] = 0


#    at_points = grid.diffusion_mesh_alpha.enriched_to_alpha[grid.interface_alphatheta.points_index]

#    S = A_alpha.copy()
    r = 1
    eps = 1e6
    FPT = 1e-7
#    qr = (param.c_alphatheta(T)+param.c_alphagamma(T))/2*np.ones(grid.diffusion_mesh_alpha.n)
    qr = param.c_alphatheta(T)*np.ones(grid.diffusion_mesh_alpha.n)
#    qr = param.c_alphagamma(T)*np.ones(grid.diffusion_mesh_alpha.n)
    while eps >= FPT:
        S_BI_cp_V = np.zeros(4*(grid.interface_alphatheta.n_edges))
        f_cp_BI = np.zeros(grid.diffusion_mesh_alpha.n)
        l = 1
        for el in grid.interface_alphatheta.edges_index:
            L = la.norm(grid.enriched_mesh.points[el[0]] - grid.enriched_mesh.points[el[1]])
            el_a = grid.diffusion_mesh_alpha.enriched_to_alpha[el]
            ### Not lumped reactoin
#            S_BI_cp_V[4*(l-1):4*(l-1)+4] = L/12*np.array([3*qr[el_a[0]] + qr[el_a[1]], qr[el_a[0]] +   qr[el_a[1]],
#                                                            qr[el_a[0]] + qr[el_a[1]], qr[el_a[0]] + 3*qr[el_a[1]]])
#            f_cp_BI[el_a] = f_cp_BI[el_a] + L/6 * np.array([2*qr[el_a[0]] +   qr[el_a[1]],
#                                                              qr[el_a[0]] + 2*qr[el_a[1]] ])
            #Lumped reaction term RHS
            S_BI_cp_V[4*(l-1):4*(l-1)+4] = L/12*np.array([3*qr[el_a[0]] + qr[el_a[0]], qr[el_a[0]] +   qr[el_a[0]],
                                                            qr[el_a[1]] + qr[el_a[1]], qr[el_a[1]] + 3*qr[el_a[1]]])
            f_cp_BI[el_a] = f_cp_BI[el_a] + L/2 * np.array([qr[el_a[0]],
                                                            qr[el_a[1]] ])
            l = l + 1

            
        S_BI_cp = coo_matrix((S_BI_cp_V,(I,J)), shape = (grid.diffusion_mesh_alpha.n,grid.diffusion_mesh_alpha.n))
        S_BI_cp = S_BI_cp.tocsr()


        #Stiffness matrix BC

        Stilde = - param.K_alphatheta(T) * (param.c0_theta() * S_BI[at_points] - S_BI_cp[at_points]) / param.c0_theta()

#        S[at_points] = A_alpha[at_points] + Stilde
#
#        f[at_points] = - param.K_alphatheta(T) * param.c_alphatheta(T) * (param.c0_theta() * f_BI[at_points] - f_cp_BI[at_points]) / param.c0_theta()


#        D = lil_matrix(A_alpha.shape)
#        Stilde_matrix[at_points] = Stilde
#        for i in at_points:
#            where_Stilde_neg = np.setdiff1d(find(Stilde_matrix[i] < 0)[1], i)
#            for j_neg in where_Stilde_neg:
#                min_neg = min(Stilde_matrix[i,j_neg], Stilde_matrix[j_neg,i])
#                D[i,j_neg] = min_neg
#                D[j_neg,i] = min_neg
#            D[i,i] = - np.sum(D[i])
#        Stilde_matrix[x_points] = 0
#        D[x_points] = 0
#        Stilde_matrix = Stilde_matrix - D
        Stilde_matrix = lil_matrix((grid.diffusion_mesh_alpha.n, grid.diffusion_mesh_alpha.n))
        Stilde_matrix[at_points] = Stilde
        Stilde_matrix = Stilde_matrix.tocsr()
        Stilde_matrix_lumped = lil_matrix((grid.diffusion_mesh_alpha.n, grid.diffusion_mesh_alpha.n))
        Stilde_matrix_lumped.setdiag(csr_matrix.sum(Stilde_matrix, axis = 1))
        Stilde_matrix_lumped = Stilde_matrix_lumped.tocsr()
#            S[at_points] = A_alpha[at_points] + Stilde
        S = A_alpha + Stilde_matrix_lumped
        S[x_points] = eye(grid.diffusion_mesh_alpha.n).tolil()[x_points]

        if not(Curvature):
            f[at_points] = - param.K_alphatheta(T) * param.c_alphatheta(T) * (param.c0_theta() * f_BI[at_points] - f_cp_BI[at_points]) / param.c0_theta()
            f[x_points] = param.c_alphagamma(T)
        else:
            for i in range(0, grid.interface_alphatheta.n):
                f[at_points[i]] = - param.K_alphatheta(T) * param.c_alphatheta_GT(T, grid.enriched_mesh.curvature_at[i]) * (param.c0_theta() * f_BI[at_points[i]] - f_cp_BI[at_points[i]]) / param.c0_theta()
            for i in range(0, len(grid.triple_points_index)):
                f[x_points[i]]  = param.c_alphagamma_GT(T, -grid.enriched_mesh.curvature_ga[grid.interface_gammaalpha.triple_points_index[i]])

        eps = la.norm(S @ qr - f, np.inf)/la.norm(f, np.inf)
#        eps = la.norm(qr - qr_old, np.inf)/la.norm(qr, np.inf)

        if not(S.getformat() == 'csr') and r == 2:
            print('S is not csr, S is ',S.getformat(), 'eq MMP DR')
        qr = spla.spsolve(S,f)

        r = r+1
        if r > 1e3 or eps > 1e7:
            print('inf c_alpha eq while loop, r = ',r, '\n', qr)
            break
    return qr




def equilibrium_gamma(grid, c_alpha, T, c0_gamma, Gx_vec = None, Curvature = False):
    BP_points = np.where(grid.is_point_BP[grid.diffusion_mesh_gamma.points_index,0] == 1)[0]
    gt_points = grid.diffusion_mesh_gamma.enriched_to_gamma[grid.interface_gammatheta.points_index]
    x_points_ga = grid.diffusion_mesh_gamma.enriched_to_gamma[grid.triple_points_index]
    ga_points = grid.interface_gammaalpha.points_index
    if Gx_vec is None:
        #Calc dc_alpha/dn
        dcdxy_alpha_ga = np.zeros((grid.interface_gammaalpha.n,2))
        k=0
        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
                if r == 20:
                    print('big r')
            other_ii = new_other_ii
            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)
            if n_row > 20:
                print('n_row large ',k, xp)

            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
        dcdn_alpha_ga = np.sum(dcdxy_alpha_ga * grid.n_ga[ga_points], axis = 1)
    else:
        dcdxy_alpha_ga = Gx_vec
        dcdn_alpha_ga = np.sum(dcdxy_alpha_ga[grid.diffusion_mesh_alpha.enriched_to_alpha[ga_points]] * grid.n_ga[ga_points], axis = 1)
    ga_points = grid.diffusion_mesh_gamma.enriched_to_gamma[grid.interface_gammaalpha.points_index]
    dcdn_alpha_ga_matrix = lil_matrix((grid.diffusion_mesh_gamma.n, grid.diffusion_mesh_gamma.n))


    I_gt = np.zeros(4*(grid.interface_gammatheta.n_edges), dtype = np.int)
    J_gt = np.zeros(4*(grid.interface_gammatheta.n_edges), dtype = np.int)
    S_BI_V = np.zeros(4*(grid.interface_gammatheta.n_edges))
    f_BI_gt = np.zeros(grid.diffusion_mesh_gamma.n)
    l = 1
    for el in grid.interface_gammatheta.edges_index:
        L = la.norm(grid.enriched_mesh.points[el[0]] - grid.enriched_mesh.points[el[1]])
        I_gt[4*(l-1):4*(l-1)+4] = grid.diffusion_mesh_gamma.enriched_to_gamma[np.array([el[0], el[0], el[1], el[1]])]
        J_gt[4*(l-1):4*(l-1)+4] = grid.diffusion_mesh_gamma.enriched_to_gamma[np.array([el[0], el[1], el[0], el[1]])]
        S_BI_V[4*(l-1):4*(l-1)+4] = L/6*np.array([2, 1, 1, 2])
        l = l + 1
        f_BI_gt[grid.diffusion_mesh_gamma.enriched_to_gamma[el]] = f_BI_gt[grid.diffusion_mesh_gamma.enriched_to_gamma[el]] + L/2 * np.ones(2)
    S_BI_gt = coo_matrix((S_BI_V,(I_gt,J_gt)), shape = (grid.diffusion_mesh_gamma.n,grid.diffusion_mesh_gamma.n))
    S_BI_gt = S_BI_gt.tocsr()

    I_ga = np.zeros(4*(grid.interface_gammaalpha.n_edges), dtype = np.int)
    J_ga = np.zeros(4*(grid.interface_gammaalpha.n_edges), dtype = np.int)
    S_BI_V = np.zeros(4*(grid.interface_gammaalpha.n_edges))
    f_BI_ga = np.zeros(grid.diffusion_mesh_gamma.n)
    l = 1
    for el in grid.interface_gammaalpha.edges_index:
        L = la.norm(grid.enriched_mesh.points[el[0]] - grid.enriched_mesh.points[el[1]], 2)
        I_ga[4*(l-1):4*(l-1)+4] = grid.diffusion_mesh_gamma.enriched_to_gamma[np.array([el[0], el[0], el[1], el[1]])]
        J_ga[4*(l-1):4*(l-1)+4] = grid.diffusion_mesh_gamma.enriched_to_gamma[np.array([el[0], el[1], el[0], el[1]])]
        S_BI_V[4*(l-1):4*(l-1)+4] = L/6*np.array([2, 1, 1, 2])
        l = l + 1
        f_BI_ga[grid.diffusion_mesh_gamma.enriched_to_gamma[el]] = f_BI_ga[grid.diffusion_mesh_gamma.enriched_to_gamma[el]] + L/2 * np.ones(2)
    S_BI_ga = coo_matrix((S_BI_V,(I_ga,J_ga)), shape = (grid.diffusion_mesh_gamma.n,grid.diffusion_mesh_gamma.n))
    S_BI_ga = S_BI_ga.tocsr()


    I = np.zeros(9*grid.diffusion_mesh_gamma.n_tri, dtype = np.int)
    J = np.zeros(9*grid.diffusion_mesh_gamma.n_tri, dtype = np.int)
    A_V = np.zeros(9*grid.diffusion_mesh_gamma.n_tri)
    l = 1
    for i in grid.diffusion_mesh_gamma.triangles_index:
        trii = grid.enriched_mesh.triangles_points_index[i]
        tri = grid.enriched_mesh.triangles[i]
        trii_gamma_coord = grid.diffusion_mesh_gamma.enriched_to_gamma[trii]
        Delta = 2*grid.enriched_mesh.area_triangle[i]

        beta0 = (tri[1,1] - tri[2,1])/Delta
        beta1 = (tri[2,1] - tri[0,1])/Delta
        beta2 = (tri[0,1] - tri[1,1])/Delta
        gamma0 = (tri[2,0] - tri[1,0])/Delta
        gamma1 = (tri[0,0] - tri[2,0])/Delta
        gamma2 = (tri[1,0] - tri[0,0])/Delta

        I[9*(l-1):9*(l-1)+9] = np.hstack((np.ones(3, dtype = np.int) * trii_gamma_coord[0], np.ones(3, dtype = np.int) * trii_gamma_coord[1], np.ones(3, dtype = np.int) * trii_gamma_coord[2]))
        J[9*(l-1):9*(l-1)+9] = np.hstack((trii_gamma_coord, trii_gamma_coord, trii_gamma_coord))



        A_V[9*(l-1):9*(l-1)+9] = - param.D_gamma(T) * np.abs(Delta) / 2 * np.array([(beta0**2 + gamma0**2)          , beta0 * beta1 + gamma0 * gamma1, beta0 * beta2 + gamma0 * gamma2,
                                                                                     beta1 * beta0 + gamma1 * gamma0, beta1**2 + gamma1**2           , beta1 * beta2 + gamma1 * gamma2,
                                                                                     beta2 * beta0 + gamma2 * gamma0, beta2 * beta1 + gamma2 * gamma1, beta2**2 + gamma2**2           ])
        l = l+1

    A_gamma = coo_matrix((A_V,(I,J)), shape = (grid.diffusion_mesh_gamma.n,grid.diffusion_mesh_gamma.n))
    A_gamma = A_gamma.tocsr()

##   #Dirichlet
#    A_gamma[BP_points] = eye(grid.diffusion_mesh_gamma.n).tolil()[BP_points]
#    A_gamma[ga_points] = eye(grid.diffusion_mesh_gamma.n).tolil()[ga_points]
#    A_gamma[gt_points] = eye(grid.diffusion_mesh_gamma.n).tolil()[gt_points]
#    A_gamma = A_gamma.tocsr()
#
#    f[BP_points] = c0_gamma
#    f[ga_points] = param.c_gammaalpha(T)
#    f[gt_points] = param.c_gammatheta(T)
#    f[x_points_ga ] = (param.c_gammaalpha(T) + param.c_gammatheta(T))/2
#
#    qr = spla.spsolve(A_gamma,f)

    [setcsrrow2id(A_gamma, BP_points[i]) for i in range(0, len(BP_points))]
#    [setcsrrow2id(A_gamma, x_points_ga[i]) for i in range(0, len(x_points_ga))]
    f = np.zeros(grid.diffusion_mesh_gamma.n)
    f[BP_points] = c0_gamma
#    f[x_points_ga] = param.c_gammatheta(T) - param.K_alphatheta(T) / param.K_gammatheta(T) * (param.c_alphatheta(T) - param.c_alphagamma(T))

    dcdn_alpha_ga_matrix[ga_points, ga_points] = dcdn_alpha_ga * f_BI_ga[ga_points]
#    print('fixed0')
    dcdn_alpha_ga_matrix = dcdn_alpha_ga_matrix.tocsr()

#    S = A_gamma.copy()
    r = 1
    eps = 1e6
    FPT = 1e-7
    qr = (param.c_gammatheta(T)+param.c_gammaalpha(T)+c0_gamma)/3*np.ones(grid.diffusion_mesh_gamma.n)
    qr = param.c_gammaalpha(T)*np.ones(grid.diffusion_mesh_gamma.n)
    while eps >= FPT:
        S_BI_cp_V = np.zeros(4*(grid.interface_gammatheta.n_edges))
        f_cp_BI_gt = np.zeros(grid.diffusion_mesh_gamma.n)
        l = 1
        for el in grid.interface_gammatheta.edges_index:
            L = la.norm(grid.enriched_mesh.points[el[0]] - grid.enriched_mesh.points[el[1]], 2)
            el_g = grid.diffusion_mesh_gamma.enriched_to_gamma[el]
            ### Not lumped reaction term
#            S_BI_cp_V[4*(l-1):4*(l-1)+4] = L/12*np.array([3*qr[el_g[0]] + qr[el_g[1]], qr[el_g[0]] +   qr[el_g[1]],
#                                                            qr[el_g[0]] + qr[el_g[1]], qr[el_g[0]] + 3*qr[el_g[1]]])
#            f_cp_BI_gt[el_g] = f_cp_BI_gt[el_g] + L/6 * np.array([2*qr[el_g[0]] +   qr[el_g[1]],
#                                                                    qr[el_g[0]] + 2*qr[el_g[1]] ])
            
            #Lumped reaction term 
            S_BI_cp_V[4*(l-1):4*(l-1)+4] = L/12*np.array([3*qr[el_g[0]] + qr[el_g[0]], qr[el_g[0]] +   qr[el_g[0]],
                                                            qr[el_g[1]] + qr[el_g[1]], qr[el_g[1]] + 3*qr[el_g[1]]])
            f_cp_BI_gt[el_g] = f_cp_BI_gt[el_g] + L/2 * np.array([qr[el_g[0]],
                                                                  qr[el_g[1]] ])
    
            l = l + 1
            
        S_BI_cp_gt = coo_matrix((S_BI_cp_V,(I_gt,J_gt)), shape = (grid.diffusion_mesh_gamma.n,grid.diffusion_mesh_gamma.n))
        S_BI_cp_gt = S_BI_cp_gt.tocsr()

        S_BI_cp_V = np.zeros(4*(grid.interface_gammaalpha.n_edges))
        f_cp_BI_ga = np.zeros(grid.diffusion_mesh_gamma.n)
        l = 1
        for el in grid.interface_gammaalpha.edges_index:
            L = la.norm(grid.enriched_mesh.points[el[0]] - grid.enriched_mesh.points[el[1]], 2)
            el_g = grid.diffusion_mesh_gamma.enriched_to_gamma[el]
            #Not lumped reaction term
#            S_BI_cp_V[4*(l-1):4*(l-1)+4] = L/12*np.array([3*qr[el_g[0]] + qr[el_g[1]], qr[el_g[0]] +   qr[el_g[1]],
#                                                            qr[el_g[0]] + qr[el_g[1]], qr[el_g[0]] + 3*qr[el_g[1]]])
#            f_cp_BI_ga[el_g] = f_cp_BI_ga[el_g] + L/6 * np.array([2*qr[el_g[0]] +   qr[el_g[1]],
#                                                                    qr[el_g[0]] + 2*qr[el_g[1]] ])
            #Lumped reaction term
            f_cp_BI_ga[el_g] = f_cp_BI_ga[el_g] + L/2 * np.array([qr[el_g[0]],
                                                                  qr[el_g[1]] ])
            S_BI_cp_V[4*(l-1):4*(l-1)+4] = L/12*np.array([3*qr[el_g[0]] + qr[el_g[0]], qr[el_g[0]] +   qr[el_g[0]],
                                                            qr[el_g[1]] + qr[el_g[1]], qr[el_g[1]] + 3*qr[el_g[1]]])
            l = l + 1
        S_BI_cp_ga = coo_matrix((S_BI_cp_V,(I_ga,J_ga)), shape = (grid.diffusion_mesh_gamma.n,grid.diffusion_mesh_gamma.n))
        S_BI_cp_ga = S_BI_cp_ga.tocsr()

    	  #Stiffness matrix BC
        Stilde_gt     =                                                     - param.K_gammatheta(T) * (param.c0_theta()      * S_BI_gt[gt_points] - S_BI_cp_gt[gt_points]) / param.c0_theta()
        if not(Curvature):
            Stilde_ga = (param.D_alpha(T) * dcdn_alpha_ga_matrix[ga_points] - param.K_gammaalpha(T) * (param.c_alphagamma(T) * S_BI_ga[ga_points] - S_BI_cp_ga[ga_points]))/ param.c_alphagamma(T)
        else:
            Stilde_ga = lil_matrix((grid.interface_gammaalpha.n, grid.diffusion_mesh_gamma.n))
            for i in range(0, grid.interface_gammaalpha.n):
                Stilde_ga[i] = (param.D_alpha(T) * dcdn_alpha_ga_matrix[ga_points[i]] - param.K_gammaalpha(T) * (param.c_alphagamma_GT(T, grid.enriched_mesh.curvature_ga[i]) * S_BI_ga[ga_points[i]] - S_BI_cp_ga[ga_points[i]]))/ param.c_alphagamma_GT(T, grid.enriched_mesh.curvature_ga[i])



        Stilde_matrix = lil_matrix((grid.diffusion_mesh_gamma.n, grid.diffusion_mesh_gamma.n))
        Stilde_matrix[ga_points] = Stilde_ga
        Stilde_matrix[gt_points] = Stilde_gt
        Stilde_matrix[x_points_ga] = Stilde_ga[grid.interface_gammaalpha.triple_points_index] + Stilde_gt[grid.interface_gammatheta.triple_points_index]

        Stilde_matrix = Stilde_matrix.tocsr()
        Stilde_matrix_lumped = lil_matrix((grid.diffusion_mesh_gamma.n, grid.diffusion_mesh_gamma.n))
        Stilde_matrix_lumped.setdiag(csr_matrix.sum(Stilde_matrix, axis = 1))
        Stilde_matrix_lumped = Stilde_matrix_lumped.tocsr()
        S = A_gamma + Stilde_matrix_lumped
        [setcsrrow2id(S, x_points_ga[i]) for i in range(0, len(x_points_ga))]

#        f[ga_points] =  - param.K_gammaalpha(T) * param.c_gammaalpha(T) * (param.c_alphagamma(T) - qI_ga) / param.c_alphagamma(T)
#        f[gt_points] =  - param.K_gammatheta(T) * param.c_gammatheta(T) * (param.c0_theta()      - qI_gt) / param.c0_theta()
        if not(Curvature):
            f[ga_points] = - param.K_gammaalpha(T) * param.c_gammaalpha(T) * (param.c_alphagamma(T) * f_BI_ga[ga_points] - f_cp_BI_ga[ga_points]) / param.c_alphagamma(T)
            f[gt_points] = - param.K_gammatheta(T) * param.c_gammatheta(T) * (param.c0_theta()      * f_BI_gt[gt_points] - f_cp_BI_gt[gt_points]) / param.c0_theta()
            f[x_points_ga] = param.c_gammatheta(T) - param.K_alphatheta(T) / param.K_gammatheta(T) * (param.c_alphatheta(T) - param.c_alphagamma(T))
        else:
            for i in range(0, grid.interface_gammaalpha.n):
                f[ga_points[i]] = - param.K_gammaalpha(T) * param.c_gammaalpha_GT(T, grid.enriched_mesh.curvature_ga[i]) * (param.c_alphagamma_GT(T, grid.enriched_mesh.curvature_ga[i]) * f_BI_ga[ga_points[i]] - f_cp_BI_ga[ga_points[i]]) / param.c_alphagamma_GT(T, grid.enriched_mesh.curvature_ga[i])
            for i in range(0, grid.interface_gammatheta.n):
                f[gt_points[i]] = - param.K_gammatheta(T) * param.c_gammatheta_GT(T, grid.enriched_mesh.curvature_gt[i]) * (param.c0_theta()      * f_BI_gt[gt_points[i]] - f_cp_BI_gt[gt_points[i]]) / param.c0_theta()
            for i in range(0, len(grid.triple_points_index)):
                x_i_ga = grid.interface_gammaalpha.triple_points_index[i]
                x_i_gt = grid.interface_gammatheta.triple_points_index[i]
                x_i_at = grid.interface_alphatheta.triple_points_index[i]
                f[x_points_ga[i]] = param.c_gammatheta_GT(T, grid.enriched_mesh.curvature_gt[x_i_gt]) - param.K_alphatheta(T) / param.K_gammatheta(T) * (param.c_alphatheta_GT(T, grid.enriched_mesh.curvature_at[x_i_at]) - param.c_alphagamma_GT(T, grid.enriched_mesh.curvature_ga[x_i_ga]))

        eps = la.norm(S @ qr - f, np.inf)/la.norm(f, np.inf)

        if not(S.getformat() == 'csr') and r == 2:
            print('S is not csr, S is ',S.getformat(), 'eq MMP DR')
        qr = spla.spsolve(S,f)

        r = r+1
        if r > 1e3 or eps > 1e7:
            print('inf c_gamma eq while loop, r = ',r, '\n', qr)
            break
    return qr, dcdn_alpha_ga
