Reference documentation for deal.II version Git 3fdf7d9 2017-04-24 17:39:03 -0400
|Table of contents|
Since this is the first tutorial program, let us comment first on how this tutorial and the rest of the deal.II documentation is supposed to work. The documentation for deal.II comes essentially at three different levels:
Let's come back to the tutorial, since you are looking at the first program (or "step") of it. Each tutorial program is subdivided into the following sections:
The tutorials are not only meant to be static documentation, but you should play with them. To this end, go to the
examples/step-1 directory (or whatever the number of the tutorial is that you're interested in) and type
The first command sets up the files that describe which include files this tutorial program depends on, how to compile it and how to run it. This command should find the installed deal.II libraries as well that were generated when you compiled and installed everything as described in the deal.II ReadMe file. If this command should fail to find the deal.II library, then you need to provide the path to the installation using the command
The second of the commands above compiles the sources into an executable, while the last one executes it (strictly speaking,
make run will also compile the code if the executable doesn't exist yet, so you could have skipped the second command if you wanted). This is all that's needed to run the code and produce the output that is discussed in the "Results" section of the tutorial programs. This sequence needs to be repeated in all of the tutorial directories you want to play with.
When learning the library, you need to play with it and see what happens. To this end, open the
examples/step-1/step-1.cc source file with your favorite editor and modify it in some way, save it and run it as above. A few suggestions for possibly modifications are given at the end of the results section of this program, where we also provide a few links to other useful pieces of information.
This and several of the other tutorial programs are also discussed and demonstrated in Wolfgang Bangerth's video lectures on deal.II and computational science. In particular, you can see the steps he executes to run this and other programs, and you will get a much better idea of the tools that can be used to work with deal.II. In particular, lectures 2 and 4 give an overview of deal.II and of the building blocks of any finite element code. (See also video lecture 2, video lecture 4.)
If you are not yet familiar with using Linux and running things on the command line, you may be interested in watching lectures 2.9 and 2.91. (See also video lecture 2.9, video lecture 2.91.) These give overviews over the command line and on what happens when compiling programs, respectively.
Note that deal.II is actively developed, and in the course of this development we occasionally rename or deprecate functions or classes that are still referenced in these video lectures. For example, the step-1 code shown in video lecture 5 uses a class HyperShellBoundary which was replaced with SphericalManifold class later on. However, in spite of such kind of details the content of the lectures is still very relevant.
Let's come back to step-1, the current program. In this first example, we don't actually do very much, but show two techniques: what is the syntax to generate triangulation objects, and some elements of simple loops over all cells. We create two grids, one which is a regularly refined square (not very exciting, but a common starting grid for some problems), and one more geometric attempt: a ring-shaped domain, which is refined towards the inner edge. Through this, you will get to know three things every finite element program will have to have somewhere: An object of type Triangulation for the mesh; a call to the GridGenerator functions to generate a mesh; and loops over all cells that involve iterators (iterators are a generalization of pointers and are frequently used in the C++ standard library; in the context of deal.II, the Iterators on mesh-like containers module talks about them).
The program is otherwise small enough that it doesn't need a whole lot of introduction.
If you are reading through this tutorial program, chances are that you are interested in continuing to use deal.II for your own projects. Thus, you are about to embark on an exercise in programming using a large-scale scientific computing library. Unless you are already an experienced user of large-scale programming methods, this may be new territory for you — with all the new rules that go along with it such as the fact that you will have to deal with code written by others, that you may have to think about documenting your own code because you may not remember what exactly it is doing a year down the road (or because others will be using it as well), or coming up with ways to test that your program is doing the right thing. None of this is something that we typically train mathematicians, engineers, or scientists in but that is important when you start writing software of more than a few hundred lines. Remember: Producing software is not the same as just writing code.
To make your life easier on this journey let us point to three resources that are worthwhile browsing through before you start any large-scale programming:
As a general recommendation: If you expect to spend more than a few days writing software in the future, do yourself the favor of learning tools that can make your life more productive, in particular debuggers and integrated development environments. (See also video lecture 7, video lecture 8, video lecture 8.01, video lecture 25.) You will find that you will get the time spent learning these tools back severalfold soon by being more productive! Several of the video lectures referenced above show how to use tools such as integrated development environments or debuggers.
The most fundamental class in the library is the Triangulation class, which is declared here:
We need the following two includes for loops over cells and/or faces:
Here are some functions to generate standard grids:
We would like to use faces and cells which are not straight lines, or bi-linear quads, so we import some classes which predefine some manifold descriptions:
Output of grids in various graphics formats:
This is needed for C++ output:
And this for the declarations of the `sqrt' and `fabs' functions:
The final step in importing deal.II is this: All deal.II functions and classes are in a namespace
dealii, to make sure they don't clash with symbols from other libraries you may want to use in conjunction with deal.II. One could use these functions and classes by prefixing every use of these names by
::, but that would quickly become cumbersome and annoying. Rather, we simply import the entire deal.II namespace for general use:
In the following, first function, we simply use the unit square as domain and produce a globally refined grid from it.
The first thing to do is to define an object for a triangulation of a two-dimensional domain:
Here and in many following cases, the string "<2>" after a class name indicates that this is an object that shall work in two space dimensions. Likewise, there are versions of the triangulation class that are working in one ("<1>") and three ("<3>") space dimensions. The way this works is through some template magic that we will investigate in some more detail in later example programs; there, we will also see how to write programs in an essentially dimension independent way.
Next, we want to fill the triangulation with a single cell for a square domain. The triangulation is the refined four times, to yield \(4^4=256\) cells in total:
Now we want to write a graphical representation of the mesh to an output file. The GridOut class of deal.II can do that in a number of different output formats; here, we choose encapsulated postscript (eps) format:
The grid in the following, second function is slightly more complicated in that we use a ring domain and refine the result once globally.
We start again by defining an object for a triangulation of a two-dimensional domain:
We then fill it with a ring domain. The center of the ring shall be the point (1,0), and inner and outer radius shall be 0.5 and 1. The number of circumferential cells could be adjusted automatically by this function, but we choose to set it explicitly to 10 as the last argument:
By default, the triangulation assumes that all boundaries are straight lines, and all cells are bi-linear quads or tri-linear hexes, and that they are defined by the cells of the coarse grid (which we just created). Unless we do something special, when new points need to be introduced; the domain is assumed to be delineated by the straight lines of the coarse mesh, and new points will simply be in the middle of the surrounding ones. Here, however, we know that the domain is curved, and we would like to have the Triangulation place new points according to the underlying geometry. Fortunately, some good soul implemented an object which describes a spherical domain, of which the ring is a section; it only needs the center of the ring and automatically figures out how to instruct the Triangulation where to place the new points. The way this works in deal.II is that you tag parts of the triangulation you want to be curved with a number that is usually referred to as "manifold indicator" and then tell the triangulation to use a particular "manifold object" for all places with this manifold indicator. How exactly this works is not important at this point (you can read up on it in step-53 and Manifold description for triangulations). Here, for simplicity, we will choose the manifold id to be zero. By default, all cells and faces of the Triangulation have their manifold_id set to numbers::invalid_manifold_id, which is the default if you want a manifold that produces straight edges, but you can change this number for individual cells and faces. In that case, the curved manifold thus associated with number zero will not apply to those parts with a non-zero manifold indicator, but other manifold description objects can be associated with those non-zero indicators. If no manifold description is associated with a particular manifold indicator, a manifold that produces straight edges is implied. (Manifold indicators are a slightly complicated topic; if you're confused about what exactly is happening here, you may want to look at the glossary entry on this topic".)
In order to demonstrate how to write a loop over all cells, we will refine the grid in five steps towards the inner circle of the domain:
Next, we need an iterator that points to a cell and which we will move over all active cells one by one. In a sense, you can think of a triangulation as a collection of cells. If it was an array, you would just get a pointer that you move from one array element to the next. In triangulations, cells aren't stored as an array, so simple pointers do not work, but one can generalize pointers to iterators (see this wikipedia link for more information). We will then get an iterator to the first cell and iterate over all of the cells until we hit the last one.
The second important piece is that we only need the active cells. Active cells are those that are not further refined, and the only ones that can be marked for further refinement, obviously. deal.II provides iterator categories that allow us to iterate over all cells (including the parent cells of active ones) or only over the active cells. Because we want the latter, we need to choose Triangulation::active_cell_iterator as data type.
Finally, by convention, we almost always use the names
endc for the iterator pointing to the present cell and to the "one-past-the-end" iterator. This is, in a sense a misnomer, because the object is not really a "cell": it is an iterator/pointer to a cell. We should really have started to call these objects
cell_iterator when deal.II started in 1998, but it is what it is.
After declaring the iterator variable, the loop over all cells is then rather trivial, and looks like any loop involving pointers instead of iterators:
Next, we want to loop over all vertices of the cells. Since we are in 2d, we know that each cell has exactly four vertices. However, instead of penning down a 4 in the loop bound, we make a first attempt at writing it in a dimension-independent way by which we find out about the number of vertices of a cell. Using the GeometryInfo class, we will later have an easier time getting the program to also run in 3d: we only have to change all occurrences of
<3>, and do not have to audit our code for the hidden appearance of magic numbers like a 4 that needs to be replaced by an 8:
If this cell is at the inner boundary, then at least one of its vertices must sit on the inner ring and therefore have a radial distance from the center of exactly 0.5, up to floating point accuracy. Compute this distance, and if we have found a vertex with this property flag this cell for later refinement. We can then also break the loop over all vertices and move on to the next cell.
Now that we have marked all the cells that we want refined, we let the triangulation actually do this refinement. The function that does so owes its long name to the fact that one can also mark cells for coarsening, and the function does coarsening and refinement all at once:
Finally, after these five iterations of refinement, we want to again write the resulting mesh to a file, again in eps format. This works just as above:
At this point, all objects created in this function will be destroyed in reverse order. Unfortunately, we defined the manifold object after the triangulation, which still has a pointer to it and the library will produce an error if the manifold object is destroyed before the triangulation. We therefore have to release it, which can be done as follows. Note that this sets the manifold object used for part "0" of the domain back to a default object, over which the triangulation has full control.
An alternative to doing so, and one that is frequently more convenient, would have been to declare the manifold object before the triangulation object. In that case, the triangulation would have let lose of the manifold object upon its destruction, and everything would have been fine.
Finally, the main function. There isn't much to do here, only to call the two subfunctions, which produce the two grids.
Running the program produces graphics of two grids (grid-1.eps and grid-2.eps). They look like this:
The left one, well, is not very exciting. The right one is — at least — unconventional.
While the second mesh is entirely artificial and made-up, and certainly not very practical in applications, to everyone's surprise it has found its way into the literature: see the paper by M. Mu titled "PDE.MART: A network-based problem-solving environment", ACM Trans. Math. Software, vol. 31, pp. 508-531, 2005. Apparently it is good for some things at least.
This program obviously does not have a whole lot of functionality, but in particular the
second_grid function has a bunch of places where you can play with it. For example, you could modify the criterion by which we decide which cells to refine. An example would be to change the condition to this:
This would refine all cells for which the \(y\)-coordinate of the cell's center is greater than zero (the
TriaAccessor::center function that we call by dereferencing the
cell iterator returns a Point<2> object; subscripting
 would give the \(x\)-coordinate, subscripting
 the \(y\)-coordinate). By looking at the functions that TriaAccessor provides, you can also use more complicated criteria for refinement.
Another possibility would be to generate meshes of entirely different geometries altogether. While for complex geometries there is no way around using meshes obtained from mesh generators, there is a good number of geometries for which deal.II can create meshes using the functions in the GridGenerator namespace. Take a look at what it provides and see how it could be used in a program like this.
We also discuss a variety of other ways to create and manipulate meshes to step-49.
We close with a comment about modifying or writing programs with deal.II in general. When you start working with tutorial programs or your own applications, you will find that mistakes happen: your program will contain code that either aborts the program right away or bugs that simply lead to wrong results. In either case, you will find it extremely helpful to know how to work with a debugger: you may get by for a while by just putting debug output into your program, compiling it, and running it, but ultimately finding bugs with a debugger is much faster, much more convenient, and more reliable because you don't have to recompile the program all the time and because you can inspect the values of variables and how they change.
Rather than postponing learning how to use a debugger till you really can't see any other way to find a bug, here's the one piece of advice we will provide in this program: learn how to use a debugger as soon as possible. It will be time well invested. The deal.II Frequently Asked Questions (FAQ) page linked to from the top-level deal.II webpage also provides a good number of hints on debugging deal.II programs.