###
# do_time_step2D.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
from Grid import cGrid2D
import parameters as param

import scipy.linalg as la
from scipy.sparse import lil_matrix
from scipy.sparse import coo_matrix
from scipy.sparse import csr_matrix
from scipy.sparse import diags
from scipy.sparse import vstack, hstack, bmat
from scipy.sparse import find
import scipy.sparse.linalg as spla
import scipy.linalg as la

from copy import deepcopy

import pickle

def interpolate(trii, tri, p, f):
        p0 = tri[0]
        p1 = tri[1]
        p2 = tri[2]
        Delta = ((p1[0] - p0[0])*(p2[1] - p0[1]) - (p1[1] - p0[1])*(p2[0] - p0[0]))

        b0 = (p1[1] - p2[1])/Delta
        b1 = (p2[1] - p0[1])/Delta
        b2 = (p0[1] - p1[1])/Delta

        g0 = (p2[0] - p1[0])/Delta
        g1 = (p0[0] - p2[0])/Delta
        g2 = (p1[0] - p0[0])/Delta

        a0 = 1 - b0*p0[0] - g0*p0[1]
        a1 = 1 - b1*p1[0] - g1*p1[1]
        a2 = 1 - b2*p2[0] - g2*p2[1]

        psi = np.zeros(3)
        psi[0] = a0 + b0*p[0] + g0*p[1]
        psi[1] = a1 + b1*p[0] + g1*p[1]
        psi[2] = a2 + b2*p[0] + g2*p[1]

        if any(psi > 1) or any(psi < 0):
            where_bad1, = np.where(psi > 1)
            where_bad0, = np.where(psi < 0)
            if any(psi[where_bad1]-1 > 1e-12) or any(np.abs(psi[where_bad0]) > 1e-12):
                print('interpolate in do_time_step2D outside of this triangle, ',trii , 'psi ', psi)

        return f[trii] @ psi

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 do_time_step2D(c_gamma, c_alpha, dcdn_alpha_ga, c0_theta, grid, t, t_old, t_end, k, ext, Reinit = False, Cooling_rate = 0, Curvature = False):
    T_old = param.T_func(t, Cooling_rate)
    #Calculate interface velocities
    vn_gammaalpha, vn_gammatheta, vn_alphatheta = calc_v_vn(grid, c_gamma, c_alpha, dcdn_alpha_ga, param.c0_theta(), T_old, Curvature)   #, v_gammaalpha, v_gammatheta, v_alphatheta
    if grid.interface_gammaalpha.n > 0:
        print('max/min vn_ga ', np.max(vn_gammaalpha), np.min(vn_gammaalpha))
    else:
        print('ga dissolved')
    if grid.interface_gammatheta.n > 1:
        print('max/min vn_gt ', np.max(vn_gammatheta), np.min(vn_gammatheta))
    else:
        print('gt dissolved')
    print('max/min vn_at ', np.max(vn_alphatheta), np.min(vn_alphatheta))

    ### Define interface indices. Interface indices for T_E as T_alpha, T_gamma
    ga_points = grid.interface_gammaalpha.points_index
    ga_points_gamma_coord = grid.diffusion_mesh_gamma.enriched_to_gamma[ga_points]
    ga_points_alpha_coord = grid.diffusion_mesh_alpha.enriched_to_alpha[ga_points]
    gt_points = grid.interface_gammatheta.points_index
    gt_points_gamma_coord = grid.diffusion_mesh_gamma.enriched_to_gamma[gt_points]
    at_points = grid.interface_alphatheta.points_index
    at_points_alpha_coord = grid.diffusion_mesh_alpha.enriched_to_alpha[at_points]



    vn_t = np.array([vn_gammatheta, vn_alphatheta])

    ### Extend interface velocity over T_E
    vn_ex_ga, vn_ex_t = Laplace_vn(grid, vn_gammaalpha, vn_t)

    vn_ex_ga = vn_ex_ga[:,None] * grid.background_mesh.n_ga
    vn_ex_t  = vn_ex_t[:,None]  * grid.background_mesh.n_t

    abs_vn = max(np.max(np.abs(vn_ex_ga)), np.max(np.abs(vn_ex_t)))
    ### Caclucate distance of triple points mesh
    min_distance_xi3 = 10*max(grid.xN-grid.x0, grid.yN-grid.y0)
    for ii3 in grid.triple_points_index:
        xi3 = grid.enriched_mesh.points[ii3]
        for jj3 in np.setdiff1d(grid.triple_points_index, ii3):
            xj3 = grid.enriched_mesh.points[jj3]
            if la.norm(xi3 - xj3) < min_distance_xi3:
                min_distance_xi3 = la.norm(xi3 - xj3)



    CFL0 = 1/5

    dt = CFL0 * grid.background_mesh.h_min / abs_vn
    if t+dt > t_end:
        print('t+dt bigger than t_end, set dt = t_end - t. This is not configured with triple dissolve moments. If triples close and this happens, take bigger t_end')
        dt = t_end - t



    copy_ga = grid.level_set_gammaalpha_enriched.phi_background.copy()
    copy_t  = grid.level_set_alphatheta_enriched.phi_background.copy()

    #UPDATE PHI
    grid.level_set_gammaalpha_enriched.update(grid, dt, vn_ex_ga, 1)
    grid.level_set_alphatheta_enriched.update(grid, dt, vn_ex_t , 2)

    copy_ga_1dt = grid.level_set_gammaalpha_enriched.phi_background.copy()
    copy_t_1dt  = grid.level_set_alphatheta_enriched.phi_background.copy()

    new_x3, shiftbendcut_tri_i0 = grid.find_triple(grid.level_set_gammaalpha_enriched.phi_background, grid.level_set_alphatheta_enriched.phi_background)
    print('CFL ',CFL0, 'which tri ', shiftbendcut_tri_i0[:,0])

    ### Caclucate distance of triple points new mesh with 1 dt
    new_min_distance_xi3 = 10*max(grid.xN-grid.x0, grid.yN-grid.y0)
    i = 0
    for xi3 in new_x3:
        for xj3 in new_x3[np.setdiff1d(np.arange(0,len(new_x3)), i)]:
            if la.norm(xi3 - xj3) < new_min_distance_xi3:
                new_min_distance_xi3 = la.norm(xi3 - xj3)
        i = i+1

    patch_bendcut_tri_i = lil_matrix((len(shiftbendcut_tri_i0), grid.background_mesh.n_tri), dtype = np.int)
    i = 0
    for tri_i in shiftbendcut_tri_i0[:,0]:
        trii = grid.background_mesh.triangles_points_index[tri_i]
        patch_tri_i = np.unique(grid.background_mesh.point_in_tri[trii].nonzero()[1].flatten())
        patch_bendcut_tri_i[i, patch_tri_i] = 1
        i = i+1
    sum_patch_bendcut = lil_matrix.sum(patch_bendcut_tri_i, axis = 0).A1

    ### Check if there is overlap of the triangles containing the triple points and the triangles neighbouring those
    if any(sum_patch_bendcut > 1) and np.unique(shiftbendcut_tri_i0[:,0]).shape[0] == shiftbendcut_tri_i0[:,0].shape[0]: # new_min_distance_xi3 < 2 * grid.background_mesh.h_mean:
        CFL = CFL0 * 2
        grid.level_set_gammaalpha_enriched.phi_background = copy_ga
        grid.level_set_alphatheta_enriched.phi_background = copy_t

        dt = CFL * grid.background_mesh.h_min / abs_vn

        #UPDATE PHI
        grid.level_set_gammaalpha_enriched.update(grid, dt, vn_ex_ga, 1)
        grid.level_set_alphatheta_enriched.update(grid, dt, vn_ex_t , 2)

        new_x3, shiftbendcut_tri_i2 = grid.find_triple(grid.level_set_gammaalpha_enriched.phi_background, grid.level_set_alphatheta_enriched.phi_background)
        print('CFL ', CFL, 'which tri ', shiftbendcut_tri_i2[:,0])

        ### Caclucate distance of triple points new mesh with 2 dt
        new_min_distance_xi3 = 10*max(grid.xN-grid.x0, grid.yN-grid.y0)
        i = 0
        for xi3 in new_x3:
            for xj3 in new_x3[np.setdiff1d(np.arange(0,len(new_x3)), i)]:
                if la.norm(xi3 - xj3) < new_min_distance_xi3:
                    new_min_distance_xi3 = la.norm(xi3 - xj3)
            i = i+1
        patch_bendcut_tri_i = lil_matrix((len(shiftbendcut_tri_i2), grid.background_mesh.n_tri), dtype = np.int)
        i = 0
        for tri_i in shiftbendcut_tri_i2[:,0]:
            trii = grid.background_mesh.triangles_points_index[tri_i]
            patch_tri_i = np.unique(grid.background_mesh.point_in_tri[trii].nonzero()[1].flatten())
            patch_bendcut_tri_i[i, patch_tri_i] = 1
            i = i+1
        sum_patch_bendcut = lil_matrix.sum(patch_bendcut_tri_i, axis = 0).A1

        ### Check again if there is overlap of the triangles containing the triple points and the triangles neighbouring those
        ### If still, do 4 dt time step
        if shiftbendcut_tri_i2.shape[0] == shiftbendcut_tri_i0[:,0].shape[0]:
            CFL = CFL * 2
            grid.level_set_gammaalpha_enriched.phi_background = copy_ga
            grid.level_set_alphatheta_enriched.phi_background = copy_t

            dt = CFL * grid.background_mesh.h_min / abs_vn

            #UPDATE PHI
            grid.level_set_gammaalpha_enriched.update(grid, dt, vn_ex_ga, 1)
            grid.level_set_alphatheta_enriched.update(grid, dt, vn_ex_t , 2)

            new_x3, shiftbendcut_tri_i4 = grid.find_triple(grid.level_set_gammaalpha_enriched.phi_background, grid.level_set_alphatheta_enriched.phi_background)
#            ax2.plot(new_x3[:,0], new_x3[:,1], 'or')
#            print('new_x3 CFL .8 \n', new_x3)
            print('CFL ', CFL, 'which tri ', shiftbendcut_tri_i4[:,0])
            patch_bendcut_tri_i4 = lil_matrix((len(shiftbendcut_tri_i4), grid.background_mesh.n_tri), dtype = np.int)
            i = 0
            ### Check again if there is overlap of the triangles containing the triple points and the triangles neighbouring those
            ### If still, go back to 1 dt and see if meshing is working correctly
            for tri_i in shiftbendcut_tri_i4[:,0]:
                trii = grid.background_mesh.triangles_points_index[tri_i]
                patch_tri_i = np.unique(grid.background_mesh.point_in_tri[trii].nonzero()[1].flatten())
                patch_bendcut_tri_i4[i, patch_tri_i] = 1
                i = i+1
            sum_patch_bendcut = lil_matrix.sum(patch_bendcut_tri_i, axis = 0).A1
            if patch_bendcut_tri_i4.shape[0] == shiftbendcut_tri_i0[:,0].shape[0]:
                CFL = CFL0
                grid.level_set_gammaalpha_enriched.phi_background = copy_ga_1dt
                grid.level_set_alphatheta_enriched.phi_background = copy_t_1dt

    else:
        CFL = CFL0
    dt = CFL * grid.background_mesh.h_min / abs_vn
    if t+dt > t_end:
        dt = t_end - t
    print('t ', t,'dt ', dt, 'CFL ', CFL)

    ### REINITIALIZE
    if Reinit:
        grid.update_Tau_all(1, 1)
        ### ga_reinit/at_reinit are booleans saying if reinitialization of phi^ga and/or phi^t is done
        ga_reinit = grid.level_set_gammaalpha_enriched.reinitialize(grid, 1, Force_reinit = False, plot_grad = True, t_reinit = k)
        at_reinit = grid.level_set_alphatheta_enriched.reinitialize(grid, 2)

    ### phi on enriched mesh is reset to background mesh
    grid.level_set_gammaalpha_enriched.phi = grid.level_set_gammaalpha_enriched.phi_background.copy()
    grid.level_set_alphatheta_enriched.phi = grid.level_set_alphatheta_enriched.phi_background.copy()

    ### MAKE ENRICHED MESH
    grid.update(1)

    if Curvature:
        grid.enriched_mesh.curvature_ga = np.zeros(grid.interface_gammaalpha.n)
        ga_points = grid.interface_gammaalpha.points_index
        for i in range(0, grid.interface_gammaalpha.n):
            grid.enriched_mesh.curvature_ga[i] = -grid.level_set_gammaalpha_enriched.curvature_enr(grid, ga_points[i], 1)

        grid.enriched_mesh.curvature_gt = np.zeros(grid.interface_gammatheta.n)
        gt_points = grid.interface_gammatheta.points_index
        for i in range(0, grid.interface_gammatheta.n):
            grid.enriched_mesh.curvature_gt[i] = -grid.level_set_alphatheta_enriched.curvature_enr(grid, gt_points[i], 2)

        grid.enriched_mesh.curvature_at = np.zeros(grid.interface_alphatheta.n)
        at_points = grid.interface_alphatheta.points_index
        for i in range(0, grid.interface_alphatheta.n):
            grid.enriched_mesh.curvature_at[i] = -grid.level_set_alphatheta_enriched.curvature_enr(grid, at_points[i], 2)



    t = t+dt
    T = param.T_func(t, Cooling_rate)

    ga_points = grid.interface_gammaalpha.points_index
    ga_points_gamma_coord = grid.diffusion_mesh_gamma.enriched_to_gamma[ga_points]
    gt_points = grid.interface_gammatheta.points_index
    at_points = grid.interface_alphatheta.points_index

    ### Extending old concentrations to new mesh
    if ext == 'Shift' or ext == 'Interpolate':
        #SHIFT POINTS
        cexn_alpha = c_alpha[grid.I_diff_alpha]
        vn_mesh_alpha = np.zeros((grid.diffusion_mesh_alpha.n,2))
        vn_mesh_alpha = np.array([(grid.diffusion_mesh_alpha.points[i] - grid.diffusion_points_old_alpha[grid.I_diff_alpha[i]])/dt for i in range(0,grid.diffusion_mesh_alpha.n)])
        cexn_gamma = c_gamma[grid.I_diff_gamma]
        vn_mesh_gamma = np.zeros((grid.diffusion_mesh_gamma.n,2))
        vn_mesh_gamma = np.array([(grid.diffusion_mesh_gamma.points[i] - grid.diffusion_points_old_gamma[grid.I_diff_gamma[i]])/dt for i in range(0,grid.diffusion_mesh_gamma.n)])

        for i in range(0, grid.enriched_mesh.n):
            xi = grid.enriched_mesh.points[i]
            i_gamma = grid.diffusion_mesh_gamma.enriched_to_gamma[i]
            i_alpha = grid.diffusion_mesh_alpha.enriched_to_alpha[i]
            if any(np.array([i_gamma, i_alpha]) != -1): #else this point is on ghost interface gammaalpha (in theta) or it's a triple point (is allowed/must be to pointing to old triple
                if grid.has_point_been_shifted[i, 1] == 1: #1 gammaalpha 2 alphatheta 3 gammatheta 4 triple 11 gammaalpha 'ghost'
                    I_i = grid.diffusion_points_old_gamma[grid.I_diff_gamma[i_gamma]]
                    I_enr_ii = grid.diffusion_points_index_old_gamma[grid.I_diff_gamma[i_gamma]]
                    patch3_old = grid.enriched_three_neighbours_old[I_enr_ii].nonzero()[1]
                    patch3_ga  = patch3_old[np.where(grid.has_point_been_shifted_old[patch3_old,1]   == 1)[0]]
                    patch3_3   = patch3_old[np.where(grid.has_point_been_shifted_old[patch3_old,1]%4 == 0)[0]]
                    patch3_ga  = np.union1d(patch3_ga, patch3_3)
                    edges_ii = np.zeros((0,2), dtype = np.int)
                    for x_j in patch3_ga:
                        patch_x_j = np.setdiff1d(grid.enriched_trii_old[grid.point_in_tri_old[x_j].nonzero()[1]].flatten(), x_j)
                        other_ii_ga = patch_x_j[np.where(grid.has_point_been_shifted_old[patch_x_j,1]   == 1)[0]]
                        other_ii_3  = patch_x_j[np.where(grid.has_point_been_shifted_old[patch_x_j,1]%4 == 0)[0]]
                        other_ii_ga = np.union1d(other_ii_ga, other_ii_3)
                        for x_j_other in other_ii_ga:
                            edges_ii = np.vstack((edges_ii, np.sort(np.hstack((x_j, x_j_other))) ))
                    edge_ii_proj = np.hstack((grid.enriched_mesh.n+1, grid.enriched_mesh.n+2))
                    x_proj = np.array([-10+grid.x0, -10+grid.y0])
                    edges_ii = cGrid2D.cGrid.unique_rows(edges_ii)
                    for edge_ii in edges_ii:
                        edge_i = grid.enriched_points_old[edge_ii]
                        i_projected = la.solve(np.array([[edge_i[0,0] - edge_i[1,0], edge_i[0,1] - edge_i[1,1]],
                                                         [edge_i[1,1] - edge_i[0,1], edge_i[0,0] - edge_i[1,0]]]),
                                              -np.array([-xi[0]  * (edge_i[0,0] - edge_i[1,0]) - xi[1]  * (edge_i[0,1] - edge_i[1,1]),
                                                         -edge_i[1,1] * (edge_i[0,0] - edge_i[1,0]) + edge_i[1,0] * (edge_i[0,1] - edge_i[1,1])]))
                        where_furtherst = np.argmax(la.norm(edge_i - i_projected, 2, axis = 1))
                        distance_xp_furthest = la.norm(i_projected - edge_i[where_furtherst], 2)
                        length_edge = la.norm(edge_i[0] - edge_i[1], 2)
                        if distance_xp_furthest < length_edge:
                            if la.norm(x_proj - xi, 2) > la.norm(i_projected - xi, 2):
                                x_proj = i_projected
                                edge_ii_proj = edge_ii
                                where_furthest_proj = where_furtherst
                                distance_xp_furthest_x_proj = distance_xp_furthest
                                length_edge_x_proj = length_edge
                    #if found projected point is closer to xi than I_i, closest interface point, set concentration of xi to interpolated value on the edge
                    if la.norm(x_proj - xi, 2) < la.norm(I_i - xi, 2):
                        cexn_gamma[i_gamma] = (1 - distance_xp_furthest_x_proj / length_edge_x_proj) * c_gamma[grid.diffusion_enriched_to_gamma_old[edge_ii_proj[where_furthest_proj]]] + distance_xp_furthest_x_proj / length_edge_x_proj * c_gamma[grid.diffusion_enriched_to_gamma_old[edge_ii_proj[1-where_furthest_proj]]]
                        vn_mesh_gamma[i_gamma] = (xi - x_proj)/dt
                        cexn_alpha[i_alpha] = (1 - distance_xp_furthest_x_proj / length_edge_x_proj) * c_alpha[grid.diffusion_enriched_to_alpha_old[edge_ii_proj[where_furthest_proj]]] + distance_xp_furthest_x_proj / length_edge_x_proj * c_alpha[grid.diffusion_enriched_to_alpha_old[edge_ii_proj[1-where_furthest_proj]]]
                        vn_mesh_alpha[i_alpha] = (xi - x_proj)/dt
                        vector_proj = x_proj - xi
                        edge_i_proj = grid.enriched_points_old[edge_ii_proj]
                        vector_edge_proj = edge_i_proj[0] - edge_i_proj[1]
                        if np.abs(vector_proj @ vector_edge_proj) > 1e-14:
                            print('projection0 is not perp, i', i)
                            print('inproduct = ', np.abs(vector_proj @ vector_edge_proj))
                        if distance_xp_furthest_x_proj / length_edge_x_proj > 1 or distance_xp_furthest_x_proj / length_edge_x_proj < 0:
                            print('interpolation factors in extending c incorrect. i =', i)
                    #else put it on the concentration value of the closest interface point I_i, as done in shift
                    else:
                        pass
                elif grid.has_point_been_shifted[i, 1] == 3: #gammatheta
                    I_i = grid.diffusion_points_old_gamma[grid.I_diff_gamma[i_gamma]]
                    I_enr_ii = grid.diffusion_points_index_old_gamma[grid.I_diff_gamma[i_gamma]]
                    patch3_old = grid.enriched_three_neighbours_old[I_enr_ii].nonzero()[1]
                    patch3_gt  = patch3_old[np.where(grid.has_point_been_shifted_old[patch3_old,1]   == 3)[0]]
                    patch3_3   = patch3_old[np.where(grid.has_point_been_shifted_old[patch3_old,1]%4 == 0)[0]]
                    patch3_gt  = np.union1d(patch3_gt, patch3_3)
                    edges_ii = np.zeros((0,2), dtype = np.int)
                    for x_j in patch3_gt:
                        patch_x_j = np.setdiff1d(grid.enriched_trii_old[grid.point_in_tri_old[x_j].nonzero()[1]].flatten(), x_j)
                        other_ii_gt = patch_x_j[np.where(grid.has_point_been_shifted_old[patch_x_j,1]   == 3)[0]]
                        other_ii_3  = patch_x_j[np.where(grid.has_point_been_shifted_old[patch_x_j,1]%4 == 0)[0]]
                        other_ii_gt = np.union1d(other_ii_gt, other_ii_3)
                        for x_j_other in other_ii_gt:
                            edges_ii = np.vstack((edges_ii, np.sort(np.hstack((x_j, x_j_other))) ))
                    edge_ii_proj = np.hstack((grid.enriched_mesh.n+1, grid.enriched_mesh.n+2))
                    x_proj = np.array([-10+grid.x0, -10+grid.y0])
                    edges_ii = cGrid2D.cGrid.unique_rows(edges_ii)
                    for edge_ii in edges_ii:
                        edge_i = grid.enriched_points_old[edge_ii]
                        i_projected = la.solve(np.array([[edge_i[0,0] - edge_i[1,0], edge_i[0,1] - edge_i[1,1]],
                                                  [edge_i[1,1] - edge_i[0,1], edge_i[0,0] - edge_i[1,0]]]),
                                              -np.array([-xi[0]  * (edge_i[0,0] - edge_i[1,0]) - xi[1]  * (edge_i[0,1] - edge_i[1,1]),
                                                    -edge_i[1,1] * (edge_i[0,0] - edge_i[1,0]) + edge_i[1,0] * (edge_i[0,1] - edge_i[1,1])]))
                        where_furtherst = np.argmax(la.norm(edge_i - i_projected, 2, axis = 1))
                        distance_xp_furthest = la.norm(i_projected - edge_i[where_furtherst], 2)
                        length_edge = la.norm(edge_i[0] - edge_i[1], 2)
                        if distance_xp_furthest < length_edge:
                            if la.norm(x_proj - xi, 2) > la.norm(i_projected - xi, 2):
                                x_proj = i_projected
                                edge_ii_proj = edge_ii
                                where_furthest_proj = where_furtherst
                                distance_xp_furthest_x_proj = distance_xp_furthest
                                length_edge_x_proj = length_edge
                    #if found projected point is closer to xi than I_i, closest interface point, set concentration of xi to interpolated value on the edge
                    if la.norm(x_proj - xi, 2) < la.norm(I_i - xi):
                        cexn_gamma[i_gamma] = (1 - distance_xp_furthest_x_proj / length_edge_x_proj) * c_gamma[grid.diffusion_enriched_to_gamma_old[edge_ii_proj[where_furthest_proj]]] + distance_xp_furthest_x_proj / length_edge_x_proj * c_gamma[grid.diffusion_enriched_to_gamma_old[edge_ii_proj[1-where_furthest_proj]]]
                        vn_mesh_gamma[i_gamma] = (xi - x_proj)/dt
                        vector_proj = x_proj - xi
                        edge_i_proj = grid.enriched_points_old[edge_ii_proj]
                        vector_edge_proj = edge_i_proj[0] - edge_i_proj[1]
                        if np.abs(vector_proj @ vector_edge_proj) > 1e-14:
                            print('projection0 is not perp, i', i)
                            print('inproduct = ', np.abs(vector_proj @ vector_edge_proj))
                        if distance_xp_furthest_x_proj / length_edge_x_proj > 1 or distance_xp_furthest_x_proj / length_edge_x_proj < 0:
                            print('interpolation factors in extending c incorrect. i =', i)
                    #else put it on the concentration value of the closest interface point I_i, as done in shift
                    else:
                        pass
                elif grid.has_point_been_shifted[i, 1] == 2: #alphatheta
                    I_i = grid.diffusion_points_old_alpha[grid.I_diff_alpha[i_alpha]]
                    I_enr_ii = grid.diffusion_points_index_old_alpha[grid.I_diff_alpha[i_alpha]]
                    patch3_old = grid.enriched_three_neighbours_old[I_enr_ii].nonzero()[1]
                    patch3_at  = patch3_old[np.where(grid.has_point_been_shifted_old[patch3_old,1]   == 2)[0]]
                    patch3_3   = patch3_old[np.where(grid.has_point_been_shifted_old[patch3_old,1]%4 == 0)[0]]
                    patch3_at  = np.union1d(patch3_at, patch3_3)
                    edges_ii = np.zeros((0,2), dtype = np.int)
                    for x_j in patch3_at:
                        patch_x_j = np.setdiff1d(grid.enriched_trii_old[grid.point_in_tri_old[x_j].nonzero()[1]].flatten(), x_j)
                        other_ii_at = patch_x_j[np.where(grid.has_point_been_shifted_old[patch_x_j,1]   == 2)[0]]
                        other_ii_3  = patch_x_j[np.where(grid.has_point_been_shifted_old[patch_x_j,1]%4 == 0)[0]]
                        other_ii_at = np.union1d(other_ii_at, other_ii_3)
                        for x_j_other in other_ii_at:
                            edges_ii = np.vstack((edges_ii, np.sort(np.hstack((x_j, x_j_other))) ))
                    edge_ii_proj = np.hstack((grid.enriched_mesh.n+1, grid.enriched_mesh.n+2))
                    x_proj = np.array([-10+grid.x0, -10+grid.y0])
                    edges_ii = cGrid2D.cGrid.unique_rows(edges_ii)
                    for edge_ii in edges_ii:
                        edge_i = grid.enriched_points_old[edge_ii]
                        i_projected = la.solve(np.array([[edge_i[0,0] - edge_i[1,0], edge_i[0,1] - edge_i[1,1]],
                                                  [edge_i[1,1] - edge_i[0,1], edge_i[0,0] - edge_i[1,0]]]),
                                              -np.array([-xi[0]  * (edge_i[0,0] - edge_i[1,0]) - xi[1]  * (edge_i[0,1] - edge_i[1,1]),
                                                    -edge_i[1,1] * (edge_i[0,0] - edge_i[1,0]) + edge_i[1,0] * (edge_i[0,1] - edge_i[1,1])]))
                        where_furtherst = np.argmax(la.norm(edge_i - i_projected, 2, axis = 1))
                        distance_xp_furthest = la.norm(i_projected - edge_i[where_furtherst], 2)
                        length_edge = la.norm(edge_i[0] - edge_i[1], 2)
                        if distance_xp_furthest < length_edge:
                            if la.norm(x_proj - xi, 2) > la.norm(i_projected - xi, 2):
                                x_proj = i_projected
                                edge_ii_proj = edge_ii
                                where_furthest_proj = where_furtherst
                                distance_xp_furthest_x_proj = distance_xp_furthest
                                length_edge_x_proj = length_edge
                    #if found projected point is closer to xi than I_i, closest interface point, set concentration of xi to interpolated value on the edge
                    if la.norm(x_proj - xi, 2) < la.norm(I_i - xi):
                        cexn_alpha[i_alpha] = (1 - distance_xp_furthest_x_proj / length_edge_x_proj) * c_alpha[grid.diffusion_enriched_to_alpha_old[edge_ii_proj[where_furthest_proj]]] + distance_xp_furthest_x_proj / length_edge_x_proj * c_alpha[grid.diffusion_enriched_to_alpha_old[edge_ii_proj[1-where_furthest_proj]]]
                        vn_mesh_alpha[i_alpha] = (xi - x_proj)/dt
                        vector_proj = x_proj - xi
                        edge_i_proj = grid.enriched_points_old[edge_ii_proj]
                        vector_edge_proj = edge_i_proj[0] - edge_i_proj[1]
                        if np.abs(vector_proj @ vector_edge_proj) > 1e-14:
                            print('projection0 is not perp, i', i)
                            print('inproduct = ', np.abs(vector_proj @ vector_edge_proj))
                        if distance_xp_furthest_x_proj / length_edge_x_proj > 1 or distance_xp_furthest_x_proj / length_edge_x_proj < 0:
                            print('interpolation factors in extending c incorrect. i =', i)
                    #else put it on the concentration value of the closest interface point I_i. as done in shift
                    else:
                        pass

        if ext == 'Interpolate':
            ext_interpolate = np.zeros((0,5), dtype = np.int)
            gamma_background = grid.diffusion_mesh_gamma.points_index[np.where(grid.diffusion_mesh_gamma.points_index < grid.N)]
            where_new_gamma, = np.nonzero(grid.has_point_been_shifted[gamma_background,0] - grid.has_point_been_shifted_old[gamma_background,0])
            alpha_background = grid.diffusion_mesh_alpha.points_index[np.where(grid.diffusion_mesh_alpha.points_index < grid.N)]
            where_new_alpha, = np.nonzero(grid.has_point_been_shifted[alpha_background,0] - grid.has_point_been_shifted_old[alpha_background,0])            
            for i in where_new_gamma:
                #'new' points where concentration is 'unknown' (not on interface( are only points that are from the background mesh (index < n(N))
                if grid.diffusion_mesh_gamma.points_index[i] < grid.N:
                    pi = grid.diffusion_mesh_gamma.points[i]
                    enr_i = grid.diffusion_mesh_gamma.points_index[i]
                    pi_old = pi - vn_ex_ga[enr_i] * dt
                    for tri_i in grid.point_in_tri_old[enr_i].nonzero()[1]:
                        try:
                            trii = grid.enriched_trii_old[tri_i]
                            trii_gamma = grid.diffusion_enriched_to_gamma_old[trii]
                            tri = grid.enriched_points_old[trii]
                            lambda_barycentric = la.solve( np.array([[tri[0,0], tri[1,0], tri[2,0]],
                                                                     [tri[0,1], tri[1,1], tri[2,1]],
                                                                     [1       , 1       , 1       ] ]), np.array([ pi_old[0], pi_old[1], 1]))
                            if all(lambda_barycentric <= 1) and all(lambda_barycentric >= 0):
                                ext_interpolate = np.vstack((ext_interpolate, np.hstack([enr_i, i, trii])))
                                cexn_gamma[i] = interpolate(trii_gamma, tri, pi_old, c_gamma)
                                vn_mesh_gamma[i] = (pi-pi_old)/dt

                        except:
                            pass
            ext_interpolate2 = np.zeros((0,5), dtype = np.int)
            for i in where_new_alpha:
                #'new' points where concentration is 'unknown' (not on interface( are only points that are from the background mesh (index < n(N))
                if grid.diffusion_mesh_alpha.points_index[i] < grid.N:
                    pi = grid.diffusion_mesh_alpha.points[i]
                    enr_i = grid.diffusion_mesh_alpha.points_index[i]
                    pi_old = pi - vn_ex_t[enr_i] * dt
                    for tri_i in grid.point_in_tri_old[enr_i].nonzero()[1]:
                        try:
                            trii = grid.enriched_trii_old[tri_i]
                            trii_alpha = grid.diffusion_enriched_to_alpha_old[trii]
                            tri = grid.enriched_points_old[trii]
                            lambda_barycentric = la.solve( np.array([[tri[0,0], tri[1,0], tri[2,0]],
                                                                     [tri[0,1], tri[1,1], tri[2,1]],
                                                                     [1       , 1       , 1       ] ]), np.array([ pi_old[0], pi_old[1], 1]))
                            if all(lambda_barycentric <= 1) and all(lambda_barycentric >= 0):
                                ext_interpolate2 = np.vstack((ext_interpolate2, np.hstack([enr_i, i, trii])))
                                cexn_alpha[i] = interpolate(trii_alpha, tri, pi_old, c_alpha)
                                vn_mesh_alpha[i] = (pi-pi_old)/dt
                        except:
                            pass



    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_alpha = 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)
    Diff_V = np.zeros(9*grid.diffusion_mesh_alpha.n_tri)
    Conv_V = np.zeros(9*grid.diffusion_mesh_alpha.n_tri)
    M_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]

        beta = np.zeros(3)
        beta[0] = (tri[1,1] - tri[2,1])/Delta
        beta[1] = (tri[2,1] - tri[0,1])/Delta
        beta[2] = (tri[0,1] - tri[1,1])/Delta
        gamma = np.zeros(3)
        gamma[0] = (tri[2,0] - tri[1,0])/Delta
        gamma[1] = (tri[0,0] - tri[2,0])/Delta
        gamma[2] = (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))

            #Diffussive part
        Diff_V[9*(l-1):9*(l-1)+9] = - param.D_alpha(T) * np.abs(Delta) * ((beta.reshape(3,1) @ beta.reshape(1,3)) + (gamma.reshape(3,1) @ gamma.reshape(1,3))).reshape(9)

            #Convective part
        Conv_V[9*(l-1):9*(l-1)+9] = np.abs(Delta) / 24 *  np.array([2*(vn_mesh_alpha[trii_alpha_coord[0],0] * beta[0] + vn_mesh_alpha[trii_alpha_coord[0],1] * gamma[0]) + vn_mesh_alpha[trii_alpha_coord[1],0] * beta[0] + vn_mesh_alpha[trii_alpha_coord[1],1] * gamma[0] + vn_mesh_alpha[trii_alpha_coord[2],0] * beta[0] + vn_mesh_alpha[trii_alpha_coord[2],1] * gamma[0],
                                                                        2*(vn_mesh_alpha[trii_alpha_coord[0],0] * beta[1] + vn_mesh_alpha[trii_alpha_coord[0],1] * gamma[1]) + vn_mesh_alpha[trii_alpha_coord[1],0] * beta[1] + vn_mesh_alpha[trii_alpha_coord[1],1] * gamma[1] + vn_mesh_alpha[trii_alpha_coord[2],0] * beta[1] + vn_mesh_alpha[trii_alpha_coord[2],1] * gamma[1],
                                                                        2*(vn_mesh_alpha[trii_alpha_coord[0],0] * beta[2] + vn_mesh_alpha[trii_alpha_coord[0],1] * gamma[2]) + vn_mesh_alpha[trii_alpha_coord[1],0] * beta[2] + vn_mesh_alpha[trii_alpha_coord[1],1] * gamma[2] + vn_mesh_alpha[trii_alpha_coord[2],0] * beta[2] + vn_mesh_alpha[trii_alpha_coord[2],1] * gamma[2],
                                                                        2*(vn_mesh_alpha[trii_alpha_coord[1],0] * beta[0] + vn_mesh_alpha[trii_alpha_coord[1],1] * gamma[0]) + vn_mesh_alpha[trii_alpha_coord[0],0] * beta[0] + vn_mesh_alpha[trii_alpha_coord[0],1] * gamma[0] + vn_mesh_alpha[trii_alpha_coord[2],0] * beta[0] + vn_mesh_alpha[trii_alpha_coord[2],1] * gamma[0],
                                                                        2*(vn_mesh_alpha[trii_alpha_coord[1],0] * beta[1] + vn_mesh_alpha[trii_alpha_coord[1],1] * gamma[1]) + vn_mesh_alpha[trii_alpha_coord[0],0] * beta[1] + vn_mesh_alpha[trii_alpha_coord[0],1] * gamma[1] + vn_mesh_alpha[trii_alpha_coord[2],0] * beta[1] + vn_mesh_alpha[trii_alpha_coord[2],1] * gamma[1],
                                                                        2*(vn_mesh_alpha[trii_alpha_coord[1],0] * beta[2] + vn_mesh_alpha[trii_alpha_coord[1],1] * gamma[2]) + vn_mesh_alpha[trii_alpha_coord[0],0] * beta[2] + vn_mesh_alpha[trii_alpha_coord[0],1] * gamma[2] + vn_mesh_alpha[trii_alpha_coord[2],0] * beta[2] + vn_mesh_alpha[trii_alpha_coord[2],1] * gamma[2],
                                                                        2*(vn_mesh_alpha[trii_alpha_coord[2],0] * beta[0] + vn_mesh_alpha[trii_alpha_coord[2],1] * gamma[0]) + vn_mesh_alpha[trii_alpha_coord[0],0] * beta[0] + vn_mesh_alpha[trii_alpha_coord[0],1] * gamma[0] + vn_mesh_alpha[trii_alpha_coord[1],0] * beta[0] + vn_mesh_alpha[trii_alpha_coord[1],1] * gamma[0],
                                                                        2*(vn_mesh_alpha[trii_alpha_coord[2],0] * beta[1] + vn_mesh_alpha[trii_alpha_coord[2],1] * gamma[1]) + vn_mesh_alpha[trii_alpha_coord[0],0] * beta[1] + vn_mesh_alpha[trii_alpha_coord[0],1] * gamma[1] + vn_mesh_alpha[trii_alpha_coord[1],0] * beta[1] + vn_mesh_alpha[trii_alpha_coord[1],1] * gamma[1],
                                                                        2*(vn_mesh_alpha[trii_alpha_coord[2],0] * beta[2] + vn_mesh_alpha[trii_alpha_coord[2],1] * gamma[2]) + vn_mesh_alpha[trii_alpha_coord[0],0] * beta[2] + vn_mesh_alpha[trii_alpha_coord[0],1] * gamma[2] + vn_mesh_alpha[trii_alpha_coord[1],0] * beta[2] + vn_mesh_alpha[trii_alpha_coord[1],1] * gamma[2]])
            #Mass matrix
        M_V[9*(l-1):9*(l-1)+9]    = np.abs(Delta) / 24 *  (np.ones((3,3)) + np.eye(3)).reshape(9)

        l = l+1

    Diff_alpha = coo_matrix((Diff_V,(I,J)), shape = (grid.diffusion_mesh_alpha.n,grid.diffusion_mesh_alpha.n))
    Diff_alpha = Diff_alpha.tocsr()


    Conv_alpha = coo_matrix((Conv_V,(I,J)), shape = (grid.diffusion_mesh_alpha.n,grid.diffusion_mesh_alpha.n))
    Conv_alpha = Conv_alpha.tocsr()

    M_alpha = coo_matrix((M_V,(I,J)), shape = (grid.diffusion_mesh_alpha.n,grid.diffusion_mesh_alpha.n))

    M_alpha = M_alpha.tocsr()
    M_L_alpha = diags(np.asarray(csr_matrix.sum(M_alpha, axis = 1)).ravel(), format = 'csr')

    Mcexn = M_alpha @ cexn_alpha
    M_Lcexn = M_L_alpha @ cexn_alpha

    I = np.zeros(4*(grid.interface_alphatheta.n_edges), dtype = np.int)
    J = np.zeros(4*(grid.interface_alphatheta.n_edges), dtype = np.int)
    R_V = np.zeros(4*(grid.interface_alphatheta.n_edges))
    f_R = 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]
        I[4*(l-1):4*(l-1)+4] = np.array([el_a[0], el_a[0], el_a[1], el_a[1]])
        J[4*(l-1):4*(l-1)+4] = np.array([el_a[0], el_a[1], el_a[0], el_a[1]])
        R_V[4*(l-1):4*(l-1)+4] = L/6*np.array([2, 1, 1, 2])
        l = l + 1
        f_R[grid.diffusion_mesh_alpha.enriched_to_alpha[el]] = f_R[grid.diffusion_mesh_alpha.enriched_to_alpha[el]] + L/2 * np.ones(2)
    R = coo_matrix((R_V,(I,J)), shape = (grid.diffusion_mesh_alpha.n,grid.diffusion_mesh_alpha.n))
    R = R.tocsr()

    #f(t^n+1,c^n+1)
    f_alpha = np.zeros(grid.diffusion_mesh_alpha.n)
    if not(Curvature):
        f_alpha[ag_points] = - param.c_alphagamma(T)
    else:
        for i in range(0, grid.interface_gammaalpha.n):
            f_alpha[ag_points[i]] = - param.c_alphagamma_GT(T, -grid.enriched_mesh.curvature_ga[i])


    ### Fluc Correction
    i = 0
    neighbours = lil_matrix((grid.diffusion_mesh_alpha.n, grid.diffusion_mesh_alpha.n), dtype = bool)
    for j in grid.diffusion_mesh_alpha.points_index:
        other_ii_all = np.unique(grid.enriched_mesh.triangles_points_index[grid.enriched_mesh.point_in_tri[j].nonzero()[1]].flatten())
        #alpha
        other_ii_alpha = grid.diffusion_mesh_alpha.enriched_to_alpha[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) ])]
        neighbours[i, other_ii_alpha] = True
        i = i+1    

    #Antidifussive part of diffusion matrix
    Diff_plus_alpha = lil_matrix((grid.diffusion_mesh_alpha.n, grid.diffusion_mesh_alpha.n))
    where_neg = find(Diff_alpha < 0)
    Diff_plus_alpha[where_neg[0], where_neg[1]] = Diff_alpha[where_neg[0], where_neg[1]]
    Diff_plus_alpha.setdiag(0)
    Diff_plus_alpha.setdiag(-lil_matrix.sum(Diff_plus_alpha, axis = 1))
    Diff_plus_alpha = Diff_plus_alpha.tocsr()
    Diff_tilde_alpha = Diff_alpha - Diff_plus_alpha

    #Artificial diffusion matrix D to limit Conv_alpha
    D = lil_matrix((grid.diffusion_mesh_alpha.n, grid.diffusion_mesh_alpha.n))
    for i in np.unique(Conv_alpha.indices):
        where_Conv_neg = np.setdiff1d(find(Conv_alpha[i] < 0)[1], i)
        for j_neg in where_Conv_neg:
            min_neg = max(-Conv_alpha[i,j_neg], 0, -Conv_alpha[j_neg,i])
            D[i,j_neg] = min_neg
            D[j_neg,i] = min_neg
    D.setdiag(-lil_matrix.sum(D, axis = 1))
    D = D.tocsr()
    Conv_tilde_alpha = Conv_alpha + D

    r = 1
    eps = 1e6
    FPT = 1e-6
    qr_tilde = cexn_alpha.copy()
    qr_tilde[at_points] = param.c_alphatheta(T)
    qr_tilde[ag_points] = param.c_gammaalpha(T)
    qr_H_alpha = cexn_alpha.copy()
    qr_H_alpha[at_points] = param.c_alphatheta(T)
    qr_H_alpha[ag_points] = param.c_gammaalpha(T)
    dqrdt_tilde = np.ones(grid.diffusion_mesh_alpha.n)
    M_alpha = M_alpha.tocoo()
    while eps >= FPT:
        R_cp_V = np.zeros(4*(grid.interface_alphatheta.n_edges))
        f_R_cp = 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]], 2)
            el_a = grid.diffusion_mesh_alpha.enriched_to_alpha[el]
#            ### Not lumped Reaction parts
#            R_cp_V[4*(l-1):4*(l-1)+4] = L/12*np.array([3*qr_tilde[el_a[0]] + qr_tilde[el_a[1]], qr_tilde[el_a[0]] +   qr_tilde[el_a[1]],
#                                                         qr_tilde[el_a[0]] + qr_tilde[el_a[1]], qr_tilde[el_a[0]] + 3*qr_tilde[el_a[1]]])
#            f_R_cp[el_a] = f_R_cp[el_a] + L/6 * np.array([2*qr_tilde[el_a[0]] +   qr_tilde[el_a[1]],
#                                                            qr_tilde[el_a[0]] + 2*qr_tilde[el_a[1]] ])
            ### Lumped Reaction parts
            R_cp_V[4*(l-1):4*(l-1)+4] = L/12*np.array([3*qr_tilde[el_a[0]] + qr_tilde[el_a[0]], qr_tilde[el_a[0]] +   qr_tilde[el_a[0]],
                                                         qr_tilde[el_a[1]] + qr_tilde[el_a[1]], qr_tilde[el_a[1]] + 3*qr_tilde[el_a[1]]])
            f_R_cp[el_a] = f_R_cp[el_a] + L/2 * np.array([qr_tilde[el_a[0]],
                                                          qr_tilde[el_a[1]] ])
    
            l = l + 1

            
        R_cp = coo_matrix((R_cp_V,(I,J)), shape = (grid.diffusion_mesh_alpha.n, grid.diffusion_mesh_alpha.n))
        R_cp = R_cp.tocsr()

        #Reaction matrix BC
        React_alpha      = - param.K_alphatheta(T) * (param.c0_theta() * R[at_points] - R_cp[at_points]) / param.c0_theta()
        React_alpha_matrix = lil_matrix((grid.diffusion_mesh_alpha.n, grid.diffusion_mesh_alpha.n))
        React_alpha_matrix[at_points] = React_alpha
        React_alpha_matrix = React_alpha_matrix.tocsr()
        React_alpha_matrix_lumped = diags(np.asarray(csr_matrix.sum(React_alpha_matrix, axis = 1)).ravel(), format = 'csr')

        if not(Curvature):
            f_alpha[at_points] = - param.K_alphatheta(T) * param.c_alphatheta(T) * (param.c0_theta() * f_R[at_points] - f_R_cp[at_points]) / param.c0_theta()
            f_alpha[x_points_alpha] = - param.c_alphagamma(T)
        else:
            for i in range(0, grid.interface_alphatheta.n):
                f_alpha[at_points[i]] = - param.K_alphatheta(T) * param.c_alphatheta_GT(T, grid.enriched_mesh.curvature_at[i]) * (param.c0_theta() * f_R[at_points[i]] - f_R_cp[at_points[i]]) / param.c0_theta()
            for i in range(0, len(grid.triple_points_index)):
                x_i = grid.interface_gammaalpha.triple_points_index[i]
                f_alpha[i] = - param.c_alphagamma_GT(T, -grid.enriched_mesh.curvature_ga[x_i])
        B_tilde_alpha = (M_L_alpha - dt*(Diff_tilde_alpha + Conv_tilde_alpha + React_alpha_matrix_lumped)).tocsr() #React_alpha_matrix_lumped
        ### Diriclet BC
        [setcsrrow2id(B_tilde_alpha, ag_points[i]) for i in range(0, len(ag_points))]
        g_alpha = M_Lcexn - dt*f_alpha
        if not(Curvature):
            g_alpha[ag_points] = param.c_alphagamma(T)
        else:
            for i in range(0, grid.interface_gammaalpha.n):
                g_alpha[ag_points[i]] = param.c_alphagamma_GT(T, -grid.enriched_mesh.curvature_ga[i])

        eps_qr_tilde = la.norm(B_tilde_alpha @ qr_tilde - g_alpha, np.inf) / la.norm(g_alpha, np.inf)

        qr_tilde = spla.spsolve(B_tilde_alpha, g_alpha)


        #Approximate d/dt u_i
        Su = (Conv_tilde_alpha + Diff_tilde_alpha + React_alpha_matrix_lumped) @ qr_tilde - f_alpha
        Su[ag_points] = 0

        eps_dqrdt_tilde = la.norm(M_L_alpha @ dqrdt_tilde - Su, np.inf) / la.norm(Su, np.inf)

        dqrdt_tilde = spla.spsolve(M_L_alpha, Su)#tilde
#        #Approximate d/dt u_i, less diffusive way
#        Su = (Conv_alpha + Diff_alpha + React_alpha_matrix_lumped) @ qr_tilde - f_alpha
#        Su[ag_points] = 0
#        dqrdt_tilde = spla.spsolve(M_alpha.tocsr(), Su)



        R_cp_V = np.zeros(4*(grid.interface_alphatheta.n_edges))
        f_R_cp = 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]], 2)
            el_a = grid.diffusion_mesh_alpha.enriched_to_alpha[el]
            ### Not lumped Reaction parts
#            R_cp_V[4*(l-1):4*(l-1)+4] = L/12*np.array([3*qr_H_alpha[el_a[0]] + qr_H_alpha[el_a[1]], qr_H_alpha[el_a[0]] +   qr_H_alpha[el_a[1]],
#                                                         qr_H_alpha[el_a[0]] + qr_H_alpha[el_a[1]], qr_H_alpha[el_a[0]] + 3*qr_H_alpha[el_a[1]]])
#            f_R_cp[el_a] = f_R_cp[el_a] + L/6 * np.array([2*qr_H_alpha[el_a[0]] +   qr_H_alpha[el_a[1]],
#                                                            qr_H_alpha[el_a[0]] + 2*qr_H_alpha[el_a[1]] ])

            ### Lumped Reaction parts
            R_cp_V[4*(l-1):4*(l-1)+4] = L/12*np.array([3*qr_H_alpha[el_a[0]] + qr_H_alpha[el_a[0]], qr_H_alpha[el_a[0]] +   qr_H_alpha[el_a[0]],
                                                         qr_H_alpha[el_a[1]] + qr_H_alpha[el_a[1]], qr_H_alpha[el_a[1]] + 3*qr_H_alpha[el_a[1]]])
            f_R_cp[el_a] = f_R_cp[el_a] + L/2 * np.array([qr_H_alpha[el_a[0]],
                                                          qr_H_alpha[el_a[1]] ])
    
            l = l + 1

            
        R_cp = coo_matrix((R_cp_V,(I,J)), shape = (grid.diffusion_mesh_alpha.n, grid.diffusion_mesh_alpha.n))
        R_cp = R_cp.tocsr()

        #Reaction matrix BC
        React_alpha      = - param.K_alphatheta(T) * (param.c0_theta() * R[at_points] - R_cp[at_points]) / param.c0_theta()
        React_alpha_matrix = lil_matrix((grid.diffusion_mesh_alpha.n, grid.diffusion_mesh_alpha.n))
        React_alpha_matrix[at_points] = React_alpha
        React_alpha_matrix = React_alpha_matrix.tocsr()
        React_alpha_matrix_lumped = diags(np.asarray(csr_matrix.sum(React_alpha_matrix, axis = 1)).ravel(), format = 'csr')


        B_alpha = (M_alpha.tocsr() - dt*(Diff_alpha + Conv_alpha + React_alpha_matrix_lumped)).tocsr() #React_alpha_matrix_lumped
        [setcsrrow2id(B_alpha, ag_points[i]) for i in range(0, len(ag_points))]
        g_alpha = Mcexn - dt*f_alpha
        if not(Curvature):
            g_alpha[ag_points] = param.c_alphagamma(T)
        else:
            for i in range(0, grid.interface_gammaalpha.n):
                g_alpha[ag_points[i]] = param.c_alphagamma_GT(T, -grid.enriched_mesh.curvature_ga[i])

        eps_qr_H = la.norm(B_alpha @ qr_H_alpha - g_alpha, np.inf) / la.norm(g_alpha, np.inf)

        qr_H_alpha = spla.spsolve(B_alpha, g_alpha)

        eps = eps_qr_H
        r=r+1
        print('eps alpha %.3e' % eps, 'qr_tilde, dqrtildedt %.3e %.3e' %(eps_qr_tilde ,eps_dqrdt_tilde))

    ### Flux Correction
    #flux f_ij = (mij d/dt + d_ij + l+_ij)(u_i-u_j) c^n+1 approx bij qr_tilde, l+ should be zero
    flux_M_alpha_V = np.zeros(M_alpha.nnz)
    for ij in range(0, M_alpha.nnz):
        i = M_alpha.row[ij]
        j = M_alpha.col[ij]
        flux_M_alpha_V[ij] = M_alpha.data[ij] * (dqrdt_tilde[i] - dqrdt_tilde[j])
    DDiff_plus_alpha = (D-Diff_plus_alpha).tocoo()
    flux_DDiff_plus_alpha_V = np.zeros(DDiff_plus_alpha.nnz)
    for ij in range(0, DDiff_plus_alpha.nnz):
        i = DDiff_plus_alpha.row[ij]
        j = DDiff_plus_alpha.col[ij]
        flux_DDiff_plus_alpha_V[ij] = DDiff_plus_alpha.data[ij] * (qr_tilde[i] - qr_tilde[j])
    flux = (coo_matrix((flux_M_alpha_V,(M_alpha.row,M_alpha.col)), shape = (grid.diffusion_mesh_alpha.n,grid.diffusion_mesh_alpha.n)).tocsr() +
            coo_matrix((flux_DDiff_plus_alpha_V,(DDiff_plus_alpha.row,DDiff_plus_alpha.col)), shape = (grid.diffusion_mesh_alpha.n,grid.diffusion_mesh_alpha.n)).tocsr()).tolil()

    flux.setdiag(0)
    flux.setdiag(lil_matrix.sum(flux, axis = 1))
    flux = flux.tocoo()

    for ij in range(0, flux.nnz):
        i = flux.row[ij]
        j = flux.col[ij]
        if flux.data[ij] * (qr_tilde[j] - qr_tilde[i]) > 0:
            flux.data[ij] = 0
    flux = flux.tocsr()



    qr_max = -np.ones(grid.diffusion_mesh_alpha.n)
    qr_min = -np.ones(grid.diffusion_mesh_alpha.n)
    for i in range(0, grid.diffusion_mesh_alpha.n):
        qr_max[i] = np.max(qr_tilde[np.asarray(neighbours[i].todense()).reshape(grid.diffusion_mesh_alpha.n)])
        qr_min[i] = np.min(qr_tilde[np.asarray(neighbours[i].todense()).reshape(grid.diffusion_mesh_alpha.n)])

    if not(Curvature):
        qr_min[at_points] = np.maximum(qr_min[at_points], param.c_alphatheta(T))
    else:
        for i in range(0, grid.interface_alphatheta.n):
            qr_min[at_points[i]] = np.maximum(qr_min[at_points[i]], param.c_alphatheta_GT(T, grid.enriched_mesh.curvature_at[i]))

    R_min = np.ones(grid.diffusion_mesh_alpha.n)
    R_plus = np.ones(grid.diffusion_mesh_alpha.n)
    nonzero = flux.nonzero()
    nonzero_i = np.setdiff1d(np.unique(nonzero[0]), ag_points)
    for i in nonzero_i:
        row_flux_pos = (flux[i] + np.abs(flux[i]))/2
        row_flux_neg = (flux[i] - np.abs(flux[i]))/2
        P_plus = csr_matrix.sum(row_flux_pos) - row_flux_pos[0,i]
        P_min  = csr_matrix.sum(row_flux_neg) - row_flux_neg[0,i]

        Q_plus = M_L_alpha[i,i] / dt * (qr_max[i] - qr_tilde[i])
        Q_min  = M_L_alpha[i,i] / dt * (qr_min[i] - qr_tilde[i])

        if P_plus != 0:
            R_plus[i] = np.min([1, Q_plus/P_plus])
        if P_min != 0:
            R_min[i]  = np.min([1, Q_min/P_min])


    coeff = lil_matrix((grid.diffusion_mesh_alpha.n, grid.diffusion_mesh_alpha.n))
    for ki in range(0,len(nonzero[0])):
        i = nonzero[0][ki]
        j = nonzero[1][ki]
        if flux[i,j] > 0:
            coeff[i,j] = np.min([R_plus[i], R_min[j]])
        elif flux[i,j] < 0:
            coeff[i,j] = np.min([R_min[i], R_plus[j]])
    coeff = coeff.tocsr()
    fbar = coeff.multiply(flux)
    fbar = fbar.tolil()
    fbar.setdiag(0)
    fbar = lil_matrix.sum(fbar, axis = 1)
    fbar[ag_points] = 0
    fbar = np.squeeze(np.asarray(fbar))

    qr = spla.spsolve(M_L_alpha, M_L_alpha @ qr_tilde + dt * fbar)



    c_alpha = qr.copy()
    
    
    ### Calculate c_gamma
    c_gamma = np.zeros(0)
    dcdn_alpha_ga = np.zeros(0)
    if grid.diffusion_mesh_gamma.n > 0:
        #Calc dc_alpha/dn
        ga_points = grid.interface_gammaalpha.points_index
        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]

        dcdxy_alpha_ga = np.zeros((grid.interface_gammaalpha.n,2))
        ki=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
            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[ki] = B/hp
            ki=ki+1
        dcdn_alpha_ga = np.sum(dcdxy_alpha_ga * grid.n_ga[ga_points], axis = 1)

        ga_points = grid.diffusion_mesh_gamma.enriched_to_gamma[grid.interface_gammaalpha.points_index]
        dcdn_alpha_matrix = lil_matrix((grid.diffusion_mesh_gamma.n, grid.diffusion_mesh_gamma.n))


        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)
        Diff_V = np.zeros(9*grid.diffusion_mesh_gamma.n_tri)
        Conv_V = np.zeros(9*grid.diffusion_mesh_gamma.n_tri)
        M_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]

            beta = np.zeros(3)
            beta[0] = (tri[1,1] - tri[2,1])/Delta
            beta[1] = (tri[2,1] - tri[0,1])/Delta
            beta[2] = (tri[0,1] - tri[1,1])/Delta
            gamma = np.zeros(3)
            gamma[0] = (tri[2,0] - tri[1,0])/Delta
            gamma[1] = (tri[0,0] - tri[2,0])/Delta
            gamma[2] = (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))

            #Diffussive part
            Diff_V[9*(l-1):9*(l-1)+9]      = - param.D_gamma(T) * np.abs(Delta) * ((beta.reshape(3,1) @ beta.reshape(1,3)) + (gamma.reshape(3,1) @ gamma.reshape(1,3))).reshape(9)

            #Convective part
            Conv_V[9*(l-1):9*(l-1)+9] = np.abs(Delta) / 24 *  np.array([2*(vn_mesh_gamma[trii_gamma_coord[0],0] * beta[0] + vn_mesh_gamma[trii_gamma_coord[0],1] * gamma[0]) + vn_mesh_gamma[trii_gamma_coord[1],0] * beta[0] + vn_mesh_gamma[trii_gamma_coord[1],1] * gamma[0] + vn_mesh_gamma[trii_gamma_coord[2],0] * beta[0] + vn_mesh_gamma[trii_gamma_coord[2],1] * gamma[0],
                                                                            2*(vn_mesh_gamma[trii_gamma_coord[0],0] * beta[1] + vn_mesh_gamma[trii_gamma_coord[0],1] * gamma[1]) + vn_mesh_gamma[trii_gamma_coord[1],0] * beta[1] + vn_mesh_gamma[trii_gamma_coord[1],1] * gamma[1] + vn_mesh_gamma[trii_gamma_coord[2],0] * beta[1] + vn_mesh_gamma[trii_gamma_coord[2],1] * gamma[1],
                                                                            2*(vn_mesh_gamma[trii_gamma_coord[0],0] * beta[2] + vn_mesh_gamma[trii_gamma_coord[0],1] * gamma[2]) + vn_mesh_gamma[trii_gamma_coord[1],0] * beta[2] + vn_mesh_gamma[trii_gamma_coord[1],1] * gamma[2] + vn_mesh_gamma[trii_gamma_coord[2],0] * beta[2] + vn_mesh_gamma[trii_gamma_coord[2],1] * gamma[2],
                                                                            2*(vn_mesh_gamma[trii_gamma_coord[1],0] * beta[0] + vn_mesh_gamma[trii_gamma_coord[1],1] * gamma[0]) + vn_mesh_gamma[trii_gamma_coord[0],0] * beta[0] + vn_mesh_gamma[trii_gamma_coord[0],1] * gamma[0] + vn_mesh_gamma[trii_gamma_coord[2],0] * beta[0] + vn_mesh_gamma[trii_gamma_coord[2],1] * gamma[0],
                                                                            2*(vn_mesh_gamma[trii_gamma_coord[1],0] * beta[1] + vn_mesh_gamma[trii_gamma_coord[1],1] * gamma[1]) + vn_mesh_gamma[trii_gamma_coord[0],0] * beta[1] + vn_mesh_gamma[trii_gamma_coord[0],1] * gamma[1] + vn_mesh_gamma[trii_gamma_coord[2],0] * beta[1] + vn_mesh_gamma[trii_gamma_coord[2],1] * gamma[1],
                                                                            2*(vn_mesh_gamma[trii_gamma_coord[1],0] * beta[2] + vn_mesh_gamma[trii_gamma_coord[1],1] * gamma[2]) + vn_mesh_gamma[trii_gamma_coord[0],0] * beta[2] + vn_mesh_gamma[trii_gamma_coord[0],1] * gamma[2] + vn_mesh_gamma[trii_gamma_coord[2],0] * beta[2] + vn_mesh_gamma[trii_gamma_coord[2],1] * gamma[2],
                                                                            2*(vn_mesh_gamma[trii_gamma_coord[2],0] * beta[0] + vn_mesh_gamma[trii_gamma_coord[2],1] * gamma[0]) + vn_mesh_gamma[trii_gamma_coord[0],0] * beta[0] + vn_mesh_gamma[trii_gamma_coord[0],1] * gamma[0] + vn_mesh_gamma[trii_gamma_coord[1],0] * beta[0] + vn_mesh_gamma[trii_gamma_coord[1],1] * gamma[0],
                                                                            2*(vn_mesh_gamma[trii_gamma_coord[2],0] * beta[1] + vn_mesh_gamma[trii_gamma_coord[2],1] * gamma[1]) + vn_mesh_gamma[trii_gamma_coord[0],0] * beta[1] + vn_mesh_gamma[trii_gamma_coord[0],1] * gamma[1] + vn_mesh_gamma[trii_gamma_coord[1],0] * beta[1] + vn_mesh_gamma[trii_gamma_coord[1],1] * gamma[1],
                                                                            2*(vn_mesh_gamma[trii_gamma_coord[2],0] * beta[2] + vn_mesh_gamma[trii_gamma_coord[2],1] * gamma[2]) + vn_mesh_gamma[trii_gamma_coord[0],0] * beta[2] + vn_mesh_gamma[trii_gamma_coord[0],1] * gamma[2] + vn_mesh_gamma[trii_gamma_coord[1],0] * beta[2] + vn_mesh_gamma[trii_gamma_coord[1],1] * gamma[2]])
            #Mass matrix
            M_V[9*(l-1):9*(l-1)+9]      = np.abs(Delta) / 24 *  (np.ones((3,3)) + np.eye(3)).reshape(9)


            l = l+1

        Diff_gamma = coo_matrix((Diff_V,(I,J)), shape = (grid.diffusion_mesh_gamma.n,grid.diffusion_mesh_gamma.n))
        Diff_gamma = Diff_gamma.tocsr()

        Conv_gamma = coo_matrix((Conv_V,(I,J)), shape = (grid.diffusion_mesh_gamma.n,grid.diffusion_mesh_gamma.n))
        Conv_gamma = Conv_gamma.tocsr()

        M_gamma = coo_matrix((M_V,(I,J)), shape = (grid.diffusion_mesh_gamma.n,grid.diffusion_mesh_gamma.n))
        M_gamma = M_gamma.tocsr()
        M_L_gamma = diags(np.asarray(csr_matrix.sum(M_gamma, axis = 1)).ravel(), format = 'csr')

        Mcexn = M_gamma @ cexn_gamma
        M_Lcexn = M_L_gamma @ cexn_gamma

        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)
        R_V = np.zeros(4*(grid.interface_gammatheta.n_edges))
        f_R_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]])]
            R_V[4*(l-1):4*(l-1)+4] = L/6*np.array([2, 1, 1, 2])
            l = l + 1
            f_R_gt[grid.diffusion_mesh_gamma.enriched_to_gamma[el]] = f_R_gt[grid.diffusion_mesh_gamma.enriched_to_gamma[el]] + L/2 * np.ones(2)
        R_gt = coo_matrix((R_V,(I_gt,J_gt)), shape = (grid.diffusion_mesh_gamma.n,grid.diffusion_mesh_gamma.n))
        R_gt = R_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)
        R_V = np.zeros(4*(grid.interface_gammaalpha.n_edges))
        f_R_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]])
            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]])]
            R_V[4*(l-1):4*(l-1)+4] = L/6*np.array([2, 1, 1, 2])
            l = l + 1
            f_R_ga[grid.diffusion_mesh_gamma.enriched_to_gamma[el]] = f_R_ga[grid.diffusion_mesh_gamma.enriched_to_gamma[el]] + L/2 * np.ones(2)
        R_ga = coo_matrix((R_V,(I_ga,J_ga)), shape = (grid.diffusion_mesh_gamma.n,grid.diffusion_mesh_gamma.n))
        R_ga = R_ga.tocsr()

        f_gamma = np.zeros(grid.diffusion_mesh_gamma.n)
        if not(Curvature):
            f_gamma[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, 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]
                f_gamma[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_gt[x_i_gt]) - param.c_alphagamma_GT(T, grid.enriched_mesh.curvature_ga[x_i_ga])))
                
        ### Flux Correction
        i = 0
        neighbours = lil_matrix((grid.diffusion_mesh_gamma.n, grid.diffusion_mesh_gamma.n), dtype = bool)
        for j in grid.diffusion_mesh_gamma.points_index:
            other_ii_all = np.unique(grid.enriched_mesh.triangles_points_index[grid.enriched_mesh.point_in_tri[j].nonzero()[1]].flatten())
            other_ii_gamma = np.zeros(0,dtype = np.int)
            for ii in other_ii_all:
                #gamma
                if grid.level_set_gammaalpha_enriched.phi[ii] <= 0 and grid.level_set_alphatheta_enriched.phi[ii] <= 0:
                    other_ii_gamma = np.append(other_ii_gamma, ii)
            other_ii_gamma = grid.diffusion_mesh_gamma.enriched_to_gamma[other_ii_gamma]
            neighbours[i, other_ii_gamma] = True
            i = i+1

        #Antidifussive part of diffusion matrix
        Diff_plus_gamma = lil_matrix((grid.diffusion_mesh_gamma.n, grid.diffusion_mesh_gamma.n))
        where_neg = find(Diff_gamma < 0)
        Diff_plus_gamma[where_neg[0], where_neg[1]] = Diff_gamma[where_neg[0], where_neg[1]]
        Diff_plus_gamma.setdiag(0)
        Diff_plus_gamma.setdiag(-lil_matrix.sum(Diff_plus_gamma, axis = 1))
        Diff_plus_gamma = Diff_plus_gamma.tocsr()
        Diff_tilde_gamma = Diff_gamma - Diff_plus_gamma

        #Artificial diffusion matrix D to limit Conv_gamma
        D = lil_matrix((grid.diffusion_mesh_gamma.n, grid.diffusion_mesh_gamma.n))
        for i in np.unique(Conv_gamma.indices):
            where_Conv_neg = np.setdiff1d(find(Conv_gamma[i] < 0)[1], i)
            for j_neg in where_Conv_neg:
                min_neg = max(-Conv_gamma[i,j_neg], 0, -Conv_gamma[j_neg,i])
                D[i,j_neg] = min_neg
                D[j_neg,i] = min_neg
        D.setdiag(-lil_matrix.sum(D, axis = 1))
        D = D.tocsr()
        Conv_tilde_gamma = Conv_gamma + D

        dcdn_alpha_matrix[ga_points, ga_points] = dcdn_alpha_ga * f_R_ga[ga_points]
        dcdn_alpha_matrix = dcdn_alpha_matrix.tocsr()

        r = 1
        eps = 1e6
        FPT = 1e-6
        qr = cexn_gamma.copy()
        qr_tilde = cexn_gamma.copy()
        qr_tilde[ga_points] = param.c_gammaalpha(T)
        qr_tilde[x_points_ga] = (param.c_gammatheta(T) - param.K_alphatheta(T) / param.K_gammatheta(T) * (param.c_alphatheta(T) - param.c_alphagamma(T)))
        qr_H_gamma = cexn_gamma.copy()
        qr_H_gamma[ga_points] = param.c_gammaalpha(T)
        qr_H_gamma[x_points_ga] = (param.c_gammatheta(T) - param.K_alphatheta(T) / param.K_gammatheta(T) * (param.c_alphatheta(T) - param.c_alphagamma(T)))

        dqrdt_tilde = np.ones(grid.diffusion_mesh_gamma.n)
        dqrdt_tilde = np.zeros(grid.diffusion_mesh_gamma.n)
        M_gamma = M_gamma.tocoo()
        while eps >= FPT and r < 21:
            R_cp_V = np.zeros(4*(grid.interface_gammatheta.n_edges))
            f_R_cp_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]
#                ### Reaction not lumped
#                R_cp_V[4*(l-1):4*(l-1)+4] = L/12*np.array([3*qr_tilde[el_g[0]] + qr_tilde[el_g[1]], qr_tilde[el_g[0]] +   qr_tilde[el_g[1]],
#                                                             qr_tilde[el_g[0]] + qr_tilde[el_g[1]], qr_tilde[el_g[0]] + 3*qr_tilde[el_g[1]]])
#                f_R_cp_gt[el_g] = f_R_cp_gt[el_g] + L/6 * np.array([2*qr_tilde[el_g[0]] +   qr_tilde[el_g[1]],
#                                                                      qr_tilde[el_g[0]] + 2*qr_tilde[el_g[1]] ])
                ### Reaction lumped
                R_cp_V[4*(l-1):4*(l-1)+4] = L/12*np.array([3*qr_tilde[el_g[0]] + qr_tilde[el_g[0]], qr_tilde[el_g[0]] +   qr_tilde[el_g[0]],
                                                             qr_tilde[el_g[1]] + qr_tilde[el_g[1]], qr_tilde[el_g[1]] + 3*qr_tilde[el_g[1]]])
                f_R_cp_gt[el_g] = f_R_cp_gt[el_g] + L/2 * np.array([qr_tilde[el_g[0]],
                                                                    qr_tilde[el_g[1]] ])
                l = l + 1
    
                
            R_cp_gt = coo_matrix((R_cp_V,(I_gt,J_gt)), shape = (grid.diffusion_mesh_gamma.n, grid.diffusion_mesh_gamma.n))
            R_cp_gt = R_cp_gt.tocsr()

            R_cp_V = np.zeros(4*(grid.interface_gammaalpha.n_edges))
            f_R_cp_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]
#                ### Reaction not lumped
#                R_cp_V[4*(l-1):4*(l-1)+4] = L/12*np.array([3*qr_tilde[el_g[0]] + qr_tilde[el_g[1]], qr_tilde[el_g[0]] +   qr_tilde[el_g[1]],
#                                                             qr_tilde[el_g[0]] + qr_tilde[el_g[1]], qr_tilde[el_g[0]] + 3*qr_tilde[el_g[1]]])
#                f_R_cp_ga[el_g] = f_R_cp_ga[el_g] + L/6 * np.array([2*qr_tilde[el_g[0]] +   qr_tilde[el_g[1]],
#                                                                      qr_tilde[el_g[0]] + 2*qr_tilde[el_g[1]] ])
                ### Reaction lumped
                R_cp_V[4*(l-1):4*(l-1)+4] = L/12*np.array([3*qr_tilde[el_g[0]] + qr_tilde[el_g[0]], qr_tilde[el_g[0]] +   qr_tilde[el_g[0]],
                                                             qr_tilde[el_g[1]] + qr_tilde[el_g[1]], qr_tilde[el_g[1]] + 3*qr_tilde[el_g[1]]])
                f_R_cp_ga[el_g] = f_R_cp_ga[el_g] + L/2 * np.array([qr_tilde[el_g[0]],
                                                                    qr_tilde[el_g[1]] ])
    
                l = l + 1
    
                
            R_cp_ga = coo_matrix((R_cp_V,(I_ga,J_ga)), shape = (grid.diffusion_mesh_gamma.n,grid.diffusion_mesh_gamma.n))
            R_cp_ga = R_cp_ga.tocsr()

            #Reaction matrix BC
            if not(Curvature):
                React_gamma_ga      = (param.D_alpha(T) * dcdn_alpha_matrix[ga_points] - param.K_gammaalpha(T) * (param.c_alphagamma(T) * R_ga[ga_points] - R_cp_ga[ga_points]))/param.c_alphagamma(T)
            else:
                React_gamma_ga = lil_matrix((grid.interface_gammaalpha.n, grid.diffusion_mesh_gamma.n))
                for i in range(0, grid.interface_gammaalpha.n):
                    React_gamma_ga[i] = (param.D_alpha(T) * dcdn_alpha_matrix[ga_points[i]] - param.K_gammaalpha(T) * (param.c_alphagamma_GT(T, grid.enriched_mesh.curvature_ga[i]) * R_ga[ga_points[i]] - R_cp_ga[ga_points[i]]))/param.c_alphagamma_GT(T, grid.enriched_mesh.curvature_ga[i])
            React_gamma_ga_matrix = lil_matrix((grid.diffusion_mesh_gamma.n, grid.diffusion_mesh_gamma.n))
            React_gamma_ga_matrix[ga_points] = React_gamma_ga
            React_gamma_ga_matrix = React_gamma_ga_matrix.tocsr()

            React_gamma_gt      = - param.K_gammatheta(T) * (param.c0_theta() * R_gt[gt_points] - R_cp_gt[gt_points]) / param.c0_theta()
            React_gamma_gt_matrix = lil_matrix((grid.diffusion_mesh_gamma.n, grid.diffusion_mesh_gamma.n))
            React_gamma_gt_matrix[gt_points] = React_gamma_gt
            React_gamma_gt_matrix = React_gamma_gt_matrix.tocsr()

            React_gamma_matrix = React_gamma_gt_matrix + React_gamma_ga_matrix
            React_gamma_matrix_lumped = diags(np.asarray(csr_matrix.sum(React_gamma_matrix, axis = 1)).ravel(), format = 'csr')    

            if not(Curvature):
                f_gamma[ga_points]   = - param.K_gammaalpha(T) * param.c_gammaalpha(T) * (param.c_alphagamma(T) * f_R_ga[ga_points] - f_R_cp_ga[ga_points]) / param.c_alphagamma(T)
                f_gamma[gt_points]   = - param.K_gammatheta(T) * param.c_gammatheta(T) * (param.c0_theta() * f_R_gt[gt_points] - f_R_cp_gt[gt_points]) / param.c0_theta()
                f_gamma[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_gamma[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_R_ga[ga_points[i]] - f_R_cp_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_gamma[gt_points[i]]   = - param.K_gammatheta(T) * param.c_gammatheta_GT(T, grid.enriched_mesh.curvature_gt[i]) * (param.c0_theta() * f_R_gt[gt_points[i]] - f_R_cp_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]
                    f_gamma[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_gt[x_i_gt]) - param.c_alphagamma_GT(T, grid.enriched_mesh.curvature_ga[x_i_ga])))



            B_tilde_gamma = (M_L_gamma - dt*(Diff_tilde_gamma + Conv_tilde_gamma + React_gamma_matrix_lumped)).tocsr() #React_gamma_matrix_lumped
            ### Dirichlet
            [setcsrrow2id(B_tilde_gamma, x_points_ga[i]) for i in range(0, len(x_points_ga))]

            g_gamma = M_Lcexn - dt*f_gamma
            if not(Curvature):
                g_gamma[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, 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]
                    g_gamma[x_points_ga] = 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_gt[x_i_gt]) - param.c_alphagamma_GT(T, grid.enriched_mesh.curvature_ga[x_i_ga]))

            eps_qr_tilde = la.norm(B_tilde_gamma @ qr_tilde - g_gamma, np.inf) / la.norm(g_gamma, np.inf)

            qr_tilde = spla.spsolve(B_tilde_gamma, g_gamma)
            #Approximate d/dt u_i, easy way
            Su = (Conv_tilde_gamma + Diff_tilde_gamma + React_gamma_matrix_lumped) @ qr_tilde - f_gamma
            Su[x_points_ga] = 0
            eps_dqrdt_tilde = la.norm(M_L_gamma @ dqrdt_tilde - Su, np.inf) / la.norm(Su, np.inf)

            dqrdt_tilde = spla.spsolve(M_L_gamma, Su)
    #        #Approximate d/dt u_i, less diffusive way
    #        Su = (Conv_gamma + Diff_gamma + React_gamma_matrix_lumped) @ qr_tilde - f_gamma
    #        Su[x_points_ga] = 0
    #        dqrdt_tilde = spla.spsolve(M_gamma.tocsr(), Su)

            R_cp_V = np.zeros(4*(grid.interface_gammatheta.n_edges))
            f_R_cp_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]
#                ### Reaction not lumped                
#                R_cp_V[4*(l-1):4*(l-1)+4] = L/12*np.array([3*qr_H_gamma[el_g[0]] + qr_H_gamma[el_g[1]], qr_H_gamma[el_g[0]] +   qr_H_gamma[el_g[1]],
#                                                             qr_H_gamma[el_g[0]] + qr_H_gamma[el_g[1]], qr_H_gamma[el_g[0]] + 3*qr_H_gamma[el_g[1]]])
#                f_R_cp_gt[el_g] = f_R_cp_gt[el_g] + L/6 * np.array([2*qr_H_gamma[el_g[0]] +   qr_H_gamma[el_g[1]],
#                                                                      qr_H_gamma[el_g[0]] + 2*qr_H_gamma[el_g[1]] ])
                ### Reaction lumped
                R_cp_V[4*(l-1):4*(l-1)+4] = L/12*np.array([3*qr_H_gamma[el_g[0]] + qr_H_gamma[el_g[0]], qr_H_gamma[el_g[0]] +   qr_H_gamma[el_g[0]],
                                                             qr_H_gamma[el_g[1]] + qr_H_gamma[el_g[1]], qr_H_gamma[el_g[1]] + 3*qr_H_gamma[el_g[1]]])
                l = l + 1
                
                f_R_cp_gt[el_g] = f_R_cp_gt[el_g] + L/2 * np.array([qr_H_gamma[el_g[0]],
                                                                    qr_H_gamma[el_g[1]] ])
            R_cp_gt = coo_matrix((R_cp_V,(I_gt,J_gt)), shape = (grid.diffusion_mesh_gamma.n, grid.diffusion_mesh_gamma.n))
            R_cp_gt = R_cp_gt.tocsr()

            R_cp_V = np.zeros(4*(grid.interface_gammaalpha.n_edges))
            f_R_cp_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]
#                ### Reaction not lumped                
#                R_cp_V[4*(l-1):4*(l-1)+4] = L/12*np.array([3*qr_H_gamma[el_g[0]] + qr_H_gamma[el_g[1]], qr_H_gamma[el_g[0]] +   qr_H_gamma[el_g[1]],
#                                                             qr_H_gamma[el_g[0]] + qr_H_gamma[el_g[1]], qr_H_gamma[el_g[0]] + 3*qr_H_gamma[el_g[1]]])
#                f_R_cp_ga[el_g] = f_R_cp_ga[el_g] + L/6 * np.array([2*qr_H_gamma[el_g[0]] +   qr_H_gamma[el_g[1]],
#                                                                      qr_H_gamma[el_g[0]] + 2*qr_H_gamma[el_g[1]] ])
                ### Reaction lumped
                R_cp_V[4*(l-1):4*(l-1)+4] = L/12*np.array([3*qr_H_gamma[el_g[0]] + qr_H_gamma[el_g[0]], qr_H_gamma[el_g[0]] +   qr_H_gamma[el_g[0]],
                                                             qr_H_gamma[el_g[1]] + qr_H_gamma[el_g[1]], qr_H_gamma[el_g[1]] + 3*qr_H_gamma[el_g[1]]])
                f_R_cp_ga[el_g] = f_R_cp_ga[el_g] + L/2 * np.array([qr_H_gamma[el_g[0]],
                                                                    qr_H_gamma[el_g[1]] ])
                l = l + 1
    
                
            R_cp_ga = coo_matrix((R_cp_V,(I_ga,J_ga)), shape = (grid.diffusion_mesh_gamma.n,grid.diffusion_mesh_gamma.n))
            R_cp_ga = R_cp_ga.tocsr()

            #Reaction matrix BC
            if not(Curvature):
                React_gamma_ga = (param.D_alpha(T) * dcdn_alpha_matrix[ga_points] - param.K_gammaalpha(T) * (param.c_alphagamma(T) * R_ga[ga_points] - R_cp_ga[ga_points]))/param.c_alphagamma(T)
            else:
                React_gamma_ga = lil_matrix((grid.interface_gammaalpha.n, grid.diffusion_mesh_gamma.n))
                for i in np.arange(0, grid.interface_gammaalpha.n):
                    React_gamma_ga[i] = (param.D_alpha(T) * dcdn_alpha_matrix[ga_points[i]] - param.K_gammaalpha(T) * (param.c_alphagamma_GT(T, grid.enriched_mesh.curvature_ga[i]) * R_ga[ga_points[i]] - R_cp_ga[ga_points[i]]))/param.c_alphagamma_GT(T, grid.enriched_mesh.curvature_ga[i])
            React_gamma_ga_matrix = lil_matrix((grid.diffusion_mesh_gamma.n, grid.diffusion_mesh_gamma.n))
            React_gamma_ga_matrix[ga_points] = React_gamma_ga
            React_gamma_ga_matrix = React_gamma_ga_matrix.tocsr()


            React_gamma_gt = - param.K_gammatheta(T) * (param.c0_theta() * R_gt[gt_points] - R_cp_gt[gt_points]) / param.c0_theta()
            React_gamma_gt_matrix = lil_matrix((grid.diffusion_mesh_gamma.n, grid.diffusion_mesh_gamma.n))
            React_gamma_gt_matrix[gt_points] = React_gamma_gt
            React_gamma_gt_matrix = React_gamma_gt_matrix.tocsr()

            React_gamma_matrix = React_gamma_gt_matrix + React_gamma_ga_matrix
            React_gamma_matrix_lumped = diags(np.asarray(csr_matrix.sum(React_gamma_matrix, axis = 1)).ravel(), format = 'csr')

            if not(Curvature):
                f_gamma[ga_points]   = - param.K_gammaalpha(T) * param.c_gammaalpha(T) * (param.c_alphagamma(T) * f_R_ga[ga_points] - f_R_cp_ga[ga_points]) / param.c_alphagamma(T)
                f_gamma[gt_points]   = - param.K_gammatheta(T) * param.c_gammatheta(T) * (param.c0_theta() * f_R_gt[gt_points] - f_R_cp_gt[gt_points]) / param.c0_theta()
                f_gamma[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_gamma[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_R_ga[ga_points[i]] - f_R_cp_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_gamma[gt_points[i]]   = - param.K_gammatheta(T) * param.c_gammatheta_GT(T, grid.enriched_mesh.curvature_gt[i]) * (param.c0_theta() * f_R_gt[gt_points[i]] - f_R_cp_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]
                    f_gamma[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_gt[x_i_gt]) - param.c_alphagamma_GT(T, grid.enriched_mesh.curvature_ga[x_i_ga])))


            B_gamma = (M_gamma.tocsr() - dt*(Diff_gamma + Conv_gamma + React_gamma_matrix_lumped)).tocsr() #React_gamma_matrix_lumped
            [setcsrrow2id(B_gamma, x_points_ga[i]) for i in range(0, len(x_points_ga))]
            g_gamma = Mcexn - dt*f_gamma
            if not(Curvature):
                g_gamma[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, 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]
                    g_gamma[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_gt[x_i_gt]) - param.c_alphagamma_GT(T, grid.enriched_mesh.curvature_ga[x_i_ga]))

            eps_qr_H = la.norm(B_gamma @ qr_H_gamma - g_gamma, np.inf) / la.norm(g_gamma, np.inf)

            qr_H_gamma = spla.spsolve(B_gamma, g_gamma)


            eps = max(eps_qr_tilde, eps_dqrdt_tilde, eps_qr_H)
            eps = eps_qr_H
            print('eps gamma %.3e' % eps, 'qr_tilde, dqrtildedt %.3e %.3e' %(eps_qr_tilde ,eps_dqrdt_tilde))
            r=r+1
        plot_c_gamma = False
        if r > 20 and eps > FPT:
            print('Stopped gamma reaction iteration')

        ### Flux Reaction
        #flux f_ij = (mij d/dt + d_ij + l+_ij)(u_i-u_j) c^n+1 approx bij qr_tilde, l+ should be zero
        flux_M_gamma_V = np.zeros(M_gamma.nnz)
        for ij in range(0, M_gamma.nnz):
            i = M_gamma.row[ij]
            j = M_gamma.col[ij]
            flux_M_gamma_V[ij] = M_gamma.data[ij] * (dqrdt_tilde[i] - dqrdt_tilde[j])
        DDiff_plus_gamma = (D-Diff_plus_gamma).tocoo()
        flux_DDiff_plus_gamma_V = np.zeros(DDiff_plus_gamma.nnz)
        for ij in range(0, DDiff_plus_gamma.nnz):
            i = DDiff_plus_gamma.row[ij]
            j = DDiff_plus_gamma.col[ij]
            flux_DDiff_plus_gamma_V[ij] = DDiff_plus_gamma.data[ij] * (qr_tilde[i] - qr_tilde[j])
        flux = (coo_matrix((flux_M_gamma_V,(M_gamma.row,M_gamma.col)), shape = (grid.diffusion_mesh_gamma.n,grid.diffusion_mesh_gamma.n)).tocsr() +
                coo_matrix((flux_DDiff_plus_gamma_V,(DDiff_plus_gamma.row,DDiff_plus_gamma.col)), shape = (grid.diffusion_mesh_gamma.n,grid.diffusion_mesh_gamma.n)).tocsr()).tolil()

        flux.setdiag(0)
        flux.setdiag(lil_matrix.sum(flux, axis = 1))
        flux = flux.tocoo()


        for ij in range(0, flux.nnz):
            i = flux.row[ij]
            j = flux.col[ij]
            if flux.data[ij] * (qr_tilde[j] - qr_tilde[i]) > 0:
                flux.data[ij] = 0
        flux = flux.tocsr()




        qr_max = -np.ones(grid.diffusion_mesh_gamma.n)
        qr_min = -np.ones(grid.diffusion_mesh_gamma.n)
        for i in range(0, grid.diffusion_mesh_gamma.n):
            qr_max[i] = np.max(qr_tilde[np.asarray(neighbours[i].todense()).reshape(grid.diffusion_mesh_gamma.n)])
            qr_min[i] = np.min(qr_tilde[np.asarray(neighbours[i].todense()).reshape(grid.diffusion_mesh_gamma.n)])

        if not(Curvature):
            qr_min[gt_points] = np.maximum(qr_min[gt_points], param.c_gammatheta(T))
            qr_min[ga_points] = np.maximum(qr_min[ga_points], param.c_gammaalpha(T))
            qr_min[x_points_ga] = np.maximum(qr_min[x_points_ga], np.maximum(param.c_gammatheta(T), param.c_gammaalpha(T)))
        else:
            for i in range(0, grid.interface_gammaalpha.n):
                qr_min[ga_points[i]] = np.maximum(qr_min[ga_points[i]], param.c_gammaalpha_GT(T, grid.enriched_mesh.curvature_ga[i]))
            for i in range(0, grid.interface_gammatheta.n):
                qr_min[gt_points[i]] = np.maximum(qr_min[gt_points[i]], param.c_gammatheta_GT(T, grid.enriched_mesh.curvature_gt[i]))
            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]
                qr_min[x_points_ga[i]] = np.maximum(qr_min[x_points_ga[i]], np.maximum(param.c_gammatheta_GT(T, grid.enriched_mesh.curvature_gt[x_i_gt]), param.c_gammaalpha_GT(T, grid.enriched_mesh.curvature_ga[x_i_ga])) )

        R_min = np.ones(grid.diffusion_mesh_gamma.n)
        R_plus = np.ones(grid.diffusion_mesh_gamma.n)
        nonzero = flux.nonzero()
        nonzero_i = np.setdiff1d(np.unique(nonzero[0]), x_points_ga)
        for i in nonzero_i:
            row_flux_pos = (flux[i] + np.abs(flux[i]))/2
            row_flux_neg = (flux[i] - np.abs(flux[i]))/2
            P_plus = csr_matrix.sum(row_flux_pos) - row_flux_pos[0,i]
            P_min  = csr_matrix.sum(row_flux_neg) - row_flux_neg[0,i]

            Q_plus = M_L_gamma[i,i] / dt * (qr_max[i] - qr_tilde[i])
            Q_min  = M_L_gamma[i,i] / dt * (qr_min[i] - qr_tilde[i])

            if P_plus != 0:
                R_plus[i] = np.min([1, Q_plus/P_plus])
            if P_min != 0:
                R_min[i]  = np.min([1, Q_min/P_min])


        coeff = lil_matrix((grid.diffusion_mesh_gamma.n, grid.diffusion_mesh_gamma.n))
        for k in range(0,len(nonzero[0])):
            i = nonzero[0][k]
            j = nonzero[1][k]
            if flux[i,j] > 0:
                coeff[i,j] = np.min([R_plus[i], R_min[j]])
            elif flux[i,j] < 0:
                coeff[i,j] = np.min([R_min[i], R_plus[j]])

        coeff = coeff.tocsr()
        fbar = coeff.multiply(flux)
        fbar = fbar.tolil()
        fbar.setdiag(0)
        fbar = lil_matrix.sum(fbar, axis = 1)
        fbar[x_points_ga] = 0
        fbar = np.squeeze(np.asarray(fbar))

        qr = spla.spsolve(M_L_gamma, M_L_gamma @ qr_tilde + dt * fbar)

        c_gamma = qr.copy()
    return c_gamma, c_alpha, dcdn_alpha_ga ,t






def calc_v_vn(grid, c_gamma, c_alpha, dcdn_alpha_ga, c0_theta, T_old, Curvature):
    #Calc dc_alpha/dn
    ga_points = grid.interface_gammaalpha.points_index
    ga_points_gamma_coord = grid.diffusion_mesh_gamma.enriched_to_gamma[ga_points]
    gt_points = grid.interface_gammatheta.points_index
    gt_points_gamma_coord = grid.diffusion_mesh_gamma.enriched_to_gamma[gt_points]
    at_points = grid.interface_alphatheta.points_index #in backgr version no triple points for at (already in gt)
    at_points_alpha_coord = grid.diffusion_mesh_alpha.enriched_to_alpha[at_points]

    if not(Curvature):
        vn_gammaalpha = -param.D_alpha(T_old)/param.c_alphagamma(T_old) * dcdn_alpha_ga + param.K_gammaalpha(T_old)/param.c_alphagamma(T_old)*(param.c_gammaalpha(T_old) - c_gamma[ga_points_gamma_coord])
        vn_gammatheta = (param.K_gammatheta(T_old) * (param.c_gammatheta(T_old) - c_gamma[gt_points_gamma_coord])) /param.c0_theta()
        vn_alphatheta = (param.K_alphatheta(T_old) * (param.c_alphatheta(T_old) - c_alpha[at_points_alpha_coord])) /param.c0_theta()
    else:
        vn_gammaalpha = np.zeros(grid.interface_gammaalpha.n)
        vn_gammatheta = np.zeros(grid.interface_gammatheta.n)
        vn_alphatheta = np.zeros(grid.interface_alphatheta.n)
        for i in range(0, grid.interface_gammaalpha.n):
            vn_gammaalpha[i] = -param.D_alpha(T_old)/param.c_alphagamma_GT(T_old, grid.enriched_mesh.curvature_ga[i]) * dcdn_alpha_ga[i] + param.K_gammaalpha(T_old)/param.c_alphagamma(T_old)*(param.c_gammaalpha_GT(T_old, grid.enriched_mesh.curvature_ga[i]) - c_gamma[ga_points_gamma_coord[i]])
        for i in range(0, grid.interface_gammatheta.n):
            vn_gammatheta[i] = (param.K_gammatheta(T_old) * (param.c_gammatheta_GT(T_old, grid.enriched_mesh.curvature_gt[i]) - c_gamma[gt_points_gamma_coord[i]])) /param.c0_theta()
        for i in range(0, grid.interface_alphatheta.n):
            vn_alphatheta[i] = (param.K_alphatheta(T_old) * (param.c_alphatheta_GT(T_old, grid.enriched_mesh.curvature_at[i]) - c_alpha[at_points_alpha_coord[i]])) /param.c0_theta()

    return vn_gammaalpha, vn_gammatheta, vn_alphatheta#, v_gammaalpha, v_gammatheta, v_alphatheta


def Laplace_vn(grid, vn_ga, vn_t):
    #vn_t should be given stacked with gt first, then at
    ga_points = grid.interface_gammaalpha.points_index
    ga_ghost_points = grid.interface_gammaalpha.ghost_points_index
    at_points = grid.interface_alphatheta.points_index
    gt_points = grid.interface_gammatheta.points_index
    t_points  = np.union1d(gt_points, at_points)
    I = np.zeros(9*grid.enriched_mesh.n_tri, dtype = np.int)
    J = np.zeros(9*grid.enriched_mesh.n_tri, dtype = np.int)
    A_V = np.zeros(9*grid.enriched_mesh.n_tri)
    l = 1
    for i in range(0,grid.enriched_mesh.n_tri):
        trii = grid.enriched_mesh.triangles_points_index[i]
        tri = grid.enriched_mesh.triangles[i]
        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[0], np.ones(3, dtype = np.int) * trii[1], np.ones(3, dtype = np.int) * trii[2]))
        J[9*(l-1):9*(l-1)+9] = np.hstack((trii, trii, trii))



        A_V[9*(l-1):9*(l-1)+9] = - 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 = coo_matrix((A_V,(I,J)), shape = (grid.enriched_mesh.n,grid.enriched_mesh.n))
    A = A.tocsr()


    A_ga = A.copy()

    [setcsrrow2id(A_ga, ga_points[i]) for i in range(0, len(ga_points))]
    f_ga = np.zeros(grid.enriched_mesh.n)
    f_ga[ga_points] = vn_ga

#    ###ghost extension only Dirichlet on t1/2 point vn = 0
#    setcsrrow2id(A_ga, grid.interface_gammaalpha.ghost_points_index[grid.interface_gammaalpha.ghost_index_max_d_triple])
#    f_ga[grid.interface_gammaalpha.ghost_points_index[grid.interface_gammaalpha.ghost_index_max_d_triple]] = 0

#     ###ghost extension over whole ga_ghost, with max dist/distance
#    [setcsrrow2id(A_ga, ga_ghost_points[i]) for i in range(0, len(ga_ghost_points))]
#
#    for i in range(0,len(ga_ghost_points)):
#        L_x = [la.norm(grid.interface_gammaalpha.ghost_points[i] - grid.enriched_mesh.points[grid.triple_points_index[j]], 2) for j in range(0, len(grid.triple_points_index))]
#        L = np.max(L_x)
#        where_L = np.argmax(L_x)
#        f_ga[ga_ghost_points[i]] = (1 - L/grid.interface_gammaalpha.ghost_max_d_triple) * vn_ga[grid.interface_gammaalpha.triple_points_index[where_L]]

#    ###ghost extension over whole ga_ghost with fitted vn_ga by paramaterisation with v = 0 on max distance ghost point
#    [setcsrrow2id(A_ga, ga_ghost_points[i]) for i in range(0, len(ga_ghost_points))]
#    t_12 = grid.interface_gammaalpha.parametrisation_ghost[grid.interface_gammaalpha.ghost_index_max_d_triple]
#    vn_ga_ghost_fit = la.solve(np.array([[1, 0, 0, 0, 0, 0],
#                                         [1, 1, 1, 1, 1, 1],
#                                         [1, t_12, t_12**2, t_12**3, t_12**4, t_12**5],
#                                         [0, 1, 0, 0, 0, 0],
#                                         [0, 1, 2, 3, 4, 5],
#                                         [0, 1, 2*t_12, 3*t_12**2, 4*t_12**3, 5*t_12**4]]),
#                                np.array([vn_ga[grid.interface_gammaalpha.triple_points_index[0]], #v_t0
#                                          vn_ga[grid.interface_gammaalpha.triple_points_index[1]], #v_t1
#                                          0, #v_t_12
#                                          0, #dv/dt|x0
#                                          0, #dv/dt|x0
#                                          0])) #dv/dt|t_12
#
#
#    vn_ga_ghost_parametrisation = np.sum(vn_ga_ghost_fit * np.vstack((grid.interface_gammaalpha.parametrisation_ghost**0, grid.interface_gammaalpha.parametrisation_ghost**1, grid.interface_gammaalpha.parametrisation_ghost**2, grid.interface_gammaalpha.parametrisation_ghost**3, grid.interface_gammaalpha.parametrisation_ghost**4, grid.interface_gammaalpha.parametrisation_ghost**5)).T, axis = 1)
#    f_ga[grid.interface_gammaalpha.ghost_points_index] = vn_ga_ghost_parametrisation

    ### Check if there is ghost interface and if 1 or 2 particles and if there is fixed point
    if grid.interface_gammatheta.n > 0 and grid.interface_gammaalpha.ghost_edges_index.shape[0] > 0 and grid.two_particles == None and grid.background_mesh.i_ghost_fixed != -1:
        ###ghost extension over whole ga_ghost with fitted vn_ga by paramaterisation with v = 0 on fixed point
        [setcsrrow2id(A_ga, ga_ghost_points[i]) for i in range(0, len(ga_ghost_points))]

#        edges = grid.enriched_mesh.points[grid.interface_gammaalpha.ghost_edges_index]
#        l_ga_g = np.sum(la.norm(edges[:,0]-edges[:,1],2,axis = 1))
        l_ga_g = grid.interface_gammaalpha.parametrisation_ghost[0,grid.interface_gammaalpha.ghost_triple_points_index[1]]

#        #fit 2nd order polynomial to vn_ga in order to approximate dvn_ga/dt in triple points
#        x = np.zeros((3,2))
#        t = np.zeros(3)
#        xii = -np.ones(3, dtype = np.int)
#        xii_ga = -np.ones(3, dtype = np.int)
#        x_i = grid.triple_points_index[0]
#        x[0] = grid.enriched_mesh.points[x_i]
#        xii[0] = x_i
#        xii_ga[0], = np.where(grid.interface_gammaalpha.points_index == x_i)
#        h_sum = 0
#        for i in range(1,3):
#            patch_x_i = np.setdiff1d(grid.enriched_mesh.triangles_points_index[grid.enriched_mesh.point_in_tri[x_i].nonzero()[1]].flatten(), xii)
#            where_ga, = np.where(grid.has_point_been_shifted[patch_x_i,1] == 1)
#            x_i = patch_x_i[where_ga[0]]
#            x[i] = grid.enriched_mesh.points[x_i]
#            xii[i] = x_i
#            xii_ga[i], = np.where(grid.interface_gammaalpha.points_index == x_i)
#            h_sum = h_sum - la.norm(x[i] - x[i-1], 2)
#            t[i] = h_sum
#        T = np.vstack((np.ones(3), t, t**2)).T
#        dvdt0 = la.solve(T, vn_ga[xii_ga])[1]
##        dvdt0 = 0
#
#        x = np.zeros((3,2))
#        t = np.zeros(3)
#        xii = -np.ones(3, dtype = np.int)
#        xii_ga = -np.ones(3, dtype = np.int)
#        x_i = grid.triple_points_index[1]
#        x[0] = grid.enriched_mesh.points[x_i]
#        xii[0] = x_i
#        xii_ga[0], = np.where(grid.interface_gammaalpha.points_index == x_i)
#        h_sum = l_ga_g
#        t[0] = h_sum
#        for i in range(1,3):
#            patch_x_i = np.setdiff1d(grid.enriched_mesh.triangles_points_index[grid.enriched_mesh.point_in_tri[x_i].nonzero()[1]].flatten(), xii)
#            where_ga, = np.where(grid.has_point_been_shifted[patch_x_i,1] == 1)
#            x_i = patch_x_i[where_ga[0]]
#            x[i] = grid.enriched_mesh.points[x_i]
#            xii[i] = x_i
#            xii_ga[i], = np.where(grid.interface_gammaalpha.points_index == x_i)
#            h_sum = h_sum + la.norm(x[i] - x[i-1], 2)
#            t[i] = h_sum
#        T = np.vstack((np.ones(3), t, t**2)).T
#        solve = la.solve(T, vn_ga[xii_ga])
#        dvdt1 = solve[1] + 2 * solve[2] * t[0]
##        dvdt1 = 0

        ### Fit 1nd order polynomial to vn_ga in order to approximate dvn_ga/dt in triple points. 
        ### First order will ensure the two line segments through the triple point will stay almost a straight line
        x = np.zeros((2,2))
        t = np.zeros(2)
        xii = -np.ones(2, dtype = np.int)
        xii_ga = -np.ones(2, dtype = np.int)
        x_i = grid.triple_points_index[0]
        x[0] = grid.enriched_mesh.points[x_i]
        xii[0] = x_i
        xii_ga[0], = np.where(grid.interface_gammaalpha.points_index == x_i)
        h_sum = 0
        for i in range(1,2):
            patch_x_i = np.setdiff1d(grid.enriched_mesh.triangles_points_index[grid.enriched_mesh.point_in_tri[x_i].nonzero()[1]].flatten(), xii)
            where_ga, = np.where(grid.has_point_been_shifted[patch_x_i,1] == 1)
            x_i = patch_x_i[where_ga[0]]
            x[i] = grid.enriched_mesh.points[x_i]
            xii[i] = x_i
            xii_ga[i], = np.where(grid.interface_gammaalpha.points_index == x_i)
            h_sum = h_sum - la.norm(x[i] - x[i-1], 2)
            t[i] = h_sum
        T = np.vstack((np.ones(2), t)).T
        dvdt0_lin = la.solve(T, vn_ga[xii_ga])[1]

        x = np.zeros((2,2))
        t = np.zeros(2)
        xii = -np.ones(2, dtype = np.int)
        xii_ga = -np.ones(2, dtype = np.int)
        x_i = grid.triple_points_index[1]
        x[0] = grid.enriched_mesh.points[x_i]
        xii[0] = x_i
        xii_ga[0], = np.where(grid.interface_gammaalpha.points_index == x_i)
        h_sum = l_ga_g
        t[0] = h_sum
        for i in range(1,2):
            patch_x_i = np.setdiff1d(grid.enriched_mesh.triangles_points_index[grid.enriched_mesh.point_in_tri[x_i].nonzero()[1]].flatten(), xii)
            where_ga, = np.where(grid.has_point_been_shifted[patch_x_i,1] == 1)
            x_i = patch_x_i[where_ga[0]]
            x[i] = grid.enriched_mesh.points[x_i]
            xii[i] = x_i
            xii_ga[i], = np.where(grid.interface_gammaalpha.points_index == x_i)
            h_sum = h_sum + la.norm(x[i] - x[i-1], 2)
            t[i] = h_sum
        T = np.vstack((np.ones(2), t)).T
        solve = la.solve(T, vn_ga[xii_ga])
        dvdt1_lin = solve[1]

        dvdt0 = dvdt0_lin/10
        dvdt1 = dvdt1_lin/10

        t_12 = grid.background_mesh.t_parametrisation_ghost_fixed[0]
        vn_ga_ghost_fit = la.solve(np.array([[1, 0, 0, 0, 0, 0],
                                             [1, l_ga_g, l_ga_g**2, l_ga_g**3, l_ga_g**4, l_ga_g**5],
                                             [1, t_12  , t_12**2  , t_12**3  , t_12**4  , t_12**5],
                                             [0, 1, 0, 0, 0, 0],
                                             [0, 1, 2*l_ga_g, 3*l_ga_g**2, 4*l_ga_g**3, 5*l_ga_g**4],
                                             [0, 1, 2*t_12, 3*t_12**2, 4*t_12**3, 5*t_12**4]]),
                                    np.array([vn_ga[grid.interface_gammaalpha.triple_points_index[0]], #v_t0
                                              vn_ga[grid.interface_gammaalpha.triple_points_index[1]], #v_t1
                                              0, #v_t_12
                                              dvdt0, #dv/dt|t0
                                              dvdt1, #dv/dt|t1
                                              0])) #dv/dt|t_12


        vn_ga_ghost_parametrisation = np.sum(vn_ga_ghost_fit * np.vstack((grid.interface_gammaalpha.parametrisation_ghost**0, grid.interface_gammaalpha.parametrisation_ghost**1, grid.interface_gammaalpha.parametrisation_ghost**2, grid.interface_gammaalpha.parametrisation_ghost**3, grid.interface_gammaalpha.parametrisation_ghost**4, grid.interface_gammaalpha.parametrisation_ghost**5)).T, axis = 1)
        f_ga[grid.interface_gammaalpha.ghost_points_index] = vn_ga_ghost_parametrisation
    elif grid.two_particles == 1 and grid.interface_gammatheta.n > 0 and grid.interface_gammaalpha.ghost_edges_index.shape[0] > 0:
        ### ghost extension over whole ga_ghost with fitted vn_ga by paramaterisation with v = 0 on fixed point (ellipse[0], ellipse[1] - b_ga)
        [setcsrrow2id(A_ga, ga_ghost_points[i]) for i in range(0, len(ga_ghost_points))]
        for i3 in range(0,int(len(grid.triple_points_index)/2)):
            parametrisation = grid.interface_gammaalpha.parametrisation_ghost[i3]
            where_param, = np.where(grid.interface_gammaalpha.ghost_points_index_which_param == i3)
            parametrisation = parametrisation[where_param]

            l_ga_g = np.max(parametrisation)

            x_triple = np.zeros((2,2))

            #fit 1nd order polynomial to vn_ga in order to approximate dvn_ga/dt in triple points. First order will ensure the two line segments through the triple point will stay almost a straight line
            x = np.zeros((2,2))
            t = np.zeros(2)
            xii = -np.ones(2, dtype = np.int)
            xii_ga = -np.ones(2, dtype = np.int)
            x_i = grid.interface_gammaalpha.ghost_points_index[grid.interface_gammaalpha.parametrisation_triple_points_index[i3,0]]
            x[0] = grid.enriched_mesh.points[x_i]
            x_triple[0] = x[0]
            xii[0] = x_i
            xii_ga[0], = np.where(grid.interface_gammaalpha.points_index == x_i)
            h_sum = 0
            for i in range(1,2):
                patch_x_i = np.setdiff1d(grid.enriched_mesh.triangles_points_index[grid.enriched_mesh.point_in_tri[x_i].nonzero()[1]].flatten(), xii)
                where_ga, = np.where(grid.has_point_been_shifted[patch_x_i,1] == 1)
                x_i = patch_x_i[where_ga[0]]
                x[i] = grid.enriched_mesh.points[x_i]
                xii[i] = x_i
                xii_ga[i], = np.where(grid.interface_gammaalpha.points_index == x_i)
                h_sum = h_sum - la.norm(x[i] - x[i-1], 2)
                t[i] = h_sum
            T = np.vstack((np.ones(2), t)).T
            dvdt0_lin = la.solve(T, vn_ga[xii_ga])[1]

            x = np.zeros((2,2))
            t = np.zeros(2)
            xii = -np.ones(2, dtype = np.int)
            xii_ga = -np.ones(2, dtype = np.int)
            x_i = grid.interface_gammaalpha.ghost_points_index[grid.interface_gammaalpha.parametrisation_triple_points_index[i3,1]]
            x[0] = grid.enriched_mesh.points[x_i]
            x_triple[1] = x[0]
            xii[0] = x_i
            xii_ga[0], = np.where(grid.interface_gammaalpha.points_index == x_i)
            h_sum = l_ga_g
            t[0] = h_sum
            for i in range(1,2):
                patch_x_i = np.setdiff1d(grid.enriched_mesh.triangles_points_index[grid.enriched_mesh.point_in_tri[x_i].nonzero()[1]].flatten(), xii)
                where_ga, = np.where(grid.has_point_been_shifted[patch_x_i,1] == 1)
                x_i = patch_x_i[where_ga[0]]
                x[i] = grid.enriched_mesh.points[x_i]
                xii[i] = x_i
                xii_ga[i], = np.where(grid.interface_gammaalpha.points_index == x_i)
                h_sum = h_sum + la.norm(x[i] - x[i-1], 2)
                t[i] = h_sum
            T = np.vstack((np.ones(2), t)).T
            solve = la.solve(T, vn_ga[xii_ga])
            dvdt1_lin = solve[1]

            dvdt0 = dvdt0_lin/10
            dvdt1 = dvdt1_lin/10

            res = l_ga_g
            for i in where_param:
                if np.abs(grid.interface_gammaalpha.parametrisation_ghost[i3,i] - l_ga_g / 2) < res:
                    t_12 = grid.interface_gammaalpha.parametrisation_ghost[i3, i]
                    res = np.abs(grid.interface_gammaalpha.parametrisation_ghost[i3,i] - l_ga_g / 2)

            x3_0, = np.where(grid.interface_gammaalpha.points_index == grid.interface_gammaalpha.ghost_points_index[grid.interface_gammaalpha.parametrisation_triple_points_index[i3,0]])[0]
            x3_1, = np.where(grid.interface_gammaalpha.points_index == grid.interface_gammaalpha.ghost_points_index[grid.interface_gammaalpha.parametrisation_triple_points_index[i3,1]])[0]

            vn_ga_ghost_fit = la.solve(np.array([[1, 0, 0, 0, 0, 0],
                                                 [1, l_ga_g, l_ga_g**2, l_ga_g**3, l_ga_g**4, l_ga_g**5],
                                                 [1, t_12  , t_12**2  , t_12**3  , t_12**4  , t_12**5],
                                                 [0, 1, 0, 0, 0, 0],
                                                 [0, 1, 2*l_ga_g, 3*l_ga_g**2, 4*l_ga_g**3, 5*l_ga_g**4],
                                                 [0, 1, 2*t_12, 3*t_12**2, 4*t_12**3, 5*t_12**4]]),
                                        np.array([vn_ga[x3_0], #v_t0
                                                  vn_ga[x3_1], #v_t1
                                                  0, #v_t_12
                                                  dvdt0, #dv/dt|t0
                                                  dvdt1, #dv/dt|t1
                                                  0])) #dv/dt|t_12


            vn_ga_ghost_parametrisation = np.sum(vn_ga_ghost_fit * np.vstack((parametrisation**0, parametrisation, parametrisation**2, parametrisation**3, parametrisation**4, parametrisation**5)).T, axis = 1)
            f_ga[grid.interface_gammaalpha.ghost_points_index[where_param]] = vn_ga_ghost_parametrisation

#    f_ga = np.append(f_ga, 0)

#    F_ga = np.append(grid.background_mesh.Tau_fixed_ghost_point, np.zeros(grid.enriched_mesh.n - grid.N))
#    F_ga = np.zeros(grid.enriched_mesh.n)
#    F_ga[grid.background_mesh.i_ghost_fixed] = 1

#    System_ga = bmat([[A_ga, F_ga[:,None]], [F_ga[:,None].T, None]], format = 'csr')

#    vn_ex_ga = spla.spsolve(System_ga, f_ga)[0:grid.enriched_mesh.n]


    vn_ex_ga = spla.spsolve(A_ga, f_ga)
    if grid.background_mesh.i_ghost_fixed != -1 and grid.two_particles == None:
        vn_ex_ga[grid.background_mesh.i_ghost_fixed] = 0


#    #with multiple ghost interfaces (multiple particles)
#    [setcsrrow2id(A_ga, ga_ghost_points[i]) for i in range(0, len(ga_ghost_points))]
#    for i3 in range(0, int(len(grid.triple_points_index) / 2)):
#        triple_i3 = grid.interface_gammaalpha.triple_points_index[grid.interface_gammaalpha.ghost_points_index_which_param[grid.interface_gammaalpha.triple_points_index] == i3]
#        vn_ga_ghost_fit = la.solve(np.array([[1, 0, 0, 0],
#                                             [1, 1, 1, 1],
#                                             [0, 1, 0, 0],
#                                             [0, 1, 2, 3]]),
#                                    np.array([vn_ga[grid.interface_gammaalpha.triple_points_index[triple_i3[0]]], #v_t0
#                                              vn_ga[grid.interface_gammaalpha.triple_points_index[triple_i3[1]]], #v_t1
#                                              0, #dv/dt|x0
#                                              0])) #dv/dt|x0
#
#        vn_ga_ghost_parametrisation = np.sum(vn_ga_ghost_fit * np.vstack((grid.interface_gammaalpha.parametrisation_ghost[i3]**0, grid.interface_gammaalpha.parametrisation_ghost[i3]**1, grid.interface_gammaalpha.parametrisation_ghost[i3]**2, grid.interface_gammaalpha.parametrisation_ghost[i3]**3)).T, axis = 1)
#        f_ga[grid.interface_gammaalpha.ghost_points_index[grid.interface_gammaalpha.ghost_points_index_which_param == i3]] = vn_ga_ghost_parametrisation



    A_t = A.copy()
    ### Diriclet
    [setcsrrow2id(A_t, t_points[i]) for i in range(0, len(t_points))]
    f_t = np.zeros(grid.enriched_mesh.n)
    f_t[gt_points] = vn_t[0]
    f_t[at_points] = vn_t[1]
    
    vn_ex_t = spla.spsolve(A_t, f_t)

    vn_ex_ga_background = np.zeros(grid.background_mesh.n)
    vn_ex_t_background = np.zeros(grid.background_mesh.n)

    #Mapping back to background mesh
    for i in range(0, grid.background_mesh.n):
        pi = grid.background_mesh.points[i]
        patch_i = i
        update_count = 0
        has_update = 0
        if grid.has_point_been_shifted[i,0] == -1 or i == grid.background_mesh.i_ghost_fixed:
            vn_ex_ga_background[i] = vn_ex_ga[i]
            vn_ex_t_background[i] = vn_ex_t[i]
            has_update = 1
        elif grid.is_point_BP[i,0] == 1:
            pi_E = grid.enriched_mesh.points[i]
            if grid.has_point_been_shifted[i,1]%10 == 1:
                pass
            elif grid.has_point_been_shifted[i,1] == 2 or grid.has_point_been_shifted[i,1] == 3:
                pass
            else:
                print('shifted point over outer boundary not ga(ghost) or gt/at. i, has_been_shift1', i, grid.has_point_been_shifted[i,1])
            patch_i = np.setdiff1d(grid.background_mesh.triangles_points_index[grid.background_mesh.point_in_tri[i].nonzero()[1]].flatten(), i)
            where_BP,  = np.where(grid.is_point_BP[patch_i, 1] == grid.is_point_BP[i,1])
            where_corner, = np.where(grid.is_point_BP[patch_i, 1] == 5)
            where_BP = np.union1d(where_BP, where_corner)
            for xjj in patch_i[where_BP]:
                edgeii = np.hstack((i, xjj))
                edge = grid.background_mesh.points[edgeii]
                #Check if point is shifted over right/left outer boundary or upper/lower
                if grid.is_point_BP[i,1] == 2 or grid.is_point_BP[i,1] == 4: #2/4 right/left resp.
                    t = (edge[0, 1] - pi_E[1]) / (edge[0,1] - edge[1,1])
                elif grid.is_point_BP[i,1] == 1 or grid.is_point_BP[i,1] == 3: #1/3 upper/lower resp.
                    t = (edge[0, 0] - pi_E[0]) / (edge[0,0] - edge[1,0])
                if t >= 0 and t <= 1:
                    vn_ex_ga_background[i] = (1-t)*vn_ex_ga[edgeii[0]] + t*vn_ex_ga[edgeii[1]]
                    vn_ex_t_background[i]  = (1-t)*vn_ex_t[edgeii[0]] + t*vn_ex_t[edgeii[1]]
                    has_update = 1
        while has_update != 1:
            patch_tri_i = grid.enriched_mesh.point_in_tri[patch_i].nonzero()[1]
            j = 0
            while has_update == 0 and j < len(patch_tri_i):
                tri_i = patch_tri_i[j]
                trii = grid.enriched_mesh.triangles_points_index[tri_i]
                tri = grid.enriched_mesh.triangles[tri_i]
                lambda_barycentric = la.solve( np.array([[tri[0,0], tri[1,0], tri[2,0]],
                                                         [tri[0,1], tri[1,1], tri[2,1]],
                                                         [1       , 1       , 1       ] ]), np.array([ pi[0], pi[1], 1]))
                if all(lambda_barycentric <= 1) and all(lambda_barycentric >= 0):
                    vn_ex_ga_background[i] = interpolate(trii, tri, pi, vn_ex_ga)
                    vn_ex_t_background[i] = interpolate(trii, tri, pi, vn_ex_t)
                    has_update = 1
                j = j+1
            patch_trii = grid.enriched_mesh.triangles_points_index[patch_tri_i]
            if has_update != 1 and isinstance(patch_i, int):
                tmin = 10
                j = 0
                while has_update == 0 and j < len(patch_trii):
                    trii = patch_trii[j]
                    other_ii = np.setdiff1d(trii, i)
                    for edgeii in np.array([[i, other_ii[0]], [i, other_ii[1]]]):
                        edge = grid.enriched_mesh.points[edgeii]
                        t = (edge[0] - pi) / (edge[0] - edge[1])
                        if np.abs(t[0]-t[1]) < tmin:
                            tmin = np.abs(t[0]-t[1])
                        if all(t >= 0) and all(t <= 1) and np.abs(t[0]-t[1]) < 1e-10:
                            t = np.mean(t)
                            vn_ex_ga_background[i] = (1-t)*vn_ex_ga[edgeii[0]] + t*vn_ex_ga[edgeii[1]]
                            vn_ex_t_background[i]  = (1-t)*vn_ex_t[edgeii[0]] + t*vn_ex_t[edgeii[1]]
                            has_update = 1
                    j = j+1
            patch_i = np.setdiff1d(patch_trii.flatten(), i)
            update_count = update_count + 1
            if update_count == 5:
                print('mapping v^ex back to background takes (too?) big patches to find background point i in enriched mesh', i, pi)
    return vn_ex_ga_background, vn_ex_t_background

