Skip to content
Snippets Groups Projects
Select Git revision
  • 5bc09c02f6d01a28ea394433ee00dcf3529e29df
  • main default protected
2 results

main.py

Blame
  • Code owners
    Assign users and groups as approvers for specific file changes. Learn more.
    main.py 26.68 KiB
    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 = 10
        
        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())