### Strange error with boundary conditions

269
views
0
12 months ago by
Hi,
I have modified the example on how to implement multiple boundary conditions from the Fenics tutorial (p. 95-) and I am getting a strange error: The following code implements  $-\Delta u=0$Δu=0  on the unit square with the following boundary conditions:  $u\left(x,y\right)=x$u(x,y)=x on the boundary. The analytical solution of this is a nice sloped plane. I have chosen to split up the boundary in four parts (because I actually want to have different bcs later on), so the code looks like this:
from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
from fenics import *

resol=4

mesh = RectangleMesh(Point(0,0), Point(1,1), 2**resol, 2**resol)
V = FunctionSpace(mesh, 'P', 1)

f = Constant(0)

boundary_markers = FacetFunction("size_t", mesh)

class BoundaryY0(SubDomain):
tol = 1E-14
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 0, self.tol)
class BoundaryX1(SubDomain):
tol = 1E-14
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 1, self.tol)
class BoundaryY1(SubDomain):
tol = 1E-14
def inside(self, x, on_boundary):
return on_boundary and near(x[1], 1, self.tol)

class BoundaryX0(SubDomain):
tol = 1E-14
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0, self.tol)

bx0 = BoundaryX0()
bx0.mark(boundary_markers, 0)
by0 = BoundaryY0()
by0.mark(boundary_markers, 1)
bx1 = BoundaryX1()
bx1.mark(boundary_markers, 2)
by1 = BoundaryY1()
by1.mark(boundary_markers, 3)

boundary_conditions = {0: {'Dirichlet': Constant(0.0)}, 1: {'Dirichlet':   Expression("x[0]", degree=2)}, 2: {'Dirichlet': Constant(1)}, 3: {'Dirichlet':   Expression("x[0]", degree=2)}}

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)

L = f*v*dx # = 0

uSol = Function(V)
solve(a == L, uSol, bcs)
plot(uSol)
interactive()​
I have more or less exactly followed the tutorial and only adapted it to my case. Now the result is completely wrong:
I'd like to upload a plot of the obtained solution but I can't see how, so: The solution is exactly zero everywhere but on the boundary, where the bc is enforced (so the function is discontinuous).

Now if I change the code only very slightly by numbering the boundary conditions differently, i.e. skipping 0 and marking BoundaryX0 with a "4" instead, everything works (just change the following part in the code):
bx0 = BoundaryX0()
bx0.mark(boundary_markers, 4)
by0 = BoundaryY0()
by0.mark(boundary_markers, 1)
bx1 = BoundaryX1()
bx1.mark(boundary_markers, 2)
by1 = BoundaryY1()
by1.mark(boundary_markers, 3)

boundary_conditions = {4: {'Dirichlet': Constant(0.0)}, 1: {'Dirichlet':   Expression("x[0]", degree=2)}, 2: {'Dirichlet': Constant(1)}, 3: {'Dirichlet':   Expression("x[0]", degree=2)}}​
The solution now looks exactly like it should, i.e. it's a sloped plane. Can anyone explain the reason for this? (I'm actually interested in more complicated bcs, with Neumann instead of Dirichlet on two of the four boundaries and here I have similar problems but I couldn't fix them by re-numbering the bcs so I need some better understanding of the mechanics here)
Community: FEniCS Project

3
12 months ago by
When you initialise a new FacetFunction of type sizeT it marks all the facets with 0..

When you then create the corresponding DirichletBC for the region marked 0. you end up applying the boundary condition to all the facets in the domain.

This first boundary condition therefore sets all interior degrees of freedom to zero. When you apply the subsequent boundary conditions they change some of the DoF on the boundaries back to the values you prescribed.
I feel like this is something that is worth mentioning in the tutorial. Thanks!
written 12 months ago by CalabiYau