\section{Nonlinear problem}
\label{ss:nonlinear problem}
We define the nonlinear problem
\begin{eqnarray}
-\Delta u + u^4 &=& f~~~\mbox{in } \Omega \subset \mathbb{R}^{dim}\\
u &=& g~~~ \mbox{on } \partial\Omega.
\end{eqnarray}
We choose the functions $f$ and $g$ so that the exact solution again is $u(x)=e^{-10x^2}$. This leads to 
\begin{eqnarray}
f(x)&=& -\left(400-20dow\right)e^{-10x^2} + \left(e^{-10x^2}\right)^4\\
g(x)&=&e^{-10x^2},
\end{eqnarray}
with $dow$ the world dimension.

We linearize the problem using the Newton method. First, we define an
initial guess $u_0$ of the solution which is $0$ for the first
adaptation loop iteration. In later iterations we can use the solution
of the last iteration interpolated to the current mesh as initial
guess. In each Newton step, a correction $d$ for the solution of the
last step is computed by solving
\begin{eqnarray}
\label{eq:newton step}
DF(u_n)(d) = F(u_n)
\end{eqnarray}
for $d$, where $F(u):=-\Delta u + u^4 - f$ and
\begin{eqnarray}
DF(u_n)(d) &=& \lim_{h \rightarrow 0} \frac{F(u_n+hd)-F(u_n)}{h}\\
&=& \lim_{h \rightarrow 0} \frac{-\Delta u_n - h \Delta d + \Delta u_n}{h} + \lim_{h \rightarrow 0}\frac{(u_n+hd)^4-u_n^4}{h}\\
&=&  -\Delta d + 4 u_n^3 d
\end{eqnarray}
the directional derivative of $F$ at $u_n$ along $d$. 

Then the solution is updated:
\begin{equation}
u_{n+1} := u_n - d.
\end{equation}
We repeat this procedure until $||d||_{L^2} < tol$ with $tol$ a given tolerance for the Newton method. 

In our example, equation (\ref{eq:newton step}) reads:
\begin{equation}
-\Delta d + 4 u_n^3 d = - \Delta u_n + u_n^4 - f.
\end{equation}

\begin{figure}
\center
\includegraphics[width=0.7\textwidth]{fig/nonlin_iteration.pdf}
\caption{Solve step of the nonlinear problem.}
\label{f:nonlin iteration}
\end{figure}
In Figure \ref{f:nonlin iteration}, the Newton method is illustrated.

\begin{table}
\center
\begin{tabular}{|cc|c|c|c|}
\hline 
& & {\bf 1d} & {\bf 2d} & {\bf 3d} \\
\hline 
{\bf source code} & \tt src/ & \multicolumn{3}{|c|}{\tt nonlin.cc} \\
\hline 
{\bf parameter file} & \tt init/ & \tt nonlin.dat.1d & \tt nonlin.dat.2d & \tt nonlin.dat.3d \\
\hline 
{\bf macro file} & \tt macro/ & \tt macro\_big.stand.1d & \tt macro\_big.stand.2d & \tt macro\_big.stand.3d \\
\hline 
{\bf output files} & \tt output/ & \multicolumn{3}{|c|}{\tt nonlin.mesh, nonlin.dat} \\
\hline 
\end{tabular}
\caption{Files of the {\tt nonlin} example.}
\end{table}

\subsection{Source code}
The main idea is to realize the Newton method as implementation of the {\it ProblemIterationInterface}, which replaces the standard iteration in the adaptation loop. The Newton method has to know the problem which has to be solved, as well as the implementation of one Newton step. Estimation and adaptation is done by the problem object. The assemble step and solve step are delegated to a Newton-step object. 

Now, we describe the code step by step. The function $g$ is defined like in the previous examples. In the following, we define a zero-function which is later used to implement the Dirichlet boundary condition for the Newton-step implementation (at domain boundaries no correction has to be done). The function \verb+F+ implements the right hand side function $f$.

\begin{lstlisting}{}
class Zero : public AbstractFunction<double, WorldVector<double> >
{
public:
  double operator()(const WorldVector<double>& x) const 
  {
    return 0.0;
  }
};

class F : public AbstractFunction<double, WorldVector<double> >
{
public:
  /// Constructor
  F(int degree)
    : AbstractFunction<double, WorldVector<double> >(degree)
  {}

  /// Implementation of AbstractFunction::operator().
  double operator()(const WorldVector<double>& x) const 
  {
    int dow = x.getSize();
    double r2 = x * x;
    double ux = exp(-10.0 * r2);
    double ux4 = ux * ux * ux * ux;
    return ux4 -(400.0 * r2 - 20.0 * dow) * ux;
  }
};
\end{lstlisting}

The class \verb+X3+ implements the function $u^3(x)$ used within the Newton step.

\begin{lstlisting}{}
class X3 : public AbstractFunction<double, double>
{
public:
  X3() : AbstractFunction<double, double>(3) {}

  /// Implementation of AbstractFunction::operator().
  double operator()(const double& x) const 
  {
    return x * x * x;
  }
};
\end{lstlisting}

In the following, we define an interface which has to be implemented by the Newton-step object. 

\begin{lstlisting}{}
class NewtonStepInterface
{
public:
  virtual void initNewtonStep(AdaptInfo *adaptInfo) = 0;
  virtual void exitNewtonStep(AdaptInfo *adaptInfo) = 0;
  virtual void assembleNewtonStep(AdaptInfo *adaptInfo, Flag flag) = 0;
  virtual void solveNewtonStep(AdaptInfo *adaptInfo) = 0;
  virtual DOFVector<double> *getCorrection() = 0;
};
\end{lstlisting}

The \verb+initNewtonStep+ method is called before each Newton step, the method \verb+exitNewtonStep+ after each Newton step. \verb+assembleNewtonStep+ assembles the linear system of equations needed for the next step, \verb+solveNewtonStep+ solves this system of equations. The solution is the correction $d$. The method \verb+getCorrection+ returns a pointer to the vector storing the correction.

Now, the Newton method will be implemented. Actually, the class \verb+NewtonMethod+ replaces the whole iteration in the adaptation loop, including mesh adaptation and error estimation. The Newton method, which is a loop over Newton steps, is one part of this iteration.

\begin{lstlisting}{}
class NewtonMethod : public ProblemIterationInterface
{
public:
  NewtonMethod(const char *name,
               ProblemScal *problem,
               NewtonStepInterface *step)
    : problemNonlin(problem),
      newtonStep(step),
      newtonTolerance(1e-8),
      newtonMaxIter(100)
  {
    GET_PARAMETER(0, std::string(name) + "->tolerance", "%f", 
                  &newtonTolerance);
    GET_PARAMETER(0, std::string(name) + "->max iteration", "%d", 
                  &newtonMaxIter);
    solution = problemNonlin->getSolution();
    correction = newtonStep->getCorrection();
  }
\end{lstlisting}

In the constructor, pointers to the nonlinear problem and to the Newton-step object are stored to the class members \verb+problemNonlin+ and \verb+newtonStep+. Furthermore, the parameters \verb+newtonTolerance+ and \verb+newtonMaxIter+ are initialized, and pointers to the nonlinear solution and to the correction vector are stored. 

The following methods define one iteration in the adaptation loop.

\begin{lstlisting}{}
  void beginIteration(AdaptInfo *adaptInfo) 
  {
    FUNCNAME("NewtonMethod::beginIteration()");
    MSG("\n");
    MSG("begin of iteration %d\n", adaptInfo->getSpaceIteration()+1);
    MSG("=============================\n");
  }

  Flag oneIteration(AdaptInfo *adaptInfo, Flag toDo = FULL_ITERATION)
  {
    Flag flag = 0, markFlag = 0;

    if (toDo.isSet(MARK)) markFlag = problemNonlin->markElements(adaptInfo);
    if (toDo.isSet(ADAPT) && markFlag.isSet(MESH_REFINED))
      flag = problemNonlin->refineMesh(adaptInfo);
    if (toDo.isSet(ADAPT) && markFlag.isSet(MESH_COARSENED))
      flag |= problemNonlin->coarsenMesh(adaptInfo);

    if (toDo.isSet(SOLVE)) {
      newtonStep->initNewtonStep(adaptInfo);
      int newtonIteration = 0;
      double res = 0.0;
      do {    
	newtonIteration++;
	newtonStep->assembleNewtonStep(adaptInfo, flag);    
	newtonStep->solveNewtonStep(adaptInfo);
	res = correction->L2Norm();
	*solution -= *correction;
	MSG("newton iteration %d: residual %f (tol: %f)\n",
	    newtonIteration, res, newtonTolerance);
      } while ((res > newtonTolerance) && (newtonIteration < newtonMaxIter));    

      newtonStep->exitNewtonStep(adaptInfo);
    }

    if (toDo.isSet(ESTIMATE)) problemNonlin->estimate(adaptInfo);
    return flag;
  }
  
  void endIteration(AdaptInfo *adaptInfo) 
  {
    FUNCNAME("NewtonMethod::endIteration()");
    MSG("\n");
    MSG("end of iteration number: %d\n", adaptInfo->getSpaceIteration()+1);
    MSG("=============================\n");
  }
\end{lstlisting}

The methods \verb+beginIteration+ and \verb+endIteration+ only print some information to the standard output. In \verb+oneIteration+, the iteration, including the loop over the Newton steps, is defined.

Finally, the methods \verb+getNumProblems+ and \verb+getProblem+ are implemented to complete the \\\verb+ProblemIterationInterface+, and the private class members are defined.

\begin{lstlisting}{}
  int getNumProblems() { return 1; };

  ProblemStatBase *getProblem(int number = 0) 
  {
    FUNCNAME("NewtonMethod::getProblem()");
    if (number == 0) return problemNonlin;
    ERROR_EXIT("invalid problem number\n");
    return NULL;
  }

private:
  ProblemScal *problemNonlin;
  NewtonStepInterface *newtonStep;
  double newtonTolerance;
  int newtonMaxIter;
  DOFVector<double> *solution;
  DOFVector<double> *correction;
};
\end{lstlisting}

The class \verb+Nonlin+ implements both, the nonlinear problem and the Newton-step. Since the Newton step is accessed over an own interface, it is always clear, in which role a \verb+Nonlin+ instance is called by the Newton method. 

\begin{lstlisting}{}
class Nonlin : public ProblemScal,
	       public NewtonStepInterface
{
public:
  Nonlin(const char *name)
    : ProblemScal(name)
  {}
\end{lstlisting}

In the constructor, the base class constructor of \verb+ProblemScal+ is called and the name is given to it.

In the initialization, the base class intialization is called, the correction vector is created and initialized, and Dirichlet boundary conditions are created.

\begin{lstlisting}{}
  void initialize(Flag initFlag,
		  ProblemScal *adoptProblem = NULL,
		  Flag adoptFlag = INIT_NOTHING)
  {
    ProblemScal::initialize(initFlag, adoptProblem, adoptFlag);
    correction = new DOFVector<double>(this->getFESpace(), "old solution");
    correction->set(0.0);

    dirichletZero = new DirichletBC(1, &zero, feSpace_);
    dirichletG = new DirichletBC(1, &g, feSpace_);

    solution->getBoundaryManager()->addBoundaryCondition(dirichletG);
    systemMatrix->getBoundaryManager()->
      addBoundaryCondition(dirichletZero);
    rhs->getBoundaryManager()->addBoundaryCondition(dirichletZero);
    correction->getBoundaryManager()->addBoundaryCondition(dirichletZero);
  };
\end{lstlisting}

To the solution of the nonlinear problem the function $g$ is applied as Dirichlet boundary function. The system matrix, the correction vector and the right hand side vector build the system for the Newton step. Since no correction has to be done at the domain boundaries, zero Dirichlet conditions are applied to them.  

In the destructor, the allocated memory is freed.

\begin{lstlisting}{}
  ~Nonlin()
  { 
    delete correction; 
    delete dirichletZero;
    delete dirichletG;
  }
\end{lstlisting}

Now, we implement the Newton step functionality. First, in \verb+initNewtonStep+, we fill the solution vector with boundary values. This will not be done automatically because we let the \verb+solution_+ pointer point to \verb+correction+. The address of \verb+solution_+ is stored in \verb+tmp+. After the Newton method is finished, the \verb+solution_+ pointer is reset to its original value in \verb+exitNewtonStep+. 

\begin{lstlisting}{}
  void initNewtonStep(AdaptInfo *adaptInfo) 
  {
    solution->getBoundaryManager()->initVector(solution);
    TraverseStack stack;
    ElInfo *elInfo = stack.traverseFirst(mesh_, -1, 
					 Mesh::CALL_LEAF_EL | 
					 Mesh::FILL_COORDS |
					 Mesh::FILL_BOUND);
    while (elInfo) {
      solution->getBoundaryManager()->fillBoundaryConditions(elInfo, solution);
      elInfo = stack.traverseNext(elInfo);
    }
    solution->getBoundaryManager()->exitVector(solution);

    tmp = solution;
    solution = correction;
  }

  void exitNewtonStep(AdaptInfo *adaptInfo) 
  {
    solution = tmp;
  }
\end{lstlisting}

The implementation of \verb+assembleNewtonStep+ and \verb+solveNewtonStep+ just delegates the calls to the base class implementations in \verb+ProblemScal+.

\begin{lstlisting}{}
  void assembleNewtonStep(AdaptInfo *adaptInfo, Flag flag) 
  {
    ProblemScal::buildAfterCoarsen(adaptInfo, flag);
  }

  void solveNewtonStep(AdaptInfo *adaptInfo) 
  {
    ProblemScal::solve(adaptInfo);
  }
\end{lstlisting}

Finally, the \verb+getCorrection+ method is implemented and private class members are defined.

\begin{lstlisting}{}
  DOFVector<double> *getCorrection() 
  { 
    return correction; 
  }

private:
  DOFVector<double> *correction;
  DOFVector<double> *tmp;
  Zero zero;
  G g;
  DirichletBC *dirichletZero;
  DirichletBC *dirichletG;
};
\end{lstlisting}

Now, we start with the main program.

\begin{lstlisting}{}
int main(int argc, char* argv[])
{
  FUNCNAME("main");

  TEST_EXIT(argc == 2)("usage: nonlin initfile\n");

  Parameters::init(false, argv[1]);

  Nonlin nonlin("nonlin");
  nonlin.initialize(INIT_ALL);

  AdaptInfo *adaptInfo = new AdaptInfo("nonlin->adapt", 1);

  NewtonMethod newtonMethod("nonlin->newton", &nonlin, &nonlin);

  AdaptStationary *adapt = new AdaptStationary("nonlin->adapt",
					       &newtonMethod,
					       adaptInfo);
\end{lstlisting}

An instance of class \verb+NewtonMethod+ was created with a pointer to
a \verb+Nonlin+ object as nonlinear problem and as Newton-step
implementation. Instead of the nonlinear problem, now, the object
\verb+newtonMethod+ is handed to the adaptation loop as implementation
of \verb+ProblemIterationInterface+.

We have to add operators representing the Newton step equation
$-\Delta d + 4u_n^3d=-\Delta u_n + u_n^4 -f$ as well as operators
representing the nonlinear problem $-\Delta u_n + u_n^4 = f$. When the
operators are given to the problem, one can determine an assemble
factor (second argument of \verb+addMatrixOperator+ and
\verb+addVectorOperator+) as well as an estimation factor (third
argument) for each operator. So, we can manage both equations in one
problem instance. Note that in the Newton step we solve for $d$. $u_n$
there is known from the last Newton step. The term $u_n^4$ was
implemented as $u_n^3 v$, where $v$ for the Newton step equation is
equal to $d$ and in the nonlinear problem equation equal to $u_n$. So,
the corresponding operator can be used in both equations just with
different factors. In Figure \ref{f:nonlin operators}, the operator
factors for the assemble and the estimate step are shown.
\begin{figure}
\center
\includegraphics[width=0.5\textwidth]{fig/nonlin_operators.pdf}
\caption{Operator factors for the assemble step and for the estimate step.}
\label{f:nonlin operators}
\end{figure}

\begin{lstlisting}{}
  // ===== create operators =====
  double four = 4.0;
  double one = 1.0;
  double zero = 0.0;
  double minusOne = -1.0;

  Operator *nonlinOperator0 = new Operator(nonlin.getFESpace());
  nonlinOperator0->setUhOld(nonlin.getSolution());
  nonlinOperator0->addZeroOrderTerm(new VecAtQP_ZOT(nonlin.getSolution(), 
						    new X3));

  nonlin.addMatrixOperator(nonlinOperator0, &four, &one);
  nonlin.addVectorOperator(nonlinOperator0, &one, &zero);

  Operator *nonlinOperator2 = new Operator(nonlin.getFESpace());
  nonlinOperator2->setUhOld(nonlin.getSolution());
  nonlinOperator2->addSecondOrderTerm(new Laplace_SOT);

  nonlin.addMatrixOperator(nonlinOperator2, &one, &one);
  nonlin.addVectorOperator(nonlinOperator2, &one, &zero);

  int degree = nonlin.getFESpace()->getBasisFcts()->getDegree();

  Operator* rhsFunctionOperator = new Operator(nonlin.getFESpace());
  rhsFunctionOperator->addZeroOrderTerm(new CoordsAtQP_ZOT(new F(degree)));
  
  nonlin.addVectorOperator(rhsFunctionOperator, &minusOne, &one);
\end{lstlisting}

Finally, the adaptation loop is started and after it is fnished the result is written.

\begin{lstlisting}{}
  adapt->adapt();
  nonlin.writeFiles(adaptInfo, true);
}
\end{lstlisting}

\subsection{Parameter file}
The used parameter file \verb+nonlin.dat.2d+ looks like:

\begin{lstlisting}{}
dimension of world:   2

nonlinMesh->macro file name:      ./macro/macro_big.stand.2d
nonlinMesh->global refinements:   0

nonlin->mesh:  nonlinMesh

nonlin->dim:                  2
nonlin->polynomial degree:    1

nonlin->newton->tolerance:     1e-8
nonlin->newton->max iteration: 100

nonlin->solver:                 cg
nonlin->solver->max iteration:  1000
nonlin->solver->tolerance:      1.e-8
nonlin->solver->left precon:    diag

nonlin->estimator:              residual
nonlin->estimator->C0:          0.1
nonlin->estimator->C1:          0.1

nonlin->marker->strategy:       2
nonlin->marker->MSGamma:        0.5

nonlin->adapt->tolerance:       1e-1
nonlin->adapt->max iteration:   100

nonlin->output->filename:       output/nonlin
nonlin->output->AMDiS format:   1
nonlin->output->AMDiS mesh ext: .mesh
nonlin->output->AMDiS data ext: .dat
\end{lstlisting}

Here, as macro file \verb+macro_big.stand.2d+ is used, which is described in Section \ref{s:nonlin macro}. The parameters \verb+nonlin->newton->tolerance+ and \verb+nonlin->newton->max iteration+ determine the tolerance of the Newton solver and the maximal number of Newton steps within each iteration of the adaptation loop. 

\subsection{Macro file}
\label{s:nonlin macro}
The used macro file \verb+macro_big.stand.2d+ is very similar to the macro file of the first example described in Section \ref{s:ellipt macro}. Only the domain was changed from $\Omega=[0,1]^2$ to $\Omega=[-1, 1]^2$ by adapting the vertex coordinates correspondingly.

\subsection{Output}

\begin{figure}
\center
\begin{tabular}{cc}
\includegraphics[width=0.4\textwidth]{fig/nonlin_solution.jpg}&
\includegraphics[width=0.4\textwidth]{fig/nonlin_mesh.jpg}\\
solution & mesh
\end{tabular}
\caption{Solution and final mesh of the nonlinear problem}
\label{f:nonlin mesh and solution}
\end{figure}

In Figure \ref{f:nonlin mesh and solution}, the solution and the final mesh written after the adaptation loop are visualized. The solution is shown as height field where the values are interpreted as z-coordinates.