Select Git revision
hjsc authored
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
main.py 22.72 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)
# Defining \Gamma
class boundaryN(SubDomain):
#\Gamma_{mini}:
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) )
#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)))
#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_{large}:
#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)))
#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())
H11t = Function(Vs)
H12t = Function(Vs)
H22t = Function(Vs)
#Compute the noise free 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 noisy data to the mesh for solving the inverse problem
H11 = project(H11t,Vs1)
H12 = project(H12t,Vs1)
H22 = project(H22t,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.0,4.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('AdaptMini.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('AdaptMini.mat',{'xy2':xy2,
'AdaptMedRec':A.vector().get_local(),
'AdaptTrue': sigt.vector().get_local()})
np.savez('AdaptMini.npz',xy2=xy2,A=A.vector().get_local(),sigt=sigt.vector().get_local())