Hello All,
Recently, I encounter an accuracy problem for the Poisson equation with subdomains. Here I use the question in this post (Solving PDEs in Python - <br> The FEniCS Tutorial Volume I) for test and explanation. the unit square is divided by two subdomains. and the top and bottom boundaries are biased at certain voltages by DirichletBC (top: 1V, bottom: 0V). the permittivity for upper and lower materials is 6 and 2, respectively. the Poisson equation is sourceless. Theoretically, the potential line at y=0.5 is constantly 0.75V.
I find that the result for potential at y=0.5 is sensitive to the way generating mesh. For example, if I use the built-in mesh “UnitSquareMesh(10,10)” to generate the mesh, the potential at y=0.5 is exactly 0.75.
x potential
0.4 0.75
0.5 0.75
0.6 0.75
However, if I use below code to generate mesh,
geometry=Rectangle(Point(0,0),Point(1,1))
mesh=generate_mesh(geometry,10,'cgal')
the result has significant difference to the exact result.
x potential
0.4 0.7299
0.5 0.7463
0.6 0.7379
I am noticed that in fenics, the accuracy for vertices in the mesh is 10^(−15), while the accuracy for non-vertex points is 10^(−2) since interpolation is used for non-vertex points. Thus I was trying to move the points in the mesh to the points that would be calculated.
# mesh
geometry=Rectangle(Point(0,0),Point(1,1))
mesh=generate_mesh(geometry,10,'cgal')
vtkfile=File('data_mesh_before.pvd')
vtkfile<<mesh
calc_x=np.linspace(0.4,0.6,3)
# moving points in mesh
tree=BoundingBoxTree()
tree.build(mesh)
cell_index=mesh.cells()
vertex_index=mesh.coordinates()
for x in calc_x:
p_flake=Point(x,0.5)
first=tree.compute_first_entity_collision(p_flake)
vertex_index[cell_index[first][0]]=[p_flake.x(),p_flake.y()]
vtkfile=File('data_mesh_after.pvd')
vtkfile<<mesh
The mesh now can include the points to be calculated. But it still doesn’t work for the potential at y=0.5.
x potential
0.4 0.7267
0.5 0.7475
0.6 0.7477
I will extend this 2D problem to 3D thus I prefer to use something like “generate_mesh(geometry,10,‘cgal’)” to generate mesh. May I ask if anyone can shed some light on this issue? Thanks. Below I attach the code without moving points in mesh.
from fenics import *
from mshr import *
import numpy as np
# mesh
geometry=Rectangle(Point(0,0),Point(1,1))
mesh=generate_mesh(geometry,10,'cgal')
#mesh=UnitSquareMesh(10,10)
vtkfile=File('data_mesh.pvd')
vtkfile<<mesh
tol=1E-14
# boundary conditions
def Gate1(x,on_boundary):
return on_boundary and near(x[1],0,tol)
def Gate2(x,on_boundary):
return on_boundary and near(x[1],1,tol)
class hBN(SubDomain):
def inside(self,x,on_boundary):
return between(x[1],(0.0-tol,0.5+tol))
class SiO2(SubDomain):
def inside(self,x,on_boundary):
return between(x[1],(0.5-tol,1.0+tol))
hbn=hBN()
sio2=SiO2()
# define subdomain markers and integration measure
markers=MeshFunction('size_t',mesh,mesh.topology().dim())
dx=Measure('dx',domain=mesh,subdomain_data=markers)
hbn.mark(markers,0)
sio2.mark(markers,1)
# define permittivity for different domains
class Permittivity(UserExpression):
def __init__(self,markers,**kwargs):
super().__init__(markers,**kwargs)
self.markers=markers
def eval_cell(self,values,x,cell):
if self.markers[cell.index]==0:
values[0]=2 # hBN
elif self.markers[cell.index]==1:
values[0]=6 # SiO2
def value_shape(self):
return ()
er=Permittivity(markers,degree=0)
# initial voltages for gate 1 and 3
# Define variational problem
V=FunctionSpace(mesh,'CG',1)
u=TrialFunction(V)
v=TestFunction(V)
aint=er*dot(grad(u),grad(v))*dx
Lint=Constant(0)*v*dx
bc1=DirichletBC(V,0,Gate1)
bc2=DirichletBC(V,1,Gate2)
bcs=[bc1,bc2]
A,b=assemble_system(aint,Lint,bcs)
u=Function(V)
solve(A,u.vector(),b)
vtkfile=File('potential.pvd')
vtkfile<<u
# calculate potential profiles
calc_x=np.linspace(0.4,0.6,3)
points=[(x_,0.5) for x_ in calc_x]
p_line=np.array([[point[0],u(point)] for point in points])
np.savetxt('data_p_line.txt',p_line)