### Neumann boundary condition not fulfilled

122

views

0

Hi,

I'm solving the Poisson problem on the unit square with the following boundary conditions: u = 0 on the "left" boundary (i.e. x=0), u(x,y) = x+y on the "lower right" boundary (y=0 and (x=1, 0 <= y <= 1/2)) and Neumann-0 boundary conditions on the rest ("upper right"). Problem is that the solution does not really enforce the Neumann boundary conditions. See the following working example

I'm solving the Poisson problem on the unit square with the following boundary conditions: u = 0 on the "left" boundary (i.e. x=0), u(x,y) = x+y on the "lower right" boundary (y=0 and (x=1, 0 <= y <= 1/2)) and Neumann-0 boundary conditions on the rest ("upper right"). Problem is that the solution does not really enforce the Neumann boundary conditions. See the following working example

```
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
from fenics import *
resol=8
mesh = RectangleMesh(Point(0,0), Point(1,1), 2**resol, 2**resol)
V = FunctionSpace(mesh, 'P', 1)
f = Expression("pow(x[0]-0.4, 2) + pow(x[1]-0.3, 2) <= 0.05 ? -10 : 0.0", degree=5)
boundary_markers = FacetFunction("size_t", mesh)
def boundaryclassifier(x):
tol = 1E-14
if near(x[0], 0, tol):
return "DirichletLeft"
elif near(x[1], 0, tol) or (near(x[0], 1, tol) and x[1] <= 0.5):
return "DirichletRightBottom"
elif near(x[1], 1, tol) or (near(x[0], 1, tol) and x[1] > 0.5):
return "NeumannRightTop"
else:
return None
class BoundaryLeft(SubDomain):
tol = 1E-14
def inside(self, x, on_boundary):
return on_boundary and boundaryclassifier(x) == "DirichletLeft"
class BoundaryRB(SubDomain):
tol = 1E-14
def inside(self, x, on_boundary):
return on_boundary and boundaryclassifier(x) == "DirichletRightBottom"
class BoundaryN(SubDomain):
tol = 1E-14
def inside(self, x, on_boundary):
return on_boundary and boundaryclassifier(x) == "NeumannRightTop"
bl = BoundaryLeft()
bl.mark(boundary_markers, 1)
brb = BoundaryRB()
brb.mark(boundary_markers, 2)
bn = BoundaryN()
bn.mark(boundary_markers, 3)
boundary_conditions = {1: {'Dirichlet': Constant(0.0)}, 2: {'Dirichlet': Expression("x[0]+x[1]", degree=5)}, 3: {'Neumann': Constant(0.0)}}
bcs = []
for i in boundary_conditions:
if 'Dirichlet' in boundary_conditions[i]:
bc = DirichletBC(V, boundary_conditions[i]['Dirichlet'], boundary_markers, i)
bcs.append(bc)
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v))*dx
L = f*v*dx
uSol = Function(V)
solve(a == L, uSol, bcs)
plot(uSol)
boundary_conditions2 = {1: {'Dirichlet': Constant(0.0)}, 2: {'Dirichlet': Constant(0.0)}, 3: {'Neumann': Constant(0.0)}}
ds_surf = Measure("ds", domain=mesh, subdomain_data=boundary_markers)
n = FacetNormal(mesh)
print(assemble(uSol*ds_surf(1)))
print(assemble(dot(grad(uSol),n)*ds_surf(3)))
```

The last two lines are the diagnostic I am interested in: The second-to-last line evaluates the boundary integral $\int_{\Gamma_{left}}u$∫_{Γleƒ t}`u`over the Dirichlet boundary and returns 0 (as it should). The last line evaluates $\int_{\Gamma_{Neumann}}\partial_nu$∫_{ΓNeumann}∂_{n}`u`but returns a number which is significantly not zero (something like 0.05). What's the reason for this bad compliance with the boundary condition at this point? Is my grid ( $2^8\times2^8$2^{8}×2^{8}) too coarse? But even $2^9\times2^9$2^{9}×2^{9}only reduces this slightly.
Community: FEniCS Project

### 1 Answer

4

There is interpolation error associated with the P1 function space whose first derivative is constant over cells. Using the quadratic P2 space, the error reduces to ~ -0.0038.

Please login to add an answer/comment or follow this question.