from numpy import array, linspace, sqrt, linalg, where
from numpy import logical_not, logical_or, logical_and
import scipy.signal as signal
import scipy.io
import ufl



def sig(Vs):
    """ Conductivity function """
    # Initialization
    xy = Vs.tabulate_dof_coordinates().reshape((-1,2))
    
    N = xy[:,0].size
    
    #Test case 1:
    #a = 2.0
    #R = 0.8
    #sig = 1 + np.where(np.power(xy[:,0],2)+np.power(xy[:,1],2) <= np.power(R,2),np.exp(a-np.divide(a,1-np.divide(np.power(xy[:,0],2)+np.power(xy[:,1],2),np.power(R,2)))),0)
    
    #Test case 2:
    sig = 1 + np.where(np.power(xy[:,0]+1/2,2)+np.power(xy[:,1],2) <= np.power(0.3,2),1,0)  + np.where(np.power(xy[:,0],2)+np.power(xy[:,1]+1/2,2) <= np.power(0.1,2),1,0) + np.where(np.power(xy[:,0]-1/2,2)+np.power(xy[:,1]-1/2,2) <= np.power(0.1,2),1,0)

    
    sigsqrt = np.sqrt(sig)
    
    return sig,sigsqrt

   


if __name__ == '__main__':
    import sys
    # Initialize numpy, IO, matplotlib
    import numpy as np
    import scipy.io as io
    #import matplotlib
    #matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    from matplotlib import ticker
    #import matplotlib
    
    plt.ion()
    
    # Load 
    #import cmap as cmap
    from dolfin import __version__ as DOLFINversion
    from dolfin import TrialFunction, TestFunction, FunctionSpace, MeshFunction
    from dolfin import project, Point, triangle, MixedElement, SubDomain, Measure
    from dolfin import inner, dot, grad, dx, ds, VectorFunctionSpace, PETScLUSolver, as_backend_type
    from dolfin import Function, assemble, Expression, parameters, VectorFunctionSpace
    from dolfin import DirichletBC, as_matrix, interpolate, as_vector, UserExpression, errornorm, norm
    from dolfin import MeshFunction, cells, solve, DOLFIN_EPS, near, Constant, FiniteElement
    from mshr import generate_mesh, Circle
    
    # Load dolfin plotting
    from dolfin import plot as dolfinplot, File as dolfinFile

    def plotVs(vec,**kwargs):
        """ dolfin.plot-wrapper """
        fn = Function(Vs,vec)
        return dolfinplot(fn,**kwargs)
        
    def plotVs1(vec,**kwargs):
        """ dolfin.plot-wrapper """
        fn = Function(Vs1,vec)
        return dolfinplot(fn,**kwargs)
        
        
    # Define Function spaces
    
    def VSpace(mesh, degree=1): # endeligt underrum af H¹_{\diamond} funktioner i H¹ der integrerer til 0 på randen - det svarer til H¹(\Omega) x \mathbb{R}
        E1 = FiniteElement('P', triangle, degree)
        E0 = FiniteElement('R', triangle, 0)

        return FunctionSpace(mesh, MixedElement(E1,E0))

    def VsSpace(mesh, degree=1): # Corresponding to H^1
        E1 = FiniteElement('CG', triangle, degree)
        return FunctionSpace(mesh, E1)
    
    def VqSpace(mesh, degree=1): # Used for the gradients
        return VectorFunctionSpace(mesh, 'CG', degree)
        
    def Vector(V):
        return Function(V).vector()
    
    
    # Define how to calculate the gradients
    def Solver(Op):
        s = PETScLUSolver(as_backend_type(Op),'mumps') # Constructs the linear operator Ks for the  linear system Ks u = f using LU factorization, where the method 'numps' is used
        return s
    
    def GradientSolver(Vq,Vs): #Calculate derivatives on the quadrature points
        """
        Based on:
        https://fenicsproject.org/qa/1425/derivatives-at-the-quadrature-points/
        """
        uq = TrialFunction(Vq)
        vq = TestFunction(Vq)
        M = assemble(inner(uq,vq)*dx)
        femSolver = Solver(M)
    
        u = TrialFunction(Vs)
        P = assemble(inner(vq,grad(u))*dx)

        def GradSolver(uvec):
            gv = Vector(Vq)
            g = P*uvec
            femSolver.solve(gv, g)
            dx = Vector(Vs)
            dy = Vector(Vs)
            dx[:] = gv[0::2].copy()
            dy[:] = gv[1::2].copy()
            return dx,dy

        return GradSolver
    
    # Define the mesh for the unit disk:
    
    def UnitCircleMesh(n):
        C = Circle(Point(0,0),1)
        return generate_mesh(C,n)
        
    #cmaps = cmap.twilights()

    # ------------------------------    
    # Setup mesh and FEM-spaces
    
    parameters['allow_extrapolation'] = True
    
    #Mesh to generate the power density data:
    #Ms = 150
    Ms = 200 #For N_{medium} as in Table 2
    #Ms = 250 #For N_{large} as in Table 2
    
    m = UnitCircleMesh(Ms)
    V = VSpace(m,1) 
    Vs = VsSpace(m,1) 
    Vq = VqSpace(m,1)
    
    #Mesh to solve the inverse problem:
    #Ms2 = 100
    Ms2 = 160 #For N_{medium} as in Table 2
    #Ms2 = 200 #For N_{large} as in Table 2
    
    m2 = UnitCircleMesh(Ms2)
    #V1 = VSpace(m2,1) 
    Vs1 = VsSpace(m2,1)
    Vq1 = VqSpace(m2,1)
    
    # Gradient solver
    GradSolver = GradientSolver(Vq,Vs)
    GradSolver2 = GradientSolver(Vq1,Vs1)
    
    xy = Vs.tabulate_dof_coordinates().reshape((-1,2))
    xy2 = Vs1.tabulate_dof_coordinates().reshape((-1,2))
    
    N = xy[:,0].size
    print(N)
    N2 = xy2[:,0].size
    print(N2)
      
    # ------------------------------
    # Conductivity
    sigt = Function(Vs1)
    sigsqrtt = Vector(Vs1)
    
    sig1 = Function(Vs)
    sigsqrt1 = Function(Vs)
    
    sigt1,sigsqrtt1  = sig(Vs)
    sigt.vector()[:],sigsqrtt[:] = sig(Vs1)
    
    sig1.vector().set_local(sigt1)
    sigsqrt1.vector().set_local(sigsqrtt1)
    
    # Plot
    
    plot_settings = {
        #'levels': np.linspace(-1,2,120),
        #'levels': np.linspace(0,5,120),
        'cmap': plt.get_cmap('inferno'),
    }

    
    
    # ------------------------------    
    # Boundary conditions
    
    #Gamma_Mini:
    #f1 = Expression('cos(std::atan2(x[1],x[0]))-(2.0*sqrt(2.0))/pi', degree=2)
    #f2 = Expression('sin(std::atan2(x[1],x[0]))+(2.0*sqrt(2.0)-4.0)/pi', degree=1)
    
    #Gamma_Sma:
    #f1 = Expression('cos(std::atan2(x[1],x[0]))-2.0/pi', degree=2)
    #f2 = Expression('sin(std::atan2(x[1],x[0]))-2.0/pi', degree=1)
    
    #Gamma_SmaMed:
    #f1 = Expression('cos(std::atan2(x[1],x[0]))-(2.0*sqrt(2.0))/(3.0*pi)', degree=2)
    #f2 = Expression('sin(std::atan2(x[1],x[0]))-(2.0*sqrt(2.0)+4.0)/(3.0*pi)', degree=1)
    
    #Gamma_Med:
    f1 = Expression('cos(std::atan2(x[1],x[0]))', degree=2)
    f2 = Expression('sin(std::atan2(x[1],x[0]))-2.0/pi', degree=1)
    
    #Gamma_MedLar:
    #f1 = Expression('cos(std::atan2(x[1],x[0]))+(2.0*sqrt(2.0))/(5.0*pi)', degree=2)
    #f2 = Expression('sin(std::atan2(x[1],x[0]))-(2.0*sqrt(2.0)+4.0)/(5.0*pi)', degree=1)
    
    #Gamma_Lar:
    #f1 = Expression('cos(std::atan2(x[1],x[0]))+2.0/(3.0*pi)', degree=2)
    #f2 = Expression('sin(std::atan2(x[1],x[0]))-2.0/(3.0*pi)', degree=1)
    
    #Gamma_Huge:
    #f1 = Expression('cos(std::atan2(x[1],x[0]))+(2.0*sqrt(2.0))/(7.0*pi)', degree=2)
    #f2 = Expression('sin(std::atan2(x[1],x[0]))+(2.0*sqrt(2.0)-4.0)/(7.0*pi)', degree=1)
    
    
    
    #Gamma_MedLar:
    #class MyExpression1(UserExpression):
    #    def eval(self, value, x):
    #        if (0 < ufl.atan_2(x[1],x[0]) and ufl.atan_2(x[1],x[0]) <= ufl.atan(1)*4.0):
    #            value[0] = ufl.cos(ufl.atan_2(x[1],x[0]))+(2.0*sqrt(2.0))/(5.0*ufl.atan(1)*4.0)
    #        elif (-ufl.atan(1)*4.0 < ufl.atan_2(x[1],x[0]) and ufl.atan_2(x[1],x[0]) < -(3.0/4.0)*ufl.atan(1)*4.0):
    #            value[0] = ufl.cos(ufl.atan_2(x[1],x[0])+2.0*ufl.atan(1)*4.0)+(2.0*sqrt(2.0))/(5.0*ufl.atan(1)*4.0)
    
    #class MyExpression2(UserExpression):
    #    def eval(self, value, x):
    #        if (0 < ufl.atan_2(x[1],x[0]) and ufl.atan_2(x[1],x[0]) <= ufl.atan(1)*4.0):
    #            value[0] = ufl.sin(ufl.atan_2(x[1],x[0]))-(2.0*sqrt(2.0)+4.0)/(5.0*ufl.atan(1)*4.0)
    #        elif (-ufl.atan(1)*4.0 < ufl.atan_2(x[1],x[0]) and ufl.atan_2(x[1],x[0]) < -(3.0/4.0)*ufl.atan(1)*4.0):
    #            value[0] = ufl.sin(ufl.atan_2(x[1],x[0])+2.0*ufl.atan(1)*4.0)-(2.0*sqrt(2.0)+4.0)/(5.0*ufl.atan(1)*4.0)
                
    #f1 = MyExpression1(degree=2)
    #f2 = MyExpression2(degree=2)

    # Defining \Gamma    
    class boundaryN(SubDomain):
    
        #\Gamma_{small}:
        #def inside(self, x, on_boundary):
        #    return (on_boundary and (ufl.atan_2(x[1],x[0])<(1.0/4.0)*np.pi and ufl.atan_2(x[1],x[0])>0) )
    
        #\Gamma_{small}:
        #def inside(self, x, on_boundary):
        #    return (on_boundary and (ufl.atan_2(x[1],x[0])<(1.0/2.0)*np.pi and ufl.atan_2(x[1],x[0])>0) )
            
        #\Gamma_{small}:
        #def inside(self, x, on_boundary):
        #    return (on_boundary and (ufl.atan_2(x[1],x[0])<(3.0/4.0)*np.pi and ufl.atan_2(x[1],x[0])>0) )
        
        #\Gamma_{medium}:
        def inside(self, x, on_boundary):
            return (on_boundary and (ufl.atan_2(x[1],x[0])<np.pi and ufl.atan_2(x[1],x[0])>0) )
        
        #\Gamma_{medlar}:    
        #def inside(self, x, on_boundary):
        #    return (on_boundary and ((ufl.atan_2(x[1],x[0])<-3.0/4.0*np.pi and ufl.atan_2(x[1],x[0])>-np.pi)  or (ufl.atan_2(x[1],x[0])<np.pi and ufl.atan_2(x[1],x[0])>0)))
        
        #\Gamma_{large}:
        #def inside(self, x, on_boundary):
        #    return (on_boundary and ((ufl.atan_2(x[1],x[0])<-1.0/2.0*np.pi and ufl.atan_2(x[1],x[0])>-np.pi)  or (ufl.atan_2(x[1],x[0])<np.pi and ufl.atan_2(x[1],x[0])>0)))
        
        #\Gamma_{huge}:
        #def inside(self, x, on_boundary):
        #    return (on_boundary and ((ufl.atan_2(x[1],x[0])<-1.0/4.0*np.pi and ufl.atan_2(x[1],x[0])>-np.pi)  or (ufl.atan_2(x[1],x[0])<np.pi and ufl.atan_2(x[1],x[0])>0)))
        
        #Full:
        #def inside(self, x, on_boundary):
        #    return on_boundary
        
    
        
    bN = boundaryN()
    
    boundary_markers = MeshFunction("size_t",m,m.topology().dim()-1,0)
    boundary_markers.set_all(9999)
    bN.mark(boundary_markers,0)
    
    ds = Measure('ds', domain=m,subdomain_data=boundary_markers)
   
        
        
    def boundary2(x, on_boundary):
        return on_boundary

    
    
    
    (u1,c1) = TrialFunction(V)
    (v1,d1) = TestFunction(V)
    (u2,c2) = TrialFunction(V)
    (v2,d2) = TestFunction(V)
    
    
    #Defining and solving the variational equations
    a1 = (inner(sig1*grad(u1),grad(v1))+c1*v1+u1*d1)*dx
    a2 = (inner(sig1*grad(u2),grad(v2))+c2*v2+u2*d2)*dx 
    
    L1 = f1*v1*ds(0)
    L2 = f2*v2*ds(0)
    
    #L1 = f1*v1*ds
    #L2 = f2*v2*ds
    
    w1 = Function(V)
    w2 = Function(V)
    
    solve(a1 == L1,w1)
    solve(a2 == L2,w2)
    
    (u1,c1) = w1.split()
    (u2,c2) = w2.split()
    
    u1new = interpolate(u1,Vs)
    u2new = interpolate(u2,Vs)
    
    #Defining the gradients
    dU1 = GradSolver(u1new.vector())
    dU2 = GradSolver(u2new.vector())

    H11tmp = Vector(Vs)
    H12tmp = Vector(Vs)
    H22tmp = Vector(Vs)
    
    #Compute the noise free power density data
    H11tmp[:]  = sig1.vector()*(dU1[0]*dU1[0]+dU1[1]*dU1[1])
    H12tmp[:]  = sig1.vector()*(dU1[0]*dU2[0]+dU1[1]*dU2[1])
    H22tmp[:]  = sig1.vector()*(dU2[0]*dU2[0]+dU2[1]*dU2[1])
 
    
    H11t = Function(Vs)
    H12t = Function(Vs)
    H21t = Function(Vs)
    H22t = Function(Vs)
    
    
    #Add noise to the data by following the approach in section 5.4
    np.random.seed(50) 
    
    e1 = Vector(Vs)
    e2 = Vector(Vs)
    e3 = Vector(Vs)
    e4 = Vector(Vs)
    
    e1[:]= np.random.randn(N)
    e2[:]= np.random.randn(N)
    e3[:]= np.random.randn(N)
    e4[:]= np.random.randn(N)
    
    #Define the noise level
    alpha = 1
    
    H11t.vector()[:] = H11tmp + alpha/100.0*(e1/norm(e1))*H11tmp
    H12t.vector()[:] = H12tmp + alpha/100.0*(e2/norm(e2))*H12tmp
    H21t.vector()[:] = H12tmp + alpha/100.0*(e3/norm(e3))*H12tmp
    H22t.vector()[:] = H22tmp + alpha/100.0*(e4/norm(e4))*H22tmp
    
    S11 = Function(Vs)
    S12 = Function(Vs)
    S21 = Function(Vs)
    S22 = Function(Vs)
    
    S11.vector()[:] = sigsqrt1.vector()*dU1[0]
    S21.vector()[:] = sigsqrt1.vector()*dU1[1]
    S12.vector()[:] = sigsqrt1.vector()*dU2[0]
    S22.vector()[:] = sigsqrt1.vector()*dU2[1]
    
    #Project the noisy data to the mesh for solving the inverse problem
    H11tmp = project(H11t,Vs1)
    H12tmp = project(H12t,Vs1)
    H21tmp = project(H21t,Vs1)
    H22tmp = project(H22t,Vs1)
    
    H11 = Function(Vs1)
    H12 = Function(Vs1)
    H22 = Function(Vs1)
    
    #Define a lower bound for the eigenvalues of \tilde{H}
    lb = 1e-6
    
    lbvec = Vector(Vs1)
    
    for i in range(N2):
        H11mat = H11tmp.vector()[i]
        H12mat = (1/2)*(H12tmp.vector()[i] + H21tmp.vector()[i])
        H22mat = H22tmp.vector()[i]
        Hmat = np.array([[H11mat,H12mat],[H12mat,H22mat]])
        wh,vh = np.linalg.eig(Hmat)
        H11.vector()[i] = np.maximum(wh[0],lb)*vh[0,0]*vh[0,0] + np.maximum(wh[1],lb)*vh[0,1]*vh[0,1]
        H12.vector()[i] = np.maximum(wh[0],lb)*vh[0,0]*vh[1,0] + np.maximum(wh[1],lb)*vh[0,1]*vh[1,1]
        H22.vector()[i] = np.maximum(wh[0],lb)*vh[1,0]*vh[1,0] + np.maximum(wh[1],lb)*vh[1,1]*vh[1,1]
        if (np.maximum(wh[1],lb)==lb or np.maximum(wh[0],lb)==lb):
            lbvec[i] = 1
    
    plot_settingslb = {
        'levels': np.linspace(0,1.01,250),
        #'levels': np.linspace(0,5,120),
        'cmap': plt.get_cmap('cool'),
    }
    
    plt.figure(20).clear()
    h = plotVs1(lbvec,**plot_settingslb) # dolfinplot
    plt.gca().axis('off')
    cb=plt.colorbar(h)
    cb.ax.tick_params(labelsize=20)
    tick_locator = ticker.MaxNLocator(nbins=6)
    cb.locator = tick_locator
    cb.update_ticks()
    
    for h1 in h.collections:
        h1.set_edgecolor("face")
      
        
    plt.savefig('lb1e6.pdf',format='pdf')
    
    
    
    
    H11tlog = Vector(Vs1)
    H11tlog[:] = np.log(H11tmp.vector())
    
    H11log = Vector(Vs1)
    H11log[:] = np.log(H11.vector())

    
    plt.figure(31).clear()
    h = plotVs1(H11log,**plot_settings) # dolfinplot
    plt.gca().axis('off')
    cb=plt.colorbar(h)
    cb.ax.tick_params(labelsize=20)
    tick_locator = ticker.MaxNLocator(nbins=6)
    cb.locator = tick_locator
    cb.update_ticks()
    
    for h1 in h.collections:
        h1.set_edgecolor("face")
        
    plt.savefig('H11_1noiselb1e6cont.pdf',format='pdf')
    
    
    S11_1 = project(S11,Vs1)
    S12_1 = project(S12,Vs1)
    S21_1 = project(S21,Vs1)
    S22_1 = project(S22,Vs1)
    

    
    H11log = Vector(Vs1)
    H11log[:] = np.log(H11.vector())
    
    H12log = Vector(Vs1)
    H12log[:] = np.sign(H12.vector())*np.log(1.0+np.abs(H12.vector())/(5e-3))
    #H12log[:] = np.sign(H12.vector())*np.log(np.abs(H12.vector()))
    
    H22log = Vector(Vs1)
    H22log[:] = np.log(H22.vector())
    
    miH11 = np.amin(H11log.get_local())
    maH11 = np.amax(H11log.get_local())
    
    plot_settingsH11 = {
        'levels': np.linspace(miH11,maH11,250),
        #'levels': np.linspace(0,5,120),
        #'cmap': plt.set_cmap('RdBu'),
        'cmap': plt.get_cmap('inferno'),
    }
    
    miH12 = np.amin(H12log.get_local())
    maH12 = np.amax(H12log.get_local())
    
    plot_settingsH12 = {
        'levels': np.linspace(miH12,maH12,250),
        #'levels': np.linspace(0,5,120),
        #'cmap': plt.set_cmap('RdBu'),
        'cmap': plt.get_cmap('inferno'),
    }
    
    miH22 = np.amin(H22log.get_local())
    maH22 = np.amax(H22log.get_local())
    
    plot_settingsH22 = {
        'levels': np.linspace(miH22,maH22,250),
        #'levels': np.linspace(0,5,120),
        #'cmap': plt.set_cmap('RdBu'),
        'cmap': plt.get_cmap('inferno'),
    }
        
    
    plt.figure(30)
    h = plotVs1(H11log,**plot_settingsH11) # dolfinplot
    plt.gca().axis('off')
    cb=plt.colorbar(h)
    cb.ax.tick_params(labelsize=20)
    tick_locator = ticker.MaxNLocator(nbins=6)
    cb.locator = tick_locator
    cb.update_ticks()
    
    for h1 in h.collections:
        h1.set_edgecolor("face")
        
    plt.savefig('H11medlog.pdf',format='pdf')
    
    plt.figure(31)
    h = plotVs1(H12log,**plot_settingsH12) # dolfinplot
    plt.gca().axis('off')
    cb=plt.colorbar(h)
    cb.ax.tick_params(labelsize=20)
    tick_locator = ticker.MaxNLocator(nbins=6)
    cb.locator = tick_locator
    cb.update_ticks()
    
    for h1 in h.collections:
        h1.set_edgecolor("face")
        
    plt.savefig('H12medlog.pdf',format='pdf')
    
    #plt.figure(32).clear()
    h = plotVs1(H22log,**plot_settingsH22) # dolfinplot
    plt.gca().axis('off')
    cb=plt.colorbar(h)
    cb.ax.tick_params(labelsize=20)
    tick_locator = ticker.MaxNLocator(nbins=6)
    cb.locator = tick_locator
    cb.update_ticks()
    
    for h1 in h.collections:
        h1.set_edgecolor("face")
        
    plt.savefig('H22medlog.pdf',format='pdf')
    
    
    
    #sys.exit()
    
    S11_1 = project(S11,Vs1)
    S12_1 = project(S12,Vs1)
    S21_1 = project(S21,Vs1)
    S22_1 = project(S22,Vs1)
    
    dt = H11.vector()*H22.vector()-H12.vector()*H12.vector()

    lJac = Vector(Vs1)
    lJac[:] = np.log(dt)
    
    #lJacold = Vector(Vs1)
    #lJacold[:] = np.log(dtold)
    
    midet = np.amin(dt.get_local())
    madet = np.amax(dt.get_local())
    
    plot_settingsdet = {
        'levels': np.linspace(midet,madet,250),
        #'levels': np.linspace(0,5,120),
        'cmap': plt.get_cmap('inferno'),
    }
    
    mildet = np.amin(lJac.get_local())
    maldet = np.amax(lJac.get_local())
    
    plot_settingslJac = {
        'levels': np.linspace(mildet,maldet,250),
        #'levels': np.linspace(0,5,120),
        'cmap': plt.get_cmap('inferno'),
    }
    

    
    #Illustrating log(det(H)):
    plt.figure(3)
    h = plotVs1(lJac,**plot_settingslJac) # dolfinplot
    plt.gca().axis('off')
    cb=plt.colorbar(h)
    cb.ax.tick_params(labelsize=20)
    tick_locator = ticker.MaxNLocator(nbins=6)
    cb.locator = tick_locator
    cb.update_ticks()
    
    for h1 in h.collections:
        h1.set_edgecolor("face")
    
        
    plt.savefig('logDetCoordMed.pdf',format='pdf') 
    
    #sys.exit()

    dtest = H11.vector()*H22.vector()-H12.vector()*H12.vector()
    
    print(min(dtest))
    
    d = Vector(Vs1)
    d[:] = np.sqrt(dtest)
    
    #Compute the T matrix as in section 4.3 
    T11 = np.divide(1,np.sqrt(H11.vector()))
    T12 = np.zeros(N2)
    T21 = -np.divide(H12.vector(),np.multiply(np.sqrt(H11.vector()),d))
    T22 = np.divide(np.sqrt(H11.vector()),d)

    H1211 = Vector(Vs1)
    H1211[:] = np.divide(H12.vector(),H11.vector())
    dH1211 = GradSolver2(H1211)
    
    #Compute the vector field V21 as in equation (4.7)
    V210 = Vector(Vs1)
    V211 = Vector(Vs1)

    V210[:] = -np.multiply(np.divide(H11.vector(),d),dH1211[0])
    V211[:] = -np.multiply(np.divide(H11.vector(),d),dH1211[1])
    
    ld1 = Vector(Vs1)
    ld1[:] = np.log(np.power(d,2))
    dld1 = GradSolver2(ld1)
    
    #Compute the right hand side \mathbf{F} for the Poisson equation (4.3)
    dtheta0 = Function(Vs1)
    dtheta1 = Function(Vs1)
    dtheta0.vector()[:] = (1/2)*(-V210 + (1/2)*dld1[1])
    dtheta1.vector()[:] = (1/2)*(-V211 - (1/2)*dld1[0])
    
    #Compute the true R and theta
    R11_1 = np.multiply(S11_1.vector(),T11) + np.multiply(S12_1.vector(),T12)
    R21_1 = np.multiply(S21_1.vector(),T11) + np.multiply(S22_1.vector(),T12)
    
    theta1t = Function(Vs1)
    theta1t.vector()[:] = np.angle(R11_1+np.multiply(1j,R21_1))
    
    theta1testfun = Function(Vs1)
    theta1testfun.vector()[:] = theta1t.vector()
    
    ##The modification of \theta defined in equation (5.1) when using \Gamma_{small}:
    
    
    #Gamma_mini:
    #indBdr1tmp = np.where((np.arctan2(xy2[:,1],xy2[:,0])< np.pi/2.0) & (np.arctan2(xy2[:,1],xy2[:,0])>0.7688187211))
    
    
    #indBdr1 = np.asarray(indBdr1tmp)

    #indBdr1 = np.reshape(indBdr1,indBdr1.shape[1])

    
    ##theta1testfun.vector()[indBdr1] = theta1testfun.vector()[indBdr1] + 2*np.pi
    #theta1testfun.vector()[indBdr1] = theta1testfun.vector()[indBdr1] - 2*np.pi
    
    
    ThetFun = Function(Vs1,theta1t.vector())
    ThetFunMod = Function(Vs1,theta1testfun.vector())
    
    ##Illustration of the modification of \theta at the boundary:
    
    #plt.figure(4)
    #ax = plt.gca()
    #Nt = 100
    #ang = np.linspace(-np.pi,np.pi,Nt)
    #r = 1
    #bdryTf = [ThetFun(r*np.cos(t),r*np.sin(t)) for t in ang]
    #bdryTfmod = [ThetFunMod(r*np.cos(t),r*np.sin(t)) for t in ang]
    #ax.plot(ang, bdryTf,'b-',label=r'$\theta^c\vert_{\partial \Omega}(t)$')
    #ax.plot(ang, bdryTfmod,'r--',label=r'$\tilde{\theta}^c\vert_{\partial \Omega}(t)$')
    #ax.legend(loc=2,prop={'size': 16})
    ##plt.yticks([-3*np.pi/2,-np.pi,-np.pi/2,0,np.pi/2, np.pi],[r'-$\frac{3\pi}{2}$',r'-$\pi$',r'-$\frac{\pi}{2}$',0,r'$\frac{\pi}{2}$',r'$\pi$'])
    #plt.xticks([-np.pi,-np.pi/2,0,np.pi/2, np.pi],[r'-$\pi$',r'-$\frac{\pi}{2}$',0,r'$\frac{\pi}{2}$',r'$\pi$'])
    #plt.xlabel(r'$t$',fontsize=20)
    #ax.tick_params(axis='both', which='major', labelsize=20)
    #plt.grid(True)
    ##plt.savefig('ThetaBdrycSma2', bbox_inches="tight")
    
    #for h1 in ax.collections:
    #    h1.set_edgecolor("face")
    
        
    #plt.savefig('ThetaAdapt.pdf',format='pdf') 
    
    #sys.exit()
    
    #plot_settingsT = {
    #    'levels': np.linspace(-np.pi,np.pi,250),
    #    #'levels': np.linspace(0,5,120),
    #    'cmap': cmaps['twilight_shifted']
    #}
    
    #plt.figure(5)
    #h = plotVs1(theta1t.vector(),**plot_settings) # dolfinplot
    #plt.gca().axis('off')
    #cb=plt.colorbar(h)
    #cb.ax.tick_params(labelsize=20)
    #tick_locator = ticker.MaxNLocator(nbins=6)
    #cb.locator = tick_locator
    #cb.update_ticks()
    ##plt.savefig('TrueThetaSmaMod')
 
    bctheta1 = DirichletBC(Vs1, theta1testfun, boundary2)

    theta1 = TrialFunction(Vs1)
    vt1 = TestFunction(Vs1)
    
    #Defining and solving the variational equation for the Poisson problem in (13)
    a = inner(grad(theta1),grad(vt1))*dx 
    L = inner(as_vector([dtheta0,dtheta1]),grad(vt1))*dx
    
    theta1 = Function(Vs1)
    solve(a == L,theta1,bctheta1)
    
    print((errornorm(theta1,theta1t,'L2')/norm(theta1t,'L2'))*100)
    
    #theta1.vector()[:] = theta1testfun.vector()

    
    cos2 = Function(Vs1)
    sin2 = Function(Vs1)
    cos2t = Function(Vs1)
    sin2t = Function(Vs1)
    cos2.vector()[:] = np.cos(2*theta1.vector())
    sin2.vector()[:] = np.sin(2*theta1.vector())
    cos2t.vector()[:] = np.cos(2*theta1t.vector())
    sin2t.vector()[:] = np.sin(2*theta1t.vector())
    
    mis2 = np.amin(cos2.vector().get_local())
    mas2 = np.amax(cos2.vector().get_local())
    
    plot_settingssin2 = {
        'levels': np.linspace(mis2,mas2,250),
        'cmap': plt.get_cmap('inferno'),
    }
    
    #plt.figure(6)
    #h = plotVs1(theta1.vector(),**plot_settings) # dolfinplot
    #plt.gca().axis('off')
    #cb=plt.colorbar(h)
    #cb.ax.tick_params(labelsize=20)
    #tick_locator = ticker.MaxNLocator(nbins=6)
    #cb.locator = tick_locator
    #cb.update_ticks()
    
    #Illustrations of sin(2\theta) and its reconstructed version
    #plt.figure(6).clear()
    #h = plotVs1(sin2t.vector(),**plot_settingssin2) # dolfinplot
    #plt.gca().axis('off')
    #cb=plt.colorbar(h)
    #cb.ax.tick_params(labelsize=20)
    #tick_locator = ticker.MaxNLocator(nbins=6)
    #cb.locator = tick_locator
    #cb.update_ticks()
    
    #for h1 in h.collections:
    #    h1.set_edgecolor("face")
    #    
    #plt.savefig('TrueThetaLarSin.pdf',format='pdf')
    
    #plt.figure(7).clear()
    #h = plotVs1(sin2.vector(),**plot_settingssin2) # dolfinplot
    #plt.gca().axis('off')
    #cb=plt.colorbar(h)
    #cb.ax.tick_params(labelsize=20)
    #tick_locator = ticker.MaxNLocator(nbins=6)
    #cb.locator = tick_locator
    #cb.update_ticks()
    
    #for h1 in h.collections:
    #    h1.set_edgecolor("face")
        
    #plt.savefig('RecThetaLarSin.pdf',format='pdf')
    
    #print((errornorm(cos2,cos2t,'L2')/norm(cos2t,'L2'))*100)
    #print((errornorm(sin2,sin2t,'L2')/norm(sin2t,'L2'))*100)
    
    #Defining the vector fields V11 and V22 defined in equation (4.7)
    V11in = Vector(Vs1)
    V11in[:] = np.log(np.divide(1,np.sqrt(H11.vector())))
    V11 = GradSolver2(V11in)
    V22in = Vector(Vs1)
    V22in[:] = np.log(np.divide(np.sqrt(H11.vector()),d))
    V22 = GradSolver2(V22in)
    
    V22min = np.amin(V22in.get_local())
    V22max = np.amax(V22in.get_local())
    
    plot_settingsV22 = {
        'levels': np.linspace(V22min,V22max,250),
        'cmap': plt.set_cmap('inferno'),
    }
    
    #Illustration of \log(\sqrt(H11)/D) as in figure 8
    #plt.figure(11).clear()
    #h = plotVs1(V22in,**plot_settingsV22) # dolfinplot
    #plt.gca().axis('off')
    #cb=plt.colorbar(h)
    #cb.ax.tick_params(labelsize=20)
    #tick_locator = ticker.MaxNLocator(nbins=6)
    #cb.locator = tick_locator
    #cb.update_ticks()
    
    #Computing the vector field \mathbf{K} for the right hand side in equation (4.4)
    Fc0 = V11[0] - V22[0] + V211
    Fc1 = -V11[1] + V22[1] + V210
        
    #Computing \mathbf{G} for the right hand side in (4.4)
    dlogA0 = Function(Vs1)
    dlogA1 = Function(Vs1)
    
    dlogA0.vector()[:] = np.multiply(np.cos(2*theta1.vector()),Fc0) - np.multiply(np.sin(2*theta1.vector()),Fc1)
    dlogA1.vector()[:] = np.multiply(np.cos(2*theta1.vector()),Fc1) + np.multiply(np.sin(2*theta1.vector()),Fc0)
    
    logAt = Function(Vs1)
    logAt.vector()[:] = np.log(sigt.vector())
    
    
    bclogA = DirichletBC(Vs1, logAt, boundary2)
    
    
    logA = TrialFunction(Vs1)
    vA = TestFunction(Vs1)
    
    #Defining and solving the variational formulation of the Poisson problem (4.5)
    aA = inner(grad(logA),grad(vA))*dx 
    LA = inner(as_vector([dlogA0,dlogA1]),grad(vA))*dx

    logA = Function(Vs1)
    solve(aA == LA,logA,bclogA)
    
    A = Function(Vs1)
    A.vector()[:] = np.exp(logA.vector())
    
    plot_settings4 = {
        'levels': np.linspace(0.0,5.1,250),
        'cmap': plt.set_cmap('inferno'),
    }
    
    plot_settings5 = {
        'levels': np.linspace(0.8,5.1,250),
        'cmap': plt.set_cmap('inferno'),
    }

    
    #plt.savefig('TrueSigma2.pdf',format='pdf')
    
    #plt.figure(9)
    #h = plotVs1(sigt.vector(),**plot_settings5) # dolfinplot
    #plt.gca().axis('off')
    #cb=plt.colorbar(h)
    #cb.ax.tick_params(labelsize=20)
    #tick_locator = ticker.MaxNLocator(nbins=6)
    #cb.locator = tick_locator
    #cb.update_ticks()
    
    #for h1 in h.collections:
    #    h1.set_edgecolor("face")
    #    
    #plt.savefig('TrueSigma.pdf',format='pdf')
    
    plt.figure(11)
    h = plotVs1(A.vector(),**plot_settings5) # dolfinplot
    plt.gca().axis('off')
    #cb=plt.colorbar(h)
    #cb.ax.tick_params(labelsize=20)
    #tick_locator = ticker.MaxNLocator(nbins=6)
    #cb.locator = tick_locator
    #cb.update_ticks()
    
    for h1 in h.collections:
        h1.set_edgecolor("face")
        
    plt.savefig('MedCoord1.pdf',format='pdf')
    
    #Computing the relative L2-error
    print((errornorm(A,sigt,'L2')/norm(sigt,'L2'))*100)
    
    
    #Saving the data to a file to do illustrations in Matlab
    io.savemat('AdaptMedLar.mat',{'xy2':xy2,
    'AdaptMedLar':A.vector().get_local(),
    'AdaptTrue': sigt.vector().get_local()})
    
    np.savez('AdaptMedLar.npz',xy2=xy2,A=A.vector().get_local(),sigt=sigt.vector().get_local())