\section{Poisson equation}%
\label{S:poisson-impl}%
\idx{implementation of model problems!Poisson equation}%
\idx{Poisson equation!implementation}

In this section we describe a model implementation for the Poisson equation
\begin{alignat*}{2}
-\Delta u &= f &\qquad&\mbox{in } \Omega \subset \R^d,\\
        u &= g &&\text{on } \Gamma_d,\\
        \partial_\nu u + \alpha_r\,u &= g_n &&\text{on }
        \Gamma_n,\text{ with } \partial\Omega = \Gamma_d \dot\cup
        \Gamma_n.
\end{alignat*}
Apart from the slightly complicated boundary conditions this is the
most simple elliptic problem, but the program presents all major
ingredients for general scalar stationary problems. Also, Poisson
equations often occur as sub-problems in much more complicated
settings. Modifications needed for a nonlinear problem are presented
in \secref{S:nonlin-impl}.

\begin{figure}[htbp]
\centerline{\includegraphics[scale=0.35]{EPS/ellipt_uh}
            \quad \includegraphics[scale=0.35]{EPS/ellipt_mesh}}
\caption[Solution of the linear Poisson problem and 
         corresponding mesh]{Solution of the linear Poisson problem and 
         corresponding mesh. The pictures were produced by GRAPE.}
\label{F:ellipt}
\end{figure}

Data and parameters described below lead in 2d to the solution and mesh shown
in \figref{F:ellipt}. The implementation of the Poisson problem is split into
several major steps which are now described in detail.



\subsection{Include file and global variables}

All \ALBERTA source files must include the header file \code{alberta.h}
with all \ALBERTA type definitions, function prototypes and macro definitions:
\bv\begin{lstlisting}
#include <alberta.h>
\end{lstlisting}\ev
This is realised by including the header file \code{alberta-demo.h} which
additionally includes the header files \code{graphics.h} and
\code{geomview-graphics.h} for graphical output:
\bv\begin{lstlisting}
#include "alberta-demo.h"
\end{lstlisting}\ev
leads to:
\bv\begin{lstlisting}
#include <alberta.h>
#include "graphics.h"
#include "geomview-graphics.h"
\end{lstlisting}\ev
For the linear scalar elliptic problem we use four global pointers to
data structures holding the finite element space and components of the
linear system of equations. These are used in different subroutines where
such information cannot be passed via parameters.
\bv\begin{lstlisting}
static const FE_SPACE *fe_space;
static DOF_REAL_VEC   *u_h;
static DOF_REAL_VEC   *f_h;
static DOF_MATRIX     *matrix;
\end{lstlisting}\ev
\begin{descr}
\kitem{fe\_space} a pointer to the actually used finite element space;
                  it is initialized by the function 
                  \code{main()}, see \secref{S:ellipt_init_dofs};
\kitem{u\_h} a pointer to a DOF vector storing the coefficients of the 
             discrete solution;
                  it is initialized by the function 
                  \code{main()}
\kitem{f\_h} a pointer to a DOF vector storing the load vector;
                  it is initialized by the function 
                  \code{main()}
\kitem{matrix} a pointer to a DOF matrix storing the system matrix;
                  it is initialized by the function 
                  \code{main()}
\end{descr}
The data structure \code{FE\_SPACE} is explained in 
\secref{S:FES-data}, \code{DOF\_REAL\_VEC} in 
\secref{S:DOF_VEC}, and
\code{DOF\_MATRIX} in \secref{S:DOF_MATRIX}.  Details about DOF
administration \code{DOF\_ADMIN} can be found in \secref{S:DOF_ADMIN}
and about the data structure \code{MESH} for a finite element mesh in
\secref{S:mesh_data_structure}.

We use another set of three global variable which store information
about the boundary conditions in use:
\begin{lstlisting}
static REAL robin_alpha = -1.0;
static bool pure_neumann = false;
static BNDRY_FLAGS dirichlet_mask; /* bit-mask of Dirichlet segments */
\end{lstlisting}
\begin{descr}
  \kitem{robin\_alpha} The zero-order factor if Robin-boundary
  conditions are prescribed, see \secref{S:robin_bound}.
  %%
  \kitem{pure\_neumann} When prescribing Neumann boundary conditions on
  all parts of the boundary, then the solution is only determined up
  to an additive constant. In this case it is necessary to perform a
  mean-value ``correction'' before, e.g., computing an error in the
  $L^2$-norm.
  %%
  \kitem{dirichlet\_mask} Initialized by a call to
  \code{GET\_PARAMTER()} in the main-function. A bit-mask tagging
  boundary segments on which the discrete solution is subject to
  Dirichlet boundary conditions. See \secref{S:boundary} and
  \secref{S:dirichlet_bound}.
\end{descr}

\subsection{The main program for the Poisson equation}%
\label{S:ellipt_main}

The main program is very simple, it just includes the main steps needed to
implement any stationary problem. Special problem-dependent aspects are hidden
in other subroutines described below.
\smallskip

We first read a parameter file (indicating which data, algorithms, and
solvers should be used; the file is described below in
\secref{S:ellipt_par}). The call to \code{parse\_parameters()} is
further explained in the section \secref{S:parse_parameters} above.
The parameters fetched from the parameter file at this point in the
code are:
\begin{descr}
  \kitem{dim} Dimension of the mesh.
  %%
  \kitem{filename} The file-name for the macro-triangulation.
  %%
  \kitem{degree} The desired polynomial degree for the finite element
  triangulation (should be between $1$ and $4$).
  %%
  \kitem{n\_refine} The number of global refinements of the mesh to be
  performed before starting the simulation.
  %%
  \kitem{do\_graphics} A boolean value for disabling all graphical
  output (individual windows can be disabled separately, see below in
  \secref{S:ellipt_par}.
  %%
  \kitem{dirichlet\_bit} The number of the boundary segment where
  Dirichlet boundary conditions should be imposed. See also
  \secref{S:boundary}. Boundary segments having another number than
  \code{dirichlet\_bit} are Neumann- or Robin-boundaries.
  %%
  \kitem{robin\_alpha} If positive, the zero-order parameter for a
  Robin-boundary condition. If negative and no boundary segment is a
  Dirichlet-boundary, then the discrete right-hand side will be forced
  to obey the mean-value zero compatibility condition. See
  \secref{S:boundary_conditions}.
\end{descr}

Having fetched those basic parameters from the data file
\code{INIT/ellipt.dat} we read the macro triangulation and initialize
the mesh (the basic geometric data structure). The subdirectories
\code{Macro/} in the \code{alberta-VERSION-demo/src/*d/} directories
contain data for several sample macro triangulations. How to read and
write macro triangulation files is explained in \secref{S:macro_tria}.

Now that the domain's geometry is defined, we allocate standard
Lagrange basis functions and from them generate a finite element space
through a call to \code{get\_fe\_space()}. The mesh is globally
refined if necessary. A call to \code{graphics()} displays the initial
mesh, unless the parameter \code{do\_graphics} has been initialized to
\code{false}, in which case no graphical output at all will appear.

Afterwards, the DOF vectors \code{u\_h} and \code{f\_h}, and the DOF
matrix \code{matrix} are allocated. The vector \code{u\_h}
additionally is initialized with zeros and the function pointers for
an automatic interpolation during refinement and coarsening are
adjusted to the predefined functions in \code{fe\_space->bas\_fcts}.
The load vector \code{f\_h} and the system matrix \code{matrix} are
newly assembled on each call of \code{build()}. Thus, there is no need
for interpolation during mesh modifications or initialization.
%%
Additionally, we initialize the global variable
\code{dirichlet\_mask}, setting the bit \code{dirichlet\_bit} to mark
those parts of the boundary, which are subject to Dirichlet boundary
conditions. The variable \code{dirichlet\_mask} is later on used by
several other routines: for matrix assembly, to install Dirichlet
boundary conditions into the load vector, and during the computation
of the error estimate.

The basic algorithmic data structure \code{ADAPT\_STAT} introduced in
\secref{S:adapt_stat_in_ALBERTA} specifies the behaviour of the
adaptive finite element method for stationary problems.  A
pre-initialized data structure is accessed by the function
\code{get\_adapt\_stat()}; the most important members
(\code{adapt->tolerance}, \code{adapt->strategy}, etc.) are
automatically initialized with values from the parameter file; other
members can be also initialized by adding similar lines for these
members to the parameter file (compare \secref{S:get_adapt}).
Eventually, function pointers for the problem dependent routines have
to be set (\code{estimate}, \code{get\_el\_est}, \code{build},
\code{solve}).  Since the assemblage is done in one step after all
mesh modifications, only \code{adapt->build\_after\_coarsen} is used,
no assemblage is done before refinement or before coarsening. These
additional assemblage steps are possible and may be needed in a more
general application, for details see \secref{S:adapt_stat_in_ALBERTA}.

The adaptive procedure is started by a call of
\code{adapt\_method\_stat()}. This automatically solves the discrete
problem, computes the error estimate, and refines the mesh until the
given tolerance is met, or the maximal number of iterations is
reached, compare \secref{S:adapt_stat_in_ALBERTA}. Finally,
\code{WAIT\_REALLY} allows an inspection of the final solution by
preventing a direct program exit with closure of the graphics windows.
The \code{WAIT\_REALLY}-blocker is not necessary when using the
\code{gltools} package for the graphical output.

\bv\begin{lstlisting}
int main(int argc, char **argv)
{
  FUNCNAME("main");
  MACRO_DATA        *data;
  MESH              *mesh;
  int                n_refine = 0, dim, degree = 1, dirichlet_bit = 1;
  const BAS_FCTS    *lagrange;
  static ADAPT_STAT *adapt;
  char               filename[PATH_MAX];

  /*****************************************************************************
   * first of all, initialize the access to parameters of the init file
   ****************************************************************************/

  parse_parameters(argc, argv, "INIT/ellipt.dat");

  GET_PARAMETER(1, "mesh dimension", "%d", &dim);
  GET_PARAMETER(1, "macro file name", "%s", filename);
  GET_PARAMETER(1, "polynomial degree", "%d", &degree);
  GET_PARAMETER(1, "global refinements", "%d", &n_refine);
  GET_PARAMETER(1, "online graphics", "%B", &do_graphics);
  GET_PARAMETER(1, "dirichlet boundary", "%d", &dirichlet_bit);
  GET_PARAMETER(1, "robin factor", "%f", &robin_alpha);

  /*****************************************************************************
   *  get a mesh, and read the macro triangulation from file
   ****************************************************************************/
  data = read_macro(filename);
  mesh = GET_MESH(dim, "ALBERTA mesh", data,
                  NULL /* init_node_projection() */,
                  NULL /* init_wall_trafos() */);
  free_macro_data(data);

  init_leaf_data(mesh, sizeof(struct ellipt_leaf_data),
                 NULL /* refine_leaf_data() */,
                 NULL /* coarsen_leaf_data() */);


  /*****************************************************************************
   *  initialize the global variables shared across build(), solve()
   *  and estimate().
   ****************************************************************************/
  lagrange = get_lagrange(mesh->dim, degree);
  TEST_EXIT(lagrange, "no lagrange BAS_FCTS\n");
  fe_space = get_fe_space(mesh, lagrange->name, lagrange, 1 /* rdim */,
                          ADM_FLAGS_DFLT);

  global_refine(mesh, n_refine * mesh->dim, FILL_NOTHING);

  if (do_graphics) {
    MSG("Displaying the mesh.\n");
    graphics(mesh, NULL /* u_h */, NULL /* get_est()*/ , NULL /* u_exact() */,
             HUGE_VAL /* time */);
  }

  matrix = get_dof_matrix("A", fe_space, NULL /* col_fe_space */);
  f_h    = get_dof_real_vec("f_h", fe_space);
  u_h    = get_dof_real_vec("u_h", fe_space);
  u_h->refine_interpol = fe_space->bas_fcts->real_refine_inter;
  u_h->coarse_restrict = fe_space->bas_fcts->real_coarse_inter;
  dof_set(0.0, u_h);      /*  initialize u_h */

  if (dirichlet_bit > 0) {
    BNDRY_FLAGS_SET(dirichlet_mask, dirichlet_bit);
  }

  /*****************************************************************************
   *  init adapt structure and start adaptive method
   ****************************************************************************/
  adapt = get_adapt_stat(mesh->dim, "ellipt", "adapt", 2,
                         NULL /* ADAPT_STAT storage area, optional */);
  adapt->estimate = estimate;
  adapt->get_el_est = get_el_est;
  adapt->build_after_coarsen = build;
  adapt->solve = solve;

  adapt_method_stat(mesh, adapt);

  if (do_graphics) {
    MSG("Displaying u_h, u, (u_h-u) and the final estimate.\n");
    graphics(mesh, u_h, get_el_est, u, HUGE_VAL /* time */);
  }
  WAIT_REALLY;

  return 0;
}
\end{lstlisting}\ev

\subsection{The parameter file for the Poisson equation}%
\label{S:ellipt_par}

The following parameter file \code{INIT/ellipt.dat} is used for the
\code{ellipt.c} program in the 2d case:
\bv\begin{verbatim}
mesh dimension:       2
macro file name:      Macro/macro.amc
global refinements:   1
polynomial degree:    3

dirichlet boundary:   1 % type of the Dirichlet boundary segment,
                        % must correspond to the boundary types used
                        % used in the macro triangulation. Use a value
                        % <= 0 to disable Dirichlet boundary
                        % conditions. Neumann boundary conditions will
                        % hold for all boundary segments with a type
                        % different from the value specified here.
robin factor:        -1 % > 0: Robin b.c.

online graphics:      true % global gfx kill-switch

% graphic windows: solution, estimate, mesh, and error if size > 0
graphic windows:      400 400 400 400
% for gltools graphics you can specify the range for the values of
% discrete solution for displaying:  min max
% automatical scaling by display routine if min >= max
gltools range:  0.0 -1.0

solver:                2 % 1: BICGSTAB 2: CG 3: GMRES 4: ODIR 5: ORES
solver max iteration:  10000
solver restart:        10  %  only used for GMRES
solver tolerance:      1.e-8
solver info:           2   
solver precon:         2   % 0: no precon
                           % 1: diag precon
                           % 2: HB precon
                           % 3: BPX precon
                           % 4: SSOR, omega = 1.0, #iter = 3
                           % 5: SSOR, with control over omega and #iter
                           % 6: ILU(k)
precon ssor omega:    1.0  % for precon == 5
precon ssor iter:     1    % for precon == 5
precon ilu(k):        8    % for precon == 6

error norm:           1   % 1: H1_NORM, 2: L2_NORM
estimator C0:         0.1 % constant of element residual
estimator C1:         0.1 % constant of jump residual
estimator C2:         0.0 % constant of coarsening estimate

adapt->strategy:          1 % 0: no adaption / 1: GR / 2: MS / 3: ES / 4: GERS
adapt->tolerance:         1.e-4
adapt->MS_gamma:          0.5
adapt->MS_gamma_c:        0.1
adapt->ES_theta:          1.9
adapt->ES_theta_c:        0.2
adapt->GERS_theta_star:   0.6
adapt->GERS_nu:           0.1
adapt->GERS_theta_c:      0.1
adapt->coarsen_allowed:   1
adapt->max_iteration:     20
adapt->info:              8

WAIT: 1
\end{verbatim}\ev

The file \code{Macro/macro.amc} storing data about the macro
triangulation for $\Omega = (0,1)^d$ can be found in
\secref{S:macro_tria} for 2d and 3d. The \code{polynomial degree}
parameter selects the third order Lagrange elements. \code{dirichlet
  boundary} marks those parts of the boundary which are subject to
Dirichlet boundary conditions, see also \secref{S:boundary}. The value
of \code{dirichlet boundary} corresponds to the numbers assigned to
boundary segments in the macro-triangulation.

By \code{graphic windows}, the number and sizes of graphics output
windows are selected. This line is used by the \code{graphics()}
routine. For gltools graphics, the range of function values might be
specified (used for graph coloring and height). If no graphical output
at all is desired, then \code{online graphics} can be set to
\code{false}. Individual output windows can be disabled by setting
their size to $0$. The size is specified in units of screen pixels.

The solver for the linear system of equations is selected (here: the
conjugate gradient solver), and corresponding parameters like
preconditioner and tolerance. Some preconditioners need additional
parameters, these are specified here as well.

Parameters for the error estimator include values of different
constants and selection of the error norm to be estimated ($H^1$- or
$L^2$-norm, selection leads to multiplication with different powers of
the local mesh size in the error indicators), see
Section~\ref{S:ellipt_est}.

An error tolerance and selection of a marking strategy with
corresponding parameters are main data given to the adaptive method.
For the meaning of the individual parameters the reader is referred to
the conceptional \secref{book:S:refinement_strategies} and
\secref{book:S:coarsening_strategies} in the book-part of the manual, and
to \secref{book:S:adaptive_methods} which describes the implementation of
adaptive methods in \ALBERTA.

Finally, the \code{WAIT} parameter specifies whether the program
should wait for user interaction at additional breakpoints, whenever a
\code{WAIT} statement is executed as in the routine \code{graphics()},
for instance, in case the \code{gltools} package is not in use.
  
The solution and corresponding mesh in 2d for the above parameters are
shown in \figref{F:ellipt}.  As optimal parameter sets might differ
for different space dimensions, separate parameter files exist in
\code{1d/INIT/}, \code{2d/INIT/}, and \code{3d/INIT/}.

\subsection{Initialization of the finite element space}%
\label{S:ellipt_init_dofs}

In contrast to prior versions of \ALBERTA, finite element spaces may
be newly allocated at any time. Since this involves updating DOF
information on all elements, however, it is advisable to allocate
finite element spaces before refining a mesh, see also Sections
\ref{S:dof_access} and \ref{S:access_fe_space}. 

For the scalar elliptic problem we need one finite element space for
the discretization. In this example, we use Lagrange elements and we
initialize the degree of the elements via a parameter. The
corresponding \code{fe\_space} is initialized by
\code{get\_fe\_space()} which automatically stores at the mesh
information about the DOFs used by this finite element space.

It is possible to allocate several finite element spaces, for instance
in a mixed finite element method, compare \secref{S:access_fe_space}.

\subsection{Functions for leaf data}%
\label{S:ellipt_leaf_data}

As explained in \secref{S:leaf_data_info}, we can ``hide'' information
which is only needed on a leaf element at the pointer to the second
child.  Such information, which we use here, is the local error
indicator on an element.  For this elliptic problem we need one
\code{REAL} for storing this element indicator.

After mesh initialization by \code{GET\_MESH()} in the main program,
we have to give information about the size of leaf data to be stored
and how to transform leaf data from parent to children during
refinement and vice versa during coarsening.  The function
\code{init\_leaf\_data()} initializes the leaf data used for this
problem.  Here, leaf data is one
structure \code{struct ellipt\_leaf\_data} and no transformation
during mesh modifications is needed. The details of the
\code{LEAF\_DATA\_INFO} data structure are stated in
\secref{S:leaf_data_info}.
\bv\begin{lstlisting}
init_leaf_data(mesh, sizeof(struct ellipt_leaf_data),
               NULL /* refine_leaf_data() */,
               NULL /* coarsen_leaf_data() */);
\end{lstlisting}\ev
The error estimation is done by the library function
\code{ellipt\_est()}, see \secref{S:ellipt_est}.  
For \code{ellipt\_est()}, we need a function
which gives read and write access to the local element error, and for
the marking function of the adaptive procedure, we need a function
which returns the local error indicator, see \secref{S:adapt_stat_in_ALBERTA}.
The indicator is stored as the \code{REAL} member \code{estimate} of
\code{struct ellipt\_leaf\_data} and the function \code{rw\_el\_est()}
returns for each element a pointer to this member. The function
\code{get\_el\_est()} returns the value stored at that member for each
element.

\bv\begin{lstlisting}
struct ellipt_leaf_data
{
  REAL estimate;            /*  one real for the estimate */
};

static REAL *rw_el_est(EL *el)
{
  if (IS_LEAF_EL(el))
    return &((struct ellipt_leaf_data *)LEAF_DATA(el))->estimate;
  else
    return NULL;
}

static REAL get_el_est(EL *el)
{
  if (IS_LEAF_EL(el))
    return ((struct ellipt_leaf_data *)LEAF_DATA(el))->estimate;
  else
    return 0.0;
}
\end{lstlisting}\ev


\subsection{Data of the differential equation}%
\label{S:ellipt_data}

Data for the Poisson problem are the right hand side $f$ and
boundary values $g$.
For test purposes it is convenient to have access to an exact solution
of the problem. In this example we use the function
\[
u(x) = \mbox{e}^{-10\, |x|^2}
\] 
as exact solution, resulting in 
\[
\nabla u(x) = -20\,x\,\mbox{e}^{-10\, |x|^2}
\]
and
\[
f(x) = -\Delta u(x) = -(400\,|x|^2 - 20\,d)\,\mbox{e}^{-10\, |x|^2}.
\]
Here, $d$ denotes the space dimension, $\Omega\subset\R^d$.  The
functions \code{u()} and \code{grd\_u()} are the implementation of $u$
and $\nabla u$ and are optional (and usually not known for a general
problem). The functions \code{g()}, \code{gn()} and \code{f()} are
implementations of the boundary values and the right hand side and are
not optional. Of course, \code{g()} needs only to be implemented when
Dirichlet boundary conditions apply, likewise \code{gn()} only for
inhomogeneous Robin or Neumann boundary conditions (see
\secref{S:robin_bound} and \secref{S:neumann_bound}.

\bv\begin{lstlisting} 
#define GAUSS_SCALE 10.0

static REAL u(const REAL_D x)
{
  return exp(-GAUSS_SCALE*SCP_DOW(x,x));
}

static const REAL *grd_u(const REAL_D x, REAL_D grd)
{
  static REAL_D buffer;
  REAL          ux = exp(-GAUSS_SCALE*SCP_DOW(x,x));
  int           n;

  if (!grd) {
    grd = buffer;
  }

  for (n = 0;  n < DIM_OF_WORLD; n++)
    grd[n] = -2.0*GAUSS_SCALE*x[n]*ux;

  return grd;
}

/*******************************************************************************
 * problem data: right hand side, boundary values
 ******************************************************************************/

static REAL g(const REAL_D x) /* boundary values, not optional */
{
  return u(x);
}

static REAL gn(const REAL_D x, const REAL_D normal) /* Neumann b.c. */
{
  return robin_alpha > 0.0
    ? SCP_DOW(grd_u(x, NULL), normal) + robin_alpha * u(x)
    : SCP_DOW(grd_u(x, NULL), normal);
}

static REAL f(const REAL_D x) /* -Delta u, not optional        */
{
  REAL  r2 = SCP_DOW(x,x), ux  = exp(-GAUSS_SCALE*r2);
  return -(4.0*SQR(GAUSS_SCALE)*r2 - 2.0*GAUSS_SCALE*DIM_OF_WORLD)*ux;
}
\end{lstlisting}\ev
A common principle in the implementation of functions of the type
\code{grd\_u} is that we store the result either at the
caller-specified pointer \code{input}, if provided, or overwrite a
local static \code{buffer} on each call.

\subsection{The assemblage of the discrete system}%
\label{S:ellipt_build}

For the assemblage of the discrete system we use the tools described
in Sections \ref{S:matrix_assemblage}, \ref{S:vector_update}, and
\ref{S:dirichlet_bound}. For the matrix assemblage we have to provide
an element-wise description of the differential operator. Following
the description in Section~\ref{book:S:FE-dis-2nd} we provide the function
\code{init\_element()} for an initialization of the operator on an
element and the function \code{LALt()} for the computation of \,$\det
|D F_S| \Lambda A \Lambda^t$\, on the actual element, where $\Lambda$
is the Jacobian of the barycentric coordinates, $D F_S$ the the
Jacobian of the element parameterization, and $A$ the matrix of the
second order term. For $-\Delta$, we have $A = id$ and $\det |D F_S|
\Lambda\Lambda^t$ is the description of the complete differential
operator since no lower order terms are involved.

For passing information about the Jacobian $\Lambda$ of the
barycentric coordinates and $\det |D F_S|$ from the function
\code{init\_element()} to the function \code{LALt()} we use the data
structure \code{struct op\_data} which stores the Jacobian and the
determinant. The function \code{init\_element()} calculates the
Jacobian and the determinant by the library functions
\code{el\_grd\_lambda\_?d()} and the function \code{LALt()} uses these
values in order to compute $\det |D F_S| \Lambda \Lambda^t$. The mesh
dimension $d$ given by \code{mesh->dim} is always less than or equal
to the world dimension $n$ given by the macro \code{DIM\_OF\_WORLD},
hence we comment out irrelevant parts of the code.

Pointers to these functions and to one structure
\code{struct~op\_info} are members of a structure
\code{OPERATOR\_INFO} which is used for the initialization of a
function for the automatic assemblage of the global system matrix (see
also \exampleref{Ex:LALt} in \secref{S:matrix_assemblage} for the
access to a structure \code{matrix\_info}). For more general equations
with lower order terms, additional functions \code{Lb0}, \code{Lb1},
and/or \code{c} have to be defined at that point. This initialization
is done on the first call of the function \code{build()} which is
called by \code{adapt\_method\_stat()} during the adaptive cycle
(compare \secref{S:adapt_stat_in_ALBERTA}). 

By calling \code{dof\_compress()}, unused DOF indices are removed such
that the valid DOF indices are consecutive in their range. This
guarantees optimal performance of the BLAS1 routines used in the
iterative solvers and \code{admin->size\_used} is the dimension of the
current finite element space. This dimension is printed for
information.

On each call of \code{build()} the matrix is assembled by first
clearing the matrix using the function \code{clear\_dof\_matrix()} and
then adding element contributions by \code{update\_matrix()}. This
function will call \code{init\_element()} and \code{LALt()} on each
element.

The load vector \code{f\_h} is then initialized with zeros and the
right hand side is added by \code{L2scp\_fct\_bas()}. Finally, the
boundary conditions are installed into the load-vector, and possibly
also into the matrix in the case of Robin boundary conditions.
Dirichlet boundary values are also interpolated into the vector
\code{u\_h} for the discrete solution. If only Dirichlet boundary
conditions are desired, then the call to \code{boundary\_conditions()}
quoted below could be replaced by a less complicated call to
\code{dirichlet\_bound()}:
%%
\begin{lstlisting}
dirichlet_bound(f_h, u_h, NULL, dirichlet_mask, g);
\end{lstlisting}
%%
Analogously, if only inhomogeneous Neumann boundary conditions should
be implemented, then a call to \code{bndry\_L2scp\_fct\_bas()} could
replace the call to \code{boundary\_conditions()}. Compare Sections
\ref{S:vector_update}, \ref{S:dirichlet_bound},
\ref{S:neumann_bound}, \ref{S:robin_bound} and
\ref{S:boundary_conditions}.

\bv\begin{lstlisting}
struct op_data
{
  REAL_BD  Lambda; /*  the gradient of the barycentric coordinates */
  REAL    det;     /*  |det D F_S|                                 */
};

static
bool init_element(const EL_INFO *el_info, const QUAD *quad[3], void *ud)
{
  struct op_data *info = (struct op_data *)ud;

  /* ..._0cd: co-dimension 0 version of el_grd_lambda(dim, ...) */
  info->det = el_grd_lambda_0cd(el_info, info->Lambda);

  return false; /* not parametric */
}

static
const REAL_B *LALt(const EL_INFO *el_info, const QUAD *quad,
                   int iq, void *ud)
{
  static REAL_BB LALt;
  struct op_data *info = (struct op_data *)ud;
  int            i, j, dim = el_info->mesh->dim;

  for (i = 0; i < N_VERTICES(dim); i++) {
    LALt[i][i] = info->det*SCP_DOW(info->Lambda[i], info->Lambda[i]);
    for (j = i+1; j < N_VERTICES(dim); j++) {
      LALt[i][j] = SCP_DOW(info->Lambda[i], info->Lambda[j]);
      LALt[i][j] *= info->det;
      LALt[j][i] = LALt[i][j];
    }
  }

  return (const REAL_B *)LALt;
}

static void build(MESH *mesh, U_CHAR flag)
{
  FUNCNAME("build");
  static const EL_MATRIX_INFO *matrix_info;

  dof_compress(mesh);
  MSG("%d DOFs for %s\n", fe_space->admin->size_used, fe_space->name);

  if (!matrix_info) {
    /* information for matrix assembling (only once) */
    OPERATOR_INFO  o_info = { NULL, };
    static struct op_data user_data; /* storage for det and Lambda */

    o_info.row_fe_space    = o_info.col_fe_space = fe_space;
    o_info.init_element    = init_element;
    o_info.LALt.real       = LALt;
    o_info.LALt_pw_const   = true;        /* pw const. assemblage is faster */
    o_info.LALt_symmetric  = true;        /* symmetric assemblage is faster */
    BNDRY_FLAGS_CPY(o_info.dirichlet_bndry,
                    dirichlet_mask);          /* Dirichlet bndry conditions */
    o_info.user_data = (void *)&user_data;    /* application data */
    o_info.fill_flag = CALL_LEAF_EL|FILL_COORDS; /* only FILL_BOUND is added
                                                  * automatically.
                                                  */
    matrix_info = fill_matrix_info(&o_info, NULL);
  }

  /* assembling of matrix */
  clear_dof_matrix(matrix);
  update_matrix(matrix, matrix_info, NoTranspose);

  /* assembling of load vector */
  dof_set(0.0, f_h);
  L2scp_fct_bas(f, NULL /* quadrature */, f_h);

  /* Boundary values, the combination alpha_r < 0.0 flags automatic
   * mean-value correction iff f_h has non-zero mean-value and no
   * non-Neumann boundary conditions were detected during mesh
   * traversal.
   */
  pure_neumann =
    !boundary_conditions(matrix, f_h, u_h, NULL /* bound */,
                         dirichlet_mask,
                         g, gn,
                         robin_alpha, /* < 0: mean-value correction */
                         NULL /* wall_quad, use default */);
}
\end{lstlisting}\ev


\subsection{The solution of the discrete system}%
\label{S:ellipt_solve}

The function \code{solve()} computes the solution of the resulting
linear system. It is called by \code{adapt\_method\_stat()} (compare
\secref{S:adapt_stat_in_ALBERTA}). The system matrix for the Poisson
equation is positive definite and symmetric for non-Dirichlet DOFs.
Thus, the solution of the resulting linear system is rather easy and
we can use any preconditioned Krylov-space solver
(\code{oem\_solve\_s()}), compare \secref{S:ALBERTA_OEM_solvers}. On
the first call of \code{solve()}, the parameters for the linear solver
are initialized and stored in \code{static} variables.  For the
\textsf{OEM} solver we have to initialize the \code{solver}, the
tolerance \code{tol} for the residual, a maximal number of iterations
\code{max\_iter}, the level of information printed by the linear
solver, and the use of a preconditioner by the parameter \code{icon},
which may be \code{0} (no preconditioning), \code{1} (diagonal
preconditioning), \code{2} (hierarchical basis preconditioning),
\code{3} (BPX preconditioning), \code{4} (SSOR preconditioning, with
given \code{omega} = 1.0, \code{\#iter} = 3), \code{5} (SSOR
preconditioning, with control over \code{omega} and \code{\#iter} ),
or \code{6} (ILU(k) preconditioning). If GMRes is used, then the
dimension of the Krylov-space for the minimizing procedure is needed,
too. If ILU(k) is used, then the level $k$ is needed, too (ILU(k)
denotes the ILU-flavour described in \cite{templates:94}).

After solving the discrete system, the discrete solution (and mesh)
is displayed by calling \code{graphics()}.
\bv\begin{lstlisting}
static void solve(MESH *mesh)
{
  FUNCNAME("solve");
  static REAL tol = 1.e-8, ssor_omega = 1.0;
  static int  max_iter = 1000, info = 2, restart = 0;
  static int  ssor_iter = 1, ilu_k = 8;
  static OEM_PRECON icon = DiagPrecon;
  static OEM_SOLVER solver = NoSolver;
  const PRECON *precon;

  if (solver == NoSolver) {
    GET_PARAMETER(1, "solver", "%d", &solver);
    GET_PARAMETER(1, "solver tolerance", "%f", &tol);
    GET_PARAMETER(1, "solver precon", "%d", &icon);
    GET_PARAMETER(1, "solver max iteration", "%d", &max_iter);
    GET_PARAMETER(1, "solver info", "%d", &info);
    if (icon == __SSORPrecon) {
        GET_PARAMETER(1, "precon ssor omega", "%f", &ssor_omega);
        GET_PARAMETER(1, "precon ssor iter", "%d", &ssor_iter);
    }
    if (icon == ILUkPrecon)
        GET_PARAMETER(1, "precon ilu(k)", "%d", &ilu_k);
    if (solver == GMRes) {
      GET_PARAMETER(1, "solver restart", "%d", &restart);
    }
  }

  if (icon == ILUkPrecon)
    precon = init_oem_precon(matrix, NULL, info, ILUkPrecon, ilu_k);
  else
    precon = init_oem_precon(matrix, NULL, info, icon, ssor_omega, ssor_iter);
  oem_solve_s(matrix, NULL, f_h, u_h,
              solver, tol, precon, restart, max_iter, info);

  if (do_graphics) {
    MSG("Displaying u_h, u and (u_h-u).\n");
    graphics(mesh, u_h, NULL /* get_el_est */, u, HUGE_VAL /* time */);
  }

  return;
}
\end{lstlisting}\ev

\subsection{Error estimation}%
\label{S:ellipt_estimate}

The last ingredient missing for the adaptive procedure is a function
for an estimation of the error. For an elliptic problem with constant
coefficients in the second order term this can done by the library
function \code{ellipt\_est()} which implements the standard residual
type error estimator and is described in \secref{S:ellipt_est}.
\code{ellipt\_est()} needs a pointer to a function for writing the
local error indicators (the function \code{rw\_el\_est()} described
above in \secref{S:ellipt_leaf_data}) and a function \code{r()} for
the evaluation of the lower order terms of the element residuals at
quadrature nodes. For the Poisson equation, this function has to
return the negative value of the right hand side $f$ at that node
(which is implemented in \code{r()}). Since we only have to evaluate
the right hand side $f$, the init flag \code{r\_flag} is zero. For an
equation with lower order term involving the discrete solution or its
derivative this flag has to be \code{INIT\_UH} and/or
\code{INIT\_GRD\_UH}, if needed by \code{r()}, compare
Example~\ref{E:est-impl}. Finally, for inhomogeneous Neumann or Robin
boundary conditions we must pass a pointer to yet another function
\code{est\_gn()} to \code{ellipt\_est()} which describes the
inhomogeneity. The information about which boundaries are subject to
Dirichlet boundary conditions is provided through the bit-mask
\code{dirichlet\_mask}, which is passed to \code{ellipt\_est()},
compare \secref{S:boundary}.

The function \code{estimate()}, which is called by
\code{adapt\_method\_stat()}, first initializes parameters for the
error estimator, like the estimated norm and constants in front of the
residuals. On each call the error estimate is computed by
\code{ellipt\_est()}. The degrees for quadrature formulas are chosen
according to the degree of finite element basis functions.
Additionally, as the exact solution for our test problem is known
(defined by \code{u()} and \code{grd\_u()}), the true error between
discrete and exact solutions is calculated by the function
\code{H1\_err()} or \code{L2\_err()}, and the ratio of the true and
estimated errors is printed (which should be approximately constant).
The experimental orders of convergence of the estimated and exact
errors are calculated, which should both be, when using global
refinement with $d$ bisection refinements,
\code{fe\_space->bas\_fcts->degree} for the $H^1$ norm and
\code{fe\_space->bas\_fcts->degree+1} for the $L^2$ norm.  Finally,
the error indicators are displayed by calling \code{graphics()}.

\bv\begin{lstlisting}
static REAL r(const EL_INFO *el_info, const QUAD *quad, int iq,
              REAL uh_at_qp, const REAL_D grd_uh_at_qp)
{
  REAL_D x;

  coord_to_world(el_info, quad->lambda[iq], x);

  return -f(x);
}

static REAL est_gn(const EL_INFO *el_info,
                   const QUAD *quad,
                   int qp,
                   REAL uh_at_qp,
                   const REAL_D normal)
{
  /* we simply return gn(), exploiting the fact that the geometry cache
   * of the quadrature already contains the world-coordinates of the
   * quadrature points.
   */
  const QUAD_EL_CACHE *qelc =
    fill_quad_el_cache(el_info, quad, FILL_EL_QUAD_WORLD);

  if (robin_alpha > 0.0) {
    return gn(qelc->world[qp], normal) - robin_alpha * uh_at_qp;
  } else {
    return gn(qelc->world[qp], normal);
  }
}

#define EOC(e,eo) log(eo/MAX(e,1.0e-15))/M_LN2

static REAL estimate(MESH *mesh, ADAPT_STAT *adapt)
{
  FUNCNAME("estimate");
  static int   norm = -1;
  static REAL  C[3] = {1.0, 1.0, 0.0};
  static REAL  est_old = -1.0, err_old = -1.0;
  REAL         est, err;
  REAL_DD      A = {{0.0}};
  int          n;

  for (n = 0; n < DIM_OF_WORLD; n++) {
    A[n][n] = 1.0; /* set diagonal of A; all other elements are zero */
  }

  if (norm < 0) {
    norm = H1_NORM;
    GET_PARAMETER(1, "error norm", "%d", &norm);
    GET_PARAMETER(1, "estimator C0", "%f", &C[0]);
    GET_PARAMETER(1, "estimator C1", "%f", &C[1]);
    GET_PARAMETER(1, "estimator C2", "%f", &C[2]);
  }

  est = ellipt_est(u_h, adapt, rw_el_est, NULL /* rw_est_c() */,
                   -1 /* quad_degree */,
                   norm, C,
                   (const REAL_D *) A,
                   dirichlet_mask,
                   r, 0 /* (INIT_UH | INIT_GRD_UH), if needed by r() */,
                   est_gn, robin_alpha > 0.0 ? INIT_UH : 0);

  MSG("estimate   = %.8le", est);
  if (est_old >= 0)
    print_msg(", EOC: %.2lf\n", EOC(est,est_old));
  else
    print_msg("\n");
  est_old = est;

  if (norm == L2_NORM)
    err = L2_err(u, u_h, NULL /* quad */,
                 false /* relative error*/,
                 pure_neumann /* mean-value adjust */,
                 NULL /* rw_err_el()*/, NULL /* max_err_el2 */);
  else
    err = H1_err(grd_u, u_h, NULL /* quad */,
                 false /* relative error */,
                 NULL /* rw_err_el()*/, NULL /* max_err_el2 */);

  MSG("||u-uh||%s = %.8le", norm == L2_NORM ? "L2" : "H1", err);
  if (err_old >= 0)
    print_msg(", EOC: %.2lf\n", EOC(err,err_old));
  else
    print_msg("\n");
  err_old = err;
  MSG("||u-uh||%s/estimate = %.2lf\n", norm == L2_NORM ? "L2" : "H1",
      err/MAX(est,1.e-15));

  if (do_graphics) {
    MSG("Displaying the estimate.\n");
    graphics(mesh, NULL /* u_h */, get_el_est, NULL /* u_exact() */,
             HUGE_VAL /* time */);
  }

  return adapt->err_sum;
}
\end{lstlisting}\ev

%%% Local Variables: 
%%% mode: latex
%%% TeX-master: "alberta-man.tex"
%%% End: