Integral constraint to replace Neumann Boundary Condition

variational formulation
3

Thanks, for the quick reply. Your solution makes sense but it only works if the problem with Dirichlet BCs is solved first. If I comment out the “boring DirichletBC” section and run the revised code it still does not converge. Any other suggestions?

Here is the revised code based on your suggestion:

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
Fl0 = dot(grad(c), grad(v))*dx

n = FacetNormal(mesh)
Fl = dot(grad(c), grad(v))*dx - dot(grad(c),n)*v*ds

# 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(Fl0 == 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()