The step-14 tutorial program

Table of contents
  1. Introduction
  2. The commented program
  1. Results
  2. The plain program

Introduction

The maths

The Heidelberg group of Professor Rolf Rannacher, to which the three main authors of the deal.II library belonged during their PhD time and partly also afterwards, has been involved with adaptivity and error estimation for finite element discretizations since the mid-90ies. The main achievement is the development of error estimates for arbitrary functionals of the solution, and of optimal mesh refinement for its computation.

We will not discuss the derivation of these concepts in too great detail, but will implement the main ideas in the present example program. For a thorough introduction into the general idea, we refer to the seminal work of Becker and Rannacher [BR95],[BR96r], and the overview article of the same authors in Acta Numerica [BR01]; the first introduces the concept of error estimation and adaptivity for general functional output for the Laplace equation, while the second gives many examples of applications of these concepts to a large number of other, more complicated equations. For applications to individual types of equations, see also the publications by Becker [Bec95], [Bec98], Kanschat [Kan96], [FK97], Suttmeier [Sut96], [RS97], [RS98c], [RS99], Bangerth [BR99b], [Ban00w], [BR01a], [Ban02], and Hartmann [Har02], [HH01], [HH01b]. All of these works, from the original introduction by Becker and Rannacher to individual contributions to particular equations, have later been summarized in a book by Bangerth and Rannacher that covers all of these topics, see [BR03].

The basic idea is the following: in applications, one is not usually interested in the solution per se, but rather in certain aspects of it. For example, in simulations of flow problems, one may want to know the lift or drag of a body immersed in the fluid; it is this quantity that we want to know to best accuracy, and whether the rest of the solution of the describing equations is well resolved is not of primary interest. Likewise, in elasticity one might want to know about values of the stress at certain points to guess whether maximal load values of joints are safe, for example. Or, in radiative transfer problems, mean flux intensities are of interest.

In all the cases just listed, it is the evaluation of a functional $J(u)$ of the solution which we are interested in, rather than the values of $u$ everywhere. Since the exact solution $u$ is not available, but only its numerical approximation $u_h$, it is sensible to ask whether the computed value $J(u_h)$ is within certain limits of the exact value $J(u)$, i.e. we want to bound the error with respect to this functional, $J(u)-J(u_h)$.

For simplicity of exposition, we henceforth assume that both the quantity of interest $J$, as well as the equation are linear, and we will in particular show the derivation for the Laplace equation with homogeneous Dirichlet boundary conditions, although the concept is much more general. For this general case, we refer to the references listed above. The goal is to obtain bounds on the error, $J(e)=J(u)-J(u_h)$. For this, let us denote by $z$ the solution of a dual problem, defined as follows:

\[ a(\varphi,z) = J(\varphi) \qquad \forall \varphi, \]

where $a(\cdot,\cdot)$ is the bilinear form associated with the differential equation, and the test functions are chosen from the corresponding solution space. Then, taking as special test function $\varphi=e$ the error, we have that

\[ J(e) = a(e,z) \]

and we can, by Galerkin orthogonality, rewrite this as

\[ J(e) = a(e,z-\varphi_h) \]

for all possible functions $\varphi_h$ from the discrete test space.

Concretely, for Laplace's equation, the error identity reads

\[ J(e) = (\nabla e, \nabla(z-\varphi_h)). \]

For reasons that we will not explain, we do not want to use this formula as is, but rather split the scalar products into terms on all cells, and integrate by parts on each of them:

\begin{eqnarray*} J(e) &=& \sum_K (\nabla (u-u_h), \nabla (z-\varphi_h))_K \\ &=& \sum_K (-\Delta (u-u_h), z-\varphi_h)_K + (\partial_n (u-u_h), z-z_h)_{\partial K}. \end{eqnarray*}

Next we use that $-\Delta u=f$, and that $\partial_n u$ is a quantity that is continuous almost everywhere, so the terms involving $\partial_n u$ on one cell cancels with that on its neighbor, where the normal vector has the opposite sign. At the boundary of the domain, where there is no neighbor cell with which this term could cancel, the weight $z-\varphi_h$ can be chosen as zero, since $z$ has zero boundary values, and $\varphi_h$ can be chosen to have the same.

Thus, we have

\begin{eqnarray*} J(e) &=& \sum_K (f+\Delta u_h, z-\varphi_h)_K - (\partial_n u_h, z-\varphi_h)_{\partial K\backslash \partial\Omega}. \end{eqnarray*}

In a final step, note that when taking the normal derivative of $u_h$, we mean the value of this quantity as taken from this side of the cell (for the usual Lagrange elements, derivatives are not continuous across edges). We then rewrite the above formula by exchanging half of the edge integral of cell $K$ with the neighbor cell $K'$, to obtain

\begin{eqnarray*} J(e) &=& \sum_K (f+\Delta u_h, z-\varphi_h)_K - \frac 12 (\partial_n u_h|_K + \partial_{n'} u_h|_{K'}, z-\varphi_h)_{\partial K\backslash \partial\Omega}. \end{eqnarray*}

Using that for the normal vectors $n'=-n$ holds, we define the jump of the normal derivative by

\[ [\partial_n u_h] := \partial_n u_h|_K + \partial_{n'} u_h|_{K'} = \partial_n u_h|_K - \partial_n u_h|_{K'}, \]

and get the final form after setting the discrete function $\varphi_h$, which is by now still arbitrary, to the point interpolation of the dual solution, $\varphi_h=I_h z$:

\begin{eqnarray*} J(e) &=& \sum_K (f+\Delta u_h, z-I_h z)_K - \frac 12 ([\partial_n u_h], z-I_h z)_{\partial K\backslash \partial\Omega}. \end{eqnarray*}

With this, we have obtained an exact representation of the error of the finite element discretization with respect to arbitrary (linear) functionals $J(\cdot)$. Its structure is a weighted form of a residual estimator, as both $f+\Delta u_h$ and $[\partial_n u_h]$ are cell and edge residuals that vanish on the exact solution, and $z-I_h z$ are weights indicating how important the residuals on a certain cell is for the evaluation of the given functional. Furthermore, it is a cell-wise quantity, so we can use it as a mesh refinement criterion. The question, is: how to evaluate it? After all, the evaluation requires knowledge of the dual solution $z$, which carries the information about the quantity we want to know to best accuracy.

In some, very special cases, this dual solution is known. For example, if the functional $J(\cdot)$ is the point evaluation, $J(\varphi)=\varphi(x_0)$, then the dual solution has to satisfy

\[ -\Delta z = \delta(x-x_0), \]

with the Dirac delta function on the right hand side, and the dual solution is the Green's function with respect to the point $x_0$. For simple geometries, this function is analytically known, and we could insert it into the error representation formula.

However, we do not want to restrict ourselves to such special cases. Rather, we will compute the dual solution numerically, and approximate $z$ by some numerically obtained $\tilde z$. We note that it is not sufficient to compute this approximation $\tilde z$ using the same method as used for the primal solution $u_h$, since then $\tilde z-I_h \tilde z=0$, and the overall error estimate would be zero. Rather, the approximation $\tilde z$ has to be from a larger space than the primal finite element space. There are various ways to obtain such an approximation (see the cited literature), and we will choose to compute it with a higher order finite element space. While this is certainly not the most efficient way, it is simple since we already have all we need to do that in place, and it also allows for simple experimenting. For more efficient methods, again refer to the given literature, in particular [BR95], [BR03].

With this, we end the discussion of the mathematical side of this program and turn to the actual implementation.

The software

The step-14 example program builds heavily on the techniques already used in the step-13 program. Its implementation of the dual weighted residual error estimator explained above is done by deriving a second class, properly called DualSolver, from the Solver base class, and having a class (WeightedResidual) that joins the two again and controls the solution of the primal and dual problem, and then uses both to compute the error indicator for mesh refinement.

The program continues the modular concept of the previous example, by implementing the dual functional, describing quantity of interest, by an abstract base class, and providing two different functionals which implement this interface. Adding a different quantity of interest is thus simple.

One of the more fundamental differences is the handling of data. A common case is that you develop a program that solves a certain equation, and test it with different right hand sides, different domains, different coefficients and boundary values, etc. Usually, these have to match, so that exact solutions are known, or that their combination makes sense at all.

We demonstrate a way how this can be achieved in a simple, yet very flexible way. We will put everything that belongs to a certain setup into one class, and provide a little C++ mortar around it, so that entire setups (domains, coefficients, right hand sides, etc.) can be exchanged by only changing something in one place.

Going this way a little further, we have also centralized all the other parameters that describe how the program is to work in one place, such as the order of the finite element, the maximal number of degrees of freedom, the evaluation objects that shall be executed on the computed solutions, and so on. This allows for simpler configuration of the program, and we will show in a later program how to use a library class that can handle setting these parameters by reading an input file. The general aim is to reduce the places within a program where one may have to look when wanting to change some parameter, as it has turned out in practice that one forgets where they are as programs grow. Furthermore, putting all options describing what the program does in a certain run into a file (that can be stored with the results) helps repeatability of results more than if the various flags were set somewhere in the program, where their exact values are forgotten after the next change to this place.

Unfortunately, the program has become rather long. While this admittedly reduces its usefulness as an example program, we think that it is a very good starting point for development of a program for other kinds of problems, involving different equations than the Laplace equation treated here. Furthermore, it shows everything that we can show you about our way of a posteriori error estimation, and its structure should make it simple for you to adjust this method to other problems, other functionals, other geometries, coefficients, etc.

The author believes that the present program is his masterpiece among the example programs, regarding the mathematical complexity, as well as the simplicity to add extensions. If you use this program as a basis for your own programs, we would kindly like to ask you to state this fact and the name of the author of the example program, Wolfgang Bangerth, in publications that arise from that, of your program consists in a considerable part of the example program.

Bibliography

[Ban00w]
Wolfgang Bangerth.
Mesh adaptivity and error control for a finite element approximation of the elastic wave equation.
In Alfredo Bermudez, Dolores Gomez, Christophe Hazard, Patrick Joly, and Jean E. Roberts, editors, Proceedings of the Fifth International Conference on Mathematical and Numerical Aspects of Wave Propagation (Waves2000), Santiago de Compostela, Spain, 2000, pages 725–729. SIAM, 2000.

[Ban02]
Wolfgang Bangerth.
Adaptive Finite Element Methods for the Identification of Distributed Coefficient in Partial Differential Equations.
PhD thesis, University of Heidelberg, 2002.

[BR99b]
Wolfgang Bangerth and Rolf Rannacher.
Finite element approximation of the acoustic wave equation: Error control and mesh adaptation.
East–West J. Numer. Math., 7(4):263–282, 1999.

[BR03]
Wolfgang Bangerth and Rolf Rannacher.
Adaptive Finite Element Methods for Differential Equations.
Birkhäuser Verlag, Basel, 2003.

[BR01a]
Wolfgang Bangerth and Rolf Rannacher.
Adaptive finite element techniques for the acoustic wave equation.
J. Comput. Acoustics, 9(2):575–591, 2001.

[BR01]
Roland Becker and Rolf Rannacher.
An optimal control approach to error estimation and mesh adaptation in finite element methods.
Acta Numerica, 10:1–102, 2001.

[Bec95]
Roland Becker.
An Adaptive Finite Element Method for the Incompressible Navier-Stokes Equations on Time-dependent Domains.
Dissertation, Universität Heidelberg, 1995.

[Bec98]
Roland Becker.
Weighted error estimators for finite element approximations of the incompressible Navier-Stokes equations.
Preprint 98-20, SFB 359, Universität Heidelberg, 1998.

[BR96r]
Roland Becker and Rolf Rannacher.
A feed-back approach to error control in finite element methods: Basic analysis and examples.
East–West J. Numer. Math., 4:237–264, 1996.

[BR95]
Roland Becker and Rolf Rannacher.
Weighted a posteriori error control in FE methods.
In H. G. Bock et al., ed.s, ENUMATH 95, pages 621–637, Paris, September 1998. World Scientific Publ., Singapure.
in [enumath97].

[enumath97]
Hans Georg Bock, Franco Brezzi, Roland Glowinsky, Guido Kanschat, Yuri A. Kuznetsov, Jacques Periaux, and Rolf Rannacher, editors.
ENUMATH 97, Proceedings of the 2nd European Conference on Numerical Mathematics and Advanced Applications, Singapore, 1998. World Scientific.

[FK97]
Christian Führer and Guido Kanschat.
A posteriori error control in radiative transfer.
Computing, 58(4):317–334, 1997.

[Har02]
Ralf Hartmann.
Adaptive Finite Element Methods for the Compressible Euler Equations.
PhD thesis, University of Heidelberg, 2002.

[HH01]
Ralf Hartmann and Paul Houston.
Adaptive discontinuous Galerkin finite element methods for nonlinear hyperbolic conservation laws.
SIAM J. Sci. Comput., 24 (2002), pp. 979-1004.

[HH01b]
Ralf Hartmann and Paul Houston.
Adaptive discontinuous Galerkin finite element methods for the compressible Euler equations.
J. Comput. Phys. 183 (2002), pp. 508-532.

[Kan96]
Guido Kanschat.
Parallel and Adaptive Galerkin Methods for Radiative Transfer Problems.
Dissertation, Universität Heidelberg, 1996.

[RS97]
Rolf Rannacher and Franz-Theo Suttmeier.
A feed-back approach to error control in finite element methods: Application to linear elasticity.
Comp. Mech., 19(5):434–446, 1997.

[RS98c]
Rolf Rannacher and Franz-Theo Suttmeier.
A posteriori error control in finite element methods via duality techniques: Application to perfect plasticity.
Comp. Mech., 21(2):123–133, 1998.

[RS99]
Rolf Rannacher and Franz-Theo Suttmeier.
A posteriori error control and mesh adaptation for finite element models in elasticity and elasto-plasticity.
Comput. Methods Appl. Mech. Engrg., pages 333–361, 1999.

[Sut96]
Franz-Theo Suttmeier.
Adaptive Finite Element Approximation of Problems in Elasto-Plasticity Theory.
Dissertation, Universität Heidelberg, 1996.

The commented program

Start out with well known things...

 #include <base/quadrature_lib.h>
 #include <base/function.h>
 #include <base/logstream.h>
 #include <base/thread_management.h>
 #include <lac/vector.h>
 #include <lac/full_matrix.h>
 #include <lac/sparse_matrix.h>
 #include <lac/solver_cg.h>
 #include <lac/precondition.h>
 #include <grid/tria.h>
 #include <grid/grid_generator.h>
 #include <grid/grid_out.h>
 #include <grid/tria_accessor.h>
 #include <grid/tria_iterator.h>
 #include <grid/grid_refinement.h>
 #include <dofs/dof_handler.h>
 #include <dofs/dof_constraints.h>
 #include <dofs/dof_accessor.h>
 #include <dofs/dof_tools.h>
 #include <fe/fe_q.h>
 #include <fe/fe_values.h>
 #include <fe/fe_tools.h>
 #include <numerics/vectors.h>
 #include <numerics/matrices.h>
 #include <numerics/data_out.h>
 #include <numerics/error_estimator.h>
 
 #include <iostream>
 #include <fstream>
 #include <list>
 #include <algorithm>
 #include <numeric>
 #include <sstream>

The last step is as in all previous programs:

 using namespace dealii;

Evaluating the solution

As mentioned in the introduction, significant parts of the program have simply been taken over from the step-13 example program. We therefore only comment on those things that are new.

First, the framework for evaluation of solutions is unchanged, i.e. the base class is the same, and the class to evaluate the solution at a grid point is unchanged:

 namespace Evaluation
 {

The EvaluationBase class

   template <int dim>
   class EvaluationBase 
   {
     public:
       virtual ~EvaluationBase ();
 
       void set_refinement_cycle (const unsigned int refinement_cycle);
       
       virtual void operator () (const DoFHandler<dim> &dof_handler,
                                 const Vector<double>  &solution) const = 0;
     protected:
       unsigned int refinement_cycle;
   };
 
 
   template <int dim>
   EvaluationBase<dim>::~EvaluationBase ()
   {}
   
 
   
   template <int dim>
   void
   EvaluationBase<dim>::set_refinement_cycle (const unsigned int step)
   {
     refinement_cycle = step;
   }

The PointValueEvaluation class

   template <int dim>
   class PointValueEvaluation : public EvaluationBase<dim>
   {
     public:
       PointValueEvaluation (const Point<dim>   &evaluation_point);
       
       virtual void operator () (const DoFHandler<dim> &dof_handler,
                                 const Vector<double>  &solution) const;
       
       DeclException1 (ExcEvaluationPointNotFound,
                       Point<dim>,
                       << "The evaluation point " << arg1
                       << " was not found among the vertices of the present grid.");
     private:
       const Point<dim>  evaluation_point;
   };
 
 
   template <int dim>
   PointValueEvaluation<dim>::
   PointValueEvaluation (const Point<dim>   &evaluation_point)
                   :
                   evaluation_point (evaluation_point)
   {}
   
 
 
   template <int dim>
   void
   PointValueEvaluation<dim>::
   operator () (const DoFHandler<dim> &dof_handler,
                const Vector<double>  &solution) const 
   {
     double point_value = 1e20;
 
     typename DoFHandler<dim>::active_cell_iterator
       cell = dof_handler.begin_active(),
       endc = dof_handler.end();
     bool evaluation_point_found = false;
     for (; (cell!=endc) && !evaluation_point_found; ++cell)
       for (unsigned int vertex=0;
            vertex<GeometryInfo<dim>::vertices_per_cell;
            ++vertex)
         if (cell->vertex(vertex).distance (evaluation_point)
             <
             cell->diameter() * 1e-8)
           {
             point_value = solution(cell->vertex_dof_index(vertex,0));
 
             evaluation_point_found = true;
             break;
           };
 
     AssertThrow (evaluation_point_found,
                  ExcEvaluationPointNotFound(evaluation_point));
 
     std::cout << "   Point value=" << point_value
               << std::endl;
   }

The PointXDerivativeEvaluation class

Besides the class implementing the evaluation of the solution at one point, we here provide one which evaluates the gradient at a grid point. Since in general the gradient of a finite element function is not continuous at a vertex, we have to be a little bit more careful here. What we do is to loop over all cells, even if we have found the point already on one cell, and use the mean value of the gradient at the vertex taken from all adjacent cells.

Given the interface of the PointValueEvaluation class, the declaration of this class provides little surprise, and neither does the constructor:

   template <int dim>
   class PointXDerivativeEvaluation : public EvaluationBase<dim>
   {
     public:
       PointXDerivativeEvaluation (const Point<dim>   &evaluation_point);
       
       virtual void operator () (const DoFHandler<dim> &dof_handler,
                                 const Vector<double>  &solution) const;
       
       DeclException1 (ExcEvaluationPointNotFound,
                       Point<dim>,
                       << "The evaluation point " << arg1
                       << " was not found among the vertices of the present grid.");
     private:
       const Point<dim>  evaluation_point;
   };
 
 
   template <int dim>
   PointXDerivativeEvaluation<dim>::
   PointXDerivativeEvaluation (const Point<dim>   &evaluation_point)
                   :
                   evaluation_point (evaluation_point)
   {}

The more interesting things happen inside the function doing the actual evaluation:

   template <int dim>
   void
   PointXDerivativeEvaluation<dim>::
   operator () (const DoFHandler<dim> &dof_handler,
                const Vector<double>  &solution) const 
   {

This time initialize the return value with something useful, since we will have to add up a number of contributions and take the mean value afterwards...

     double point_derivative = 0;

...then have some objects of which the meaning wil become clear below...

     QTrapez<dim>  vertex_quadrature;
     FEValues<dim> fe_values (dof_handler.get_fe(),
                              vertex_quadrature,
                              update_gradients | update_quadrature_points);
     std::vector<Tensor<1,dim> >
       solution_gradients (vertex_quadrature.size());

...and next loop over all cells and their vertices, and count how often the vertex has been found:

     typename DoFHandler<dim>::active_cell_iterator
       cell = dof_handler.begin_active(),
       endc = dof_handler.end();
     unsigned int evaluation_point_hits = 0;
     for (; cell!=endc; ++cell)
       for (unsigned int vertex=0;
            vertex<GeometryInfo<dim>::vertices_per_cell;
            ++vertex)
         if (cell->vertex(vertex) == evaluation_point)
           {

Things are now no more as simple, since we can't get the gradient of the finite element field as before, where we simply had to pick one degree of freedom at a vertex.

Rather, we have to evaluate the finite element field on this cell, and at a certain point. As you know, evaluating finite element fields at certain points is done through the FEValues class, so we use that. The question is: the FEValues object needs to be a given a quadrature formula and can then compute the values of finite element quantities at the quadrature points. Here, we don't want to do quadrature, we simply want to specify some points!

Nevertheless, the same way is chosen: use a special quadrature rule with points at the vertices, since these are what we are interested in. The appropriate rule is the trapezoidal rule, so that is the reason why we used that one above.

Thus: initialize the FEValues object on this cell,

             fe_values.reinit (cell);

and extract the gradients of the solution vector at the vertices:

             fe_values.get_function_grads (solution,
                                           solution_gradients);

Now we have the gradients at all vertices, so pick out that one which belongs to the evaluation point (note that the order of vertices is not necessarily the same as that of the quadrature points):

             unsigned int q_point = 0;
             for (; q_point<solution_gradients.size(); ++q_point)
               if (fe_values.quadrature_point(q_point) ==
                   evaluation_point)
                 break;

Check that the evaluation point was indeed found,

             Assert (q_point < solution_gradients.size(),
                     ExcInternalError());

and if so take the x-derivative of the gradient there as the value which we are interested in, and increase the counter indicating how often we have added to that variable:

             point_derivative += solution_gradients[q_point][0];
             ++evaluation_point_hits;

Finally break out of the innermost loop iterating over the vertices of the present cell, since if we have found the evaluation point at one vertex it cannot be at a following vertex as well:

             break;
           };

Now we have looped over all cells and vertices, so check whether the point was found:

     AssertThrow (evaluation_point_hits > 0,
                  ExcEvaluationPointNotFound(evaluation_point));

We have simply summed up the contributions of all adjacent cells, so we still have to compute the mean value. Once this is done, report the status:

     point_derivative /= evaluation_point_hits;
     std::cout << "   Point x-derivative=" << point_derivative
               << std::endl;
   }

The GridOutput class

Since this program has a more difficult structure (it computed a dual solution in addition to a primal one), writing out the solution is no more done by an evaluation object since we want to write both solutions at once into one file, and that requires some more information than available to the evaluation classes.

However, we also want to look at the grids generated. This again can be done with one such class. Its structure is analog to the SolutionOutput class of the previous example program, so we do not discuss it here in more detail. Furthermore, everything that is used here has already been used in previous example programs.

   template <int dim>
   class GridOutput : public EvaluationBase<dim>
   {
     public:
       GridOutput (const std::string &output_name_base);
       
       virtual void operator () (const DoFHandler<dim> &dof_handler,
                                 const Vector<double>  &solution) const;
     private:
       const std::string output_name_base;
   };
 
 
   template <int dim>
   GridOutput<dim>::
   GridOutput (const std::string &output_name_base)
                   :
                   output_name_base (output_name_base)
   {}
   
 
   template <int dim>
   void
   GridOutput<dim>::operator () (const DoFHandler<dim> &dof_handler,
                                 const Vector<double>  &/ *solution* /) const
   {
     std::ostringstream filename;
     filename << output_name_base << "-"
              << this->refinement_cycle
              << ".eps"
              << std::ends;
     
     std::ofstream out (filename.str().c_str());
     GridOut().write_eps (dof_handler.get_tria(), out);
   }
 }

The Laplace solver classes

Next are the actual solver classes. Again, we discuss only the differences to the previous program.

 namespace LaplaceSolver
 {

Before everything else, forward-declare one class that we will have later, since we will want to make it a friend of some of the classes that follow, which requires the class to be known:

   template <int dim> class WeightedResidual;

The Laplace solver base class

This class is almost unchanged, with the exception that it declares two more functions: output_solution will be used to generate output files from the actual solutions computed by derived classes, and the set_refinement_cycle function by which the testing framework sets the number of the refinement cycle to a local variable in this class; this number is later used to generate filenames for the solution output.

   template <int dim>
   class Base
   {
     public:
       Base (Triangulation<dim> &coarse_grid);
       virtual ~Base ();
 
       virtual void solve_problem () = 0;
       virtual void postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const = 0;
       virtual void refine_grid () = 0;
       virtual unsigned int n_dofs () const = 0;
 
       virtual void set_refinement_cycle (const unsigned int cycle);
 
       virtual void output_solution () const = 0;
       
     protected:
       const SmartPointer<Triangulation<dim> > triangulation;
 
       unsigned int refinement_cycle;
   };
 
 
   template <int dim>
   Base<dim>::Base (Triangulation<dim> &coarse_grid)
                   :
                   triangulation (&coarse_grid)
   {}
 
 
   template <int dim>
   Base<dim>::~Base () 
   {}
 
 
 
   template <int dim>
   void
   Base<dim>::set_refinement_cycle (const unsigned int cycle)
   {
     refinement_cycle = cycle;
   }

The Laplace Solver class

Likewise, the Solver class is entirely unchanged and will thus not be discussed.

   template <int dim>
   class Solver : public virtual Base<dim>
   {
     public:
       Solver (Triangulation<dim>       &triangulation,
               const FiniteElement<dim> &fe,
               const Quadrature<dim>    &quadrature,
               const Quadrature<dim-1>  &face_quadrature,              
               const Function<dim>      &boundary_values);
       virtual
       ~Solver ();
 
       virtual
       void
       solve_problem ();
 
       virtual
       void
       postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const;
 
       virtual
       unsigned int
       n_dofs () const;
       
     protected:
       const SmartPointer<const FiniteElement<dim> >  fe;
       const SmartPointer<const Quadrature<dim> >     quadrature;
       const SmartPointer<const Quadrature<dim-1> >   face_quadrature;      
       DoFHandler<dim>                                dof_handler;
       Vector<double>                                 solution;
       const SmartPointer<const Function<dim> >       boundary_values;
 
       virtual void assemble_rhs (Vector<double> &rhs) const = 0;
     
     private:
       struct LinearSystem
       {
           LinearSystem (const DoFHandler<dim> &dof_handler);
 
           void solve (Vector<double> &solution) const;
         
           ConstraintMatrix     hanging_node_constraints;
           SparsityPattern      sparsity_pattern;
           SparseMatrix<double> matrix;
           Vector<double>       rhs;
       };
 
       void
       assemble_linear_system (LinearSystem &linear_system);
 
       void
       assemble_matrix (LinearSystem                                         &linear_system,
                        const typename DoFHandler<dim>::active_cell_iterator &begin_cell,
                        const typename DoFHandler<dim>::active_cell_iterator &end_cell,
                        Threads::ThreadMutex                                 &mutex) const;
   };
 
 
 
   template <int dim>
   Solver<dim>::Solver (Triangulation<dim>       &triangulation,
                        const FiniteElement<dim> &fe,
                        const Quadrature<dim>    &quadrature,
                        const Quadrature<dim-1>  &face_quadrature,
                        const Function<dim>      &boundary_values)
                   :
                   Base<dim> (triangulation),
                   fe (&fe),
                   quadrature (&quadrature),
                   face_quadrature (&face_quadrature),    
                   dof_handler (triangulation),
                   boundary_values (&boundary_values)
   {}
 
 
   template <int dim>
   Solver<dim>::~Solver () 
   {
     dof_handler.clear ();
   }
 
 
   template <int dim>
   void
   Solver<dim>::solve_problem ()
   {
     dof_handler.distribute_dofs (*fe);
     solution.reinit (dof_handler.n_dofs());
 
     LinearSystem linear_system (dof_handler);
     assemble_linear_system (linear_system);
     linear_system.solve (solution);
   }
 
 
   template <int dim>
   void
   Solver<dim>::
   postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const
   {
     postprocessor (dof_handler, solution);
   }
 
 
   template <int dim>
   unsigned int
   Solver<dim>::n_dofs () const
   {
     return dof_handler.n_dofs();
   }
   
 
   template <int dim>
   void
   Solver<dim>::assemble_linear_system (LinearSystem &linear_system)
   {
     typedef
       typename DoFHandler<dim>::active_cell_iterator
       active_cell_iterator;
 
     const unsigned int n_threads = multithread_info.n_default_threads;
     std::vector<std::pair<active_cell_iterator,active_cell_iterator> >
       thread_ranges 
       = Threads::split_range<active_cell_iterator> (dof_handler.begin_active (),
                                                     dof_handler.end (),
                                                     n_threads);
 
     Threads::ThreadMutex mutex;
     Threads::ThreadGroup<> threads;
     for (unsigned int thread=0; thread<n_threads; ++thread)
       threads += Threads::spawn (*this, &Solver<dim>::assemble_matrix)
                  (linear_system,
                   thread_ranges[thread].first,
                   thread_ranges[thread].second,
                   mutex);
 
     assemble_rhs (linear_system.rhs);
     linear_system.hanging_node_constraints.condense (linear_system.rhs);
 
     std::map<unsigned int,double> boundary_value_map;
     VectorTools::interpolate_boundary_values (dof_handler,
                                               0,
                                               *boundary_values,
                                               boundary_value_map);
     
     threads.join_all ();
     linear_system.hanging_node_constraints.condense (linear_system.matrix);
 
     MatrixTools::apply_boundary_values (boundary_value_map,
                                         linear_system.matrix,
                                         solution,
                                         linear_system.rhs);
   }
 
 
   template <int dim>
   void
   Solver<dim>::assemble_matrix (LinearSystem                                         &linear_system,
                                 const typename DoFHandler<dim>::active_cell_iterator &begin_cell,
                                 const typename DoFHandler<dim>::active_cell_iterator &end_cell,
                                 Threads::ThreadMutex                                 &mutex) const
   {
     FEValues<dim> fe_values (*fe, *quadrature, 
                              update_gradients | update_JxW_values);
 
     const unsigned int   dofs_per_cell = fe->dofs_per_cell;
     const unsigned int   n_q_points    = quadrature->size();
 
     FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
 
     std::vector<unsigned int> local_dof_indices (dofs_per_cell);
 
     for (typename DoFHandler<dim>::active_cell_iterator cell=begin_cell;
          cell!=end_cell; ++cell)
       {
         cell_matrix = 0;
 
         fe_values.reinit (cell);
 
         for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
           for (unsigned int i=0; i<dofs_per_cell; ++i)
             for (unsigned int j=0; j<dofs_per_cell; ++j)
               cell_matrix(i,j) += (fe_values.shape_grad(i,q_point) *
                                    fe_values.shape_grad(j,q_point) *
                                    fe_values.JxW(q_point));
 
 
         cell->get_dof_indices (local_dof_indices);
         Threads::ThreadMutex::ScopedLock lock (mutex);
         for (unsigned int i=0; i<dofs_per_cell; ++i)
           for (unsigned int j=0; j<dofs_per_cell; ++j)
             linear_system.matrix.add (local_dof_indices[i],
                                       local_dof_indices[j],
                                       cell_matrix(i,j));
       };
   }
 
 
   template <int dim>
   Solver<dim>::LinearSystem::
   LinearSystem (const DoFHandler<dim> &dof_handler)
   {
     hanging_node_constraints.clear ();
 
     void (*mhnc_p) (const DoFHandler<dim> &,
                     ConstraintMatrix      &)
       = &DoFTools::make_hanging_node_constraints;
     
     Threads::Thread<>
       mhnc_thread = Threads::spawn (mhnc_p)(dof_handler, hanging_node_constraints);
 
     sparsity_pattern.reinit (dof_handler.n_dofs(),
                              dof_handler.n_dofs(),
                              dof_handler.max_couplings_between_dofs());
     DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
 
     mhnc_thread.join ();
     hanging_node_constraints.close ();
     hanging_node_constraints.condense (sparsity_pattern);
 
     sparsity_pattern.compress();
     matrix.reinit (sparsity_pattern);
     rhs.reinit (dof_handler.n_dofs());
   }
 
 
 
   template <int dim>
   void
   Solver<dim>::LinearSystem::solve (Vector<double> &solution) const
   {
     SolverControl           solver_control (5000, 1e-12);
     SolverCG<>              cg (solver_control);
 
     PreconditionSSOR<> preconditioner;
     preconditioner.initialize(matrix, 1.2);
 
     cg.solve (matrix, solution, rhs, preconditioner);
 
     hanging_node_constraints.distribute (solution);
   }

The PrimalSolver class

The PrimalSolver class is also mostly unchanged except for overloading the functions solve_problem, n_dofs, and postprocess of the base class, and implementing the output_solution function. These overloaded functions do nothing particular besides calling the functions of the base class -- that seems superfluous, but works around a bug in a popular compiler which requires us to write such functions for the following scenario: Besides the PrimalSolver class, we will have a DualSolver, both derived from Solver. We will then have a final classes which derived from these two, which will then have two instances of the Solver class as its base classes. If we want, for example, the number of degrees of freedom of the primal solver, we would have to indicate this like so: PrimalSolver<dim>n_dofs(). However, the compiler does not accept this since the n_dofs function is actually from a base class of the PrimalSolver class, so we have to inject the name from the base to the derived class using these additional functions.

Regarding the implementation of the output_solution function, we keep the GlobalRefinement and RefinementKelly classes in this program, and they can then rely on the default implementation of this function which simply outputs the primal solution. The class implementing dual weighted error estimators will overload this function itself, to also output the dual solution.

Except for this, the class is unchanged with respect to the previous example.

   template <int dim>
   class PrimalSolver : public Solver<dim>
   {
     public:
       PrimalSolver (Triangulation<dim>       &triangulation,
                     const FiniteElement<dim> &fe,
                     const Quadrature<dim>    &quadrature,
                     const Quadrature<dim-1>  &face_quadrature,
                     const Function<dim>      &rhs_function,
                     const Function<dim>      &boundary_values);
 
       virtual
       void solve_problem ();
       
       virtual
       unsigned int n_dofs () const;
       
       virtual
       void postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const;
 
       virtual
       void output_solution () const;
       
     protected:
       const SmartPointer<const Function<dim> > rhs_function;
       virtual void assemble_rhs (Vector<double> &rhs) const;

Now, in order to work around some problems in one of the compilers this library can be compiled with, we will have to declare a class that is actually derived from the present one, as a friend (strange as that seems). The full rationale will be explained below.

       friend class WeightedResidual<dim>;
   };
 
 
   template <int dim>
   PrimalSolver<dim>::
   PrimalSolver (Triangulation<dim>       &triangulation,
                 const FiniteElement<dim> &fe,
                 const Quadrature<dim>    &quadrature,
                 const Quadrature<dim-1>  &face_quadrature,
                 const Function<dim>      &rhs_function,
                 const Function<dim>      &boundary_values)
                   :
                   Base<dim> (triangulation),
                   Solver<dim> (triangulation, fe,
                                quadrature, face_quadrature,
                                boundary_values),
                   rhs_function (&rhs_function)
   {}
 
 
   template <int dim>
   void
   PrimalSolver<dim>::solve_problem ()
   {
     Solver<dim>::solve_problem ();
   }
 
 
 
   template <int dim>
   unsigned int
   PrimalSolver<dim>::n_dofs() const
   {
     return Solver<dim>::n_dofs();
   }
 
 
   template <int dim>
   void
   PrimalSolver<dim>::
   postprocess (const Evaluation::EvaluationBase<dim> &postprocessor) const
   {
     Solver<dim>::postprocess(postprocessor);
   }
 
 
   template <int dim>
   void
   PrimalSolver<dim>::output_solution () const
   {
     DataOut<dim> data_out;
     data_out.attach_dof_handler (this->dof_handler);
     data_out.add_data_vector (this->solution, "solution");
     data_out.build_patches ();
 
     std::ostringstream filename;
     filename << "solution-"
              << this->refinement_cycle
              << ".gnuplot"
              << std::ends;
     
     std::ofstream out (filename.str().c_str());
     data_out.write (out, DataOut<dim>::gnuplot);
   }
   
 
 
   template <int dim>
   void
   PrimalSolver<dim>::
   assemble_rhs (Vector<double> &rhs) const 
   {
     FEValues<dim> fe_values (*this->fe, *this->quadrature, 
                              update_values  | update_quadrature_points  |
                              update_JxW_values);
 
     const unsigned int   dofs_per_cell = this->fe->dofs_per_cell;
     const unsigned int   n_q_points    = this->quadrature->size();
 
     Vector<double>       cell_rhs (dofs_per_cell);
     std::vector<double>  rhs_values (n_q_points);
     std::vector<unsigned int> local_dof_indices (dofs_per_cell);
 
     typename DoFHandler<dim>::active_cell_iterator
       cell = this->dof_handler.begin_active(),
       endc = this->dof_handler.end();
     for (; cell!=endc; ++cell)
       {
         cell_rhs = 0;
 
         fe_values.reinit (cell);
 
         rhs_function->value_list (fe_values.get_quadrature_points(),
                                   rhs_values);
       
         for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
           for (unsigned int i=0; i<dofs_per_cell; ++i)
             cell_rhs(i) += (fe_values.shape_value(i,q_point) *
                             rhs_values[q_point] *
                             fe_values.JxW(q_point));
 
         cell->get_dof_indices (local_dof_indices);
         for (unsigned int i=0; i<dofs_per_cell; ++i)
           rhs(local_dof_indices[i]) += cell_rhs(i);
       }
   }

The RefinementGlobal and RefinementKelly classes

For the following two classes, the same applies as for most of the above: the class is taken from the previous example as-is:

   template <int dim>
   class RefinementGlobal : public PrimalSolver<dim>
   {
     public:
       RefinementGlobal (Triangulation<dim>       &coarse_grid,
                         const FiniteElement<dim> &fe,
                         const Quadrature<dim>    &quadrature,
                         const Quadrature<dim-1>  &face_quadrature,
                         const Function<dim>      &rhs_function,
                         const Function<dim>      &boundary_values);
 
       virtual void refine_grid ();
   };
 
 
 
   template <int dim>
   RefinementGlobal<dim>::
   RefinementGlobal (Triangulation<dim>       &coarse_grid,
                     const FiniteElement<dim> &fe,
                     const Quadrature<dim>    &quadrature,
                     const Quadrature<dim-1>  &face_quadrature,
                     const Function<dim>      &rhs_function,
                     const Function<dim>      &boundary_values)
                   :
                   Base<dim> (coarse_grid),
                   PrimalSolver<dim> (coarse_grid, fe, quadrature,
                                      face_quadrature, rhs_function,
                                      boundary_values)
   {}
 
 
 
   template <int dim>
   void
   RefinementGlobal<dim>::refine_grid ()
   {
     this->triangulation->refine_global (1);
   }
 
 
 
   template <int dim>
   class RefinementKelly : public PrimalSolver<dim>
   {
     public:
       RefinementKelly (Triangulation<dim>       &coarse_grid,
                        const FiniteElement<dim> &fe,
                        const Quadrature<dim>    &quadrature,
                        const Quadrature<dim-1>  &face_quadrature,
                        const Function<dim>      &rhs_function,
                        const Function<dim>      &boundary_values);
 
       virtual void refine_grid ();
   };
 
 
 
   template <int dim>
   RefinementKelly<dim>::
   RefinementKelly (Triangulation<dim>       &coarse_grid,
                    const FiniteElement<dim> &fe,
                    const Quadrature<dim>    &quadrature,
                    const Quadrature<dim-1>  &face_quadrature,
                    const Function<dim>      &rhs_function,
                    const Function<dim>      &boundary_values)
                   :
                   Base<dim> (coarse_grid),
                   PrimalSolver<dim> (coarse_grid, fe, quadrature,
                                      face_quadrature,
                                      rhs_function, boundary_values)
   {}
 
 
 
   template <int dim>
   void
   RefinementKelly<dim>::refine_grid ()
   {
     Vector<float> estimated_error_per_cell (this->triangulation->n_active_cells());
     KellyErrorEstimator<dim>::estimate (this->dof_handler,
                                         QGauss<dim-1>(3),
                                         typename FunctionMap<dim>::type(),
                                         this->solution,
                                         estimated_error_per_cell);
     GridRefinement::refine_and_coarsen_fixed_number (*this->triangulation,
                                                      estimated_error_per_cell,
                                                      0.3, 0.03);
     this->triangulation->execute_coarsening_and_refinement ();
   }

The RefinementWeightedKelly class

This class is a variant of the previous one, in that it allows to weight the refinement indicators we get from the library's Kelly indicator by some function. We include this class since the goal of this example program is to demonstrate automatic refinement criteria even for complex output quantities such as point values or stresses. If we did not solve a dual problem and compute the weights thereof, we would probably be tempted to give a hand-crafted weighting to the indicators to account for the fact that we are going to evaluate these quantities. This class accepts such a weighting function as argument to its constructor:

   template <int dim>
   class RefinementWeightedKelly : public PrimalSolver<dim>
   {
     public:
       RefinementWeightedKelly (Triangulation<dim>       &coarse_grid,
                                const FiniteElement<dim> &fe,
                                const Quadrature<dim>    &quadrature,
                                const Quadrature<dim-1>  &face_quadrature,
                                const Function<dim>      &rhs_function,
                                const Function<dim>      &boundary_values,
                                const Function<dim>      &weighting_function);
 
       virtual void refine_grid ();
 
     private:
       const SmartPointer<const Function<dim> > weighting_function;
   };
 
 
 
   template <int dim>
   RefinementWeightedKelly<dim>::
   RefinementWeightedKelly (Triangulation<dim>       &coarse_grid,
                            const FiniteElement<dim> &fe,
                            const Quadrature<dim>    &quadrature,
                            const Quadrature<dim-1>  &face_quadrature,
                            const Function<dim>      &rhs_function,
                            const Function<dim>      &boundary_values,
                            const Function<dim>      &weighting_function)
                   :
                   Base<dim> (coarse_grid),
                   PrimalSolver<dim> (coarse_grid, fe, quadrature,
                                      face_quadrature,
                                      rhs_function, boundary_values),
                   weighting_function (&weighting_function)
   {}

Now, here comes the main function, including the weighting:

   template <int dim>
   void
   RefinementWeightedKelly<dim>::refine_grid ()
   {

First compute some residual based error indicators for all cells by a method already implemented in the library. What exactly is computed can be read in the documentation of that class.

     Vector<float> estimated_error (this->triangulation->n_active_cells());
     KellyErrorEstimator<dim>::estimate (this->dof_handler,
                                         *this->face_quadrature,
                                         typename FunctionMap<dim>::type(),
                                         this->solution,
                                         estimated_error);

Now we are going to weight these indicators by the value of the function given to the constructor:

     typename DoFHandler<dim>::active_cell_iterator
       cell = this->dof_handler.begin_active(),
       endc = this->dof_handler.end();
     for (unsigned int cell_index=0; cell!=endc; ++cell, ++cell_index)
       estimated_error(cell_index)
         *= weighting_function->value (cell->center());
     
     GridRefinement::refine_and_coarsen_fixed_number (*this->triangulation,
                                                      estimated_error,
                                                      0.3, 0.03);
     this->triangulation->execute_coarsening_and_refinement ();
   }
 
 }

Equation data

In this example program, we work with the same data sets as in the previous one, but as it may so happen that someone wants to run the program with different boundary values and right hand side functions, or on a different grid, we show a simple technique to do exactly that. For more clarity, we furthermore pack everything that has to do with equation data into a namespace of its own.

The underlying assumption is that this is a research program, and that there we often have a number of test cases that consist of a domain, a right hand side, boundary values, possibly a specified coefficient, and a number of other parameters. They often vary all at the same time when shifting from one example to another. To make handling such sets of problem description parameters simple is the goal of the following.

Basically, the idea is this: let us have a structure for each set of data, in which we pack everything that describes a test case: here, these are two subclasses, one called BoundaryValues for the boundary values of the exact solution, and one called RightHandSide, and then a way to generate the coarse grid. Since the solution of the previous example program looked like curved ridges, we use this name here for the enclosing class. Note that the names of the two inner classes have to be the same for all enclosing test case classes, and also that we have attached the dimension template argument to the enclosing class rather than to the inner ones, to make further processing simpler. (From a language viewpoint, a namespace would be better to encapsulate these inner classes, rather than a structure. However, namespaces cannot be given as template arguments, so we use a structure to allow a second object to select from within its given argument. The enclosing structure, of course, has no member variables apart from the classes it declares, and a static function to generate the coarse mesh; it will in general never be instantiated.)

The idea is then the following (this is the right time to also take a brief look at the code below): we can generate objects for boundary values and right hand side by simply giving the name of the outer class as a template argument to a class which we call here Data::SetUp, and it then creates objects for the inner classes. In this case, to get all that characterizes the curved ridge solution, we would simply generate an instance of Data::SetUp<Data::CurvedRidge>, and everything we need to know about the solution would be static member variables and functions of that object.

This approach might seem like overkill in this case, but will become very handy once a certain set up is not only characterized by Dirichlet boundary values and a right hand side function, but in addition by material properties, Neumann values, different boundary descriptors, etc. In that case, the SetUp class might consist of a dozen or more objects, and each descriptor class (like the CurvedRidges class below) would have to provide them. Then, you will be happy to be able to change from one set of data to another by only changing the template argument to the SetUp class at one place, rather than at many.

With this framework for different test cases, we are almost finished, but one thing remains: by now we can select statically, by changing one template argument, which data set to choose. In order to be able to do that dynamically, i.e. at run time, we need a base class. This we provide in the obvious way, see below, with virtual abstract functions. It forces us to introduce a second template parameter dim which we need for the base class (which could be avoided using some template magic, but we omit that), but that's all.

Adding new testcases is now simple, you don't have to touch the framework classes, only a structure like the CurvedRidges one is needed.

 namespace Data
 {

The SetUpBase and SetUp classes

Based on the above description, the SetUpBase class then looks as follows. To allow using the SmartPointer class with this class, we derived from the Subscriptor class.

   template <int dim>
   struct SetUpBase : public Subscriptor
   {
       virtual
       const Function<dim> &  get_boundary_values () const = 0;
 
       virtual
       const Function<dim> &  get_right_hand_side () const = 0;
 
       virtual
       void create_coarse_grid (Triangulation<dim> &coarse_grid) const = 0;
   };

And now for the derived class that takes the template argument as explained above. For some reason, C++ requires us to define a constructor (which maybe empty), as otherwise a warning is generated that some data is not initialized.

Here we pack the data elements into private variables, and allow access to them through the methods of the base class.

   template <class Traits, int dim>
   struct SetUp : public SetUpBase<dim>
   {
       SetUp () {}
 
       virtual
       const Function<dim> &  get_boundary_values () const;
 
       virtual
       const Function<dim> &  get_right_hand_side () const;
       
 
       virtual
       void create_coarse_grid (Triangulation<dim> &coarse_grid) const;
 
     private:
       static const typename Traits::BoundaryValues boundary_values;
       static const typename Traits::RightHandSide  right_hand_side;
   };

We have to provide definitions for the static member variables of the above class:

   template <class Traits, int dim>
   const typename Traits::BoundaryValues  SetUp<Traits,dim>::boundary_values;
   template <class Traits, int dim>
   const typename Traits::RightHandSide   SetUp<Traits,dim>::right_hand_side;

And definitions of the member functions:

   template <class Traits, int dim>
   const Function<dim> &
   SetUp<Traits,dim>::get_boundary_values () const 
   {
     return boundary_values;
   }
 
 
   template <class Traits, int dim>
   const Function<dim> &
   SetUp<Traits,dim>::get_right_hand_side () const 
   {
     return right_hand_side;
   }
 
 
   template <class Traits, int dim>
   void
   SetUp<Traits,dim>::
   create_coarse_grid (Triangulation<dim> &coarse_grid) const 
   {
     Traits::create_coarse_grid (coarse_grid);
   }

The CurvedRidges class

The class that is used to describe the boundary values and right hand side of the curved ridge problem already used in the step-13 example program is then like so:

   template <int dim>
   struct CurvedRidges
   {
       class BoundaryValues : public Function<dim>
       {
         public:
           BoundaryValues () : Function<dim> () {}
           
           virtual double value (const Point<dim>   &p,
                                 const unsigned int  component) const;
       };
 
 
       class RightHandSide : public Function<dim>
       {
         public:
           RightHandSide () : Function<dim> () {}
           
           virtual double value (const Point<dim>   &p,
                                 const unsigned int  component) const;
       };
 
       static
       void
       create_coarse_grid (Triangulation<dim> &coarse_grid);
   };
   
     
   template <int dim>
   double
   CurvedRidges<dim>::BoundaryValues::
   value (const Point<dim>   &p,
          const unsigned int  / *component* /) const
   {
     double q = p(0);
     for (unsigned int i=1; i<dim; ++i)
       q += std::sin(10*p(i)+5*p(0)*p(0));
     const double exponential = std::exp(q);
     return exponential;
   }
 
 
 
   template <int dim>
   double
   CurvedRidges<dim>::RightHandSide::value (const Point<dim>   &p,
                                            const unsigned int  / *component* /) const
   {
     double q = p(0);
     for (unsigned int i=1; i<dim; ++i)
       q += std::sin(10*p(i)+5*p(0)*p(0));
     const double u = std::exp(q);
     double t1 = 1,
            t2 = 0,
            t3 = 0;
     for (unsigned int i=1; i<dim; ++i)
       {
         t1 += std::cos(10*p(i)+5*p(0)*p(0)) * 10 * p(0);
         t2 += 10*std::cos(10*p(i)+5*p(0)*p(0)) -
               100*std::sin(10*p(i)+5*p(0)*p(0)) * p(0)*p(0);
         t3 += 100*std::cos(10*p(i)+5*p(0)*p(0))*std::cos(10*p(i)+5*p(0)*p(0)) -
               100*std::sin(10*p(i)+5*p(0)*p(0));
       };
     t1 = t1*t1;
     
     return -u*(t1+t2+t3);
   }
 
 
   template <int dim>
   void
   CurvedRidges<dim>::
   create_coarse_grid (Triangulation<dim> &coarse_grid)
   {
     GridGenerator::hyper_cube (coarse_grid, -1, 1);
     coarse_grid.refine_global (2);
   }

The Exercise_2_3 class

This example program was written while giving practical courses for a lecture on adaptive finite element methods and duality based error estimates. For these courses, we had one exercise, which required to solve the Laplace equation with constant right hand side on a square domain with a square hole in the center, and zero boundary values. Since the implementation of the properties of this problem is so particularly simple here, lets do it. As the number of the exercise was 2.3, we take the liberty to retain this name for the class as well.

   template <int dim>
   struct Exercise_2_3
   {

We need a class to denote the boundary values of the problem. In this case, this is simple: it's the zero function, so don't even declare a class, just a typedef:

       typedef ZeroFunction<dim> BoundaryValues;

Second, a class that denotes the right hand side. Since they are constant, just subclass the corresponding class of the library and be done:

       class RightHandSide : public ConstantFunction<dim>
       {
         public:
           RightHandSide () : ConstantFunction<dim> (1.) {}
       };

Finally a function to generate the coarse grid. This is somewhat more complicated here, see immediately below.

       static
       void
       create_coarse_grid (Triangulation<dim> &coarse_grid);
   };

As stated above, the grid for this example is the square [-1,1]^2 with the square [-1/2,1/2]^2 as hole in it. We create the coarse grid as 4 times 4 cells with the middle four ones missing.

Of course, the example has an extension to 3d, but since this function cannot be written in a dimension independent way we choose not to implement this here, but rather only specialize the template for dim=2. If you compile the program for 3d, you'll get a message from the linker that this function is not implemented for 3d, and needs to be provided.

For the creation of this geometry, the library has no predefined method. In this case, the geometry is still simple enough to do the creation by hand, rather than using a mesh generator.

   template <>
   void
   Exercise_2_3<2>::
   create_coarse_grid (Triangulation<2> &coarse_grid)
   {

First define the space dimension, to allow those parts of the function that are actually dimension independent to use this variable. That makes it simpler if you later takes this as a starting point to implement the 3d version.

     const unsigned int dim = 2;

Then have a list of vertices. Here, they are 24 (5 times 5, with the middle one omitted). It is probably best to draw a sketch here. Note that we leave the number of vertices open at first, but then let the compiler compute this number afterwards. This reduces the possibility of having the dimension to large and leaving the last ones uninitialized.

     static const Point<2> vertices_1[]
       = {  Point<2> (-1.,   -1.),
              Point<2> (-1./2, -1.),
              Point<2> (0.,    -1.),
              Point<2> (+1./2, -1.),
              Point<2> (+1,    -1.),
              
              Point<2> (-1.,   -1./2.),
              Point<2> (-1./2, -1./2.),
              Point<2> (0.,    -1./2.),
              Point<2> (+1./2, -1./2.),
              Point<2> (+1,    -1./2.),
              
              Point<2> (-1.,   0.),
              Point<2> (-1./2, 0.),
              Point<2> (+1./2, 0.),
              Point<2> (+1,    0.),
              
              Point<2> (-1.,   1./2.),
              Point<2> (-1./2, 1./2.),
              Point<2> (0.,    1./2.),
              Point<2> (+1./2, 1./2.),
              Point<2> (+1,    1./2.),
              
              Point<2> (-1.,   1.),
              Point<2> (-1./2, 1.),
              Point<2> (0.,    1.),                          
              Point<2> (+1./2, 1.),
              Point<2> (+1,    1.)    };
     const unsigned int
       n_vertices = sizeof(vertices_1) / sizeof(vertices_1[0]);

From this static list of vertices, we generate an STL vector of the vertices, as this is the data type the library wants to see.

     const std::vector<Point<dim> > vertices (&vertices_1[0],
                                              &vertices_1[n_vertices]);

Next, we have to define the cells and the vertices they contain. Here, we have 8 vertices, but leave the number open and let it be computed afterwards:

     static const int cell_vertices[][GeometryInfo<dim>::vertices_per_cell]
       = {{0, 1, 5, 6},
          {1, 2, 6, 7},
          {2, 3, 7, 8},
          {3, 4, 8, 9},
          {5, 6, 10, 11},
          {8, 9, 12, 13},
          {10, 11, 14, 15},
          {12, 13, 17, 18},
          {14, 15, 19, 20},
          {15, 16, 20, 21},
          {16, 17, 21, 22},
          {17, 18, 22, 23}};
     const unsigned int
       n_cells = sizeof(cell_vertices) / sizeof(cell_vertices[0]);

Again, we generate a C++ vector type from this, but this time by looping over the cells (yes, this is boring). Additionally, we set the material indicator to zero for all the cells:

     std::vector<CellData<dim> > cells (n_cells, CellData<dim>());
     for (unsigned int i=0; i<n_cells; ++i) 
       {
         for (unsigned int j=0;
              j<GeometryInfo<dim>::vertices_per_cell;
              ++j)
           cells[i].vertices[j] = cell_vertices[i][j];
         cells[i].material_id = 0;
       };

Finally pass all this information to the library to generate a triangulation. The last parameter may be used to pass information about non-zero boundary indicators at certain faces of the triangulation to the library, but we don't want that here, so we give an empty object:

     coarse_grid.create_triangulation (vertices,
                                       cells,
                                       SubCellData());

And since we want that the evaluation point (3/4,3/4) in this example is a grid point, we refine once globally:

     coarse_grid.refine_global (1);
   }
 }

Discussion

As you have now read through this framework, you may be wondering why we have not chosen to implement the classes implementing a certain setup (like the CurvedRidges class) directly as classes derived from Data::SetUpBase. Indeed, we could have done very well so. The only reason is that then we would have to have member variables for the solution and right hand side classes in the CurvedRidges class, as well as member functions overloading the abstract functions of the base class giving access to these member variables. The SetUp class has the sole reason to relieve us from the need to reiterate these member variables and functions that would be necessary in all such classes. In some way, the template mechanism here only provides a way to have default implementations for a number of functions that depend on external quantities and can thus not be provided using normal virtual functions, at least not without the help of templates.

However, there might be good reasons to actually implement classes derived from Data::SetUpBase, for example if the solution or right hand side classes require constructors that take arguments, which the Data::SetUpBase class cannot provide. In that case, subclassing is a worthwhile strategy. Other possibilities for special cases are to derive from Data::SetUp<SomeSetUp> where SomeSetUp denotes a class, or even to explicitly specialize Data::SetUp<SomeSetUp>. The latter allows to transparently use the way the SetUp class is used for other set-ups, but with special actions taken for special arguments.

A final observation favoring the approach taken here is the following: we have found numerous times that when starting a project, the number of parameters (usually boundary values, right hand side, coarse grid, just as here) was small, and the number of test cases was small as well. One then starts out by handcoding them into a number of switch statements. Over time, projects grow, and so does the number of test cases. The number of switch statements grows with that, and their length as well, and one starts to find ways to consider impossible examples where domains, boundary values, and right hand sides do not fit together any more, and starts loosing the overview over the whole structure. Encapsulating everything belonging to a certain test case into a structure of its own has proven worthwhile for this, as it keeps everything that belongs to one test case in one place. Furthermore, it allows to put these things all in one or more files that are only devoted to test cases and their data, without having to bring their actual implementation into contact with the rest of the program.

Dual functionals

As with the other components of the program, we put everything we need to describe dual functionals into a namespace of its own, and define an abstract base class that provides the interface the class solving the dual problem needs for its work.

We will then implement two such classes, for the evaluation of a point value and of the derivative of the solution at that point. For these functionals we already have the corresponding evaluation objects, so they are comlementary.

 namespace DualFunctional
 {

The DualFunctionalBase class

First start with the base class for dual functionals. Since for linear problems the characteristics of the dual problem play a role only in the right hand side, we only need to provide for a function that assembles the right hand side for a given discretization:

   template <int dim>
   class DualFunctionalBase : public Subscriptor
   {
     public:
       virtual
       void
       assemble_rhs (const DoFHandler<dim> &dof_handler,
                     Vector<double>        &rhs) const = 0;
   };

The PointValueEvaluation class

As a first application, we consider the functional corresponding to the evaluation of the solution's value at a given point which again we assume to be a vertex. Apart from the constructor that takes and stores the evaluation point, this class consists only of the function that implements assembling the right hand side.

   template <int dim>
   class PointValueEvaluation : public DualFunctionalBase<dim>
   {
     public:
       PointValueEvaluation (const Point<dim> &evaluation_point);
 
       virtual
       void
       assemble_rhs (const DoFHandler<dim> &dof_handler,
                     Vector<double>        &rhs) const;
       
       DeclException1 (ExcEvaluationPointNotFound,
                       Point<dim>,
                       << "The evaluation point " << arg1
                       << " was not found among the vertices of the present grid.");
 
     protected:
       const Point<dim> evaluation_point;
   };
 
 
   template <int dim>
   PointValueEvaluation<dim>::
   PointValueEvaluation (const Point<dim> &evaluation_point)
                   :
                   evaluation_point (evaluation_point)
   {}

As for doing the main purpose of the class, assembling the right hand side, let us first consider what is necessary: The right hand side of the dual problem is a vector of values J(phi_i), where J is the error functional, and phi_i is the i-th shape function. Here, J is the evaluation at the point x0, i.e. J(phi_i)=phi_i(x0).

Now, we have assumed that the evaluation point is a vertex. Thus, for the usual finite elements we might be using in this program, we can take for granted that at such a point exactly one shape function is nonzero, and in particular has the value one. Thus, we set the right hand side vector to all-zeros, then seek for the shape function associated with that point and set the corresponding value of the right hand side vector to one:

   template <int dim>
   void
   PointValueEvaluation<dim>::
   assemble_rhs (const DoFHandler<dim> &dof_handler,
                 Vector<double>        &rhs) const