Upload New File

CutOffBoundaryFunsNoise/Gamma5/NoiseLevel1/main.py

+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('AdaptMedLar.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())
+
+    
+    
+    
+    
+    
+
+    
+    
+