Neumann boundary condition not fulfilled


122
views
0
4 months ago by
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

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ƒ tu  over the Dirichlet boundary and returns 0 (as it should). The last line evaluates  $\int_{\Gamma_{Neumann}}\partial_nu$ΓNeumannnu 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$28×28 ) too coarse? But even  $2^9\times2^9$29×29 only reduces this slightly.
Community: FEniCS Project

1 Answer


4
4 months ago by
pf4d  
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.

Similar posts:
Search »