Hi everyone,
so I've struggled with this problem for a few days now and can't find the answer. I've worked out a simple example to present to you.
The domain $\Omega$ is a square $[-1,1]\times [-1,1]$, divided in two subdomains $\Omega_1 = {(x,y): x^2+y^2 < 0.4}$ and $\Omega_2$ is the complementary of $\Omega_1$.
The PDE is of the type
$\operatorname{div}\sigma(x,y)\nabla u(x,y) = 0 $
where $\sigma$ is piecewise constant, equal to 1 in $\Omega_1$ and to $k=0.01$ in $\Omega_2$.
I set Dirichlet boundary conditions on the left side of the square, and Neumann conditions $\partial_n u=1$ on the right side. On the two other sides we have the Neumann conditions $\partial_n u=0$.
I use RectangleMesh for the grid so I expect the solution to be perfectly symmetric since the data is symmetric. However I get a quite large error (in terms of symmetry I mean), which decreases as the mesh gets finer, but for my application it is quite important to have a perfectly symmetric solution for any size of the grid.
My suspicion is that the asymmetry appears because the domain $\Omega_1$ is not aligned with the grid (it's a disk). The problem is I don't know what Fenics is doing with the subdomains exactly so I'm kind of stuck.
Any idea on how I could get a symmetric solution? Should I change the way the subdomains are defined? Below you can find the code. It's printing the maximum difference between the solution U and its mirror image with respect to the x-axis.
from dolfin import *
import numpy as np
import pdb
#------------------------------
Nx,Ny,lx,ly = [50,50,1.0,1.0]
k = 0.01
# mesh and function spaces
mesh = RectangleMesh(-lx,-ly,lx,ly, Nx, Ny)
V = FunctionSpace(mesh, 'CG', 1)
phi = Expression('x[1]*x[1] + x[0]*x[0] -0.4')
# Define domains and boundaries
class Omega1(SubDomain):
def inside(self, x, on_boundary):
return between(x[0], (-lx, lx)) and between(x[1], (-ly, ly)) and phi(x) < 0
class DirBd(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], -lx)
class DirNeu(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], lx)
omega1 = Omega1()
dirBd,dirNeu = [DirBd(),DirNeu()]
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
dirBd.mark( boundaries, 1)
dirNeu.mark( boundaries, 2)
domains = CellFunction("size_t", mesh)
domains.set_all(0)
omega1.mark(domains, 1)
dx = Measure('dx')[domains]
ds = Measure("ds")[boundaries]
# Solve equation
bcd = DirichletBC(V, (0.0), boundaries, 1)
u,v = [TrialFunction(V), TestFunction(V)]
a = inner(grad(u),grad(v))*dx(1) + k*inner(grad(u),grad(v))*dx(0)
R = v*ds(2)
U,A,b = [Function(V),assemble(a),assemble(R)]
bcd.apply(A)
bcd.apply(b)
solver = LUSolver(A)
solver.solve(U.vector(), b)
# Compute the maximum difference between the solution U and its mirror image U_flip
U_array = U.vector().array()
U_vec = U_array[V.dofmap().dof_to_vertex_map(mesh)]
U_mat = U_vec.reshape( Ny+1, Nx+1)
U_flip = np.flipud(U_mat)
symm_diff = U_flip - U_mat
print("Maximal symmetry error : %s " % symm_diff.max())
print("Maximal value of U : %s " % U_mat.max())
pdb.set_trace()
