Marking Neumann boundary conditions isn't working

Question

Hey,

how do I mark Neumann boundary conditions correctly? After declaring the mesh as

nx = 120
xmax = 12
mesh1D = Interval(nx, 0, xmax)

I want u to fulfil the Neumann boundary condition (=g) only to the left border, i.e. x=0. I thought

class left(SubDomain):
  def inside(self, x, on_boundary):
    tol = 1E-13   # tolerance for coordinate comparisons
    return on_boundary and x[0] < tol 

Gamma_0 = left()
# Create mesh function over cell facets
exterior_facet_domains = FacetFunction("uint", mesh1D)
exterior_facet_domains.set_all(1)
Gamma_0.mark(exterior_facet_domains, 0)

would do the job, but strangely I have an influx at both ends, the boundary conditions is also applied to the other end of the interval, and not only to the left end.

If I change the mark of Gamma_0 to any value except 0, the boundary conditions doesn't seem to be applied at all. What did I do wrong?

Here the corresponding source code:

from dolfin import *
import numpy as np

# to shut off the output
set_log_level(WARNING)

nx = 120
xmax = 12
mesh1D = Interval(nx, 0, xmax)

# definition of solving method
V = FunctionSpace(mesh1D, 'CG', 1)

T = 20      
dt = 0.25
t = dt

# Define boundary segments for Neumann conditions
class left(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-13   # tolerance for coordinate comparisons
        return on_boundary and x[0] < tol

Gamma_0 = left()

# Create mesh function over cell facets
exterior_facet_domains = FacetFunction("uint", mesh1D)
exterior_facet_domains.set_all(1)
Gamma_0.mark(exterior_facet_domains, 0)

g = Expression('2.0')

# RHS
def rhs(u):
  return - u

# initial values
u_start = Constant(0.0)

u = Function(V)
u_old = Function(V)
du = TrialFunction(V)
v = TestFunction(V)
u_old = interpolate(u_start, V)

F = (u - u_old)*v*dx + dt*dot(grad(u), grad(v))*dx - dt*rhs(u)*v*dx - g*v*ds(0) 
dF = derivative(F,u, du)
problem = NonlinearVariationalProblem(F, u, J = dF)
pdesys_newton = NonlinearVariationalSolver(problem)

#constants for ploting
x_points = []
for i in range(0,nx):
  x_points.append((xmax*i)/(1.0*nx)) 

x_array = np.array(x_points)

ymax = 0

solution = []

#--------------- Solving the Equation
while t <= T:
    #print 'time =', t

    u_array = u_old.vector().array()
    u_nodal_values = []

    for i in range(0,nx):
      u_nodal_values.append(u_array[i])

    pdesys_newton.solve()

    t += dt
    t_out += dt
    u_old.assign(u)

Answer 1

Add

ds = ds[exterior_facet_domains]

or

problem = NonlinearVariationalProblem(F, u, J = dF,
              exterior_facet_domains=exterior_facet_domains)

Comments

Thanks a lot Jan, your first suggestions helped me, although I got another question:
the boundary condition is applied to the right end of my interval, instead of the left end.
I added

class right(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-13   # tolerance for coordinate comparisons
        return on_boundary and x[0] >= xmax - DOLFIN_EPS

Gamma_r = right()
....
h = Expression('2')

as another Neumann boundary to be applied to the right end.
The result is the opposite: h is applied at x=0, and g is applied at x=xmax. Why?

Could you please as well explain what's happening, why did I have to add the line

ds = ds[exterior_facet_domains]

?

---

Please, submit a complete and short code demonstrating the problem.

ds

is a predefined measure on a whole exterior facets and

ds[exterior_facet_domains]

assigns to it markers based on users's data so that you're able to use a measure

ds(i)

on subset of ds where exterior_facet_domains take value i.

---

Sorry, my output seemed to cause the problem and switched the sides, it works, thank you a lot :-)