diff --git a/CutOffBoundaryFunsNoise/Gamma5/NoiseLevel10/main.py b/CutOffBoundaryFunsNoise/Gamma5/NoiseLevel10/main.py
new file mode 100644
index 0000000000000000000000000000000000000000..16d3a86a55c41bf58db47ac65ea28487f647f40d
--- /dev/null
+++ b/CutOffBoundaryFunsNoise/Gamma5/NoiseLevel10/main.py
@@ -0,0 +1,887 @@
+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('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())
+
+
+
+
+
+
+
+
+
+