I would like to compare the outward flux through various parts of the boundary of the exact and numerical solutions to the following boundary-value problem:
\begin{equation}
\Delta u = 0 \text{ in } \Omega, \quad \int_\Omega u \, dx = 0,
\end{equation}
\begin{equation}
\partial_n u = 0 \text{ on } \Gamma_0, \quad \partial_n u = 1 \text{ on } \Gamma_1, \quad \partial_n u = 0 \text{ on } \Gamma_2, \quad \partial_n u = -1 \text{ on } \Gamma_3,
\end{equation}
where $\Omega := [0,1]\times [0,1]$ and
\begin{equation}
\Gamma_0 := { (x,0): 0 < x < 1 },
\end{equation}
\begin{equation}
\Gamma_1 := { (1, y): 0 < y < 1 },
\end{equation}
\begin{equation}
\Gamma_2 := { (x,1): 0 < x < 1 },
\end{equation}
\begin{equation}
\Gamma_3 := { (y, 0): 0 < x < 1 }.
\end{equation}
The exact solution is $u_\text{exact}(x,y) = x - \frac{1}{2}$ and the associated outward-flux quantities are
\begin{equation}
\int_{\Gamma_0} \nabla u_\text{exact} \cdot n \, ds = 0, \ \int_{\Gamma_1} \nabla u_\text{exact} \cdot n \, ds = 1,
\end{equation}
\begin{equation}
\int_{\Gamma_2} \nabla u_\text{exact} \cdot n \, ds = 0, \ \int_{\Gamma_3} \nabla u_\text{exact} \cdot n \, ds = -1.
\end{equation}
In order to obtain the numerical solution $u$ and the associated outward-flux quantities, I use the script given at the end of this post.
The problem is that the script (at the end of this post) obtains the following values for the outward-flux quantities:
\begin{equation}
\int_{\Gamma_i} \nabla u \cdot n \, ds = -1.69655955951 * 10^{-15}, \quad i = 0, 1, 2, 3,
\end{equation}
which is correct for $i = 0, 2$ but incorrect for $i = 1, 3$.
May I ask for advice on how to fix my script in order to obtain reasonable numerical approximations to the exact values of all four outward-flux quantities?
from dolfin import *
# ------------------------
# BEGIN: class definitions
class LowerBoundaryOfUnitSquare(SubDomain):
""" This class represents and manipulates
the lower boundary of the unit square
[0,1]x[0,1]. """
def inside(self, x, on_boundary):
tol = 1E-14 # tolerance for coordinate comparisons
return on_boundary and abs(x[1]) < tol
class UpperBoundaryOfUnitSquare(SubDomain):
""" This class represents and manipulates
the upper boundary of the unit square
[0,1]x[0,1]. """
def inside(self, x, on_boundary):
tol = 1E-14 # tolerance for coordinate comparisons
return on_boundary and abs(x[1] - 1) < tol
class LeftBoundaryOfUnitSquare(SubDomain):
""" This class represents and manipulates
the left boundary of the unit square
[0,1]x[0,1]. """
def inside(self, x, on_boundary):
tol = 1E-14 # tolerance for coordinate comparisons
return on_boundary and abs(x[0]) < tol
class RightBoundaryOfUnitSquare(SubDomain):
""" This class represents and manipulates
the right boundary of the unit square
[0,1]x[0,1]. """
def inside(self, x, on_boundary):
tol = 1E-14 # tolerance for coordinate comparisons
return on_boundary and abs(x[0] - 1) < tol
# END: class definitions
# ----------------------
# --------------------
# BEGIN: Main function
nx = 64
ny = 64
mesh = UnitSquareMesh(nx, ny)
# create a mesh function over cell facets
mesh_fn_marking_bdry_parts = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
Gamma_0 = LowerBoundaryOfUnitSquare()
Gamma_0.mark(mesh_fn_marking_bdry_parts, 0)
Gamma_1 = RightBoundaryOfUnitSquare()
Gamma_1.mark(mesh_fn_marking_bdry_parts, 1)
Gamma_2 = UpperBoundaryOfUnitSquare()
Gamma_2.mark(mesh_fn_marking_bdry_parts, 2)
Gamma_3 = LeftBoundaryOfUnitSquare()
Gamma_3.mark(mesh_fn_marking_bdry_parts, 3)
# set up and solve the weak formulation of the Neumann
# boundary-value problem for the Laplacian
V = FunctionSpace(mesh, "CG", 1)
R = FunctionSpace(mesh, "R", 0)
W = V * R
(u, c) = TrialFunction(W)
(v, d) = TestFunction(W)
f = Constant("0")
g_0 = Constant("0")
g_1 = Constant("1")
g_2 = Constant("0")
g_3 = Constant("-1")
(u, c) = TrialFunction(W)
(v, d) = TestFunction(W)
a = ( inner(grad(u),grad(v)) + c*v + d*u )*dx
L = f*v*dx + g_0*v*ds(0) + g_1*v*ds(1) + g_2*v*ds(2) + g_3*v*ds(3)
A = assemble(a, exterior_facet_domains = mesh_fn_marking_bdry_parts)
b = assemble(L, exterior_facet_domains = mesh_fn_marking_bdry_parts)
w = Function(W)
solve(A, w.vector(), b, 'lu')
(u, c) = w.split()
# plot the numerically obtained solution
plot(u, title = 'The solution u')
interactive()
# compute and print the flux through the boundary
unitNormal = FacetNormal(mesh)
for i in range(0,4):
flux = dot(grad(u), unitNormal) * ds(i)(mesh)
flux = assemble(flux)
print "outward flux of u through Gamma_{0}: ".format(i), flux
# END: Main function
# --------------------------
