From cec18690bf22a03ddfd563803a4331d3a37df593 Mon Sep 17 00:00:00 2001 From: wolf Date: Thu, 6 Jan 2000 15:52:28 +0000 Subject: [PATCH] Finish step-6. git-svn-id: https://svn.dealii.org/trunk@2161 0785d39b-7218-0410-832d-ea1e28bc413d --- .../examples/step-by-step/step-6/Makefile | 2 +- .../examples/step-by-step/step-6/step-6.cc | 393 +++++++++++++++++- deal.II/examples/step-6/Makefile | 2 +- deal.II/examples/step-6/step-6.cc | 393 +++++++++++++++++- 4 files changed, 746 insertions(+), 44 deletions(-) diff --git a/deal.II/deal.II/Attic/examples/step-by-step/step-6/Makefile b/deal.II/deal.II/Attic/examples/step-by-step/step-6/Makefile index dc2798f839..8753051799 100644 --- a/deal.II/deal.II/Attic/examples/step-by-step/step-6/Makefile +++ b/deal.II/deal.II/Attic/examples/step-by-step/step-6/Makefile @@ -50,7 +50,7 @@ $(target) : $(libraries) # rule how to run the program run: $(target) - @echo ============================ Running $@ + @echo ============================ Running $< @./$(target) diff --git a/deal.II/deal.II/Attic/examples/step-by-step/step-6/step-6.cc b/deal.II/deal.II/Attic/examples/step-by-step/step-6/step-6.cc index 04af2fa3d9..2c89c63110 100644 --- a/deal.II/deal.II/Attic/examples/step-by-step/step-6/step-6.cc +++ b/deal.II/deal.II/Attic/examples/step-by-step/step-6/step-6.cc @@ -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 #include #include @@ -17,6 +19,7 @@ #include #include #include +#include #include #include #include @@ -25,19 +28,55 @@ #include #include - // 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 -#include - - //... + // 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 + + // 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 #include - + // 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 class LaplaceProblem { @@ -56,7 +95,7 @@ class LaplaceProblem FEQ1 fe; DoFHandler 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::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 void LaplaceProblem::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 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::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::estimate (dof_handler, QGauss3(), 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::refine_grid () template void LaplaceProblem::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::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 void LaplaceProblem::run () { @@ -390,7 +583,32 @@ void LaplaceProblem::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 boundary; @@ -399,8 +617,19 @@ void LaplaceProblem::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::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::EpsFlags eps_flags; eps_flags.z_scaling = 4; @@ -428,13 +664,128 @@ void LaplaceProblem::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; }; diff --git a/deal.II/examples/step-6/Makefile b/deal.II/examples/step-6/Makefile index dc2798f839..8753051799 100644 --- a/deal.II/examples/step-6/Makefile +++ b/deal.II/examples/step-6/Makefile @@ -50,7 +50,7 @@ $(target) : $(libraries) # rule how to run the program run: $(target) - @echo ============================ Running $@ + @echo ============================ Running $< @./$(target) diff --git a/deal.II/examples/step-6/step-6.cc b/deal.II/examples/step-6/step-6.cc index 04af2fa3d9..2c89c63110 100644 --- a/deal.II/examples/step-6/step-6.cc +++ b/deal.II/examples/step-6/step-6.cc @@ -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 #include #include @@ -17,6 +19,7 @@ #include #include #include +#include #include #include #include @@ -25,19 +28,55 @@ #include #include - // 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 -#include - - //... + // 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 + + // 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 #include - + // 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 class LaplaceProblem { @@ -56,7 +95,7 @@ class LaplaceProblem FEQ1 fe; DoFHandler 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::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 void LaplaceProblem::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 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::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::estimate (dof_handler, QGauss3(), 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::refine_grid () template void LaplaceProblem::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::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 void LaplaceProblem::run () { @@ -390,7 +583,32 @@ void LaplaceProblem::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 boundary; @@ -399,8 +617,19 @@ void LaplaceProblem::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::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::EpsFlags eps_flags; eps_flags.z_scaling = 4; @@ -428,13 +664,128 @@ void LaplaceProblem::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; }; -- 2.39.5