FEniCS demo: demo_neumann-poisson.py

Question

I am looking at demo_neumann-poisson.py (#13 from the Documented DOLFIN demos), which solves a Neumann boundary-value problem for the Laplace operator over the unit square. I would like to ask the following two questions:

1) After reading the verbal description of the demo and carefully studying the variational formulation over $V \times \mathbb{R}$, I am left with the impression that the function $u$ which is plotted by demo_neumann-poisson.py is the unique member of $V$ which solves

\begin{eqnarray}
-\Delta u = f - c \text{ in } \Omega, \
\partial_n u = g \text{ on } \partial \Omega, \
\int_\Omega u \, dx \, dy = 0,
\end{eqnarray}

where $\Omega := [0,1]\times [0,1]$ and

\begin{equation}
c := \underbrace{\frac{1}{|\Omega|}}_{=1} \left[ \int_\Omega f\, dx + \int_{\partial \Omega} g \, ds \right].
\end{equation}

Could you confirm if my understanding is correct?

2) How can I impose a Neumann boundary condition $g = g(x,y)$ such as

\begin{equation}
g(0, y) := 1,\ g(1,y) := -1 \text{ for all } |y| \leq 1,
\end{equation}

\begin{equation}
g(x, 0) := 0,\ g(x, 1) := 0 \text{ for all } |x| \leq 1,
\end{equation}

in a FEniCS program?

Answer 1

Hello kdma-at-dtu,

1)
unfortunately, your understanding is not quite correct. So let's try to clear things up.
The function that is plotted is the unique solution to

$$ -\Delta u = f \;\text{in}\; \Omega, \;\text{and} \;\partial_n u = g \text{ on } \partial\Omega, \int_\Omega u dx dy = 0$$

The *weak form" of this PDE has the integrals from your definition of $c$ on the right hand side.

2)
to code your example, see the section "multiple neumann, robin and Dirichlet condition" in http://fenicsproject.org/documentation/tutorial/materials.html

Comments

@multigrid202: Thanks for your answer.

About 2): I will check out the tutorial, as you recommended.

About 1): I am afraid I don't quite understand what you mean; let me explain the source of my confusion: the weak form of the boundary-value problem you wrote is

\begin{equation}
\forall v \in H^1(\Omega): \quad \int_\Omega \nabla u \cdot \nabla v \, dx = \int_\Omega fv \, dx + \int_{\partial \Omega} g (v|_{\partial \Omega}) \, ds;
\end{equation}

this is just the interpretation of your boundary-value problem in terms of generalized functions on $\Omega$, and there is no reference to an entity $c$.

If we consider to have found a solution $(u,c) \in H^1(\Omega) \times \mathbb{R}$ to the weak formulation from http://fenicsproject.org/documentation/dolfin/dev/python/demo/documented/neumann-poisson/python/documentation.html, then:

\begin{eqnarray}
d = 1, v = 0 &\Rightarrow& \int_\Omega u \, dx = 0
\end{eqnarray}
\begin{eqnarray}
d = 0, v = 1 &\Rightarrow& \int_\Omega c \, dx = \int_\Omega f \, dx + \int_{\partial \Omega} g \, ds \ \Leftrightarrow \ c = \frac{1}{|\Omega|} \left( \int_\Omega f \, dx + \int_{\partial \Omega} g \, ds \right);
\end{eqnarray}
\begin{eqnarray}
d = 0, \forall v \in C_c^\infty &\Rightarrow& -\Delta u = f - c;
\end{eqnarray}

finally, applying $-\Delta u = f - c$ to the weak formulation from http://fenicsproject.org/documentation/dolfin/dev/python/demo/documented/neumann-poisson/python/documentation.html leads to
\begin{equation}
\forall v \in H^1(\Omega): \quad \int_{\partial \Omega} ( \partial_n u - g ) (v|_{\partial \Omega})\, ds = 0 \ \Rightarrow \ \partial_n u = g \text{ on } \partial \Omega.
\end{equation}

These computations seem to imply that the function $u$ which is plotted by demo_neumann-poisson.py is the unique member of V which solves the boundary-value problem I listed in my first post.

Could you kindly let me know your thoughts on this?

---

Hello again :)

Your explanation makes total sense so I actually have to thank you for clearing things up for me.

The reason I got confused was the following: In order for the original Problem (the one I wrote), $c$ must be zero for it to have any solution. I got confused because in that case, your formulation would have coincided with mine, but would have been an over-complication. In the case that $c\neq 0$, the original formulation wouldn't have made any sense (from a mathematical point of view) so I didn't really get why one would try to repair a (mathematically) meaningless problem.

Do you (or anybody else) happen to know any physical motivation for the "trick" used in the demo? Basically $c$ seems to be the unique constant shift that makes the original problem solvable. But this could have also been accomplished using other (non-constant) shifts on the right hand side.

As an example:

• The above BVP comes up when modeling axial dispalcement of a 1D bar with neumann conditions. If $c\neq 0$, then the Problem is also meaningless from a physical point of view and introducing $c$ simply solves the wrong problem

So to conclude, I am sorry for actually providing more confusion than answers, but would be happy to hear your thoughts on this matter.

---

@multigrid202: I'm glad we were able to clear things up!

As for your follow-up question, my purely mathematical sense is telling me that it makes sense to solve the Neumann boundary-value problem for the Laplace operator only if $f, g$ are compatible with each other.

I agree with you that the introduction of $c$ is more of trick than anything else; in other words, it seems to be the computational experimentalist's responsibility to only consider compatible $f, g$. Would anyone else care to weigh in on this?