## "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,
15 \usepackage{hyperref}
16 \usepackage{makeidx}
17
18 \usepackage{graphicx}
19
20 \DeclareGraphicsExtensions{.pdf, .png, .jpg}
21
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}
28
29 %The theorems
30 \theorembodyfont{\upshape}
32 \newtheorem{exc}{Exercise}[chapter]
33 \newtheorem{rem}[exc]{Remark}
34 \newtheorem{lst}{Listing}
35 \newtheorem{warn}[exc]{Warning}
36
38
39 \title{The Distributed and Unified Numerics Environment (\Dune{}) Grid
40   Interface HOWTO}
41
42 \input{config.inc}
43
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$}
52
53 \date{Version \version, \today}
54
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
60 \includegraphics[width=0.32\textwidth]{img/alucube3d}
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 }
83
84 \begin{document}
85
86 \maketitle
87
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!
102
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}
107
108 \tableofcontents
109
110
111 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
112 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
113 \chapter{Introduction}
114 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
115 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
116
117 \section{What is \texorpdfstring{\Dune{}}{DUNE} anyway?}
118
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}
143
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
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}
171
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).
180
182
183 The source code of the \Dune{} framework can be
184 downloaded from the web page.  To get started, it is easiest to
186 \texttt{dune-common}, \texttt{dune-grid} and \texttt{dune-grid-howto}.
188 %
189 \begin{center}
191 \end{center}
192 %
193
195 anonymous Git.  For further information, please see the web page.
196
197 \section{Installation}
198
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}
204
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.
212
213 In order to build the \Dune{} framework, you need a standards
214 compliant C++ compiler.  We tested compiling with GNU \texttt{g++} in
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.
219
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
229 \Dune{} modules in \texttt{dune-home} with a basic configuration.
230
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}
257
258 For information on how to build and configure the respective grids,
259 please see the \Dune{} web page.
260
261 \section{Code documentation}
262
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)
269 \begin{center}
270 \url{http://www.dune-project.org/doxygen/master/index.html}
271 \end{center}
272
273 \section{Licence}
274
275 The \Dune{} library and headers are licensed under version 2 of the
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.
282
283 The exact wording of the exception reads as follows:
284
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
296 \end{quote}
297
298 %\section{Contributing to DUNE}
299
300
301
302
303 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
304 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
305 \chapter{Getting started}
306 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
307 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
308
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.
312
314
315 Let us start with a replacement of the famous \emph{hello world}
316 program given below.
317
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}
322
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
334 functions which print useful information about a grid.
335
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.
346
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.
358
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:
362
363 \begin{lst}[Output of gettingstarted] \mbox{}
364
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}
372
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}.
391
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}
397
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}
405
406 \section{Traversing a grid --- A first look at the grid interface}
407
408 After looking at very first simple example we are now ready to go on
409 to a more complicated one. Here it is:
410
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}
416
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}.
425
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.
431
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.
436
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
446 different kinds of entities and which also provide different patterns
447 of access.
448
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.
453
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.
462
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.
473
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}
486
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.
491
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.
506
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).
522
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.
526
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.
538
539 The following listing shows the output of the program.
540
541 \begin{lst}[Output of traversal] \mbox{}
542
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)
550
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)
562
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
566
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}
574
575 \begin{rem} Define the end iterator for efficiency.
576 \end{rem}
577
578 \begin{exc} Play with different dimensions, codimension
579   (\lstinline!YaspGrid! supports all codimenions) and refinements.
580 \end{exc}
581
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}
586
587
588 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
589 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
590 \chapter{The \texorpdfstring{\Dune{}}{DUNE} grid interface}
591 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
592 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
593
594
595 \section{Grid definition}
596
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.
603
604 \minisec{Reference elements}
605
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.
612
613 \minisec{Dimension and world dimension}
614
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.
618
619 \minisec{Faces, entities and codimension}
620
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$.
631
632 \minisec{Conformity}
633
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}.
639
640 \minisec{Element types}
641
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.
646
647 \minisec{Hierarchically nested grids, macro grid}
648
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.
662
663 \minisec{Leaf grid}
664
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.
677
678 \minisec{Refinement rules}
679
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.
693
694 \minisec{Summary}
695
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}
709
710
711
712 \section{Concepts}
713
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.
722
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.
730
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.
740
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.
761
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}
771
772
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}
781
782
783 \subsection{Common types}
784
785 Some types in the grid interface do not depend on a specific model,
786 i.~e.~they are shared by all implementations.
787
788 \minisec{\href{https://www.dune-project.org/doxygen/master/group__GeometryReferenceElements.html}{Dune::ReferenceElement}}
789
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.
793
794 \minisec{\href{http://www.dune-project.org/doxygen/master/classDune_1_1GeometryType.html}{Dune::GeometryType}}
795
796 defines names for the reference elements.
797
798 \minisec{\href{http://www.dune-project.org/doxygen/master/classDune_1_1CollectiveCommunication.html}{Dune::CollectiveCommunication}}
799
800 defines an interface to global communication operations in a portable
801 and transparent way. In particular also for sequential grids.
802
803
804 \subsection{Concepts of the \texorpdfstring{\Dune{}}{DUNE} grid interface}
805
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).
810
811 \minisec{\href{http://www.dune-project.org/doxygen/master/classDune_1_1Grid.html}{Grid}}
812
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.
819
820 \minisec{\href{http://www.dune-project.org/doxygen/master/classDune_1_1GridView.html}{GridView}}
821
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.
827
828 \minisec{\href{http://www.dune-project.org/doxygen/master/classDune_1_1Entity.html}{Entity}}
829
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.
835
836 \minisec{\href{http://www.dune-project.org/doxygen/master/classDune_1_1Geometry.html}{Geometry}}
837
838     Geometry
839     encapsulates the geometric part of an entity by mapping local
840     coordinates in a reference element to world coordinates.
841
842 \minisec{\href{http://www.dune-project.org/doxygen/master/classDune_1_1EntityPointer.html}{EntityPointer}}
843
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.
848
849 \minisec{\href{http://www.dune-project.org/doxygen/master/classDune_1_1EntityIterator.html}{EntityIterator}}
850
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.
856
857 \minisec{\href{http://www.dune-project.org/doxygen/master/classDune_1_1IntersectionIterator.html}{IntersectionIterator}}
858
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.
871
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}}
873
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}).
880
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}}
882
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}).
890
891
892 \section{Propagation of type information}
893
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.
905
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.
920
921 \chapter{Constructing grid objects}
922
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.
926
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.
936
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.
944
945 \section{Creating Structured Grids}
946
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.
953
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);
959
960     Dune::FieldVector<double,dim> upper(1.0);
961
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!:
967
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$:
981
982 \begin{lstlisting}
983  std::array<std::vector<ctype,2> > coords;
984  coords[0] = {-1., 1.};
985  coords[1] = {-1., 1.};
986
987  YaspGrid<dim, TensorProductCoordinates<ctype,dim> > grid(coords);
988 \end{lstlisting}
989
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.
998
999 \section{Reading Unstructured Grids from Files}
1000
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.
1006
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.
1023 and write
1024 \begin{lstlisting}
1026 \end{lstlisting}
1027
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%
1031 installed.  Then
1032 \begin{lstlisting}
1034 \end{lstlisting}
1035 reads the grid in \lstinline!filename! into the FooGrid.
1036
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!.
1044
1045
1046 \section{The \texorpdfstring{\Dune}{Dune} Grid Format (DGF)}
1047
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}}.
1052
1053 Here we only give a short introduction.
1054 Basically one can choose the grid manager during the make process of
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 \\
1065 \lstinline!ALBERTAGRID,ALUGRID_CUBE,ALUGRID_SIMPLEX,ALUGRID_CONFORM,ONEDGRID,UGGRID!,
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.
1092
1093 \section{The Grid Factory}
1094 \label{sec:grid_factory}
1095
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.
1105
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:
1112
1113 \begin{enumerate}
1114 \item Create a GridFactory Object
1115
1116 Get a new GridFactory object by calling
1117
1118 \begin{lstlisting}[basicstyle=\small]
1119     GridFactory< FooGrid > factory;
1120 \end{lstlisting}
1121
1122 \item Enter the Vertices
1123
1124 Insert the grid vertices by calling
1125
1126 \begin{lstlisting}[basicstyle=\small]
1127     factory.insertVertex(const FieldVector<ctype,dimworld>& position);
1128 \end{lstlisting}
1129
1130 for each vertex. The order of insertion determines the level- and leaf indices of your level 0 vertices.
1131
1132 \item Enter the elements
1133
1134 For each element call
1135
1136 \begin{lstlisting}[basicstyle=\small]
1137     factory.insertElement(Dune::GeometryType type,
1138         const std::vector<int>& vertices);
1139 \end{lstlisting}
1140
1141 The parameters are
1142
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}
1148
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.
1151
1152 \item Parametrized Domains
1153
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.
1155
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
1158
1159 \begin{lstlisting}[basicstyle=\small]
1160     template <int dim, int dimworld> Dune::BoundarySegment
1161 \end{lstlisting}
1162
1163 This is an abstract base class which requires you to overload the method
1164
1165 \begin{lstlisting}[basicstyle=\small]
1166     virtual FieldVector< double, dimworld >
1167          operator() (const FieldVector< double, dim-1 > &local)
1168 \end{lstlisting}
1169
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
1171
1172 \begin{lstlisting}[basicstyle=\small]
1173     factory.insertBoundarySegment(const std::vector<int>& vertices,
1174                                   const BoundarySegment<dim,dimworld> *
1175                                      boundarySegment = NULL);
1176 \end{lstlisting}
1177
1178 Control over the allocated objects is taken from you, and the grid object will take care of their destruction.
1179
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}).
1185
1186 \item Finish construction
1187
1188 To finish off the construction of the FooGrid object call
1189
1190 \begin{lstlisting}[basicstyle=\small]
1191     FooGrid* grid = factory.createGrid();
1192 \end{lstlisting}
1193
1194 This time it is you who gets full responsibility for the allocated object.
1195 \end{enumerate}
1196
1197 \section{Attaching Data to a New Grid}
1198 \label{sec:import_data_for_new_grids}
1199
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.
1205
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.
1218
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.
1225
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.
1239
1240
1241 \section{Example: The \texorpdfstring{\lstinline{UnitCube}}{UnitCube} class}
1242
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$.
1246
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.
1257
1258 The definition of the general template is as follows.
1259
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}
1265
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.
1274
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.
1282
1283 For historic reasons, there are also specializations for
1284 \lstinline!ALUGrid! and \lstinline!AlbertaGrid!.
1285
1286 \minisec{YaspGrid}
1287
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.
1293
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}
1299
1300 \minisec{AlbertaGrid}
1301
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!.
1313
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}
1319
1320 \minisec{ALUGrid}
1321
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.
1332
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}
1338
1339
1340
1341 %\section{Capabilities}
1342
1343 % \chapter{Reference elements}
1344
1345
1348
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.
1352
1353 \section{Numerical integration}
1354
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}.$$
1359
1360
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.
1368
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$.
1380
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.
1394
1395 In the following we show how the formula (\ref{Eq:IntegrateDomain})
1396 can be realised within \Dune.
1397
1398
1399 \section{Functors}
1400
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.
1406
1407
1408 \begin{lst}[dune-grid-howto/functors.hh] Here are some examples for functors.
1409
1410 \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
1411 numberstyle=\tiny, numbersep=5pt]{../functors.hh}
1412 \end{lst}
1413
1414
1415 \section{Integration over a single element}
1416
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.
1421
1422 \begin{lst}[dune-grid-howto/integrateentity.hh] \mbox{}
1423
1424 \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
1425 numberstyle=\tiny, numbersep=5pt]{../integrateentity.hh}
1426 \end{lst}
1427
1428 Line \ref{ieh:qr} extracts a reference to a
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.
1435
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$.
1441
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}).
1447
1448 \section{Integration with global error estimation}
1449
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.
1467
1468 \begin{lst}[dune-grid-howto/integration.cc] \mbox{}
1469
1470 \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
1471 numberstyle=\tiny, numbersep=5pt]{../integration.cc}
1472 \end{lst}
1473
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:
1477
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}
1492
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$.
1496
1497
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}
1502
1503 \begin{exc} Try different grid implementations and dimensions and
1504   compare the run-time.
1505 \end{exc}
1506
1507 \begin{exc} Try different integrands $f$ and look at the development
1508   of the (estimated) error in the integral.
1509 \end{exc}
1510
1511 \chapter{Attaching user data to a grid}
1512
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.
1520
1521 \section{Mappers}\label{ch:mappers}
1522
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.
1527
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}
1540
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.
1545
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}
1567
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.
1571
1572 \minisec{Index based mappers}
1573
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.
1582
1583 \minisec{Id based mappers}
1584
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$.
1594
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.
1599
1600 \section{Visualization of discrete functions}
1601
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.
1606
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:
1609
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}
1615
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.
1627
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.
1636
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.
1646
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.
1656
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.
1660
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}
1666
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]!.
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}
1682
1683 Finally the following listing shows the main program that calls the
1684 two functions just discussed:
1685
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}
1691
1692 \section{Cell centered finite volumes}
1693 \label{Sec:CellCenteredFV}
1694
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.
1698
1699 We will solve the linear hyperbolic PDE
1700
1701 \frac{\partial c}{\partial t} + \nabla\cdot (uc) = 0 \ \ \text{ in
1702   $\Omega\times T$}
1703 \label{Eq:TransportEquation}
1704
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
1712 c(x,0) = c_0(x) \ \ x\in\Omega
1713
1714 and the boundary condition
1715
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
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$.
1721
1722 \subsection{Numerical Scheme}
1723
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.
1727
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
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
1740 Using integration by parts we arrive at
1741
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
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
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
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
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
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
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
1782 Inserting this into (\ref{Eq:DiscreteForm}) and solving for
1783 $C_i^{n+1}$ we obtain
1784
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
1792 One can show that the scheme is stable provided the following condition holds:
1793
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
1802 When we rewrite \ref{Eq:IterativeForm1} in the form
1803
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
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.
1815
1816
1817 \subsection{Implementation}
1818
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.
1822
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}
1828
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.
1836
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}
1842
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}
1854
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}
1867
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.
1873
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
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.
1893
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.
1898
1899 Finally, line \ref{evh:updc} computes the new concentration by adding
1900 the scaled update to the current concentration.
1901
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.
1905
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}
1911
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.
1919
1920 Finally, the main program:
1921
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}
1927
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}.
1935
1936 \section{A FEM example: The Poisson equation}
1937 \label{Sec:FEMPoisson}
1938
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}
1944
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.
1946
1947 The starting point of the Finite Element Method is the variational formulation of \ref{eq:fempoisson}, which is obtained by partial integration:
1948
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
1952
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}
1962
1963 As we want to use linear finite elements we choose our test function space to be
1964
1965  V_h = \{ u\in C(\bar{\Omega} )\ \Big| \ u\big|_{\Delta}\in\mathcal{P}_1 (\Delta ) \ \forall\Delta\in\mathcal{T}\}
1966
1967
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.
1969
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}
1975
1976 to evaluate the basis functions and their gradients.
1977
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
1980  A_{ij}=a(\phi_i ,\phi_j ) \ \ b_i=\ell (\phi_i )
1981
1982
1983 The coefficients $u_1,...,u_N$ are then obtained by solving $Au=b$.
1984
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$
1986
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}
1995
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
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
2003
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
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
2009
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.
2011
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}.
2013
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}
2041
2042 This is possible as---using the nodal basis---the coefficients match the value of the numerical solution at the corresponding node.
2043
2044 \subsection{Implementation}
2045
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!
2047
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:
2049
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}
2055
2056 And now the actual FEM code:
2057
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}
2063
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.
2065
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}.
2067
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}.
2069
2070 When you visualize your results, you should get something like figure \ref{Fig:FEM1} or \ref{Fig:FEM2}!
2071
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}
2078
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}
2085
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}
2092
2093 \begin{exc}
2094  Modify the code in order to make it handle Neumann boundary conditions too!
2095 \end{exc}
2096
2097
2098
2100
2102
2104
2105 In this section we describe briefly the adaptive multigrid integration
2106 algorithm presented in \cite{Deuflhard93}.
2107
2108 \minisec{Global error estimation}
2109
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}).
2113
2114 \minisec{Local error estimation}
2115
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
2121 \bar{\epsilon}(\omega) = |I_f^p(\omega)-I_f^q(\omega)|
2122
2123 is an estimator for the local error on the element $\omega$.
2124
2125 \minisec{Refinement strategy}
2126
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
2144 \epsilon^+(\omega) = \frac{\bar{\epsilon}(\omega)^2}{\bar{\epsilon}(\omega^-)}
2145
2146 is an estimate for the local error on $\omega^+$,
2147 i.~e.~$\bar{\epsilon}(\omega^+)$.
2148
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.
2153
2154 By combining the two approaches we get the threshold value $\kappa$
2155 actually used in the code:
2156
2157 \kappa = \min\left(\max\limits_{\omega} \epsilon^+(\omega), \frac12
2158   \max\limits_{\omega} \bar{\epsilon}(\omega) \right).
2159
2160
2161 \minisec{Algorithm}
2162
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}
2175
2176 \subsection{Implementation of the algorithm}
2177
2178 The algorithm above is realized in the following code.
2179
2181 \nopagebreak
2182 \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
2184 \end{lst}
2185
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
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
2203
2204 \begin{figure}
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.}
2211 \end{figure}
2212
2213 Figure \ref{Fig:AdaptiveIntegration} shows two grids generated by the
2215
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}
2220
2221
2222
2223 \section{Adaptive cell centered finite volumes}
2224
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}.
2232
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}
2249
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.
2255
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}
2269
2270 The complete mesh adaptation is done by the function
2271 \lstinline!finitevolumeadapt! in the following listing:
2272
2274 \nopagebreak
2275 \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
2277 \end{lst}
2278
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.
2290
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}.
2299
2300 Here is the new main program with an adapted \lstinline!timeloop!:
2301
2303 \nopagebreak
2304 \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
2306 \end{lst}
2307
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.
2319
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}
2324
2325 \chapter{Parallelism}
2326
2327 \section{\texorpdfstring{\Dune{}}{DUNE} Data Decomposition Model}
2328
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.
2337
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}
2354
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$.
2368
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}
2390
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$}
2397
2398  \put(12,31){$c=0$}
2399  \put(45,31){$c=1$}
2400  \put(80,31){$c=2$}
2401
2402  \put(12,1){$c=0$}
2403  \put(45,1){$c=1$}
2404  \put(80,1){$c=2$}
2405 \end{overpic}
2406
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}
2415
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.
2434
2435
2436 \section{Communication Interfaces}
2437
2438 This section explains how the exchange of data between the partitions
2439 in different processes is organized in a flexible and portable way.
2440
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$.
2457
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:
2467
2468
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}
2474
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.
2498
2499 \section{Parallel finite volume scheme}
2500
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}.
2505
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}
2511
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}.
2518
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}.
2528
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.
2537
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.
2545
2546 Finally, we need a new main program, which is in the following listing:
2547
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}
2553
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
2566
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.
2573
2574 % \chapter{Input and Output}
2575
2576 % \section{Visualization with Grape}
2577
2578
2579
2580 % \section{Visualization with VTK}
2581
2582 % \section{Dune portable grid format}
2583
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
2588
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}}}
2591
2592 % in the dune-grid documentation.
2593
2594 % \section{Amiramesh output}
2595
2596
2597 % \chapter{Iterative Solver Template Library}
2598
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}
2604
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}
2610
2612 % \nopagebreak
2613 % \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left,
2615 % \end{lst}
2616
2617
2618 % \chapter{Outlook}
2619
2620
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}
2626
2627 \includegraphics[width=0.32\textwidth]{img/alberta3d}\
2628 \includegraphics[width=0.32\textwidth]{img/alusimplex3d}\
2629 \includegraphics[width=0.32\textwidth]{img/alucube3d}
2630
2631
2632 \includegraphics[width=0.32\textwidth]{img/ugsimplex3d}\
2633 \includegraphics[width=0.32\textwidth]{img/ugcube3d}\
2634 \includegraphics[width=0.32\textwidth]{img/iso}
2635
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}
2647
2648
2649
2650 \bibliographystyle{plain}
2651 \bibliography{grid-howto.bib}
2652
2653 \end{document}
`