Table of contents | |
---|---|
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. The latter is certainly not very useful and is probably only rarely used in numerical analysis for PDEs (although, to everyone's surprise, it has actually 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 :-), but looks nice and illustrates how loops over cells are written and some of the things you can do with cells.
The most fundamental class in the library is the Triangulation class, which is declared here:
#include <grid/tria.h>
We need the following two includes for loops over cells and/or faces:
#include <grid/tria_accessor.h> #include <grid/tria_iterator.h>
Here are some functions to generate standard grids:
#include <grid/grid_generator.h>
We would like to use boundaries which are not straight lines, so we import some classes which predefine some boundary descriptions:
#include <grid/tria_boundary_lib.h>
Output of grids in various graphics formats:
#include <grid/grid_out.h>
This is needed for C++ output:
#include <fstream>
And this for the declarations of the `sqrt' and `fabs' functions:
#include <cmath>
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:
using namespace dealii;
In the following, first function, we simply use the unit square as domain and produce a globally refined grid from it.
void first_grid ()
{
The first thing to do is to define an object for a triangulation of a two-dimensional domain:
Triangulation<2> triangulation;
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:
GridGenerator::hyper_cube (triangulation); triangulation.refine_global (4);
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:
std::ofstream out ("grid-1.eps");
GridOut grid_out;
grid_out.write_eps (triangulation, out);
}
The grid in the following, second function is slightly more complicated in that we use a ring domain and refine the result once globally.
void second_grid ()
{
We start again by defining an object for a triangulation of a two-dimensional domain:
Triangulation<2> triangulation;
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 explicitely to 10 as the last argument:
const Point<2> center (1,0); const double inner_radius = 0.5, outer_radius = 1.0; GridGenerator::hyper_shell (triangulation, center, inner_radius, outer_radius, 10);
By default, the triangulation assumes that all boundaries are straight and given by the cells of the coarse grid (which we just created). It uses this information when cells at the boundary are refined and new points need to be introduced on the boundary; if the boundary is assumed to be straight, then new points will simply be in the middle of the surrounding ones.
Here, however, we would like to have a curved boundary. Fortunately, some good soul implemented an object which describes the boundary of a ring domain; it only needs the center of the ring and automatically figures out the inner and outer radius when needed. Note that we associate this boundary object with that part of the boundary that has the "boundary indicator" zero. By default, all boundary parts have this number, but you can change this number for some parts of the boundary. In that case, the curved boundary thus associated with number zero will not apply on those parts with a non-zero boundary indicator, but other boundary description objects can be associated with those non-zero indicators. If no boundary description is associated with a particular boundary indicator, a straight boundary is implied.
const HyperShellBoundary<2> boundary_description(center);
triangulation.set_boundary (0, boundary_description);
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:
for (unsigned int step=0; step<5; ++step) {
Next, we need an iterator which points to a cell and which we will move over all active cells one by one (active cells are those that are not further refined, and the only ones that can be marked for further refinement, obviously). By convention, we almost always use the names cell
and endc
for the iterator pointing to the present cell and to the one-past-the-end
iterator:
Triangulation<2>::active_cell_iterator cell = triangulation.begin_active(), endc = triangulation.end();
The loop over all cells is then rather trivial, and looks like any loop involving pointers instead of iterators:
for (; cell!=endc; ++cell)
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 <2>
to <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:
for (unsigned int v=0; v < GeometryInfo<2>::vertices_per_cell; ++v) {
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.
const double distance_from_center = center.distance (cell->vertex(v)); if (std::fabs(distance_from_center - inner_radius) < 1e-10) { cell->set_refine_flag (); break; } }
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:
triangulation.execute_coarsening_and_refinement (); }
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:
std::ofstream out ("grid-2.eps");
GridOut grid_out;
grid_out.write_eps (triangulation, out);
At this point, all objects created in this function will be destroyed in reverse order. Unfortunately, we defined the boundary object after the triangulation, which still has a pointer to it and the library will produce an error if the boundary object is destroyed before the triangulation. We therefore have to release it, which can be done as follows. Note that this sets the boundary object used for part "0" of the boundary back to a default object, over which the triangulation has full control.
triangulation.set_boundary (0);
An alternative to doing so, and one that is frequently more convenient, would have been to declare the boundary object before the triangulation object. In that case, the triangulation would have let lose of the boundary 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.
int main ()
{
first_grid ();
second_grid ();
}
The program has, after having been run, produced two grids, which look like this:
![]() | ![]() |
The left one, well, is not very exciting. The right one is — at least — unconventional.
(If you are looking at a locally installed deal.II version, then the program can be found at /home /hazelsct /repositories /deal.ii /examples /step-1 /step-1.cc . Otherwise, this is only the path on some remote server.)
/ * Subversion Id: @ref step_1 "step-1".cc 19994 2009-10-25 02:47:44Z bangerth * * Copyright (C) 1999, 2000, 2001, 2002, 2003, 2005, 2006, 2007, 2009 by the deal.II authors * * This file is subject to QPL and may not be distributed * without copyright and license information. Please refer * to the file deal.II/doc/license.html for the text and * further information on this license. * /
#include <grid/tria.h> #include <grid/tria_accessor.h> #include <grid/tria_iterator.h> #include <grid/grid_generator.h> #include <grid/tria_boundary_lib.h> #include <grid/grid_out.h> #include <fstream> #include <cmath> using namespace dealii;
void first_grid () { Triangulation<2> triangulation; GridGenerator::hyper_cube (triangulation); triangulation.refine_global (4); std::ofstream out ("grid-1.eps"); GridOut grid_out; grid_out.write_eps (triangulation, out); }
void second_grid () { Triangulation<2> triangulation; const Point<2> center (1,0); const double inner_radius = 0.5, outer_radius = 1.0; GridGenerator::hyper_shell (triangulation, center, inner_radius, outer_radius, 10); const HyperShellBoundary<2> boundary_description(center); triangulation.set_boundary (0, boundary_description); for (unsigned int step=0; step<5; ++step) { Triangulation<2>::active_cell_iterator cell = triangulation.begin_active(), endc = triangulation.end(); for (; cell!=endc; ++cell) for (unsigned int v=0; v < GeometryInfo<2>::vertices_per_cell; ++v) { const double distance_from_center = center.distance (cell->vertex(v)); if (std::fabs(distance_from_center - inner_radius) < 1e-10) { cell->set_refine_flag (); break; } } triangulation.execute_coarsening_and_refinement (); } std::ofstream out ("grid-2.eps"); GridOut grid_out; grid_out.write_eps (triangulation, out); triangulation.set_boundary (0); }
int main ()
{
first_grid ();
second_grid ();
}