From 3c22202f92f477e849e1463d6801e5827d0d5c1d Mon Sep 17 00:00:00 2001 From: wolf Date: Thu, 26 Jan 2006 06:52:55 +0000 Subject: [PATCH] Review the documentation. git-svn-id: https://svn.dealii.org/trunk@12173 0785d39b-7218-0410-832d-ea1e28bc413d --- .../step-3.data/results.html | 24 +- deal.II/examples/step-3/step-3.cc | 880 +++++++++--------- 2 files changed, 486 insertions(+), 418 deletions(-) diff --git a/deal.II/doc/tutorial/chapter-2.step-by-step/step-3.data/results.html b/deal.II/doc/tutorial/chapter-2.step-by-step/step-3.data/results.html index 39219e4a43..836212f226 100644 --- a/deal.II/doc/tutorial/chapter-2.step-by-step/step-3.data/results.html +++ b/deal.II/doc/tutorial/chapter-2.step-by-step/step-3.data/results.html @@ -8,12 +8,11 @@ Number of active cells: 1024 Total number of cells: 1365 Number of degrees of freedom: 1089 DEAL:cg::Starting value 0.121094 -DEAL:cg::Starting value 0.323075 -DEAL:cg::Convergence step 47 value 5.33692e-13 +DEAL:cg::Convergence step 48 value 5.33692e-13 The first three lines is what we wrote to cout. The last -three lines were generated without our interaction by the CG +two lines were generated without our intervention by the CG solver. The first two lines state the residual at the start of the iteration, while the last line tells us that the solver needed 47 iterations to bring the norm of the residual to 5.3e-13, i.e. below @@ -29,6 +28,25 @@ Apart from the output shown above, the program generated the file viewed as follows: invoke GNUPLOT and enter the following sequence of commands at its prompt:

+examples/step-3> gnuplot
+
+        G N U P L O T
+        Version 3.7 patchlevel 3
+        last modified Thu Dec 12 13:00:00 GMT 2002
+        System: Linux 2.6.11.4-21.10-default
+
+        Copyright(C) 1986 - 1993, 1998 - 2002
+        Thomas Williams, Colin Kelley and many others
+
+        Type `help` to access the on-line reference manual
+        The gnuplot FAQ is available from
+        http://www.gnuplot.info/gnuplot-faq.html
+
+        Send comments and requests for help to 
+        Send bugs, suggestions and mods to 
+
+
+Terminal type set to 'x11'
 gnuplot> set data style lines
 gnuplot> splot "solution.gpl"
 
diff --git a/deal.II/examples/step-3/step-3.cc b/deal.II/examples/step-3/step-3.cc index 071dcc9350..3d9991a818 100644 --- a/deal.II/examples/step-3/step-3.cc +++ b/deal.II/examples/step-3/step-3.cc @@ -4,35 +4,42 @@ /* $Id$ */ /* Version: $Name$ */ /* */ -/* Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004 by the deal.II authors */ +/* Copyright (C) 1999, 2000, 2001, 2002, 2003, 2004, 2006 by the deal.II authors */ /* */ /* This file is subject to QPL and may not be distributed */ /* without copyright and license information. Please refer */ /* to the file deal.II/doc/license.html for the text and */ /* further information on this license. */ + + // @sect3{Many new include files} + // These include files are already // known to you. They declare the // classes which handle - // triangulations and enumerate the - // degrees of freedom. + // triangulations and enumeration of + // degrees of freedom: #include #include // And this is the file in which the - // functions are declared which - // create grids. + // functions are declared that + // create grids: #include - // The next three files contain - // classes which are needed for loops - // over all cells and to get the - // information from the cell objects. + // The next three files contain classes which + // are needed for loops over all cells and to + // get the information from the cell + // objects. The first two have been used + // before to get geometric information from + // cells; the last one is new and provides + // information about the degrees of freedom + // local to a cell: #include #include #include - // In this file are the finite - // element descriptions. + // In this file contains the description of + // the Lagrange interpolation finite element: #include // And this file is needed for the @@ -45,7 +52,7 @@ // assembling the matrix using // quadrature on each cell. The // classes declared in them will be - // explained below. + // explained below: #include #include @@ -56,25 +63,21 @@ #include #include - // These include files are for the - // linear algebra which we employ to - // solve the system of equations - // arising from the finite element - // discretization of the Laplace - // equation. We will use vectors and - // full matrices for assembling the - // system of equations locally on - // each cell, and transfer the - // results into a sparse matrix. We - // will then use a Conjugate Gradient - // solver to solve the problem, for - // which we need a preconditioner (in - // this program, we use the identity - // preconditioner which does nothing, - // but we need to include the file - // anyway), and a class which - // provides the solver with some - // memory for temporary vectors. + // We're now almost to the end. The second to + // last group of include files is for the + // linear algebra which we employ to solve + // the system of equations arising from the + // finite element discretization of the + // Laplace equation. We will use vectors and + // full matrices for assembling the system of + // equations locally on each cell, and + // transfer the results into a sparse + // matrix. We will then use a Conjugate + // Gradient solver to solve the problem, for + // which we need a preconditioner (in this + // program, we use the identity + // preconditioner which does nothing, but we + // need to include the file anyway): #include #include #include @@ -88,35 +91,39 @@ #include - // Instead of the procedural - // programming of previous examples, - // we encapsulate everything into a - // class for this program. The class - // consists of functions which do - // certain aspects of a finite - // element program, a `main' function - // which controls what is done first - // and what is done next, and a list - // of member variables. + // @sect3{The ``LaplaceProblem'' class} + + // Instead of the procedural programming of + // previous examples, we encapsulate + // everything into a class for this + // program. The class consists of functions + // which each perform certain aspects of a + // finite element program, a `main' function + // which controls what is done first and what + // is done next, and a list of member + // variables. + + // The public part of the class is rather + // short: it has a constructor and a function + // `run' that is called from the outside and + // acts as something like the `main' + // function: it coordinates which operations + // of this class shall be run in which + // order. Everything else in the class, + // i.e. all the functions that actually do + // anything, are in the private section of + // the class: class LaplaceProblem { public: - // This is the constructor: LaplaceProblem (); - // And the top-level function, - // which is called from the - // outside to start the whole - // program (see the `main' - // function at the bottom of this - // file): void run (); - // Then there are some member - // functions that mostly do what - // their names suggest. Since - // they do not need to be called - // from outside, they are made + // Then there are the member functions + // that mostly do what their names + // suggest. Since they do not need to be + // called from outside, they are made // private to this class. private: void make_grid_and_dofs (); @@ -124,10 +131,10 @@ class LaplaceProblem void solve (); void output_results () const; - // And then we have the member + // And finally we have some member // variables. There are variables // describing the triangulation - // and the numbering of the + // and the global numbering of the // degrees of freedom (we will // specify the exact polynomial // degree of the finite element @@ -152,27 +159,35 @@ class LaplaceProblem Vector system_rhs; }; - - // Here comes the constructor. It - // does not much more than first to - // specify that we want bi-linear - // elements (denoted by the parameter - // to the finite element object, - // which specifies the polynomial - // degree), and to associate the - // dof_handler variable to the - // triangulation we use. All the - // other member variables of the - // LaplaceProblem class have a - // default constructor which does all - // we want. + // @sect4{LaplaceProblem::LaplaceProblem} + + // Here comes the constructor. It does not + // much more than first to specify that we + // want bi-linear elements (denoted by the + // parameter to the finite element object, + // which indicates the polynomial degree), + // and to associate the dof_handler variable + // to the triangulation we use. (Note that + // the triangulation isn't set up with a mesh + // at all at the present time, but the + // DoFHandler doesn't care: it only wants to + // know which triangulation it will be + // associated with, and it only starts to + // care about an actual mesh once you try to + // distribute degree of freedom on the mesh + // using the distribute_dofs() function.) All + // the other member variables of the + // LaplaceProblem class have a default + // constructor which does all we want. LaplaceProblem::LaplaceProblem () : fe (1), dof_handler (triangulation) {} - // Now, the first thing we've got to + // @sect4{LaplaceProblem::make_grid_and_dofs} + + // Now, the first thing we've got to // do is to generate the // triangulation on which we would // like to do our computation and @@ -194,48 +209,43 @@ void LaplaceProblem::make_grid_and_dofs () // total of 1024. GridGenerator::hyper_cube (triangulation, -1, 1); triangulation.refine_global (5); - // Unsure that 1024 is the correct - // number? Let's see: - // n_active_cells return the number - // of terminal cells. By terminal - // we mean the cells on the finest - // grid. + // Unsure that 1024 is the correct number? + // Let's see: n_active_cells return the + // number of active cells: std::cout << "Number of active cells: " << triangulation.n_active_cells() << std::endl; - // We stress the adjective - // `terminal' or `active', since - // there are more cells, namely the - // parent cells of the finest - // cells, their parents, etc, up to - // the one cell which made up the - // initial grid. Of course, on the - // next coarser level, the number - // of cells is one quarter that of - // the cells on the finest level, - // i.e. 256, then 64, 16, 4, and - // 1. We can get the total number - // of cells like this: + // Here, by active we mean the cells on the + // finest level, i.e. cells that aren't + // refined any further. We stress the + // adjective `active', since there are more + // cells, namely the parent cells of the + // finest cells, their parents, etc, up to + // the one cell which made up the initial + // grid. Of course, on the next coarser + // level, the number of cells is one + // quarter that of the cells on the finest + // level, i.e. 256, then 64, 16, 4, and + // 1. We can get the total number of cells + // like this: std::cout << "Total number of cells: " << triangulation.n_cells() << std::endl; // Note the distinction between // n_active_cells() and n_cells(). - // Next we enumerate all the - // degrees of freedom. This is done - // by using the distribute_dofs - // function, as we have seen in - // previous examples. Since we use - // the FEQ1 class, i.e. bilinear - // elements, this associates one - // degree of freedom with each - // vertex. + // Next we enumerate all the degrees of + // freedom. This is done by using the + // distribute_dofs function, as we have + // seen in the step-2 example. Since we use + // the ``FE_Q'' class with a polynomial + // degree of 1, i.e. bilinear elements, + // this associates one degree of freedom + // with each vertex. While we're at + // generating output, let us also take a + // look at how many degrees of freedom are + // generated: dof_handler.distribute_dofs (fe); - - // Now that we have the degrees of - // freedom, we can take a look at - // how many there are: std::cout << "Number of degrees of freedom: " << dof_handler.n_dofs() << std::endl; @@ -244,13 +254,16 @@ void LaplaceProblem::make_grid_and_dofs () // 32 grid, the number of DoFs // should be 33 times 33, or 1089. - // As we have seen in the previous - // example, we set up a sparse - // matrix for the system matrix and - // tag those entries that might be - // nonzero. Since that has already - // been done, we won't discuss the - // next few lines: + // As we have seen in the previous example, + // we set up a sparsity pattern for the + // system matrix and tag those entries that + // might be nonzero. Compared to what we + // did in step-2, the only change is that + // instead of giving a magically obtained + // maximal number of nonzero entries per + // row, we now use a function in the + // ``DoFHandler'' class that can compute + // this number for us: sparsity_pattern.reinit (dof_handler.n_dofs(), dof_handler.n_dofs(), dof_handler.max_couplings_between_dofs()); @@ -262,7 +275,7 @@ void LaplaceProblem::make_grid_and_dofs () // `compress' has been called, you // can't add nonzero entries // anymore; the sparsity pattern is - // `sealed', so to say), and we can + // `sealed', so to say), we can // initialize the matrix itself // with it. Note that the // SparsityPattern object does @@ -296,6 +309,8 @@ void LaplaceProblem::make_grid_and_dofs () system_rhs.reinit (dof_handler.n_dofs()); } + // @sect4{LaplaceProblem::assemble_system} + // Now comes the difficult part: // assembling matrices and @@ -306,132 +321,178 @@ void LaplaceProblem::make_grid_and_dofs () // other things in the functions // above and below. // - // The general way to assemble - // matrices and vectors is to loop - // over all cells, and on each cell - // compute the contribution of that - // cell to the global matrix and - // right hand side by quadrature. The - // idea now is that since we only - // need the finite element shape - // functions on the quadrature points - // of each cell, we don't need the - // shape functions of the finite - // element themselves any - // more. Therefore, we won't deal - // with the finite element object - // `fe' (which was of type FEQ1), but - // with another object which only - // provides us with the values, - // gradients, etc of the shape - // functions at the quadrature - // points. The objects which do this - // are of type FEValues. + // The general way to assemble matrices and + // vectors is to loop over all cells, and on + // each cell compute the contribution of that + // cell to the global matrix and right hand + // side by quadrature. The point to realize + // now is that we need the values of the + // shape functions at the locations of + // quadrature points on the real + // cell. However, both the finite element + // shape functions as well as the quadrature + // points are only defined on the unit + // cell. They are therefore of little help to + // us, and we will in fact hardly ever query + // information about finite element shape + // functions or quadrature points from these + // objects directly. + // + // Rather, what is required is a way to map + // this data from the unit cell to the real + // cell. Classes that can do that are derived + // from the Mapping class, though one again + // often does not have to deal with them + // directly: many functions in the library + // can take a mapping object as argument, but + // when it is omitted they simply resort to + // the standard bilinear Q1 mapping. We will + // go this route, and not bother with it for + // the moment (we come back to this in + // step-10, step-11, and step-12). + // + // So what we now have is a collection of + // three classes to deal with: finite + // element, quadrature, and mapping + // objects. That's too much, so there is one + // type of class that orchestrates + // information exchange between these three: + // the ``FEValues'' class. If given one + // instance of each three of these objects, + // it will be able to provide you with + // information about values and gradients of + // shape functions at quadrature points on a + // real cell. + // + // Using all this, we will assemble the + // linear system for this problem in the + // following function: void LaplaceProblem::assemble_system () { - // Ok, let's start: we need a - // quadrature formula for the - // evaluation of the integrals on - // each cell. Let's take a Gauss - // formula with two quadrature - // points in each direction, i.e. a - // total of four points since we - // are in 2D. This quadrature - // formula integrates polynomials - // of degrees up to three exactly - // (in 1D). Since the integrands in - // the matrix entries are quadratic - // (in 1D), this is sufficient. The - // same holds for 2D. + // Ok, let's start: we need a quadrature + // formula for the evaluation of the + // integrals on each cell. Let's take a + // Gauss formula with two quadrature points + // in each direction, i.e. a total of four + // points since we are in 2D. This + // quadrature formula integrates + // polynomials of degrees up to three + // exactly (in 1D). It is easy to check + // that this is sufficient for the present + // problem: QGauss<2> quadrature_formula(2); - // And we initialize the object - // which we have briefly talked - // about above. It needs to be told - // which the finite element is that - // we want to use, the quadrature - // points and their - // weights. Finally, we have to - // tell it what we want it to - // compute on each cell: we need - // the values of the shape - // functions at the quadrature - // points, their gradients, and - // also the weights of the - // quadrature points and the + // And we initialize the object which we + // have briefly talked about above. It + // needs to be told which finite element we + // want to use, and the quadrature points + // and their weights (jointly described by + // a Quadrature object). As mentioned, we + // use the implied Q1 mapping, rather than + // specifying one ourselves + // explicitly. Finally, we have to tell it + // what we want it to compute on each cell: + // we need the values of the shape + // functions at the quadrature points (for + // the right hand side (f,phi)), their + // gradients (for the matrix entries (grad + // phi_i, grad phi_j)), and also the + // weights of the quadrature points and the // determinants of the Jacobian - // transformations from the unit - // cell to the real cells. The - // values of the shape functions - // computed by specifying - // update_values; the gradients are - // done alike, using - // update_gradients. The - // determinants of the Jacobians - // and the weights are always used - // together, so only the products - // (Jacobians times weights, or - // short JxW) are computed; since - // we also need them, we have to - // list them as well. The advantage - // of this proceeding is that we - // calculate only what we - // need. This optimatizes the - // process of solving: + // transformations from the unit cell to + // the real cells. + // + // This list of what kind of information we + // actually need is given as a bitwise + // connection of flags as the third + // argument to the constructor of + // ``FEValues''. Since these values have to + // be recomputed, or updated, every time we + // go to a new cell, all of these flags + // start with the prefix ``update_'' and + // then indicate what it actually is that + // we want updated. The flag to give if we + // want the values of the shape functions + // computed is ``update_values''; for the + // gradients it is + // ``update_gradients''. The determinants + // of the Jacobians and the quadrature + // weights are always used together, so + // only the products (Jacobians times + // weights, or short ``JxW'') are computed; + // since we need them, we have to list + // ``update_JxW_values'' as well: FEValues<2> fe_values (fe, quadrature_formula, - UpdateFlags(update_values | - update_gradients | - update_JxW_values)); - - // For use further down below, we - // define two short cuts for the - // number of degrees of freedom on - // each cell (since we are in 2D - // and degrees of freedom are - // associated with vertices only, - // this number is four). We also - // define an abbreviation for the - // number of quadrature points - // (here that should be nine). In - // general, it is a good idea to - // use their symbolic names instead - // of hard-coding these number even - // if you know them, since you may - // want to change the quadrature - // formula and/or finite element at - // some time; the program will just - // work with these changes, without - // the need to change the matrix - // assemblage. + update_values | update_gradients | update_JxW_values); + // The advantage of this proceeding is that + // we can specify what kind of information + // we actually need on each cell. It is + // easily understandable that this approach + // can significant speed up finite element + // computations, compared to approaches + // where everything, including second + // derivatives, normal vectors to cells, + // etc are computed on each cell, + // regardless whether they are needed or + // not. + + // For use further down below, we define + // two short cuts for values that will be + // used very frequently. First, an + // abbreviation for the number of degrees + // of freedom on each cell (since we are in + // 2D and degrees of freedom are associated + // with vertices only, this number is four, + // but we rather want to write the + // definition of this variable in a way + // that does not preclude us from later + // choosing a different finite element that + // has a different number of degrees of + // freedom per cell, or work in a different + // space dimension). + // + // Secondly, we also define an abbreviation + // for the number of quadrature points + // (here that should be four). In general, + // it is a good idea to use their symbolic + // names instead of hard-coding these + // number even if you know them, since you + // may want to change the quadrature + // formula and/or finite element at some + // time; the program will just work with + // these changes, without the need to + // change anything in this function. // - // The shortcuts, finally, are only - // defined to make the following - // loops a bit more readable. You - // will see them in many places in - // larger programs, and - // `dofs_per_cell' and `n_q_points' - // are more or less standard names - // for these purposes. + // The shortcuts, finally, are only defined + // to make the following loops a bit more + // readable. You will see them in many + // places in larger programs, and + // `dofs_per_cell' and `n_q_points' are + // more or less by convention the standard + // names for these purposes: const unsigned int dofs_per_cell = fe.dofs_per_cell; const unsigned int n_q_points = quadrature_formula.n_quadrature_points; - // Now, we said that we wanted to - // assemble the global matrix and - // vector cell-by-cell. We could - // write the results directly into - // the global matrix, but this is - // not very efficient since access - // to the elements of a sparse - // matrix is slow. Rather, we first - // compute the contribution of each - // ell in a small matrix with the - // degrees of freedom on the - // present cell, and only transfer - // them to the global matrix when - // the computations are finished - // for this cell. We do the same - // for the right hand side vector, - // although access times are not so - // problematic for them. + // Now, we said that we wanted to assemble + // the global matrix and vector + // cell-by-cell. We could write the results + // directly into the global matrix, but + // this is not very efficient since access + // to the elements of a sparse matrix is + // slow. Rather, we first compute the + // contribution of each cell in a small + // matrix with the degrees of freedom on + // the present cell, and only transfer them + // to the global matrix when the + // computations are finished for this + // cell. We do the same for the right hand + // side vector. So let's first allocate + // these objects (these being local + // objects, all degrees of freedom are + // coupling with all others, and we should + // use a full matrix object rather than a + // sparse one for the local operations; + // everything will be transferred to a + // global sparse matrix later on): FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); Vector cell_rhs (dofs_per_cell); @@ -445,72 +506,70 @@ void LaplaceProblem::assemble_system () // we transfer the result into the // global matrix, we have to know // the global numbers of the - // degrees of freedom. When we get + // degrees of freedom. When we query // them, we need a scratch // (temporary) array for these // numbers: std::vector local_dof_indices (dofs_per_cell); - // Now for the loop over all - // cells. You have seen before how - // this works, so this should be - // familiar to you: - DoFHandler<2>::active_cell_iterator cell = dof_handler.begin_active(), - endc = dof_handler.end(); + // Now for the loop over all cells. We have + // seen before how this works, so this + // should be familiar including the + // conventional names for these variables: + DoFHandler<2>::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); for (; cell!=endc; ++cell) { - // We are on one cell, and we - // would like the values and - // gradients of the shape - // functions be computed, as - // well as the determinants of - // the Jacobian matrices of the - // mapping between unit cell - // and true cell, at the - // quadrature points. Since all - // these values depend on the - // geometry of the cell, we - // have to have the FEValues - // object re-compute them on - // each cell: + // We are now sitting on one cell, and + // we would like the values and + // gradients of the shape functions be + // computed, as well as the + // determinants of the Jacobian + // matrices of the mapping between unit + // cell and true cell, at the + // quadrature points. Since all these + // values depend on the geometry of the + // cell, we have to have the FEValues + // object re-compute them on each cell: fe_values.reinit (cell); - // Reset the values of the - // contributions of this cell - // to global matrix and global - // right hand side to zero, - // before we fill them. + // Next, reset the local cell's + // contributions contributions to + // global matrix and global right hand + // side to zero, before we fill them: cell_matrix = 0; cell_rhs = 0; - // Assemble the matrix: For the - // Laplace problem, the matrix - // on each cell is the integral - // over the gradients of shape - // function i and j. Since we - // do not integrate, but rather - // use quadrature, this is the - // sum over all quadrature - // points of the integrands - // times the determinant of the - // Jacobian matrix at the - // quadrature point times the - // weight of this quadrature - // point. You can get the - // gradient of shape function i - // at quadrature point q_point - // by using + // Then finally assemble the matrix: + // For the Laplace problem, the matrix + // on each cell is the integral over + // the gradients of shape function i + // and j. Since we do not integrate, + // but rather use quadrature, this is + // the sum over all quadrature points + // of the integrands times the + // determinant of the Jacobian matrix + // at the quadrature point times the + // weight of this quadrature point. You + // can get the gradient of shape + // function i at quadrature point + // q_point by using // fe_values.shape_grad(i,q_point); - // this gradient is a - // 2-dimensional vector (in - // fact it is of type - // Tensor<1,dim>, with here - // dim=2) and the product of - // two such vectors is the - // scalar product, i.e. the - // product of the two - // shape_grad function calls is - // the dot product. + // this gradient is a 2-dimensional + // vector (in fact it is of type + // Tensor<1,dim>, with here dim=2) and + // the product of two such vectors is + // the scalar product, i.e. the product + // of the two shape_grad function calls + // is the dot product. This is in turn + // multiplied by the Jacobian + // determinant and the quadrature point + // weight (that one gets together by + // the call to + // ``fe_values.JxW''). Finally, this is + // repeated for all shape functions + // phi_i and phi_j: for (unsigned int i=0; i boundary_values; VectorTools::interpolate_boundary_values (dof_handler, 0, @@ -686,7 +732,9 @@ void LaplaceProblem::assemble_system () } - // The following function simply + // @sect4{LaplaceProblem::solve} + + // The following function simply // solves the discretized // equation. As the system is quite a // large one for direct solvers such @@ -697,42 +745,42 @@ void LaplaceProblem::assemble_system () // variables here (only 1089) is a // very small number for finite // element computations, where - // 100.000 is a more usual number; - // for this number of variables, + // 100.000 is a more usual number. + // For this number of variables, // direct methods are no longer // usable and you are forced to use // methods like CG. void LaplaceProblem::solve () { - // We need to tell the algorithm - // where to stop. This is done by - // using a SolverControl object, - // and as stopping criterion we - // say: maximally 1000 iterations - // (which is far more than is - // needed for 1089 variables; see - // the results section to find out - // how many were really used), and - // stop if the norm of the residual - // is below 1e-12. In practice, the - // latter criterion will be the one - // which stops the iteration. + // First, we need to have an object that + // knows how to tell the CG algorithm when + // to stop. This is done by using a + // ``SolverControl'' object, and as + // stopping criterion we say: stop after a + // maximum of 1000 iterations (which is far + // more than is needed for 1089 variables; + // see the results section to find out how + // many were really used), and stop if the + // norm of the residual is below 1e-12. In + // practice, the latter criterion will be + // the one which stops the iteration: SolverControl solver_control (1000, 1e-12); - // Then we need the solver - // itself. The template parameters - // here are the matrix type and the - // type of the vectors, but the - // empty angle brackets indicate - // that we simply take the default - // arguments. + // Then we need the solver itself. The + // template parameters to the ``SolverCG'' + // class are the matrix type and the type + // of the vectors, but the empty angle + // brackets indicate that we simply take + // the default arguments (which are + // ``SparseMatrix'' and + // ``Vector''): SolverCG<> cg (solver_control); - // Now solve the system of - // equations. The CG solver takes a - // preconditioner, but we don't - // want to use one, so we tell it - // to use the identity operation as - // preconditioner. + // Now solve the system of equations. The + // CG solver takes a preconditioner as its + // fourth argument. We don't feel ready to + // delve into this yet, so we tell it to + // use the identity operation as + // preconditioner: cg.solve (system_matrix, solution, system_rhs, PreconditionIdentity()); // Now that the solver has done its @@ -742,6 +790,8 @@ void LaplaceProblem::solve () } + // @sect4{LaplaceProblem::output_results} + // The last part of a typical finite // element program is to output the // results and maybe do some @@ -757,23 +807,20 @@ void LaplaceProblem::output_results () const // To write the output to a file, // we need an object which knows // about output formats and the - // like. This is the DataOut class, + // like. This is the ``DataOut'' class, // and we need an object of that // type: DataOut<2> data_out; - // Now we have to tell it where to - // take the values from which it - // shall write. We tell it which - // DoFHandler object to use, and we - // add the solution vector (and the - // name by which it shall be - // written to disk) to the list of - // data that is to be written. If - // we had more than one vector - // which we would like to look at - // in the output (for example right - // hand sides, errors per cell, - // etc) we would add them as well: + // Now we have to tell it where to take the + // values from which it shall write. We + // tell it which ``DoFHandler'' object to + // use, and the solution vector (and + // the name by which the solution variable + // shall appear in the output file). If + // we had more than one vector which we + // would like to look at in the output (for + // example right hand sides, errors per + // cell, etc) we would add them as well: data_out.attach_dof_handler (dof_handler); data_out.add_data_vector (solution, "solution"); // After the DataOut object knows @@ -783,9 +830,9 @@ void LaplaceProblem::output_results () const // handle. The reason is that we // have separated the frontend // (which knows about how to treat - // DoFHandler objects and data + // ``DoFHandler'' objects and data // vectors) from the back end (which - // knows several output formats) + // knows many different output formats) // and use an intermediate data // format to transfer data from the // front- to the backend. The data @@ -806,15 +853,16 @@ void LaplaceProblem::output_results () const } - // The following function is the main - // function which calls all the other - // functions of the LaplaceProblem - // class. The order in which this is - // done resembles the order in which - // most finite element programs - // work. Since the names are mostly - // self-explanatory, there is not - // much to comment about: + // @sect4{LaplaceProblem::run} + + // Finally, the last function of this class + // is the main function which calls all the + // other functions of the ``LaplaceProblem'' + // class. The order in which this is done + // resembles the order in which most finite + // element programs work. Since the names are + // mostly self-explanatory, there is not much + // to comment about: void LaplaceProblem::run () { make_grid_and_dofs (); @@ -823,7 +871,8 @@ void LaplaceProblem::run () output_results (); } - + + // @sect3{The ``main'' function} // This is the main function of the // program. Since the concept of a @@ -832,11 +881,12 @@ void LaplaceProblem::run () // programming, it often does not // much more than creating an object // of the top-level class and calling - // it principle function. This is - // what is done here as well. + // its principle function. This is + // what is done here as well: int main () { LaplaceProblem laplace_problem; laplace_problem.run (); + return 0; } -- 2.39.5