]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Finish step-6.
authorwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 6 Jan 2000 15:52:28 +0000 (15:52 +0000)
committerwolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 6 Jan 2000 15:52:28 +0000 (15:52 +0000)
git-svn-id: https://svn.dealii.org/trunk@2161 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/Attic/examples/step-by-step/step-6/Makefile
deal.II/deal.II/Attic/examples/step-by-step/step-6/step-6.cc
deal.II/examples/step-6/Makefile
deal.II/examples/step-6/step-6.cc

index dc2798f83951ccc2f42089e835bfb57050d1cd1a..8753051799cf996c1cb4b25beca741c301fc978d 100644 (file)
@@ -50,7 +50,7 @@ $(target) : $(libraries)
 
 # rule how to run the program
 run: $(target)
-       @echo ============================ Running $@
+       @echo ============================ Running $<
        @./$(target)
 
 
index 04af2fa3d95b34e70df8047d250474e5daeea7dc..2c89c63110fbf6cb13551bcf66423df63272591d 100644 (file)
@@ -1,8 +1,10 @@
 /* $Id$ */
 /* Author: Wolfgang Bangerth, University of Heidelberg, 2000 */
 
-
-                                // still unfinished
+                                // The first few files have already
+                                // been covered in previous examples
+                                // and will thus not be further
+                                // commented on.
 #include <base/quadrature_lib.h>
 #include <base/function.h>
 #include <base/logstream.h>
@@ -17,6 +19,7 @@
 #include <grid/grid_generator.h>
 #include <grid/tria_accessor.h>
 #include <grid/tria_iterator.h>
+#include <grid/tria_boundary_lib.h>
 #include <dofs/dof_accessor.h>
 #include <dofs/dof_tools.h>
 #include <fe/fe_lib.lagrange.h>
 #include <numerics/matrices.h>
 #include <numerics/data_out.h>
 
-                                // out statt in
+                                // We will not read the grid from a
+                                // file as in the previous example,
+                                // but generate it using a function
+                                // of the library. However, we will
+                                // want to write out the locally
+                                // refined grids in each step, so we
+                                // need the following include file
+                                // instead of ``grid_in.h'':
 #include <grid/grid_out.h>
 
-#include <grid/tria_boundary_lib.h>
-
-                                //...
+                                // When using locally refined grids,
+                                // we will get so-called ``hanging
+                                // nodes''. However, the standard
+                                // finite element methods assumes
+                                // that the discrete solution spaces
+                                // be continuous, so we need to make
+                                // sure that the degrees of freedom
+                                // on hanging nodes conform to some
+                                // constraints such that the global
+                                // solution is continuous. The
+                                // following file contains a class
+                                // which is used to handle these
+                                // constraints:
 #include <dofs/dof_constraints.h>
+
+                                // Finally, we would like to use a
+                                // simple way to adaptively refine
+                                // the grid. While in general,
+                                // adaptivity is very
+                                // problem-specific, the error
+                                // indicator in the following file
+                                // often yields quite nicely adapted
+                                // grids for a wide class of
+                                // problems.
 #include <numerics/error_estimator.h>
 
 #include <fstream>
 
 
-
+                                // The main class is again almost
+                                // unchanged. Two additions, however,
+                                // are made: we have added the
+                                // ``refine'' function, which is used
+                                // to adaptively refine the grid
+                                // (instead of the global refinement
+                                // in the previous examples), and a
+                                // variable which will hold the
+                                // constraints associated to the
+                                // hanging nodes.
 template <int dim>
 class LaplaceProblem 
 {
@@ -56,7 +95,7 @@ class LaplaceProblem
     FEQ1<dim>            fe;
     DoFHandler<dim>      dof_handler;
 
-                                    // This is the only addition to
+                                    // This is the new variable in
                                     // the main class. We need an
                                     // object which holds a list of
                                     // the constraints from the
@@ -343,22 +382,164 @@ void LaplaceProblem<dim>::solve ()
 
                                 // Instead of global refinement, we
                                 // now use a slightly more elaborate
-                                // scheme.
-                                // ...
+                                // scheme. We will use the
+                                // ``KellyErrorEstimator'' class
+                                // which implements an error
+                                // estimator for the Laplace
+                                // equation; it can in principle
+                                // handle variable coefficients, but
+                                // we will not use these advanced
+                                // features, but rather use its most
+                                // simple form since we are not
+                                // interested in quantitative results
+                                // but only in a quick way to
+                                // generate locally refined grids.
+                                //
+                                // Although the error estimator
+                                // derived by Kelly et al. was
+                                // originally developed for Laplace's
+                                // equation, we have found that it is
+                                // also well suited to quickly
+                                // generate locally refined grids for
+                                // a wide class of
+                                // problems. Basically, it looks at
+                                // the jumps of the gradients of the
+                                // solution over the faces of cells
+                                // (which is a measure for the second
+                                // derivatives) and scales it by the
+                                // size of the cell. It is therefore
+                                // a measure for the local smoothness
+                                // of the solution at the place of
+                                // each cell and it is thus
+                                // understandable that it yields
+                                // reasonable grids also for
+                                // hyperbolic transport problems or
+                                // the wave equation as well,
+                                // although these grids are certainly
+                                // suboptimal compared to approaches
+                                // specially tailored to the
+                                // problem. This error estimator may
+                                // therefore be understood as a quick
+                                // way to test an adaptive program.
 template <int dim>
 void LaplaceProblem<dim>::refine_grid ()
 {
+                                  // The output of the error
+                                  // estimator class is an error
+                                  // indicator for each cell. We
+                                  // therefore need a vector with as
+                                  // many elements as there are
+                                  // active cells. Since accuracy is
+                                  // not that important here, the
+                                  // data type for the error values
+                                  // on each cell is ``float''
+                                  // instead of ``double''.
   Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
 
+                                  // Next, the error estimator can
+                                  // handle Neumann boundary
+                                  // conditions. For this, it needs
+                                  // to know which parts of the
+                                  // boundary have Neumann boundary
+                                  // conditions and teh respective
+                                  // boundary values there. This
+                                  // information is mediated by a map
+                                  // in which the keys are the
+                                  // boundary part numbers and the
+                                  // values are pointers to the
+                                  // boundary value functions. We
+                                  // create such a map, but since we
+                                  // do not use Neumann boundary
+                                  // conditions, the map will not
+                                  // contain entries.
   KellyErrorEstimator<dim>::FunctionMap neumann_boundary;
+
+                                  // Now we call the error
+                                  // estimator. The parameters should
+                                  // be clear apart from the
+                                  // quadrature formula: as said
+                                  // above, the jump of the gradients
+                                  // of the solution across the faces
+                                  // of a cell are considered. They
+                                  // are integrated along the face,
+                                  // but as usual in finite element
+                                  // programs the integration is done
+                                  // using quadrature. Since the
+                                  // error estimator class can't know
+                                  // itself which quadrature formula
+                                  // might be appropriate, we have to
+                                  // pass one to the function. Note
+                                  // that since the quadrature has to
+                                  // take place along faces, the
+                                  // dimension of the quadrature
+                                  // formula is ``dim-1'' rather then
+                                  // ``dim''.
   KellyErrorEstimator<dim>::estimate (dof_handler,
                                      QGauss3<dim-1>(),
                                      neumann_boundary,
                                      solution,
                                      estimated_error_per_cell);
 
+                                  // The above function returned one
+                                  // error indicator value for each
+                                  // cell in the
+                                  // ``estimated_error_per_cell''
+                                  // array. Refinement is now done as
+                                  // follows: refine those 30 per
+                                  // cent of the cells with the
+                                  // highest error values, and
+                                  // coarsen the 3 per cent of cells
+                                  // with the lowest values.
+                                  //
+                                  // One can easily verify that if
+                                  // the second number were zero,
+                                  // this would approximately result
+                                  // in a doubling of cells in each
+                                  // step in two space dimensions,
+                                  // since for each of the 30 per
+                                  // cent of cells four new would be
+                                  // replaced. In practice, some more
+                                  // cells are usually produced since
+                                  // it is disallowed that a cell is
+                                  // refined twice while the neighbor
+                                  // cell is not refined; in that
+                                  // case, the neighbor cell would be
+                                  // refined as well.
+                                  //
+                                  // In many applications, the number
+                                  // of cells to be coarsened would
+                                  // be set to something larger than
+                                  // only three per cent. A non-zero
+                                  // value is useful especially if
+                                  // for some reason the initial
+                                  // (coarse) grid is already rather
+                                  // refined. In that case, it might
+                                  // be necessary to refine it in
+                                  // some regions, while coarsening
+                                  // in some other regions is
+                                  // useful. In our case here, the
+                                  // initial grid is very coarse, so
+                                  // coarsening is only necessary in
+                                  // a few regions where
+                                  // over-refinement may have taken
+                                  // place. Thus a small, non-zero
+                                  // value is appropriate here.
   triangulation.refine_and_coarsen_fixed_number (estimated_error_per_cell,
                                                 0.3, 0.03);
+
+                                  // After the previous function has
+                                  // exited, some cells are flagged
+                                  // for refinement, and some other
+                                  // for coarsening. The refinement
+                                  // or coarsening itself is not
+                                  // performed by now, however, since
+                                  // there are many cases where
+                                  // further modifications of these
+                                  // flags is useful. Here, we don't
+                                  // want to do any such thing, so we
+                                  // can tell the triangulation to
+                                  // perform the actions for which
+                                  // the cells are flagged.
   triangulation.execute_coarsening_and_refinement ();
 };
 
@@ -367,7 +548,18 @@ void LaplaceProblem<dim>::refine_grid ()
 template <int dim>
 void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
 {
-                                  // ...
+                                  // We want to write the grid in
+                                  // each cycle. Here is another way
+                                  // to quickly produce a filename
+                                  // based on the cycle number. It
+                                  // assumes that the numbers `0'
+                                  // through `9' are represented
+                                  // consecutively in the character
+                                  // set (which is the case in all
+                                  // known character sets). However,
+                                  // this will only work if the cycle
+                                  // number is less than ten, which
+                                  // we check by an assertion.
   string filename = "grid-";
   filename += ('0' + cycle);
   Assert (cycle < 10, ExcInternalError());
@@ -375,12 +567,13 @@ void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
   filename += ".eps";
   ofstream output (filename.c_str());
 
+                                  // Using this filename, we write
+                                  // each grid as a postscript file.
   GridOut grid_out;
   grid_out.write_eps (triangulation, output);
 };
 
 
-
 template <int dim>
 void LaplaceProblem<dim>::run () 
 {
@@ -390,7 +583,32 @@ void LaplaceProblem<dim>::run ()
 
       if (cycle == 0)
        {
-                                          //...
+                                          // Instead of reading the
+                                          // grid from a file on disk
+                                          // as in the previous
+                                          // example, we now again
+                                          // create it using a
+                                          // library function. The
+                                          // domain is again a
+                                          // circle, which is why we
+                                          // have to provide a
+                                          // suitable boundary object
+                                          // as well.
+                                          //
+                                          // You will notice by
+                                          // looking at the coarse
+                                          // grid that it is of
+                                          // inferior quality than
+                                          // the one which we read
+                                          // from the file in the
+                                          // previous example: the
+                                          // cells are less equally
+                                          // formed. However, using
+                                          // the library function
+                                          // this program works in
+                                          // any space dimension,
+                                          // which was not the case
+                                          // before.
          GridGenerator::hyper_ball (triangulation);
 
          static const HyperBallBoundary<dim> boundary;
@@ -399,8 +617,19 @@ void LaplaceProblem<dim>::run ()
          triangulation.refine_global (1);
        }
       else
-                                        // ...
-       refine_grid ();
+                                        // In case this is not the
+                                        // first cycle, we want to
+                                        // refine the grid. Unlike
+                                        // the global refinement
+                                        // employed in the last
+                                        // example, we now use the
+                                        // adaptive procedure
+                                        // described in the function
+                                        // which we now call:
+       {
+         refine_grid ();
+       };
+      
 
       cout << "   Number of active cells: "
           << triangulation.n_active_cells()
@@ -412,7 +641,14 @@ void LaplaceProblem<dim>::run ()
       output_results (cycle);
     };
 
-                                  // ...
+                                  // The solution on the final grid
+                                  // is now written to a file. As
+                                  // already done in one of the
+                                  // previous examples, we use the
+                                  // EPS format for output, and to
+                                  // obtain a reasonable view on the
+                                  // solution, we rescale the z-axis
+                                  // by a factor of four.
   DataOut<dim>::EpsFlags eps_flags;
   eps_flags.z_scaling = 4;
   
@@ -428,13 +664,128 @@ void LaplaceProblem<dim>::run ()
 };
 
     
-
+                                // The main function is unaltered in
+                                // its functionality against the
+                                // previous example, but we have
+                                // taken a step of additional
+                                // caution. Sometimes, something goes
+                                // wrong (such as insufficient disk
+                                // space upon writing an output file,
+                                // not enough memory when trying to
+                                // allocate a vector or a matrix, or
+                                // if we can't read from or write to
+                                // a file for whatever reason), and
+                                // in these cases the library will
+                                // throw exceptions. Since they do
+                                // not constitute programming errors,
+                                // these exceptions also are not
+                                // switched off in optimized mode, in
+                                // constrast to the ``Assert'' macro
+                                // which we have used to test against
+                                // programming errors. If uncought,
+                                // these exceptions propagate the
+                                // call tree up to the ``main''
+                                // function, and if they are not
+                                // caught there either, the program
+                                // is aborted. In many cases, like if
+                                // there is not enough memory or disk
+                                // space, we can't do anything but we
+                                // can at least print some text
+                                // trying to explain the reason why
+                                // the program failed. A way to do so
+                                // is shown in the following. It is
+                                // certainly useful to write any
+                                // larger program in this way, and
+                                // you can do so by more or less
+                                // copying this function apart from
+                                // the ``try'' block which contains
+                                // the code that constitutes the
+                                // actual functionality.
 int main () 
 {
-  deallog.depth_console (0);
 
-  LaplaceProblem<2> laplace_problem_2d;
-  laplace_problem_2d.run ();
-  
+                                  // The general idea behind the
+                                  // layout of this function is as
+                                  // follows: let's try to run the
+                                  // program as we did before...
+  try
+    {
+      deallog.depth_console (0);
+
+      LaplaceProblem<2> laplace_problem_2d;
+      laplace_problem_2d.run ();
+    }
+                                  // ...and if this should fail, try
+                                  // to gather as much information as
+                                  // possible. Specifically, if the
+                                  // exception that was thrown is an
+                                  // object of a class that is
+                                  // derived from the C++ standard
+                                  // class ``exception'', then we can
+                                  // use the ``what'' member function
+                                  // to get a string which describes
+                                  // the reason why the exception was
+                                  // thrown. 
+                                  //
+                                  // The deal.II exception classes
+                                  // are all derived from the
+                                  // standard class, and in
+                                  // particular, the ``exc.what()''
+                                  // function will return
+                                  // approximately the same string as
+                                  // would be generated if the
+                                  // exception was thrown using the
+                                  // ``Assert'' macro. You have seen
+                                  // the output of such an exception
+                                  // in the previous example, and you
+                                  // then know that it contains the
+                                  // file and line number of where
+                                  // the exception occured, and some
+                                  // other information. This is also
+                                  // what would be printed in the
+                                  // following.
+  catch (exception &exc)
+    {
+      cerr << endl << endl
+          << "----------------------------------------------------"
+          << endl;
+      cerr << "Exception on processing: " << exc.what() << endl
+          << "Aborting!" << endl
+          << "----------------------------------------------------"
+          << endl;
+                                      // We can't do much more than
+                                      // printing as much information
+                                      // as we can get to, so abort
+                                      // with error:
+      return 1;
+    }
+                                  // If the exception that was thrown
+                                  // somewhere was not an object of a
+                                  // class derived from the standard
+                                  // ``exception'' class, then we
+                                  // can't do anything at all. We
+                                  // then simply print an error
+                                  // message and exit.
+  catch (...) 
+    {
+      cerr << endl << endl
+          << "----------------------------------------------------"
+          << endl;
+      cerr << "Unknown exception!" << endl
+          << "Aborting!" << endl
+          << "----------------------------------------------------"
+          << endl;
+      return 1;
+    };
+
+                                  // If we got to this point, there
+                                  // was no exception which
+                                  // propagated up to the main
+                                  // functino (maybe there were some,
+                                  // but they were caught somewhere
+                                  // in the program or the
+                                  // library). Therefore, the program
+                                  // performed as was expected and we
+                                  // can return without error.
   return 0;
 };
index dc2798f83951ccc2f42089e835bfb57050d1cd1a..8753051799cf996c1cb4b25beca741c301fc978d 100644 (file)
@@ -50,7 +50,7 @@ $(target) : $(libraries)
 
 # rule how to run the program
 run: $(target)
-       @echo ============================ Running $@
+       @echo ============================ Running $<
        @./$(target)
 
 
index 04af2fa3d95b34e70df8047d250474e5daeea7dc..2c89c63110fbf6cb13551bcf66423df63272591d 100644 (file)
@@ -1,8 +1,10 @@
 /* $Id$ */
 /* Author: Wolfgang Bangerth, University of Heidelberg, 2000 */
 
-
-                                // still unfinished
+                                // The first few files have already
+                                // been covered in previous examples
+                                // and will thus not be further
+                                // commented on.
 #include <base/quadrature_lib.h>
 #include <base/function.h>
 #include <base/logstream.h>
@@ -17,6 +19,7 @@
 #include <grid/grid_generator.h>
 #include <grid/tria_accessor.h>
 #include <grid/tria_iterator.h>
+#include <grid/tria_boundary_lib.h>
 #include <dofs/dof_accessor.h>
 #include <dofs/dof_tools.h>
 #include <fe/fe_lib.lagrange.h>
 #include <numerics/matrices.h>
 #include <numerics/data_out.h>
 
-                                // out statt in
+                                // We will not read the grid from a
+                                // file as in the previous example,
+                                // but generate it using a function
+                                // of the library. However, we will
+                                // want to write out the locally
+                                // refined grids in each step, so we
+                                // need the following include file
+                                // instead of ``grid_in.h'':
 #include <grid/grid_out.h>
 
-#include <grid/tria_boundary_lib.h>
-
-                                //...
+                                // When using locally refined grids,
+                                // we will get so-called ``hanging
+                                // nodes''. However, the standard
+                                // finite element methods assumes
+                                // that the discrete solution spaces
+                                // be continuous, so we need to make
+                                // sure that the degrees of freedom
+                                // on hanging nodes conform to some
+                                // constraints such that the global
+                                // solution is continuous. The
+                                // following file contains a class
+                                // which is used to handle these
+                                // constraints:
 #include <dofs/dof_constraints.h>
+
+                                // Finally, we would like to use a
+                                // simple way to adaptively refine
+                                // the grid. While in general,
+                                // adaptivity is very
+                                // problem-specific, the error
+                                // indicator in the following file
+                                // often yields quite nicely adapted
+                                // grids for a wide class of
+                                // problems.
 #include <numerics/error_estimator.h>
 
 #include <fstream>
 
 
-
+                                // The main class is again almost
+                                // unchanged. Two additions, however,
+                                // are made: we have added the
+                                // ``refine'' function, which is used
+                                // to adaptively refine the grid
+                                // (instead of the global refinement
+                                // in the previous examples), and a
+                                // variable which will hold the
+                                // constraints associated to the
+                                // hanging nodes.
 template <int dim>
 class LaplaceProblem 
 {
@@ -56,7 +95,7 @@ class LaplaceProblem
     FEQ1<dim>            fe;
     DoFHandler<dim>      dof_handler;
 
-                                    // This is the only addition to
+                                    // This is the new variable in
                                     // the main class. We need an
                                     // object which holds a list of
                                     // the constraints from the
@@ -343,22 +382,164 @@ void LaplaceProblem<dim>::solve ()
 
                                 // Instead of global refinement, we
                                 // now use a slightly more elaborate
-                                // scheme.
-                                // ...
+                                // scheme. We will use the
+                                // ``KellyErrorEstimator'' class
+                                // which implements an error
+                                // estimator for the Laplace
+                                // equation; it can in principle
+                                // handle variable coefficients, but
+                                // we will not use these advanced
+                                // features, but rather use its most
+                                // simple form since we are not
+                                // interested in quantitative results
+                                // but only in a quick way to
+                                // generate locally refined grids.
+                                //
+                                // Although the error estimator
+                                // derived by Kelly et al. was
+                                // originally developed for Laplace's
+                                // equation, we have found that it is
+                                // also well suited to quickly
+                                // generate locally refined grids for
+                                // a wide class of
+                                // problems. Basically, it looks at
+                                // the jumps of the gradients of the
+                                // solution over the faces of cells
+                                // (which is a measure for the second
+                                // derivatives) and scales it by the
+                                // size of the cell. It is therefore
+                                // a measure for the local smoothness
+                                // of the solution at the place of
+                                // each cell and it is thus
+                                // understandable that it yields
+                                // reasonable grids also for
+                                // hyperbolic transport problems or
+                                // the wave equation as well,
+                                // although these grids are certainly
+                                // suboptimal compared to approaches
+                                // specially tailored to the
+                                // problem. This error estimator may
+                                // therefore be understood as a quick
+                                // way to test an adaptive program.
 template <int dim>
 void LaplaceProblem<dim>::refine_grid ()
 {
+                                  // The output of the error
+                                  // estimator class is an error
+                                  // indicator for each cell. We
+                                  // therefore need a vector with as
+                                  // many elements as there are
+                                  // active cells. Since accuracy is
+                                  // not that important here, the
+                                  // data type for the error values
+                                  // on each cell is ``float''
+                                  // instead of ``double''.
   Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
 
+                                  // Next, the error estimator can
+                                  // handle Neumann boundary
+                                  // conditions. For this, it needs
+                                  // to know which parts of the
+                                  // boundary have Neumann boundary
+                                  // conditions and teh respective
+                                  // boundary values there. This
+                                  // information is mediated by a map
+                                  // in which the keys are the
+                                  // boundary part numbers and the
+                                  // values are pointers to the
+                                  // boundary value functions. We
+                                  // create such a map, but since we
+                                  // do not use Neumann boundary
+                                  // conditions, the map will not
+                                  // contain entries.
   KellyErrorEstimator<dim>::FunctionMap neumann_boundary;
+
+                                  // Now we call the error
+                                  // estimator. The parameters should
+                                  // be clear apart from the
+                                  // quadrature formula: as said
+                                  // above, the jump of the gradients
+                                  // of the solution across the faces
+                                  // of a cell are considered. They
+                                  // are integrated along the face,
+                                  // but as usual in finite element
+                                  // programs the integration is done
+                                  // using quadrature. Since the
+                                  // error estimator class can't know
+                                  // itself which quadrature formula
+                                  // might be appropriate, we have to
+                                  // pass one to the function. Note
+                                  // that since the quadrature has to
+                                  // take place along faces, the
+                                  // dimension of the quadrature
+                                  // formula is ``dim-1'' rather then
+                                  // ``dim''.
   KellyErrorEstimator<dim>::estimate (dof_handler,
                                      QGauss3<dim-1>(),
                                      neumann_boundary,
                                      solution,
                                      estimated_error_per_cell);
 
+                                  // The above function returned one
+                                  // error indicator value for each
+                                  // cell in the
+                                  // ``estimated_error_per_cell''
+                                  // array. Refinement is now done as
+                                  // follows: refine those 30 per
+                                  // cent of the cells with the
+                                  // highest error values, and
+                                  // coarsen the 3 per cent of cells
+                                  // with the lowest values.
+                                  //
+                                  // One can easily verify that if
+                                  // the second number were zero,
+                                  // this would approximately result
+                                  // in a doubling of cells in each
+                                  // step in two space dimensions,
+                                  // since for each of the 30 per
+                                  // cent of cells four new would be
+                                  // replaced. In practice, some more
+                                  // cells are usually produced since
+                                  // it is disallowed that a cell is
+                                  // refined twice while the neighbor
+                                  // cell is not refined; in that
+                                  // case, the neighbor cell would be
+                                  // refined as well.
+                                  //
+                                  // In many applications, the number
+                                  // of cells to be coarsened would
+                                  // be set to something larger than
+                                  // only three per cent. A non-zero
+                                  // value is useful especially if
+                                  // for some reason the initial
+                                  // (coarse) grid is already rather
+                                  // refined. In that case, it might
+                                  // be necessary to refine it in
+                                  // some regions, while coarsening
+                                  // in some other regions is
+                                  // useful. In our case here, the
+                                  // initial grid is very coarse, so
+                                  // coarsening is only necessary in
+                                  // a few regions where
+                                  // over-refinement may have taken
+                                  // place. Thus a small, non-zero
+                                  // value is appropriate here.
   triangulation.refine_and_coarsen_fixed_number (estimated_error_per_cell,
                                                 0.3, 0.03);
+
+                                  // After the previous function has
+                                  // exited, some cells are flagged
+                                  // for refinement, and some other
+                                  // for coarsening. The refinement
+                                  // or coarsening itself is not
+                                  // performed by now, however, since
+                                  // there are many cases where
+                                  // further modifications of these
+                                  // flags is useful. Here, we don't
+                                  // want to do any such thing, so we
+                                  // can tell the triangulation to
+                                  // perform the actions for which
+                                  // the cells are flagged.
   triangulation.execute_coarsening_and_refinement ();
 };
 
@@ -367,7 +548,18 @@ void LaplaceProblem<dim>::refine_grid ()
 template <int dim>
 void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
 {
-                                  // ...
+                                  // We want to write the grid in
+                                  // each cycle. Here is another way
+                                  // to quickly produce a filename
+                                  // based on the cycle number. It
+                                  // assumes that the numbers `0'
+                                  // through `9' are represented
+                                  // consecutively in the character
+                                  // set (which is the case in all
+                                  // known character sets). However,
+                                  // this will only work if the cycle
+                                  // number is less than ten, which
+                                  // we check by an assertion.
   string filename = "grid-";
   filename += ('0' + cycle);
   Assert (cycle < 10, ExcInternalError());
@@ -375,12 +567,13 @@ void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
   filename += ".eps";
   ofstream output (filename.c_str());
 
+                                  // Using this filename, we write
+                                  // each grid as a postscript file.
   GridOut grid_out;
   grid_out.write_eps (triangulation, output);
 };
 
 
-
 template <int dim>
 void LaplaceProblem<dim>::run () 
 {
@@ -390,7 +583,32 @@ void LaplaceProblem<dim>::run ()
 
       if (cycle == 0)
        {
-                                          //...
+                                          // Instead of reading the
+                                          // grid from a file on disk
+                                          // as in the previous
+                                          // example, we now again
+                                          // create it using a
+                                          // library function. The
+                                          // domain is again a
+                                          // circle, which is why we
+                                          // have to provide a
+                                          // suitable boundary object
+                                          // as well.
+                                          //
+                                          // You will notice by
+                                          // looking at the coarse
+                                          // grid that it is of
+                                          // inferior quality than
+                                          // the one which we read
+                                          // from the file in the
+                                          // previous example: the
+                                          // cells are less equally
+                                          // formed. However, using
+                                          // the library function
+                                          // this program works in
+                                          // any space dimension,
+                                          // which was not the case
+                                          // before.
          GridGenerator::hyper_ball (triangulation);
 
          static const HyperBallBoundary<dim> boundary;
@@ -399,8 +617,19 @@ void LaplaceProblem<dim>::run ()
          triangulation.refine_global (1);
        }
       else
-                                        // ...
-       refine_grid ();
+                                        // In case this is not the
+                                        // first cycle, we want to
+                                        // refine the grid. Unlike
+                                        // the global refinement
+                                        // employed in the last
+                                        // example, we now use the
+                                        // adaptive procedure
+                                        // described in the function
+                                        // which we now call:
+       {
+         refine_grid ();
+       };
+      
 
       cout << "   Number of active cells: "
           << triangulation.n_active_cells()
@@ -412,7 +641,14 @@ void LaplaceProblem<dim>::run ()
       output_results (cycle);
     };
 
-                                  // ...
+                                  // The solution on the final grid
+                                  // is now written to a file. As
+                                  // already done in one of the
+                                  // previous examples, we use the
+                                  // EPS format for output, and to
+                                  // obtain a reasonable view on the
+                                  // solution, we rescale the z-axis
+                                  // by a factor of four.
   DataOut<dim>::EpsFlags eps_flags;
   eps_flags.z_scaling = 4;
   
@@ -428,13 +664,128 @@ void LaplaceProblem<dim>::run ()
 };
 
     
-
+                                // The main function is unaltered in
+                                // its functionality against the
+                                // previous example, but we have
+                                // taken a step of additional
+                                // caution. Sometimes, something goes
+                                // wrong (such as insufficient disk
+                                // space upon writing an output file,
+                                // not enough memory when trying to
+                                // allocate a vector or a matrix, or
+                                // if we can't read from or write to
+                                // a file for whatever reason), and
+                                // in these cases the library will
+                                // throw exceptions. Since they do
+                                // not constitute programming errors,
+                                // these exceptions also are not
+                                // switched off in optimized mode, in
+                                // constrast to the ``Assert'' macro
+                                // which we have used to test against
+                                // programming errors. If uncought,
+                                // these exceptions propagate the
+                                // call tree up to the ``main''
+                                // function, and if they are not
+                                // caught there either, the program
+                                // is aborted. In many cases, like if
+                                // there is not enough memory or disk
+                                // space, we can't do anything but we
+                                // can at least print some text
+                                // trying to explain the reason why
+                                // the program failed. A way to do so
+                                // is shown in the following. It is
+                                // certainly useful to write any
+                                // larger program in this way, and
+                                // you can do so by more or less
+                                // copying this function apart from
+                                // the ``try'' block which contains
+                                // the code that constitutes the
+                                // actual functionality.
 int main () 
 {
-  deallog.depth_console (0);
 
-  LaplaceProblem<2> laplace_problem_2d;
-  laplace_problem_2d.run ();
-  
+                                  // The general idea behind the
+                                  // layout of this function is as
+                                  // follows: let's try to run the
+                                  // program as we did before...
+  try
+    {
+      deallog.depth_console (0);
+
+      LaplaceProblem<2> laplace_problem_2d;
+      laplace_problem_2d.run ();
+    }
+                                  // ...and if this should fail, try
+                                  // to gather as much information as
+                                  // possible. Specifically, if the
+                                  // exception that was thrown is an
+                                  // object of a class that is
+                                  // derived from the C++ standard
+                                  // class ``exception'', then we can
+                                  // use the ``what'' member function
+                                  // to get a string which describes
+                                  // the reason why the exception was
+                                  // thrown. 
+                                  //
+                                  // The deal.II exception classes
+                                  // are all derived from the
+                                  // standard class, and in
+                                  // particular, the ``exc.what()''
+                                  // function will return
+                                  // approximately the same string as
+                                  // would be generated if the
+                                  // exception was thrown using the
+                                  // ``Assert'' macro. You have seen
+                                  // the output of such an exception
+                                  // in the previous example, and you
+                                  // then know that it contains the
+                                  // file and line number of where
+                                  // the exception occured, and some
+                                  // other information. This is also
+                                  // what would be printed in the
+                                  // following.
+  catch (exception &exc)
+    {
+      cerr << endl << endl
+          << "----------------------------------------------------"
+          << endl;
+      cerr << "Exception on processing: " << exc.what() << endl
+          << "Aborting!" << endl
+          << "----------------------------------------------------"
+          << endl;
+                                      // We can't do much more than
+                                      // printing as much information
+                                      // as we can get to, so abort
+                                      // with error:
+      return 1;
+    }
+                                  // If the exception that was thrown
+                                  // somewhere was not an object of a
+                                  // class derived from the standard
+                                  // ``exception'' class, then we
+                                  // can't do anything at all. We
+                                  // then simply print an error
+                                  // message and exit.
+  catch (...) 
+    {
+      cerr << endl << endl
+          << "----------------------------------------------------"
+          << endl;
+      cerr << "Unknown exception!" << endl
+          << "Aborting!" << endl
+          << "----------------------------------------------------"
+          << endl;
+      return 1;
+    };
+
+                                  // If we got to this point, there
+                                  // was no exception which
+                                  // propagated up to the main
+                                  // functino (maybe there were some,
+                                  // but they were caught somewhere
+                                  // in the program or the
+                                  // library). Therefore, the program
+                                  // performed as was expected and we
+                                  // can return without error.
   return 0;
 };

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.