\section{Nonlinear reaction--diffusion equation}%
\label{S:nonlin-impl}%
\idx{implementation of model problems!nonlinear reaction--diffusion equation}%
\idx{nonlinear reaction--diffusion equation!implementation}

In this section, we discuss the implementation of a stationary,
nonlinear problem. Due to the nonlinearity, the computation of the
discrete solution is more complex. The solver for the nonlinear
reaction--diffusion equation and the solver for Poisson equation,
described in Section~\ref{S:poisson-impl}, thus mainly differ in the
routines \code{build()} and \code{solve()}.

Here we describe the solution by a Newton method, which involves the
assemblage and solution of a linear system in each iteration. Hence,
we do not split the assemble and solve routines in \code{build()} and
\code{solve()} as in the solver for the Poisson equation (compare
Sections~\ref{S:ellipt_build} and \ref{S:ellipt_solve}), but only set
Dirichlet boundary values for the initial guess in \code{build()} and
solve the nonlinear equation (including the assemblage of linearized
systems) in \code{solve()}.  The actual solution process is
implemented by several subroutines in the separate file
\code{nlsolve.c}, see Sections~\ref{S:nonlin_build} and
\ref{S:nonlin_solve}.

Additionally we describe a simple way to handle different problem data easily,
see Sections \ref{S:nonlin_org} and \ref{S:nonlin_data}.
\medskip

We consider the following nonlinear reaction--diffusion equation:
\begin{subequations}\label{E:reac-diff}
\begin{alignat}{2}
  -k \Delta u + \sigma \, u^4  &= f + \sigma \,u_{ext}^4&
   \qquad &\mbox{in }\Omega \subset \R^d, \\
  u &= g  & &\mbox{on }\partial\Omega.
\end{alignat}
\end{subequations}
For $\Omega\subset \R^2$, this equation models the heat transport in a thin
plate $\Omega$ which radiates heat and is heated by an external heat source
$f$. Here, $k$ is the constant heat conductivity, $\sigma$ the
Stefan--Boltzmann constant, $g$ the temperature at the edges of the plate and
$u_{ext}$ the temperature of the surrounding space (absolute temperature in
${}^\circ\!K$).

The solver is applied to following data:
\begin{itemize}
\item[$\bullet$]
For testing the solver we again use the `exponential peak'
\[
   u(x) = e^{-10\,|x|^2}, \qquad x \in \Omega = (-1,1)^d,\, k=1,\,
   \sigma=1,\,   u_{ext} = 0.
\]
\item[$\bullet$] In general (due to the nonlinearity), the problem is
  not uniquely solvable; depending on the initial guess for
  the nonlinear solver at least two discrete solutions can be obtained
  by using data
\[
 \Omega = (0,1)^d,\, k=1,\,\sigma = 1,\, f\equiv 1,\, g\equiv 0,
\,u_{ext} = 0.
\]
and the interpolant of 
\[
  u_0(x) = 4^d\,U_0 \prod_{i=1}^{d} x_i (1 - x_i)
 \qquad \mbox{with } U_0 \in [-5.0, 1.0].
\]
as initial guess for the discrete solution on the coarsest grid.
\item[$\bullet$]
The last application now addresses a physical problem in 2d
with following data:
\begin{align*}
\Omega = (-1,1)^2,\, k=2,\,\sigma = 5.67\mbox{e-8},
\, g\equiv 300,\, u_{ext}= 273,\,
  f(x) =  
    \begin{cases}
      150, & \mbox{if } x \in (-\frac12,\frac12)^2,\\
        0, & \mbox{otherwise.}
    \end{cases}
\end{align*}
\end{itemize}

\subsection{Program organization and header file}\label{S:nonlin_org}

The implementation is split into three source files:
\begin{descr}
\kitem{nonlin.c} main program with all subroutines for the adaptive procedure;
                 initializes DOFs, leaf data and problem dependent data in 
                 \code{main()}
                 and the \code{solve()} routine calls the nonlinear solver;
\kitem{nlprob.c} definition of problem dependent data;
\kitem{nlsolve.c} implementation of the nonlinear solver.
\end{descr}
Data structures used in all source files, and prototypes of functions
are defined in the header file \code{nonlin.h}, which includes the
\code{alberta.h} header file on the first line. This file is included
by all three source files.
\bv\begin{verbatim}
typedef struct prob_data PROB_DATA;
struct prob_data
{
  MACRO_DATA   *data;
  REAL          k, sigma;

  REAL        (*g)(const REAL_D x);
  REAL        (*f)(const REAL_D x);

  REAL        (*u0)(const REAL_D x);

  REAL        (*u)(const REAL_D x);
  const REAL *(*grd_u)(const REAL_D x, REAL_D input);
};

/*--- file nlprob.c --------------------------------------------------------*/
const PROB_DATA *init_problem(MESH *mesh);

/*--- file nlsolve.c -------------------------------------------------------*/
int nlsolve(DOF_REAL_VEC *, REAL, REAL, REAL (*)(const REAL_D));
\end{verbatim}\ev
The data structure \code{PROB\_DATA} yields following information:
\begin{descr}
\kitem{data} pointer to a macro triangulation object;
\kitem{k} diffusion coefficient (constant heat conductivity);
\kitem{sigma} reaction coefficient (Stefan--Boltzmann constant);
\kitem{g} pointer to a function for evaluating boundary values;
\kitem{f} pointer to a function for evaluating the right-hand side
          ($f$ + $\sigma\,u_{ext}^4$);
\kitem{u0} pointer to a function for evaluating an initial guess for the
           discrete solution on the macro triangulation, if not \nil;
\kitem{u} pointer to a function for evaluating the true solution, if not \nil
          (only for test purpose); 
\kitem{grd\_u} pointer to a function for evaluating the gradient of the 
          true solution, if not \nil (only for test purpose).
\end{descr}

The function \code{init\_problem()} initializes problem data, like
boundary values, right hand side, etc. which is stored in a
\code{PROB\_DATA} structure and reads data of the macro triangulation
for the actual problem.  The function \code{nlsolve()} implements the
nonlinear solver by a Newton method including the assemblage and
solution of the linearized sub-problems.

\subsection{Global variables}
In the main source file for the nonlinear solver \code{nonlin.c}
we use the following global variables:
\bv\begin{verbatim}
#include "nonlin.h"

#include "alberta-demo.h"          /* proto-types for support functions */

static bool do_graphics = true;    /* global graphics switch          */

static const FE_SPACE  *fe_space;  /* initialized by init_dof_admin() */
static DOF_REAL_VEC    *u_h;       /* initialized by build()          */
static const PROB_DATA *prob_data; /* initialized by main()           */
static BNDRY_FLAGS dirichlet_mask; /* bit-mask for Dirichlet segments */
\end{verbatim}\ev
As in the solver for the linear Poisson equation, we have a pointer to
the used \code{fe\_space} and the discrete solution \code{u\_h}.  In
this file, we do not need a pointer to a \code{DOF\_MATRIX} for
storing the system matrix and a pointer to a \code{DOF\_REAL\_VEC} for
storing the right hand side. The system matrix and right hand side are
handled by the nonlinear solver \code{nlsolve()}, implemented in
\code{nlsolve.c}. Data about the problem is handled via the
\code{prob\_data} pointer. The variable \code{dirichlet\_mask} marks
those segments on which Dirichlet boundary conditions are imposed, see
\secref{S:boundary}. It is initialized by the \code{main()} function.

\subsection{The main program for the nonlinear reaction--diffusion equation}%
\label{S:nonlin_main}

The main program is very similar to the main program of the Poisson
problem described in Section~\ref{S:ellipt_main}. 

After initializing the access to the parameter file and processing
command-line parameters (see \secref{S:parse_parameters}), the mesh
with the used leaf data is initialized, problem dependent data,
including the macro triangulation, are initialized by
\code{init\_problem(mesh)} (see Section~\ref{S:nonlin_data}), a finite
element space is allocated, the structure for the adaptive method is
filled, and finally the adaptive method is started.
%%
\bv\begin{verbatim} int main(int argc, char **argv)
 { FUNCNAME("main"); MESH *mesh; const BAS_FCTS *lagrange; ADAPT_STAT
 *adapt; int dim, degree = 1, n_refine;

  /*****************************************************************************
   * first of all, initialize the access to parameters of the init file
   ****************************************************************************/
  parse_parameters(argc, argv, "INIT/nonlin.dat");

  GET_PARAMETER(1, "global refinements", "%d", &n_refine);
  GET_PARAMETER(1, "polynomial degree", "%d", &degree);
  GET_PARAMETER(1, "mesh dimension", "%d", &dim);
  GET_PARAMETER(1, "online graphics", "%d", &do_graphics);

  BNDRY_FLAGS_ALL(dirichlet_mask); /* Only Dirichlet b.c. supported here */
  /*****************************************************************************
   * init problem dependent data and read macro triangulation
   ****************************************************************************/

  prob_data = init_problem();

  /*****************************************************************************
   *  get a mesh with DOFs and leaf data                               
   ****************************************************************************/

  mesh = GET_MESH(dim,"Nonlinear problem mesh", prob_data->data, NULL, NULL);

  free_macro_data(prob_data->data);

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

  lagrange = get_lagrange(mesh->dim, degree);
  TEST_EXIT(lagrange, "no lagrange BAS_FCTS\n");

  fe_space = get_fe_space(mesh, lagrange->name, lagrange, 1, ADM_FLAGS_DFLT);

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

  /*****************************************************************************
   *  init adapt structure and start adaptive method
   ****************************************************************************/
  adapt = get_adapt_stat(dim, "nonlin", "adapt", 1, NULL);
  adapt->estimate = estimate;
  adapt->get_el_est = get_el_est;
  adapt->build_after_coarsen = build;
  adapt->solve = solve;

  adapt_method_stat(mesh, adapt);

  WAIT_REALLY;

  return 0;
}
\end{verbatim}\ev

\subsection{Initialization of leaf data}%
\label{S:nonlin_dofs_leaf}

The functions for initializing leaf data (\code{init\_leaf\_data()}),
and for accessing leaf data (\code{rw\_el\_est()},
\code{get\_el\_est()}) are exactly the same as in the solver for the
linear Poisson equation, compare \secref{S:ellipt_leaf_data}.

\subsection{The build routine}%
\label{S:nonlin_build}

As mentioned above, inside the build routine we only access one vector
for storing the discrete solution. On the coarsest grid, the
discrete solution is initialized with zeros, or by interpolating
the function \code{prob\_data->u0}, which implements an initial 
guess for the discrete solution. On a refined grid we do not
initialize the discrete solution again. Here, we use the discrete
solution from the previous step, which is interpolated during
mesh modifications, as an initial guess.

In each adaptive cycle, Dirichlet boundary values are set for
the discrete solution. This ensures $u_0 \in g_h + \Xc_h$ for
the initial guess of the Newton method.
\bv\begin{verbatim}
static void build(MESH *mesh, U_CHAR flag)
{
  FUNCNAME("build");

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

  if (!u_h)                 /*  access and initialize discrete solution    */
  {
    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;
    if (prob_data->u0)
      interpol(prob_data->u0, u_h);
    else
      dof_set(0.0, u_h);
  }

  /* set boundary values */
  dirichlet_bound(u_h, NULL, NULL, dirichlet_mask, prob_data->g);

  return;
}
\end{verbatim}\ev

\subsection{The solve routine}%
\label{S:nonlin_solve}

The \code{solve()} routine solves the nonlinear equation by calling
the function \code{nlsolve()} which is implemented in \code{nlsolve.c}
and described below in \secref{S:nonlin_solver}. 
After solving the discrete problem, the 
new discrete solution and true error is displayed via the 
\code{graphics()} routine. The true error can be computed only for
the first application, where the true solution is known
(\code{prob\_data->u()} and \code{prob\_data->grd\_u()} are not \nil).
\bv\begin{verbatim}
static void solve(MESH *mesh)
{
  nlsolve(u_h, prob_data->k, prob_data->sigma, prob_data->f, dirichlet_mask);

  if (do_graphics) {
    graphics(mesh, u_h, NULL, prob_data->u, HUGE_VAL /* time */);
  }

  return;
}
\end{verbatim}\ev

\subsection{The estimator for the nonlinear problem}%
\label{S:nonlin_est}

In comparison to the Poisson program, the function \code{r()} which implements
the lower order term in the element residual changes due to the term $\sigma
u^4$ in the differential operator, compare \secref{S:ellipt_est}.  The right
hand side $f + \sigma u_{ext}^4$ is already implemented in the function
\code{prob\_data->f()}.

In the function \code{estimate()} we have to initialize the diagonal
of \code{A} with the heat conductivity \code{prob\_data->k} and for
the function \code{r()} we need the values of $u_h$ at the quadrature
node, thus \code{r\_flag = INIT\_UH} is set. The initialization of
parameters for the estimator is the same as in
\secref{S:ellipt_estimate}. Finally, the error indicator is displayed
by \code{graphics()}.

\bv\begin{verbatim}
static REAL r(const EL_INFO *el_info, const QUAD *quad, int iq, REAL uh_iq,
              const REAL_D grd_uh_iq)
{
  REAL_D      x;
  REAL        uhx2 = SQR(uh_iq);

  coord_to_world(el_info, quad->lambda[iq], x);
  return(prob_data->sigma*uhx2*uhx2 - (*prob_data->f)(x));
}

#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     degree, norm = -1;
  static REAL    C[3] = {1.0, 1.0, 0.0};
  static REAL    est, est_old = -1.0, err = -1.0, err_old = -1.0;
  static REAL    r_flag = INIT_UH;
  REAL_DD        A = {{0.0}};
  int            n;

  for (n = 0; n < DIM_OF_WORLD; n++)
    A[n][n] = prob_data->k;  /* set diagonal of A; other elements are zero  */

  if (norm < 0)
  {
    norm = H1_NORM;
    GET_PARAMETER(1, "error norm", "%d", &norm);
    GET_PARAMETER(1, "estimator C0", "%f", C);
    GET_PARAMETER(1, "estimator C1", "%f", C+1);
    GET_PARAMETER(1, "estimator C2", "%f", C+2);
  }
  degree = 2*u_h->fe_space->bas_fcts->degree;
  est = ellipt_est(u_h, adapt, rw_el_est, NULL, degree, norm, C,
                   (const REAL_D *) A, r, r_flag);

  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  &&  prob_data->u)
    err = L2_err(prob_data->u, u_h, NULL, 0, NULL, NULL);
  else if (norm == H1_NORM  &&  prob_data->grd_u)
    err = H1_err(prob_data->grd_u, u_h, NULL, 0, NULL, NULL);

  if (err >= 0)
  {
    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) {
    graphics(mesh, NULL, get_el_est, NULL, HUGE_VAL /* time */);
  }

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

\subsection{Initialization of problem dependent data}%
\label{S:nonlin_data}

The file \code{nlprob.c} contains all problem dependent data. On the
first line, \code{nonlin.h} is included and then two variables for
storing the values of the heat conductivity and the Stefan --Boltzmann
constant are declared. These values  are used by several functions:
\bv\begin{verbatim}
#include "nonlin.h"
static REAL  k = 1.0, sigma = 1.0;
\end{verbatim}\ev

The following functions are used in the first example for testing the
nonlinear solver (\code{problem number: 0}):
\bv\begin{verbatim}
static REAL u_0(const REAL_D x)
{
  REAL   x2 = SCP_DOW(x,x);
  return(exp(-10.0*x2));
}

static const REAL *grd_u_0(const REAL_D x, REAL_D input)
{
  static REAL_D buffer = {};
  REAL         *grd = input ? input : buffer;
  REAL          ux = exp(-10.0*SCP_DOW(x,x));
  int           n;

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

  return(grd);
}

static REAL f_0(const REAL_D x)
{
  REAL  r2 = SCP_DOW(x,x), ux  = exp(-10.0*r2), ux4 = ux*ux*ux*ux;
  return(sigma*ux4 - k*(400.0*r2 - 20.0*DIM_OF_WORLD)*ux);
}
\end{verbatim}\ev

For the computation of a stable and an unstable (but non-physical)
solution, depending on the initial choice of the discrete solution, the
following functions are used, which also use a global variable
\code{U0}.  Such an unstable solution in 3d is shown in
Figure~\ref{F:nl-sol-grid3d}.  
Data is given as follows (\code{problem number: 1}):
\bv\begin{verbatim}
static REAL U0 = 0.0;

static REAL g_1(const REAL_D x)
{
#if DIM_OF_WORLD == 1
  return(4.0*U0*x[0]*(1.0-x[0]));
#endif
#if DIM_OF_WORLD == 2
  return(16.0*U0*x[0]*(1.0-x[0])*x[1]*(1.0-x[1]));
#endif
#if DIM_OF_WORLD == 3
  return(64.0*U0*x[0]*(1.0-x[0])*x[1]*(1.0-x[1])*x[2]*(1.0-x[2]));
#endif
}

static REAL f_1(const REAL_D x)
{
  return(1.0);
}
\end{verbatim}\ev

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.48\hsize]{EPS/nl-sol3d}\hfill%
\includegraphics[width=0.48\hsize]{EPS/nl-grid3d}
\end{center}
\caption[Graph of the unstable solution, nonlinear reaction-diffusion
problem]{Graph of the unstable solution with corresponding mesh of the
  nonlinear reaction-diffusion problem in 3d on the clipping plane
  $z=0.5$.  The pictures were produced by the gltools.}
\label{F:nl-sol-grid3d} 
\end{figure}

The last example needs functions for boundary data and right hand side
and variables for the temperature at the edges, and $\sigma\,u_{ext}^4$.
A solution to this problem is depicted in Figure~\ref{F:nl-sol-grid2d}
and problem data is (\code{problem number: 2}):
\bv\begin{verbatim}
static REAL g2 = 300.0, sigma_uext4 = 0.0;
static REAL g_2(const REAL_D x)
{
  return(g2);
}

static REAL f_2(const REAL_D x)
{
  if (x[0] >= -0.25  &&  x[0] <= 0.25  &&  x[1] >= -0.25  &&  x[1] <= 0.25)
    return(150.0 + sigma_uext4);
  else
    return(sigma_uext4);
}
\end{verbatim}\ev
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.48\hsize]{EPS/nl-sol2d}\hfill%
\includegraphics[width=0.48\hsize]{EPS/nl-grid2d}
\end{center}
\caption[Graph of the physical solution, nonlinear
reaction-diffusion problem]{Graph of the solution to the physical
  problem with corresponding mesh of the nonlinear reaction-diffusion
  problem in 2d. The pictures were produced by the gltools.}
\label{F:nl-sol-grid2d} \end{figure}

Depending on the chosen problem via the parameter \code{problem
number}, the function \code{init\_problem()} initializes the entries
of a \code{PROB\_DATA} structure, adjusts the corresponding function
pointers, reads the macro triangulation, and returns a pointer to the
filled \code{PROB\_DATA} structure. Information stored in
\code{PROB\_DATA} is then used in the \code{build()} and the
\code{nlsolve()} routines.
\bv\begin{verbatim}
const PROB_DATA *init_problem(MESH *mesh)
{
  FUNCNAME("init_problem");
  static PROB_DATA prob_data;
  int              pn = 2;

  GET_PARAMETER(1, "problem number", "%d", &pn);
  switch (pn)
  {
  case 0:   /*---  problem with known true solution  -----------------------*/
    k = 1.0;
    sigma = 1.0;

    prob_data.g = u_0;
    prob_data.f = f_0;

    prob_data.u = u_0;
    prob_data.grd_u = grd_u_0;

    prob_data.data = read_macro("Macro/macro-big.amc");
    break;
  case 1:   /*---  problem for computing a stable and an unstable sol.  ----*/
    k = 1.0;
    sigma = 1.0;

    prob_data.g = g_1;
    prob_data.f = f_1;

    prob_data.u0 = g_1;
    GET_PARAMETER(1, "U0", "%f", &U0);

    prob_data.data = read_macro("Macro/macro.amc");
    break;
  case 2:   /*---  physical problem  ---------------------------------------*/
    k = 2.0;
    sigma = 5.67e-8;
    sigma_uext4 = sigma*273*273*273*273;

    prob_data.g = g_2;
    prob_data.f = f_2;
    prob_data.data = read_macro("Macro/macro-big.amc");
    break;
  default:
    ERROR_EXIT("no problem defined with problem no. %d\n", pn);
  }
  prob_data.k = k;
  prob_data.sigma = sigma;

  return &prob_data;
}
\end{verbatim}\ev

\subsection{The parameter file for the nonlinear reaction--diffusion equation}%
\label{S:nonlin_param}

The following parameter file \code{INIT/nonlin.dat} is read by \code{main()} 
for 2d.
\bv\begin{verbatim}
mesh dimension:       2
problem number:       2
global refinements:   1
polynomial degree:    2

online graphics:      false

U0:     -5.0          % height of initial guess for Problem 1

% graphic windows: solution, estimate, mesh, and error if size > 0
graphic windows:         500 500 0 0
% 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:  1.0 0.0

newton tolerance: 1.e-6     % tolerance for Newton
newton max. iter: 50        % maximal number of iterations of Newton
newton info:      6         % information level of Newton
newton restart:   10        % number of iterations for step size control

linear solver max iteration:  1000
linear solver restart:        10  %  only used for GMRES
linear solver tolerance:      1.e-8
linear solver info:           0   
linear solver precon:         2   % 0: no precon 1: diag precon
                                  % 2: HB precon 3: BPX precon

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:        2   % 0: no adaption 1: GR 2: MS 3: ES 4:GERS
adapt->tolerance:       1.e-2
adapt->MS_gamma:        0.5
adapt->max_iteration:   15
adapt->info:            4

WAIT: 1
\end{verbatim}\ev

Besides the parameters for the Newton solver and the height of the
initial guess \code{U0} in Problem 1, the file is very similar to the 
parameter file \code{ellipt.dat} for the Poisson problem,
compare Section~\ref{S:ellipt_par}. As mentioned above, additional
parameters may be defined or overwritten by command line
arguments, see Section~\ref{S:nonlin_main}.


\subsection{Implementation of the nonlinear solver}%
\label{S:nonlin_solver}

In this section, we now describe the solution of the nonlinear problem
which differs most from the solver for the Poisson equation.  It is
the last module missing for the adaptive solver.  We use the abstract
Newton methods of Section~\ref{S:nls} for solving
\[
u_h \in g_h + \Xc_h: \qquad F(u_h) = 0 \qquad \mbox{in } \Xc_h^*,
\]
where $g_h \in X_h$ is an approximation to boundary data $g$.  Using
the classical Newton method, we start with an initial guess $u_0\in
g_h+\Xc_h$, where Dirichlet boundary values are set in the
\code{build()} routine (compare Section~\ref{S:nonlin_build}).  For $m
\geq 0$ we compute
\[
d_m \in \Xc_h: \qquad DF(u_m) d_m = F(u_m) \qquad \mbox{in } \Xc_h^*
\]
and set
\[
u_{m+1} = u_m - d_m
\]
until some suitable norm $\|d_m\|$ or $\|F(u_{m+1})\|$ is sufficiently
small. Since the correction $d_m$ satisfies $d_m\in\Xc_h$, all Newton
iterates $u_m$ satisfy $u_m \in g_h+\Xc_h$, $m \geq 0$.
Newton methods with step size control solve similar defect equations and
perform similar update steps, compare \secref{S:nls}.

For $v\in g_h+\Xc_h$ the functional $F(v)\in\Xc_h^*$ 
of the nonlinear reaction--diffusion equation is defined by
\begin{equation}\label{E:F(v)}
\ldual{F(v)}{\varphi_j}{\Xc_h^* \times \Xc_h} =
 \int_\Omega k \nabla \varphi_j \nabla v  + \sigma\, \varphi_j\, v^4\,dx-
 \int_\Omega (f  + u_{ext}^4) \varphi_j \,dx
\qquad \mbox{for all } \varphi_j \in \Xc_h,
\end{equation}
and the Frechet derivative $DF(v)$ of $F$ is given by
\begin{equation}\label{E:DF(v)}
\ldual{DF(v)\,\varphi_i}{\varphi_j}{\Xc_h^* \times \Xc_h} =
 \int_\Omega k \nabla \varphi_j \nabla \varphi_i + 
  4 \sigma\, v^3\, \varphi_j\, \varphi_i \,dx 
\qquad \mbox{for all } \varphi_i,\, \varphi_j\in \Xc_h.
\end{equation}

The Newton solvers need a function for assembling the right hand side
vector of the discrete system \mathref{E:F(v)}, and the system matrix
of the linearized equation \mathref{E:DF(v)} for some given $v$ in
$X_h$. The system matrix is always symmetric. It is positive definite,
if $v \geq 0$, and is then solved by the conjugate gradient method.
For $v\not\geq 0$ BiCGStab is used. We choose the $H^1$ semi--norm
as problem dependent norm $\|.\|$.

\subsubsection{Problem dependent data structures for assembling and solving}

Similar to the assemblage of the system matrix for the Poisson problem,
we define a data structure \code{struct op\_info} in order to
pass information to the routines which describe the differential operator.
In the assembling of the linearized system around a given 
finite element function $v$ we additionally need the 
diffusion coefficient $k$ and reaction coefficient $\sigma$.
In general, $v$ is not constant on the elements, thus we have
to compute the zero order term by numerical quadrature on each
element. For this we need access to the used quadrature for this
term, and a vector storing the values of $v$ for all quadrature nodes.
\bv\begin{verbatim}
struct op_info
{
  REAL_BD  Lambda;          /*  the gradient of the barycentric coordinates */
  REAL    det;              /*  |det D F_S|                                 */

  REAL    k, sigma;         /*  diffusion and reaction coefficient          */

  const QUAD_FAST *quad_fast;  /*  quad_fast for the zero order term        */
  const REAL      *v_qp;       /*  v at all quadrature nodes of quad_fast   */
};
\end{verbatim}\ev

The general Newton solvers pass data about the actual problem by 
\code{void} pointers to the problem dependent routines. Information that
is used by these routines are collected in the data structure
\code{NEWTON\_DATA}
\bv\begin{verbatim}
typedef struct newton_data NEWTON_DATA;
struct newton_data
{
  const FE_SPACE   *fe_space;     /* used finite element space              */
  BNDRY_FLAGS      dirichlet_mask;

  REAL       k;                   /* diffusion coefficient                  */
  REAL       sigma;               /* reaction coefficient                   */
  REAL       (*f)(const REAL_D);  /* for evaluation f + sigma u_ext^4       */

  DOF_MATRIX *DF;                 /* pointer to system matrix               */

/*---  parameters for the linear solver  -----------------------------------*/
  OEM_SOLVER solver;              /* used solver: CG (v >= 0) else BiCGStab */
  REAL       tolerance;
  REAL	     ssor_omega;
  int        max_iter;
  int	     ssor_iter;
  int        ilu_k;
  int        restart;
  int        info;
  OEM_PRECON icon;
  const PRECON *precon;
};
\end{verbatim}\ev
All entries of this structure besides \code{solver} are initialized
in the function \code{nlsolve()}. The entry \code{solver} is set
every time the linearized matrix is assembled.

\subsubsection{The assembling routine}

\newcommand{\Nb}{N\vphantom{\Nc}}
Denote by $\{\varphi_0,\dots,\varphi_{\Nc}\}$ the basis of $\Xc_h$, by
$\{\varphi_0,\dots,\varphi_{\Nb}\}$ the basis of $X_h$. Let $\boldsymbol{A}$
be the stiffness matrix, i.e.
\[
A_{ij} = 
\begin{cases}
\int_\Omega k \nabla \varphi_j \nabla \varphi_i\,dx &
i = 0,\dots,\Nc,\, j = 0,\dots,\Nb,\\
\delta_{ij} & i = \Nc+1,\dots,\Nb,\, j = 0,\dots,\Nb,
\end{cases}
\]
and $\boldsymbol{M} = \boldsymbol{M}(v)$ the mass matrix, i.e.
\[
M_{ij} = 
\begin{cases}
\int_\Omega \sigma\, v^3\, \varphi_j\, \varphi_i\,dx &
 i = 0,\dots,\Nc,\, j = 0,\dots,\Nb,\\
0 & i = \Nc+1,\dots,\Nb,\, j = 0,\dots,\Nb.
\end{cases}
\]
The system matrix $\boldsymbol{L}$, representing $DF(v)$,
of the linearized equation is then given as
\[
\boldsymbol{L} = \boldsymbol{A} + 4 \boldsymbol{M}.
\]
The right hand side vector $\boldsymbol{F}$, representing $F(v)$ is
for all non--Dirichlet DOFs $j$ given by
\begin{align}\label{E:ApMv}
F_j &= 
\int_\Omega k \nabla v \nabla \varphi_j  + \sigma\, v^4\, \varphi_j \,dx
- \int_\Omega (f + \sigma u_{ext}^4) \varphi_j  \,dx\notag \\
&= (\boldsymbol{A}\,\boldsymbol{v} +  \boldsymbol{M}\,\boldsymbol{v})_j
- \int_\Omega (f + \sigma u_{ext}^4) \varphi_j  \,dx,
\end{align}
where $\boldsymbol{v}$ denotes the coefficient vector of $v$. Thus, we
want to use information assembled into $\boldsymbol{A}$ and
$\boldsymbol{M}$ for both system matrix and right hand side vector.

Unfortunately, this can not be done \emph{after} assembling
$\boldsymbol{A} + 4\, \boldsymbol{M}$ into the system matrix
$\boldsymbol{L}$ due to the different scaling of $\boldsymbol{M}$ in
the system matrix (factor $4$) and right hand side (factor $1$).
Storing both matrices $\boldsymbol{A}$ \emph{and} $\boldsymbol{M}$ is
too costly, since matrices are the objects in finite element codes
which need most memory.

The solution to this problem comes from the observation, that
\mathref{E:ApMv} holds also element--wise for the element contributions
of the right hand side and element matrices $\boldsymbol{A}_S$ and 
$\boldsymbol{M}_S$ when replacing $\boldsymbol{v}$ by the local 
coefficient vector $\boldsymbol{v}_S$.
Hence, on elements $S$ we compute the element contributions of 
$\boldsymbol{A}_S$ and $\boldsymbol{M}_S$, add them to the system
matrix, and use them and the local coefficient vector
$\boldsymbol{v}_S$ for adding the right hand side contribution
to the load vector.

The resulting assembling routine is more complicated in comparison to
the very simple routine used for the linear Poisson problem. On the
other hand, using \ALBERTA routines for the computation of element
matrices, extracting local coefficient vectors, and boundary
information, the routine is still rather easy to implement. The
implementation still does not depend on the actually used set of local
basis functions.

The function \code{update()} which is now described in detail, can be
seen as an example for the very flexible implementation of rather
complex nonlinear and time dependent problems which often show the
same structure (compare the implementation of the assembling routine
for the time dependent heat equation, Section~\ref{S:heat_build}). It
demonstrates the functionality and flexibility of the \ALBERTA tools:
the assemblage of complex problems is still quite easy, whereas the
resulting code is quite efficient.

Similar to the linear Poisson solver, we provide a function
\code{LALt()} for the second order term. Besides the additional
scaling by the heat conductivity $k$, it is exactly the same as for the
Poisson problem.  For the nonlinear reaction--diffusion equation we
also need a function \code{c()} for the zero order term. This term is
assembled using element-wise quadrature and thus needs information
about the function $v$ used in the linearization at all quadrature
nodes.  Information for \code{LALt()} and \code{c()} is stored in the
data structure \code{struct op\_info}, see above.  The members of this
structure are initialized during mesh traversal in \code{update()}.
\bv\begin{verbatim}
static const REAL_B *LALt(const EL_INFO *el_info,
			  const QUAD *quad, 
			  int iq, void *ud)
{
  struct op_data *info = (struct op_data *)ud;
  REAL           fac = info->k*info->det;
  int            i, j, k, dim = el_info->mesh->dim;
  static REAL_BB LALt;

  for (i = 0; i <= dim; i++) {
    for (j = i; j <= dim; j++) {
      for (LALt[i][j] = k = 0; k < DIM_OF_WORLD; k++)
	LALt[i][j] += info->Lambda[i][k]*info->Lambda[j][k];
      LALt[i][j] *= fac;
      LALt[j][i] = LALt[i][j];
    }
  }
  return (const REAL_B *)LALt;
}

static REAL c(const EL_INFO *el_info, const QUAD *quad, int iq, void *ud)
{
  struct op_data *info = (struct op_data *)ud;
  REAL  v3;

  DEBUG_TEST_EXIT(info->quad_fast->quad == quad, "quads differ\n");

  v3 = info->v_qp[iq]*info->v_qp[iq]*info->v_qp[iq];

  return(info->sigma*info->det*v3);
}
\end{verbatim}\ev

As mentioned above, we use a general Newton solver and a pointer to the 
\code{update()} routine is adjusted inside the function \code{nlsolve()} 
in the data structure for this solver. Such a solver does not have any
information about the actual problem, nor information about the
\ALBERTA data structures for storing DOF vectors and matrices. This is
also reflected in the arguments of \code{update()}:
\bv\begin{verbatim}
static void update(void *ud, int dim, const REAL *v, int up_DF, REAL *F);
\end{verbatim}\ev
Here, \code{dim} is the dimension of the discrete nonlinear problem,
\code{v} is a \emph{vector} storing the coefficients of the 
finite element function which is used for the linearization, \code{up\_DF}
is a flag indicating whether $DF(v)$ should be assembled or not.
If \code{F} is not \nil, then $F(v)$ should be assembled and
stored in the \emph{vector} \code{F}. Information about the
\ALBERTA finite element space, a pointer to a DOF matrix, etc. can
be passed to \code{update()} by the \code{ud} pointer. The declaration
\bv\begin{verbatim}
  NEWTON_DATA    *data = (NEWTON_DATA *)ud;
\end{verbatim}\ev
converts the \code{void *} pointer \code{ud} into a pointer \code{data} to
a structure \code{NEWTON\_DATA} which gives access to all 
information, used for the assembling (see above). This structure
is initialized in \code{nlsolve()} before starting the Newton method.

The \code{update()} routine contains three parts: an
initialization of the assembling functions (only done on the first
call), a conversion of the vectors that are arguments to the routine
into DOF vectors, and finally the assembling.

\paragraph{Initialization of the assembling functions.}

The initialization of \ALBERTA functions for the assembling is similar
to the initialization in the \code{build()} routine of the linear
Poisson equation (compare Section~\ref{S:ellipt_build}). There are
minor differences:

\begin{enumerate}
\item In addition to the assemblage of the 2nd order term
(see the function \code{LALt()}), we now have to assemble the zero
order term too (see the function \code{c()}).  The integration of the
zero order term has to be done by using an element wise quadrature
which needs the values of $v^3$ at all quadrature nodes.  The two
element matrices are computed separately. This makes it possible to
use them for the system matrix and right hand side.

\item In the solver for the Poisson problem, we have filled an
  \code{OPERATOR\_INFO} structure with information about the
  differential operator. This structure is an argument to
  \code{fill\_matrix\_info()} which returns a pointer to a structure
  \code{EL\_MATRIX\_INFO}. This pointer is used for the complete
  assemblage of the system matrix by some \ALBERTA routine. A detailed
  description of this structures and the general assemblage routines
  for matrices can be found in \secref{S:matrix_assemblage}.  Here, we
  want to use only the function for computing the element matrices.
  Thus, we only need the entries \code{el\_matrix\_fct()} and
  \code{fill\_info} of the \code{EL\_MATRIX\_INFO} structure, which
  are used to compute the element matrix (\code{fill\_info} is the
  second argument to \code{el\_matrix\_fct()}). We initialize a
  function pointer \code{fill\_a} with data pointer \code{a\_info} for
  the computation of the element matrix $\boldsymbol{A}_S$ and a
  function pointer \code{fill\_c} with data pointer \code{c\_info} for
  the computation $\boldsymbol{M}_S$.
  
  All other information inside the \code{EL\_MATRIX\_INFO} structure
  is used for the automatic assembling of element matrices into the
  system matrix by \code{update\_matrix()}. Such information can be
  ignored here, since this is now done in \code{update()}.
  
\item For the assembling of the element matrix into the system matrix
  and the element contribution of the right hand side into the load
  vector we need information about the number of local basis
  functions, \code{n\_phi}, and how to access global DOFs from the
  elements, \code{get\_dof()}. This function uses the DOF
  administration \code{admin} of the finite element space.  We also
  need information about the boundary type of the local basis
  functions, \code{get\_bound()}, and for the computation of the
  values of $v$ at quadrature nodes, we have to extract the local
  coefficient vector from the global one, \code{get\_v\_loc()}. These
  functions and the number of local basis functions can be accessed
  via the \code{bas\_fcts} inside the \code{data->fe\_space}
  structure. The used \code{admin} is the \code{admin} structure in
  \code{data->fe\_space}. For details about these functions we refer
  to Sections \ref{S:basfct_data}, \ref{S:dof_access}, and
  \ref{book:S:eval_fe}.
\end{enumerate}

\paragraph{Conversion of the vectors into DOF vectors.}

The input vector \code{v} of \code{update()} is a vector storing the
coefficients of the function used for the linearization. It is not a
DOF vector, but \ALBERTA routines for extracting a local coefficient
vector need a DOF vector. Thus, we have to ``convert'' \code{v} into
some DOF vector \code{dof\_v}. This is done by calls
%%
\bv\begin{lstlisting}
  init_dof_real_vec_skel(dof_v, "v", data->fe_space);
  distribute_to_dof_real_vec_skel(dof_v, v);
\end{lstlisting}\ev
%%
We refer the reader to \secref{S:dof_vec_skel} for a more detailed
discussion.

In the same way we have to convert \code{F} to a DOF vector \code{dof\_F}
if \code{F} is not \nil.

\paragraph{The assemblage of the linearized system.}

If the system matrix has to be assembled, then the
DOF matrix \code{data->DF} is cleared and we check which solver
can be used for solving the linearized equation. 

If the right hand side has to be assembled, then this vector 
is initialized with values
\[
- \int_\Omega (f + \sigma u_{ext}^4) \varphi_j  \,dx.
\]

For the assemblage of the element contributions we use the
non--recursive mesh traversal routines. On each element we access
the local coefficient vector \code{v\_loc}, the global DOFs \code{dof}
and boundary types \code{bound} of the local basis functions.

Next, we initialize the Jacobian of the barycentric coordinates and
compute the values of $v$ at the quadrature node by \code{uh\_at\_qp()}.
Hence $v^3$ can easily be calculated in \code{c()} at all quadrature
nodes. Routines for evaluating finite element functions and their
derivatives are described in detail in \secref{S:eval}.

Now, all members of \code{struct op\_info} are initialized, and
we compute the element matrices $\boldsymbol{A}_S$ by the function
\code{fill\_a()} and $\boldsymbol{M}_S$ by the function \code{fill\_c()}.

These contributions are added to the system matrix if \code{up\_DF} is
not zero.  Finally, the right hand side contributions for all non
Dirichlet DOFs are computed, and zero Dirichlet boundary values are
set for Dirichlet DOFs, if \code{F} is not \nil.

The following sources code listing quotes the entire \code{update()}
sub-routine:
%%
\bv\begin{verbatim}
static void update(void *ud, int dim, const REAL *v, bool up_DF, REAL *F)
{
  /* Some quantities remembered across calls. Think of this routine
   * like being a "library function" ... The stuff is re-initialized
   * whenever the finite element space changes. We use fe_space->admin
   * to check for changes in the finite element space because
   * DOF_ADMIN's are persisitent within ALBERTA, while fe-space are
   * not.
   */
  static EL_MATRIX_INFO elmi2, elmi0;
  static const DOF_ADMIN *admin = NULL;
  static struct op_data op_data[1]; /* storage for det and Lambda */

  /* Remaining (non-static) variables. */
  const BAS_FCTS  *bas_fcts = NULL;
  int             n_phi;
  int             mesh_dim;
  NEWTON_DATA    *data = (NEWTON_DATA *)ud;
  FLAGS          fill_flag;
  DOF_REAL_VEC   dof_v[1];
  DOF_REAL_VEC   dof_F[1];

/*--------------------------------------------------------------------------*/
/*  init functions for assembling DF(v) and F(v)                            */
/*--------------------------------------------------------------------------*/

  bas_fcts = data->fe_space->bas_fcts;
  n_phi    = bas_fcts->n_bas_fcts;
  mesh_dim = bas_fcts->dim;

  if (admin != data->fe_space->admin) {
    OPERATOR_INFO o_info2 = { NULL, }, o_info0 = { NULL, };
    const QUAD    *quad;

    admin     = data->fe_space->admin;

    quad     = get_quadrature(mesh_dim, 2*bas_fcts->degree-2);
    o_info2.row_fe_space   = data->fe_space;
    o_info2.quad[2]        = quad;
    o_info2.LALt.real      = LALt;
    o_info2.LALt_pw_const  = true;
    o_info2.LALt_symmetric = true;
    o_info2.user_data      = op_data;

    fill_matrix_info(&o_info2, &elmi2);

    o_info0.row_fe_space   = data->fe_space;
    o_info0.quad[0]        = quad;
    o_info0.c.real         = c;
    o_info0.c_pw_const     = false;
    o_info0.user_data      = op_data;

    fill_matrix_info(&o_info0, &elmi0);

    op_data->quad_fast = get_quad_fast(bas_fcts, quad, INIT_PHI);
  }

/*--------------------------------------------------------------------------*/
/*  make a DOF vector from input vector v_vec                               */
/*--------------------------------------------------------------------------*/
  init_dof_real_vec_skel(dof_v, "v", data->fe_space);
  distribute_to_dof_real_vec_skel(dof_v, v);

/*--------------------------------------------------------------------------*/
/*  make a DOF vector from F, if not NULL                                    */
/*--------------------------------------------------------------------------*/
  if (F) {
    init_dof_real_vec_skel(dof_F, "F(v)", data->fe_space);
    distribute_to_dof_real_vec_skel(dof_F, F);
  }

/*--------------------------------------------------------------------------*/
/*  and now assemble DF(v) and/or F(v)                                      */
/*--------------------------------------------------------------------------*/
  op_data->k     = data->k;
  op_data->sigma = data->sigma;

  if (up_DF)
  {
/*--- if v_vec[i] >= 0 for all i => matrix is positive definite (p=1)  ----*/
    data->solver = dof_min(dof_v) >= 0 ? CG : BiCGStab;
    clear_dof_matrix(data->DF);
  }

  if (F)
  {
    dof_set(0.0, dof_F);			//!! Seggi
    L2scp_fct_bas(data->f, op_data->quad_fast->quad, dof_F);
    dof_scal(-1.0, dof_F);
  }

  fill_flag = CALL_LEAF_EL|FILL_COORDS|FILL_BOUND;
  TRAVERSE_FIRST(data->fe_space->mesh, -1, fill_flag) {
    const EL_REAL_VEC  *v_loc;
    const EL_DOF_VEC   *dof;
    const EL_BNDRY_VEC *bndry_bits;
    EL_SCHAR_VEC bound[n_phi];
    const EL_MATRIX *elmat2, *elmat0;

    v_loc      = fill_el_real_vec(NULL, el_info->el, dof_v);
    dof        = get_dof_indices(NULL, data->fe_space, el_info->el);
    bndry_bits = get_bound(NULL, bas_fcts, el_info);

/*--------------------------------------------------------------------------*/
/* initialization of values used by LALt and c                              */
/*--------------------------------------------------------------------------*/

    op_data->det = el_grd_lambda_0cd(el_info, op_data->Lambda);

    op_data->v_qp = uh_at_qp(NULL, op_data->quad_fast, v_loc);

    elmat2 = elmi2.el_matrix_fct(el_info, elmi2.fill_info);

    elmat0 = elmi0.el_matrix_fct(el_info, elmi0.fill_info);

    /* Translate the geometric boundary classification into
     * Dirichlet/Neumann/Interior boundary condition
     * interpretation. Inside the loop over the mesh-elements we need
     * only to care about Dirichlet boundary conditions.
     */
    dirichlet_map(bound, bndry_bits, data->dirichlet_mask);

    if (up_DF)  /*---  add element contribution to matrix DF(v)  ----------*/
    {
/*--------------------------------------------------------------------------*/
/*  add a(phi_i,phi_j) + 4*m(u^3*phi_i,phi_j) to matrix                     */
/*--------------------------------------------------------------------------*/
      add_element_matrix(data->DF, 1.0, elmat2, NoTranspose, dof, dof, bound);
      add_element_matrix(data->DF, 4.0, elmat0, NoTranspose, dof, dof, bound);
    }

    if (F)     /*---  add element contribution to F(v)  --------------------*/
    {
      int i;
/*--------------------------------------------------------------------------*/
/*  F(v) += a(v, phi_i) + m(v^4, phi_i)                                     */
/*--------------------------------------------------------------------------*/
      bi_mat_el_vec(1.0, elmat2, 1.0, elmat0, v_loc, 1.0, dof_F, dof, bound);
      for (i = 0; i < n_phi; i++) {
	if (bound->vec[i] >= DIRICHLET) {
	  F[dof->vec[i]] = 0.0;  /*--- zero Dirichlet boundary conditions! -*/
	}
      }
    }
    
  } TRAVERSE_NEXT();

  /* Record that the boundary conditions are built into the matrix, needed
   * e.g. by the hierarchical preconditioners.
   */
  BNDRY_FLAGS_CPY(data->DF->dirichlet_bndry, data->dirichlet_mask);
\end{verbatim}\ev

\subsubsection{The linear sub--solver}

For the solution of the linearized problem we use the \code{oem\_solve\_s()}
function, which is also used in the solver for the linear Poisson equation
(compare Section~\ref{S:ellipt_solve}). Similar to the \code{update()}
function, we have to convert the right hand side vector \code{F}
and the solution vector \code{d} to DOF vectors. Information about
the system matrix and parameters for the solver are passed by \code{ud}.
The member \code{data->solver} is initialized in \code{update()}.

\bv\begin{verbatim}
static int solve(void *ud, int dim, const REAL *F, REAL *d)
{
  NEWTON_DATA   *data = (NEWTON_DATA *)ud;
  int           iter;
  DOF_REAL_VEC  dof_F[1];
  DOF_REAL_VEC  dof_d[1];
  
/*--------------------------------------------------------------------------*/
/*  make DOF vectors from F and d                                           */
/*--------------------------------------------------------------------------*/
  init_dof_real_vec_skel(dof_F, "F", data->fe_space);
  distribute_to_dof_real_vec_skel(dof_F, F);
  
  init_dof_real_vec_skel(dof_d, "d", data->fe_space);
  distribute_to_dof_real_vec_skel(dof_d, d);

  
  if (data->icon == ILUkPrecon)
	  data->precon = init_oem_precon(data->DF, NULL, data->info, ILUkPrecon, data->ilu_k);
  else
	  data->precon = init_oem_precon(data->DF, NULL, data->info, data->icon,
                                 data->ssor_omega, data->ssor_iter);
  iter = oem_solve_s(data->DF, NULL, dof_F, dof_d, data->solver, 
		     data->tolerance, data->precon, data->restart, 
		     data->max_iter, data->info);

  return iter;
}
\end{verbatim}\ev

\subsubsection{The computation of the $H^1$ semi norm}

The $H^1$ semi norm can easily be calculated by converting the input
vector \code{v} into a DOF-vector and then calling the \ALBERTA routine
\code{H1\_norm\_uh()} (compare Section~\ref{S:eval_norm}).
\bv\begin{verbatim}
static REAL norm(void *ud, int dim, const REAL *v)
{
  NEWTON_DATA   *data = (NEWTON_DATA *)ud;
  DOF_REAL_VEC  dof_v[1]; /* = {NULL, NULL, "v"};*/
  
  init_dof_real_vec_skel(dof_v, "v", data->fe_space);
  distribute_to_dof_real_vec_skel(dof_v, v);

  return H1_norm_uh(NULL, dof_v);
}
\end{verbatim}\ev

\subsubsection{The nonlinear solver}

The function \code{nlsolve()} initializes the structure \code{NEWTON\_DATA}
with problem dependent information. Here, we have to allocate a
DOF matrix for storing the system matrix (only on the first call),
and initialize parameters for the linear sub--solver and problem dependent
data (like heat conductivity $k$, etc.)

The structure \code{NLS\_DATA} is filled with information for the
general Newton solver (the problem dependent routines \code{update()},
\code{solve()}, and \code{norm()} described above). All these routines use
the same structure \code{NEWTON\_DATA} for problem dependent
information.

The dimension of the discrete equation is 
\bv\begin{verbatim}
dim = u0->fe_space->admin->size_used;
\end{verbatim}\ev
where \code{u0} is a pointer to a DOF vector storing the initial
guess.  Note, that after the call to \code{dof\_compress()} in the
\code{build()} routine, \code{dim} holds the true dimension of the
discrete equation.  Without a \code{dof\_compress()} there may be
holes in DOF vectors, and \verb|u0->fe_space->admin->size_used|
bigger than the last \emph{used} index, and again \code{dim} is the
dimension of the discrete equation for the Newton solver. The
\ALBERTA routines do not operate on unused indices, whereas the
Newton solvers do operate on unused indices too, because they do not
know about used and unused indices. In this situation, all unused
DOFs would have to be cleared for the initial solution \code{u0} by
\bv\begin{verbatim} 
FOR_ALL_FREE_DOFS(u0->fe_space->admin, u0->vec[dof] = 0.0);
\end{verbatim}\ev
The same applies to the vector storing the right hand side in 
\code{update()}. The \code{dof\_set()} function only initializes
used indices.

Finally, we reallocate the workspace used by the Newton solvers
(compare Section~\ref{S:nls}) and start the Newton method.
\bv\begin{verbatim}
int nlsolve(DOF_REAL_VEC *u0, REAL k, REAL sigma, REAL (*f)(const REAL_D),
	    const BNDRY_FLAGS dirichlet_mask)
{
  FUNCNAME("nlsolve");
  static NEWTON_DATA  data =
    { NULL, { 0, }, 0, 0, NULL, NULL, CG, 1.e-8, 1.0, 1000, 1, 8, 0, 2, 0, NULL };
  static NLS_DATA     nls_data;
  int                 iter, dim = u0->fe_space->admin->size_used;

  if (!data.fe_space)
  {
/*--------------------------------------------------------------------------*/
/*--  init parameters for newton  ------------------------------------------*/
/*--------------------------------------------------------------------------*/
    nls_data.update      = update; 
    nls_data.update_data = &data;
    nls_data.solve       = solve;
    nls_data.solve_data  = &data;
    nls_data.norm        = norm;
    nls_data.norm_data   = &data;

    nls_data.tolerance = 1.e-4;
    GET_PARAMETER(1, "newton tolerance", "%e", &nls_data.tolerance);
    nls_data.max_iter = 50;
    GET_PARAMETER(1, "newton max. iter", "%d", &nls_data.max_iter);
    nls_data.info = 8;
    GET_PARAMETER(1, "newton info", "%d", &nls_data.info);
    nls_data.restart = 0;
    GET_PARAMETER(1, "newton restart", "%d", &nls_data.restart);
    

/*--------------------------------------------------------------------------*/
/*--  init data for update and solve  --------------------------------------*/
/*--------------------------------------------------------------------------*/

    data.fe_space = u0->fe_space;
    data.DF       = get_dof_matrix("DF(v)", u0->fe_space, NULL);

    data.tolerance = 1.e-2*nls_data.tolerance;
    GET_PARAMETER(1, "linear solver tolerance", "%f", &data.tolerance);
    GET_PARAMETER(1, "linear solver max iteration", "%d", &data.max_iter);
    GET_PARAMETER(1, "linear solver info", "%d", &data.info);
    GET_PARAMETER(1, "linear solver precon", "%d", &data.icon);
    if (data.icon == __SSORPrecon) {
    	GET_PARAMETER(1, "linear precon ssor omega", "%f", &data.ssor_omega);
    	GET_PARAMETER(1, "linear precon ssor iter", "%d", &data.ssor_iter);
    }
    if (data.icon == ILUkPrecon)
        GET_PARAMETER(1, "linear precon ilu(k)", "%d", &data.ilu_k);
    GET_PARAMETER(1, "linear solver restart", "%d", &data.restart);
  }
  TEST_EXIT(data.fe_space == u0->fe_space, "can't change f.e. spaces\n");
  BNDRY_FLAGS_CPY(data.dirichlet_mask, dirichlet_mask);

/*--------------------------------------------------------------------------*/
/*--  init problem dependent parameters  -----------------------------------*/
/*--------------------------------------------------------------------------*/
  data.k     = k;
  data.sigma = sigma;
  data.f     = f;

/*--------------------------------------------------------------------------*/
/*--  enlarge workspace used by newton(_fs), and solve by Newton  ----------*/
/*--------------------------------------------------------------------------*/
  if (nls_data.restart)
  {
    nls_data.ws = REALLOC_WORKSPACE(nls_data.ws, 4*dim*sizeof(REAL));
    iter = nls_newton_fs(&nls_data, dim, u0->vec);
  }
  else
  {
    nls_data.ws = REALLOC_WORKSPACE(nls_data.ws, 2*dim*sizeof(REAL));
    iter = nls_newton(&nls_data, dim, u0->vec);
  }

  return iter;
}
\end{verbatim}\ev

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