* further information on this license.
*/
+ // @sect3{Include files}
+
// The most fundamental class in the
// library is the ``Triangulation''
// class, which is declared here:
// `sqrt' and `fabs' functions:
#include <cmath>
+ // @sect3{Creating the first mesh}
- // In the following function, we
+ // In the following, first function, we
// simply use the unit square as
// domain and produce a globally
// refined grid from it.
void first_grid ()
{
- // Define an object for a
+ // 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
- // version working in one ("<1>")
- // and three ("<3>") space
- // dimensions, or for all
- // dimensions. We will see such
- // constructs in later examples,
- // where we show how to program
- // dimension independently.
- // (At present, only one through
- // three space dimensions are
- // supported, but that is not a
- // restriction. In case someone
- // would like to implement four
- // dimensional finite elements, for
- // example for general relativity,
- // this would be a straightforward
- // thing.)
+ // 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.
- // Fill it with a square
+ // 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);
-
- // Refine all cells four times, to
- // yield 4^4=256 cells in total
triangulation.refine_global (4);
- // Now we want to write it to some
- // output, here in postscript
- // format
+ // 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 function
- // is slightly more complicated in
- // that we use a ring domain and
- // refine the result once globally
+ // @sect3{Creating the second mesh}
+
+ // 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 ()
{
- // Define an object for a
- // triangulation of a
- // two-dimensional domain
+ // We start again by defining an
+ // object for a triangulation of a
+ // two-dimensional domain:
Triangulation<2> triangulation;
- // 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 as the
- // last argument
+ // 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;
// assumes that all boundaries are
// straight and given by the cells
// of the coarse grid (which we
- // just created). Here, however, we
- // would like to have a curved
+ // 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
// radius when needed. Note that we
// associate this boundary object
// with that part of the boundary
- // that has the "boundary number"
- // zero. By default, all boundary
- // parts have this number, but you
- // might want to change this number
- // for some parts, and then the
+ // 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
- // there.
+ // 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);
- // Now, just for the purpose of
- // demonstration and for no
- // particular reason, we will
- // refine the grid in five steps
- // towards the inner circle of the
- // domain:
+ // 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)
{
- // Get 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
- Triangulation<2>::active_cell_iterator cell, endc;
- cell = triangulation.begin_active();
- endc = triangulation.end();
-
- // Now loop over all cells...
+ // 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)
- // ...and over all vertices
- // of the cells. Note the
+ // 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 faces of a
- // cell
- for (unsigned int ivertex=0;
- ivertex < GeometryInfo<2>::vertices_per_cell;
- ++ivertex)
+ // 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 vertex=0;
+ vertex < GeometryInfo<2>::vertices_per_cell;
+ ++vertex)
{
// If this cell is at the
// inner boundary, then
- // at least one of its vertices
- // must have a radial
- // distance from the center
- // of 0.5
- const Point<2> vector_to_center
- = (cell->vertex(ivertex) - center);
+ // 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
- = std::sqrt(vector_to_center.square());
+ = center.distance (cell->vertex(vertex));
if (std::fabs(distance_from_center - inner_radius) < 1e-10)
{
- // Ok, this is one of
- // the cells we were
- // looking for. Flag
- // it for refinement
- // and go to the next
- // cell by breaking
- // the loop over all
- // vertices
cell->set_refine_flag ();
break;
};
};
- // Refine the cells which we
- // have marked
+ // 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 ();
};
- // Now we want to write it to some
- // output, here in postscript
- // format
+ // 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);
// 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.
}
- // Main function. Only call the two
- // subfunctions, which produce the
- // two grids.
+ // @sect3{Creating the second mesh}
+
+ // 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 ();