From: wolf Date: Mon, 20 Dec 1999 14:54:08 +0000 (+0000) Subject: Mostly finish step-3 doc. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e95f4aaa88f75b7699538647649e757781528597;p=dealii-svn.git Mostly finish step-3 doc. git-svn-id: https://svn.dealii.org/trunk@2089 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/Attic/examples/step-by-step/step-3/step-3.cc b/deal.II/deal.II/Attic/examples/step-by-step/step-3/step-3.cc index e9e5732286..eb81074355 100644 --- a/deal.II/deal.II/Attic/examples/step-by-step/step-3/step-3.cc +++ b/deal.II/deal.II/Attic/examples/step-by-step/step-3/step-3.cc @@ -49,57 +49,344 @@ class LaplaceProblem }; + // Here comes the constructor. It + // does not much more than 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. LaplaceProblem::LaplaceProblem () : dof_handler (triangulation) {}; - + // Now, the first thing we've got to + // do is to generate the + // triangulation on which we would + // like to do our computation and + // number each vertex with a degree + // of freedom. We have seen this in + // the previous examples before. Then + // we have to set up space for the + // system matrix and right hand side + // of the discretized problem. This + // is what this function does: void LaplaceProblem::make_grid_and_dofs () { + // First create the grid and refine + // all cells five times. Since the + // initial grid (which is the + // square [-1,1]x[-1,1]) consists + // of only one cell, the final grid + // has 32 times 32 cells, for a + // 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. + cout << "Number of active cells: " + << triangulation.n_active_cells() + << 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 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: + cout << "Total number of cells: " + << triangulation.n_cells() + << 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. dof_handler.distribute_dofs (fe); + // Now that we have the degrees of + // freedom, we can take a look at + // how many there are: + cout << "Number of degrees of freedom: " + << dof_handler.n_dofs() + << endl; + // There should be one DoF for each + // vertex. Since we have a 32 times + // 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: sparsity_pattern.reinit (dof_handler.n_dofs(), dof_handler.n_dofs(), dof_handler.max_couplings_between_dofs()); DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern); sparsity_pattern.compress(); + // Now the sparsity pattern is + // built and fixed (after + // `compress' has been called, you + // can't add nonzero entries + // anymore; the sparsity pattern is + // `sealed', so to say), and we can + // initialize the matrix itself + // with it. Note that the + // SparseMatrixStruct object does + // not hold the values of the + // matrix, it only stores the + // places where entries are. The + // entries are themselves stored in + // objects of type SparseMatrix, of + // which our variable system_matrix + // is one. + // + // The distinction between sparsity + // pattern and matrix was made to + // allow several matrices to use + // the same sparsity pattern. This + // may not seem relevant, but when + // you consider the size which + // matrices can have, and that it + // may take some time to build the + // sparsity pattern, this becomes + // important in large-scale + // problems. system_matrix.reinit (sparsity_pattern); + // The last thing to do in this + // function is to set the sizes of + // the right hand side vector and + // the solution vector to the right + // values: solution.reinit (dof_handler.n_dofs()); system_rhs.reinit (dof_handler.n_dofs()); }; - + // Now comes the difficult part: + // assembling matrices and + // vectors. In fact, this is not + // overly difficult, but it is + // something that the library can't + // do for you as for most of the + // 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. 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 three quadrature + // points in each direction, i.e. a + // total of nine points since we + // are in 2D: QGauss3<2> quadrature_formula; + // 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 + // determinants of the Jacobian + // transformations from the unit + // cell to the real cells. The + // values of the shape functions + // are always computed, so we don't + // have to list them. 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: FEValues<2> fe_values (fe, quadrature_formula, UpdateFlags(update_gradients | update_JxW_values)); - const unsigned int n_q_points = quadrature_formula.n_quadrature_points; + // 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. + // + // 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. 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 copmutations are finished + // for this cell. We do the same + // for the right hand side vector, + // although access times are not so + // problematic for them. FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); Vector cell_rhs (dofs_per_cell); - + + // When assembling the + // contributions of each cell, we + // do this with the local numbering + // of the degrees of freedom + // (i.e. the number running from + // zero through + // dofs_per_cell-1). However, when + // we transfer the result into the + // global matrix, we have to know + // the global numbers of the + // degrees of freedom. When we get + // them, we need a scratch array + // for these numbers: vector local_dof_indices (dofs_per_cell); + // Now for th 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(); 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: 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. cell_matrix.clear (); cell_rhs.clear (); - + + // 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. for (unsigned int i=0; iget_dof_indices (local_dof_indices); + + // Then again loop over all + // shape functions i and j and + // transfer the local elements + // to the global matrix. The + // global numbers can be + // obtained using + // local_dof_indices[i]: for (unsigned int i=0; i boundary_values; VectorTools::interpolate_boundary_values (dof_handler, - 0, - ZeroFunction<2>(), - boundary_values); + 0, + ZeroFunction<2>(), + boundary_values); + // Now that we got the list of + // boundary DoFs and their + // respective boundary values, + // let's use them to modify the + // system of equations + // accordingly. This is done by the + // following function call: MatrixTools<2>::apply_boundary_values (boundary_values, system_matrix, solution, @@ -136,30 +558,151 @@ void LaplaceProblem::assemble_system () }; - + // The following function simply + // solves the discretized + // equation. As the system is quite a + // large one for direct solvers such + // as Gauss elimination or LU + // decomposition, we use a Conjugate + // Gradient algorithm. You should + // remember that the number of + // 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, + // 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. SolverControl solver_control (1000, 1e-12); + // Furthermore, the CG algorithm + // needs some space for temporary + // vectors. Rather than allocating + // it on the stack or heap itself, + // it relies on helper objects, + // which can sometimes do a better + // job at this. The + // PrimitiveVectorMemory class is + // such a helper class which the + // solver can ask for memory. The + // angle brackets indicate that + // this class really takes a + // template parameter (here the + // data type of the vectors we + // use), which however has a + // default value, which is + // appropriate here. PrimitiveVectorMemory<> vector_memory; + // Then we need the solver + // itself. The template parameters + // here are the matrix type and the + // type of the vectors. They + // default to the ones we use here. SolverCG<> cg (solver_control, vector_memory); - + + // 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. cg.solve (system_matrix, solution, system_rhs, PreconditionIdentity()); + // Now that the solver has done its + // job, the solution variable + // contains the nodal values of the + // solution function. }; + // The last part of a typical finite + // element program is to output the + // results and maybe do some + // postprocessing (for example + // compute the maximal stress values + // at the boundary, or the average + // flux across the outflow, etc). We + // have no such postprocessing here, + // but we would like to write the + // solution to a file. void LaplaceProblem::output_results () { + // To write the output to a file, + // we need an object which knows + // about output formats and the + // 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: data_out.attach_dof_handler (dof_handler); data_out.add_data_vector (solution, "solution"); + // After the DataOut object knows + // which data it is to work on, we + // have to tell it to process them + // into something the backends can + // handle. The reason is that we + // have separated the frontend + // (which knows about how to treat + // DoFHandler objects and data + // vectors) from the backend (which + // knows several output formats) + // and use an intermediate data + // format to transfer data from the + // front- to the backend. The data + // is transformed into this + // intermediate format by the + // following function: data_out.build_patches (); - + + // Now we have everything in place + // for the actual output. Just open + // a file and write the data into + // it, using GNUPLOT format (there + // are other functions which write + // their data in postscript, AVS, + // GMV, or some other format): ofstream output ("solution.gpl"); data_out.write_gnuplot (output); }; + // 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: void LaplaceProblem::run () { make_grid_and_dofs(); @@ -170,6 +713,15 @@ void LaplaceProblem::run () + // This is the main function of the + // program. Since the concept of a + // main function is mostly a remnant + // from the pre-object era in C/C++ + // 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. int main () { LaplaceProblem laplace_problem; diff --git a/deal.II/examples/step-3/step-3.cc b/deal.II/examples/step-3/step-3.cc index e9e5732286..eb81074355 100644 --- a/deal.II/examples/step-3/step-3.cc +++ b/deal.II/examples/step-3/step-3.cc @@ -49,57 +49,344 @@ class LaplaceProblem }; + // Here comes the constructor. It + // does not much more than 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. LaplaceProblem::LaplaceProblem () : dof_handler (triangulation) {}; - + // Now, the first thing we've got to + // do is to generate the + // triangulation on which we would + // like to do our computation and + // number each vertex with a degree + // of freedom. We have seen this in + // the previous examples before. Then + // we have to set up space for the + // system matrix and right hand side + // of the discretized problem. This + // is what this function does: void LaplaceProblem::make_grid_and_dofs () { + // First create the grid and refine + // all cells five times. Since the + // initial grid (which is the + // square [-1,1]x[-1,1]) consists + // of only one cell, the final grid + // has 32 times 32 cells, for a + // 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. + cout << "Number of active cells: " + << triangulation.n_active_cells() + << 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 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: + cout << "Total number of cells: " + << triangulation.n_cells() + << 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. dof_handler.distribute_dofs (fe); + // Now that we have the degrees of + // freedom, we can take a look at + // how many there are: + cout << "Number of degrees of freedom: " + << dof_handler.n_dofs() + << endl; + // There should be one DoF for each + // vertex. Since we have a 32 times + // 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: sparsity_pattern.reinit (dof_handler.n_dofs(), dof_handler.n_dofs(), dof_handler.max_couplings_between_dofs()); DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern); sparsity_pattern.compress(); + // Now the sparsity pattern is + // built and fixed (after + // `compress' has been called, you + // can't add nonzero entries + // anymore; the sparsity pattern is + // `sealed', so to say), and we can + // initialize the matrix itself + // with it. Note that the + // SparseMatrixStruct object does + // not hold the values of the + // matrix, it only stores the + // places where entries are. The + // entries are themselves stored in + // objects of type SparseMatrix, of + // which our variable system_matrix + // is one. + // + // The distinction between sparsity + // pattern and matrix was made to + // allow several matrices to use + // the same sparsity pattern. This + // may not seem relevant, but when + // you consider the size which + // matrices can have, and that it + // may take some time to build the + // sparsity pattern, this becomes + // important in large-scale + // problems. system_matrix.reinit (sparsity_pattern); + // The last thing to do in this + // function is to set the sizes of + // the right hand side vector and + // the solution vector to the right + // values: solution.reinit (dof_handler.n_dofs()); system_rhs.reinit (dof_handler.n_dofs()); }; - + // Now comes the difficult part: + // assembling matrices and + // vectors. In fact, this is not + // overly difficult, but it is + // something that the library can't + // do for you as for most of the + // 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. 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 three quadrature + // points in each direction, i.e. a + // total of nine points since we + // are in 2D: QGauss3<2> quadrature_formula; + // 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 + // determinants of the Jacobian + // transformations from the unit + // cell to the real cells. The + // values of the shape functions + // are always computed, so we don't + // have to list them. 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: FEValues<2> fe_values (fe, quadrature_formula, UpdateFlags(update_gradients | update_JxW_values)); - const unsigned int n_q_points = quadrature_formula.n_quadrature_points; + // 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. + // + // 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. 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 copmutations are finished + // for this cell. We do the same + // for the right hand side vector, + // although access times are not so + // problematic for them. FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); Vector cell_rhs (dofs_per_cell); - + + // When assembling the + // contributions of each cell, we + // do this with the local numbering + // of the degrees of freedom + // (i.e. the number running from + // zero through + // dofs_per_cell-1). However, when + // we transfer the result into the + // global matrix, we have to know + // the global numbers of the + // degrees of freedom. When we get + // them, we need a scratch array + // for these numbers: vector local_dof_indices (dofs_per_cell); + // Now for th 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(); 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: 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. cell_matrix.clear (); cell_rhs.clear (); - + + // 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. for (unsigned int i=0; iget_dof_indices (local_dof_indices); + + // Then again loop over all + // shape functions i and j and + // transfer the local elements + // to the global matrix. The + // global numbers can be + // obtained using + // local_dof_indices[i]: for (unsigned int i=0; i boundary_values; VectorTools::interpolate_boundary_values (dof_handler, - 0, - ZeroFunction<2>(), - boundary_values); + 0, + ZeroFunction<2>(), + boundary_values); + // Now that we got the list of + // boundary DoFs and their + // respective boundary values, + // let's use them to modify the + // system of equations + // accordingly. This is done by the + // following function call: MatrixTools<2>::apply_boundary_values (boundary_values, system_matrix, solution, @@ -136,30 +558,151 @@ void LaplaceProblem::assemble_system () }; - + // The following function simply + // solves the discretized + // equation. As the system is quite a + // large one for direct solvers such + // as Gauss elimination or LU + // decomposition, we use a Conjugate + // Gradient algorithm. You should + // remember that the number of + // 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, + // 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. SolverControl solver_control (1000, 1e-12); + // Furthermore, the CG algorithm + // needs some space for temporary + // vectors. Rather than allocating + // it on the stack or heap itself, + // it relies on helper objects, + // which can sometimes do a better + // job at this. The + // PrimitiveVectorMemory class is + // such a helper class which the + // solver can ask for memory. The + // angle brackets indicate that + // this class really takes a + // template parameter (here the + // data type of the vectors we + // use), which however has a + // default value, which is + // appropriate here. PrimitiveVectorMemory<> vector_memory; + // Then we need the solver + // itself. The template parameters + // here are the matrix type and the + // type of the vectors. They + // default to the ones we use here. SolverCG<> cg (solver_control, vector_memory); - + + // 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. cg.solve (system_matrix, solution, system_rhs, PreconditionIdentity()); + // Now that the solver has done its + // job, the solution variable + // contains the nodal values of the + // solution function. }; + // The last part of a typical finite + // element program is to output the + // results and maybe do some + // postprocessing (for example + // compute the maximal stress values + // at the boundary, or the average + // flux across the outflow, etc). We + // have no such postprocessing here, + // but we would like to write the + // solution to a file. void LaplaceProblem::output_results () { + // To write the output to a file, + // we need an object which knows + // about output formats and the + // 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: data_out.attach_dof_handler (dof_handler); data_out.add_data_vector (solution, "solution"); + // After the DataOut object knows + // which data it is to work on, we + // have to tell it to process them + // into something the backends can + // handle. The reason is that we + // have separated the frontend + // (which knows about how to treat + // DoFHandler objects and data + // vectors) from the backend (which + // knows several output formats) + // and use an intermediate data + // format to transfer data from the + // front- to the backend. The data + // is transformed into this + // intermediate format by the + // following function: data_out.build_patches (); - + + // Now we have everything in place + // for the actual output. Just open + // a file and write the data into + // it, using GNUPLOT format (there + // are other functions which write + // their data in postscript, AVS, + // GMV, or some other format): ofstream output ("solution.gpl"); data_out.write_gnuplot (output); }; + // 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: void LaplaceProblem::run () { make_grid_and_dofs(); @@ -170,6 +713,15 @@ void LaplaceProblem::run () + // This is the main function of the + // program. Since the concept of a + // main function is mostly a remnant + // from the pre-object era in C/C++ + // 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. int main () { LaplaceProblem laplace_problem;