How can I set the mean value for the solution of a PDE with neumann boundary value

Question

I want to solve the PDE

$$ -\nabla (a \nabla y) = 0 \text{ on } \Omega $$
$$ a\frac{\partial y}{\partial \nu} = g $$

Now the problem is, that the solution for this PDE is if and only if unique, if there is a requirement for the mean value. I want the mean value for the solution to be zero.

Now if I use

y = y - assemble(y*dx)

the mean value of y is close to zero.
Is this the right way to do this? It feels strange..

Code

phi_init = Constant(0)
phi = interpolate(phi_init, V)
g_h = Expression("(near(x[0], 1) && ! near(x[1], 0))|| near(x[1], 1) ? k0 : k1", k0=-0.5, k1= 0.5, degree= 1)

a_phi = a_1 * (1 + phi)/2 + a_2 * (1 - phi)/2

Define variational problem

y = TrialFunction(V)
v = TestFunction(V)
a = a_phi(inner(grad(y), grad(v)))dx
L = g_hvds

Compute solution

y = Function(V)
solve(a == L, y)
#y = y - assemble(y*dx)

Answer 1

What you have written is a Neumann boundary condition. A Robin boundary condition would be

$$ a \frac{\partial y}{\partial n} + b y = g, $$

and the assembled matrix would not be singular for $b>0$.

See the singular poisson demo for solving the singular Neumann problem. Solutions are unique only up to a scalar constant, so you can fix the solution to have any given mean value.

Comments

Sorry - I meant Neumann instead of Robin (I fixed that)

so this is the right way to set the mean value for the solution to zero?

Code

y = Function(V)
solve(a == L, y)
y = y - assemble(y*dx)

---

Yes. Once you found one solution, you obtain the solution with mean value zero that way.

Note though that solutions only exists when the right hand side vanishes on the constants.

---

I find a problem.

After y = y - assemble(y*dx), I can't compute error norm, like errornorm(y_e, y), where y_e is the exact solution.

When I run errornorm(y_e, y), I got error message "*** Reason: expected a GenericVector or GenericFunction."