Initial commit

ContBC/cmap/__init__.py

+from .cmap import *

ContBC/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))
+    
+    #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.pyplot as plt
+    from matplotlib import ticker
+    plt.ion()
+    
+    # Load 
+    import cmap as cmap
+    from dolfin import __version__ as DOLFINversion
+    from dolfin import TrialFunction, TestFunction, FunctionSpace
+    from dolfin import project, Point, triangle
+    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 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
+    
+    #Mesh to generate the power density data:
+    Ms = 200
+    m = UnitCircleMesh(Ms)
+    Vs = VsSpace(m,1) 
+    Vq = VqSpace(m,1)
+    
+    #Mesh to solve the inverse problem:
+    Ms2 = 160
+    m2 = UnitCircleMesh(Ms2)
+    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
+    N2 = xy2[:,0].size
+      
+    # ------------------------------
+    # 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.set_cmap('RdBu'),
+    }
+
+    
+    
+    # ------------------------------    
+    # Boundary conditions for the different sizes of \Gamma
+    
+    #Gamma_{small}:
+    #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]) < (1.0/4.0)*ufl.atan(1)*4.0:
+    #            value[0] = ufl.cos(8*ufl.atan_2(x[1],x[0]))-1
+    #        else :
+    #            value[0] = 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]) < (1.0/4.0)*ufl.atan(1)*4.0:
+    #            value[0] = ufl.sin(8*ufl.atan_2(x[1],x[0]))
+    #        else :
+    #            value[0] = 0
+    
+    #Gamma_{medium}:
+    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(2*ufl.atan_2(x[1],x[0]))-1
+            else :
+                value[0] = 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(2*ufl.atan_2(x[1],x[0]))
+            else :
+                value[0] = 0
+    
+    #Gamma_{large]: 
+    #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((8.0/7.0)*ufl.atan_2(x[1],x[0]))-1
+    #        elif (-ufl.atan(1)*4.0 < ufl.atan_2(x[1],x[0]) and ufl.atan_2(x[1],x[0]) < -(1.0/4.0)*ufl.atan(1)*4.0):
+    #            value[0] = ufl.cos((8.0/7.0)*(ufl.atan_2(x[1],x[0])-1.5*ufl.atan(1)*4.0))-1
+    #        else :
+    #            value[0] = 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((8.0/7.0)*ufl.atan_2(x[1],x[0]))
+    #        elif (-ufl.atan(1)*4.0 < ufl.atan_2(x[1],x[0]) and ufl.atan_2(x[1],x[0]) < -(1.0/4.0)*ufl.atan(1)*4.0):
+    #            value[0] = ufl.sin((8.0/7.0)*(ufl.atan_2(x[1],x[0])-1.5*ufl.atan(1)*4.0))
+    #        else :
+    #            value[0] = 0
+    
+    u_D1 = MyExpression1(degree=2)
+    u_D2 = MyExpression2(degree=2)
+        
+    def boundary(x, on_boundary):
+        return on_boundary
+    
+    bc1 = DirichletBC(Vs, u_D1, boundary)
+    bc2 = DirichletBC(Vs, u_D2, boundary)
+
+    u1 = Function(Vs)
+    v1 = TestFunction(Vs)
+    u2 = Function(Vs)
+    v2 = TestFunction(Vs)
+    
+    
+    #Defining and solving the variational equations
+    a1 = inner(sig1*grad(u1),grad(v1))*dx
+    a2 = inner(sig1*grad(u2),grad(v2))*dx 
+    
+    solve(a1 == 0,u1,bc1)
+    solve(a2 == 0,u2,bc2)
+    
+    #Defining the gradients
+    dU1 = GradSolver(u1.vector())
+    dU2 = GradSolver(u2.vector())
+
+    H11t = Function(Vs)
+    H12t = Function(Vs)
+    H22t = Function(Vs)
+    
+    #Compute the power density data
+    H11t.vector()[:]  = sig1.vector()*(dU1[0]*dU1[0]+dU1[1]*dU1[1])
+    H12t.vector()[:]  = sig1.vector()*(dU1[0]*dU2[0]+dU1[1]*dU2[1])
+    H22t.vector()[:]  = sig1.vector()*(dU2[0]*dU2[0]+dU2[1]*dU2[1])
+    
+    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 data to the mesh for solving the inverse problem
+    H11 = project(H11t,Vs1)
+    H12 = project(H12t,Vs1)
+    H22 = project(H22t,Vs1)
+    
+    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)
+    
+    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.set_cmap('RdBu'),
+    }
+    
+    mildet = np.amin(lJac.get_local())
+    maldet = np.amax(lJac.get_local())
+    
+    plot_settingsldet = {
+        'levels': np.linspace(mildet,maldet,250),
+        #'levels': np.linspace(0,5,120),
+        #'cmap': plt.set_cmap('RdBu'),
+        'cmap': plt.get_cmap('bwr'),
+    }
+    
+    miu1 = np.amin(u1.vector().get_local())
+    mau1 = np.amax(u1.vector().get_local())
+    
+    plot_settingsu1 = {
+        'levels': np.linspace(miu1,mau1,250),
+        #'levels': np.linspace(0,5,120),
+        #'cmap': plt.set_cmap('RdBu'),
+        'cmap': plt.get_cmap('bwr'),
+    }
+    
+    miu2 = np.amin(u2.vector().get_local())
+    mau2 = np.amax(u2.vector().get_local())
+    
+    plot_settingsu2 = {
+        'levels': np.linspace(miu2,mau2,250),
+        #'levels': np.linspace(0,5,120),
+        #'cmap': plt.set_cmap('RdBu'),
+        'cmap': plt.get_cmap('bwr'),
+    }
+    
+    #Illustrating the solutions u1 and u2
+    plt.figure(1).clear()
+    h = plotVs(u1.vector(),**plot_settingsu1) # 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('u1medc1')
+    
+    plt.figure(2).clear()
+    h = plotVs(u2.vector(),**plot_settingsu2) # 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('u2medc1')
+    
+    #Illustrating log(det(H)):
+    plt.figure(3).clear()
+    h = plotVs1(lJac,**plot_settingsldet) # 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('logDetMedc1') 
+
+    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 (17)
+    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 (13)
+    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 when using \Gamma_{small} defined in equation (20):
+    
+    #indBdr1tmp = np.where((np.arctan2(xy2[:,1],xy2[:,0])<np.pi) & (np.arctan2(xy2[:,1],xy2[:,0])>np.pi/8.0+0.068))
+    
+    #indBdr1 = np.asarray(indBdr1tmp)
+
+    #indBdr1 = np.reshape(indBdr1,indBdr1.shape[1])
+
+    
+    #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).clear()
+    #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")
+    
+    plot_settingsT = {
+        'levels': np.linspace(-np.pi,0,250),
+        #'levels': np.linspace(0,5,120),
+        'cmap': cmaps['twilight_shifted']
+    }
+    
+    #plt.figure(5).clear()
+    #h = plotVs1(theta1t.vector(),**plot_settingsT) # 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, boundary)
+    
+    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)
+    
+    miT = np.amin(theta1.vector().get_local())
+    maT = np.amax(theta1.vector().get_local())
+    
+    plot_settingstest = {
+        'levels': np.linspace(miT,maT,250),
+        'cmap': plt.get_cmap('bwr'),
+    }
+    
+    #Illustrating the true vs. the reconstructed \theta-function
+    plt.figure(6).clear()
+    h = plotVs1(theta1t.vector(),**plot_settingstest) # 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('TrueThetaLarc2')
+    
+    plt.figure(7).clear()
+    h = plotVs1(theta1.vector(),**plot_settingstest) # 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('RecThetaLarc2')
+    
+    #Compute the relative L2-error for \theta
+    print((errornorm(theta1,theta1t,'L2')/norm(theta1t,'L2'))*100)
+    
+    #Calculating the relative errors for the reconstruction of \cos(2\theta) and \sin(2\theta) in the case of \Gamma_{small}:
+    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())
+    
+    mic2 = np.amin(cos2.vector().get_local())
+    mac2 = np.amax(cos2.vector().get_local())
+    
+    plot_settingscos2 = {
+        'levels': np.linspace(mic2,mac2,250),
+        'cmap': plt.get_cmap('bwr'),
+    }
+    
+    #Illustrations of cos(2\theta) and its reconstructed version in the case of \Gamma_{small}
+    #plt.figure(8).clear()
+    #h = plotVs1(cos2t.vector(),**plot_settingscos2) # 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('Truecos2ThetaMedc2')
+    
+    #plt.figure(9).clear()
+    #h = plotVs1(cos2.vector(),**plot_settingscos2) # 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('Reccos2ThetaMedc2')
+    
+    #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 (17)
+    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)
+    
+    #Computing the vector field \mathbf{K} for the right hand side in equation (14)
+    K0 = V11[0] - V22[0] + V211
+    K1 = -V11[1] + V22[1] + V210
+        
+    #Computing \mathbf{G} for the right hand side in (14)
+    dlogA0 = Function(Vs1)
+    dlogA1 = Function(Vs1)
+    
+    dlogA0.vector()[:] = np.multiply(np.cos(2*theta1.vector()),K0) - np.multiply(np.sin(2*theta1.vector()),K1)
+    dlogA1.vector()[:] = np.multiply(np.cos(2*theta1.vector()),K1) + np.multiply(np.sin(2*theta1.vector()),K0)
+    
+    logAt = Function(Vs1)
+    logAt.vector()[:] = np.log(sigt.vector())
+    
+    
+    bclogA = DirichletBC(Vs1, logAt, boundary)
+    
+    
+    logA = TrialFunction(Vs1)
+    vA = TestFunction(Vs1)
+    
+    #Defining and solving the variational formulation of the Poisson problem (15)
+    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_settings5 = {
+        'levels': np.linspace(0.2,2.0001,250),
+        'cmap': plt.get_cmap('bwr'),
+    }
+    
+    #Illustrating the true and reconstructed conductivity such as in figure 12
+    plt.figure(10).clear()
+    h = plotVs1(sigt.vector(),**plot_settings5)
+    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('TrueSigma2')
+    
+    plt.figure(11).clear()
+    h = plotVs1(A.vector(),**plot_settings5)
+    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('RecSigmaLarc2')
+    
+    #Computing the relative L2-error
+    print((errornorm(A,sigt,'L2')/norm(sigt,'L2'))*100)
+    
+    
+    
+    
+    
+
+    
+    
+

ContBC/modules.sh

+#!/bin/bash
+
+# Load modules
+#module load FEniCS/2017.2.0-with-petsc-and-slepc
+module load FEniCS/2018.1.0-with-petsc-and-slepc-and-scotch-and-newmpi
+module load scipy/0.19.1-python-3.6.2
+module load matplotlib/2.0.2-python-3.6.2
+module load ffmpeg/3.4
+module load numba/0.35.0-python-3.6.2

ContBCNoise/cmap/__init__.py

+from .cmap import *

ContBCNoise/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))
+    
+    #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.pyplot as plt
+    from matplotlib import ticker
+    plt.ion()
+    
+    # Load 
+    import cmap as cmap
+    from dolfin import __version__ as DOLFINversion
+    from dolfin import TrialFunction, TestFunction, FunctionSpace
+    from dolfin import project, Point, triangle
+    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 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
+    
+    #Mesh to generate the power density data:
+    Ms = 200
+    m = UnitCircleMesh(Ms)
+    Vs = VsSpace(m,1) 
+    Vq = VqSpace(m,1)
+    
+    #Mesh to solve the inverse problem:
+    Ms2 = 160
+    m2 = UnitCircleMesh(Ms2)
+    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
+    N2 = xy2[:,0].size
+      
+    # ------------------------------
+    # 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 = {
+        'cmap': plt.set_cmap('bwr'),
+    }
+
+    
+    
+    # ------------------------------    
+    # Boundary conditions for \Gamma_{medium}
+
+    
+    #Gamma_{medium}:
+    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(2*ufl.atan_2(x[1],x[0]))-1
+            else :
+                value[0] = 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(2*ufl.atan_2(x[1],x[0]))
+            else :
+                value[0] = 0
+
+    
+    u_D1 = MyExpression1(degree=2)
+    u_D2 = MyExpression2(degree=2)
+        
+    def boundary(x, on_boundary):
+        return on_boundary
+    
+    bc1 = DirichletBC(Vs, u_D1, boundary)
+    bc2 = DirichletBC(Vs, u_D2, boundary)
+
+    u1 = Function(Vs)
+    v1 = TestFunction(Vs)
+    u2 = Function(Vs)
+    v2 = TestFunction(Vs)
+    
+    
+    #Defining and solving the variational equations
+    a1 = inner(sig1*grad(u1),grad(v1))*dx
+    a2 = inner(sig1*grad(u2),grad(v2))*dx 
+    
+    solve(a1 == 0,u1,bc1)
+    solve(a2 == 0,u2,bc2)
+    
+    #Defining the gradients
+    dU1 = GradSolver(u1.vector())
+    dU2 = GradSolver(u2.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
+    
+    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]
+    
+    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)
+    
+    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.set_cmap('RdBu'),
+    }
+    
+    mildet = np.amin(lJac.get_local())
+    maldet = np.amax(lJac.get_local())
+    
+    plot_settingsldet = {
+        'levels': np.linspace(mildet,maldet,250),
+        #'levels': np.linspace(0,5,120),
+        #'cmap': plt.set_cmap('RdBu'),
+        'cmap': plt.get_cmap('bwr'),
+    }
+    
+    miu1 = np.amin(u1.vector().get_local())
+    mau1 = np.amax(u1.vector().get_local())
+    
+    plot_settingsu1 = {
+        'levels': np.linspace(miu1,mau1,250),
+        #'levels': np.linspace(0,5,120),
+        #'cmap': plt.set_cmap('RdBu'),
+        'cmap': plt.get_cmap('bwr'),
+    }
+    
+    miu2 = np.amin(u2.vector().get_local())
+    mau2 = np.amax(u2.vector().get_local())
+    
+    plot_settingsu2 = {
+        'levels': np.linspace(miu2,mau2,250),
+        #'levels': np.linspace(0,5,120),
+        #'cmap': plt.set_cmap('RdBu'),
+        'cmap': plt.get_cmap('bwr'),
+    }
+    
+    #Illustrating the solutions u1 and u2
+    plt.figure(1).clear()
+    h = plotVs(u1.vector(),**plot_settingsu1) # 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('u1medc1')
+    
+    plt.figure(2).clear()
+    h = plotVs(u2.vector(),**plot_settingsu2) # 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('u2medc1')
+    
+    #Illustrating log(det(H)):
+    plt.figure(3).clear()
+    h = plotVs1(lJac,**plot_settingsldet) # 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('logDetMedc1') 
+
+    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 (17)
+    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 (13)
+    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))
+
+    bctheta1 = DirichletBC(Vs1, theta1t, boundary)
+    
+    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)
+    
+    miT = np.amin(theta1.vector().get_local())
+    maT = np.amax(theta1.vector().get_local())
+    
+    plot_settingstest = {
+        'levels': np.linspace(miT,maT,250),
+        'cmap': plt.get_cmap('bwr'),
+    }
+    
+    #Illustrating the true vs. the reconstructed \theta-function
+    plt.figure(6).clear()
+    h = plotVs1(theta1t.vector(),**plot_settingstest) # 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('TrueThetaLarc2')
+    
+    plt.figure(7).clear()
+    h = plotVs1(theta1.vector(),**plot_settingstest) # 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('RecThetaLarc2')
+
+    #Compute the relative L2-error for \theta
+    print((errornorm(theta1,theta1t,'L2')/norm(theta1t,'L2'))*100)
+    
+    
+    #Defining the vector fields V11 and V22 defined in equation (17)
+    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)
+    
+    #Computing the vector field \mathbf{K} for the right hand side in equation (14)
+    K0 = V11[0] - V22[0] + V211
+    K1 = -V11[1] + V22[1] + V210
+        
+    #Computing \mathbf{G} for the right hand side in (14)
+    dlogA0 = Function(Vs1)
+    dlogA1 = Function(Vs1)
+    
+    dlogA0.vector()[:] = np.multiply(np.cos(2*theta1.vector()),K0) - np.multiply(np.sin(2*theta1.vector()),K1)
+    dlogA1.vector()[:] = np.multiply(np.cos(2*theta1.vector()),K1) + np.multiply(np.sin(2*theta1.vector()),K0)
+    
+    logAt = Function(Vs1)
+    logAt.vector()[:] = np.log(sigt.vector())
+    
+    
+    bclogA = DirichletBC(Vs1, logAt, boundary)
+    
+    
+    logA = TrialFunction(Vs1)
+    vA = TestFunction(Vs1)
+    
+    #Defining and solving the variational formulation of the Poisson problem (15)
+    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_settings5 = {
+        'levels': np.linspace(0,2.1,250),
+        'cmap': plt.get_cmap('bwr'),
+    }
+    
+    #Illustrating the true and reconstructed conductivity such as in figure 12
+    plt.figure(9).clear()
+    h = plotVs1(sigt.vector(),**plot_settings5)
+    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('TrueSigma2')
+    
+    plt.figure(10).clear()
+    h = plotVs1(A.vector(),**plot_settings5)
+    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('RecSigmaLarc2')
+    
+    #Computing the relative L2-error
+    print((errornorm(A,sigt,'L2')/norm(sigt,'L2'))*100)
+    
+    
+    
+    
+    
+
+    
+    
+

ContBCNoise/modules.sh

+#!/bin/bash
+
+# Load modules
+#module load FEniCS/2017.2.0-with-petsc-and-slepc
+module load FEniCS/2018.1.0-with-petsc-and-slepc-and-scotch-and-newmpi
+module load scipy/0.19.1-python-3.6.2
+module load matplotlib/2.0.2-python-3.6.2
+module load ffmpeg/3.4
+module load numba/0.35.0-python-3.6.2

DiscontBC/cmap/__init__.py

+from .cmap import *

DiscontBC/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))
+    
+    #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.pyplot as plt
+    from matplotlib import ticker
+    plt.ion()
+    
+    # Load 
+    import cmap as cmap
+    from dolfin import __version__ as DOLFINversion
+    from dolfin import TrialFunction, TestFunction, FunctionSpace
+    from dolfin import project, Point, triangle
+    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 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 = 200
+    m = UnitCircleMesh(Ms)
+    Vs = VsSpace(m,1) 
+    Vq = VqSpace(m,1)
+    
+    #Mesh to solve the inverse problem:
+    Ms2 = 160
+    m2 = UnitCircleMesh(Ms2)
+    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
+    N2 = xy2[:,0].size
+      
+    # ------------------------------
+    # 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.set_cmap('RdBu'),
+    }
+
+    
+    
+    # ------------------------------    
+    #Load source terms that are generated by solving the Laplace equation (18) semi-analytically in Matlab using code by Nuutti Hyvönen
+    source1 = Function(Vs)
+    source2 = Function(Vs)
+    
+    #For \Gamma_{large}
+    #source1tmp = io.loadmat('sxyf1larTest3.mat')['sxy']
+    #source2tmp = io.loadmat('sxyf2larTest3.mat')['sxy']
+    
+    #For \Gamma_{medium}
+    source1tmp = io.loadmat('sxyf1medTest3.mat')['sxy']
+    source2tmp = io.loadmat('sxyf2medTest3.mat')['sxy']
+    
+    #For \Gamma_{small}
+    #source1tmp = io.loadmat('sxyf1miniTest3.mat')['sxy']
+    #source2tmp = io.loadmat('sxyf2miniTest3.mat')['sxy']
+    
+    
+    
+    source1tmp = np.array(source1tmp)
+    source2tmp = np.array(source2tmp)
+
+    source1.vector().set_local(source1tmp.flatten())
+    source2.vector().set_local(source2tmp.flatten())
+        
+    def boundary(x, on_boundary):
+        return on_boundary
+    
+    bc1 = DirichletBC(Vs, Constant(0.0), boundary)
+    bc2 = DirichletBC(Vs, Constant(0.0), boundary)
+    
+    
+    w1 = TrialFunction(Vs)
+    v1 = TestFunction(Vs)
+    w2 = TrialFunction(Vs)
+    v2 = TestFunction(Vs)
+
+    a1 = inner(sig1*grad(w1),grad(v1))*dx
+    a2 = inner(sig1*grad(w2),grad(v2))*dx
+    
+    l1 = -inner(sig1*grad(source1),grad(v1))*dx
+    l2 = -inner(sig1*grad(source2),grad(v2))*dx
+    
+    w1 = Function(Vs)
+    w2 = Function(Vs)
+    solve(a1 == l1,w1,bc1)
+    solve(a2 == l2,w2,bc2)
+    
+    u1 = Function(Vs)
+    u2 = Function(Vs)
+    
+    u1.vector()[:] = w1.vector() + source1.vector()
+    u2.vector()[:] = w2.vector() + source2.vector()
+    
+    #Defining the gradients
+    dU1 = GradSolver(u1.vector())
+    dU2 = GradSolver(u2.vector())
+
+    H11t = Function(Vs)
+    H12t = Function(Vs)
+    H22t = Function(Vs)
+    
+    #Compute the power density data
+    H11t.vector()[:]  = sig1.vector()*(dU1[0]*dU1[0]+dU1[1]*dU1[1])
+    H12t.vector()[:]  = sig1.vector()*(dU1[0]*dU2[0]+dU1[1]*dU2[1])
+    H22t.vector()[:]  = sig1.vector()*(dU2[0]*dU2[0]+dU2[1]*dU2[1])
+    
+    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 data to the mesh for solving the inverse problem
+    H11 = project(H11t,Vs1)
+    H12 = project(H12t,Vs1)
+    H22 = project(H22t,Vs1)
+    
+    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)
+    
+    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.set_cmap('RdBu'),
+    }
+    
+    mildet = np.amin(lJac.get_local())
+    maldet = np.amax(lJac.get_local())
+    
+    plot_settingsldet = {
+        'levels': np.linspace(mildet,maldet,250),
+        #'levels': np.linspace(0,5,120),
+        #'cmap': plt.set_cmap('RdBu'),
+        'cmap': plt.get_cmap('bwr'),
+    }
+    
+    miu1 = np.amin(u1.vector().get_local())
+    mau1 = np.amax(u1.vector().get_local())
+    
+    plot_settingsu1 = {
+        'levels': np.linspace(miu1,mau1,250),
+        #'levels': np.linspace(0,5,120),
+        #'cmap': plt.set_cmap('RdBu'),
+        'cmap': plt.get_cmap('bwr'),
+    }
+    
+    miu2 = np.amin(u2.vector().get_local())
+    mau2 = np.amax(u2.vector().get_local())
+    
+    plot_settingsu2 = {
+        'levels': np.linspace(miu2,mau2,250),
+        #'levels': np.linspace(0,5,120),
+        #'cmap': plt.set_cmap('RdBu'),
+        'cmap': plt.get_cmap('bwr'),
+    }
+    
+    #Illustrating the solutions u1 and u2
+    plt.figure(1).clear()
+    h = plotVs(u1.vector(),**plot_settingsu1) # 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('u1medc1')
+    
+    plt.figure(2).clear()
+    h = plotVs(u2.vector(),**plot_settingsu2) # 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('u2medc1')
+    
+    #Illustrating log(det(H)):
+    plt.figure(3).clear()
+    h = plotVs1(lJac,**plot_settingsldet) # 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('logDetMedc1') 
+
+    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 (17)
+    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 (13)
+    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 when using \Gamma_{small} defined in equation (20):
+    
+    #indBdr1tmp = np.where((np.arctan2(xy2[:,1],xy2[:,0])<np.pi) & (np.arctan2(xy2[:,1],xy2[:,0])>np.pi/4.0+0.001))
+    
+    #indBdr1 = np.asarray(indBdr1tmp)
+
+    #indBdr1 = np.reshape(indBdr1,indBdr1.shape[1])
+
+    
+    #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).clear()
+    #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")
+    
+    plot_settingsT = {
+        'levels': np.linspace(-np.pi,0,250),
+        #'levels': np.linspace(0,5,120),
+        'cmap': cmaps['twilight_shifted']
+    }
+    
+    #plt.figure(5).clear()
+    #h = plotVs1(theta1t.vector(),**plot_settingsT) # 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, boundary)
+    
+    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)
+    
+    miT = np.amin(theta1.vector().get_local())
+    maT = np.amax(theta1.vector().get_local())
+    
+    plot_settingstest = {
+        'levels': np.linspace(miT,maT,250),
+        'cmap': plt.get_cmap('bwr'),
+    }
+    
+    #Illustrating the true vs. the reconstructed \theta-function
+    plt.figure(6).clear()
+    h = plotVs1(theta1t.vector(),**plot_settingstest) # 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('TrueThetaLarc2')
+    
+    plt.figure(7).clear()
+    h = plotVs1(theta1.vector(),**plot_settingstest) # 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('RecThetaLarc2')
+    
+    #Compute the relative L2-error for \theta
+    print((errornorm(theta1,theta1t,'L2')/norm(theta1t,'L2'))*100)
+    
+    #Calculating the relative errors for the reconstruction of \cos(2\theta) and \sin(2\theta) in the case of \Gamma_{small}:
+    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())
+    
+    mic2 = np.amin(cos2.vector().get_local())
+    mac2 = np.amax(cos2.vector().get_local())
+    
+    plot_settingscos2 = {
+        'levels': np.linspace(mic2,mac2,250),
+        'cmap': plt.get_cmap('bwr'),
+    }
+    
+    #Illustrations of cos(2\theta) and its reconstructed version in the case of \Gamma_{small}
+    #plt.figure(8).clear()
+    #h = plotVs1(cos2t.vector(),**plot_settingscos2) # 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('Truecos2ThetaMedc2')
+    
+    #plt.figure(9).clear()
+    #h = plotVs1(cos2.vector(),**plot_settingscos2) # 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('Reccos2ThetaMedc2')
+    
+    #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 (17)
+    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)
+    
+    #Computing the vector field \mathbf{K} for the right hand side in equation (14)
+    K0 = V11[0] - V22[0] + V211
+    K1 = -V11[1] + V22[1] + V210
+        
+    #Computing \mathbf{G} for the right hand side in (14)
+    dlogA0 = Function(Vs1)
+    dlogA1 = Function(Vs1)
+    
+    dlogA0.vector()[:] = np.multiply(np.cos(2*theta1.vector()),K0) - np.multiply(np.sin(2*theta1.vector()),K1)
+    dlogA1.vector()[:] = np.multiply(np.cos(2*theta1.vector()),K1) + np.multiply(np.sin(2*theta1.vector()),K0)
+    
+    logAt = Function(Vs1)
+    logAt.vector()[:] = np.log(sigt.vector())
+    
+    
+    bclogA = DirichletBC(Vs1, logAt, boundary)
+    
+    
+    logA = TrialFunction(Vs1)
+    vA = TestFunction(Vs1)
+    
+    #Defining and solving the variational formulation of the Poisson problem (15)
+    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())
+    
+    miA = np.amin(A.vector().get_local())
+    maA = np.amax(A.vector().get_local())
+    
+    plot_settings4 = {
+        'levels': np.linspace(0,2.0001,250),
+        'cmap': plt.get_cmap('RdBu'),
+    }
+    
+    plot_settings5 = {
+        'levels': np.linspace(0.2,2.0001,250),
+        'cmap': plt.get_cmap('bwr'),
+    }
+    
+    #Illustrating the true and reconstructed conductivity such as in figure 12
+    plt.figure(10).clear()
+    h = plotVs1(sigt.vector(),**plot_settings5)
+    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('TrueSigma2')
+    
+    plt.figure(11).clear()
+    h = plotVs1(A.vector(),**plot_settings5)
+    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('RecSigmaLarc2')
+    
+    #Computing the relative L2-error
+    print((errornorm(A,sigt,'L2')/norm(sigt,'L2'))*100)
+    
+    
+    
+    
+    
+
+    
+    
+