# Integral constraint to replace Neumann Boundary Condition

**URL:** <https://fenicsproject.discourse.group/t/integral-constraint-to-replace-neumann-boundary-condition/2576>\
**Category:** variational formulation\
**Created:** [March 6, 2020, 4:30pm UTC](https://fenicsproject.discourse.group/t/integral-constraint-to-replace-neumann-boundary-condition/2576 "2020-03-06T16:30:17Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![scbarton](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/scbarton/32/858_2.png) [@scbarton](https://fenicsproject.discourse.group/u/scbarton)\
**Post date:** [March 6, 2020, 4:30pm UTC](https://fenicsproject.discourse.group/t/integral-constraint-to-replace-neumann-boundary-condition/2576/1 "2020-03-06T16:30:17Z")

</div>

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](https://fenicsproject.org/qa/5518/solve-a-pde-with-an-integral-term/) and [here](https://fenicsproject.discourse.group/t/how-to-impose-integral-constraint-on-pde/590) 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!

```auto
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()

```

---

<div class="post-metadata">

**Author:** ![kamensky](https://avatars.discourse-cdn.com/v4/letter/k/e95f7d/32.png) [@kamensky](https://fenicsproject.discourse.group/u/kamensky)\
**Post date:** [March 6, 2020, 4:45pm UTC](https://fenicsproject.discourse.group/t/integral-constraint-to-replace-neumann-boundary-condition/2576/2 "2020-03-06T16:45:44Z")

</div>

You can get rid of the implicit Neumann BC at x=0 by including the boundary term from integration by parts in your weak form:

```python
# Laplace Term
n = FacetNormal(mesh)
Fl = dot(grad(c), grad(v))*dx - dot(grad(c),n)*v*ds

```

That runs and gives the expected solution of c(x) = 0.5 + x.

---

<div class="post-metadata">

**Author:** ![scbarton](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/scbarton/32/858_2.png) [@scbarton](https://fenicsproject.discourse.group/u/scbarton)\
**Post date:** [March 6, 2020, 5:18pm UTC](https://fenicsproject.discourse.group/t/integral-constraint-to-replace-neumann-boundary-condition/2576/3 "2020-03-06T17:18:43Z")

</div>

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:

```auto
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()

```

---

<div class="post-metadata">

**Author:** ![kamensky](https://avatars.discourse-cdn.com/v4/letter/k/e95f7d/32.png) [@kamensky](https://fenicsproject.discourse.group/u/kamensky)\
**Post date:** [March 6, 2020, 5:37pm UTC](https://fenicsproject.discourse.group/t/integral-constraint-to-replace-neumann-boundary-condition/2576/4 "2020-03-06T17:37:28Z")

</div>

Okay, looking a little closer, one thing to notice is that the Lagrange multiplier `m` for the constraint is absent from the variational form `F`, so the constraint has no way to influence the PDE problem (weighted against `v`). Using

```python
# Conservation term
Fc = (c - Constant(C0)) * r * dx + m*v*ds

```

produces the expected solution in this case, although that looks weird to me, because it doesn’t derive from a Lagrangian. (I would expect something like `derivative(m*(c-Constant(C0))*dx,u)`, which would have `+ m*v*dx`, although that enforces the constraint through a volume term, which leads to a different solution.)

EDIT: Also, if you want `m` to have the interpretation of a flux, you can drop the boundary term from `Fl`, although the solution for `c` stays the same.

---

<div class="post-metadata">

**Author:** ![scbarton](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/scbarton/32/858_2.png) [@scbarton](https://fenicsproject.discourse.group/u/scbarton)\
**Post date:** [March 9, 2020, 12:48pm UTC](https://fenicsproject.discourse.group/t/integral-constraint-to-replace-neumann-boundary-condition/2576/5 "2020-03-09T12:48:46Z")

</div>

Yes, thanks that works. I thought I could get away without calculating `m`, since I could with the Dirichlet BCs.

---

<div class="post-metadata">

**Author:** ![giannoko](https://avatars.discourse-cdn.com/v4/letter/g/7c8e57/32.png) [@giannoko](https://fenicsproject.discourse.group/u/giannoko)\
**Post date:** [April 19, 2021, 1:17pm UTC](https://fenicsproject.discourse.group/t/integral-constraint-to-replace-neumann-boundary-condition/2576/6 "2021-04-19T13:17:26Z")

</div>

Hello,

I have the same problem. I want to solve the Navier-Stokes equation in a 2D rectangular domain. I would like to have a fully-developed flow at the inlet so I want to impose a constraint regarding the total volumetric flow rate. I followed the steps proposed here but I get errors. The unknowns of the problem are the velocity U(Ux,Uy), pressure P and I also defined a “Real” element for flow rate Q. At the wall I impose no-slip condition and at the inlet/outlet I include the term from the integration by parts. The code is converged for Dirichlet BCs at the inlet. Here is a sample of my code:

```
*# Define function spaces *
V = VectorElement('CG' ,mesh.ufl_cell(),2)
P = FiniteElement('CG' ,mesh.ufl_cell(),1)
Q = FiniteElement('Real',mesh.ufl_cell(),0)
VPQ = FunctionSpace(mesh,MixedElement([V,P,Q]))
phi,eta,ksi = TestFunctions(VPQ)

       	 
*# Set datum pressure*
def PinPoint(x, on_boundary ):
	tol = 1.e-4
	pointPres = near( x[0], 0.995833 , tol) and near(x[1], 0.005 , tol) 

	if (pointPres is True):
		return (pointPres)

*#BCs*
 Pref = Constant(0.0)
ns2 = Constant((0.0,0.0))

vFluidWall = DirichletBC(VPQ.sub(0), ns2 ,boundaries, wall ) 

bcp = DirichletBC(VPQ.sub(1),Pref,PinPoint,'pointwise')

bcs = [bcp , vFluidWall]
bcs_ = [bcp , vFluidWall]

*# Create functions*
W = Function(VPQ)
[bc.apply(W.vector()) for bc in bcs]
u,p,q = split(W) 

*#	stresses*
    def tau(X):
	nu = Constant(0.1)
	return nu*(grad(X)+grad(X).T)

def cauchy(X,P):
	return -P*Identity(len(X)) + tau(X)

#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
# *weak forms*
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-

dx = Measure("dx")(subdomain_data=subdomains)
ds = Measure("ds")(subdomain_data=boundaries) # outer boundary
dS = Measure("dS")(subdomain_data=boundaries) # inner boundary

n = FacetNormal(mesh)
Diff = inner( cauchy(u,p) , grad(phi))*dx() 
outFBC = - inner( dot( cauchy(u,p),n ),phi )*ds(outlet) 
inFBC = - inner( dot( cauchy(u,p),n ),phi )*ds(inlet)

Cont = eta*div(u)*dx()
Conv = inner(grad(u)*u, phi)*dx()

*# Conservation term*
Co = 0.01
Fc = (u - Constant(C0)) * ksi * dx + q*phi*ds(inlet) 

Res = Diff + Cont + Conv + outFBC + inFBC + Fc
Jac = derivative(Res,W)
Corr = Function(VPQ)

*# Solve*
NewtonRaphson(Res,Jac,bcs_,W,Corr,rank)

```

Thank you in advance