]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Rewrite part of the text to make it smoother and lass sloppy.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 24 Jan 2006 20:22:56 +0000 (20:22 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Tue, 24 Jan 2006 20:22:56 +0000 (20:22 +0000)
git-svn-id: https://svn.dealii.org/trunk@12146 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-1/step-1.cc

index a00af50774e66eeec48fcb3333f821abee8cdb92..ec93ec69436963e93c1bb56be863e397a3f87b62 100644 (file)
@@ -9,6 +9,8 @@
  * 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);
@@ -86,26 +92,30 @@ void first_grid ()
 
 
 
-                                 // 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;
@@ -116,8 +126,17 @@ void second_grid ()
                                    // 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
@@ -127,82 +146,140 @@ void second_grid ()
                                    // 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);
@@ -226,13 +303,26 @@ void second_grid ()
                                    // 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 ();

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.