I’d like to replace a Neumann boundary condition with an integral constraint. A simple example is the Laplace equation \Delta c=0 on the unit interval.

I’d like to solve for c with a Dirichlet condition on the right (say c(1)=1.5) and setting the average value of c to 1.0. This should give the solution c=0.5 + x.

I’ve found a couple of examples here and here that I’ve tried to follow. My idea was to add an equation:

\int_x{(c-C_0)\ r\ dx}=0

using an additional test function, r. Based on the examples above, I used a `Real`

finite element for the additional test function. Sample code is given below that does not converge. I believe it is because of the implicit Neumann boundary condition at x=0. Any suggestions would be much appreciated!

```
from dolfin import *
import matplotlib.pyplot as pl
# Parameters
C0= 1.0 # conserved mean concentration
# Mesh
mesh = IntervalMesh(10,0,1)
# Function space
P2 = FiniteElement('Lagrange', mesh.ufl_cell(), 2)
R = FiniteElement('Real', mesh.ufl_cell(), 0)
element = MixedElement([P2, R])
V = FunctionSpace(mesh, element)
# Functions
u = Function(V)
u.assign(Constant([C0,C0]))
c,m = split(u)
v, r = TestFunctions(V)
f = Constant(0)
# Conservation term
Fc = (c - Constant(C0)) * r * dx
# Laplace Term
Fl = dot(grad(c), grad(v))*dx
# Combined terms
F= Fc+Fl
# Boundary Conditions
def left_boundary(x, on_boundary):
return on_boundary and near(x[0],0,1e-6)
def right_boundary(x, on_boundary):
return on_boundary and near(x[0],1.0,1e-6)
bc0=DirichletBC(V.sub(0), Constant(0.5), left_boundary)
bc1=DirichletBC(V.sub(0), Constant(1.5), right_boundary)
bc=[bc0, bc1]
# Solution with boring DirichletBC
solve(Fl == 0, u, bc) # this works
plot(c)
pl.show()
# Solution with integral constraint
solve(F == 0, u, bc1) # this does not converge
plot(c)
pl.show()
```