Constrain volume integral

I would like to solve the following:

∇(u) = 0
∫u = Q*vol

Solution:
u = Q

Here is what I tried:

from fenics import *
mesh = UnitIntervalMesh(10)
vol = assemble(Constant(1)*dx(domain=mesh))
Q = 1
V = FunctionSpace(mesh, "P", 1)
v = TestFunction(V)
u = Function(V)
F = inner(grad(u)[0],v) * dx + (u - Q/vol) * dx
solve(F == 0, u)

However I get an error:

ArityMismatch: Adding expressions with non-matching form arguments () vs ('v_0',).

How to formulate this problem correctly?

The volume of a mesh can be computed by

vol = assemble(Constant(1) * dx)
1 Like

Thanks! I edited this into my original post.

I would suggest having a look at this Poisson with only Neumann conditions demo
Which shows you how to use lagrange multipliers in the case Q=0. It is very easy to add a non-zero Q to this constraint

1 Like
from fenics import *
mesh = UnitIntervalMesh(10)
vol = assemble(Constant(1)*dx(domain=mesh))
cell = mesh.ufl_cell()
fe_rho      = FiniteElement("P", cell, 1)
fe_lagrange = FiniteElement("R", cell, 0)
fe = MixedElement([fe_rho, fe_lagrange])
V = FunctionSpace(mesh, fe)
Qbar = project(Constant(1/vol), V.sub(0).collapse())
v_rho, v_l = TestFunctions(V)
u = Function(V)
rho, _ = split(u)
F = inner(grad(rho)[0], v_rho) * dx + inner(rho - Qbar, v_l) * dx
solve(F == 0, u)

Seems to work. What looks a bit strange to me is that along with the “test scalar” l comes a new unknown rho, _ = split(u) that I just ignore. Does this confuse solvers? Seems bad for condition numbers etc?

Add it to the system as done in the Poisson demo. Otherwise you end up with a Zero entry.
As you have a linear system, i suggest trying to solve it with an LU solver first, as it should then become clear that you Need an additional entry.

Thanks! The following works for me:

from fenics import *
import numpy as np 

mesh = UnitIntervalMesh(4)

vol = assemble(Constant(1)*dx(domain=mesh))
cell = mesh.ufl_cell()
fe_rho      = FiniteElement("P", cell, 1)
fe_lagrange = FiniteElement("R", cell, 0)
fe = MixedElement([fe_rho, fe_lagrange])
V = FunctionSpace(mesh, fe)
Qbar = project(Constant(1/vol), V.sub(0).collapse())
v_rho, v_l = TestFunctions(V)
u = TrialFunction(V)
rho, l = split(u)

class LeftBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], 0)

facedim = mesh.topology().dim() - 1
boundary_markers = MeshFunction("size_t", mesh, facedim, 0)
LeftBoundary().mark(boundary_markers, 1)

ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)

F = inner(grad(rho)[0], v_rho) * dx + inner(rho - Qbar, v_l) * dx + inner(rho -l, v_rho) * ds(1)
sol = Function(V)
A = assemble(lhs(F)).array()
print(f"A = {A}")
print(f"cond(A) = {np.linalg.cond(A)}") # better way to compute cond without going to dense matrix?
solve(lhs(F) == rhs(F), sol)
rho, l = split(sol)
plot(rho)

It feels slightly weird to make up the boundary condition rho = l though.