%----------------------------------------------------
%\section{Non-homogeneous Neumann boundary conditions for the Laplace operator}
%----------------------------------------------------
\label{sec-neumann-laplace}
\pbindex{Poisson}%
\cindex{boundary condition!Neumann}%

In this chapter we study how to solve a ill-posed problem
with a solution defined up to a constant.

\subsubsection*{Formulation}
  Let $\Omega$ be a bounded open and simply connected subset of $\mathbb{R}^d$, 
  $d=1,2$ or $3$.
  Let $f \in L^{2}(\Omega)$ and
  $g \in H^{\frac{1}{2}}(\partial \Omega)$ satisfying
  the following compatibility condition:
  \[
	\int_\Omega f \, {\rm d}x
	+
	\int_{\partial\Omega} g \, {\rm d}s
	=
	0
  \]
  The problem writes:\\
  $(P_5)_h$: {\it find $u$, defined in $\Omega$ such that:}
  \begin{eqnarray*}
      -\Delta u &=& f \ {\rm in}\ \Omega \\
      \Frac{\partial u}{ \partial n} &=& g \ {\rm on}\ \partial \Omega
  \end{eqnarray*}
\cindex{method!conjugate gradient algorithm}%
\cindex{method!minres algorithm}%
  Since this problem only involves the derivatives of $u$,
  it is defined up to a constant:
  it is then clear that its solution is never unique~\citep[p.~11]{GirRav-1986}.
  A discrete version of this problem could be solved iteratively
  by the conjugate gradient or the MINRES algorithm~\citep{PaiSau-1975}.
  In order to solve it by a direct method,
  we get round this difficulty by seeking~$u$ in the following space
  \[	
	V = \{ v\in H^1(\Omega);\ \  b(v,1) = 0 \}
  \]
  where
  \[
	b(v,\mu) = \mu \, \int_\Omega v \, {\rm d}x,
		\ \ \forall v\in L^2(\Omega),
		\forall \mu\in \mathbb{R}
  \]
  The variational formulation of this problem writes:

  $(VF_5)$: {\it find $u \in V$ such that:}
  $$
    a(u,v) = \ell(v), \ \forall v \in V
  $$
  where
  \begin{eqnarray*}
    a(u,v) &=& \int_\Omega \nabla u . \nabla v \, \mathrm{d}x \\
    \ell(v) &=& \int_\Omega f v \, \mathrm{d}x \\
           + \int_{\partial \Omega} g v \, \mathrm{d}s
  \end{eqnarray*}
  Since the direct discretization of the space $V$ is not
  an obvious task,
  the constraint $b(u,1) = 0$ is enforced by a Lagrange multiplier
\cindex{Lagrange!multiplier}%
  $\lambda\in \mathbb{R}$.
  Let us introduce the Lagrangian, defined
  for all $v\in H^1(\Omega)$ and $\mu\in\mathbb{R}$ by:
  \[
    L(v,\mu) = \Frac{1}{2} a(v,v) + b(v,\mu) - \ell(v)
  \]
  The saddle point $(u,\lambda)\in H^1(\Omega)\times\mathbb{R}$
  of this Lagrangian is characterized as the unique solution
  of:
  \begin{eqnarray*}
	a(u,v) + b(v,\lambda)
		&=& \ell(v), \ \ \forall v\in H^ 1(\Omega) \\
	b(u,\mu) \phantom{+ b(v,\lambda)}
	        &=& 0,    \ \ \ \ \  \forall \mu\in\mathbb{R}
  \end{eqnarray*}
  It is clear that if $(u,\lambda)$ is solution of this problem,
  then $u\in V$ and $u$ is a solution of $(VF_5)$.
  Conversely, let $u\in V$ the solution of $(VF_5)$.
  Choosing $v=v_0$ where $v_0(x)=1$, $\forall x\in\Omega$
  leads to 
  \mbox{$
     \lambda\,{\rm meas}(\Omega)=\ell(v_0)
  $}.
  From the definition of $\ell(.)$ and the compatibility 
  condition between the data $f$ and $g$, we get $\lambda=0$.
  Note that the saddle point problem extends to the case
  when $f$ and $g$ does not satisfies the compatibility condition,
  and in that case $\lambda=\ell(v_0)/{\rm meas}(\Omega)$.
  
\subsubsection*{Approximation}
\cindex{Lagrange!interpolation}
  As usual, we introduce a mesh ${\cal T}_h$ of $\Omega$
  and the finite dimensional space $X_h$:
  $$
      X_h = \{ v \in H^1(\Omega); \
          v_{/K} \in P_k, \
          \forall K \in {\cal T}_h \}
  $$
  The approximate problem writes:\\
  {\it $(VF_5)_h$: find $(u_h,\lambda_h)\in X_h\times \mathbb{R}$ such that:}
  \begin{eqnarray*}
	a(u_h,v) + b(v,\lambda_h)
		&=& \ell(v), \ \ \forall v\in X_h \\
	b(u_h,\mu) \phantom{+ b(v,\lambda)}
	        &=& 0,    \ \ \ \ \  \forall \mu\in\mathbb{R}
  \end{eqnarray*}

% ------------------------------
\myexamplelicense{neumann-laplace.cc}
% ------------------------------

\subsubsection*{Comments}
  Let $\Omega \subset \mathbb{R}^d$, $d=1,2,3$.
  We choose $f(x) = 1$ and $g(x) = -1/(2d)$.
  This example is convenient, since the exact solution is known:
  $$
	u(x) = - \Frac{1}{12} + \Frac{1}{2d} \sum_{i=1}^d x_i(1-x_i)
  $$
  The code looks like the previous ones.
  Let us comment the changes.
  The discrete bilinear form $b$ is computed as $b\in X_h$
  that interprets as a linear application from $X_h$ to $\mathbb{R}$:
  $b(v_h)=m(v_h,1)$. Thus $b$ is computed as
\begin{lstlisting}[numbers=none,frame=none]
  field b = integrate(v);
\end{lstlisting} 
  Let 
  \[
     \mathcal{A}
	=
     \left( \begin{array}{cc} 
        {\tt a} & {\tt trans(b)} \\ 
        {\tt b}  & 0 
     \end{array} \right)
	, \ \ \ 
     \mathcal{U}
     	=
     \left( \begin{array}{c} 
        {\tt uh} \\ 
        {\tt lambda} 
     \end{array} \right)
	, \ \ \ 
     \mathcal{B}
     	=
     \left( \begin{array}{c} 
        {\tt lh} \\
        0
     \end{array} \right)
  \]
  The problem admits the following matrix form:
  \[
     \mathcal{A} \ \mathcal{U} = \mathcal{B}
  \]
\cindex{matrix!concatenation}
  The matrix that represents the bilinear fo
  and its right-hand side are assembled as:
\begin{lstlisting}[numbers=none,frame=none]
  form  A = {{   a,      b},
             { trans(b), 0}};
  field Bh = {   lh,     0};
\end{lstlisting}
Both the pairs $\mathcal{U}=(u_h,\lambda)$
and $\mathcal{B}=(b,0)$ belong to the
vectorial space~$X_h\times \mathbb{R}$.
and $\mathcal{B}=(b,0)$ belong to the
Then, the variable \code{Uh} could be declared as:
\begin{lstlisting}[numbers=none,frame=none]
  field Uh (Bh.get_space(), 0); 
\end{lstlisting}
%
\cindex{matrix!singular}%
\cindex{matrix!indefinite}%
\findex{ldlt}%
Note that the matrix $\mathcal{A}$ is symmetric and non-singular, but indefinite~: 
it admits eigenvalues that are either strictly positive or strictly negative.
While the Choleski factorization is not possible, its variant the $LDL^T$ one
is performed, thanks to the \code{ldlt} function:
\begin{lstlisting}[numbers=none,frame=none]
  A.set_symmetry(true);
  problem p (A);
  p.solve (Bh, Uh);
\end{lstlisting}
Then, the \code{uh} field is extracted as the first
component of the the \code{Uh} one:
\begin{lstlisting}[numbers=none,frame=none]
  dout << Uh[0];
\end{lstlisting}
%% \findex{catchmark}%
%%   Finally, the statement
%% \begin{lstlisting}[numbers=none,frame=none]
%%   dout << catchmark("u") << uh
%%        << catchmark("lambda") << lambda << endl;
%% \end{lstlisting}
%%   writes the solution $(u_h,\lambda)$.
%%   The \code{catchmark} function writes marks together with the solution in the output stream.
%%   These marks are suitable for a future reading with the same format, as:
%% \begin{lstlisting}[numbers=none,frame=none]
%%   din  >> catchmark("u") >> uh
%%        >> catchmark("lambda") >> lambda;
%% \end{lstlisting}
%%   This is useful for post-treatment, visualization and error analysis.

\subsubsection*{How to run the program}

  As usual, enter:
\begin{verbatim}
  make neumann-laplace
  mkgeo_grid -t 10 > square.geo
  ./neumann-laplace square P1 | field -
\end{verbatim}