# -*- coding: utf-8 -*-
###########################################################################
################        Yassine CutFEM           ##########################
###########################################################################
import getfem as gf
import numpy as np
import matplotlib.pyplot as plt

## Parameters
NX = 8                          # Mesh parameter.
Nb_NX = 4
Tab_NX = []
Tab_E_L2 = []
Tab_E_H1 = []
for i in range(Nb_NX):
    NX=NX*2
    print("The value of NX is")
    print(NX)
    Id_Omega_h = 1
    Id_Gamma_h = 2
    Id_Edges = 3
    # Parameter du cercle domaine de calcul
    x0 = 0.
    y0 = 0.
    r = 0.34
    
    Coef_gamma = 5.0*NX
    
    # Create a simple cartesian mesh
    m = gf.Mesh('triangles grid', np.arange(-.5,.5+1./NX,1./NX),
                np.arange(-.5,.5+1./NX,1./NX))
    
    # Create a MeshFem for u and rhs fields of dimension 1 (i.e. a scalar field)
    mfu   = gf.MeshFem(m, 1)
    ls = gf.LevelSet(m,2,'(x-{X0})*(x-{X0}) + (y-{Y0})*(y-{Y0}) - {R}*{R}'.format(X0=x0, Y0=y0, R=r))
    mls = gf.MeshLevelSet(m)
    mls.add(ls)
    mls.adapt()
    mls.cut_mesh().export_to_pos('MLS.pos')
    # assign the P1 fem to all convexes of the both MeshFem
    mfu.set_fem(gf.Fem('FEM_PK(2,1)'))
    
    #  Integration method used
    mim = gf.MeshIm('levelset',mls,'inside', gf.Integ('IM_TRIANGLE(6)'))
    mim_bound = gf.MeshIm('levelset',mls,'boundary', gf.Integ('IM_TRIANGLE(6)'))
    mim.set_integ(6)
    #mfu.dof_from_im(mim)
    
    # Reperage des domaines 
    Tab_CNV_Ind_Omega = mim.convex_index()
    m.set_region(Id_Omega_h, m.faces_from_cvid(Tab_CNV_Ind_Omega))
    #print("Omega_h")
    #print(m.region(Id_Omega_h))
    
    # Set the region Gamma_h
    m.set_region(Id_Gamma_h, m.outer_faces(Tab_CNV_Ind_Omega))
    #print("Gamma_h")
    #print(m.region(Id_Gamma_h))
    
    # Set the list of internal edges for the Ghost penalty
    Tab_band_Gamma = mim_bound.convex_index()
    Tab_all_ed_band = m.faces_from_cvid(Tab_band_Gamma)
    m.set_region(Id_Edges,Tab_all_ed_band)
    m.set_region(6,m.outer_faces(Tab_CNV_Ind_Omega))
    m.region_subtract(Id_Edges, 6)
    #print("the liste of inter edges in the band")
    #print(m.region(Id_Edges))

    # Model
    md = gf.Model('real')
    
    # Main unknown
    md.add_fem_variable('u', mfu)
    
    # Interpolate the exact solution (Assuming mfu is a Lagrange fem)
    Uexpr = ('(np.sqrt(np.square(x-%f)+np.square(y-%f))**3-%f**3)/9.0'%(x0,y0,r));
    Uex = mfu.eval(Uexpr,globals(),locals())
    
    # Laplacian term on u
    md.add_Laplacian_brick(mim, 'u')
    
    # Volumic source term
    Fexpr = ('-np.sqrt(np.square(x-%f)+np.square(y-%f))'%(x0,y0));
    F = mfu.eval(Fexpr,globals(),locals())
    md.add_initialized_fem_data('f', mfu, F)
    md.add_source_term_brick(mim, 'u', 'f')
    
    
    # Add the first Boundary condition on Gamma_h
    Var_Normal_Grad_u = "Grad_u.Normal"
    Test_of_u = "Test_u"
    md.add_linear_term(mim_bound, "-({S})*({K})".format(S=Var_Normal_Grad_u, K=Test_of_u), -1)

    #Add the second Boudary term on Gamma
    var_u = "u"
    Normal_test_of_u = "Grad_Test_u.Normal"
    md.add_linear_term(mim_bound, "-({L})*({G})".format(L=var_u, G=Normal_test_of_u), -1)
    
    # Add the 3th condition with penalization
    #g = mfu.eval('0')
    #md.add_initialized_fem_data('g', mfu, g)
    md.add_initialized_data('g', [0.0])
    md.add_Dirichlet_condition_with_penalization(mim_bound, 'u', Coef_gamma, -1, 'g', mfu)
    
    # Add the 4th Boundary condition Ghost Penalty
    md.add_initialized_data('Sigma_h', [0.1/NX])
    jump = "((Grad_u-Interpolate(Grad_u,neighbour_elt)).Normal)"
    test_jump = "((Grad_Test_u-Interpolate(Grad_Test_u,neighbour_elt)).Normal)"
    md.add_linear_term(mim_bound, "Sigma_h*({M}).({N})".format(M=jump, N=test_jump), Id_Edges)
    #md.add_Dirichlet_condition_with_penalization(mim_bound,
    #md.add_linear_term(mim_bound, 'Sigma_h*(Xfem_plus(Grad_u)-Xfem_minus(Grad_u)).(Xfem_plus(Grad_Test_u)-Xfem_minus(Grad_Test_u))', Id_Edges)
    
    # Assembly of the linear system and solve.
    md.solve()
    
    # Main unknown
    Uap = md.variable('u')
    
    # Slice pour la restriction de la solution et le calcul de l'erreur
    sl = gf.Slice(('comp',('ball', +1, [x0, y0], r)), m, 5)
    # Restriction of the solution to the physical domain
    sl.export_to_vtk('App_Solution.vtk', 'ascii', mfu, Uap)
    sl.export_to_vtk('Exc_Solution.vtk', 'ascii', mfu, Uex)
    sl.export_to_vtk('Error-Ex-App.vtk', 'ascii', mfu, Uap-Uex)
    # The approximate solution on the fictititious domain
    mfu.export_to_pos('Solution.pos', Uap,'Approximate solution')
    mfu.export_to_pos('ExSolution.pos', Uex,'Exact solution')
    mfu.export_to_pos('Poisson-Exact-App.pos', Uex,'Exact solution',
                                        Uap,'Computed solution')
    
    L2error = gf.compute(mfu, Uap-Uex, 'L2 norm', mim)
    H1error = gf.compute(mfu, Uap-Uex, 'H1 norm', mim)
    print('Error in L2 norm : ', L2error)
    print('Error in H1 norm : ', H1error)
                                           
    # Sauvegarde des résultats pour le calcul de l'ordre de convergence 
    Tab_NX.append(np.log(1./NX))
    Tab_E_L2.append(np.log(L2error))
    Tab_E_H1.append(np.log(H1error))

# plot the convergence rate 
plt.plot(Tab_NX, Tab_E_H1, 'v-', Tab_NX, Tab_E_L2, 'o-')
plt.legend(["H1 error", "L2 error"])
plt.axes()
plt.xlabel("Log of the meshsize")
plt.ylabel('Log of the error')
plt.show()

# Calcul d'ordre de convergence pour les deux normes
Order_L2 =[]
Order_H1 =[]
for i in range(1,Nb_NX):
    OL2 = (Tab_E_L2[i]-Tab_E_L2[i-1])/(Tab_NX[i]-Tab_NX[i-1])
    OH1 = (Tab_E_H1[i]- Tab_E_H1[i-1])/(Tab_NX[i]-Tab_NX[i-1])
    Order_L2.append(OL2)
    Order_H1.append(OH1)
print("L'ordre de convergence selon la norme L2")
print(Order_L2)
print("L'ordre de convergence selon la norme H1")
print(Order_H1)