\label{sec-surface}
\cindex{geometry!surface}%
\cindex{method!level set}%
This chapter deals with equations defined on a closed hypersurface. 
We present three different numerical methods: the direct resolution
of the problem on an explicit surface mesh generated independently
of \Rheolef, the direct resolution on a surface mesh generated
by \Rheolef\  from a volume mesh, and finally a level set type method
based on a volume mesh in an $h$-narrow band containing the surface.
This last method allows one to define hybrid operators between surface and
volume-based finite element fields.
These methods are demonstrated on two model problems and two
different surfaces.

% ---------------------------------------------
%\subsubsection*{Problems statement}
% ---------------------------------------------

Let us consider a closed surface $\Gamma \in \mathbb{R}^d$, $d=2$ or $3$ and $\Gamma$
is a connected $C^2$ surface of dimension $d-1$ with $\partial\Gamma=0$.
We first consider the following problem:\\
\mbox{} \ \ $(P1)$ {\it find $u$, defined on $\Gamma$ such that:}
  \begin{eqnarray}
      u -\Delta_s u &=& f \ {\rm on}\ \Gamma
      \label{eq-helmholtz-s}
  \end{eqnarray}
where $f \in L^2(\Gamma)$.
\cindex{operator!Laplace-Beltrami}%
\cindex{operator!Helmholtz-Beltrami}%
For all function $u$ defined on $\Gamma$,
$\Delta_s$ denotes the Laplace-Beltrami operator:
\[
	\Delta_s u = {\rm div}_s (\nabla_s u)
\]
where $\nabla_s$ and ${\rm div}_s$
are the tangential derivative and the surface divergence
along $\Gamma$, defined respectively,
for all scalar field $\varphi$ and vector field ${\bf v}$ by:
\begin{eqnarray*}
	\nabla_s \varphi &=& (I - \mathbf{n} \otimes \mathbf{n})\, \nabla \varphi \\
	{\rm div}_s \, {\bf v} &=& (I - \mathbf{n} \otimes \mathbf{n}) : \nabla {\bf v}
\end{eqnarray*}
Here, $\mathbf{n}$ denotes a unit normal on $\Gamma$.

We also consider the following variant of this problem:\\
\mbox{} \ \ $(P2)$ {\it find $u$, defined on $\Gamma$ such that:}
  \begin{eqnarray}
      -\Delta_s u &=& f \ {\rm on}\ \Gamma
      \label{eq-laplace-s}
  \end{eqnarray}
This second problem is similar to the first one: the Helmholtz
operator $I-\Delta_s$ has been replaced by the Laplace-Beltrami one $-\Delta_s$. 
In that case, the solution is defined up to a constant:
if $u$ is a solution, then $u+c$ is also a solution for any constant $c\in\mathbb{R}$.
Thus, we refers to $(P1)$ as the Helmholtz-Beltrami problem 
and to $(P2)$ as the Laplace-Beltrami one. 

% -------------------------------------------------
\subsection{Approximation on an explicit surface mesh}
% -------------------------------------------------
% -----------------------------------------
\subsubsection{The Helmholtz-Beltrami problem}
% -----------------------------------------
\cindex{Green formula}%
Tanks to the surface Green formula (see appendix~\ref{sec-green-surface}),
the variational formulation of problem $(P1)$ writes:\\
\mbox{}\ \ $(VF1)$: {\it find $u \in H^1(\Gamma)$ such that:}
  $$
    a(u,v) = l(v), \ \forall v \in H^1(\Gamma)
  $$
  where for all $u,v\in H^1(\Gamma)$,
  \begin{eqnarray*}
    a(u,v) &=&  \int_\Gamma \left( u \, v + \nabla_s u . \nabla_s v \right) \, {\rm d}s \\
    l(v) &=& \int_\Gamma f \, v \, {\rm d}s
  \end{eqnarray*}

Let $k\geq 1$ and consider
a $k$-th order curved surface finite element mesh $\Gamma_h$ of $\Gamma$.
We define the space $W_h$:
$$
	W_h = \left\lbrace v_h \in H^1(\Gamma_h) ; v_{|_S} \in P_k, \forall S \in \Gamma_h \right\rbrace
$$
The approximate problem writes:\\
\mbox{}\ \ {\it $(VF1)_h$: find $u_h\in W_h$ such that:}
  \begin{eqnarray*}
        a(u_h,v_h) &=& l(v_h), \ \ \forall v_h\in W_h \\
  \end{eqnarray*}

% -----------------------------------------
\myexamplelicense{helmholtz_s.cc}
% -----------------------------------------

% -----------------------------------------
\subsubsection*{Comments}
% -----------------------------------------
The problem involves the Helmholtz operator
and thus, the code is similar to
\exindex{neumann-nh.cc}%
\reffile{neumann-nh.cc}
presented page~\pageref{neumann-nh.cc}.
\findex{integrate}%
\cindex{form!{$\nabla_s u.\nabla_s v+uv$}}%
Let us comments the only differences:
\begin{lstlisting}[numbers=none,frame=none]
  form  a  = integrate (u*v + dot(grad_s(u),grad_s(v)));
\end{lstlisting}
The form refers to the \code{grad_s} operator
instead of the \code{grad} one, since only 
the coordinates related to the surface are involved.
\begin{lstlisting}[numbers=none,frame=none]
  field lh = integrate (f(d)*v);
\end{lstlisting}
The right-hand-side does not involve any boundary term,
since the surface $\Gamma$ is closed: the boundary domain $\partial \Gamma=\emptyset$.
%
As test problem, the surface $\Gamma$ is the unit circle when
$d=2$ and the unit sphere when $d=3$.
\cindex{benchmark!Dziuk-Elliott-Heine on a sphere}%
The data $f$ has been chosen as in \citet[p.~17]{DecDziEllHei-2009}.
This choice is convenient since the exact solution is known.
Recall that the spherical coordinates $(\rho,\theta,\phi)$
are defined from the Cartesian ones $(x_0,x_1,x_2)$ by:
\cindex{coordinate system!spherical}%
\cindex{geometry!sphere}%
\[
    \rho = \sqrt{x_0^2 + x_1^2 + x_2^2}	
	,\ \ 
    \phi = \arccos\left(x_2/\rho\right)
	,\ \ 
    \theta = \left\{
	      \begin{array}{ll}
                \arccos\left(x_0/\sqrt{x_0^2+x_1^2}\right) & \mbox{ when } x_1 \geq 0 \\
                2\pi - \arccos\left(x_0/\sqrt{x_0^2+x_1^2}\right) & \mbox{ otherwise }
              \end{array}
             \right.
\]
% -----------------------------------------
\myexamplelicense{sphere.icc}
% -----------------------------------------

% -------------------------------
\subsubsection*{How to run the program}
% -------------------------------
The program compile as usual:
\begin{verbatim}
  make helmholtz_s
\end{verbatim}
A mesh of a circle is generated by:
\pindex{mkgeo_ball}%
\pindexopt{mkgeo_ball}{-s}%
\pindexopt{mkgeo_ball}{-e}%
\cindex{geometry!circle}%
\begin{verbatim}
  mkgeo_ball -s -e 100 > circle.geo
  geo circle -gnuplot
\end{verbatim}
The \code{mkgeo_ball} is a convenient script that
\pindex{gmsh}%
generates a mesh with the {\tt gmsh} mesh generator.
Then, the problem resolution writes:
\begin{verbatim}
  ./helmholtz_s circle P1 > circle.field
  field circle.field 
  field circle.field -elevation
\end{verbatim}
The tridimensional case is similar:
\pindexopt{mkgeo_ball}{-t}%
\pindexopt{geo}{-stereo}%
\pindexopt{field}{-stereo}%
\begin{verbatim}
  mkgeo_ball -s -t 10 > sphere.geo
  geo sphere.geo -stereo
  ./helmholtz_s sphere.geo P1 > sphere.field
  field sphere.field
  field sphere.field -stereo 
\end{verbatim}
The solution is represented on Fig~.\ref{fig-sphere-s-geo}.left.

\begin{figure}[htb]
     %\begin{center}
       \begin{tabular}{ccc}
	  %TODO: le maillage correspond a n=5 et la solution a n=10 !
          % sur le maillage n=10 on ne voit pas les aretes courbes...
          \includegraphics[height=6.5cm]{sphere_s_P1_geo.png}   &
          \includegraphics[height=6.5cm]{sphere_s_P1_field.png} & \\
          \includegraphics[height=6.5cm]{sphere_s_P3_geo.pdf}   &
          \includegraphics[height=6.5cm]{sphere_s_P3_field.png} &
	  \stereoglasses
       \end{tabular}
     %\end{center}
     \caption{Helmholtz-Beltrami problem:
	high-order curved surface mesh
	and its corresponding isoparametric solution:
	(top) $order=1$; (bottom) $order=3$.}
     \label{fig-sphere-s-geo}
\end{figure}
Higher-order isoparametric finite elements can be considered for the
curved geometry:
\cindex{geometry!curved}%
\apindex{isoparametric}%
\pindexopt{geo}{-subdivide}%
\begin{verbatim}
  mkgeo_ball -s -e 30 -order 3 > circle-P3.geo
  geo circle-P3.geo -subdivide 10
\end{verbatim}
Observe the curved edges (see Fig~.\ref{fig-sphere-s-geo}).
The \code{-subdivide} option allows a graphical representation of the
curved edges by subdividing each edge in ten linear parts, since graphical
softwares are not yet able to represent curved elements.
The computation with the $P_3$ isoparametric approximation writes:
\begin{verbatim}
  ./helmholtz_s circle-P3 P3 > circle-P3.field
  field circle-P3.field -elevation -gnuplot
\end{verbatim}
Note that both the curved geometry and the finite element are second order.
The tridimensional counterpart writes simply:
\begin{verbatim}
  mkgeo_ball -s -t 10 -order 3 > sphere-P3.geo
  geo sphere-P3.geo -gnuplot
  ./helmholtz_s sphere-P3 P3 > sphere-P3.field
  field sphere-P3.field
  field sphere-P3.field -stereo 
\end{verbatim}
The solution is represented on Fig~.\ref{fig-sphere-s-geo}).right-bottom.
The graphical representation is not yet able to represent the high-order 
approximation: each elements is subdivided and a piecewise linear representation is
used in each sub-elements.

\cindex{error analysis}%
\cindex{convergence!error!versus mesh}%
\cindex{convergence!error!versus polynomial degree}%
\begin{figure}[htb]
  \begin{center}
    \begin{tabular}{cc}
      \includegraphics{cvge-helmholtz-s-sphere-l2.pdf} &
      \includegraphics{cvge-helmholtz-s-sphere-linf.pdf} \\
      \includegraphics{cvge-helmholtz-s-sphere-h1.pdf} &
    \end{tabular}
  \end{center}
  \caption{Curved non-polynomial surface:
	error analysis in $L^2$, $L^\infty$ and $H^1$ norms.}
  \label{fig-surface-err}
\end{figure}
\myexamplenoinput{helmholtz_s_error.cc}%
Since the exact solution is known, the error can be computed:
this is done by the program \code{helmholtz_s_error.cc}.
This file is not presented here, as it is similar to some others
examples, but can be founded in the \Rheolef\  example directory.
Figure~\ref{fig-surface-err} plots the error in various norms
versus element size for different isoparametric approximations.

\clearpage
% -----------------------------------------
\subsubsection{The Laplace-Beltrami problem}
% -----------------------------------------

This problem has been introduced in~\eqref{eq-laplace-s}, page~\pageref{eq-laplace-s}.
While the treatment of the Helmholtz-Beltrami problem was similar to the
Helmholtz problem with Neumann boundary conditions, 
here, the treatment of the Laplace-Beltrami problem is similar to the
Laplace problem with Neumann boundary conditions:
see section~\ref{sec-neumann-laplace}, page~\pageref{sec-neumann-laplace}.
Note that for both problems, the solution is defined up to a constant.
Thus, the linear problem has a singular matrix.
The \reffile{laplace_s.cc} code is similar
to the \reffile{neumann-laplace.cc} one, as presented in section~\ref{sec-neumann-laplace}.
The only change lies one the definition of the right-hand side.

% -----------------------------------------
\myexamplelicense{laplace_s.cc}

\myexamplelicense{torus.icc}
% -----------------------------------------

As test problem, the surface $\Gamma$ is the a torus when $d=3$.
\cindex{benchmark!Olshanskii-Reusken-Grande on a torus}%
The data $f$ has been chosen as in~\citet[p.~3355]{OlsReuGra-2009}.
This choice is convenient since the exact solution is known.
\cindex{coordinate system!torus}%
Let $R$ and $r$ denotes the large and small torus radii, respectively.
The torus coordinates $(\rho,\theta,\phi)$ are defined linked
to the Cartesian ones by:
\[
   \left( \begin{array}{c}
     x_0 \\ x_1 \\ x_2
   \end{array} \right)
   =
   R
   \left( \begin{array}{c}
     \cos(\phi) \\ \sin(\phi) \\ 0
   \end{array} \right)
   +
   \rho
   \left( \begin{array}{r}
     \cos(\phi) \cos(\theta) \\ \sin(\phi) \cos(\theta)  \\ \sin(\theta)
   \end{array} \right)
\]
Here $\rho$ is the distance from the point to the circle in the $x_0x_1$ plane
around $0$ with radius $R$, $\theta$ is the angle from the positive $(x_0,x_1,0)$
to $x_0$ and $\phi$ is the angle from the positive $x_0$ axis to $(x_0,x_1,0)$.

% -------------------------------------
\subsubsection*{How to run the program ?}
% -------------------------------------
\begin{figure}[htb]
     %\begin{center}
       \begin{tabular}{ccc}
          \includegraphics[width=6.5cm]{torus_s_P1_geo.png}   &
          \includegraphics[width=6.5cm]{torus_s_P1_field.png} & \\
          \includegraphics[width=6.5cm]{torus_s_P2_geo.pdf}   &
          \includegraphics[width=6.5cm]{torus_s_P2_field.png} &
	  \stereoglasses
       \end{tabular}
     %\end{center}
     \caption{Laplace-Beltrami problem on a torus:
	high-order curved surface mesh
	and its corresponding isoparametric solution:
	(top) $order=1$; (bottom) $order=2$.}
     \label{fig-torus-s-geo}
\end{figure}
\fiindex{\filesuffix{.mshcad} gmsh geometry}%
\pindex{gmsh}%
\cindex{geometry!torus}%
The surface mesh of the torus is generated by:
\begin{verbatim}
    gmsh -2 torus.mshcad -format msh2 -o torus.msh
    msh2geo torus.msh > torus.geo
    geo torus.geo -stereo
\end{verbatim}
\myexamplenoinput{torus.mshcad}%
The \reffile{torus.mshcad} is not presented here:
it can be founded in the \Rheolef\  example directory.
Then, the computation and visualization writes:
\begin{verbatim}
    make laplace_s
    ./laplace_s torus.geo P1 > torus.field
    field torus.field
    field torus.field -stereo 
\end{verbatim}
For a higher-order approximation:
% TODO P3: 
% gmsh -2 -order 3 torus.mshcad -format msh2 -o torus-P3.msh
% msh2geo torus-P3.msh > torus-P3.geo
% 	fatal{0}(geo_seq_upgrade.cc,429): invalid permutation
% =>  bug in gmsh: extrude and high order ?
\pindexopt{geo}{-gnuplot}%
\begin{verbatim}
    gmsh -2 -order 2 torus.mshcad -format msh2 -o torus-P2.msh
    msh2geo torus-P2.msh > torus-P2.geo
    geo torus-P2.geo -gnuplot
    ./laplace_s torus-P2.geo P2 > torus-P2.field
    field torus-P2.field -stereo
\end{verbatim}
The solution is represented on Fig.~\ref{fig-torus-s-geo}.
By editing \reffile{torus.mshcad} and changing the
density of discretization, we can improve the approximate solution
and converge to the exact solution.
Due to a bug~\citep{gmsh-bug-curved-high-order} in the current gmsh version~2.5.1
the convergence is not optimal ${\cal O}(h^k)$ for higher values of $k$.

\clearpage
% -------------------------------------------------------------------
\subsection{Building a surface mesh from a level set function}
% -------------------------------------------------------------------
\cindex{method!level set}%
The previous method is limited to not-too-complex surface $\Gamma$,
that can be described by a regular finite element surface mesh $\Gamma_h$.
When the surface change, as in a time-dependent process, complex
change of topology often occurs and the mesh $\Gamma_h$ can degenerate or 
be too complex to be efficiently meshed.
In that case, the surface is described implicitly as
the zero isosurface, or zero {\em level set}, of a function:
\[
    \Gamma = \{ x\in \Lambda; \ \ \phi(x) = 0 \}
\]
where $\Lambda\subset\mathbb{R}^d$ is a bounding box of the surface $\Gamma$.

The following code automatically generates
the mesh $\Gamma_h$ of the surface described by the zero isosurface
of a discrete $\phi_h\in X_h$ level set function:
\[
    \Gamma_h = \{ x\in \Lambda; \ \ \phi_h(x) = 0 \}
\]
where $X_h$ is a piecewise affine functional space over a mesh
${\cal T}_h$ of $\Lambda$:
\[
      X_h = \{ \varphi \in L^2(\Lambda) \cap C^0(\Lambda); \ 
          \varphi_{/K} \in P_1, \ 
          \forall K \in {\cal T}_h \}
\]
The polynomial approximation is actually limited here to first
order: building higher order curved finite element surface meshes
from a level set function is planed for the future versions of \Rheolef.

Finally, a computation, as performed
in the previous paragraph can be done using $\Gamma_h$.
We also point out the limitations of this approach.

% -----------------------------------------
\myexamplelicense{level_set_sphere.cc}
% -----------------------------------------

% -----------------------------------------
\subsubsection*{Comments}
% -----------------------------------------
\findex{level_set}%
\clindex{level_set_option}%
All the difficult work of building the intersection mesh $\Gamma_h$,
defined as the zero level set of the $\phi_h$ function,
is performed by the \code{level_set} function:
\begin{lstlisting}[numbers=none,frame=none]
  geo gamma = level_set (phi_h, opts);
\end{lstlisting}
When $d=3$, intersected tetrahedra leads to either triangular or quadrangular faces.
By default, quadrangular faces are split into two triangles.
An optional \code{-tq} program flag allows one to conserve quadrangles in the surface mesh:
it set the \code{split_to_triangle} optional field to false.

% -----------------------------------------
\subsubsection*{How to run the program ?}
% -----------------------------------------
\begin{figure}[htb]
  \begin{tabular}{ccc}
    \includegraphics[height=5.5cm]{level-set-circle-geo.pdf} &
    \includegraphics[height=5.5cm]{level-set-circle-field.png} & \\
    \includegraphics[height=5.5cm]{level-set-sphere-geo.png} &
    \includegraphics[height=5.5cm]{level-set-sphere-field.png} &\\
    \includegraphics[height=5.0cm]{level-set-torus-geo.png} &
    \includegraphics[height=5.0cm]{level-set-torus-field.png} &
    \stereoglasses
  \end{tabular}
  \caption{Building an explicit surface mesh from level set:
	(top) circle;
        (center) sphere;
        (bottom) torus.
  }
  \label{fig-intersection-sphere}
\end{figure}

After the compilation, generates the mesh
of a bounding box $\Lambda=[-2,2]^d$ of the surface and run the program:
\begin{verbatim}
  make level_set_sphere
  mkgeo_grid -t 20 -a -2 -b 2 -c -2 -d 2 > square2.geo
  ./level_set_sphere square2.geo > circle.geo
  geo circle.geo -gnuplot
\end{verbatim}
The computation of the previous paragraph can be reused:
\begin{verbatim}
  ./helmholtz_s circle.geo P1 | field -paraview -
\end{verbatim}
Note that, while the bounding box mesh was uniform,
the intersected mesh could present arbitrarily small edge length
(see also Fig.~\ref{fig-intersection-sphere}):
\begin{verbatim}
  geo -min-element-measure circle.geo
  geo -max-element-measure circle.geo
\end{verbatim}

Let us turn to the $d=3$ case: 
\begin{verbatim}
  mkgeo_grid -T 20 -a -2 -b 2 -c -2 -d 2 -f -2 -g 2 > cube2.geo
  ./level_set_sphere cube2.geo | geo -upgrade - > sphere.geo
  geo sphere.geo -stereo
  ./helmholtz_s sphere.geo P1 | field -
\end{verbatim}
While the bounding box mesh was uniform,
the triangular elements obtained by intersecting the 3D bounding box
mesh with the level set function can present arbitrarily irregular sizes 
and shapes (see also Fig.~\ref{fig-intersection-sphere}):
\begin{verbatim}
  geo -min-element-measure -max-element-measure sphere.geo
\end{verbatim}
Nevertheless, and surprisingly, 
\citet{OlsReuXu-2012} recently showed that
the finite element method converges on these irregular
intersected families of meshes.

%---------------------------------------
\myexamplenoinput{level_set_torus.cc}
%---------------------------------------
This approach can be extended to the Laplace-Beltrami problem
on a torus:
\begin{verbatim}
  sed -e 's/sphere/torus/' < level_set_sphere.cc > level_set_torus.cc
  make level_set_torus
  ./level_set_torus cube2.geo | geo -upgrade - > torus.geo
  geo torus.geo -stereo
  ./laplace_s torus.geo P1 | field -
\end{verbatim}
Note that the intersected mesh is also irregular:
\begin{verbatim}
  geo -min-element-measure -max-element-measure torus.geo
\end{verbatim}

\clearpage
% ---------------------------------------------------------
\subsection{The banded level set method}
% ---------------------------------------------------------
\cindex{method!level set!banded}%
The banded level set method presents the advantages of the
two previous methods without their drawback: it applies to
very general geometries, as described by a level set funtion,
and stronger convergence properties, as usual finite element methods.
The previous intersection mesh
can be circumvented by enlarging the surface $\Gamma_h$
to a band $\beta_h$ containing all the intersected elements
of ${\cal T}_h$ (see ~\citealp{OlsReuGra-2009,Abo-2010-m2r,Dic-2011-m2r}):
\[
  \beta_h = \{ K \in {\cal T}_h; K\cap\Gamma_h \neq \emptyset \}
\]
Then, we introduce $B_h$ the piecewise affine functional space over $\beta_h$:
\[
      B_h = \{ v \in L^2(\beta_h) \cap C^0(\beta_h); \ 
          v_{/K} \in P_1, \ 
          \forall K \in {\cal T}_h \}
\]
The problem is extended from $\Gamma_h$ to $\beta_h$ as:\\
\mbox{}\ \ $(VF)_h$: {\it find $u_h \in B_h$ such that:}
  $$
    a(u_h,v_h) = l(v_h), \ \forall v_h \in B_h
  $$
where, for all $u,v\in B_h$,
  \begin{eqnarray*}
    a(u,v) &=&  \int_{\Gamma_h} \left( u \, v + \nabla_s u . \nabla_s v\right) \, {\rm d}s \\
    l(v) &=& \int_{\Gamma_h} f \, v \, {\rm d}s
  \end{eqnarray*}
for all $u_h,v_h\in B_h$.
Note that while $u_h$ and $v_h$ are defined over $\beta_h$,
the summations in the variational formulations are restricted 
only to $\Gamma_h\subset \beta_h$.

% -----------------------------------------
\myexamplelicense{helmholtz_band_iterative.cc}
% -----------------------------------------

% -----------------------------------------
\subsubsection*{Comments}
% -----------------------------------------
\clindex{band}%
The band is build directly from the level set function as:
\begin{lstlisting}[numbers=none,frame=none]
  band gamma_h (phi_h);
\end{lstlisting}
The band structure is a small class that groups
the surface mesh $\Gamma_h$, available as \code{gamma_h.level_set()},
and the $\beta_h$ mesh, available as \code{gamma_h.band()}.
It also manages some correspondence between both meshes.
Then, the space of piecewise affine functions over the band is introduced:
\begin{lstlisting}[numbers=none,frame=none]
  space Bh (gamma_h.band(), "P1");
\end{lstlisting}
\cindex{function!\code{integrate}!on a band}%
Next, the bilinear form is computed by using the \code{integrate} function,
with the band \code{gamma_h} as a domain-like argument:
\begin{lstlisting}[numbers=none,frame=none]
  form a = integrate (gamma_h, u*v + dot(grad_s(u),grad_s(v)));
\end{lstlisting}
The right-hand side also admits the \code{gamma_h} argument:
\begin{lstlisting}[numbers=none,frame=none]
  field lh = integrate (gamma_h, f(d)*v);
\end{lstlisting}
Recall that summations for both forms and right-hand side
will be performed on $\Gamma_h$, represented
by \code{gamma_h.level_set()}, while the approximate functional space is $B_h$. 
\cindex{method!minres algorithm}%
\cindex{matrix!singular}%
Due to this summation on $\Gamma_h$ instead of $\beta_h$,
the matrix of the system is singular~\citep{OlsReuGra-2009,OlsReu-2010,Abo-2010-m2r}
and the MINRES algorithm has been chosen to solve the linear system:
\begin{lstlisting}[numbers=none,frame=none]
  minres (a.uu(), uh.set_u(), lh.u(), eye(), sopt);
\end{lstlisting}
\findex{diag}%
\clindex{eye}%
\cindex{matrix!diagonal}
\cindex{matrix!identity}
The \code{eye()} argument represents here the identity preconditioner, i.e. no preconditioner
at all.
It has few influence of the convergence properties of the matrix and
could be replaced by another simple one:
the diagonal of the matrix \code{diag(a.uu())}
without sensible gain of performance:
\begin{lstlisting}[numbers=none,frame=none]
  minres (a.uu(), uh.set_u(), lh.u(), diag(a.uu()), sopt);
\end{lstlisting}
See the reference manual for more about \code{minres}, 
e.g. on the \Rheolef\  web site or as unix manual
\clindex{reference manual}%
\clindex{man}%
\begin{verbatim}
  man minres
\end{verbatim}

%Finally, the $\beta_h$ meshe is saved:
%it will be required for the post-treatment of the solution.

% -----------------------------------------
\subsubsection*{How to run the program}
% -----------------------------------------
\begin{figure}[htb]
  \mbox{}\hspace{-1cm}
  \begin{tabular}{ccc}
    \includegraphics[height=6.0cm]{banded-level-set-circle-geo.png} &
    \includegraphics[height=6.0cm]{banded-level-set-circle-field.png} & \\
    \includegraphics[height=7.0cm]{banded-level-set-sphere-geo.png} &
    \includegraphics[height=6.0cm]{banded-level-set-sphere-field.png} & \\
    \includegraphics[height=4.5cm]{banded-level-set-torus-geo.png} &
    \includegraphics[height=4.5cm]{banded-level-set-torus-field.png} &
    \stereoglasses
  \end{tabular}
  \caption{The banded level set method:
	(top) circle;
        (center) sphere;
        (bottom) torus.
  }
  \label{fig-band-sphere}
\end{figure}

The compilation and run writes:
\begin{verbatim}
  make helmholtz_band_iterative
  mkgeo_grid -T 20 -a -2 -b 2 -c -2 -d 2 -f -2 -g 2 > cube-20.geo
  ./helmholtz_band_iterative cube-20.geo > sphere-band.field
\end{verbatim}
The run generates also two meshes (see Fig.~\ref{fig-band-sphere}):
the intersection mesh and the band around it.
The solution is here defined on this band: this extension has
no interpretation in terms of the initial problem and can be
restricted to the intersection mesh for visualization purpose:
\begin{verbatim}
  make proj_band
  ./proj_band < sphere-band.field | field -
\end{verbatim}
The \reffile{proj_band.cc} is presented below.
The run generates also the $\Gamma_h$ mesh (see Fig.~\ref{fig-band-sphere}),
required for the visualization.
%
The two-dimensional case is obtained simply by replacing the 3D bounding box by a 2D one:
\pindexopt{field}{-elevation}%
\pindexopt{field}{-bw}%
\pindexopt{field}{-stereo}%
\begin{verbatim}
  mkgeo_grid -t 20 -a -2 -b 2 -c -2 -d 2 > square-20.geo
  ./helmholtz_band_iterative square-20.geo > circle-band.field
  ./proj_band < circle-band.field | field -paraview -
  ./proj_band < circle-band.field | field -elevation -bw -stereo -
\end{verbatim}
% -----------------------------------------
\myexamplelicense{proj_band.cc}
% -----------------------------------------

% -------------------------------------------------------------------
\subsection{Improving the banded level set method with a direct solver}
% -------------------------------------------------------------------
The iterative algorithm previously used for solving the linear
system is not optimal:
for 3D problems on a surface, the bidimensionnal connectivity 
of the sparse matrix suggests that a direct sparse factorization
would be much more efficient.

Recall that $\phi_h=0$ on $\Gamma_h$.
Thus, if $u_h \in B_h$ is solution of the problem, then
$u_h+\alpha \phi_{h|\beta_h}\in B_h$ is also solution
for any $\alpha\in \mathbb{R}$,
where $\phi_{h|\beta_h}\in B_h$ denotes the restriction
of the level set function $\phi_h\in X_h$ on the band $\beta_h$.
Thus there is multiplicity of solutions and the matrix of the problem is singular.
The direct resolution is still possible on a modified linear system with additional
constraints in order to recover the unicity of the solution.
We impose the constraint that the solution $u_h$ should be
othogonal to $\phi_{h|\beta_h}\in B_h$.
\begin{figure}[htb]
  \begin{center}
    \begin{tabular}{cc}
      \includegraphics[height=7.5cm]{banded-level-set-cc-geo.pdf} &
    \end{tabular}
  \end{center}
  \caption{The banded level set method:
	the band is composed of several connected components. 
  }
  \label{fig-band-cc}
\end{figure}
\cindex{mesh!connected components}%
In some special cases, the band is composed of several
connected components (see Fig.~\ref{fig-band-cc}):
this appends when a vertex of the bounding box mesh belongs to $\Gamma_h$.
In that case, the constraint should be expressed on each connected component.
Fig.~\ref{fig-band-cc} shows also the case when a full side of an element is
included in $\Gamma_h$: such an element of the band is called {\em isolated}.

% -----------------------------------------
\myexamplelicense{helmholtz_band.cc}
% -----------------------------------------

% -----------------------------------------
\subsubsection*{Comments}
% -----------------------------------------
The management of the special sides and vertices that are fully included in $\Gamma_h$
is perfomed by:
\begin{lstlisting}[numbers=none,frame=none]
  Bh.block ("isolated");
  Bh.unblock ("zero");
\end{lstlisting}

The addition of linear constraints is similar
to the \reffile{neumann-laplace.cc} code, as
presented in section~\ref{sec-neumann-laplace}:
\begin{lstlisting}[numbers=none,frame=none]
  form A = {{ a, trans(b)},
            { b,    0    }};
\end{lstlisting}
Here \code{b} is a \code{vector<field>},
i.e. a vector of linear constraints, one per connected component
of the band $\beta_h$.

% -----------------------------------------
\subsubsection*{How to run the program}
% -----------------------------------------

The commands are similar to the previous iterative implementation,
just replacing \code{helmholtz_band_iterative} by \code{helmholtz_band}.

This approach could be also adapted to the Laplace-Beltrami
problem on the torus.

% -----------------------------------------
\myexamplelicense{laplace_band.cc}
% -----------------------------------------

% -----------------------------------------
\subsubsection*{Comments}
% -----------------------------------------
The code is similar to the previous one \code{helmholtz_band.cc}.
Since the solution is defined up to a constant,
an additional linear constraint has to be inserted:
\[
	\int_{\Gamma_h} u_h \, {\rm d}x = 0
\]
This writes:
\begin{lstlisting}[numbers=none,frame=none]
  field c  = integrate (gamma_h, v);
  form  A = {{    a,  trans(b), c },
             {    b,     0,     0 },
             { trans(c), 0,     0 }};
\end{lstlisting}

% -----------------------------------------
\subsubsection*{How to run the program}
% -----------------------------------------
\begin{verbatim}
  make laplace_band
  mkgeo_grid -T 20 -a -2 -b 2 -c -2 -d 2 -f -2 -g 2 > cube-20.geo
  ./laplace_band cube-20.geo > torus-band.field
  ./proj_band < torus-band.field | field -stereo -
\end{verbatim}
% do no more work:
% geo cube-20.band.geo -stereo -cut
The solution is represented on Fig.~\ref{fig-band-sphere}.bottom.

\vfill