DOLFIN

Imposing strong Dirichlet boundary condition for discontinuous Galerkin method in one dimension

Asked by Xiaoxian Liu

Hi all,

I ran into trouble when I was trying to solve a simple poisson equation "-Laplacian(u) = f" using discontinuous galerkin method. An example code in Python is copied below:
"...
mesh = UnitInterval(2**6)
nu = FacetNormal(mesh)
h = CellSize(mesh)
h_avg = (h('+') + h('-'))/2.0
V = FunctionSpace(mesh, 'DG', 1)
class Boundary(SubDomain):
 def inside(self, x, on_boundary):
                return on_boundary
bc = DirichletBC(V, Constant(0), Boundary(), 'geometric')

alpha=4.0
u = TrialFunction(V)
v = TestFunction(V)
a = dot(grad(u), grad(v)) *dx\
 - dot(jump(u,nu), avg(grad(v))) *dS\
 - dot(avg(grad(u)), jump(v,nu)) *dS\
 + alpha/h_avg * dot(jump(u,nu), jump(v,nu))*dS
f = Expression('pow(x[0], 2) + pow(x[0]-1,2)')
L = f*v*dx
u = Function(V)
solve(a == L, u, bc)
..."

The error message is :
*** Error: Unable to determine if given point is on facet.
*** Reason: Not implemented for given facet dimension.
*** Where: This error was encountered inside DirichletBC.cpp.

I guess it's not allowed to impose strong Dirichlet boundary condition for DG in one dimension? For this problem, I can impose the boundary function weakly. However, if the Dirichlet boundary is not homogeneous, say "0" on the left end and "1" on the right end, I will have to impose some sort of strong boundary condition. Can it be done or avoided at all in 1D?

Thank you very much!
Xiaoxian

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Jan Blechta
Solved:
Last query:
Last reply:
Revision history for this message
Jan Blechta (blechta) said : 		#1

method='pointwise' works in 1D.
on_boundary does not work with pointwise.
So use following. Works with development version as well as 1.0.0.

class Boundary(SubDomain):
  def inside(self, x, on_boundary):
                return x[0]<DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS

bc = DirichletBC(V, Constant(0), Boundary(), 'pointwise')

Revision history for this message
Xiaoxian Liu (liuxiaox) said : 		#2

Thanks for the reply, Jan! That did work in this problem.
However, I have a further related question: what if we have the following boundary condition:
"u=0 at x=0
u=1 at x=3"

I tried the following lines to generate the boundary condition:

" ...
V = FunctionSpace(mesh, 'DG', 1)
class LeftBoundary(SubDomain):
 def inside(self, x, on_boundary):
  return abs(x[0])<DOLFIN_EPS
class RightBoundary(SubDomain):
 def inside(self, x, on_boundary):
  return abs(x[0]-1)<DOLFIN_EPS

boundary_parts.set_all(2)
left = LeftBoundary()
left.mark(boundary_parts, 0)
right = RightBoundary()
right.mark(boundary_parts, 1)

u0 = Constant(0.0); u1 = Constant(3.0)
bcs = [DirichletBC(V, u0, boundary_parts, 0, 'pointwise'), DirichletBC(V, u0, boundary_parts, 1, 'pointwise')]
..."""

When I tried to recover the boundary values like below:
"""...
bcs[0].get_boundary_values()
"""

The program was terminated and exited Python. Can you help me?

Revision history for this message
Best Jan Blechta (blechta) said : 		#4

> I tried the following lines to generate the boundary condition:
>
> " ...
> V = FunctionSpace(mesh, 'DG', 1)
> class LeftBoundary(SubDomain):
> def inside(self, x, on_boundary):
> return abs(x[0])<DOLFIN_EPS
> class RightBoundary(SubDomain):
> def inside(self, x, on_boundary):
> return abs(x[0]-1)<DOLFIN_EPS
>
> boundary_parts.set_all(2)

I assume that boundary_parts is defined
boundary_parts = MeshFunction('uint', mesh, 0)
or
boundary_parts = MeshFunction('sizet', mesh, 0)

> left = LeftBoundary()
> left.mark(boundary_parts, 0)
> right = RightBoundary()
> right.mark(boundary_parts, 1)
>
> u0 = Constant(0.0); u1 = Constant(3.0)
> bcs = [DirichletBC(V, u0, boundary_parts, 0, 'pointwise'),
> DirichletBC(V, u0, boundary_parts, 1, 'pointwise')] ..."""

I can reproduce it: this approach with mesh function does not work on
1.0.0 (segfault) as well as trunk with error:
*** Error: Unable to computing Dirichlet boundary values, pointwise
search. *** Reason: A SubDomain is required for pointwise search.

> When I tried to recover the boundary values like below:
> """...
> bcs[0].get_boundary_values()
> """

According to documentation this is not correct call. It requires at
least one more argument. Check
http://fenicsproject.org/documentation/dolfin/1.0.0/python/programmers-reference/cpp/DirichletBC.html#dolfin.cpp.DirichletBC.get_boundary_values

> The program was terminated and exited Python. Can you help me?

I suggest writing a bug report for 1D case with boundary markers
approach. Meanwhile you can use this approach:
...
class LeftBoundary(SubDomain):
 def inside(self, x, on_boundary):
               return x<DOLFIN_EPS
class RightBoundary(SubDomain):
 def inside(self, x, on_boundary):
               return x>1.-DOLFIN_EPS
bc_left = DirichletBC(V, Constant(0), LeftBoundary(), 'pointwise')
bc_right = DirichletBC(V, Constant(3), RightBoundary(), 'pointwise')
bc = [bc_left, bc_right]
...

Revision history for this message
Xiaoxian Liu (liuxiaox) said : 		#5

Thanks Jan Blechta, that solved my question.

Revision history for this message
Charles (rodbourn) said : 		#6

Why does on_boundary not work? This solved my problem as well - but threw me for a spin :)