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

1 \documentclass[11pt,a4paper,headinclude,footinclude,DIV16,headings=normal]{scrreprt} 2 \usepackage[automark]{scrpage2} 3 \usepackage{scrhack} 4 \usepackage[ansinew]{inputenc} 5 \usepackage{amsmath} 6 \usepackage{amsfonts} 7 \usepackage{theorem} 8 \usepackage{xspace} 9 \usepackage{color} 10 \usepackage{overpic} 11 \usepackage{listings} 12 \lstset{language=C++, basicstyle=\ttfamily, 13 keywordstyle=\color{black}\bfseries, tabsize=4, stringstyle=\ttfamily, 14 commentstyle=\itshape, extendedchars=true, escapeinside={/*@}{@*/}} 15 \usepackage{hyperref} 16 \usepackage{makeidx} 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} 31 \theoremheaderfont{\sffamily\bfseries} 32 \newtheorem{exc}{Exercise}[chapter] 33 \newtheorem{rem}[exc]{Remark} 34 \newtheorem{lst}{Listing} 35 \newtheorem{warn}[exc]{Warning} 36 37 \pagestyle{scrheadings} 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 59 \includegraphics[width=0.32\textwidth]{img/adaptiveintegration_alberta2d}\hfill 60 \includegraphics[width=0.32\textwidth]{img/alucube3d} 61 %\includegraphics[width=0.32\textwidth]{img/ug2dquad-view2} 62 \\ 63 \vspace{10mm} 64 {\normalsize $^\ast$Abteilung `Simulation gro\ss er Systeme', 65 Universit\"at Stuttgart,\\ 66 Universit\"atsstr.~38, D-70569 Stuttgart, Germany}\\ 67 % 68 \bigskip 69 {\normalsize $^\dagger$Abteilung f\"ur Angewandte Mathematik, Universit\"at Freiburg,\\ 70 Hermann-Herder-Str.~10, D-79104 Freiburg, Germany}\\ 71 % 72 \bigskip 73 {\normalsize $^\mathparagraph$Institut f\"ur Numerische und Angewandte Mathematik, Universit\"at M\"unster,\\ 74 Einsteinstr.~62, D-48149 M\"unster, Germany}\\ 75 % 76 \bigskip 77 {\normalsize $^\ddagger$Institut f\"ur Mathematik II,\\ Freie Universit\"at Berlin, 78 Arnimallee 6, D-14195 Berlin, Germany}\\ 79 % 80 \bigskip 81 {\normalsize \texttt{\url{http://www.dune-project.org/}}}\\ 82 } 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 152 quadrature rules. 153 \item \texttt{dune-grid} is the grid interface and is covered in 154 this document. It defines nonconforming, hierarchically nested, 155 multi-element-type, parallel grids in arbitrary space dimensions. 156 Graphical output with several packages is available, e.~g.~file 157 output to VTK (parallel XML format for 158 unstructured grids). The graphics package Grape, \cite{Grape} has 159 been integrated in interactive mode. 160 \item \texttt{dune-istl} -- \textit{Iterative Solver Template 161 Library.} Provides generic sparse matrix/vector classes and a 162 variety of solvers based on these classes. A special feature is the 163 use of templates to exploit the recursive block structure of finite 164 element matrices at compile time. Available solvers include Krylov 165 methods, (block-) incomplete decompositions and aggregation-based 166 algebraic multigrid. 167 \item \texttt{dune-localfunctions} -- \textit{Library of local base 168 functions.} Provides classes for base functions on reference elements 169 from which global discrete function spaces can be constructed. 170 \end{itemize} 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 181 \section{Download} 182 183 The source code of the \Dune{} framework can be 184 downloaded from the web page. To get started, it is easiest to 185 download the latest stable version of the tarballs of 186 \texttt{dune-common}, \texttt{dune-grid} and \texttt{dune-grid-howto}. 187 These are available on the \Dune{} download page: 188 % 189 \begin{center} 190 \url{http://www.dune-project.org/download.html} 191 \end{center} 192 % 193 194 Alternatively, you can download the latest development version via 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 215 version 4.4 or newer and \texttt{clang} in version 3.4 or newer. Recent versions 216 of Intel \texttt{icc} ($\geq 15$) should work, too. For older versions of 217 \texttt{icc}, like $14.0.3$, several issues are known and usage of these 218 is discouraged. 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 228 you downloaded if necessary. This should configure and build all 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) 268 or its latest version is available at 269 \begin{center} 270 \url{http://www.dune-project.org/doxygen/master/index.html} 271 \end{center} 272 273 \section{Licence} 274 275 The \Dune{} library and headers are licensed under version 2 of the 276 GNU General Public License% 277 \footnote{\url{http://www.gnu.org/licenses/gpl-2.0.html}}, with a special 278 exception for linking and compiling against \Dune{}, the so-called 279 ``runtime exception.'' The license is intended to be similar to the 280 GNU Lesser General Public License, which by itself isn't suitable for 281 a C++ template library. 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 295 Public License. 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 313 \section{Creating your first grid} 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 333 \lstinline!dune/grid/common/gridinfo.hh! loads the headers of some 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 445 knows several different iterators which provide access to the 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 859 IntersectionIterator provides access to all entities of 860 codimension 0 that have an intersection of codimension 1 with a 861 given entity of codimension 0. In a conforming mesh these are the 862 face neighbors of an element. For two entities with a common 863 intersection the IntersectionIterator can be dereferenced as an Intersection 864 object which in turn provides information about the geometric location of 865 the intersection. Furthermore this Intersection class also provides 866 information about intersections of an entity 867 with 868 the internal or external boundaries. 869 The IntersectionIterator provides intersections between 870 codimension 0 entities among the same GridView. 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 1001 Unstructured grids usually require much more information than what can 1002 reasonable be provided within the program code. Instead, they are usually 1003 read from special files. A large variety of different file formats for 1004 finite element grids exists, and \Dune provides support for some of them. 1005 Again, no interface specification exists for file readers in \Dune. 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. 1022 To read such a file into a FooGrid, include \lstinline!dune/grid/io/file/gmshreader.hh! 1023 and write 1024 \begin{lstlisting} 1025 FooGrid* grid = GmshReader<GridType>::read(filename); 1026 \end{lstlisting} 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% 1030 \footnote{\url{http://amira.com/downloads/patch-archive/patches412/81.html}} 1031 installed. Then 1032 \begin{lstlisting} 1033 FooGrid* grid = AmiraMeshReader<GridType>::read(filename); 1034 \end{lstlisting} 1035 reads the grid in \lstinline!filename! into the FooGrid. 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 1055 your program: 1056 \lstinline!make GRIDTYPE=MYGRID GRIDDIM=d GRIDWORLD=w! 1057 Including \lstinline!config.h! will then 1058 introduce the namespace \lstinline!GridSelector! into the 1059 \lstinline!Dune! namespace. This contains the typedef 1060 \lstinline!GridType! which is the type of the grid. Furthermore the 1061 required header files for the grid manager are included. 1062 Through the module \Dune{}-grid the following grid managers can be used 1063 (for \lstinline!MYGROD! in the example above): 1064 \\ 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 1346 \chapter{Quadrature rules} 1347 \label{sec:quadrature} 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 1429 \lstinline!Dune::QuadratureRule! from the 1430 \lstinline!Dune::QuadratureRules! singleton which is a container 1431 containing quadrature rules for all the different reference element 1432 types and different orders of approximation. A rule of order $p$ is 1433 guaranteed to integrate polynomials of order $p$ exactly, up to 1434 floating point errors. 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]!. 1676 \item Use \lstinline!addVertexData! instead of \lstinline!addCellData! 1677 on the \lstinline!Dune::VTKWriter!. 1678 \item Pass \lstinline!polynomialOrder=1! instead of 0 as 1679 the second last argument of \lstinline!grape.displayVector!. This argument is the polynomial 1680 degree of the approximation. 1681 \end{itemize} 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 \begin{equation} 1701 \frac{\partial c}{\partial t} + \nabla\cdot (uc) = 0 \ \ \text{ in 1702 $\Omega\times T$} 1703 \label{Eq:TransportEquation} 1704 \end{equation} 1705 where $\Omega\subset\mathbb{R}^d$ is a domain, $T=(0,t_{\text{end}})$ 1706 is a time interval, 1707 $c:\Omega\times T\to\mathbb{R}$ is the unknown concentration and 1708 $u:\Omega\times T\to\mathbb{R}^d$ is a given velocity field. We 1709 require that the velocity field is divergence free for all times. 1710 The equation is subject to the initial condition 1711 \begin{equation} 1712 c(x,0) = c_0(x) \ \ x\in\Omega 1713 \end{equation} 1714 and the boundary condition 1715 \begin{equation} 1716 c(x,t) = b(x,t) \ \ t>0, x\in\Gamma_{\text{in}}(t)=\{y\in\partial\Omega\ 1717 |\ u(y,t)\cdot\nu(y)<0\}. 1718 \end{equation} 1719 Here $\nu(x)$ is the unit outer normal at a point $y\in\partial\Omega$ 1720 and $\Gamma_{\text{in}}(t)$ is the inflow boundary at time $t$. 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 \begin{equation} 1734 \int\limits_{\omega_i}\int\limits_{t_n}^{t_{n+1}}\frac{\partial 1735 c}{\partial t}\ dt\,dx + 1736 \int\limits_{\omega_i}\int\limits_{t_n}^{t_{n+1}} \nabla\cdot (uc) \ 1737 dt\,dx = 0 \ \ \forall i. 1738 \label{Eq:TransportEquationIntegrated} 1739 \end{equation} 1740 Using integration by parts we arrive at 1741 \begin{equation} 1742 \int\limits_{\omega_i} c(x,t_{n+1})\ dx - \int\limits_{\omega_i} 1743 c(x,t_{n})\ dx + \int\limits_{t_n}^{t_{n+1}} 1744 \int\limits_{\partial\omega_i} c u\cdot\nu\ ds\,dt = 0 \ \ \forall i. 1745 \end{equation} 1746 Now we approximate $c$ by a cell-wise constant function $C$, where 1747 $C_i^n$ denotes the value in cell $\omega_i$ at time $t_n$. 1748 Moreover we subdivide the boundary $\partial\omega_i$ into 1749 facets $\gamma_{ij}$ which are either intersections with 1750 other cells $\partial\omega_i\cap\partial\omega_j$, 1751 or intersections with the boundary 1752 $\partial\omega_i\cap\partial\Omega$. 1753 Evaluation of the fluxes 1754 at time level $t_n$ leads to the following equation for the unknown 1755 cell values at $t_{n+1}$: 1756 \begin{equation} 1757 C_i^{n+1}|\omega_i| - C_i^{n}|\omega_i| + \sum\limits_{\gamma_{ij}} 1758 \phi(C_i^n, C_j^n, u_{ij}^n\cdot\nu_{ij}; \gamma_{ij}, t_n ) 1759 |\gamma_{ij}| \Delta t^n = 0 \ \ \forall i, 1760 \label{Eq:DiscreteForm} 1761 \end{equation} 1762 where $\Delta t^n=t_{n+1}-t_n$, $u_{ij}^n$ is the velocity on the facet $\gamma_{ij}$ at time 1763 $t_n$, $\nu_{ij}$ is the unit outer normal of the facet $\gamma_{ij}$ and $\phi$ is the flux function defined as 1764 \begin{equation} 1765 \phi(C_i^n, C_j^n, u_{ij}^n\cdot\nu_{ij}; \gamma_{ij}, t_n ) = 1766 \left\{\begin{array}{ll} 1767 b(\gamma_{ij})\, u_{ij}^n\cdot\nu_{ij} & \gamma_{ij} \subset \Gamma_{\text{in}}(t) \\ 1768 C_j^n\, u_{ij}^n\cdot\nu_{ij}& 1769 \gamma_{ij}=\partial\omega_i\cap\partial\omega_j\,\wedge\,u_{ij}^n\cdot\nu_{ij} < 0\\ 1770 C_i^n\, u_{ij}^n\cdot\nu_{ij}& u_{ij}^n\cdot\nu_{ij} \geq 0 1771 \end{array}\right.. 1772 \end{equation} 1773 Here $b(\gamma_{ij})$ denotes evaluation of the boundary condition on 1774 an inflow facet $\gamma_{ij}$. If we formally set 1775 $C_j^n=b(\gamma_{ij})$ on an inflow facet $\gamma_{ij} \subset 1776 \Gamma_{\text{in}}(t)$ we can derive the following shorthand notation 1777 for the flux function: 1778 \begin{equation} 1779 \phi(C_i^n, C_j^n, u_{ij}^n\cdot\nu_{ij}; \gamma_{ij}, t_n ) = C_i^n 1780 \max(0,u_{ij}^n\cdot\nu_{ij}) - C_j^n \max(0,-u_{ij}^n\cdot\nu_{ij}). 1781 \end{equation} 1782 Inserting this into (\ref{Eq:DiscreteForm}) and solving for 1783 $C_i^{n+1}$ we obtain 1784 \begin{equation} 1785 C_i^{n+1} = C_i^n\left(1-\Delta t^n \sum\limits_{\gamma_{ij}} 1786 \frac{|\gamma_{ij}|}{|\omega_i|} 1787 \max(0,u_{ij}^n\cdot\nu_{ij})\right) + \Delta t^n 1788 \sum\limits_{\gamma_{ij}} C_j^n \frac{|\gamma_{ij}|}{|\omega_i|} 1789 \max(0,-u_{ij}^n\cdot\nu_{ij}) \ \ \forall i. 1790 \label{Eq:IterativeForm1} 1791 \end{equation} 1792 One can show that the scheme is stable provided the following condition holds: 1793 \begin{equation} 1794 \forall i:\ 1-\Delta t^n \sum\limits_{\gamma_{ij}} 1795 \frac{|\gamma_{ij}|}{|\omega_i|} 1796 \max(0,u_{ij}^n\cdot\nu_{ij})\geq 0 \ \Leftrightarrow\ \Delta 1797 t^n\leq \min_i \left(\sum\limits_{\gamma_{ij}} 1798 \frac{|\gamma_{ij}|}{|\omega_i|} 1799 \max(0,u_{ij}^n\cdot\nu_{ij})\right)^{-1}. 1800 \label{Eq:TimeStepControl} 1801 \end{equation} 1802 When we rewrite \ref{Eq:IterativeForm1} in the form 1803 \begin{equation} 1804 C_i^{n+1} = C_i^n - \Delta t^n 1805 \underbrace{\sum\limits_{\gamma_{ij}} \frac{|\gamma_{ij}|}{|\omega_i|} 1806 \left( C_i^n \max(0,u_{ij}^n\cdot\nu_{ij}) + 1807 C_j^n \max(0,-u_{ij}^n\cdot\nu_{ij})\right)}_{\delta_i} \ \ \forall i 1808 \label{Eq:IterativeForm2} 1809 \end{equation} 1810 then it becomes clear that the optimum time step $\Delta t^n$ and the 1811 update $\delta_i$ for each cell can be computed in a single iteration 1812 over the grid. The computation $C^{n+1} = C^n - \Delta t^n \delta$ can then 1813 be realized with a simple vector update. In this form, the algorithm 1814 can also be parallelized in a straightforward way. 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 1883 \lstinline!boundary()! returns true (line \ref{evh:bndry}), see also 1884 Figure~\ref{Fig:IsIt}. An intersection $\gamma_{ij}$ is 1885 described by several mappings: (i) from a reference element of the 1886 intersection (with a dimension equal to the grid's dimension minus 1) 1887 to the reference elements of the two elements $\omega_i$ and 1888 $\omega_j$ and (ii) from a reference element of the intersection to 1889 the global coordinate system (with the world dimension). If an 1890 intersection is with another element then the \lstinline!outside()! 1891 method returns an \lstinline!EntityPointer! to an entity of 1892 codimension 0. 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 \begin{equation} 1949 \underbrace{\int\limits_{\Omega}\nabla\ u\cdot\nabla v dx}_{ =:a(u,v)} 1950 = \underbrace{\int\limits_{\Omega} fv dx}_{=:\ell (v)} \ \ v\in V_h 1951 \end{equation} 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 \begin{equation} 1965 V_h = \{ u\in C(\bar{\Omega} )\ \Big| \ u\big|_{\Delta}\in\mathcal{P}_1 (\Delta ) \ \forall\Delta\in\mathcal{T}\} 1966 \end{equation} 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 \begin{equation} 1980 A_{ij}=a(\phi_i ,\phi_j ) \ \ b_i=\ell (\phi_i ) 1981 \end{equation} 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 \begin{equation} 1998 \label{equ:computea} 1999 \Rightarrow \ A_{ij} = \sum\limits_{\Delta\in\tau} \sum\limits_{k=1}^{p_1} 2000 \omega_k(J_g^{-T}\hat{\nabla}\hat{\phi_i} )(\xi_k ) 2001 \cdot (J_g^{-T}\hat{\nabla}\hat{\phi_j} )(\xi_k )\mu (\xi_k ) 2002 \end{equation} 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 \begin{equation} 2006 \label{equ:computeb} 2007 b_i = \sum\limits_{\Delta\in\tau}\sum\limits_{k=1}^{p_2}\omega_k f(g(\xi_k ))\phi_i (g(\xi_k ))\mu(\xi_k ) 2008 \end{equation} 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 2099 \chapter{Adaptivity} 2100 2101 \section{Adaptive integration} 2102 2103 \subsection{Adaptive multigrid integration} 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 \begin{equation} 2121 \bar{\epsilon}(\omega) = |I_f^p(\omega)-I_f^q(\omega)| 2122 \end{equation} 2123 is an estimator for the local error on the element $\omega$. 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 \begin{equation} 2144 \epsilon^+(\omega) = \frac{\bar{\epsilon}(\omega)^2}{\bar{\epsilon}(\omega^-)} 2145 \end{equation} 2146 is an estimate for the local error on $\omega^+$, 2147 i.~e.~$\bar{\epsilon}(\omega^+)$. 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 \begin{equation} 2157 \kappa = \min\left(\max\limits_{\omega} \epsilon^+(\omega), \frac12 2158 \max\limits_{\omega} \bar{\epsilon}(\omega) \right). 2159 \end{equation} 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 2180 \begin{lst}[File dune-grid-howto/adaptiveintegration.cc] \mbox{} 2181 \nopagebreak 2182 \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left, 2183 numberstyle=\tiny, numbersep=5pt]{../adaptiveintegration.cc} 2184 \end{lst} 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 2197 \lstinline!preAdapt()!, \lstinline!adapt()! and 2198 \lstinline!postAdapt()! will be explained with the next example. Note 2199 the flexibility of this algorithm: It runs in any space dimension on 2200 any kind of grid and different integration orders can easily be 2201 incorporated. And that with just about 100 lines of code including 2202 comments. 2203 2204 \begin{figure} 2205 \includegraphics[width=0.48\textwidth]{img/adaptiveintegration_alberta2d}\hfill 2206 \includegraphics[width=0.48\textwidth]{img/adaptiveintegration_ug3d} 2207 \caption{Two and three-dimensional grids generated by the adaptive 2208 integration algorithm applied to the needle pulse. Left grid is 2209 generated using Alberta, right grid is generated using UG.} 2210 \label{Fig:AdaptiveIntegration} 2211 \end{figure} 2212 2213 Figure \ref{Fig:AdaptiveIntegration} shows two grids generated by the 2214 adaptive integration algorithm. 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 2273 \begin{lst}[File dune-grid-howto/finitevolumeadapt.hh] \mbox{} 2274 \nopagebreak 2275 \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left, 2276 numberstyle=\tiny, numbersep=5pt]{../finitevolumeadapt.hh} 2277 \end{lst} 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 2302 \begin{lst}[File dune-grid-howto/adativefinitevolume.cc] \mbox{} 2303 \nopagebreak 2304 \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left, 2305 numberstyle=\tiny, numbersep=5pt]{../adaptivefinitevolume.cc} 2306 \end{lst} 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 2565 communicaton overhead increases. 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 2611 % \begin{lst}[File dune-grid-howto/groundwateradapt.hh] \mbox{} 2612 % \nopagebreak 2613 % \lstinputlisting[basicstyle=\ttfamily\scriptsize,numbers=left, 2614 % numberstyle=\tiny, numbersep=5pt]{../groundwateradapt.hh} 2615 % \end{lst} 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}