"Fossies" - the Fresh Open Source Software Archive

Member "dune-grid-howto-2.7.0/doc/grid-howto.tex" (20 Dec 2019, 121110 Bytes) of package /linux/misc/dune/dune-grid-howto-2.7.0.tar.gz:

As a special service "Fossies" has tried to format the requested source page into HTML format using (guessed) TeX and LaTeX source code syntax highlighting (style: standard) with prefixed line numbers. Alternatively you can here view or download the uninterpreted source code file.

    1 \documentclass[11pt,a4paper,headinclude,footinclude,DIV16,headings=normal]{scrreprt}
    2 \usepackage[automark]{scrpage2}
    3 \usepackage{scrhack}
    4 \usepackage[ansinew]{inputenc}
    5 \usepackage{amsmath}
    6 \usepackage{amsfonts}
    7 \usepackage{theorem}
    8 \usepackage{xspace}
    9 \usepackage{color}
   10 \usepackage{overpic}
   11 \usepackage{listings}
   12 \lstset{language=C++, basicstyle=\ttfamily,
   13   keywordstyle=\color{black}\bfseries, tabsize=4, stringstyle=\ttfamily,
   14   commentstyle=\itshape, extendedchars=true, escapeinside={/*@}{@*/}}
   15 \usepackage{hyperref}
   16 \usepackage{makeidx}
   18 \usepackage{graphicx}
   20 \DeclareGraphicsExtensions{.pdf, .png, .jpg}
   22 \newcommand{\C}{\mathbb{C}}
   23 \newcommand{\R}{\mathbb{R}}
   24 \newcommand{\N}{\mathbb{N}}
   25 \newcommand{\Z}{\mathbb{Z}}
   26 \newcommand{\Q}{\mathbb{Q}}
   27 \newcommand{\Dune}{{\sffamily\bfseries DUNE}\xspace}
   29 %The theorems
   30 \theorembodyfont{\upshape}
   31 \theoremheaderfont{\sffamily\bfseries}
   32 \newtheorem{exc}{Exercise}[chapter]
   33 \newtheorem{rem}[exc]{Remark}
   34 \newtheorem{lst}{Listing}
   35 \newtheorem{warn}[exc]{Warning}
   37 \pagestyle{scrheadings}
   39 \title{The Distributed and Unified Numerics Environment (\Dune{}) Grid
   40   Interface HOWTO}
   42 \input{config.inc}
   44 \author{Peter Bastian$^\ast$ \and
   45 Markus Blatt$^\ast$ \and
   46 Andreas Dedner$^\dagger$ \and
   47 Christian Engwer$^\ast$ \and
   48 Robert Kl\"ofkorn$^\dagger$ \and
   49 Martin Nolte$^\dagger$ \and
   50 Mario Ohlberger$^\mathparagraph$ \and
   51 Oliver Sander$^\ddagger$}
   53 \date{Version \version, \today}
   55 \publishers{%
   56 \vspace{10mm}
   57 %\includegraphics[width=0.32\textwidth]{img/alberta2d-view2}\hfill
   58 \includegraphics[width=0.32\textwidth]{img/ug2dtri-view2}\hfill
   59 \includegraphics[width=0.32\textwidth]{img/adaptiveintegration_alberta2d}\hfill
   60 \includegraphics[width=0.32\textwidth]{img/alucube3d}
   61 %\includegraphics[width=0.32\textwidth]{img/ug2dquad-view2}
   62 \\
   63 \vspace{10mm}
   64 {\normalsize $^\ast$Abteilung `Simulation gro\ss er Systeme',
   65 Universit\"at Stuttgart,\\
   66 Universit\"atsstr.~38, D-70569 Stuttgart, Germany}\\
   67 %
   68 \bigskip
   69 {\normalsize $^\dagger$Abteilung f\"ur Angewandte Mathematik, Universit\"at Freiburg,\\
   70 Hermann-Herder-Str.~10, D-79104 Freiburg, Germany}\\
   71 %
   72 \bigskip
   73 {\normalsize $^\mathparagraph$Institut f\"ur Numerische und Angewandte Mathematik, Universit\"at M\"unster,\\
   74 Einsteinstr.~62, D-48149 M\"unster, Germany}\\
   75 %
   76 \bigskip
   77 {\normalsize $^\ddagger$Institut f\"ur Mathematik II,\\ Freie Universit\"at Berlin,
   78 Arnimallee 6, D-14195 Berlin, Germany}\\
   79 %
   80 \bigskip
   81 {\normalsize \texttt{\url{http://www.dune-project.org/}}}\\
   82 }
   84 \begin{document}
   86 \maketitle
   88 \begin{abstract}
   89 This document gives an introduction to the Distributed and Unified
   90 Numerics Environment (\Dune). \Dune{} is a template library for the
   91 numerical solution of partial differential equations. It is based on
   92 the following principles: i) Separation of data structures and
   93 algorithms by abstract interfaces, ii) Efficient implementation of these
   94 interfaces using generic programming techniques (templates) in C++ and
   95 iii) Reuse of existing finite element packages with a large body of
   96 functionality. This introduction covers only the abstract grid interface
   97 of \Dune{} which is currently the most developed part. However, part of
   98 \Dune{} are also the Iterative Solver Template Library (ISTL, providing a
   99 large variety of solvers for sparse linear systems) and a flexible class
  100 hierarchy for finite element methods. These will be described in
  101 subsequent documents. Now have fun!
  103 \vspace*{\fill}
  104 Thanks to Martin Drohmann for adapting this howto to version 1.2 of the
  105 DUNE grid interface.
  106 \end{abstract}
  108 \tableofcontents
  111 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  112 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  113 \chapter{Introduction}
  114 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  115 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  117 \section{What is \texorpdfstring{\Dune{}}{DUNE} anyway?}
  119 \Dune{} is a software framework for the numerical solution of partial
  120 differential equations with grid-based methods. It is based on the
  121 following main principles:
  122 \begin{itemize}
  123 \item \textit{Separation of data structures and
  124 algorithms by abstract interfaces.} This provides more functionality
  125 with less code and also ensures maintainability and
  126 extendability of the framework.
  127 \item \textit{Efficient implementation of these
  128 interfaces using generic programming techniques}. Static polymorphism
  129 allows the compiler to do more optimizations, in particular function
  130 inlining, which in turn allows the interface to have very small
  131 functions (implemented by one or few machine instructions) without a
  132 severe performance penalty. In essence the algorithms are parametrized
  133 with a particular data structure and the interface is removed at
  134 compile time. Thus the resulting code is as efficient as if it would
  135 have been written for the special case.
  136 \item \textit{Reuse of existing finite element packages with a large body of
  137 functionality.} In particular the finite element codes UG, \cite{ug},
  138 Alberta, \cite{Alberta}, and ALU3d, \cite{ALU3d}, have been
  139 adapted to the \Dune{} framework. Thus, parallel and adaptive meshes with
  140 multiple element types and refinement rules are available. All these
  141 packages can be linked together in one executable.
  142 \end{itemize}
  144 The framework consists of a number of modules which are downloadable
  145 as separate packages.  The current core modules are:
  146 \begin{itemize}
  147 \item \texttt{dune-common} contains the basic classes used by all
  148   \Dune{}-modules.  It provides some infrastructural classes for
  149   debugging and exception handling as well as a library to handle
  150   dense matrices and vectors.
  151 \item \texttt{dune-geometry} provides element geometry classes and
  152   quadrature rules.
  153 \item \texttt{dune-grid} is the grid interface and is covered in
  154   this document. It defines nonconforming, hierarchically nested,
  155   multi-element-type, parallel grids in arbitrary space dimensions.
  156   Graphical output with several packages is available, e.~g.~file
  157   output to VTK (parallel XML format for
  158   unstructured grids). The graphics package Grape, \cite{Grape} has
  159   been integrated in interactive mode.
  160 \item \texttt{dune-istl} -- \textit{Iterative Solver Template
  161     Library.} Provides generic sparse matrix/vector classes and a
  162   variety of solvers based on these classes. A special feature is the
  163   use of templates to exploit the recursive block structure of finite
  164   element matrices at compile time. Available solvers include Krylov
  165   methods, (block-) incomplete decompositions and aggregation-based
  166   algebraic multigrid.
  167  \item \texttt{dune-localfunctions} -- \textit{Library of local base
  168   functions.} Provides classes for base functions on reference elements
  169   from which global discrete function spaces can be constructed.
  170 \end{itemize}
  172 Before starting to work with \Dune{} you might want to update your
  173 knowledge about C++ and templates in particular. For that you should
  174 have the bible, \cite{Stroustrup}, at your desk. A good introduction,
  175 besides its age, is still the book by Barton and Nackman,
  176 \cite{BN}. The definitive guide to template programming is
  177 \cite{VandervoordeJosuttis}. A very useful compilation of template
  178 programming tricks with application to scientific computing is given
  179 in \cite{Veldhui99} (if you can't find it on the web, contact us).
  181 \section{Download}
  183 The source code of the \Dune{} framework can be
  184 downloaded from the web page.  To get started, it is easiest to
  185 download the latest stable version of the tarballs of
  186 \texttt{dune-common}, \texttt{dune-grid} and \texttt{dune-grid-howto}.
  187 These are available on the \Dune{} download page:
  188 %
  189 \begin{center}
  190 \url{http://www.dune-project.org/download.html}
  191 \end{center}
  192 %
  194 Alternatively, you can download the latest development version via
  195 anonymous Git.  For further information, please see the web page.
  197 \section{Installation}
  199 The official installation instructions are available on the web page
  200 %
  201 \begin{center}
  202 \url{http://www.dune-project.org/doc/installation-notes.html}
  203 \end{center}
  205 Obviously, we do not want to copy all this information because it might
  206 get outdated and inconsistent then. To make this document
  207 self-contained, we describe only how to install \Dune{} from the
  208 tarballs.  If you prefer to use the version from Git, see the web page
  209 for further information.  Moreover, we assume that you use a UNIX
  210 system (usually Linux or OS X). If you have the Redmond system then
  211 you have to explore new horizons.
  213 In order to build the \Dune{} framework, you need a standards
  214 compliant C++ compiler.  We tested compiling with GNU \texttt{g++} in
  215 version 4.4 or newer and \texttt{clang} in version 3.4 or newer. Recent versions
  216 of Intel \texttt{icc} ($\geq 15$) should work, too. For older versions of
  217 \texttt{icc}, like $14.0.3$, several issues are known and usage of these
  218 is discouraged.
  220 Now extract the tarballs of \texttt{dune-common},
  221 \texttt{dune-geometry}, \texttt{dune-istl}, \texttt{dune-grid},
  222 and \texttt{dune-grid\-howto} into a common directory, say
  223 \texttt{dune-home}.  Change to this directory and call
  224 \begin{lstlisting}
  225 > dune-common-2.3/bin/dunecontrol all
  226 \end{lstlisting}
  227 Replace ``\texttt{2.3}'' by the actual version number of the package
  228 you downloaded if necessary.  This should configure and build all
  229 \Dune{} modules in \texttt{dune-home} with a basic configuration.
  231 For many of the examples in this howto you need adaptive grids or the
  232 parallel features of \Dune{}.  To use adaptive grids, you need to
  233 install one of the external grid packages which \Dune{} provides
  234 interfaces for, for instance Alberta, UG and ALUGrid.
  235 \begin{itemize}
  236 \item Alberta -- \url{http://www.alberta-fem.de/}
  237 \item UG -- \url{http://www.iwr.uni-heidelberg.de/frame/iwrwikiequipment/software/ug}
  238 \item \Dune{}-ALUGrid -- \url{http://users.dune-project.org/projects/dune-alugrid}
  239 \end{itemize}
  240 To use the parallel code of \Dune{}, you need an implementation of the
  241 Message Passing Interface (MPI), for example MPICH or Open MPI.  For the
  242 \Dune{} build system to find these libraries, the \texttt{configure}
  243 scripts of the particular \Dune{} modules must be passed the locations
  244 of the respective installations.  The \texttt{dunecontrol} script
  245 facilitates to pass options to the \texttt{configure} via a
  246 configuration file.  Such a configuration file might look like this:
  247 \begin{lstlisting}
  248 CONFIGURE_FLAGS="--with-alberta=/path/to/alberta "\
  249 "--with-ug=/path/to/ug --enable-parallel"
  250 MAKE_FLAGS="-j 2"
  251 \end{lstlisting}
  252 If this is saved under the name \texttt{dunecontrol.opts}, you
  253 can tell \texttt{dunecontrol} to consider the file by calling
  254 \begin{lstlisting}
  255 > dune-common-2.3/bin/dunecontrol --opts=dunecontrol.opts all
  256 \end{lstlisting}
  258 For information on how to build and configure the respective grids,
  259 please see the \Dune{} web page.
  261 \section{Code documentation}
  263 Documentation of the files and classes in \Dune{} is provided in code and
  264 can be extracted using the
  265 Doxygen\footnote{\url{http://www.doxygen.org/}}
  266 software available elsewhere. The code documentation can either be built
  267 locally on your machine (in HTML and other formats, e.\,g.~\LaTeX)
  268 or its latest version is available at
  269 \begin{center}
  270 \url{http://www.dune-project.org/doxygen/master/index.html}
  271 \end{center}
  273 \section{Licence}
  275 The \Dune{} library and headers are licensed under version 2 of the
  276 GNU General Public License%
  277 \footnote{\url{http://www.gnu.org/licenses/gpl-2.0.html}}, with a special
  278 exception for linking and compiling against \Dune{}, the so-called
  279 ``runtime exception.''  The license is intended to be similar to the
  280 GNU Lesser General Public License, which by itself isn't suitable for
  281 a C++ template library.
  283 The exact wording of the exception reads as follows:
  285 \begin{quote}
  286    As a special exception, you may use the \Dune{} source files as part
  287    of a software library or application without restriction.
  288    Specifically, if other files instantiate templates or use macros or
  289    inline functions from one or more of the \Dune{} source files, or you
  290    compile one or more of the \Dune{} source files and link them with
  291    other files to produce an executable, this does not by itself cause
  292    the resulting executable to be covered by the GNU General Public
  293    License.  This exception does not however invalidate any other
  294    reasons why the executable file might be covered by the GNU General
  295    Public License.
  296 \end{quote}
  298 %\section{Contributing to DUNE}
  303 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  304 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  305 \chapter{Getting started}
  306 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  307 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  309 In this section we will take a quick tour through the abstract
  310 grid interface provided by \Dune. This should give you an overview of
  311 the different classes before we go into the details.
  313 \section{Creating your first grid}
  315 Let us start with a replacement of the famous \emph{hello world}
  316 program given below.
  318 \begin{lst}[File dune-grid-howto/gettingstarted.cc] \mbox{}
  319 \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
  320 numberstyle=\tiny, numbersep=5pt]{../gettingstarted.cc}
  321 \end{lst}
  323 This program is quite simple. It starts with some includes in lines
  324 \ref{gs:inc0}-\ref{gs:inc1}. The file \lstinline!config.h! has been
  325 produced by the CMake
  326 build system. It contains the current configuration and can be used to
  327 compile different versions of your code depending on the configuration
  328 selected. It is important that this file is included before any other
  329 \Dune{} header files. The file \lstinline!dune/grid/yaspgrid.hh!
  330 includes the headers for the \lstinline!YaspGrid! class which provides a
  331 special implementation of the \Dune{} grid interface with a
  332 structured mesh of arbitrary dimension. Then
  333 \lstinline!dune/grid/common/gridinfo.hh!  loads the headers of some
  334 functions which print useful information about a grid.
  336 Since the dimension will be used as a template parameter in many
  337 places below we define it as a constant in line number \ref{gs:dim}.
  338 The \lstinline!YaspGrid! class template takes its dimension as a
  339 template parameter.  Technically, \lstinline!YaspGrid! has a second template
  340 parameter, that defines how it stores coordinates. In this example, we
  341 use an equidistant grid, which allows us to stay with the default value
  342 \lstinline!EquidistantCoordinates!. For ease of writing we
  343 define in line \ref{gs:gridtype} the type \lstinline!GridType! using
  344 the selected value for the dimension. All identifiers of the \Dune{}
  345 framework are within the \lstinline!Dune! namespace.
  347 Lines \ref{gs:par0}-\ref{gs:par1} prepare the arguments for the
  348 construction of a \lstinline!YaspGrid! object. The first argument use
  349 the class template \lstinline!FieldVector<T,n>! which is a vector with
  350 \lstinline!n!  components of type \lstinline!T!. You can either assign
  351 the same value to all components in the constructor (as is done here)
  352 or you could use \lstinline!operator[]! to assign values to individual
  353 components. The variable \lstinline!length! defines the lengths of
  354 the cube. The variable \lstinline!elements! defines the number of
  355 cells or elements to be used in the respective dimension of the grid.
  356 Finally in line \ref{gs:grid} we are now able to
  357 instantiate the \lstinline!YaspGrid!  object.
  359 The only thing we do with the grid in this little example is printing
  360 some information about it. After successfully running the executable
  361 \lstinline!gettingstarted! you should see an output like this:
  363 \begin{lst}[Output of gettingstarted] \mbox{}
  365 \begin{lstlisting}[basicstyle=\ttfamily\scriptsize]
  366 => Dune::YaspGrid<3, Dune::EquidistantCoordinates<double, 3> > (dim=3, dimworld=3)
  367 level 0 codim[0]=64 codim[1]=240 codim[2]=300 codim[3]=125
  368 leaf    codim[0]=64 codim[1]=240 codim[2]=300 codim[3]=125
  369 leaf dim=3 geomTypes=((cube, 3)[0]=64,(cube, 2)[1]=240,(simplex, 1)[2]=300,(simplex, 0)[3]=125)
  370 \end{lstlisting}
  371 \end{lst}
  373 The first line tells you that you are looking at an \lstinline!YaspGrid!
  374 object of the given dimensions. The \Dune{} grid interface supports
  375 unstructured, locally refined, logically nested grids. The coarsest
  376 grid is called level-0-grid or macro grid. Elements can be
  377 individually refined into a number of smaller elements. Each element
  378 of the macro grid and all its descendants obtained from refinement
  379 form a tree structure. All elements at depth $n$ of a refinement tree
  380 form the level-$n$-grid. All elements that are leaves of a refinement
  381 tree together form the so-called leaf grid. The second line of the
  382 output tells us that this grid object consists only of a single level
  383 (level $0$) while the next line tells us that that level 0 coincides
  384 also with the leaf grid in this case. Each line reports about the
  385 number of grid entities which make up the grid. We see that there are
  386 64 elements (codimension 0), 240 faces (codimension 1), 300 edges
  387 (codimension 2) and 125 vertices (codimension 3) in the grid. The last
  388 line reports on the different types of entities making up the grid. In
  389 this case all entities are of type \emph{cube}, as a line and a point
  390 are both \emph{cube} and \emph{simplex}.
  392 \begin{exc} Try to play around with different grid sizes by assigning
  393   different values to the \lstinline!element! parameter. You can also change
  394   the dimension of the grid by varying \lstinline!dim!. Don't be
  395   modest. Also try dimensions 4 and 5!
  396 \end{exc}
  398 \begin{exc}
  399   The static methods \lstinline!Dune::gridlevellist! and
  400   \lstinline!Dune::gridleaflist! produce a very detailed output of the grid's
  401   elements on a specified grid level. Change the code and print out this
  402   information for the leaf grid or a grid on lower level. Try to understand the
  403   output.
  404 \end{exc}
  406 \section{Traversing a grid --- A first look at the grid interface}
  408 After looking at very first simple example we are now ready to go on
  409 to a more complicated one. Here it is:
  411 \begin{lst}[File dune-grid-howto/traversal.cc] \mbox{}
  412 \nopagebreak
  413 \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
  414 numberstyle=\tiny, numbersep=5pt]{../traversal.cc}
  415 \end{lst}
  417 The \lstinline!main! function near the end of the listing is pretty
  418 similar to the previous one except that we use a 2d grid for the unit
  419 square that just consists of one cell. In line \ref{tc:refine} this
  420 cell is refined once using the standard method of grid refinement of
  421 the implementation. Here, the cell is refined into four smaller cells.
  422 The main work is done in a call to the function \lstinline!traversal!
  423 in line \ref{tc:call}.  This function is given in lines
  424 \ref{tc:tra0}-\ref{tc:tra1}.
  426 The function \lstinline!traversal! is a function template that is
  427 parametrized by a class \lstinline!G! that is assumed to implement
  428 the \Dune{} grid interface.  Thus, it will work on \textit{any} grid
  429 available in \Dune{} without any changes. We now go into the details
  430 of this function.
  432 The algorithm should work in any dimension so we extract the grid's
  433 dimension in line \ref{tc:dim}. Next, each \Dune{} grid defines a type
  434 that it uses to represent positions. This type is extracted in line
  435 \ref{tc:ct} for later use.
  437 A grid is considered to be a container of ``entities'' which are
  438 abstractions for geometric objects like vertices, edges,
  439 quadrilaterals, tetrahedra, and so on. This is very similar to the
  440 standard template library (STL), see e.~g.~\cite{Stroustrup},
  441 which is part of any C++ system.
  442 A key difference is, however, that there is not just one type of entity but
  443 several. As in the STL the elements of any container can be accessed
  444 with iterators which are generalized pointers. Again, a \Dune{} grid
  445 knows several different iterators which provide access to the
  446 different kinds of entities and which also provide different patterns
  447 of access.
  449 As we usually do not want to use the entire hierarchy of the grid, we first
  450 define a view on that part of the grid we are interested in. This can be a
  451 level or the leaf part of the grid. In line \ref{tc:lfgv} a type for a
  452 \lstinline!GridView! on the leaf grid is defined.
  454 Line \ref{tc:ittype} extracts the type of an iterator from this view
  455 class. \lstinline!Codim! is a \lstinline!struct! within the grid class
  456 that takes an integer template parameter specifying the codimension
  457 over which to iterate. Within the \lstinline!Codim! structure the type
  458 \lstinline!Iterator! is defined. Since we specified codimension 0
  459 this iterator is used to iterate
  460 over the elements which are not refined any further, i.e.~which are
  461 the leaves of the refinement trees.
  463 The \lstinline!for!-loop in line \ref{tc:forel} now visits every such
  464 element. The \lstinline!begin! and \lstinline!end! on the
  465 \lstinline!LeafGridView!
  466 class deliver the first leaf element and one past the last leaf element. Note
  467 that the \lstinline!template! keyword must be used and template parameters are
  468 passed explicitly. Within the loop body in
  469 lines \ref{tc:forel0}-\ref{tc:forel1} the iterator \lstinline!it! acts
  470 like a pointer to an entity of dimension \lstinline!dim! and
  471 codimension 0. The exact type would be
  472 \lstinline!typename G::template Codim<0>::Entity! just to mention it.
  474 Please note, that from Dune 2.4 on, C++11 range based for statements will
  475 be used to iterate over entities of a grid. To enable this feature, you will
  476 need \texttt{g++} version 4.6 or higher. There are no issues with the mentioned
  477 version of \texttt{icc} or \texttt{clang}. The entire iteration loop,
  478 will then boil down to
  479 \begin{lstlisting}[basicstyle=\ttfamily\scriptsize,numbers=left,
  480 numberstyle=\tiny, numbersep=5pt]
  481  for (auto&& e : elements(leafView))
  482  {
  483    // the iteration loop
  484  }
  485 \end{lstlisting}
  487 There is no more need of getting the iterator type by hand. Note that you
  488 should always write \lstinline!auto&&! !In contrast to the above procedure you don't
  489 get an \lstinline!EntityPointer!, but a reference to the \lstinline!Entity!
  490 itself.
  492 An important part of an entity is its geometrical shape and position.
  493 All geometrical information is factored out into a sub-object that can
  494 be accessed via the \lstinline!geometry()!  method. The geometry
  495 object is in general a mapping from a $d$-dimensional polyhedral
  496 reference element to $w$ dimensional space. Here we have $d=$
  497 \lstinline!G::dimension! and $w=$ \lstinline!G::dimensionworld!. This
  498 mapping is also called the ``local to global'' mapping.  The
  499 corresponding reference element has a certain type which is extracted
  500 in line \ref{tc:reftype}. Since the reference elements are polyhedra
  501 they consist of a finite number of corners. The images of the corners
  502 under the local to global map can be accessed via the
  503 \lstinline!corner(int n)! method. Line \ref{tc:print} prints the geometry type
  504 and the position of the first corner of the element. Then line
  505 \ref{tc:count} just counts the number of elements visited.
  507 Suppose now that we wanted to iterate over the vertices of the leaf
  508 grid instead of the elements. Now vertices have the codimension
  509 \lstinline!dim! in a \lstinline!dim!-dimensional grid and a
  510 corresponding iterator is provided by each grid class. It is extracted
  511 in line \ref{tc:vertit} for later use. The \lstinline!for!-loop
  512 starting in line \ref{tc:forve} is very similar to the first one
  513 except that it now uses the \lstinline!VertexLeafIterator!.  As you
  514 can see the different entities can be accessed with the same methods.
  515 We will see later that codimensions 0 and \lstinline!dim! are
  516 specializations with an extended interface compared to all other
  517 codimensions. You can also access the codimensions between 0 and
  518 \lstinline!dim!. However, currently not all implementations of the
  519 grid interface support these intermediate codimensions (though this
  520 does not restrict the implementation of finite element methods with
  521 degrees of freedom associated to, say, faces).
  523 Again, above task, can be done with a range based for statement by
  524 using the function \lstinline!vertices! instead of \lstinline!elements!
  525 in the above code snippet.
  527 Finally, we show in lines \ref{tc:level0}-\ref{tc:level1} how the
  528 hierarchic structure of the mesh can be accessed. To that end a
  529 \lstinline!LevelGridView! is used. It provides via an \lstinline!Iterator!
  530 access to all entities of a given codimension (here 0) on a given grid level.
  531 The coarsest
  532 grid level (the initial macro grid) has number zero and the number of
  533 the finest grid level is returned by the \lstinline!maxLevel()! method
  534 of the grid.  The methods \lstinline!begin()! and \lstinline!end()!
  535 on the view deliver iterators to the first and one-past-the-last
  536 entity of a given grid level supplied as an integer argument to these
  537 methods.
  539 The following listing shows the output of the program.
  541 \begin{lst}[Output of traversal] \mbox{}
  543 \begin{lstlisting}[basicstyle=\ttfamily\scriptsize]
  544 *** Traverse codim 0 leaves
  545 visiting leaf (cube, 2) with first vertex at -1 -1
  546 visiting leaf (cube, 2) with first vertex at 0 -1
  547 visiting leaf (cube, 2) with first vertex at -1 0
  548 visiting leaf (cube, 2) with first vertex at 0 0
  549 there are/is 4 leaf element(s)
  551 *** Traverse codim 2 leaves
  552 visiting (cube, 0) at -1 -1
  553 visiting (cube, 0) at 0 -1
  554 visiting (cube, 0) at 1 -1
  555 visiting (cube, 0) at -1 0
  556 visiting (cube, 0) at 0 0
  557 visiting (cube, 0) at 1 0
  558 visiting (cube, 0) at -1 1
  559 visiting (cube, 0) at 0 1
  560 visiting (cube, 0) at 1 1
  561 there are/is 9 leaf vertices(s)
  563 *** Traverse codim 0 level-wise
  564 visiting (cube, 2) with first vertex at -1 -1
  565 there are/is 1 element(s) on level 0
  567 visiting (cube, 2) with first vertex at -1 -1
  568 visiting (cube, 2) with first vertex at 0 -1
  569 visiting (cube, 2) with first vertex at -1 0
  570 visiting (cube, 2) with first vertex at 0 0
  571 there are/is 4 element(s) on level 1
  572 \end{lstlisting}
  573 \end{lst}
  575 \begin{rem} Define the end iterator for efficiency.
  576 \end{rem}
  578 \begin{exc} Play with different dimensions, codimension
  579   (\lstinline!YaspGrid! supports all codimenions) and refinements.
  580 \end{exc}
  582 \begin{exc} The method \lstinline!corners()! of the geometry returns
  583   the number of corners of an entity. Modify the code such that the
  584   positions of all corners are printed.
  585 \end{exc}
  588 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  589 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  590 \chapter{The \texorpdfstring{\Dune{}}{DUNE} grid interface}
  591 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  592 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  595 \section{Grid definition}
  597 There is a great variety of grids: conforming and non-conforming
  598 grids, single-element-type and multiple-element-type grids, locally
  599 and globally refined grids, nested and non-nested grids,
  600 bisection-type grids, red-green-type grids, sparse grids and so on. In
  601 this section we describe in some detail the type of grids that are
  602 covered by the \Dune{} grid interface.
  604 \minisec{Reference elements}
  606 A computational grid is a non-overlapping subdivision of a domain
  607 $\Omega\subset\R^w$ into elements of ``simple'' shape. Here ``simple''
  608 means that the element can be represented as the image of a reference
  609 element under a transformation. A reference element is a convex
  610 polytope, which is a bounded intersection of a finite set of
  611 half-spaces.
  613 \minisec{Dimension and world dimension}
  615 A grid has a dimension $d$ which is the dimensionality of
  616 its reference elements. Clearly we have $d\leq w$. In the case $d<w$ the grid
  617 discretizes a $d$-dimensional manifold.
  619 \minisec{Faces, entities and codimension}
  621 The intersection of a $d$-dimensional convex polytope (in
  622 $d$-dimensional space) with a
  623 tangent plane is called a face (note that there are faces of
  624 dimensionality $0,\ldots,d-1$). Consequently, a face of a grid element
  625 is defined as the image of a face of its reference element under the
  626 transformation. The elements and faces of elements of a grid are
  627 called its entities. An entity is said to be of codimension $c$ if it
  628 is a $d-c$-dimensional object. Thus the elements of the grid are
  629 entities of codimension 0, facets of an element have codimension 1,
  630 edges have codimension $d-1$ and vertices have codimension $d$.
  632 \minisec{Conformity}
  634 Computational grids come in a variety of flavours: A
  635 {conforming} grid is one where the intersection of two
  636 elements is either empty or a face of each of the two elements.
  637 Grids where the intersection of two elements may have an
  638 arbitrary shape are called {nonconforming}.
  640 \minisec{Element types}
  642 A {simplicial} grid is one where the reference elements are
  643 simplices. In a {multi-element-type} grid a finite number of
  644 different reference elements are allowed. The \Dune{} grid interface
  645 can represent conforming as well as non-conforming grids.
  647 \minisec{Hierarchically nested grids, macro grid}
  649 A {hierarchically nested} grid consists of a collection of $J+1$
  650 grids that are subdivisions of nested domains $$\Omega=\Omega_0 \supseteq \Omega_1 \supseteq
  651 \ldots \supseteq \Omega_J.$$ Note that only $\Omega_0$ is required to
  652 be identical to $\Omega$. If $\Omega_0=\Omega_1=\ldots=\Omega_J$ the
  653 grid is {globally refined}, otherwise it is {locally refined}.
  654 The grid that discretizes $\Omega_0$ is called the macro grid and its
  655 elements the macro elements. The
  656 grid for $\Omega_{l+1}$ is obtained from the grid
  657 for $\Omega_l$ by possibly subdividing each of its elements into
  658 smaller elements. Thus, each element of the macro grid and the
  659 elements that are obtained from refining it form a tree structure. The
  660 grid discretizing $\Omega_l$ with $0\leq l \leq J$ is called the level-$l$-grid and its
  661 elements are obtained from an $l$-fold refinement of some macro elements.
  663 \minisec{Leaf grid}
  665 Due to the nestedness of the domains we can partition the domain
  666 $\Omega$ into $$\Omega = \Omega_J \cup \bigcup_{l=0}^{J-1}
  667 \Omega_l\setminus\Omega_{l+1}.$$ As a consequence of the hierarchical
  668 construction a computational grid discretizing $\Omega$ can be
  669 obtained by taking the elements of the level-$J$-grid plus
  670 the elements of the level-$J-1$-grid in the region
  671 $\Omega_{J-1}\setminus\Omega_{J}$ plus the elements of the level-$J-2$-grid in the region
  672 $\Omega_{J-2}\setminus\Omega_{J-1}$ and so on plus the elements of the level-$0$-grid in the region
  673 $\Omega_{0}\setminus\Omega_{1}$. The grid resulting from this
  674 procedure is called the leaf grid
  675 because it is formed by the leaf elements of the trees emanating at
  676 the macro elements.
  678 \minisec{Refinement rules}
  680 There is a variety of ways how to hierarchically refine a grid. The
  681 refinement is called conforming if the leaf grid is always a
  682 conforming grid, otherwise the refinement is called
  683 non-conforming. Note that the grid on each level $l$ might be
  684 conforming while the leaf grid is not.
  685 There are also many ways how to subdivide an individual element into
  686 smaller elements. Bisection always subdivides elements into two
  687 smaller elements, thus the resulting data structure is a binary tree
  688 (independent of the dimension of the grid). Bisection is sometimes
  689 called ``green'' refinement. The so-called ``red'' refinement is the
  690 subdivision of an element into $2^d$ smaller elements, which is most
  691 obvious for cube elements. In many practical situation anisotropic
  692 refinement, i.~e.~refinement in a preferred direction, may be required.
  694 \minisec{Summary}
  696 The \Dune{} grid interface is able to represent grids with the
  697 following properties:
  698 \begin{itemize}
  699 \item Arbitrary dimension.
  700 \item Entities of all codimensions.
  701 \item Any kind of reference elements (you could define the icosahedron
  702   as a reference element if you wish).
  703 \item Conforming and non-conforming grids.
  704 \item Grids are always hierarchically nested.
  705 \item Any type of refinement rules.
  706 \item Conforming and non-conforming refinement.
  707 \item Parallel, distributed grids.
  708 \end{itemize}
  712 \section{Concepts}
  714 Generic algorithms are based on concepts. A concept is a kind of
  715 ``generalized'' class with a well defined set of members.
  716 Imagine a function template that takes a type \lstinline!T!
  717 as template argument. All the members of \lstinline!T!, i.e.
  718 methods, enumerations, data (rarely) and nested
  719 classes  used by the function template form the concept. From that
  720 definition it is clear that the concept does not necessarily exist as
  721 program text.
  723 A class that implements a concept is called a
  724 \textit{model} of the concept. E.\,g.~in the standard template library (STL)
  725 the class \lstinline!std::vector<int>! is a model of the concept
  726 ``container''. If all instances of a class template are a model of
  727 a given concept we can also say that the class template is a model of
  728 the concept. In that sense \lstinline!std::vector! is also a model of
  729 container.
  731 In standard OO language a concept would be formulated as
  732 an abstract base class and all the models would be implemented as
  733 derived classes. However, for reasons of efficiency we do not want to
  734 use dynamic polymorphism. Moreover, concepts are more powerful because
  735 the models of a concept can use different types, e.\,g.~as return types of
  736 methods. As an example consider the STL where the begin method on a
  737 vector of int returns \lstinline!std::vector<int>::iterator! and on a
  738 list of int it returns \lstinline!std::list<int>::iterator! which may
  739 be completely different types.
  741 Concepts are difficult to describe when they do not exist as concrete
  742 entities (classes or class templates) in a program. The STL way of
  743 specifying concepts is to describe the members \lstinline!X::foo()! of
  744 some arbitrary model named \lstinline!X!. Since this decription of the
  745 concept is not processed by the compiler it can get inconsistent and
  746 there is no way to check conformity of a model to the interface. As a
  747 consequence, strange error messages from the compiler may be the
  748 result (well C++ compilers can always produce strange error messages).
  749 There are two ways to improve the situation:
  750 \begin{itemize}
  751 \item \textit{Engines:} A class template is defined that wraps the
  752   model (which is the template parameter) and forwards all member
  753   function calls to it. In addition all the nested types and
  754   enumerations of the model are copied into the wrapper class.
  755   The model can be seen as an engine that powers the wrapper class,
  756   hence the name. Generic
  757   algorithms are written in terms of the wrapper class. Thus the
  758   wrapper class encapsulates the concept and it can be ensured
  759   formally by the compiler that
  760   all members of the concept are implemented.
  762 \item \textit{Barton-Nackman trick:} This is a refinement of the
  763   engine approach where the models are derived from the wrapper class
  764   template in addition. Thus static polymorphism is combined
  765   with a traditional class hierarchy, see \cite{Veldhui99,BN}.
  766   However, the
  767   Barton-Nackman trick gets rather involved when the derived classes
  768   depend on additional template parameters and several types are related
  769   with each other. That is why it is not used at all places in \Dune.
  770 \end{itemize}
  773 The \Dune{} grid interface now consists of a \textit{set of related concepts}.
  774 Either the engine or the Barton-Nackman approach are used to clearly
  775 define the concepts. In order to avoid any inconsistencies we refer as
  776 much as possible to the doxygen-generated documentation. For an
  777 overview of the grid interface see the web page
  778 \begin{center}
  779 \url{http://www.dune-project.org/doxygen/master/group__Grid.html}.
  780 \end{center}
  783 \subsection{Common types}
  785 Some types in the grid interface do not depend on a specific model,
  786 i.~e.~they are shared by all implementations.
  788 \minisec{\href{https://www.dune-project.org/doxygen/master/group__GeometryReferenceElements.html}{Dune::ReferenceElement}}
  790 describes the topology and geometry of standard entities. Any given
  791 entity of the grid can be completely specified by a reference element
  792 and a map from this reference element to world coordinate space.
  794 \minisec{\href{http://www.dune-project.org/doxygen/master/classDune_1_1GeometryType.html}{Dune::GeometryType}}
  796 defines names for the reference elements.
  798 \minisec{\href{http://www.dune-project.org/doxygen/master/classDune_1_1CollectiveCommunication.html}{Dune::CollectiveCommunication}}
  800 defines an interface to global communication operations in a portable
  801 and transparent way. In particular also for sequential grids.
  804 \subsection{Concepts of the \texorpdfstring{\Dune{}}{DUNE} grid interface}
  806 In the following a short description of each concept in the \Dune{}
  807 grid interface is given. For the details click on the link that leads
  808 you to the documentation of the corresponding wrapper class template
  809 (in the engine sense).
  811 \minisec{\href{http://www.dune-project.org/doxygen/master/classDune_1_1Grid.html}{Grid}}
  813     The grid
  814     is a container of entities that allows to access these entities
  815     and that knows the number of its entities. You create instances of
  816     a grid class in your applications, while objects of the other
  817     classes are typically aggregated in the grid class and accessed via
  818     iterators.
  820 \minisec{\href{http://www.dune-project.org/doxygen/master/classDune_1_1GridView.html}{GridView}}
  822     The GridView gives a view on a level or the leaf part of a grid. It
  823     provides iterators for access to the elements of this view and a
  824     communication method for parallel computations. Alternatively, a
  825     LevelIterator of a LeafIterator  can be directly accessed from a grid.
  826     These iterator types are described below.
  828 \minisec{\href{http://www.dune-project.org/doxygen/master/classDune_1_1Entity.html}{Entity}}
  830     The entity class encapsulates
  831     the topological part of an entity, i.\,e. its hierarchical
  832     construction from subentities and the relation to other
  833     entities. Entities cannot be created, copied or modified by the
  834     user. They can only be read-accessed through immutable iterators.
  836 \minisec{\href{http://www.dune-project.org/doxygen/master/classDune_1_1Geometry.html}{Geometry}}
  838     Geometry
  839     encapsulates the geometric part of an entity by mapping local
  840     coordinates in a reference element to world coordinates.
  842 \minisec{\href{http://www.dune-project.org/doxygen/master/classDune_1_1EntityPointer.html}{EntityPointer}}
  844     EntityPointer  is a
  845     dereferenceable type that delivers a reference to an
  846     entity. Moreover it is immutable, i.\,e. the referenced entity can
  847     not be modified.
  849 \minisec{\href{http://www.dune-project.org/doxygen/master/classDune_1_1EntityIterator.html}{EntityIterator}}
  851     EntityIterator is an
  852     immutable iterator that provides access to an entity. It can by
  853     incremented to visit all entities of a given
  854     codimension of a GridView. An EntityPointer is assignable
  855     from an Iterator.
  857 \minisec{\href{http://www.dune-project.org/doxygen/master/classDune_1_1IntersectionIterator.html}{IntersectionIterator}}
  859     IntersectionIterator provides access to all entities of
  860     codimension 0 that have an intersection of codimension 1 with a
  861     given entity of codimension 0. In a conforming mesh these are the
  862     face neighbors of an element. For two entities with a common
  863     intersection the IntersectionIterator can be dereferenced as an Intersection
  864     object which in turn provides information about the geometric location of
  865     the intersection.  Furthermore this Intersection class also provides
  866     information about intersections of an entity
  867     with
  868     the internal or external boundaries.
  869     The IntersectionIterator provides intersections between
  870     codimension 0 entities among the same GridView.
  872 \minisec{\href{http://www.dune-project.org/doxygen/master/classDune_1_1IndexSet.html}{LevelIndexSet}, \href{http://www.dune-project.org/doxygen/master/classDune_1_1IndexSet.html}{LeafIndexSet}}
  874     LevelIndexSet and LeafIndexSet, which are both models of
  875     Dune::IndexSet, are used to attach any kind of user-defined data to
  876     (subsets of) entities of the grid. This data is supposed to be
  877     stored in one-dimensional arrays for reasons of efficiency. An IndexSet is
  878     usually not used directly but through a Mapper (c.f.~chapter
  879     \ref{ch:mappers}).
  881 \minisec{\href{http://www.dune-project.org/doxygen/master/classDune_1_1IdSet.html}{LocalIdSet}, \href{http://www.dune-project.org/doxygen/master/classDune_1_1IdSet.html}{GlobalIdSet}}
  883     LocalIdSet and GlobalIdSet which are both models of Dune::IdSet
  884     are used to save user data during a grid refinement phase and
  885     during dynamic load balancing in the parallel case. The LocalIdSet is
  886     unique for all entities on the current partition, whereas the GlobalIdSet
  887     gives a unique mapping over all grid partitions. An IdSet is usually not used
  888     directly but through a Mapper (c.f.~chapter
  889     \ref{ch:mappers}).
  892 \section{Propagation of type information}
  894 The types making up one grid implementation cannot be mixed with the
  895 types making up another grid implementation. Say, we have two
  896 implementations of the grid interface \lstinline!XGrid! and
  897 \lstinline!YGrid!. Each implementation provides a LevelIterator
  898 class, named \lstinline!XLevelIterator! and
  899 \lstinline!YLevelIterator! (in fact, these are class templates because
  900 they are parametrized by the codimension and other
  901 parameters). Although these types implement the same interface they
  902 are distinct classes that are not related in any way for the
  903 compiler. As in the Standard Template Library strange error messages
  904 may occur if you try to mix these types.
  906 In order to avoid these problems the related types of an
  907 implementation are provided as public types by most classes of an
  908 implementation. E.\,g., in order to extract the
  909 \lstinline!XLevelIterator! (for codimension 0) from the
  910 \lstinline!XGrid! class you would write
  911 \begin{lstlisting}
  912 XGrid::template Codim<0>::LevelIterator
  913 \end{lstlisting}
  914 Because most of the types are parametrized by certain parameters like
  915 dimension, codimension or partition type simple typedefs (as in the
  916 STL) are not sufficient here. The types are rather placed in a
  917 struct template, named \lstinline!Codim! here, where the template
  918 parameters of the struct are those of the type. This concept may even
  919 be applied recursively.
  921 \chapter{Constructing grid objects}
  923 So far we have talked about the grid interface and how you can access and
  924 manipulate grids.  This chapter will show you how you create grids in the
  925 first place.  There are several ways to do this.
  927 The central idea of \Dune is that all grid implementations behave equally
  928 and conform to the same interface.  However, this concept fails when it
  929 comes to constructing grid objects, because grid implementations differ too
  930 much to make one construction method work for all.  For example, for an
  931 unstructured grid you have to specify all vertex positions, whereas for a
  932 structured grid this would be a waste of time.  On the other hand, for
  933 a structured grid you may need to give the bounding box which, for an
  934 unstructured grid, is not necessary.  In practice, these differences do
  935 not pose serious problems.
  937 In this chapter, creating a grid always means creating a grid with only
  938 a single level.  Such grid is alternatively called a {\em coarse grid}
  939 or a {\em macro grid}.  There is currently no functionality in \Dune to
  940 set up hierarchical grids directly.  The underlying assumption is that
  941 the user will create a coarse grid first and then generate a hierarchy
  942 using refinement.  Despite the name (and grid implementations permitting),
  943 the coarse grid can of course be as large and fine as desired.
  945 \section{Creating Structured Grids}
  947 Creating structured grids is comparatively simple, as little information
  948 needs to be provided.  In general, for uniform structured grids, the grid
  949 dimension, bounding box, and number of elements in each direction suffices.
  950 Such information can be given directly with the constructor of the grid
  951 object. \Dune does not currently specify the signature of grid constructors,
  952 and hence they are all slightly different for different grid implementations.
  954 As already seen int he introduction, the code to construct a sequential \
  955 \lstinline!YaspGrid! with 10 cells in each direction is
  956 \begin{lstlisting}
  957     std::array<int,dim> n;
  958     std::fill(n.begin(). n.end(), 10);
  960     Dune::FieldVector<double,dim> upper(1.0);
  962     YaspGrid<dim> grid(upper, n);
  963 \end{lstlisting}
  964 Note that for an equidistant grid, you do not have to specify the lower left
  965 corner as \lstinline!YaspGrid! hardwires it to zero. This limitation can be
  966 overcome by using another feature of \lstinline!YaspGrid!:
  968 \lstinline!YaspGrid! can also be used a tensorproduct grid. A tensorproduct grid
  969 is defined by a set of coordinates $\{ x_{i,0}, \dots ,x_{i,n_i}\}$ for each
  970 direction $i$. The vertices of the grid are then given by the tuples
  971 \begin{displaymath}
  972  \{
  973   (x_{0,i_0},\dots ,x_{d-1,i_{d-1}})\ \
  974    with\ \  0\leq i_j\leq n_j \ \ \forall j\}
  975 \end{displaymath}
  976 Such tensorproduct grid is a structured non-equidistant grid. Having a non-zero
  977 lower left corner is possible by specifying coordinates accordingly. \lstinline!YaspGrid!s
  978 tensorproduct features are enabled by using \lstinline!TensorProductCoordinates<ctype,dim>!
  979 as the second template paramter. Coordinates are given as as a \lstinline!std::array<std::vector<ctype>,dim>!.
  980 See the following example which initializes a grid is one cell, that covers $[-1,1]^2$:
  982 \begin{lstlisting}
  983  std::array<std::vector<ctype,2> > coords;
  984  coords[0] = {-1., 1.};
  985  coords[1] = {-1., 1.};
  987  YaspGrid<dim, TensorProductCoordinates<ctype,dim> > grid(coords);
  988 \end{lstlisting}
  990 Another unstructured is the one-dimensional \lstinline!OneDGrid!, which has a constructor
  991 \begin{lstlisting}
  992     OneDGrid grid(10,    // number of elements
  993                   0.0,   // left domain boundary
  994                   1.0    // right domain boundary
  995                   );
  996 \end{lstlisting}
  997 for uniform grids.
  999 \section{Reading Unstructured Grids from Files}
 1001 Unstructured grids usually require much more information than what can
 1002 reasonable be provided within the program code.  Instead, they are usually
 1003 read from special files.  A large variety of different file formats for
 1004 finite element grids exists, and \Dune provides support for some of them.
 1005 Again, no interface specification exists for file readers in \Dune.
 1007 Arguably the most important file format currently supported by \Dune is
 1008 the \lstinline!gmsh! format.  \lstinline!Gmsh!\footnote{\url{http://geuz.org/gmsh/}}
 1009 is an open-source geometry modeler and grid generator.  It allows to define
 1010 geometries using a boundary representation (interactively and via its
 1011 own modeling language resulting in \lstinline!.geo!-files),
 1012 creates simplicial grids in 2d and 3d (using
 1013 \lstinline!tetgen! or \lstinline!netgen!)
 1014 and stores them in files ending in \lstinline!.msh!. Current precompiled
 1015 releases $\geq$ 2.4.2 of \lstinline!Gmsh! have \lstinline!OpenCascade!, an open-source CAD
 1016 kernel, as built-in geometry modeler. Thus these releases are able to import
 1017 CAD geometries, e.\,g.~from IGES or STEP files, and to generate meshes for them to
 1018 be subsequently used in \Dune. Be aware that most versions of Gmsh available in
 1019 the package repositories of your operating system still lack this functionality.
 1020 Further, Gmsh and the Gmsh reader of \Dune support second order
 1021 boundary segments thus providing a rudimentary support for curved boundaries.
 1022 To read such a file into a FooGrid, include \lstinline!dune/grid/io/file/gmshreader.hh!
 1023 and write
 1024 \begin{lstlisting}
 1025     FooGrid* grid = GmshReader<GridType>::read(filename);
 1026 \end{lstlisting}
 1028 A second format is AmiraMesh, which is the native format of the Amira.\footnote{\url{http://www.amira.com/}}
 1029 To read AmiraMesh files you need to have the library libamiramesh%
 1030 \footnote{\url{http://amira.com/downloads/patch-archive/patches412/81.html}}
 1031 installed.  Then
 1032 \begin{lstlisting}
 1033     FooGrid* grid = AmiraMeshReader<GridType>::read(filename);
 1034 \end{lstlisting}
 1035 reads the grid in \lstinline!filename! into the FooGrid.
 1037 Further available formats are StarCD and the native Alberta format.
 1038 See the class documentation of dune-grid for an up-to-date list.
 1039 Demo grids for each format can be found in \lstinline!dune-grid/doc/grids!.
 1040 They exist for documentation and also as example grids for the unit
 1041 tests of the file readers.  The unit tests should not hardwire the path
 1042 to the example grids.  Instead, the path should be provided in the preprocessor
 1043 variable \lstinline!DUNE_GRID_EXAMPLE_GRIDS_PATH!.
 1046 \section{The \texorpdfstring{\Dune}{Dune} Grid Format (DGF)}
 1048 Dune has its own macro grid format, the \underline{D}une \underline{G}rid \underline{F}ormat.
 1049 A detailed description of the DGF and how to use it can be found on the
 1050 homepage of \Dune%
 1051 \footnote{\url{http://www.dune-project.org/doxygen/master/classDune_1_1DuneGridFormatParser.html}}.
 1053 Here we only give a short introduction.
 1054 Basically one can choose the grid manager during the make process of
 1055 your program:
 1056 \lstinline!make GRIDTYPE=MYGRID GRIDDIM=d GRIDWORLD=w!
 1057 Including \lstinline!config.h! will then
 1058 introduce the namespace \lstinline!GridSelector! into the
 1059 \lstinline!Dune! namespace. This contains the typedef
 1060 \lstinline!GridType! which is the type of the grid. Furthermore the
 1061 required header files for the grid manager are included.
 1062 Through the module \Dune{}-grid the following grid managers can be used
 1063 (for \lstinline!MYGROD! in the example above):
 1064 \\
 1066 \\
 1067 and \lstinline!YASPGRID!.
 1068 \\
 1069 The following example shows how an
 1070 instance of the defined grid is generated. Given a DGF file, for example
 1071 \lstinline!unitcube2.dgf!, a grid pointer is created as follows
 1072 \begin{lstlisting}[basicstyle=\ttfamily\scriptsize]
 1073 Dune::GridPtr<Dune::GridSelector::GridType> gridPtr( "unitcube2.dgf" );
 1074 \end{lstlisting}
 1075 The grid is accessed by dereferencing the grid pointer.
 1076 \begin{lstlisting}[basicstyle=\ttfamily\scriptsize]
 1077 GridType& grid = *gridPtr;
 1078 \end{lstlisting}
 1079 To change the grid one simply has to re-compile the code using the following make command.
 1080 \begin{lstlisting}[basicstyle=\ttfamily\scriptsize]
 1081 make GRIDDIM=2 GRIDTYPE=ALBERTAGRID integration
 1082 \end{lstlisting}
 1083 This will compile the application \texttt{integration} with grid type \lstinline!ALBERTAGRID! and grid dimension $2$.
 1084 Note that before the re-compilation works,
 1085 the corresponding object file has to be removed.
 1086 If \lstinline!WORLDDIM! is not
 1087 provide then \lstinline!WORLDDIM=GRIDDIM! is assumed.
 1088 To use some grid manager by default, i.e., without providing the grid type
 1089 during the make process, \lstinline!GRIDTYPE! and
 1090 \lstinline!GRIDDIM,WORLDDIM! can be set in \lstinline!Makefile.am!. It is
 1091 then still possible to change the default during the make process.
 1093 \section{The Grid Factory}
 1094 \label{sec:grid_factory}
 1096 While there is currently no convention on what a file reader should look like,
 1097 there is a formally specified low-level interface for the construction of
 1098 unstructured coarse grids.  This interface, which goes by the name of
 1099 GridFactory, provides methods to, e.\,g.~insert vertices and elements one by one.
 1100 It is the basis of the file readers described in the previous section.
 1101 The main reason why you may want to program the GridFactory directly is
 1102 when writing your own grid readers.  However, in some cases it may also be
 1103 most convenient to be able to specify your coarse grid entirely in the
 1104 C++ code.  You can do that using the GridFactory.
 1106 The GridFactory is programmed as a factory class (hence the name).
 1107 You default-construct an object of the factory class, provide it with all
 1108 necessary information, and it will create and hand over a grid for you.
 1109 In the following we will describe the use of the GridFactory in more detail.
 1110 Say you are interested in creating a new grid of type \lstinline!FooGrid!.  Then
 1111 you proceed as follows:
 1113 \begin{enumerate}
 1114 \item Create a GridFactory Object
 1116 Get a new GridFactory object by calling
 1118 \begin{lstlisting}[basicstyle=\small]
 1119     GridFactory< FooGrid > factory;
 1120 \end{lstlisting}
 1122 \item Enter the Vertices
 1124 Insert the grid vertices by calling
 1126 \begin{lstlisting}[basicstyle=\small]
 1127     factory.insertVertex(const FieldVector<ctype,dimworld>& position);
 1128 \end{lstlisting}
 1130 for each vertex. The order of insertion determines the level- and leaf indices of your level 0 vertices.
 1132 \item Enter the elements
 1134 For each element call
 1136 \begin{lstlisting}[basicstyle=\small]
 1137     factory.insertElement(Dune::GeometryType type,
 1138         const std::vector<int>& vertices);
 1139 \end{lstlisting}
 1141 The parameters are
 1143 \begin{itemize}
 1144     \item type - The element type. The grid implementation is expected to throw
 1145     an exception if an element type that cannot be handled is encountered.
 1146     \item vertices - The indices of the vertices of this element.
 1147 \end{itemize}
 1149 The numbering of the vertices of each element is expected to follow the \Dune
 1150 conventions. Refer to the page on reference elements for the details.
 1152 \item Parametrized Domains
 1154 \lstinline!FooGrid! may support parametrized domains. That means that you can provide a smooth description of your grid boundary. The actual grid may always be piecewise linear; however, as you refine, the grid will approach your prescribed boundary. You don't have to do this. If you do not specify the boundary geometry it is left to the grid implementation.
 1156 In order to create curved boundary segments, for each segment you have to write a class which implements the correct geometry. These classes are then handed over to the
 1157 factory. Boundary segment implementations must be derived from
 1159 \begin{lstlisting}[basicstyle=\small]
 1160     template <int dim, int dimworld> Dune::BoundarySegment
 1161 \end{lstlisting}
 1163 This is an abstract base class which requires you to overload the method
 1165 \begin{lstlisting}[basicstyle=\small]
 1166     virtual FieldVector< double, dimworld >
 1167          operator() (const FieldVector< double, dim-1 > &local)
 1168 \end{lstlisting}
 1170 This methods must compute the world coordinates from the local ones on the boundary segment. Give these classes to your grid factory by calling
 1172 \begin{lstlisting}[basicstyle=\small]
 1173     factory.insertBoundarySegment(const std::vector<int>& vertices,
 1174                                   const BoundarySegment<dim,dimworld> *
 1175                                      boundarySegment = NULL);
 1176 \end{lstlisting}
 1178 Control over the allocated objects is taken from you, and the grid object will take care of their destruction.
 1180 Note that you can call \lstinline!insertBoundarySegment! with only the first argument.
 1181 In that case, the boundary geometry is left to the grid implementation.  However,
 1182 the boundary segments get ordered in the way you inserted them.  This may be
 1183 helpful if you have data attached to your coarse grid boundary
 1184 (see Sec.~\ref{sec:import_data_for_new_grids}).
 1186 \item Finish construction
 1188 To finish off the construction of the FooGrid object call
 1190 \begin{lstlisting}[basicstyle=\small]
 1191     FooGrid* grid = factory.createGrid();
 1192 \end{lstlisting}
 1194 This time it is you who gets full responsibility for the allocated object.
 1195 \end{enumerate}
 1197 \section{Attaching Data to a New Grid}
 1198 \label{sec:import_data_for_new_grids}
 1200 In many cases there is data attached to new grids.  This data may be initial
 1201 values, spatially distributed material parameters, boundary conditions, etc.
 1202 It is associated to elements or vertices, or the boundary segments of the
 1203 coarse grid.  The data may be available in a separate data file or even
 1204 included in the same file with the grid.
 1206  The connection with the grid in
 1207 the grid file is usually made implicitly.  For example, vertex data is ordered
 1208 in the same order as the vertices itself.  Hence the grid-reading process must
 1209 guarantee that vertices and elements are not reordered during grid creation.
 1210 More specifically, \Dune guarantees the following: {\em the level and leaf
 1211 indices of zero-level vertices and elements are defined by the order in which they
 1212 were inserted into the grid factory}.  Note that this does not mean that
 1213 the vertices and elements are traversed in this order by the Level- and
 1214 LeafIterators.  What matters are the indices.  Note also that no such
 1215 promise is made concerning edges, faces and the like.  Hence it is currently
 1216 not possible to read edge and face data along with a grid without some
 1217 trickery.
 1219 It is also possible to attach data to boundary segments of the coarse grids.
 1220 For this, the method \lstinline!Intersection::boundaryId! (which should
 1221 really be called \lstinline!boundaryIndex!) returns an index when called
 1222 for a boundary intersection.  If the boundary intersection is on level zero
 1223 the index is consecutive and zero-starting.  For all other boundary intersections
 1224 it is the index of the zero-level ancestor boundary segment of the intersection.
 1226 If you have a list of data associated to certain boundary segments of your
 1227 coarse grid, you need some control on how the boundary ids are set.  Remember
 1228 from Sec.\,\ref{sec:grid_factory} that you can create a grid without mentioning
 1229 the boundary at all.  If you do that, the boundary ids are set automatically
 1230 by the grid implementation and the exact order is implementation-specific.
 1231 If you set boundary segments explicitly using the \lstinline!insertBoundarySegment!
 1232 method, then {\em the boundary segments are numbered in the order of their
 1233 insertion}.  If you do not set all boundary segments the remaining ones
 1234 get automatic, implementation-specific ids.  This is why the second argument
 1235 of \lstinline!insertBoundarySegment! is optional: you may want to
 1236 influence the ordering of the boundary segments, but leave the boundary
 1237 geometry to the grid implementation.  Calling \lstinline!insertBoundarySegment!
 1238 with a single argument allows you to do just this.
 1241 \section{Example: The \texorpdfstring{\lstinline{UnitCube}}{UnitCube} class}
 1243 In this chapter we give example code that shows how the different available grid
 1244 classes are instantiated. We create grids for the unit
 1245 cube $\Omega=(0,1)^d$ in various dimensions $d$.
 1247 Not all grid classes have the same interface for instantiation.
 1248 Unstructured grids are created using the \lstinline!GridFactory! class,
 1249 but for structured grids there is more variation.  In order make the
 1250 examples in later chapters easier to write we want to have a class template
 1251 \lstinline!UnitCube! that we parametrize with a type \lstinline!T! and
 1252 an integer parameter \lstinline!variant!. \lstinline!T! should be
 1253 one of the available grid types and \lstinline!variant! can be used to
 1254 generate different grids (e.\,g.~triangular or quadrilateral) for the
 1255 same type \lstinline!T!. The advantage of the \lstinline!UnitCube!
 1256 template is that the instantiation is hidden from the user.
 1258 The definition of the general template is as follows.
 1260 \begin{lst}[File dune-grid-howto/unitcube.hh] \mbox{}
 1261 \nopagebreak
 1262 \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
 1263 numberstyle=\tiny, numbersep=5pt, breaklines=true]{../unitcube.hh}
 1264 \end{lst}
 1266 This is a default implementation that uses the utility class
 1267 \lstinline!StructuredGridFactory! (from the header \lstinline!dune-grid/dune/grid/utility/structuredgridfactory.hh!)
 1268 to create grids for the unit cube.  The \lstinline!StructuredGridFactory!
 1269 uses the \lstinline!GridFactory! class (Section~\ref{sec:grid_factory})
 1270 internally to create structured simplicial and hexahedral grids.
 1271 Depending on the template parameter \lstinline!variant!, a
 1272 hexahedral (\lstinline!variant==1!) or simplicial (\lstinline!variant==2!)
 1273 grid is created.
 1275 The \lstinline!GridFactory! class is a required part of the grid interface
 1276 for all unstructured grids.  Hence the default implementation of \lstinline!UnitCube!
 1277 should work for all unstructured grids, namely \lstinline!UGGrid!,
 1278 \lstinline!OneDGrid!, \lstinline!ALUGrid!, and \lstinline!AlbertaGrid!.
 1279 The construction of structured grid objects is currently not standardized.
 1280 Therefore \lstinline!UnitCube! is specialized for each structured grid
 1281 type.  We now look at each specialization in turn.
 1283 For historic reasons, there are also specializations for
 1284 \lstinline!ALUGrid! and \lstinline!AlbertaGrid!.
 1286 \minisec{YaspGrid}
 1288 The following listing instantiates a \lstinline!YaspGrid! object. The
 1289 \lstinline!variant! parameter specifies the number of elements in each
 1290 direction of the cube. In the parallel case all available processes
 1291 are used and the overlap is set to one element. Periodicity is not
 1292 used.
 1294 \begin{lst}[File dune-grid-howto/unitcube\_yaspgrid.hh] \mbox{}
 1295 \nopagebreak
 1296 \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
 1297 numberstyle=\tiny, numbersep=5pt]{../unitcube_yaspgrid.hh}
 1298 \end{lst}
 1300 \minisec{AlbertaGrid}
 1302 The following listing contains specializations of the
 1303 \lstinline!UnitCube! template for Alberta in two and three
 1304 dimensions. When using Alberta versions less than $2.0$ the \Dune{} framework
 1305 has to be configured with a dimension (\lstinline!--with-alberta-dim=2!,
 1306 \lstinline!--with-alberta-world-dim=2!) and only this dimension can then be
 1307 used.
 1308 The dimension from the configure run is available in the macro
 1309 \lstinline!ALBERTA_DIM! and \lstinline!ALBERTA_WORLD_DIM!
 1310 in the file \lstinline!config.h! (see
 1311 next section). The \lstinline!variant! parameter must be 1.
 1312 The grid factory concept is used by the base class \lstinline!BasicUnitCube!.
 1314 \begin{lst}[File dune-grid-howto/unitcube\_albertagrid.hh] \mbox{}
 1315 \nopagebreak
 1316 \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
 1317 numberstyle=\tiny, numbersep=5pt]{../unitcube_albertagrid.hh}
 1318 \end{lst}
 1320 \minisec{ALUGrid}
 1322 The next listing shows the instantiation of \lstinline!ALUGrid!\ for
 1323 simplices and cubes.
 1324 The ALUGrid implementation supports either simplicial grids, i.e.\
 1325 tetrahedral or triangular grids, and hexahedral grids and the
 1326 element type has to be chosen at compile-time. This is done by choosing
 1327 either \lstinline!Dune::cube! or \lstinline!Dune::simplex! as the third
 1328 template argument for \lstinline!ALUGrid!.
 1329 The \lstinline!variant!\ parameter must be 1.
 1330 As in the default implementation, grid objects are set up with help of the
 1331 \lstinline!StructuredGridFactory!\ class.
 1333 \begin{lst}[File dune-grid-howto/unitcube\_alugrid.hh] \mbox{}
 1334 \nopagebreak
 1335 \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
 1336 numberstyle=\tiny, numbersep=5pt, breaklines=true]{../unitcube_alugrid.hh}
 1337 \end{lst}
 1341 %\section{Capabilities}
 1343 % \chapter{Reference elements}
 1346 \chapter{Quadrature rules}
 1347 \label{sec:quadrature}
 1349 In this chapter we explore how an integral $$\int\limits_{\Omega} f(x)\ dx$$
 1350 over some function $f:\Omega\to\mathbb{R}$ can be computed numerically
 1351 using a \Dune{} grid object.
 1353 \section{Numerical integration}
 1355 Assume first the simpler task that $\Delta$ is a reference element
 1356 and that we want to
 1357 compute the integral over some function $\hat{f}:\Delta\to\mathbb{R}$
 1358 over the reference element:$$\int\limits_{\Delta} \hat{f}(\hat{x})\ d\hat{x}.$$
 1361 A quadrature rule is a formula that approximates integrals of
 1362 functions over a reference element $\Delta$. In general it has the form
 1363 $$\int\limits_{\Delta} \hat{f}(\hat{x})\ d\hat{x} = \sum_{i=1}^n
 1364 \hat{f}(\xi_i) w_i + \text{error}.$$
 1365 The positions $\xi_i$ and weight factors $w_i$ are dependent on the
 1366 type of reference element and the number of quadrature points $n$ is
 1367 related to the error.
 1369 Using the transformation formula for integrals we can now compute
 1370 integrals over domains $\omega\subseteq\Omega$ that are mapped from a
 1371 reference element, i.~e.~$\omega=\{x\in\Omega\ |\
 1372 x=g(\hat{x}), \hat{x}\in\Delta\}$, by some function $g:\Delta\to\Omega$:
 1373 \begin{equation}
 1374 \int\limits_{\Omega} f(x) = \int\limits_{\Delta} f(g(\hat{x}))\mu(\hat{x})\
 1375 d\hat{x} = \sum_{i=1}^n f(g(\xi_i))\mu(\xi_i)w_i + \text{error}.
 1376 \label{Eq:IntegrateEntity}
 1377 \end{equation}
 1378 Here $\mu(\hat{x}) = \sqrt{|\det J^T(\hat{x})J(\hat{x})|}$ is the
 1379 integration element and $J(\hat{x})$ the Jacobian matrix of the map $g$.
 1381 The integral over the whole domain $\Omega$ requires a grid
 1382 $\overline{\Omega}=\bigcup_k \overline{\omega}_k$. Using
 1383 (\ref{Eq:IntegrateEntity}) on each element we obtain finally
 1384 \begin{equation}
 1385 \int\limits_{\Omega} f(x)\ dx = \sum\limits_{k} \sum_{i=1}^{n_k}
 1386 f(g^k(\xi^k_i))\mu^k(\xi^k_i)w^k_i + \sum\limits_{k} \text{error}^k.
 1387 \label{Eq:IntegrateDomain}
 1388 \end{equation}
 1389 Note that each element $\omega_k$ may in principle have its own
 1390 reference element which means that quadrature points and weights as
 1391 well as the transformation and integration element may depend on
 1392 $k$. The total error is a sum of the errors on the individual
 1393 elements.
 1395 In the following we show how the formula (\ref{Eq:IntegrateDomain})
 1396 can be realised within \Dune.
 1399 \section{Functors}
 1401 The function $f$ is represented as a functor, i.~e.~a class having an
 1402 \lstinline!operator()! with appropriate arguments. A point
 1403 $x\in\Omega$ is represented by an object of type
 1404 \lstinline!FieldVector<ct,dim>! where \lstinline!ct! is the type for
 1405 each component of the vector and \lstinline!dim! is its dimension.
 1408 \begin{lst}[dune-grid-howto/functors.hh] Here are some examples for functors.
 1410 \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
 1411 numberstyle=\tiny, numbersep=5pt]{../functors.hh}
 1412 \end{lst}
 1415 \section{Integration over a single element}
 1417 The function \lstinline!integrateentity! in the following listing
 1418 computes the integral over a single element of the mesh with a
 1419 quadrature rule of given order.
 1420 This relates directly to formula (\ref{Eq:IntegrateEntity}) above.
 1422 \begin{lst}[dune-grid-howto/integrateentity.hh] \mbox{}
 1424 \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
 1425 numberstyle=\tiny, numbersep=5pt]{../integrateentity.hh}
 1426 \end{lst}
 1428 Line \ref{ieh:qr} extracts a reference to a
 1429 \lstinline!Dune::QuadratureRule!  from the
 1430 \lstinline!Dune::QuadratureRules! singleton which is a container
 1431 containing quadrature rules for all the different reference element
 1432 types and different orders of approximation. A rule of order $p$ is
 1433 guaranteed to integrate polynomials of order $p$ exactly, up to
 1434 floating point errors.
 1436 Both classes are
 1437 parametrized by dimension and the basic type used for the coordinate
 1438 positions. \lstinline!Dune::QuadratureRule! in turn is a container of
 1439 \lstinline!Dune::QuadraturePoint! supplying positions $\xi_i$ and
 1440 weights $w_i$.
 1442 Line \ref{ieh:for} shows the loop over all quadrature points in the
 1443 quadrature rules. For each quadrature point $i$ the function value at
 1444 the transformed position (line \ref{ieh:fval}), the weight (line
 1445 \ref{ieh:weight}) and the integration element (line \ref{ieh:detjac})
 1446 are computed and summed (line \ref{ieh:result}).
 1448 \section{Integration with global error estimation}
 1450 In the listing below function \lstinline!uniformintegration!
 1451 computes the integral over the whole domain via formula
 1452 (\ref{Eq:IntegrateDomain}) and in addition provides an estimate of the
 1453 error. This is done as follows. Let $I_c$ be the value of the numerically
 1454 computed integral on some grid and let $I_f$ be the value of the
 1455 numerically computed integral on a grid where each element has been
 1456 refined. Then
 1457 \begin{equation}
 1458 \label{Eq:GlobalError}
 1459 E \approx |I_f-I_c|
 1460 \end{equation}
 1461 is an estimate for the error. If
 1462 the refinement is such that every element is bisected in every
 1463 coordinate direction, the function to be integrated is sufficiently
 1464 smooth and the order of the quadrature rule is $p$,
 1465 then the error should be reduced by a factor of $(1/2)^{p+1}$ after
 1466 each mesh refinement.
 1468 \begin{lst}[dune-grid-howto/integration.cc] \mbox{}
 1470 \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
 1471 numberstyle=\tiny, numbersep=5pt]{../integration.cc}
 1472 \end{lst}
 1474 Running the executable \lstinline!integration! on a YaspGrid in two
 1475 space dimensions with a quadrature rule of order
 1476 one the following output is obtained:
 1478 \begin{lstlisting}[basicstyle=\ttfamily\scriptsize]
 1479 elements=       1 integral=1.000000000000e+00 error=1.000000000000e+100
 1480 elements=       4 integral=6.674772311008e-01 error=3.325227688992e-01
 1481 elements=      16 integral=6.283027311366e-01 error=3.917449996419e-02
 1482 elements=      64 integral=6.192294777551e-01 error=9.073253381426e-03
 1483 elements=     256 integral=6.170056966109e-01 error=2.223781144285e-03
 1484 elements=    1024 integral=6.164524949226e-01 error=5.532016882082e-04
 1485 elements=    4096 integral=6.163143653145e-01 error=1.381296081435e-04
 1486 elements=   16384 integral=6.162798435779e-01 error=3.452173662133e-05
 1487 elements=   65536 integral=6.162712138101e-01 error=8.629767731416e-06
 1488 elements=  262144 integral=6.162690564098e-01 error=2.157400356695e-06
 1489 elements= 1048576 integral=6.162685170623e-01 error=5.393474630244e-07
 1490 elements= 4194304 integral=6.162683822257e-01 error=1.348366243104e-07
 1491 \end{lstlisting}
 1493 The ratio of the errors on two subsequent grids nicely approaches the
 1494 value $1/4$ as the grid is refined, which corresponds to an error
 1495 decay of order $p+1$.
 1498 \begin{exc} Try different quadrature orders. For that just change the
 1499   last argument of the call to \lstinline!integrateentity! in line
 1500   \ref{ic:call} in file \lstinline!integration.cc!.
 1501 \end{exc}
 1503 \begin{exc} Try different grid implementations and dimensions and
 1504   compare the run-time.
 1505 \end{exc}
 1507 \begin{exc} Try different integrands $f$ and look at the development
 1508   of the (estimated) error in the integral.
 1509 \end{exc}
 1511 \chapter{Attaching user data to a grid}
 1513 In most useful applications there will be the need to associate
 1514 user-defined data with certain entities of a grid. The standard
 1515 example are, of course, the degrees of freedom of a finite element
 1516 function. But it could be as simple as a boolean value that indicates
 1517 whether an entity has already been visited by some algorithm or
 1518 not. In this chapter we will show with some examples how arbitrary
 1519 user data can be attached to a grid.
 1521 \section{Mappers}\label{ch:mappers}
 1523 The general situation is that a user wants to store some arbitrary
 1524 data with a subset of the entities of a grid. Remember that entities
 1525 are all the vertices, edges, faces, elements, etc., on all the levels
 1526 of a grid.
 1528 An important design decision in the \Dune{} grid interface was that
 1529 user-defined data is stored in user space. This has a number of
 1530 implications:
 1531 \begin{itemize}
 1532 \item \Dune{} grid objects do not need to know anything about the user
 1533   data.
 1534 \item Data structures used in the implementation of a \Dune{} grid do
 1535   not have to be extensible.
 1536 \item Types representing the user data can be arbitrary.
 1537 \item The user is responsible for possibly reorganizing the data when
 1538   a grid is modified (i.~e.~refined, coarsened, load balanced).
 1539 \end{itemize}
 1541 Since efficiency is important in scientific computing the second
 1542 important design decision was that user data is stored in arrays (or
 1543 random access containers) and that the data is accessed via an
 1544 index. The set of indices starts at zero and is consecutive.
 1546 Let us assume that the set of all entities in the grid is $E$ and
 1547 that $E^\prime\subseteq E$ is the subset of entities for which data is
 1548 to be stored. E.\,g.~this could be all the vertices in the leaf grid in
 1549 the case of $P_1$ finite elements. Then the access from grid entities
 1550 to user data is a two stage process: A so-called \textit{mapper}
 1551 provides a map
 1552 \begin{equation}
 1553 m : E^\prime \to I_{E^\prime}
 1554 \end{equation}
 1555 where $I_{E^\prime}=\{0,\ldots,|E^\prime|-1\}\subset\mathbb{N}$ is the consecutive and
 1556 zero-starting index set associated to the entity set. The user data
 1557 $D(E^\prime)=\{d_e\ |\ e\in E^\prime\}$ is stored in an array, which
 1558 is another map
 1559 \begin{equation}
 1560 a : I_{E^\prime} \to D(E^\prime).
 1561 \end{equation}
 1562 In order to get the data $d_e\in D(E^\prime)$ associated to entity
 1563 $e\in E^\prime$ we therefore have to evaluate the two maps:
 1564 \begin{equation}
 1565 d_e = a(m(e)) .
 1566 \end{equation}
 1568 \Dune{} provides different implementations of mappers that differ in
 1569 functionality and cost (with respect to storage and
 1570 run-time). Basically there are two different kinds of mappers.
 1572 \minisec{Index based mappers}
 1574 An index-based mapper is allocated for a grid and can be used as long
 1575 as the grid is not changed (i.e. refined, coarsened or load
 1576 balanced). The implementation of these mappers is based on a
 1577 \lstinline!Dune::IndexSet! and evaluation of the map $m$ is typically
 1578 of $O(1)$ complexity with a very small constant.
 1579 Index-based mappers are only available for restricted (but
 1580 usually sufficient) entity sets. They will be used in the examples
 1581 shown below.
 1583 \minisec{Id based mappers}
 1585 Id-based mappers can also be used while a grid changes, i.e.~it is
 1586 ensured that the map $m$ can still be evaluated for all entities $e$
 1587 that are still in the grid after modification. For that it
 1588 has to be implemented on the basis of a \lstinline!Dune::IdSet!. This may be
 1589 relatively slow because the data type used for ids is usually not an
 1590 \lstinline!int! and the non-consecutive ids require more complicated search data
 1591 structures (typically a map). Evaluation of the map $m$ therefore
 1592 typically costs $O(\log |E^\prime|)$ . On the other hand, id-based
 1593 mappers are not restricted to specific entity sets $E^\prime$.
 1595 In adaptive applications one would use an index-based mapper to do in
 1596 the calculations on a certain grid and only in the adaption phase an
 1597 id-based mapper would be used to transfer the required data
 1598 (e.~g.~only the finite element solution) from one grid to the next grid.
 1600 \section{Visualization of discrete functions}
 1602 Let us use mappers to evaluate a function $f:\Omega\to\mathbb{R}$ for
 1603 certain entities and store the values in a vector. Then, in order to
 1604 do something useful, we use the vector to produce a graphical
 1605 visualization of the function.
 1607 The first example evaluates the function at the centers of all
 1608 elements of the leaf grid and stores this value. Here is the listing:
 1610 \begin{lst}[File dune-grid-howto/elementdata.hh] \mbox{}
 1611 \nopagebreak
 1612 \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
 1613 numberstyle=\tiny, numbersep=5pt]{../elementdata.hh}
 1614 \end{lst}
 1616 The class template
 1617 \lstinline!Dune::LeafMultipleCodimMultipleGeomTypeMapper!  provides an
 1618 index-based mapper where the entities in the subset $E^\prime$ are all
 1619 leaf entities and can further be selected depending on the codimension
 1620 and the geometry type. To that end the constructor takes a function object
 1621 called with a geometry type and the grid dimension and returning a Boolean
 1622 value. When the function object returns true for a combination of geometry
 1623 type and grid dimension then all leaf entities with that geometry type
 1624 will be in the subset $E^\prime$. The mapper object is
 1625 constructed in line \ref{edh:mapper}. A similar mapper is available
 1626 also for the entities of a grid level.
 1628 The data vector is allocated in line \ref{edh:c}. Here we use a
 1629 \lstinline!std::vector<double>!. The \lstinline!size()! method of the
 1630 mapper returns the number of entities in the set $E^\prime$. Instead
 1631 of the STL vector one can use any other type with an
 1632 \lstinline!operator[]!, even built-in arrays (however, built-in arrays
 1633 will not work in this example because the VTK output
 1634 below requires a container with a
 1635 \lstinline!size()! method.
 1637 Now the loop in lines \ref{edh:loop0}-\ref{edh:loop1} iterates through
 1638 all leaf elements. The next three statements within the loop body
 1639 compute the position of the center of the element in global
 1640 coordinates. Then the essential statement is in line \ref{edh:feval}
 1641 where the function is evaluated and the value is assigned to the
 1642 corresponding entry in the \lstinline!c! array. The evaluation of the
 1643 map $m$ is performed by \lstinline!mapper.map(*it)!  where
 1644 \lstinline!*it! is the entity which is passed as a const reference to
 1645 the mapper.
 1647 The remaining lines of code produce graphical output. Lines
 1648 \ref{edh:vtk0}-\ref{edh:vtk1} produce an output file for the
 1649 Visualization Toolkit (VTK), \cite{VTK}, in its XML format. If the
 1650 grid is distributed over several processes the
 1651 \lstinline!Dune::VTKWriter! produces one file per process and the
 1652 corresponding XML metafile. Using Paraview, \cite{Paraview}, you can
 1653 visualize these files. Lines \ref{edh:grape0}-\ref{edh:grape1} enable
 1654 online interactive visualization with the Grape, \cite{Grape},
 1655 graphics package, if it is installed on your machine.
 1657 The next list shows a function \lstinline!vertexdata! that does the
 1658 same job except that the data is associated with the vertices of the
 1659 grid.
 1661 \begin{lst}[File dune-grid-howto/vertexdata.hh] \mbox{}
 1662 \nopagebreak
 1663 \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
 1664 numberstyle=\tiny, numbersep=5pt]{../vertexdata.hh}
 1665 \end{lst}
 1667 The differences to the \lstinline!elementdata! example are the
 1668 following:
 1669 \begin{itemize}
 1670 \item In the \lstinline!P1Layout! struct the method
 1671   \lstinline!contains! returns true if \lstinline!codim==dim!.
 1672 \item Use a leaf iterator for codimension \lstinline!dim! instead of
 1673   \lstinline!0!.
 1674 \item Evaluate the function at the vertex position which is directly
 1675   available via \lstinline!it->geometry()[0]!.
 1676 \item Use \lstinline!addVertexData! instead of \lstinline!addCellData!
 1677   on the \lstinline!Dune::VTKWriter!.
 1678 \item Pass \lstinline!polynomialOrder=1! instead of 0 as
 1679 the second last argument of \lstinline!grape.displayVector!. This argument is the polynomial
 1680   degree of the approximation.
 1681 \end{itemize}
 1683 Finally the following listing shows the main program that calls the
 1684 two functions just discussed:
 1686 \begin{lst}[File dune-grid-howto/visualization.cc] \mbox{}
 1687 \nopagebreak
 1688 \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
 1689 numberstyle=\tiny, numbersep=5pt]{../visualization.cc}
 1690 \end{lst}
 1692 \section{Cell centered finite volumes}
 1693 \label{Sec:CellCenteredFV}
 1695 In this section we show a first complete example for the numerical
 1696 solution of a partial differential equation (PDE), although a very simple
 1697 one.
 1699 We will solve the linear hyperbolic PDE
 1700 \begin{equation}
 1701 \frac{\partial c}{\partial t} + \nabla\cdot (uc) = 0 \ \ \text{ in
 1702   $\Omega\times T$}
 1703 \label{Eq:TransportEquation}
 1704 \end{equation}
 1705 where $\Omega\subset\mathbb{R}^d$ is a domain, $T=(0,t_{\text{end}})$
 1706 is a time interval,
 1707 $c:\Omega\times T\to\mathbb{R}$ is the unknown concentration and
 1708 $u:\Omega\times T\to\mathbb{R}^d$ is a given velocity field. We
 1709 require that the velocity field is divergence free for all times.
 1710 The equation is subject to the initial condition
 1711 \begin{equation}
 1712 c(x,0) = c_0(x) \ \ x\in\Omega
 1713 \end{equation}
 1714 and the boundary condition
 1715 \begin{equation}
 1716 c(x,t) = b(x,t) \ \ t>0, x\in\Gamma_{\text{in}}(t)=\{y\in\partial\Omega\
 1717 |\ u(y,t)\cdot\nu(y)<0\}.
 1718 \end{equation}
 1719 Here $\nu(x)$ is the unit outer normal at a point $y\in\partial\Omega$
 1720 and $\Gamma_{\text{in}}(t)$ is the inflow boundary at time $t$.
 1722 \subsection{Numerical Scheme}
 1724 To keep the presentation simple we use a cell-centered finite volume
 1725 discretization in space, full upwind evaluation of the fluxes and an
 1726 explicit Euler scheme in time.
 1728 The grid consists of cells (elements) $\omega$ and the time interval
 1729 $T$ is discretized into discrete steps $0=t_0, t_1, \ldots, t_n,
 1730 t_{n+1},\ldots, t_N=t_{\text{end}}$. Cell centered finite volume schemes
 1731 integrate the PDE (\ref{Eq:TransportEquation}) over a cell $\omega_i$
 1732 and a time interval $(t_n,t_{n+1})$:
 1733 \begin{equation}
 1734 \int\limits_{\omega_i}\int\limits_{t_n}^{t_{n+1}}\frac{\partial
 1735   c}{\partial t}\ dt\,dx +
 1736 \int\limits_{\omega_i}\int\limits_{t_n}^{t_{n+1}} \nabla\cdot (uc) \
 1737 dt\,dx = 0 \ \ \forall i.
 1738 \label{Eq:TransportEquationIntegrated}
 1739 \end{equation}
 1740 Using integration by parts we arrive at
 1741 \begin{equation}
 1742 \int\limits_{\omega_i} c(x,t_{n+1})\ dx - \int\limits_{\omega_i}
 1743 c(x,t_{n})\ dx + \int\limits_{t_n}^{t_{n+1}}
 1744 \int\limits_{\partial\omega_i} c u\cdot\nu\ ds\,dt = 0 \ \ \forall i.
 1745 \end{equation}
 1746 Now we approximate $c$ by a cell-wise constant function $C$, where
 1747 $C_i^n$ denotes the value in cell $\omega_i$ at time $t_n$.
 1748 Moreover we subdivide the boundary $\partial\omega_i$ into
 1749 facets $\gamma_{ij}$ which are either intersections with
 1750 other cells $\partial\omega_i\cap\partial\omega_j$,
 1751 or intersections with the boundary
 1752 $\partial\omega_i\cap\partial\Omega$.
 1753 Evaluation of the fluxes
 1754 at time level $t_n$ leads to the following equation for the unknown
 1755 cell values at $t_{n+1}$:
 1756 \begin{equation}
 1757 C_i^{n+1}|\omega_i| - C_i^{n}|\omega_i| + \sum\limits_{\gamma_{ij}}
 1758 \phi(C_i^n, C_j^n, u_{ij}^n\cdot\nu_{ij}; \gamma_{ij}, t_n )
 1759 |\gamma_{ij}| \Delta t^n = 0 \ \ \forall i,
 1760 \label{Eq:DiscreteForm}
 1761 \end{equation}
 1762 where $\Delta t^n=t_{n+1}-t_n$, $u_{ij}^n$ is the velocity on the facet $\gamma_{ij}$ at time
 1763 $t_n$, $\nu_{ij}$ is the unit outer normal of the facet $\gamma_{ij}$ and $\phi$ is the flux function defined as
 1764 \begin{equation}
 1765 \phi(C_i^n, C_j^n, u_{ij}^n\cdot\nu_{ij}; \gamma_{ij}, t_n ) =
 1766 \left\{\begin{array}{ll}
 1767 b(\gamma_{ij})\, u_{ij}^n\cdot\nu_{ij} & \gamma_{ij} \subset \Gamma_{\text{in}}(t) \\
 1768 C_j^n\, u_{ij}^n\cdot\nu_{ij}&
 1769 \gamma_{ij}=\partial\omega_i\cap\partial\omega_j\,\wedge\,u_{ij}^n\cdot\nu_{ij} < 0\\
 1770 C_i^n\, u_{ij}^n\cdot\nu_{ij}& u_{ij}^n\cdot\nu_{ij} \geq 0
 1771 \end{array}\right..
 1772 \end{equation}
 1773 Here $b(\gamma_{ij})$ denotes evaluation of the boundary condition on
 1774 an inflow facet $\gamma_{ij}$. If we formally set
 1775 $C_j^n=b(\gamma_{ij})$ on an inflow facet $\gamma_{ij} \subset
 1776 \Gamma_{\text{in}}(t)$ we can derive the following shorthand notation
 1777 for the flux function:
 1778 \begin{equation}
 1779 \phi(C_i^n, C_j^n, u_{ij}^n\cdot\nu_{ij}; \gamma_{ij}, t_n ) = C_i^n
 1780 \max(0,u_{ij}^n\cdot\nu_{ij}) - C_j^n \max(0,-u_{ij}^n\cdot\nu_{ij}).
 1781 \end{equation}
 1782 Inserting this into (\ref{Eq:DiscreteForm}) and solving for
 1783 $C_i^{n+1}$ we obtain
 1784 \begin{equation}
 1785 C_i^{n+1} = C_i^n\left(1-\Delta t^n \sum\limits_{\gamma_{ij}}
 1786   \frac{|\gamma_{ij}|}{|\omega_i|}
 1787   \max(0,u_{ij}^n\cdot\nu_{ij})\right) + \Delta t^n
 1788 \sum\limits_{\gamma_{ij}} C_j^n  \frac{|\gamma_{ij}|}{|\omega_i|}
 1789 \max(0,-u_{ij}^n\cdot\nu_{ij}) \ \ \forall i.
 1790 \label{Eq:IterativeForm1}
 1791 \end{equation}
 1792 One can show that the scheme is stable provided the following condition holds:
 1793 \begin{equation}
 1794 \forall i:\ 1-\Delta t^n \sum\limits_{\gamma_{ij}}
 1795   \frac{|\gamma_{ij}|}{|\omega_i|}
 1796   \max(0,u_{ij}^n\cdot\nu_{ij})\geq 0 \ \Leftrightarrow\ \Delta
 1797   t^n\leq  \min_i \left(\sum\limits_{\gamma_{ij}}
 1798   \frac{|\gamma_{ij}|}{|\omega_i|}
 1799   \max(0,u_{ij}^n\cdot\nu_{ij})\right)^{-1}.
 1800 \label{Eq:TimeStepControl}
 1801 \end{equation}
 1802 When we rewrite \ref{Eq:IterativeForm1} in the form
 1803 \begin{equation}
 1804 C_i^{n+1} = C_i^n - \Delta t^n
 1805   \underbrace{\sum\limits_{\gamma_{ij}} \frac{|\gamma_{ij}|}{|\omega_i|}
 1806   \left( C_i^n \max(0,u_{ij}^n\cdot\nu_{ij}) +
 1807          C_j^n \max(0,-u_{ij}^n\cdot\nu_{ij})\right)}_{\delta_i} \ \ \forall i
 1808 \label{Eq:IterativeForm2}
 1809 \end{equation}
 1810 then it becomes clear that the optimum time step $\Delta t^n$ and the
 1811 update $\delta_i$ for each cell can be computed in a single iteration
 1812 over the grid. The computation $C^{n+1} = C^n - \Delta t^n \delta$ can then
 1813 be realized with a simple vector update. In this form, the algorithm
 1814 can also be parallelized in a straightforward way.
 1817 \subsection{Implementation}
 1819 First, we need to specify the problem parameters, i.e.~initial
 1820 condition, boundary condition and velocity field. This is done by the
 1821 following functions.
 1823 \begin{lst}[File dune-grid-howto/transportproblem.hh] \mbox{}
 1824 \nopagebreak
 1825 \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
 1826 numberstyle=\tiny, numbersep=5pt]{../transportproblem.hh}
 1827 \end{lst}
 1829 The initialization of the concentration vector with the initial
 1830 condition should also be straightforward now. The function
 1831 \lstinline!initialize! works on a concentration vector \lstinline!c!
 1832 that can be stored in any container type with a vector interface
 1833 (\lstinline!operator[]!, \lstinline!size()! and copy constructor are
 1834 needed). Moreover the grid and a mapper for element-wise data have to be passed
 1835 as well.
 1837 \begin{lst}[File dune-grid-howto/initialize.hh] \mbox{}
 1838 \nopagebreak
 1839 \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
 1840 numberstyle=\tiny, numbersep=5pt]{../initialize.hh}
 1841 \end{lst}
 1843 The main work is now done in the function which implements the
 1844 evolution (\ref{Eq:IterativeForm2}) with optimal time step control via
 1845 (\ref{Eq:TimeStepControl}). In addition to grid, mapper and
 1846 concentration vector the current time $t_n$ is passed and the optimum
 1847 time step $\Delta t^n$ selected by the algorithm is returned.
 1848 % Please don't enter an empty line.. otherwise tex4ht dies :-(
 1849 \begin{lst}[File dune-grid-howto/evolve.hh] \mbox{} \label{List:evolve}
 1850 \nopagebreak
 1851 \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
 1852 numberstyle=\tiny, numbersep=5pt]{../evolve.hh}
 1853 \end{lst}
 1855 \begin{figure}
 1856 \begin{center}
 1857 \begin{overpic}[width=0.48\textwidth]{img/intersection}
 1858 \put(60,38){$\omega_i$}
 1859 \put(22,55){$\omega_j$}
 1860 \put(43,46){$\gamma_{ij}$}
 1861 \put(82,52){$\gamma_{ik}$}
 1862 \end{overpic}
 1863 \end{center}
 1864 \caption{Intersection with other elements and the boundary}
 1865 \label{Fig:IsIt}
 1866 \end{figure}
 1868 Lines \ref{evh:loop0}-\ref{evh:loop1} contain the loop over all leaf
 1869 elements where the optimum $\Delta t^n$ and the cell updates
 1870 $\delta_i$ are computed. The update vector is allocated in line
 1871 \ref{evh:update}, where we assume that \lstinline!V! is a container
 1872 with copy constructor and size method.
 1874 The computation of the fluxes is done in lines
 1875 \ref{evh:flux0}-\ref{evh:flux1}. An \lstinline!IntersectionIterator!
 1876 is used to access all intersections $\gamma_{ij}$ of a leaf element
 1877 $\omega_i$. For a full documentation on the
 1878 \lstinline!Intersection!\ class, we refer to the doxygen module page on
 1879 Intersections\footnote{\url{http://www.dune-project.org/doxygen/master/classDune_1_1IntersectionIterator.html}}
 1880 An Intersection is with another element $\omega_j$ if the
 1881 \lstinline!neighbor()! method of the iterator returns true (line
 1882 \ref{evh:neighbor}) or with the external boundary if
 1883 \lstinline!boundary()! returns true (line \ref{evh:bndry}), see also
 1884 Figure~\ref{Fig:IsIt}. An intersection $\gamma_{ij}$ is
 1885 described by several mappings: (i) from a reference element of the
 1886 intersection (with a dimension equal to the grid's dimension minus 1)
 1887 to the reference elements of the two elements $\omega_i$ and
 1888 $\omega_j$ and (ii) from a reference element of the intersection to
 1889 the global coordinate system (with the world dimension). If an
 1890 intersection is with another element then the \lstinline!outside()!
 1891 method returns an \lstinline!EntityPointer! to an entity of
 1892 codimension 0.
 1894 The $\Delta t^n$ calculation is done in line \ref{evh:dt} where the
 1895 minimum over all cells is taken. Then, line \ref{evh:.99} multiplies
 1896 the optimum $\Delta t^n$ with a safety factor to avoid any instability
 1897 due to round-off errors.
 1899 Finally, line \ref{evh:updc} computes the new concentration by adding
 1900 the scaled update to the current concentration.
 1902 The function \lstinline!vtkout! in the following listing provides an
 1903 output of the grid and the solution using the Visualization Toolkit's \cite{VTK}
 1904 XML file format.
 1906 \begin{lst}[File dune-grid-howto/vtkout.hh] \mbox{}
 1907 \nopagebreak
 1908 \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
 1909 numberstyle=\tiny, numbersep=5pt]{../vtkout.hh}
 1910 \end{lst}
 1912 In addition to the snapshots that are produced at each timestep, this function
 1913 also generates a series file which stores the actual \lstinline!time! of an
 1914 evolution scheme together with the snapshots' filenames.  After executing the
 1915 shell script \lstinline!writePVD!\ on this series file, we get a
 1916 \underline{P}ara\underline{v}iew {D}ata (PVD) file with the same name as the
 1917 snapshots. This file opened with paraview then gives us a neat animation over
 1918 the time period.
 1920 Finally, the main program:
 1922 \begin{lst}[File dune-grid-howto/finitevolume.cc] \mbox{}
 1923 \nopagebreak
 1924 \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
 1925 numberstyle=\tiny, numbersep=5pt]{../finitevolume.cc}
 1926 \end{lst}
 1928 The function \lstinline!timeloop! constructs a mapper and allocates the
 1929 concentration vector with one entry per element in the leaf grid. In
 1930 line \ref{fvc:init} this vector is initialized with the initial
 1931 concentration and the loop in line \ref{fvc:loop0}-\ref{fvc:loop1}
 1932 evolves the concentration in time. Every time the current time crosses a multiple of 
 1933 \lstinline!saveStep!\,$=0.1$, the simulation result is
 1934 written to a file in line \ref{fvc:file}.
 1936 \section{A FEM example: The Poisson equation}
 1937 \label{Sec:FEMPoisson}
 1939 In this section we will put together our knowledge about the \Dune{} grid interface acquired in previous chapters to solve the Poisson equation with Dirichlet boundary conditions on the domain $\Omega = (0,1)^d$:
 1940 \begin{align} \label{eq:fempoisson}
 1941  -\Delta u & = f \ \ in \ \Omega\\
 1942 u & = 0 \ \ on \ \partial\Omega
 1943 \end{align}
 1945 The equation will be solved using P1-Finite-Elements on a simplicial grid. The implementation aims to be easy to understand and yet show the power of the \Dune{} grid interface and its generic approach.
 1947 The starting point of the Finite Element Method is the variational formulation of \ref{eq:fempoisson}, which is obtained by partial integration:
 1948 \begin{equation}
 1949 \underbrace{\int\limits_{\Omega}\nabla\ u\cdot\nabla v dx}_{ =:a(u,v)}
 1950 = \underbrace{\int\limits_{\Omega} fv dx}_{=:\ell (v)} \ \ v\in V_h
 1951 \end{equation}
 1953 Let now $\mathcal{T}$ be a conforming triangulation of the domain $\Omega$ with simplices:
 1954 \renewcommand{\labelenumi}{(\roman{enumi})}
 1955 \begin{center}
 1956 \begin{enumerate}
 1957  \item $\bigcup\limits_{\Delta\in\mathcal{T}} \overline{\Delta} = \overline{\Omega}$
 1958  \item $\Delta_i\cap\Delta_j \ i\neq j$ is an entity of
 1959 higher codimension of the elements $\Delta_i$, $\Delta_j$
 1960 \end{enumerate}
 1961 \end{center}
 1963 As we want to use linear finite elements we choose our test function space to be
 1964 \begin{equation}
 1965  V_h = \{ u\in C(\bar{\Omega} )\ \Big| \ u\big|_{\Delta}\in\mathcal{P}_1 (\Delta ) \ \forall\Delta\in\mathcal{T}\}
 1966 \end{equation}
 1968 We will not incorporate the Dirichlet boundary conditions into this function space. Instead, we will implement them in an easier way as described later on.
 1970 As a basis $\phi_1 ,...,\phi_N$ of $V_h$ we choose the nodal basis - providing us small supports and thus a sparse stiffness matrix. After transformation onto the reference element we can use the shape functions
 1971 \begin{align}
 1972 N_0(x) & = 1 - \sum\limits_{i=1}^d x_i \\
 1973 N_i(x) & = x_i \ \ i=1,...,d
 1974 \end{align}
 1976 to evaluate the basis functions and their gradients.
 1978 The numerical solution $u_h$ is a linear combination of $\phi_1,...,\phi_N$ with coefficients $u_1,...,u_N$. We assemble the stiffness Matrix $A$ and the vector $b$:
 1979 \begin{equation}
 1980  A_{ij}=a(\phi_i ,\phi_j ) \ \ b_i=\ell (\phi_i )
 1981 \end{equation}
 1983 The coefficients $u_1,...,u_N$ are then obtained by solving $Au=b$.
 1985 The integrals are transformed onto the reference element $\hat{\Delta}$ and computed with an appropriate quadrature rule. Let the transformation map be given by $g: \hat{\Delta} \to \Delta$
 1987 \begin{align}
 1988  A_{ij} & = \int\limits_{\Omega}\nabla\phi_i\cdot\nabla\phi_j\ dx\\
 1989 & = \sum\limits_{\Delta\in\mathcal{T}} \int\limits_{\Delta}\nabla\phi_i\cdot\nabla\phi_j\ dx\\
 1990 & = \sum\limits_{\Delta\in\mathcal{T}} \int\limits_{\hat{\Delta}}\nabla\phi_i (g(\hat{x}))
 1991     \cdot\nabla\phi_j (g(\hat{x}))\mu (\hat{x} )d\hat{x}\\
 1992 & = \sum\limits_{\Delta\in\mathcal{T}} \int\limits_{\hat{\Delta}} (J_g^{-T}\hat{\nabla}\hat{\phi_i} )(\hat{x} )
 1993     \cdot (J_g^{-T}\hat{\nabla}\hat{\phi_j} )(\hat{x} )\mu (\hat{x} )d\hat{x}
 1994 \end{align}
 1996 $J_g$ is the Jacobian of the map $g$ and $\mu (\hat{x}):=\sqrt{det J_g^{T} J_g}$ the Jacobian determinant. Let now $\xi_k$ be the quadrature points  of the chosen rule and $\omega_k$ the associated weights. We assume that there are $p_1$ quadrature points to evaluate:
 1997 \begin{equation}
 1998 \label{equ:computea}
 1999  \Rightarrow \ A_{ij} = \sum\limits_{\Delta\in\tau} \sum\limits_{k=1}^{p_1}
 2000 \omega_k(J_g^{-T}\hat{\nabla}\hat{\phi_i} )(\xi_k )
 2001 \cdot (J_g^{-T}\hat{\nabla}\hat{\phi_j} )(\xi_k )\mu (\xi_k )
 2002 \end{equation}
 2004 Simultaneously, the right side b is treated in the same manner. As we might want to use another quadrature rule here that better suits our function $f$, we use $p_2$ quadrature points:
 2005 \begin{equation}
 2006 \label{equ:computeb}
 2007  b_i = \sum\limits_{\Delta\in\tau}\sum\limits_{k=1}^{p_2}\omega_k f(g(\xi_k ))\phi_i (g(\xi_k ))\mu(\xi_k )
 2008 \end{equation}
 2010 In our implementation we will of course not compute the matrix entries one after another but rather iterate over all elements of the grid and update all matrix entries with a non-vanishing contribution on that element.
 2012 After assembling the matrix we implement the Dirichlet boundary conditions by overwriting the lines of the equation system associated with boundary nodes with trivial lines, see figure~\ref{Fig:dirichlet}.
 2014 \begin{figure}
 2015 \centering
 2016 \begin{tabular}{ccl}
 2017 & $\begin{array}{cccccccc}
 2018 & & &i\downarrow& & &\\
 2019 \end{array}$ & \\
 2020 $\rightarrow{i}$ &
 2021 $\left(\begin{array}{ccccccc}
 2022 \star&\star&\star&\star&\star&\star&\star\\
 2023 0&\cdots&0&1&0&\cdots&0\\
 2024 \star&\star&\star&\star&\star&\star&\star\\
 2025 \end{array}\right)$ &
 2026 $\left(\begin{array}{c}
 2027 \star \\
 2028 u_i \\
 2029 \star \\
 2030 \end{array}\right)
 2031 =
 2032 \left(\begin{array}{c}
 2033 \star \\
 2034 0 \\
 2035 \star \\
 2036 \end{array}\right)$
 2037 \end{tabular}
 2038 \caption{Lines of $A$ and $b$ are replaced by trivial lines.}
 2039 \label{Fig:dirichlet}
 2040 \end{figure}
 2042 This is possible as---using the nodal basis---the coefficients match the value of the numerical solution at the corresponding node.
 2044 \subsection{Implementation}
 2046 In this implementation we will restrict ourselves to a 2-dimensional grid. However, the code works on simplicial grids of any dimension. Try this later!
 2048 Lets first have a look at the implementation of the shape functions. This class only provides the methods to evaluate the shape functions and their gradients:
 2050 \begin{lst}[File dune-grid-howto/shapefunctions.hh] \mbox{}
 2051 \nopagebreak
 2052 \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
 2053 numberstyle=\tiny, numbersep=5pt]{../shapefunctions.hh}
 2054 \end{lst}
 2056 And now the actual FEM code:
 2058 \begin{lst}[File dune-grid-howto/finiteelements.cc] \mbox{}
 2059 \nopagebreak
 2060 \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
 2061 numberstyle=\tiny, numbersep=5pt, breaklines=true]{../finiteelements.cc}
 2062 \end{lst}
 2064 The function \lstinline!determineAdjacencyPattern()! in lines \ref{fem:adjpat1} to \ref{fem:adjpat2} does traverse the grid and stores all adjacency information in a \lstinline!std::vector< std::set<int> >!. You might wonder why this is necessary before the actual computing of the matrix entries. The reason for this is that, as data structure for the matrix $A$, we use \lstinline!BCRSMatrix! - which is specialized to hold large sparse matrices. Using this type, information about which entries do not vanish has to be known when assembling. We do give this information to the matrix from line \ref{fem:setpattern} on. Only after finishing this in line \ref{fem:endindices} we can start to fill the matrix with values.
 2066 From line \ref{fem:loop1} to \ref{fem:loop2} we have the main loop traversing the whole grid and updating the matrix entries. This does strictly follow the procedure described in previous chapters. The main calculation is done in line \ref{fem:calca} and \ref{fem:calcb} - which are one-to-one implementations of \ref{equ:computea} and \ref{equ:computeb}.
 2068 As already said above, we do directly implement Dirichlet boundaries into our matrix. This is done in lines \ref{fem:boundary1} to \ref{fem:boundary2}. We have to traverse the whole grid once again and check for each intersection of elements whether it is on the boundary. In line \ref{fem:trivialline} we overwrite the line corresponding to a node on the boundary as shown in figure \ref{Fig:dirichlet}.
 2070 When you visualize your results, you should get something like figure \ref{Fig:FEM1} or \ref{Fig:FEM2}!
 2072 \begin{figure}
 2073 \centering
 2074 \includegraphics[width=0.8\textwidth]{img/fem2d}
 2075 \caption{Solution in 2D}
 2076 \label{Fig:FEM1}
 2077 \end{figure}
 2079 \begin{figure}
 2080 \centering
 2081 \includegraphics[width=0.8\textwidth]{img/fem3d}
 2082 \caption{Solution in 3D}
 2083 \label{Fig:FEM2}
 2084 \end{figure}
 2086 \begin{exc}
 2087  Try a 3-dimensional grid! Just change the dimension in line \ref{fem:dim} and the name of the gridfile in line \ref{fem:file} to \lstinline!3dgrid.al! . You can compile the new code without reconfiguring by running
 2088 \begin{lstlisting}
 2089 make ALBERTA_DIM=3
 2090 \end{lstlisting}
 2091 \end{exc}
 2093 \begin{exc}
 2094  Modify the code in order to make it handle Neumann boundary conditions too!
 2095 \end{exc}
 2099 \chapter{Adaptivity}
 2101 \section{Adaptive integration}
 2103 \subsection{Adaptive multigrid integration}
 2105 In this section we describe briefly the adaptive multigrid integration
 2106 algorithm presented in \cite{Deuflhard93}.
 2108 \minisec{Global error estimation}
 2110 The global error can be estimated by taking the difference of the numerically
 2111 computed value for the integral on a fine and a coarse grid as given in
 2112 (\ref{Eq:GlobalError}).
 2114 \minisec{Local error estimation}
 2116 Let $I_f^p(\omega)$ and $I_f^q(\omega)$ be two integration formulas of
 2117 different orders $p>q$ for the evaluation of the integral over some
 2118 function $f$ on the element $\omega\subseteq\Omega$. If we assume that
 2119 the higher order rule is locally more accurate then
 2120 \begin{equation}
 2121 \bar{\epsilon}(\omega) = |I_f^p(\omega)-I_f^q(\omega)|
 2122 \end{equation}
 2123 is an estimator for the local error on the element $\omega$.
 2125 \minisec{Refinement strategy}
 2127 If the estimated global error is not below a user tolerance the grid
 2128 is to be refined in those places where the estimated local error is
 2129 ``high''. To be more specific, we want to achieve that each element in
 2130 the grid contributes about the same local error to the global
 2131 error. Suppose we knew the maximum local error on all the new elements
 2132 that resulted from refining the current mesh (without actually doing
 2133 so). Then it would be a good idea to refine only those elements in the
 2134 mesh where the local error is not already below that maximum local
 2135 error that will be attained anyway.
 2136 In \cite{Deuflhard93} it is shown that the local error after mesh
 2137 refinement can be effectively computed without actually doing the
 2138 refinement. Consider an element $\omega$ and its father element
 2139 $\omega^-$, i.~e.~the refinement of  $\omega^-$ resulted in
 2140 $\omega$. Moreover, assume that $\omega^+$ is a (virtual) element that
 2141 would result from a refinement of $\omega$. Then it can be shown that
 2142 under certain assumptions the quantity
 2143 \begin{equation}
 2144 \epsilon^+(\omega) = \frac{\bar{\epsilon}(\omega)^2}{\bar{\epsilon}(\omega^-)}
 2145 \end{equation}
 2146 is an estimate for the local error on $\omega^+$,
 2147 i.~e.~$\bar{\epsilon}(\omega^+)$.
 2149 Another idea to determine the refinement threshold is to look simply
 2150 at the maximum of the local errors on the current mesh and
 2151 to refine only those elements where the local error is above a certain
 2152 fraction of the maximum local error.
 2154 By combining the two approaches we get the threshold value $\kappa$
 2155 actually used in the code:
 2156 \begin{equation}
 2157 \kappa = \min\left(\max\limits_{\omega} \epsilon^+(\omega), \frac12
 2158   \max\limits_{\omega} \bar{\epsilon}(\omega) \right).
 2159 \end{equation}
 2161 \minisec{Algorithm}
 2163 The complete multigrid integration algorithm then reads as follows:
 2164 \begin{itemize}
 2165 \item Choose an initial grid.
 2166 \item Repeat the following steps
 2167 \begin{itemize}
 2168 \item Compute the value $I$ for the integral on the current grid.
 2169 \item Compute the estimate $E$ for the global error.
 2170 \item If $E<\text{tol}\cdot I$ we are done.
 2171 \item Compute the threshold $\kappa$ as defined above.
 2172 \item Refine all elements $\omega$ where $\bar{\epsilon}(\omega)\geq\kappa$.
 2173 \end{itemize}
 2174 \end{itemize}
 2176 \subsection{Implementation of the algorithm}
 2178 The algorithm above is realized in the following code.
 2180 \begin{lst}[File dune-grid-howto/adaptiveintegration.cc] \mbox{}
 2181 \nopagebreak
 2182 \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
 2183 numberstyle=\tiny, numbersep=5pt]{../adaptiveintegration.cc}
 2184 \end{lst}
 2186 The work is done in the function \lstinline!adaptiveintegration!.
 2187 Lines \ref{aic:int0}-\ref{aic:int1} compute the value of the integral
 2188 on the current mesh. After printing the result the decision whether
 2189 to continue or not is done in line \ref{aic:finish}. The extrapolation
 2190 strategy relies on the fact that every element has a father. To ensure
 2191 this, the grid is at least once refined globally in the first step
 2192 (line \ref{aic:gr}). Now the refinement threshold $\kappa$ can be
 2193 computed in lines \ref{aic:kappa0}-\ref{aic:kappa1}. Finally the last
 2194 loop in lines \ref{aic:mark0}-\ref{aic:mark1} marks elements for
 2195 refinement and lines \ref{aic:ref0}-\ref{aic:ref1} actually do the
 2196 refinement. The reason for dividing refinement into three functions
 2197 \lstinline!preAdapt()!, \lstinline!adapt()! and
 2198 \lstinline!postAdapt()! will be explained with the next example. Note
 2199 the flexibility of this algorithm: It runs in any space dimension on
 2200 any kind of grid and different integration orders can easily be
 2201 incorporated. And that with just about 100 lines of code including
 2204 \begin{figure}
 2205 \includegraphics[width=0.48\textwidth]{img/adaptiveintegration_alberta2d}\hfill
 2206 \includegraphics[width=0.48\textwidth]{img/adaptiveintegration_ug3d}
 2207 \caption{Two and three-dimensional grids generated by the adaptive
 2208   integration algorithm applied to the needle pulse. Left grid is
 2209   generated using Alberta, right grid is generated using UG.}
 2210 \label{Fig:AdaptiveIntegration}
 2211 \end{figure}
 2213 Figure \ref{Fig:AdaptiveIntegration} shows two grids generated by the
 2214 adaptive integration algorithm.
 2216 \begin{warn} The quadrature rules for prisms and pyramids are
 2217   currently only implemented for order two. Therefore adaptive
 2218   calculations with UGGrid and hexahedral elements do not work.
 2219 \end{warn}
 2223 \section{Adaptive cell centered finite volumes}
 2225 In this section we extend the example of Section
 2226 \ref{Sec:CellCenteredFV} by adaptive mesh refinement. This requires
 2227 two things: (i) a method to select cells for refinement or coarsening
 2228 (derefinement) and (ii) the transfer of a solution on a given grid to
 2229 the adapted grid. The finite volume algorithm itself has already been
 2230 implemented for adaptively refined grids in Section
 2231 \ref{Sec:CellCenteredFV}.
 2233 For the adaptive refinement and coarsening we use a very simple
 2234 heuristic strategy that works as follows:
 2235 \begin{itemize}
 2236 \item Compute global maximum and minimum of element concentrations:
 2237   $$\overline{C}=\max_i C_i,  \ \ \ \underline{C}=\min_i C_i.$$
 2238 \item As the local indicator in cell $\omega_i$ we define $$\eta_i =
 2239   \max_{\gamma_{ij}} |C_i-C_j|.$$ Here $\gamma_{ij}$ denotes
 2240   intersections with other elements in the leaf grid.
 2241 \item If for $\omega_i$ we have
 2242   $\eta_i>\overline{\text{tol}}\cdot (\overline{C}-\underline{C})$
 2243   and $\omega_i$ has not been refined more than $\overline{M}$ times
 2244   then mark $\omega_i$ and all its neighbors for refinement.
 2245 \item Mark all elements $\omega_i$ for coarsening where
 2246   $\eta_i<\underline{\text{tol}}\cdot (\overline{C}-\underline{C})$
 2247   and $\omega_i$ has been refined at least $\underline{M}$ times.
 2248 \end{itemize}
 2250 This strategy refines an element if the local gradient is ``large''
 2251 and it coarsens elements (which means it removes a previous
 2252 refinement) if the local gradient is ``small''. In addition any
 2253 element is refined at least refined $\underline{M}$ times and at most
 2254 $\overline{M}$ times.
 2256 After mesh modification the solution from the previous grid must be
 2257 transfered to the new mesh. Thereby the following situations do occur
 2258 for an element:
 2259 \begin{itemize}
 2260 \item The element is a leaf element in the new mesh and was a leaf
 2261   element in the old mesh: keep the value.
 2262 \item The element is a leaf element in the new mesh and existed in the
 2263   old mesh as a non-leaf element: Compute the cell value as an average
 2264   of the son elements in the old mesh.
 2265 \item The element is a leaf element in the new mesh and is obtained
 2266   through refining some element in the old mesh: Copy the value
 2267   from the element in the old mesh to the new mesh.
 2268 \end{itemize}
 2270 The complete mesh adaptation is done by the function
 2271 \lstinline!finitevolumeadapt! in the following listing:
 2273 \begin{lst}[File dune-grid-howto/finitevolumeadapt.hh] \mbox{}
 2274 \nopagebreak
 2275 \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
 2276 numberstyle=\tiny, numbersep=5pt]{../finitevolumeadapt.hh}
 2277 \end{lst}
 2279 The loop in lines \ref{fah:loop0}-\ref{fah:loop1} computes the
 2280 indicator values $\eta_i$ as well as the global minimum and maximum
 2281 $\overline{C},\underline{C}$. Then the next loop in lines
 2282 \ref{fah:loop2}-\ref{fah:loop3} marks the elements for refinement.
 2283 Lines \ref{fah:loop4}-\ref{fah:loop5} construct a map that stores for
 2284 each element in the mesh (on all levels) the average of the element
 2285 values in the leaf elements of the subtree of the given element. This
 2286 is accomplished by descending from the fine grid levels to the coarse
 2287 grid levels and thereby adding the value in an element to the father
 2288 element. The key into the map is the global id of an element. Thus the
 2289 value is accessible also after mesh modification.
 2291 Now the grid can really be modified in line \ref{fah:adapt} by calling the
 2292 \lstinline!adapt()! method on the grid object. The mapper is updated
 2293 to reflect the changes in the grid in line \ref{fah:update} and the
 2294 concentration vector is resized to the new size in line
 2295 \ref{fah:resize}. Then the values have to be interpolated to the new
 2296 elements in the mesh using the map and finally to be transferred to
 2297 the resized concentration vector. This is done in the loop in lines
 2298 \ref{fah:loop6}-\ref{fah:loop7}.
 2300 Here is the new main program with an adapted \lstinline!timeloop!:
 2302 \begin{lst}[File dune-grid-howto/adativefinitevolume.cc] \mbox{}
 2303 \nopagebreak
 2304 \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
 2305 numberstyle=\tiny, numbersep=5pt]{../adaptivefinitevolume.cc}
 2306 \end{lst}
 2308 The program works analogously to the non adaptive \lstinline!finitevolume!\
 2309 version from the previous chapter. The only differences are inside the
 2310 \lstinline!timeloop!\ function. During the initialization of the concentration
 2311 vector in line \ref{afv:in} and after each time step in line \ref{afv:ad} the
 2312 function \lstinline!finitevolumeadapt! is called in order to refine the grid.
 2313 The initial adaptation is repeated $\overline{M}$ times.  Note that adaptation
 2314 after each time steps is deactivated during the compiler phase for unstructured
 2315 grids with help of the \lstinline!Capabilities!\ class. This is because
 2316 structured grids do not allow a conforming refinement and are therefore
 2317 unusable for adaptive schemes. In fact, the \lstinline!adapt! method on a grid
 2318 of \lstinline!YaspGrid! e.\,g.~ results in a {\em global} grid refinement.
 2320 \begin{exc}
 2321   Compile the program with the gridtype set to \lstinline!ALUGRID_SIMPLEX!
 2322   and \lstinline!ALUGRID_CONFORM! and compare the results visually.
 2323 \end{exc}
 2325 \chapter{Parallelism}
 2327 \section{\texorpdfstring{\Dune{}}{DUNE} Data Decomposition Model}
 2329 The parallelization concept in \Dune{} follows the Single Program
 2330 Multiple Data (SPMD) data parallel programming paradigm. In this
 2331 programming model each process executes the same code but on different
 2332 data. The parallel program is parametrized
 2333 by the rank of the individual process in the set and the number of
 2334 processes $P$
 2335 involved. The processes communicate by exchanging messages, but you
 2336 will rarely have the need to bother with sending messages.
 2338 A parallel \Dune{} grid, such as YaspGrid, is a collective object which
 2339 means that all processes participating in the computations on the grid
 2340 instantiate the grid object at the same time (collectively). Each
 2341 process stores a subset of all the entities that the same program running on a
 2342 single process would have. An entity may be stored in more
 2343 than one process, in principle it may be even stored in all
 2344 processes. An entity
 2345 stored in more than one process is called a distributed entity. \Dune{}
 2346 allows quite general data decompositions but not arbitrary data
 2347 decompositions. Each entity in a process has a partition type
 2348 value assigned to it. There are five different possible partition type
 2349 values:
 2350 \begin{center}
 2351 \textit{interior}, \textit{border}, \textit{overlap},
 2352   \textit{front} and \textit{ghost}.
 2353 \end{center}
 2355 Entities of codimension 0 are restricted to the three partition types
 2356 \textit{interior}, \textit{overlap} and \textit{ghost}. Entities of
 2357 codimension greater than 0 may take all partition type values.
 2358 The codimension 0 entities with partition type \textit{interior} form a
 2359 non-overlapping decomposition of the entity set, i.e.~for each
 2360 entity of codimension 0 there is exactly one process where this entity
 2361 has partition type \textit{interior}.
 2362 Moreover, the codimension 0 leaf entities in process number $i$ form a
 2363 subdomain $\Omega_i\subseteq\Omega$ and all the $\Omega_i$, $0\leq i <
 2364 P$, form a nonoverlapping decomposition of the computational domain
 2365 $\Omega$. The leaf entities of codimension 0 in a process $i$ with
 2366 partition types \textit{interior} or \textit{overlap} together form a
 2367 subdomain $\hat{\Omega}_i\subseteq\Omega$.
 2369 Now the partition types of the entities in process $i$
 2370 with codimension greater 0 can
 2371 be determined according to the following table:
 2372 \begin{center}
 2373 \begin{tabular}{cc}
 2374 \hline
 2375 \hline
 2376 Entity located in & Partition Type value\\
 2377 \hline
 2378 $B_i=\overline{\partial\Omega_i\setminus\partial\Omega}$ &
 2379 \textit{border}\\
 2380 $\overline{\Omega_i}\setminus B_i$ & \textit{interior}\\
 2381 $F_i =
 2382 \overline{\partial\hat{\Omega}_i\setminus\partial\Omega}\setminus B_i$
 2383 & \textit{front}\\
 2384 $\overline{\hat{\Omega}_i}\setminus(B_i\cup F_i)$ & \textit{overlap}\\
 2385 Rest & \textit{ghost}\\
 2386 \hline
 2387 \hline
 2388 \end{tabular}
 2389 \end{center}
 2391 \begin{figure}
 2392 \centering
 2393 \begin{overpic}[width=\textwidth]{img/partitionsingle}
 2394  \put(12,62){$c=0$}
 2395  \put(45,62){$c=1$}
 2396  \put(80,62){$c=2$}
 2398  \put(12,31){$c=0$}
 2399  \put(45,31){$c=1$}
 2400  \put(80,31){$c=2$}
 2402  \put(12,1){$c=0$}
 2403  \put(45,1){$c=1$}
 2404  \put(80,1){$c=2$}
 2405 \end{overpic}
 2407 \caption{Color coded illustration of different data decompositions:
 2408   interior (red), border (blue), overlap (green), front (magenta) and
 2409   ghost (yellow), gray encodes entities not stored by the
 2410   process. First row shows case with interior, overlap and ghost
 2411   entities, second row shows a case with interior and overlap without
 2412   ghost and the last row shows a case with interior and ghost only.}
 2413 \label{Fig:PartitionSingle}
 2414 \end{figure}
 2416 The assignment of partition types is illustrated for three different
 2417 examples in Figure \ref{Fig:PartitionSingle}. Each example shows a
 2418 two-dimensional structured grid with $6\times 4$ elements (in
 2419 gray). The entities stored in some process $i$ are shown in color,
 2420 where color indicates the partition type as explained in the
 2421 caption. The first row shows an example where process $i$ has
 2422 codimension 0 entities of all three partition types \textit{interior},
 2423 \textit{overlap} and \textit{ghost} (leftmost picture in first
 2424 row). The corresponding assignment of partition types to entities of
 2425 codimension 1 and 2 is then shown in the middle and right most
 2426 picture. A grid implementation can choose to omit the partition type
 2427 \textit{overlap}  or \textit{ghost} or both, but not
 2428 \textit{interior}. The middle row shows an example where an
 2429 \textit{interior} partition is extended by an \textit{overlap} and no
 2430 \textit{ghost} elements are present. This is the model used in
 2431 YaspGrid. The last row shows an example where the \textit{interior} partition
 2432 is extended by one row of  \textit{ghost} cells. This is the model
 2433 used in UGGrid and ALUGrid.
 2436 \section{Communication Interfaces}
 2438 This section explains how the exchange of data between the partitions
 2439 in different processes is organized in a flexible and portable way.
 2441 The abstract situation is that data has to be sent from a copy of a
 2442 distributed entity in a process to one or more copies of the same
 2443 entity in other processes. Usually data has to be sent not only for
 2444 one entity but for many entities at a time, thus it is more efficient
 2445 pack all data that goes to the same destination process into a single
 2446 message. All entities for which data has to be sent or received form a
 2447 so-called \textit{communication interface}. As an example let us define the set
 2448 $X_{i,j}^c$ as the set of all entities of codimension $c$ in process $i$
 2449   with partition type \textit{interior} or \textit{border} that have
 2450 a copy in process $j$ with any partition type. Then in the
 2451 communication step process $i$ will send one message to any other
 2452 process $j$ when $X_{i,j}^c\neq\emptyset$. The message contains some data for
 2453 every entity in  $X_{i,j}^c$. Since all processes participate in the
 2454 communication step, process $i$ will receive data from a process $j$
 2455 whenever $X_{j,i}^c\neq\emptyset$. This data corresponds to entities
 2456 in process $i$ that have a copy in $X_{j,i}^c$.
 2458 A \Dune{} grid offers a selection of predefined interfaces. The example
 2459 above would use the parameter \lstinline!InteriorBorder_All_Interface!
 2460 in the communication function. After the
 2461 selection of the interface it remains to specify the data to be sent
 2462 per entity and how the data should be processed at the receiving
 2463 end. Since the data is in user space the user has to write a small
 2464 class that encapsulates the processing of the data at the sending and
 2465 receiving end. The following listing shows an example for a so-called
 2466 data handle:
 2469 \begin{lst}[File dune-grid-howto/parfvdatahandle.hh] \mbox{}
 2470 \nopagebreak
 2471 \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
 2472 numberstyle=\tiny, numbersep=5pt]{../parfvdatahandle.hh}
 2473 \end{lst}
 2475 Every instance of the \lstinline!VectorExchange! class template
 2476 conforms to the data handle concept. It defines a type
 2477 \lstinline!DataType! which is the type of objects that are exchanged
 2478 in the messages between the processes. The method \lstinline!contains!
 2479 should return true for all codimensions that participate in the data
 2480 exchange. Method \lstinline!fixedsize! should return true when, for
 2481 the given codimension, the same number of data items per entity is
 2482 sent. If \lstinline!fixedsize! returns false the method
 2483 \lstinline!size! is called for each entity in order to ask
 2484 for the number of items of type \lstinline!DataType! that are to be
 2485 sent for the given entity. Note that this information has only to be
 2486 given at the sender side. Then the method \lstinline!gather! is called
 2487 for each entity in a communication interface on the sender side in
 2488 order to pack the data for this entity into the message buffer. The
 2489 message buffer itself is realized as an output stream that accepts data
 2490 of type \lstinline!DataType!. After exchanging the data via message
 2491 passing the \lstinline!scatter! method is called for each entity at
 2492 the receiving end. Here the data is read from the message buffer and
 2493 stored in the user's data structures. The message buffer is realized
 2494 as an input stream delivering items of type \lstinline!DataType!. In
 2495 the \lstinline!scatter! method it is up to the user how the data is to
 2496 be processed, e.\,g.~one can simply overwrite (as is done here), add or
 2497 compute a maximum.
 2499 \section{Parallel finite volume scheme}
 2501 In this section we parallelize the (nonadaptive!) cell centered finite volume
 2502 scheme. Essentially only the \lstinline!evolve! method has to be
 2503 parallelized. The following listing shows the parallel version of this
 2504 method. Compare this with listing \ref{List:evolve} on page \pageref{List:evolve}.
 2506 \begin{lst}[File dune-grid-howto/parevolve.hh] \mbox{}
 2507 \nopagebreak
 2508 \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
 2509 numberstyle=\tiny, numbersep=5pt]{../parevolve.hh}
 2510 \end{lst}
 2512 The first difference to the sequential version is in line
 2513 \ref{peh:assert} where it is checked that the grid provides an overlap
 2514 of at least one element. The overlap may be either of partition type
 2515 \textit{overlap} or \textit{ghost}. The finite volume scheme itself
 2516 only computes the updates for the elements with partition type
 2517 \textit{interior}.
 2519 In order to iterate over entities with a specific partiton type the
 2520 leaf and level iterators can be parametrized by an additional argument
 2521 \lstinline!PartitionIteratorType! as shown in line \ref{peh:pit}. If
 2522 the argument \lstinline!All_Partition! is given then all entities are
 2523 processed, regardless of their partition type. This is also the
 2524 default behavior of the level and leaf iterators. If the partition
 2525 iterator type is specified explicitely in an iterator the same
 2526 argument has also to be specified in the begin and end methods on the
 2527 grid as shown in lines \ref{peh:end}-\ref{peh:begin}.
 2529 The next change is in line \ref{peh:inter} where the computation of
 2530 the optimum stable time step is restricted to elements of partition
 2531 type \textit{interior} because only those elements have all neighboring
 2532 elements locally available. Next, the global minimum of the time steps
 2533 sizes determined in each process is taken in line \ref{peh:min}. For
 2534 collective communication each grid returns a collective communication
 2535 object with its \lstinline!comm()! method which allows to compute
 2536 global minima and maxima, sums, broadcasts and other functions.
 2538 Finally the updates computed on the \textit{interior} cells in each
 2539 process have to be sent to all copies of the respective entities in
 2540 the other processes. This is done in lines
 2541 \ref{peh:dist0}-\ref{peh:dist1} using the data handle described above.
 2542 The \lstinline!communicate! method on the grid uses the data handle to
 2543 assemble the message buffers, exchanges the data and writes the data
 2544 into the user's data structures.
 2546 Finally, we need a new main program, which is in the following listing:
 2548 \begin{lst}[File dune-grid-howto/parfinitevolume.cc] \mbox{}
 2549 \nopagebreak
 2550 \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
 2551 numberstyle=\tiny, numbersep=5pt]{../parfinitevolume.cc}
 2552 \end{lst}
 2554 A difference to the sequential program can be found in line
 2555 \ref{pfc:rank0} where the printing of the data of the current time
 2556 step is restricted to the process with rank 0.
 2557 \lstinline!YaspGrid! does not support dynamical load balancing and therefore
 2558 needs to start with a sufficiently fine grid that allows a reasonable partition
 2559 where each processes gets a non-empty part of grid. This is why we do not use
 2560 DGF Files in the parallel example and initialize the grid by the UnitCube class
 2561 instead. For \lstinline!YaspGrid! this allows an easy selection of the grid's
 2562 initial coarseness through the second template argument of the
 2563 \lstinline!UnitCube!. This argument should be chosen sufficiently high, because
 2564 after each global refinement step the overlap region grows and therefore the
 2565 communicaton overhead increases.
 2567 If you want to use a grid with support for dynamical load balancing, uncomment
 2568 one of the possible definitions for such a grid in the code and define the
 2569 macro \lstinline!LOAD_BALANCING!. In this case in line \ref{pfv:lb} the method
 2570 \lstinline!loadBalance!\ is called on the grid.
 2571 This method re-partitions the grid in a way such that on every partition there
 2572 is an equal amount of grid elements.
 2574 % \chapter{Input and Output}
 2576 % \section{Visualization with Grape}
 2580 % \section{Visualization with VTK}
 2582 % \section{Dune portable grid format}
 2584 % Dune provides a unified macro grid format called the
 2585 % \textbf{D}une \textbf{G}rid \textbf{F}ormat, short the \textbf{DGF}.
 2586 % This enables the user to use all types of grids with one and the same
 2587 % macro grid file. For the documentation of the DGF and the usage of the parser we refer to
 2589 % \href{http://www.dune-project.org/doc/doxygen/dune-grid-html/group\_\_DuneGridFormatParser.html}%
 2590 % {{\small\texttt{http://www.dune-project.org/doc/doxygen/dune-grid-html/group\_\_DuneGridFormatParser.html}}}
 2592 % in the dune-grid documentation.
 2594 % \section{Amiramesh output}
 2597 % \chapter{Iterative Solver Template Library}
 2599 % \begin{lst}[File dune-grid-howto/groundwater.cc] \mbox{}
 2600 % \nopagebreak
 2601 % \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
 2602 % numberstyle=\tiny, numbersep=5pt]{../groundwater.cc}
 2603 % \end{lst}
 2605 % \begin{lst}[File dune-grid-howto/groundwaterproblem.hh] \mbox{}
 2606 % \nopagebreak
 2607 % \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
 2608 % numberstyle=\tiny, numbersep=5pt]{../groundwaterproblem.hh}
 2609 % \end{lst}
 2611 % \begin{lst}[File dune-grid-howto/groundwateradapt.hh] \mbox{}
 2612 % \nopagebreak
 2613 % \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
 2614 % numberstyle=\tiny, numbersep=5pt]{../groundwateradapt.hh}
 2615 % \end{lst}
 2618 % \chapter{Outlook}
 2621 \begin{figure}
 2622 \centering
 2623 \includegraphics[width=0.32\textwidth]{img/alberta2d}\
 2624 \includegraphics[width=0.32\textwidth]{img/ugsimplex2d}\
 2625 \includegraphics[width=0.32\textwidth]{img/ugcube2d}
 2627 \includegraphics[width=0.32\textwidth]{img/alberta3d}\
 2628 \includegraphics[width=0.32\textwidth]{img/alusimplex3d}\
 2629 \includegraphics[width=0.32\textwidth]{img/alucube3d}
 2632 \includegraphics[width=0.32\textwidth]{img/ugsimplex3d}\
 2633 \includegraphics[width=0.32\textwidth]{img/ugcube3d}\
 2634 \includegraphics[width=0.32\textwidth]{img/iso}
 2636 \caption{Adaptive solution of an elliptic model problem with $P_1$
 2637 conforming finite elements and residual based error
 2638 estimator. Illustrates that adaptive finite element algorithm can be formulated
 2639 independent of dimension, element type and refinement scheme. From
 2640 top to bottom, left to right: Alberta (bisection, 2d), UG (red/green
 2641 on triangles), UG (red/green on quadrilaterals), Alberta (bisection,
 2642 3d), ALU (hanging nodes on tetrahedra), ALU (hanging nodes on
 2643 hexahedra), UG (red/green on tetrahedra), UG (red/green on hexahedra,
 2644 pyramids and tetrahedra), isosurfaces of solution.}
 2645 \label{Fig:OutlookP1}
 2646 \end{figure}
 2650 \bibliographystyle{plain}
 2651 \bibliography{grid-howto.bib}
 2653 \end{document}